Step by Step: Recommended Practice drought monitoring including a cloud mask (R) – Example Central America

Overview drought
In detail drought

This page provides the step-by-step instructions to apply a multi-temporal analysis of the MODIS-based Vegetation Condition Index (VCI) for large areas using R. Depending on the extent of the study area and the hardware capacities of the computer, the user might experiece difficulties applying this R script. We therefore provided an additional R script to be able the compute the VCI and monitor the vegetation conditions for large areas (>200.000km2) which can be accessed here.

For this Recommended Practice for drought monitoring we use the Enhanced Vegetation Index (EVI) based on MODIS data at 250m spatial resolution as input for generating the Vegetation Condition Index (VCI). The only MODIS product providing EVI data is MOD13Q1 (on Terra satellite) and MYD13Q1 (on Aqua Satellite). MODIS data is for example made available on the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) by USGS.

If the Normalized Difference Vegetation Index (NDVI) is preferred for calculating the VCI, the EVI can easily be replaced with the NDVI as this product is also available for download at AppEEARS. Given this case, the caption of the resulting jpg maps should be adjusted. However, in certain circumstances, the VCI based on the NDVI does not perform as well as based on the EVI. See the "In Detail" page for more information.

Data and software required for this Recommended Practice are:

  • Shapefile of area of interest (https://gadm.org/country)
  • EVI and Pixel Reliability data provided by MOD13Q1 (see next section)
  • R and R-Studio

 If not yet installed, the R software and R-Studio can be downloaded here. To follow a step by step instruction on how to download and install R and R-Studio, please click here.

 

Content

Part A: Data preparation

Part B: Processing

Part C: Results

Step 1 Download MODIS data as Geotiffs from AppEEARS

Access the AppEEARS website and create a user account (free) in case you do not already have one (https://lpdaacsvc.cr.usgs.gov/appeears/). Click on Extract and Area Sample and start a new request.

  1. Enter a request name
  2. Define your region of interest by specifying a ESRI shapefile (zip) or drawing a polygon . Country shapefiles can be downloaded at https://gadm.org/country. Note: The boundaries and names of the shapefiles do not imply official endorsement or acceptance by the United Nations.
  3. Define the time period of your data: January 2000 to today.
  4. Select the product: MOD13Q1 and the layers of interest: EVI and pixel_reliability
  5. Define the output as geotiff with geographic Projection
  6. Click on Submit

Click on Explore to check the status of the request. Click on the Download symbol once the status says Done. Select all ordered datasets and start downloading the data.

 

After the download store the files (.tif) in two folders: one folder for the Data_EVI and one for the Data_Pixel_Reliability. It is recommended to use an external hard drive to store your input and output data if possible.

Step 2 Renaming and Structuring of the data

For processing this Recommended Practice in R, it is required to first rename and second sort the data according to their DOY. Download the folder structure here and copy it both in the folders Data_EVI and Data_Pixel_Reliability.

  • Rename the EVI and pixel reliability data following the pattern: DOY_YYYY_[Original Name]. It is recommended to use Total Commander to rename multiple files (Download Link: http://www.ghisler.com/index.htm). Renaming the files is important because this pattern is used to automatize the filenames and the titles of the resulting maps in the R script.
  • Sort the EVI and Pixel Reliability data according to the DOY in the respective folders.

Create following folders in addition to Data_EVI and Data_Pixel_Reliability:

  • "shape" (store the shapefile with the country border here including .shp, .dbf, .prj, and .shx files)
  • "VCI_Maps_jpg"
  • "VCI_Maps_tif"

Step 3 Download R Code

Download the R Code here (right click and choose "save as"). The R Code is tailored to MOD13Q1 tif data.

Step 4 Adjust the script

Open R Studio and the saved R script (File -> Open File). The script needs various adjustments to fit your data. Note hereby: R uses slashes ("/") for file names and directories and not backslashes ("\").

  • In line 23: border <- shapefile("C:/User/MOD13Q1/Data/shape/gadm36_GTM_0.shp") --- replace the data path by the data path and/or name of the shapefile that you have stored on your computer
  • In line 26: path <- "C:/User/MOD13Q1/Data/Data_EVI" --- enter link to the folder where you have stored the MODIS EVI data
  • In line 29: path_c <- C:/User/MOD13Q1/Data/Data_Pixel_Reliability" --- enter link to the folder where you have stored the MODIS Pixel Reliability data
  • In line 32: path_jpg <- "C:/User/MOD13Q1/Data/VCI_Maps_Guatemala_jpg" --- enter the link to the folder where you want to store the resulting .jpg-images
  • In line 33: path_tif <- "C:/User/MOD13Q1/Data/VCI_Maps_Guatemala_jpg" --- enter the link to the folder where you want to store the resulting .tif-files
  • In line 102: plot(VCI_all[[k]], zlim=c(0,100), col=my_palette(101), main=paste("VCI"," (EVI)"," Guatemala ",doy," ",year,sep="")) --- Adjust the country and file naming according to the data you are processing (e.g. if you base your VCI on NDVI, change to (NDVI))

Step 5 Run the R code

Install the raster, rgdal and sp packages by highlighting the respective text in the code and clicking on "Run" or typing Ctrl+Enter. 

Now, run the remaining code by highlighting it and typing CRT+Enter. Depending on the size of your data, it will take quite a while (up to several hours) for the script to run - so be patient. After a while, a progress bar will be displayed in the console. Be patient, it might take a while for the progress bar to appear. It looks like this:

======

In case you get an error message "Error: cannot allocate vector of size ... Mb" it might be that you do not have enough RAM. It is recommended to use a 64-bit operating system, and not a 32-bit system.

Another problem could be that there might not be enough space on your harddisk for the temporary files. You can either delete the temporary files manually or change the temp directory to a drive with more storage capacity by typing: rasterOptions(tmpdir="D:/temp/R_raster_temp"), where "D:/temp/R_raster_temp" of course has to be replaced by the path of your choice.

Note: Delete the temporary files of your computer and set the memory limit within R (memory.limit()) according to your free working memory for improving your computer's capacity to process the data.

If you need to stop the skript, you can do so by clicking on the red Stop sign on the top right of the Console.

Access the full Code

You can find the full code on https://rpubs.com

Part C Results

The results will be stored in the respective folders as jpg (VCI_Maps_jpg) and as tif file (VCI_Maps_tif) within your main directory. The maps are named "DOY_YYYY".

The title of the jpeg maps is following the pattern "DOY VCI (EVI) YYYY", where DOY and YYYY are automatically adjusted according to the input data and (NDVI) can be adjusted manually in the script (see above). See below an example of the resulting jpeg map including cloud pixels and a map resulting from the old procedure without cloud mask:

White color inside the country shape indicates cloud covered pixels or water surfaces.