Table of Contents

Part. A

Step 0 of 6. Data aquisition and preparation

Step 1 of 6. Login or Register

Step 2 of 6. Entering download temporal and spatial coverage in appEEARS

Step 3 of 6. Downloading and preparing country or region boundary shapefile

Step 4 of 6. Appliying the country or region boundary

Step 5 of 6. Selecting MODIS bands of interest

Step 6 of 6. Conclusion of request and download of data & data.txt file

Part. B

Step 0 of 6. R, R-Studio and R-code preparations

i. R and R-studio download and installation

ii. Final preparations - R Code

Step 1 of 6. Downloading Files if required

Step 2 of 6. Starting EVI/VCI block Data preparation & EVI masking

Step 3 of 6. Applying VCI calculation to each chunk that was previously created in step 2

Step 4 of 6. VCI chunk merging and output

Step 5 of 6. The seasonality/phenology block

Step 6 of 6. Weighted linear combination of VCI output and EVI Phenology output

---

Part. A

---

---

Step 0 of 6 Data aquisition and preparation

---

Data access - MODIS

For this guide, the appEEARS platform will be used to acquire free, ready-to-use data-samples. appEEARS is an easy to use web GUI for downloading MODIS data specifically for your study area without the need for any pre-processing. appEEARS is provided by USGS

---

Step 1 of 6 Login or Register

---

Go to https://lpdaacsvc.cr.usgs.gov/appeears/ and click “Sign In” in the upper right-hand corner, “LOG IN” with your Earthdata credentials or if you do not have an Earthdata account, click on “REGISTER” and follow the instructions on the page to create one

On the top menu, select “Extract”, then “Area Sample”.

---

Step 2 of 6 Entering download temporal and spatial coverage in appEEARS

---

Enter a request name, the start and end date for your data. Keep in mind that the Vegetation Condition Index (VCI) is more reliable for more years of data (higher temporal extent). Because the final output is based on seasonal variability, please make sure that your start date is first of January and your end date marks 31st of December. If your study-area is located south of the equator, make sure the start date is first of July and the end date is 30th of June. Also, please note that the first full month of MODIS data is March 2000. If you want to select your study area manually, you can do so using the map and its tools to the right (right map box). If you want to use specific country-borders, please see step 3) as follows.

---

Step 3 of 6 Downloading and preparing country or region boundary shapefile

---

Important: If your country border consists of multiple polygons, please use the appEEARS map tools to draw a single polygon around your study area instead.
Another option is Via http://www.gadm.org/country, you can download shapefiles for every country of the world for use within appEEARS. Visit the website by opening a new tab or window and select a country from the upper drop-down menu. From the lower drop-down menu, select “Shapefile” and click on “OK”.

On the next page, you will see a preview of the shapefile. Click on “download” (below the preview) and save the zip-file. Now, open the zip-file and delete all files except the ones ending in “adm0”. All others refer to different administrative levels. In case you are interested in a certain province or district boundary, you can extract this with some GIS-operation.

---

Step 4 of 6 Appliying the country or region boundary

---

Finally, drag-and-drop the prepared zip-file into the designated area on the appEEARS website

Your area will now appear in the map. If you are notified that “The shape you selected was too complex to display and was automatically reduced [.]” - just proceed with the following steps. It will have no negative effect on data aquisition.

---

Step 5 of 6 Selecting MODIS bands of interest

---

The data required is the 16-day Enhanced Vegetation Index (EVI) composites. In the “Search for a product” box, first, enter “MOD13Q1” and select the “MOD13Q1.006” 250m vegetation indices, then select 250m 16-days EVI and 250 16-days pixel-reliability (optional, select if quality masking is wanted)

Lastly, select “GeoTiff” Output option, click the box below (“Search for a projection”) and select “Geographic” projection. If you have done everything correctly, it should look like this below:

