Installing relevant packages
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
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
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
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
path_c <- "C:/User/MOD13Q1/Data/Data_Pixel_Reliability"
dlist_c <- dir(path_c,pattern="DOY")
Showing an example of the downloaded EVI data
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.
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:
“=======”"
pb <- txtProgressBar (min=0, max=length(dlist), style=1)
setTxtProgressBar (pb, 0)
Main code for processing the EVI data, masking out clouds, calculating the VCI and writing the results
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 max and min value for each pixel
evimax <- stackApply (evi_c, rep (1, nlayers (evi)), max, na.rm=T) #calculating the max value for the layer stack for each individual pixel
evimin <- stackApply (evi_c, rep (1, nlayers (evi)), min, na.rm=T) #calculating the min value 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
VCI_all <- ((evi_c-evimin)/(evimax-evimin))*100 #calculating VCI
# define breaks and color scheme
my_palette <- colorRampPalette(c('#8B0000', '#FF4500', '#FFFF00', '#9ACD32', '#008000'))
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(VCI_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(VCI_all[[k]],
zlim=c(0,100),
col=my_palette(101), # sets the colors as defined above
main=paste("VCI"," (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 VCI on EVI data, write (EVI)
dev.off()
writeRaster(VCI_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