example_script.R 9.68 KB
Newer Older
Florian de Boissieu's avatar
Florian de Boissieu committed
1
# ===============================================================================
Florian de Boissieu's avatar
Florian de Boissieu committed
2
# biodivMapR
Florian de Boissieu's avatar
Florian de Boissieu committed
3
4
5
6
7
8
9
# ===============================================================================
# PROGRAMMERS:
#
# Jean-Baptiste FERET <jb.feret@irstea.fr>
#
# Copyright 2019/06 Jean-Baptiste FERET
#
Florian de Boissieu's avatar
Florian de Boissieu committed
10
# biodivMapR is free software: you can redistribute it and/or modify
Florian de Boissieu's avatar
Florian de Boissieu committed
11
12
13
14
15
16
17
18
19
20
21
22
23
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>
#
# ===============================================================================
24
25
library(raster)
library(biodivMapR)
Florian de Boissieu's avatar
Florian de Boissieu committed
26
27
28
29
30
################################################################################
##              DEFINE PARAMETERS FOR DATASET TO BE PROCESSED                 ##
################################################################################
# path (absolute or relative) for the image to process
# expected to be in ENVI HDR format, BIL interleaved
Florian de Boissieu's avatar
Florian de Boissieu committed
31
Input.Image.File  = system.file('extdata', 'RASTER', 'S2A_T33NUD_20180104_Subset', package = 'biodivMapR')
32
check_data(Input.Image.File)
Florian de Boissieu's avatar
Florian de Boissieu committed
33

34
35
36
37
38
# # convert the image using Convert.Raster2BIL if not in the proper format
# Input.Image.File  = raster2BIL(Raster.Path = Input.Image.File,
#                                        Sensor = 'SENTINEL_2A',
#                                        Convert.Integer = TRUE,
#                                        Output.Directory = '~/test')
Florian de Boissieu's avatar
Florian de Boissieu committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

# full path for the Mask raster corresponding to image to process
# expected to be in ENVI HDR format, 1 band, integer 8bits
# expected values in the raster: 0 = masked, 1 = selected
# set to FALSE if no mask available
Input.Mask.File   = FALSE

# relative or absolute path for the Directory where results will be stored
# For each image processed, a subdirectory will be created after its name
Output.Dir        = 'RESULTS'

# SPATIAL RESOLUTION
# resolution of spatial units for alpha and beta diversity maps (in pixels), relative to original image
# if Res.Map = 10 for images with 10 m spatial resolution, then spatial units will be 10 pixels x 10m = 100m x 100m surfaces
# rule of thumb: spatial units between 0.25 and 4 ha usually match with ground data
jbferet's avatar
jbferet committed
54
55
# too small window_size results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image
window_size       = 10
Florian de Boissieu's avatar
Florian de Boissieu committed
56
57

# PCA FILTERING: 		Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower
58
FilterPCA         = FALSE
Florian de Boissieu's avatar
Florian de Boissieu committed
59

jbferet's avatar
jbferet committed
60
61
62
63
64
# type of PCA:
# PCA: no rescaling of the data
# SPCA: rescaling of the data
TypePCA     = 'SPCA'

Florian de Boissieu's avatar
Florian de Boissieu committed
65
66
67
################################################################################
##      Check if the image format is compatible with codes (ENVI BIL)         ##
################################################################################
68
check_data(Input.Image.File)
Florian de Boissieu's avatar
Florian de Boissieu committed
69
70
71
72

################################################################################
##                    DEFINE PARAMETERS FOR METHOD                            ##
################################################################################
73
nbCPU         = 2
Florian de Boissieu's avatar
Florian de Boissieu committed
74
75
76
77
78
79
80
81
82
83
84
MaxRAM        = 0.5
nbclusters    = 50

################################################################################
##                              PROCESS IMAGES                                ##
################################################################################
# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
NDVI.Thresh = 0.5
Blue.Thresh = 500
NIR.Thresh  = 1500
print("PERFORM RADIOMETRIC FILTERING")
jbferet's avatar
jbferet committed
85
ImPathShade         = perform_radiometric_filtering(Input.Image.File,Input.Mask.File,Output.Dir,
Florian de Boissieu's avatar
Florian de Boissieu committed
86
87
88
89
90
                                                    NDVI.Thresh = NDVI.Thresh, Blue.Thresh = Blue.Thresh,
                                                    NIR.Thresh = NIR.Thresh)

# 2- Compute PCA for a random selection of pixels in the raster
print("PERFORM PCA ON RASTER")
91
PCA.Files           = perform_PCA(Input.Image.File,ImPathShade,Output.Dir,FilterPCA=FilterPCA,nbCPU=nbCPU,MaxRAM = MaxRAM)
Florian de Boissieu's avatar
Florian de Boissieu committed
92
93

# 3- Select principal components from the PCA raster
jbferet's avatar
jbferet committed
94
select_PCA_components(Input.Image.File,Output.Dir,PCA.Files,File.Open = TRUE)
Florian de Boissieu's avatar
Florian de Boissieu committed
95
96
97
98
99
100

################################################################################
##                      MAP ALPHA AND BETA DIVERSITY                          ##
################################################################################

print("MAP SPECTRAL SPECIES")
jbferet's avatar
jbferet committed
101
map_spectral_species(Input.Image.File,Output.Dir,PCA.Files,nbCPU=nbCPU,MaxRAM=MaxRAM)
Florian de Boissieu's avatar
Florian de Boissieu committed
102
103
104
105
106


print("MAP ALPHA DIVERSITY")
# Index.Alpha   = c('Shannon','Simpson')
Index.Alpha   = c('Shannon')
jbferet's avatar
jbferet committed
107
map_alpha_div(Input.Image.File,Output.Dir,window_size,nbCPU=nbCPU,MaxRAM=MaxRAM,Index.Alpha = Index.Alpha)
Florian de Boissieu's avatar
Florian de Boissieu committed
108
109

print("MAP BETA DIVERSITY")
jbferet's avatar
jbferet committed
110
map_beta_div(Input.Image.File,Output.Dir,window_size,nbCPU=nbCPU,MaxRAM=MaxRAM)
Florian de Boissieu's avatar
Florian de Boissieu committed
111
112
113
114
115
116
117
118
119
120
121

################################################################################
##          COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS                 ##
################################################################################

# location of the spectral species raster needed for validation
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
122
vect        = system.file('extdata', 'VECTOR', package = 'biodivMapR')
Florian de Boissieu's avatar
Florian de Boissieu committed
123
124
125
Shannon.All = list()

# list vector data
126
Path.Vector         = list_shp(vect)
Florian de Boissieu's avatar
Florian de Boissieu committed
127
128
129
130
Name.Vector         = tools::file_path_sans_ext(basename(Path.Vector))

# read raster data including projection
RasterStack         = stack(Path.Raster)
131
Projection.Raster   = get_projection(Path.Raster,'raster')
Florian de Boissieu's avatar
Florian de Boissieu committed
132
133

# get alpha and beta diversity indicators corresponding to shapefiles
134
Biodiv.Indicators           = diversity_from_plots(Raster = Path.Raster, Plots = Path.Vector,NbClusters = nbclusters)
Florian de Boissieu's avatar
Florian de Boissieu committed
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
# if no name
Biodiv.Indicators$Name.Plot = seq(1,length(Biodiv.Indicators$Shannon[[1]]),by = 1)
Shannon.RS                  = c(Biodiv.Indicators$Shannon)[[1]]

####################################################
# write RS indicators
####################################################
# write indicators for alpha diversity
Path.Results = paste(Output.Dir,'/',basename(Input.Image.File),'/',TypePCA,'/VALIDATION/',sep='')
dir.create(Path.Results, showWarnings = FALSE,recursive = TRUE)
write.table(Shannon.RS, file = paste(Path.Results,"ShannonIndex.csv",sep=''), 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.csv",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)
jbferet's avatar
jbferet committed
155
156
157
158
159


####################################################
# illustrate results
####################################################
160
# apply ordination using PCoA (same as done for map_beta_div)
161
library(labdsv)
162
163
164
MatBCdist = as.dist(BC_mean, diag = FALSE, upper = FALSE)
BetaPCO   = pco(MatBCdist, k = 3)

165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
# 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)) +
jbferet's avatar
jbferet committed
190
191
  geom_point(alpha=0.6) +
  scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
192
193
194
195
196
197
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)) +
jbferet's avatar
jbferet committed
198
199
  geom_point(alpha=0.6) +
  scale_color_manual(values=c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
200
201
202
203
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)