Step by Step: Recommended Practice: Drought monitoring using the Standard Vegetation Index (SVI) - R without Cloud Mask

Overview Drought SVI
In Detail Drought SVI

Note: The original version of the Recommended Practice on Drought Monitoring using the SVI did not include a cloud mask for the scenes analysed. Additionally, as this Recommended Practice was developed before the launch of AppEEARS, the MODIS data was downloaded via the USGS Earth Explorer. The data could could be downloaded as HDF format which then required a conversion to the geotiff format. This conversion was done through the Modis Reprojection Tool (MRT). Since this conversion tool may also be useful in other situtation, the original version of the Recommended Practice is made available as well. 

For this recommended practice  the Enhanced Vegetation Index (EVI) based on MODIS data at 250m spatial resolution is used as input. The only MODIS product providing EVI data is MOD13Q1 (on Terra satellite) and MYD13Q1 (on Aqua satellite). Here we will work with the EVI band of MOD13Q1. 

Content

Part A: Data preparation (hdf conversion, folder structure)

A+B+C: Data preparation (layer stacking rescaling), SVI and Visualization

Step 0.1: Download MODIS 13Q1 data

The MODIS 13Q1 data set can be accessed via different data platforms like Reverb, GloVIS, DataPool. We have downloaded via the USGS Earthexplorer. To download the data:

  1. Register at earthexplorer and login (it is free of charge).
  2. Enter your search criteria, i.e. the area and time period of interest.
  3. Choose the data set. You will find MOD13Q1 under NASA LPDAAC Collections - MODIS Vegetation Indices - MODIS MOD13Q1.
  4. Click "Results".
  5. Add all files to Bulk Download.
  6. Once you have finished adding all files to Bulk Download, click "View item Basket". 
  7. In the Item Basket choose "HDF Format" as your product. 
  8. Click "Proceed to Checkout" and "Submit Order".
  9. You will then receive your order confirmation via email including the link to download the Bulk Download Application (BDA). If you do not already have this installed, install the BDA on your computer.
  10. Open the BDA. Login with your Earthexplorer username and password. Choose your order, select your preferred download folder, and start the download.

The MODIS13Q1 data has a certain way to name the single files. The name contains relevant information like acquisition date or acquisition tile.

Example for hdf files:

MOD13Q1.A2000049.h09v07.005.2006269133152.hdf

MOD13Q1 - Product Short Name

.A2000049 - Julian Date of Acquisition (A-YYYYDDD) --> 49th day of the year 2000

.h09v07 - Tile Identifier (horizontalXXverticalYY) --> see image below.

.005 - Collection Version

. 2006269133152 - Julian Date of Production (YYYYDDDHHMMSS)

.hdf - Data Format (HDF-EOS)

Example for tif files:

MOD13Q1.A2000049.250m_16_days_EVI.tif

MOD13Q1 - Product Short Name

.A2000049 - Julian Date of Acquisition (A-YYYYDDD) --> 49th day of the year 2000

.250m_ - spatial resolution

.16_days_EVI – temporal resolution at 16 days Maximum value composite

.tif - Data Format (GeoTiff)

Step 0.2: Download country shapefiles

Country shapefiles can be downloaded at http://www.gadm.org/country. Note: The boundaries and names of the shapefiles do not imply official endorsement or acceptance by the United Nations.

Step 1: Reprojection of hdf into tif using MRT

