tutorial.R 7.77 KB
Newer Older
1
2
3
4
5
6
7
8
## ----setup, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval=FALSE
)

## ----Input / Output files------------------------------------------------
jbferet's avatar
jbferet committed
9
#  Input_Image_File  = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'biodivMapR')
10
#  
jbferet's avatar
jbferet committed
11
12
13
14
#  # Input.Image.File  = raster2BIL(Raster.Path = Input.Image.File,
#  #                                        Sensor = 'SENTINEL_2A',
#  #                                        Convert.Integer = TRUE,
#  #                                        Output.Directory = '~/test')
15
#  
jbferet's avatar
jbferet committed
16
17
18
#  Input_Mask_File   = FALSE
#  
#  Output_Dir        = 'RESULTS'
19
20

## ----Spatial resolution--------------------------------------------------
jbferet's avatar
jbferet committed
21
#  window_size = 10
22
23

## ----PCA filtering-------------------------------------------------------
24
#  FilterPCA = FALSE
25
26

## ----Computing options---------------------------------------------------
jbferet's avatar
jbferet committed
27
#  nbCPU         = 2
28
29
30
31
#  MaxRAM        = 0.5
#  nbclusters    = 50

## ----Mask non vegetated / shaded / cloudy pixels-------------------------
jbferet's avatar
jbferet committed
32
33
34
#  NDVI_Thresh = 0.5
#  Blue_Thresh = 500
#  NIR_Thresh  = 1500
35
#  print("PERFORM RADIOMETRIC FILTERING")
jbferet's avatar
jbferet committed
36
37
38
#  Input_Mask_File = perform_radiometric_filtering(Input_Image_File, Input_Mask_File, Output_Dir,
#                                              NDVI_Thresh = NDVI_Thresh, Blue_Thresh = Blue_Thresh,
#                                              NIR_Thresh = NIR_Thresh)
39
40
41

## ----PCA-----------------------------------------------------------------
#  print("PERFORM PCA ON RASTER")
jbferet's avatar
jbferet committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#  PCA_Output        = perform_PCA(Input_Image_File, Input_Mask_File, Output_Dir,
#                                 FilterPCA = TRUE, nbCPU = nbCPU,MaxRAM = MaxRAM)
#  # path for the PCA raster
#  PCA_Files         = PCA_Output$PCA_Files
#  # number of pixels used for each partition used for k-means clustering
#  Pix_Per_Partition = PCA_Output$Pix_Per_Partition
#  # number of partitions used for k-means clustering
#  nb_partitions     = PCA_Output$nb_partitions
#  # path for the updated mask
#  Input_Mask_File   = PCA_Output$MaskPath
#  # parameters of the PCA model
#  PCA_model         = PCA_Output$PCA_model
#  # definition of spectral bands to be excluded from the analysis
#  SpectralFilter    = PCA_Output$SpectralFilter
#  
57
#  print("Select PCA components for diversity estimations")
jbferet's avatar
jbferet committed
58
#  select_PCA_components(Input_Image_File, Output_Dir, PCA_Files, File_Open = TRUE)
59

60
## ----Spectral species map------------------------------------------------
61
#  print("MAP SPECTRAL SPECIES")
jbferet's avatar
jbferet committed
62
63
#  map_spectral_species(Input_Image_File, Output_Dir, PCA_Files,PCA_model, SpectralFilter, Input_Mask_File,
#                       Pix_Per_Partition, nb_partitions, nbCPU=nbCPU, MaxRAM=MaxRAM)
64
65

## ----alpha and beta diversity maps---------------------------------------
66
67
#  print("MAP ALPHA DIVERSITY")
#  # Index.Alpha   = c('Shannon','Simpson')
jbferet's avatar
jbferet committed
68
69
70
#  Index_Alpha   = c('Shannon')
#  map_alpha_div(Input_Image_File, Output_Dir, window_size,
#                nbCPU=nbCPU, MaxRAM=MaxRAM, Index_Alpha = Index_Alpha)
71
72
#  
#  print("MAP BETA DIVERSITY")
jbferet's avatar
jbferet committed
73
74
#  map_beta_div(Input_Image_File, Output_Dir, window_size, nb_partitions=nb_partitions,
#               nbCPU=nbCPU, MaxRAM=MaxRAM)
75

Florian de Boissieu's avatar
Florian de Boissieu committed
76
## ----alpha and beta diversity indices from vector layer------------------
77
78
79
80
81
82
83
#  # location of the spectral species raster needed for validation
#  TypePCA     = 'SPCA'
#  Dir.Raster  = file.path(Output.Dir,basename(Input.Image.File),TypePCA,'SpectralSpecies')
#  Name.Raster = 'SpectralSpecies'
#  Path.Raster = file.path(Dir.Raster,Name.Raster)
#  
#  # location of the directory where shapefiles used for validation are saved
Florian de Boissieu's avatar
Florian de Boissieu committed
84
#  vect        = system.file('extdata', 'VECTOR', package = 'biodivMapR')
85
86
87
#  Shannon.All = list() # ??
#  
#  # list vector data
jbferet's avatar
jbferet committed
88
#  Path.Vector         = list_shp(vect)
89
90
91
#  Name.Vector         = tools::file_path_sans_ext(basename(Path.Vector))
#  
#  # get alpha and beta diversity indicators corresponding to shapefiles
jbferet's avatar
jbferet committed
92
#  Biodiv.Indicators           = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#  # if no name
#  Biodiv.Indicators$Name.Plot = seq(1,length(Biodiv.Indicators$Shannon[[1]]),by = 1)
#  Shannon.RS                  = c(Biodiv.Indicators$Shannon)[[1]]

## ----Write validation----------------------------------------------------
#  # write RS indicators
#  ####################################################
#  # write indicators for alpha diversity
#  Path.Results = file.path(Output.Dir, basename(Input.Image.File), TypePCA, 'VALIDATION')
#  dir.create(Path.Results, showWarnings = FALSE, recursive = TRUE)
#  ShannonIndexFile <- file.path(Path.Results, "ShannonIndex.tab")
#  write.table(Shannon.RS, file = ShannonIndexFile, sep = "\t", dec = ".", na = " ",
#              row.names = Biodiv.Indicators$Name.Plot, col.names= F, quote=FALSE)
#  
#  Results =  data.frame(Name.Vector, Biodiv.Indicators$Richness, Biodiv.Indicators$Fisher,                                Biodiv.Indicators$Shannon, Biodiv.Indicators$Simpson)
#  names(Results)  = c("ID_Plot", "Species_Richness", "Fisher", "Shannon", "Simpson")
#  write.table(Results, file = paste(Path.Results,"AlphaDiversity.tab",sep=''), sep="\t", dec=".",               na=" ", row.names = F, col.names= T,quote=FALSE)
#  
#  # write indicators for beta diversity
#  BC_mean = Biodiv.Indicators$BCdiss
#  colnames(BC_mean) = rownames(BC_mean) = Biodiv.Indicators$Name.Plot
#  write.table(BC_mean, file = paste(Path.Results,"BrayCurtis.csv",sep=''), sep="\t", dec=".", na=" ", row.names = F, col.names= T,quote=FALSE)
#  

117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
## ----PCoA on Field Plots-------------------------------------------------
#  # apply ordination using PCoA (same as done for map_beta_div)
#  library(labdsv)
#  MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
#  BetaPCO   = pco(MatBCdist, k = 3)
#  

## ----plot PCoA & Shannon-------------------------------------------------
#  # very uglily assign vegetation type to polygons in shapefiles
#  nbSamples = c(6,4,7,7)
#  vg        = c('Forest high diversity', 'Forest low diversity', 'Forest medium diversity', 'low vegetation')
#  Type_Vegetation = c()
#  for (i in 1: length(nbSamples)){
#    for (j in 1:nbSamples[i]){
#      Type_Vegetation = c(Type_Vegetation,vg[i])
#    }
#  }
#  
#  # create data frame including alpha and beta diversity
#  library(ggplot2)
#  Results     =  data.frame('vgtype'=Type_Vegetation,'pco1'= BetaPCO$points[,1],'pco2'= BetaPCO$points[,2],'pco3' = BetaPCO$points[,3],'shannon'=Shannon.RS)
#  
#  # plot field data in the PCoA space, with size corresponding to shannon index
#  ggplot(Results, aes(x=pco1, y=pco2, color=vgtype,size=shannon)) +
#    geom_point(alpha=0.6) +
#    scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
#  filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA2.png')
#  ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
#         scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
#         dpi = 600, limitsize = TRUE)
#  
#  
#  ggplot(Results, aes(x=pco1, y=pco3, color=vgtype,size=shannon)) +
#    geom_point(alpha=0.6) +
#    scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
#  filename = file.path(Path.Results,'BetaDiversity_PcoA1_vs_PcoA3.png')
#  ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
#         scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
#         dpi = 600, limitsize = TRUE)
#  
#  ggplot(Results, aes(x=pco2, y=pco3, color=vgtype,size=shannon)) +
#    geom_point(alpha=0.6) +
#    scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
#  filename = file.path(Path.Results,'BetaDiversity_PcoA2_vs_PcoA3.png')
#  ggsave(filename, plot = last_plot(), device = 'png', path = NULL,
#         scale = 1, width = NA, height = NA, units = c("in", "cm", "mm"),
#         dpi = 600, limitsize = TRUE)
#