Generating Time Series MODIS-NDVI/EVI Monthly Data using R and MRT Tools

During the past weeks, I have been trying to process MODIS Terra and Aqua product (MOD13A3 and MYD13A) data to produce monthly NDVI and EVI data of Indonesia from 2000 – 2015 (Terra) and 2002 – 2015 (Aqua). MOD13A3 and MYD13A3 is a monthly vegetation indices data generated from Modis 16 days data in 1 km spatial resolution. Those products have 11 bands such as NDVI, EVI, VI Quality, Blue, Red, NIR, MIR, Vew Zenith Angle, Sun Zenith Angle, relative azimuth angle, and pixel reliability. The problem is that MODIS data is available in tiles so that to generate time series NDVI data for Indonesia, re-projection and mosaic process are needed to be done. To mosaic and resample MODIS, Modis Reprojection (MRT) Tools is available to do this tasks. But reprojection and mosaic monthly data of 16 and 13 years of MODIS Terra and Aqua will become a long and tedious tasks if they are not performed in the batch/automatic mode. Here, I will show how to simplify the process of reprojection and mosaicking of MODIS – MOD13A3 and MYD13A3 tiles.

First, we need :

1. Modis Reprojection Tools (MRT) downloadable from https://lpdaac.usgs.gov/tools/modis_reprojection_tool
2. Rstudio and R installation
3. 3rd party renamer software such PyRenamer (linux) or Bulk Utility Rename (windows)
4. R package of rts – https://cran.r-project.org/web/packages/rts/index.html
5. MODIS data of your area of interest, can be downloaded from https://lpdaac.usgs.gov/data_access/daac2disk

Steps :

1. Download MODIS data from the website and store the data in your local drive. The downloaded data will have .hdf file and the accompanying .xml file :
   

   “MOD13A3.A2015335.h27v07.006.2016007181738.hdf
   MOD13A3.A2015335.h27v07.006.2016007181738.hdf.xml”

   please note the naming convention, MOD13A3 is the type of product, A2015335 is the acquisition date in “Year” followed by “Julian Date”, and H27v07 is the tiles number.

2. Ok, next we need to group the tiles based on acquisition date because we want to have  monthly basis observation. I did this in Linux-Ubuntu environment from Terminal but I believe windows environment could do the same automation. the steps are as follows :
   a. I first removed the type of product (MOD13A3 or MYD13A3)  from the filenames using third party renamer such pyrenamer or bulk renamer, so that the name of the tiles would be :

     MOD13A3. A2015335.h27v07.006.2016007181738.hdf

   b. go to your folder where you stored your MODIS tiles in terminal
   c. run the following code :

   #list the hdf files 
   for src in *.hdf
   do 
   #get the tiles name
   name=${src%%-*} 
   #calculate the length of the tiles names
   namelength=$((${#name}-1)) 
   #get the first eight characters (acquisition date i.e A2015335)
   letter=${name:0:8} 
   #make new directory based on the acquisition date name 
   mkdir -p "$dir/$letter" 
   #copy the tiles with the same acquisition date to the folder
   mv -u "$src" "$dir/$letter" done

   after you run this code, you will get this folder structure

   ../MOD13A3/
   ../MOD13A3/A2015305
   ../MOD13A3/A2015335
   .etc

3. Once the folder and files is ready, please do the installation of MRT tools. MRT installation is easy and can be found elsewhere. The installation of MRT also need the installation of Java, so prepare it before you install MRT.
4. Once MRT is ready, copy all of your data that has been produced in the step 2 (including the folder) into the directory of ~/MRT/data
5. Open RStudio, and load the library of rts and set ~/MRT/data as the working directory
   

   library(rts)
   setwd(~/MRT/data)

6. List the directory and subdirectory within working directory
   

   dir <- list.dirs(getwd()) #please note that the main dir also get listed

7. run the following loop to project and mosaic the hdf files in your subdirectory. the following example is using loop function. But you may use parallel computing package such as foreach and DoParallel package to speed up the process. Also dont forget to change the MRTpath based on the location of the mrtmosaic and resample in your MRT local directory.
   

   #start loop started from the first subdir (2nd element of the list)
   for(i in 2:length(dir)){
   #go into subdir [i]
   setwd(dir[i])
   #list the .hdf files in subdir [i]
   a <- list.files(getwd(), pattern = “.hdf$”)
   #perform mosaic from the listed hdf files in the subdir
   b <- mosaicHDF(hdfNames = c(a), filename = paste0(basename(getwd()), “_mos”, “.hdf”), MRTpath = ‘ /home/sanjiwana/MRT/bin’)
   #list the output mosaic hdf file
   c <- list.files(getwd(), pattern = “_mos.hdf$”)
   #reproject the mosaic hdf files to geographic projection with only the first 7 bands (NDVI, EVI, VI Quality, Blue, Red, NIR and MIR) in the MOD13A3 were kept. the output format is .tif
   reprojectHDF(hdfName = c(c), filename = paste0(basename(dir[i]), “_MOD13A3_GEO”, “.tif”), MRTpath = “/home/sanjiwana/MRT/bin”, bands_subset = “1 1 1 1 1 1 1 0 0 0 0”, proj_type = “GEO”, proj_params = “0 0 0 0 0 0 0 0 0 0 0 0 0 0 0”, datum = “WGS84”, pixel_size = 0.00892857)}

8. done

   

   I hope this will help someone to process time series data of NDVI. Please bear in mind that I am not an expert on coding so that you might be able to develop more modest workflow compared to what I’ve done.

   cheers, -SAK