If you have downloaded MODIS data in hdf format, you need to first convert them into tif files before using them in R software. These are the steps:

  1. Download and install Modis Reprojection Tool (MRT) (Download: https://lpdaac.usgs.gov/tools/modis_reprojection_tool)
  2. Open the MRT
  3. First, you have to create a Parameter File, which contains file names, types of input and output data files. Follow the steps below to create the parameter file.
  • Select any one of the Input File to set parameters (File > Open Input File), which will be used later on for the batch processing.
  • Under "selected bands" select the bands you need and move the remaining bands to "available bands" by highlighting them and clicking on the arrow; in this case we are interested in NDVI and EVI.

4. After inserting the parameters (cf. images above), save your Parameter File with Save Parameter File.

Note: Do not click on "run" nor on "convert format". The creation of parameter files for all input MODIS files as well as the conversion to tif will happen in the Command Prompt (cf. step 5).

5. Automated Batch Processing

In the following you will convert all the MODIS files (in hdf format) at once into tif format (batch processing). In summary these are the three steps for the batch processing:

  • Gather all input files into one exclusive directory (containing nothing but MRT input)
  • You will then create parameter files (.prm) for each input file in the input data directory with MRTBatch.jar. It also writes out a batch script that is later used to execute the processing jobs (MRTBatch.bat) (cf. below how to).
  • Finally, you wil convert the hdf files into tiff files (cf. below how to).

Click on Start > Search Programs and Files > Command Prompt

Using the Command Prompt, you then should navigate to the "MRT bin" directory, where the "MRTBatch.jar" file lays:

  • cd folder - This command will move you to the folder that you specify. The folder must be in the directory you are currently in. For example: If you are currently at C:\Users\username\ and you enter cd desktop you will be taken to C:\Users\username\Desktop\
  • dir - This command will list all of the folders and files in the directory you are currently at (this step is optional).

Type the following command to initiate parameter file and batch script generation (place all fields in one line; i.e., do not break this statement into multiple lines).

Note: the italic words have to be replaced by your actual directory or file name (cf. explanation below).

java –jar MRTBatch.jar –d input_directory –p parameter_directory\input_parameter_file –o parameter_directory

where

input_file_directory is the directory in which all input files are placed (e.g.,

C:\Project\InputFiles);

parameter_directory is the directory where the parameter file - created using the graphical user interface (GUI) of the MODIS Tool - was saved (-p) and to which the parameter files built by this application will be written (-o); Note that the –o parameter_directory is optional. Users can omit it from the command line, which will result in a prm directory containing the new parameter files being built in the input_file_directory instead.

input_parameter_file is the parameter file that was saved from the graphical user interface (GUI) of the MODIS Tool (e.g., mrtauto_test.prm).

Hit the enter key on your keyboard. This will result in the creation of parameter files for each MODIS hdf file.

Note: If an error occurs as in the picture below, possibly your filename contains invalid characters (spaces, underscore, hyphens...) or your file name (including the path) is too long. In that case, change the file name and/or path and re-type the command with the new filenames.

6. Convert all hdf files to tif files:

Make sure you are still in the MRT bin directory and type behind ...\MRT\bin>

mrtbatch.bat (in this example: C:/Users/dall/MRT/bin/mrtbatch.bat

Hit the enter key on your keyboard. The results of the processing job will echo out in log form on the command prompt display, and when completed, files should be ready in the output directory specified in the parameter_directory or in the new prm directory in the input_file_directory (cf. step 5).

Step 2: Database and file naming 

Create folders and sort files manually and rename files with Total Commander

Total Commander is a useful tool to rename multiple files: 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.

  1. Create one folder, where you will store all data, e.g. C:/Data.
  2. For each day of the year (DOY) create one subfolder. The name of each folder must start with "DOY_", e.g. DOY_033. Of course, you do not have to use data of all the available days of the year. You can choose those which are most relevant in your region of interest.
  3. Rename the files by adding a prefix following the pattern DOY_YYYY_, e.g. 033_2001 or 001_2005, where:
    • DOY = day of year
    • YYYY = year
  4. Create another subfolder within the main data folder called "shape". Store the shapefile with the country border here.

Do these steps automatically by using a Python script

  • The script renames the Modis files automatically. An additional prefix is added, following the pattern DOY_YYYY_, e.g. 033_2001 or 001_2005, where

The DOY (Day of the year) folders are created automatically and the files are moved to the folders accordingly.

  • DOY= day of year
  • YYYY=year
  1. Download Python and install it on your computer. We recommend to download Anaconda 3, if you are new to Python (https://www.anaconda.com/download/). You can find more information about Python on our Knowledge Portal.
  2. Download the script and save it on the computer.
  3. Adapt the path to the folder containing the files (output files from the Modis Reprojection Tool), for example: “C:/Modis_data/prm”, in the script and execute it.

We provide you with two options for using the Python script:

  • You could download the Jupyter Notebook script (.ipynb) and run it in Jupyter Notebook.

     

    1. You launch Jupyter Notebook by clicking on the icon in the Windows start menu and it will open in a web browser.
    2. Then you can navigate to the saved script (Rename_CreateFolders_Sort_ModisFiles.ipynb) and open it.
    3. Now you need to change the path to the folder, where the Modis files are stored on your computer.
    4. You need to run each Jupyter Notebook cell by clicking on the run button or hitting Shift + Enter. More information can be found here.
    5. Download the Jupyter Notebook script (.ipynb) here.
  • Or you could download the raw Python script (.py) and run it in Spyder.
    1. You open Spyder by clicking on the icon in the Windows start menu.
    2. You open the script by clicking on File, Open and navigate to the saved Python script (Rename_CreateFolders_Sort_ModisFiles.py).
    3. Now you need to change the path to the folder, where the Modis files are stored on your computer.
    4. You need to execute each cell by clicking on the run button or hitting Shift + Enter.
    5. Download the raw Python script (.py) here.

 

Step 0.3: Download R Studio

If this is the first time for you to use R software, you will need to download and install the R software AND R Studio (which is an extremely useful user interface) at http://www.rstudio.com/products/rstudio/download/

Step 0.4: Download R Code

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

Step 3-8: Run the R Code

1. Open R studio.

2. Open the R script in R Studio: File - Open File.

3. Adjust the script.

Note: R uses slashes (/) for file names and not backslashes (\). R is case sensitive.

  • In line 23 border<-shapefile("F:/Data/ElSalvadorHonduras/shapefileElSalvador/SLV_adm0_ElSalvador.shp"), replace F:/Data/ElSalvadorHonduras/shapefileElSalvador/SLV_adm0_ElSalvador.shp by the link to the shapefile that you have stored on your computer.
  • Inline 26 path <- "F:/Data/ElSalvadorHonduras/EVI13Q1", replace F:/Data/ElSalvadorHonduras/EVI13Q1 by the link to folder where you have stored the MODIS data.
  • In line 89 plot(VCI_all[[k]],zlim=c(0,100), col=my_palette(101),main=paste(doy," VCI "," (NDVI) ",year,sep="")). Adjust the file naming according to the data you are processing. This is only necessary if you use other input data than EVI (e.g. NDVI).

3. Install the raster package and the rgdal package by highlighting the respective text in the code (cf. figure below) and typing CTR+Enter.

4. 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 hard disk 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.

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

The results will be stored in the respective DOY folders as jpg (quicklook) and as tif. They are named "DOY_YYYY". 

The title of the jpeg maps is following the pattern "DOY VCI (NDVI) YYYY", where DOY and YYYY are automatically adjusted according to the input data and (NDVI) can be adjusted manually in the script (cf. above). See below an example of the resulting jpeg map: