-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathDrought_SVI.Rmd
412 lines (286 loc) · 18.4 KB
/
Drought_SVI.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
<<<<<<< HEAD
---
title: ''
author: "UN-SPIDER"
date: "September 25, 2018"
output:
html_notebook: default
html_document: default
---
##Recommended Practice: Drought monitoring using the Standardized Vegetation Index (SVI) - R script
#### ToDo (data preparation):
1. Download MOD13Q1 data at https://lpdaacsvc.cr.usgs.gov/appeears/
2. Create two folders, where you will store all data, one for the EVI files and one for the Pixel Reliability files,
e.g. C:/Data_EVI and C:/Data_Pixel_Reliability.
3. In both folders for each day of the year (DOY) create one subfolder. The name of each folder must start with "DOY_", e.g. DOY_033.
4. Rename the files by addding a prefix following the pattern DOY_YYYY_, e.g. 033_2001 or 001_2005.
Total Commander is a useful tool to rename multiple files (Download Link: http://www.ghisler.com/index.htm).
Renaming the files is important to automatize the filenames and the titles of the resulting maps.
6. Sort the EVI and Pixel Reliability data according to the DOY in the respective folders
7. Create another subfolder within the EVI data folder called "shape". Store the shapefile with the country border here.
8. Adjust the capture of the resulting jpg map to fit your area of interest in line 116
****
Installing relevant packages
****
```{r eval=FALSE}
install.packages("raster") # you only have to do this once
install.packages("rgdal") # you only have to do this once
install.packages("sp") # you only have to do this once
```
Calling packages for use
```{r eval=FALSE}
library(raster)
library(rgdal)
library(sp)
```
Load borders (Download country borders as shapefiles: http://www.gadm.org/download)
**ToDo: insert link to the shapefile with the country borders**
```{r eval=FALSE}
border <- shapefile("C:/User/MOD13Q1/Data/shape/gadm36_GTM_0.shp")
```
**ToDo: enter link to the folder where you have stored the MODIS EVI data**
```{r eval=FALSE}
path <- "C:/User/MOD13Q1/Data/Data_EVI"
dlist <- dir(path,pattern="DOY")
```
**ToDo: enter link to the folder where you have stored the MODIS Pixel Reliability data**
```{r eval=FALSE}
path_c <- "C:/User/MOD13Q1/Data/Data_Pixel_Reliability"
dlist_c <- dir(path_c,pattern="DOY")
```
Showing an example of the downloaded EVI data
<center>
</center>
Showing an example of the corresponding Pixel Reliability
</center>
The pixel reliability rank includes the following:
Rank Key|Summary Quality Assurance
--------|--------
-1 |Fill/No Data
***0*** |***Good data***
***1*** |***Marginal Data***
2 |Snow/Ice
3 |Cloudy
**ToDo: enter the links to the folder where you want to store the resulting .jpg-images and .tif-files.**
```{r eval=FALSE}
path_jpg <- "C:/User/MOD13Q1/Data/VCI_Maps_Guatemala_jpg"
path_tif <- "C:/User/MOD13Q1/Data/VCI_Maps_Guatemala_jpg"
```
Creating a progress bar in the Console, wich ends at the end of the loop. The progress bar looks like this:
"=======""
```{r eval=FALSE}
pb <- txtProgressBar (min=0, max=length(dlist), style=1)
setTxtProgressBar (pb, 0)
```
****
> Main code for processing the EVI data, masking out clouds, calculating the SVI and writing the results
****
```{r eval=FALSE}
for (i in 1:length(dlist)) { # start of the outer for-loop
fold <- paste(path,dlist[i],sep="/") # the respective DOY folder of the Data_EVI folder
fls <- dir(fold,pattern=".tif") # all files that are available in the respective EVI DOY folder
flsp <-paste(fold,fls,sep="/") # all files that are available in the respective EVI DOY folder with complete path name
evistack <- stack(flsp) # creates a layer stack of all files within the EVI DOY folder
eviresize<- crop(evistack,border) # resizes the EVI layer stack to the rectangular extent of the border shapefile
evimask<-mask(eviresize,border) # masks the EVI layer stack using the border shapefile
evi<-evimask*0.0001 # rescaling of MODIS EVI data
evi[evi==-0.3]<-NA # EVI fill value(-0,3) in NA
evi[evi<(0)]<-NA # as we are only interested in vegetation valid EVI range is 0 to 1 and all EVI values smaller than 0 set to NA
fold_c <- paste(path_c,dlist_c[i],sep="/") # the respective DOY folder of the Data_Pixel_Reliability folder
fls_c <- dir(fold_c,pattern=".tif") # all files that are available in the respective Pixel Reliability DOY folder
flsp_c <-paste(fold_c,fls_c,sep="/") # all files that are available in the respective Pixel Reliability DOY folder with complete path name
cloudstack <- stack(flsp_c) # creates a layer stack of all files within the Pixel Relaibility DOY folder
cloudresize<- crop(cloudstack,border) # resizes the Pixel Reliability layer stack to the rectangular extent of the border shapefile
cloudmask<-mask(cloudresize,border) # masks the Pixel Reliability layer stack using the border shapefile
cloudmask[cloudmask==(3)]<-NA # Pixel Reliability rank 3 pixels (cloudy) set to NA
cloudmask[cloudmask==(2)]<-NA # Pixel Reliability rank 2 pixels (snow&ice) set to NA
cloudmask[cloudmask==(0)]<-1 # Pixel Reliability rank 0 pixels (good quality) set to 1
cloudmask[cloudmask>(3)]<-NA # as valid Pixel Reliability range is -1 to 3, all Pixel Reliability values >3 set to NA
# (as -1 rank pixels show value 255)
evi_c=evi*cloudmask # multiplying the EVI layer stack by the Pixel Reliability layer stack
# to get one single layer stack with applied cloud mask
# extracting mean and standard deviation for each pixel
evimean <- stackApply (evi_c, rep (1, nlayers (evi)), mean, na.rm=T) #calculating the mean for the layer stack for each individual pixel
evisd <- stackApply (evi_c, rep (1, nlayers (evi)), sd, na.rm=T) #calculating the standard deviation for the layer stack for each individual pixel
# if na.rm is TRUE NA values are ignored, otherwise an NA value in any of the arguments will cause a value of NA to be returned,
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/Extremes.html
SVI_all <- ((evi_c-evimean)/evisd) #calculating SVI
# define breaks and color scheme
stdev <- cellStats(SVI_all, stat='sd') # calculating the standard deviation for each layer of the raster data,
# returning as many values as number of input layers
stdev2 <- 2*stdev # 2 standard deviations, returning as many values as number of input layers
stdev15 <- 1.5*stdev # 1.5 standard deviations, returning as many values as number of input layers
breaksSVI <- c(-4, -stdev2, -stdev15, -stdev, stdev, stdev15, stdev2, 4) # definition of the breaks of the resulting maps based on statistics
# (standard deviation), returning as many values as number of input layers
# times 6 (i.e. 6 values for each of the 15 layers) plus 2
# (the -4 in the beginning and 4 in the end)
l <- nlayers(SVI_all)
x <- seq(2, 6*l+1, l)
colorsSVI<-c('#4C0E0D', '#E72223', '#F19069', '#F9F6C6', '#64B14B', '#04984A', '#00320E') # definition of the color scheme of the resulting maps
fold_jpg <- paste(path_jpg) # the respective folder where you want to store the resulting .jpg-images.
fold_tif <- paste(path_tif) # the respective folder where you want to store the resulting .tif-files.
for (k in 1:nlayers(SVI_all)) { # start of the inner for-loop
year <- substr(fls[k],5,8) # extracting the fifth to eigths letter of the filename, which is the year (cf. data preparation above)
doy <- substr(fls[k],1,3) # extracting the first to third letter of the filename, which is the DOY (cf. data preparation above)
# writeRaster(evi[[k]], filename=paste(fold,"/",doy,"_",year,sep=""), format="ENVI", datatype='FLT4S', overwrite=TRUE)
# in case you would like to have Envi files (Attention: note the datatype)
jpeg(filename=paste(fold_jpg,"/",doy,"_",year,".jpg",sep=""), quality = 100)
# writes the jpg maps and names the files autmatically according to the pattern DOY_YYYY
plot(SVI_all[[k]], # zlim=c(0,100),
breaks=c(-4, breaksSVI[x], 4), # sets the breaks as defined above
col=colorsSVI, # sets the colors as defined above
main=paste("SVI"," (EVI)"," sample ",doy," ",year,sep="")) # automizes the title of the plot.
# ToDo: Adjust the file naming according to the data you are processing!
# E.g. if you base your SVI on EVI data, write (EVI)
dev.off()
writeRaster(SVI_all[[k]], filename=paste(fold_tif,"/",doy,"_",year,".tif",sep=""), format="GTiff", overwrite=TRUE)
# writes the geotiff and automizes the file naming according to the pattern DOY_YYYY
} # end of the inner for-loop
setTxtProgressBar (pb, i)
} # end of the outer for-loop
=======
---
title: ''
author: "UN-SPIDER"
date: "September 25, 2018"
output:
html_notebook: default
html_document: default
---
##Recommended Practice: Drought monitoring using the Standardized Vegetation Index (SVI) - R script
#### ToDo (data preparation):
1. Download MOD13Q1 data at https://lpdaacsvc.cr.usgs.gov/appeears/
2. Create two folders, where you will store all data, one for the EVI files and one for the Pixel Reliability files,
e.g. C:/Data_EVI and C:/Data_Pixel_Reliability.
3. In both folders for each day of the year (DOY) create one subfolder. The name of each folder must start with "DOY_", e.g. DOY_033.
4. Rename the files by addding a prefix following the pattern DOY_YYYY_, e.g. 033_2001 or 001_2005.
Total Commander is a useful tool to rename multiple files (Download Link: http://www.ghisler.com/index.htm).
Renaming the files is important to automatize the filenames and the titles of the resulting maps.
6. Sort the EVI and Pixel Reliability data according to the DOY in the respective folders
7. Create another subfolder within the EVI data folder called "shape". Store the shapefile with the country border here.
8. Adjust the capture of the resulting jpg map to fit your area of interest in line 116
****
Installing relevant packages
****
```{r eval=FALSE}
install.packages("raster") # you only have to do this once
install.packages("rgdal") # you only have to do this once
install.packages("sp") # you only have to do this once
```
Calling packages for use
```{r eval=FALSE}
library(raster)
library(rgdal)
library(sp)
```
Load borders (Download country borders as shapefiles: http://www.gadm.org/download)
**ToDo: insert link to the shapefile with the country borders**
```{r eval=FALSE}
border <- shapefile("C:/User/MOD13Q1/Data/shape/gadm36_GTM_0.shp")
```
**ToDo: enter link to the folder where you have stored the MODIS EVI data**
```{r eval=FALSE}
path <- "C:/User/MOD13Q1/Data/Data_EVI"
dlist <- dir(path,pattern="DOY")
```
**ToDo: enter link to the folder where you have stored the MODIS Pixel Reliability data**
```{r eval=FALSE}
path_c <- "C:/User/MOD13Q1/Data/Data_Pixel_Reliability"
dlist_c <- dir(path_c,pattern="DOY")
```
Showing an example of the downloaded EVI data
<center>
</center>
Showing an example of the corresponding Pixel Reliability
</center>
The pixel reliability rank includes the following:
Rank Key|Summary Quality Assurance
--------|--------
-1 |Fill/No Data
***0*** |***Good data***
***1*** |***Marginal Data***
2 |Snow/Ice
3 |Cloudy
**ToDo: enter the links to the folder where you want to store the resulting .jpg-images and .tif-files.**
```{r eval=FALSE}
path_jpg <- "C:/User/MOD13Q1/Data/VCI_Maps_Guatemala_jpg"
path_tif <- "C:/User/MOD13Q1/Data/VCI_Maps_Guatemala_jpg"
```
Creating a progress bar in the Console, wich ends at the end of the loop. The progress bar looks like this:
"=======""
```{r eval=FALSE}
pb <- txtProgressBar (min=0, max=length(dlist), style=1)
setTxtProgressBar (pb, 0)
```
****
> Main code for processing the EVI data, masking out clouds, calculating the SVI and writing the results
****
```{r eval=FALSE}
for (i in 1:length(dlist)) { # start of the outer for-loop
fold <- paste(path,dlist[i],sep="/") # the respective DOY folder of the Data_EVI folder
fls <- dir(fold,pattern=".tif") # all files that are available in the respective EVI DOY folder
flsp <-paste(fold,fls,sep="/") # all files that are available in the respective EVI DOY folder with complete path name
evistack <- stack(flsp) # creates a layer stack of all files within the EVI DOY folder
eviresize<- crop(evistack,border) # resizes the EVI layer stack to the rectangular extent of the border shapefile
evimask<-mask(eviresize,border) # masks the EVI layer stack using the border shapefile
evi<-evimask*0.0001 # rescaling of MODIS EVI data
evi[evi==-0.3]<-NA # EVI fill value(-0,3) in NA
evi[evi<(0)]<-NA # as we are only interested in vegetation valid EVI range is 0 to 1 and all EVI values smaller than 0 set to NA
fold_c <- paste(path_c,dlist_c[i],sep="/") # the respective DOY folder of the Data_Pixel_Reliability folder
fls_c <- dir(fold_c,pattern=".tif") # all files that are available in the respective Pixel Reliability DOY folder
flsp_c <-paste(fold_c,fls_c,sep="/") # all files that are available in the respective Pixel Reliability DOY folder with complete path name
cloudstack <- stack(flsp_c) # creates a layer stack of all files within the Pixel Relaibility DOY folder
cloudresize<- crop(cloudstack,border) # resizes the Pixel Reliability layer stack to the rectangular extent of the border shapefile
cloudmask<-mask(cloudresize,border) # masks the Pixel Reliability layer stack using the border shapefile
cloudmask[cloudmask==(3)]<-NA # Pixel Reliability rank 3 pixels (cloudy) set to NA
cloudmask[cloudmask==(2)]<-NA # Pixel Reliability rank 2 pixels (snow&ice) set to NA
cloudmask[cloudmask==(0)]<-1 # Pixel Reliability rank 0 pixels (good quality) set to 1
cloudmask[cloudmask>(3)]<-NA # as valid Pixel Reliability range is -1 to 3, all Pixel Reliability values >3 set to NA
# (as -1 rank pixels show value 255)
evi_c=evi*cloudmask # multiplying the EVI layer stack by the Pixel Reliability layer stack
# to get one single layer stack with applied cloud mask
# extracting mean and standard deviation for each pixel
evimean <- stackApply (evi_c, rep (1, nlayers (evi)), mean, na.rm=T) #calculating the mean for the layer stack for each individual pixel
evisd <- stackApply (evi_c, rep (1, nlayers (evi)), sd, na.rm=T) #calculating the standard deviation for the layer stack for each individual pixel
# if na.rm is TRUE NA values are ignored, otherwise an NA value in any of the arguments will cause a value of NA to be returned,
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/Extremes.html
SVI_all <- ((evi_c-evimean)/evisd) #calculating SVI
# define breaks and color scheme
stdev <- cellStats(SVI_all, stat='sd') # calculating the standard deviation for each layer of the raster data,
# returning as many values as number of input layers
stdev2 <- 2*stdev # 2 standard deviations, returning as many values as number of input layers
stdev15 <- 1.5*stdev # 1.5 standard deviations, returning as many values as number of input layers
breaksSVI <- c(-4, -stdev2, -stdev15, -stdev, stdev, stdev15, stdev2, 4) # definition of the breaks of the resulting maps based on statistics
# (standard deviation), returning as many values as number of input layers
# times 6 (i.e. 6 values for each of the 15 layers) plus 2
# (the -4 in the beginning and 4 in the end)
l <- nlayers(SVI_all)
x <- seq(2, 6*l+1, l)
colorsSVI<-c('#4C0E0D', '#E72223', '#F19069', '#F9F6C6', '#64B14B', '#04984A', '#00320E') # definition of the color scheme of the resulting maps
fold_jpg <- paste(path_jpg) # the respective folder where you want to store the resulting .jpg-images.
fold_tif <- paste(path_tif) # the respective folder where you want to store the resulting .tif-files.
for (k in 1:nlayers(SVI_all)) { # start of the inner for-loop
year <- substr(fls[k],5,8) # extracting the fifth to eigths letter of the filename, which is the year (cf. data preparation above)
doy <- substr(fls[k],1,3) # extracting the first to third letter of the filename, which is the DOY (cf. data preparation above)
# writeRaster(evi[[k]], filename=paste(fold,"/",doy,"_",year,sep=""), format="ENVI", datatype='FLT4S', overwrite=TRUE)
# in case you would like to have Envi files (Attention: note the datatype)
jpeg(filename=paste(fold_jpg,"/",doy,"_",year,".jpg",sep=""), quality = 100)
# writes the jpg maps and names the files autmatically according to the pattern DOY_YYYY
plot(SVI_all[[k]], # zlim=c(0,100),
breaks=c(-4, breaksSVI[x], 4), # sets the breaks as defined above
col=colorsSVI, # sets the colors as defined above
main=paste("SVI"," (EVI)"," sample ",doy," ",year,sep="")) # automizes the title of the plot.
# ToDo: Adjust the file naming according to the data you are processing!
# E.g. if you base your SVI on EVI data, write (EVI)
dev.off()
writeRaster(SVI_all[[k]], filename=paste(fold_tif,"/",doy,"_",year,".tif",sep=""), format="GTiff", overwrite=TRUE)
# writes the geotiff and automizes the file naming according to the pattern DOY_YYYY
} # end of the inner for-loop
setTxtProgressBar (pb, i)
} # end of the outer for-loop
>>>>>>> 1903c5fe845011d166b80bef2c1c966c5cefc03c
```