You can now click “submit” your request using the button on the bottom of the page. Depending on the size of your study area and the selected time-period, the data preparation may take some hours. You can check the progress using the “Explore” tab on the top of the page. You will also get an email notification as soon as your request is finished

---

Step 6 of 6 Conclusion of request and download of data & data.txt file

---

Once your request is finished, click on the “Explore” tab, then on the download-button of your request.

Scroll down, select “All” and save the “Save Download list”. Remember the location of your download-list. Alternatively if you wish to download and save all data for analysis at a later time select “ALL” and then download to a suitable folder of your choice.

---

Part. B

---

---

Agricultural drought impact classification using weighted linear combination of EVI phenology and VCI (See “in detail” page of the Recommended Practice for more information)

---

---

Step 0 of 6 R, R-Studio and R-code preparations

---

i. R and R-studio download and installation

This practice uses R Script. Make sure to install R if it is not installed on your machine yet (visit https://cran.r-project.org/). The R-Studio interface is also used. We recommend installing it as well (visit https://www.rstudio.com/). Before running the R-Code, you will need to download the greenbrown R-package which is used to calculate phenology (more info: http://greenbrown.r-forge.r-project.org/). Visit http://greenbrown.r-forge.r-project.org/tar/, download “greenbrown_2.4.3.tar.gz” (do not extract the package) and remember the location you saved it to.

ii. Final preparations - R Code

Download the R-Code for this practice here and open it in R-Studio. Before running, you will need to add information to the following line below as seen in the image and script.

Folders for input and output of the analysis of the agricultural drought impact

In line 6 from image above, enter a path where you have enough free space for all products and temporary files. Depending on the study-area, a lot of gigabytes may be needed.

Note: To write lines in r as comments use “#” before starting the sentence, therefore if you find any comment after “#” it is not ment is not mentto be evaluated and it will not be evaluated by R

In case you already downloaded the data from appears manually beforehand, enter the path to your data, then follow the instructions in “Step 1 of 6. Downloading Files” by commenting out the lines like seen below using # as seen below

#filesl <- scan(downloadList, what='list', #sep='\n')
#for (d in 1:length(filesl)){
#  if (grepl("_VI_Quality_",filesl[d]) == F){
#    curl_download(url=filesl[d], #destfile=paste0(dataPath,"/",basename(filesl[d])#),quiet = T, handle = new_handle())
#  }
#  print(paste0('Downloading source data (Step 1 #of 6): ',round(d / length(filesl) * 100, #digits=2),'%'))
#}

In line 10, enter the exact path (including file-ending) to your download-list that you received from appEEARS. In line 13, enter the exact path (including file-ending) to the greenbrown-package you downloaded.

Enter the path to the a folder where you have enough free space to store your MODIS data and the resulting products Please do not use backslashes(\) except you are working with a MAC or Linux

dataPath <- "D:/Ukraine_workshop/Kiev_out"

Enter the path to your appEEARS download-list (including “.txt”!). If you downloaded the data from appEEARS beforehand, enter the path to the data here (all data needs to be in one folder)

downloadList <- "C:/Users/mduguru/Desktop/Ukraine_workshop/data_text/kiev-download-list.txt"

Enter the location of the greenbrown-package you downloaded according to the step-by-step guide (including “.tar.gz”!)

dir_GB <- "D:/Ukraine_workshop/greenbrown/greenbrown_2.4.3.tar.gz"

Hemisphere of your study-area: 1 = north hemisphere ; 0 = south hemisphere. We are considering the Southern Hemisphere for this workflow

hemisphere <- 0

Quality masking: replace the ‘0’ with a ‘1’ to apply MODIS Pixel Reliability data. If you are not sure whether you should enable this option, it is recommended that you find out more about Pixel Reliability. Check the “in Detail” page for more information.

cloudmask <- 1

---

Installing all needed packages for the determinition of the Hazard classes

####Installing all needed packages
if ("curl" %in% rownames(installed.packages()) == FALSE) {
  install.packages("curl")
}
if ("raster" %in% rownames(installed.packages()) == FALSE) {
  install.packages("raster")
}
if ("strucchange" %in% rownames(installed.packages()) == FALSE) {
  install.packages("strucchange")
}
if ("Kendall" %in% rownames(installed.packages()) == FALSE) {
  install.packages("Kendall")
}
if ("ncdf4" %in% rownames(installed.packages()) == FALSE) {
  install.packages("ncdf4")
}
if ("bfast" %in% rownames(installed.packages()) == FALSE) {
  install.packages("bfast")
}
if ("phenopix" %in% rownames(installed.packages()) == FALSE) {
  install.packages("phenopix")
}
if ("fields" %in% rownames(installed.packages()) == FALSE) {
  install.packages("fields")
}
if ("doParallel" %in% rownames(installed.packages()) == FALSE) {
  install.packages("doParallel")
}
if ("snow" %in% rownames(installed.packages()) == FALSE) {
  install.packages("snow")
}
if ("quantreg" %in% rownames(installed.packages()) == FALSE) {
  install.packages("quantreg")
}
if ("parallel" %in% rownames(installed.packages()) == FALSE) {
  install.packages("parallel")
}
if ("zoo" %in% rownames(installed.packages()) == FALSE) {
  install.packages("zoo")
}
if ("sandwich" %in% rownames(installed.packages()) == FALSE) {
  install.packages("sandwich")
}
#Installing greenbrown for seasonality
if ("greenbrown" %in% rownames(installed.packages()) == FALSE) {
  install.packages(dir_GB, repos = NULL, type="source")
}
if ("rgdal" %in% rownames(installed.packages()) == FALSE) {
  install.packages("rgdal")

Load packages using the the function library when required

library(phenopix)
library(Kendall)
library(ncdf4)
library(fields)
library(doParallel)
library(strucchange)
library(bfast)
library(snow)
library(quantreg)
library(parallel)
library(zoo)
library(sandwich)
library(raster)
library(curl)
library(greenbrown)
library(rgdal)

---

Step 1 of 6 Downloading Files if required

---

Download data EVI and data Pixel Reliability which was previously ordered from the Application for Extracting and Exploring Analysis Ready Samples (appEEARS). The text file previously saved contains all granules ordered (except ‘VI quality control’ data as they are not needed)

If you downloaded your data from appEEARS beforehand, comment out the following eight lines below using #

filesl <- scan(downloadList, what='list', sep='\n')
for (d in 1:length(filesl)){
  if (grepl("_VI_Quality_",filesl[d]) == F){
    curl_download(url=filesl[d], destfile=paste0(dataPath,"/",basename(filesl[d])),quiet = T, handle = new_handle())
  }
  print(paste0('Downloading source data (Step 1 of 6): ',round(d / length(filesl) * 100, digits=2),'%'))
}

Creating a temp-directory for intermediate calculations and setting memory limit to improve the computing capacity of r

dir.create(paste0(dataPath,'/temp/'))
rasterOptions(tmpdir=paste0(dataPath,'/temp/'))
rasterOptions(tolerance=1)
memory.limit(80000)

Listing up all downloaded TIF-files in the data folder

rasterData <- list.files(path=dataPath, pattern='.tif$', recursive=F, ignore.case=T, full.names=T)
rasterFiles <- basename(rasterData)

Chunk-size for processing (in Pixel)
You can reduce these numbers if you run into RAM-Problems early on: Try reducing the numbers to 1000. If you still run into RAM-problems, reduce to 500 or 250.
Only reduce the numbers if you have memory issues (“Error: cannot allocate vector of size […]”)

chwidth <- 1000
chheight <- 1000

---

Step 2 of 6 Starting EVI/VCI block Data preparation & EVI masking

---

Listing all EVI rasters and their corresponding pixel reliability data

EVIrasterData <- rasterData[grepl('EVI',rasterData)] 
EVIqc <- rasterData[grepl('pixel_reliability',rasterData)] 

Automatically adjusting chunk size for small scenes if needed

initial <- raster(EVIrasterData[1])
if (as.numeric(ncol(initial)) <= chwidth || as.numeric(nrow(initial)) <= chheight){
  chwidth <- ceiling(ncol(initial)/2)
}

Parsing all DOY and Years from the filenames

DOYs <- unique(substr(basename(EVIrasterData),38,40))
YEARs <- unique(substr(basename(EVIrasterData),34,37))

VCI: chunkwise calculation

Creating output folder

dir.create(paste0(dataPath,'/VCI'))
dir.create(paste0(dataPath,'/EVI'))

Determining chunk-shapefile

Loading example image from the downloaded data

exRST <- raster(EVIrasterData[1])

exPR<-raster(EVIqc[1])

Determining chunks

h   <- ceiling(ncol(exRST)/ceiling(ncol(exRST) / chwidth))
v   <- ceiling(nrow(exRST)/ceiling(nrow(exRST) / chheight)) 

Creating shapefile for each chunk

split     <- aggregate(exRST,fact=c(h,v))
split[]   <- 1:ncell(split)
splitshp  <- rasterToPolygons(split)
names(splitshp) <- 'shapes'
notiles <- ncell(split)

Filter values for quality masking clouded pixels: “0” and “1” in the pixel reliability bands are accepted as sufficient quality

Rank Key	Summary QA
-1	Fill/No Data
0	Good data
1	Marginal Data
2	Snow/Ice
3	Cloudy

For this analysis we will be considering good data and marginal data

goodValues <- c(0,1)

Masking clouds/snow; Splitting data into chunks

for (d in 1:length(DOYs)){
  #Filtering all relevant data for this DOY
  vrasterData <- EVIrasterData[grepl(paste0(DOYs[d],'_aid'),EVIrasterData)]
  #..and their corresponding pixel reliability data
  vQC <- EVIqc[grepl(paste0(DOYs[d],'_aid'),EVIqc)]
  #Reading years of available data
  vYear <- substr(basename(vrasterData),34,37)
  for (y in 1:length(vYear)){
    viPRE <- raster(vrasterData[y])
    #Applying quality mask to each image (if masking was activated)
    if (cloudmask == 1){
      qc <- raster(vQC[y])
      viPRE <- overlay(viPRE, qc, fun = function(x, y) {
        x[!(y %in% goodValues)] <- NA
        return(x)
      })
    }

Splitting (masked) data into Chunks

for(i in 1:ncell(split)){
      ex          <- extent(splitshp[splitshp$shapes==i,])
      exx         <- crop(viPRE,ex)
      writeRaster(exx,filename=paste0(dataPath,'/temp/',DOYs[d],'_',vYear[y],'_EVICHUNK',formatC(i, width=3, flag='0')),format='GTiff', overwrite=TRUE)
}
}
  #(Progress report)
  print(paste0('Data preparation (VCI) & masking (Step 2 of 6): ',round(d / length(DOYs)  * 100, digits=2),'%'))
}

---

Step 3 of 6 Applying VCI calculation to each chunk that was previously created in step 2

---

#List all chunks
EVIchunks <- list.files(path=paste0(dataPath,'/temp/'), pattern='_EVICHUNK', recursive=F, ignore.case=T, full.names=T)
for (d in 1:length(DOYs)){
  for (t in 1:notiles){
    #Filtering relevant chunks for this specific date
    sEVIchunks <- EVIchunks[grepl(paste0(DOYs[d],'_'),EVIchunks)] 
    sEVIchunks <- sEVIchunks[grepl(paste0('_EVICHUNK',formatC(t, width=3, flag='0')),sEVIchunks)] 
    #Listing years of data available for each DOY
    vvYear <- substr(basename(sEVIchunks),5,8)
    
    if (length(sEVIchunks) > 0){
      sT <- stack(sEVIchunks)
      
      #Removing filler-values from EVI data
      sT[sT<(-2999)] <- NA
      
      #VCI formula application
      vimax <- stackApply(sT , rep(1, nlayers (sT)), max, na.rm=T)
      vimin <- stackApply(sT , rep(1, nlayers (sT)), min, na.rm=T)
      z <- vimax - vimin
      VCI <- ((sT -vimin)/z)*100
      
      #Writing VCI-chunks for each available year
      for (y in 1:length(vvYear)){
        writeRaster(VCI[[y]],filename=paste0(dataPath,'/temp/',DOYs[d],'_',vvYear[y],'_VCICHUNK',formatC(t, width=3, flag='0')),format='GTiff', overwrite=TRUE)
      }
    }
  }
  #(Progress report)
  print(paste0('VCI processing (Step 3 of 6): ',round(d / length(DOYs)  * 100, digits=2),'%'))
}

---

Step 4 of 6 VCI chunk merging and output

---

#Listing all created chunks
VCIchunks <- list.files(path=paste0(dataPath,'/temp/'), pattern='_VCICHUNK', recursive=F, ignore.case=T, full.names=T)
#Looping through each available DOY_YEAR combination
for (y in 1:length(YEARs)){
  for (d in 1:length(DOYs)){
    #Listing relevant chunks (EVI, VCI) for this specific date
    sVCIchunks <- VCIchunks[grepl(paste0(DOYs[d],'_',YEARs[y],'_'),VCIchunks)]
    sEVIchunks <- EVIchunks[grepl(paste0(DOYs[d],'_',YEARs[y],'_'),EVIchunks)] 
    #Creating raster-list of the EVI chunks
    if (length(sEVIchunks) > 0){
      if (length(sEVIchunks) > 1){
        sMos <- list()
        for (o in 1:length(sEVIchunks)){
          sMos <- append(sMos, raster(sEVIchunks[o]))
        }
        #Mosaicking the EVI chunks
        sMos$fun = mean
        Mos <- do.call(mosaic, sMos)
        Mos <- Mos*0.0001
        #Output of mosaicked, scaled EVI image
        writeRaster(Mos,filename=paste0(dataPath,'/EVI/',DOYs[d],'_',YEARs[y],'_EVI'),format='GTiff', overwrite=TRUE)
      }else{
        #Fail-safe in case only one chunk is available
        exp <- raster(sEVIchunks[1])
        exp <- exp*0.0001
        writeRaster(exp,filename=paste0(dataPath,'/EVI/',DOYs[d],'_',YEARs[y],'_EVI'),format='GTiff', overwrite=TRUE)
      }
    }
    #Creating raster-list of the VCI chunks
    if (length(sVCIchunks) > 0){
      if (length(sVCIchunks) > 1){
        sMos <- list()
        for (o in 1:length(sVCIchunks)){
          sMos <- append(sMos, raster(sVCIchunks[o]))
        }
        #Mosaicking the VCI chunks
        sMos$fun = mean
        Mos <- do.call(mosaic, sMos)
        
        #Output: final VCI mosaic
        writeRaster(Mos,filename=paste0(dataPath,'/VCI/',DOYs[d],'_',YEARs[y],'_VCI'),format='GTiff', overwrite=TRUE)
      }else{
        #Fail-safe in case only one chunk is available
        writeRaster(raster(sVCIchunks[1]),filename=paste0(dataPath,'/VCI/',DOYs[d],'_',YEARs[y],'_VCI'),format='GTiff', overwrite=TRUE)
      }
    }
    print(paste0('EVI output & VCI output (Step 4 of 6): ',round((d+(length(DOYs)*(y-1)))/(length(DOYs)*length(YEARs))*100, digits=2),'%')) 
  }
}

Removing temp files to free space

tmp <- list.files(path=paste0(dataPath,'/temp/'), recursive=F, ignore.case=T, full.names=T)
file.remove(tmp)

---

Step 5 of 6 The seasonality/phenology block

---

Creating output folder for the seasons

dir.create(paste0(dataPath,'/season/'))

Seasonality for each year

Listing EVI data which is the base for seasonality calculation

allEVI <- list.files(paste0(dataPath,'/EVI/'), pattern=".tif$", recursive=F, full.names = T)

Creating a stack of all EVI data

sEVI <- stack(allEVI)

Passing a starting year to greenbrown for proper output filenames

stDate <- c(as.numeric(min(YEARs)),01,01)

Start multicore processing

• The seasonality-calculation will take advantage of all CPU cores. The PC might lag during the process and may not be usable for other tasks in a productive way
• If you want to keep using your PC during the process, reduce the number of used cores by one. To do this, replace max(as.numeric(Sys.getenv(‘NUMBER_OF_PROCESSORS’)),1) with max(as.numeric(Sys.getenv(‘NUMBER_OF_PROCESSORS’))-1,1) in the following line.

no_cores <- max(as.numeric(Sys.getenv('NUMBER_OF_PROCESSORS')),1)
print(paste0(Sys.time()," Seasonality calculation (this will take a while; Step 5 of 6)"))
beginCluster(no_cores)

Seasonality calculation

#Check the greenbrown-package documentation at http://greenbrown.r-forge.r-project.org/ and the "in detail" page for more information on the calculation
phenmap <- clusterR(sEVI,fun=PhenologyRaster,args=list(start = stDate, freq = 23, tsgf=NULL, approach="Deriv", check.seasonality=NULL))

Ending multi-core processing

endCluster()

Output Phenology

n <- NamesPhenologyRaster(phenmap, start=as.numeric(YEARs[1]))
for (i in 1:nlayers(phenmap)){
  writeRaster(phenmap[[i]],filename=paste0(dataPath,"/season/",n[i],".tif"),format="GTiff",overwrite=T)
}

---

Step 6 of 6 Weighted linear combination of VCI output and EVI Phenology output

---

List all finished VCI images

allVCI <- list.files(paste0(dataPath,'/VCI/'), pattern=".tif$", recursive=F, full.names = T)

Create an output-folder for the final output

dir.create(paste0(dataPath,'/weighted'))

Looping through each year(/season)

for (sy in as.numeric(min(YEARs)):as.numeric(max(YEARs))){
  
  #Fail-safe: check if seasonality data was calculated for the year
  if (file.exists(paste0(dataPath,"/season/","SOS.",sy,".tif")) && !(sy == max(YEARs) && hemisphere != 1)){
  
  #Check for hemisphere, load VCI data for the year
  if (hemisphere == 1){
  lyVCI <- allVCI[grepl(sy,allVCI)]
  dscr <- sy
  }else{
  #Load correct data if southern hemisphere was selected (seasonality shift)
  lyVCI <- c(allVCI[grepl(sy,allVCI)],allVCI[grepl(sy+1,allVCI)])
  if (sy == min(YEARs)){
  lyVCI <- lyVCI[1:23]
   }else{
  lyVCI <- lyVCI[13:35]
  }
  #Adjusting output filenames if southern hemisphere was selected
  dscr <- paste0(sy,"_",sy+1)
  }

Stacking all VCI for the year

  yVCI <- stack(lyVCI)

Loading seasonality rasters for the year

  SoS <- raster(paste0(dataPath,"/season/","SOS.",sy,".tif"))
  Peak <- raster(paste0(dataPath,"/season/","POP.",sy,".tif"))
  EoS <- raster(paste0(dataPath,"/season/","EOS.",sy,".tif"))

Creating VCI means for the three “blocks” of seasonal timing

Start of season (SOS) variables

BlockA <- stack()
  #Iterating through all 23 images of the season
  for (i in 1:23){
    #Read start-of-season dates
    s <- SoS
    #Read dates up to 16 days before peak of season
    e <- (Peak - 1)
    
    #Check if current iteration lies within the seasonal dates of each pixel
    s[s > i] <- NA
    e[e < i] <- NA
    
    #Creating a binary mask and apply to the VCI image
    m <- s*e
    v <- yVCI[[i]]*(m/m)

    #Add layer to stack
    BlockA <- stack(BlockA, v)
    
  }

Create mean of this time-block Block-A Start Of Season

BlockA <- calc(BlockA, fun=mean, na.rm=T)

Peak of season variables

See comments above

BlockB <- stack()
  for (i in 1:23){
    s <- Peak
    e <- Peak
    
    s[s > i] <- NA
    e[e < i] <- NA
    
    m <- s*e
    v <- yVCI[[i]]*(m/m)
            
    BlockB <- stack(BlockB, v)
  }

Create mean of this time-block Block-B Peak of season

BlockB <- calc(BlockB, fun=mean, na.rm=T)

End of season (EOS) variables

See comments above

BlockC <- stack()
  for (i in 1:23){
    s <- (Peak + 1)
    e <- EoS
    
    s[s > i] <- NA
    e[e < i] <- NA

    m <- s*e
    v <- yVCI[[i]]*(m/m)
       
    BlockC <- stack(BlockC, v)
  }

Create mean of this time-block Block-c EOS

 BlockC <- calc(BlockC, fun=mean, na.rm=T)

Weighting: 5-7-1 start-peak-end of season mean

Result: weighted VCI by season

WeightedStack<-stack(BlockA,BlockA,BlockA,BlockA,BlockA,
                     BlockB,BlockB,BlockB,BlockB,BlockB,BlockB,BlockB,
                     BlockC)

Create the mean of the season weighthed block stacks

WeightedVCI <- calc(WeightedStack, fun=mean,na.rm=T)

Raster mean classification

Output: five classes

cWeightedVCI <- reclassify(WeightedVCI, c(-1, 10, 4,  10, 20, 3,  20, 30, 2, 30,40, 1,40,101,0))  
writeRaster(cWeightedVCI,filename=paste0(dataPath,"/weighted/",dscr,"_droughtFiveCl.tif"),format="GTiff",overwrite=T)

Output: three classes required for the next steps of the analysis as recomended by the Technical guidance for monitoring and reporting on progress in achieving the global targets of the Sendai Framework for Disaster Risk Reduction

#(Progress report)
  dWeightedVCI <- reclassify(WeightedVCI, c(-1, 10, 2,  10, 40, 1,  40, 101, 0))  
  writeRaster(dWeightedVCI,filename=paste0(dataPath,"/weighted/",dscr,"_droughtThreeCl.tif"),format="GTiff",overwrite=T)
  
#(Progress report)
print(paste0('Weighting indices (Step 6 of 6): ',round((sy-as.numeric(min(YEARs))) (as.numeric(max(YEARs))-as.numeric(min(YEARs)))  * 100, digits=2),'%'))
}
}

Deleting all temp files

tmp <- list.files(path=paste0(dataPath,'/temp/'), recursive=F, ignore.case=T, full.names=T)
file.remove(tmp)

Finished

#Finished
print(paste0(Sys.time(),' - Finished. You will find the data in the subfolders \'VCI\' and \'weighted\' (classified drought indicators) at ', dataPath,"."))

Once you entered all necessary information, select the whole code (CTRL + a) and click on “run” (or: CTRL + Enter). Depending on your study-area, calculations will take several hours. The script will keep you updated with the current steps being performed in the console. You will also see a message in the console once it is completely finished. Sometimes execution might seem unresponsive. As long as you can see the little stop sign above the console, you can be sure the script is still running (see image below)