MODIS data

Introduction

In this chapter we discuss processing of MODIS satellite data for rice field detection and yield estimation. We use MODIS data between 1 January 2002 and 31 December 2018. We use a number of MODIS data products, including the Enhanced Vegetation Index (EVI), Gross Primary Productivity (GPP) and Evapotranspiration (ET). We use some of these products to detect areas planted under rice and other to predict yield. For all products, the general data processing workflow is described in more detail here for NDVI.

Search and download MODIS tiles

To illustrate the workflow, we show how to download and pre-process the data. But here we only do that for a few days, to limit the download and processing needed. In the next chapters, we use use pre-processed data for the entire time period.

We use the tempdir function to get the path to the temporary files directory on your computer, and we will store that in the datadir variable (you can of course change that to another location if you wish). We will organize the files in the following directory structure:

../agrins
│
└───raw
│   │   MOD*.hdf
│
└───vi
│   │   MOD*_evi.tif
│   │   MOD*_ndfi.tif
│
└───composite
│   │   filter_evi_year*_stack.tif
│   │   filter_ndfi_year*_stack.tif
│   │   year*_phenology.tif
│
└───index
    │   index*.csv

We first configure directories for storing the MODIS data.

datadir <- file.path(dirname(tempdir()), "agrins")
moddir <- file.path(datadir, "raw")
vidir <- file.path(datadir, "vi")
dir.create(moddir, recursive = TRUE, showWarnings = FALSE)
dir.create(vidir, recursive = TRUE, showWarnings = FALSE)

Next we download MODIS tiles that we need. More details can be found here.

To download the data you need to provide the username and password for your (free) EOSDIS account. If you do not have an account, you can sign up here. Here we use passwords that are stored in a file that I read in below (sorry, we cannot show you the values).

up <- readRDS("../pwds.rds")
up <- up[up$service == "EOSDIS", ]

library(terra)
library(luna)
library(agrodata)
# study area
aoi <- data_rice("zones")
# download MODIS data; smaller time window to limit the download
start <- "2010-01-01"
end <- "2010-01-07"
fmod <- getModis(product="MOD09A1", start, end, aoi, download=TRUE, path=moddir, username=up$user, password=up$pwd)
fmod
## [1] "c:/temp/agrins/raw/MOD09A1.A2009361.h21v09.061.2021149145412.hdf"
## [2] "c:/temp/agrins/raw/MOD09A1.A2010001.h21v09.061.2021150094231.hdf"
## [3] "c:/temp/agrins/raw/MOD09A1.A2009361.h21v09.006.2015198070303.hdf"
## [4] "c:/temp/agrins/raw/MOD09A1.A2010001.h21v09.006.2015198080654.hdf"
length(fmod)
## [1] 4
fmod[1:2]
## [1] "c:/temp/agrins/raw/MOD09A1.A2009361.h21v09.061.2021149145412.hdf"
## [2] "c:/temp/agrins/raw/MOD09A1.A2010001.h21v09.061.2021150094231.hdf"

Pre-processing

The study area is covered by a single MODIS tile, but the actual region is even smaller! For cropping the MODIS tiles with study area boundary, we reproject the boundary to MODIS sinusoidal projection system.

prj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
pols <- project(aoi, prj)

Create quality mask

Now we remove pixels that are affected by clouds and cloud-shadow. We follow the method described here. We specify a matrix in which each row has start and end of the quality assessment (QA) bits considered, and the values to be rejected (in the image).

from <- c(1,3,11,12)
to   <- c(2,6,11,14)
reject <- c("10,11", "1100,1101,1110,1111", "1", "000,110,111")
qmat <- cbind(from, to, reject)

Apply mask and compute vegetation indices

After masking out the bad pixels, we compute number of vegetation indices for each tiles. The pre-processing steps remain the same. First, we produce a quality mask, mask out bad quality pixels from MODIS images, clip to area of study extents and compute vegetation indices for each scene. However, in this example we show how to perform the operation for 2 tiles that can also be applied over hundreds of others.

for(i in 1:length(fmod)){
  # first we specify output filename
  evifile <- gsub(".hdf$", "_evi.tif", basename(fmod[i]))
  ndfifile <- gsub(".hdf$", "_ndfi.tif", basename(fmod[i]))
  # add output file directory
  evifile <- file.path(vidir, evifile)
  ndfifile <- file.path(vidir, ndfifile)
  if (file.exists(ndfifile) & file.exists(evifile)) next
  # read hdf file
  r <- rast(fmod[i])
  # get the bands we need. For this product
  # red=1; nir=2; blue=3; swir1=6; quality=12
  r <- r[[c(1:3, 7, 12)]]
  names(r) <- c("red", "nir", "blue", "swir2", "quality")
  # to make the process faster crop image using AOI before other calculations
  r <- crop(r, pols)
  # Generate quality mask
  quality_mask <- modis_mask(r$quality, 16, qmat)
  # remove ("mask") low quality pixels
  r <- mask(r, quality_mask)
  # Ensure all values are between 0 and 1
  r <- clamp(r, 0, 1)
  ## Vegetation indices
  # Enhanced vegetation index (EVI)
  evi <- 2.5*((r$nir - r$red)/(r$nir + (6*r$red)-(7.5*r$blue)+1))
  # Normalized difference flood index, instead of LSWI
  ndfi <- (r$red - r$swir2)/(r$nir + r$swir2)
  # save results
  writeRaster(evi, filename = evifile, overwrite = TRUE)
  writeRaster(ndfi, filename = ndfifile, overwrite = TRUE)
}

We calculated only two vegetation indices that are important for studying rice. However a number of other indices can be computed from the available bands. Sensitivity of vegetation indices to crop growth depends on number of factors including crop type, variety, management, location and season. Hatfield and Prueger, 2010 compared multiple vegetation indices to quantity characteristics of maize, soybean, wheat, and rapeseed under difference management practices using eight years of data.

Additional MODIS datasets

In the previous sections, we showed how to download and process MOD09A1 product that can be used for computing vegetation indices. However, a number of other MODIS products such as evapotranspiration (ET), Fraction of Photosynthetically Active Radiation (FPAR), Leaf Area Index (LAI), and gross primary productivity (GPP) can be of interest for this exercise. We can follow the same steps described above to download and process these additional datasests. Note that the layer names and quality mask handling will be different for each products. More information about these datasets can be found with productInfo('MOD17A2H'), productInfo('MOD15A2H'), and productInfo('MOD16A2'). Here is an example of processing GPP products (go here for understanding the quality flags used. Search for “Description: QC (quality control) flags for Psn_500M biophysical variable”

se <- matrix(c(1,1,4,5), ncol=2, byrow=TRUE)
reject <- list("1",c("01","10"))
# function to crop, mask and save gpp for the study area
getGPP <- function(i, fmod, pols, vidir) {
  # first we specify output filename
  gppfile <- gsub(".hdf$", "_gpp.tif", basename(fmod[i]))
  gppfile <- file.path(vidir, gppfile)
  # read hdf file
  r <- rast(fmod[i])
  # Crop the raster using the AOI
  r <- crop(r, pols)
  # Generate quality mask
  quality_mask <- modis_mask(r[[3]], 8, se, reject)
  # band index for this product
  r <- r[[1]]
  names(r) <- c("gpp")
  # Mask out bad quality pixels
  r <- mask(r, quality_mask)
  # save output
  writeRaster(r, filename = gppfile, overwrite = TRUE)
}
# download
getModis(product="MOD17A2H", start="2010-01-01", end ="2010-01-31",
         aoi, download=TRUE, path=moddir)
# list GPP hdf files in the folder
ffmod <- list.files(moddir, pattern = "*.hdf", full.names = TRUE, recursive=TRUE)
fmod <- grep("MOD17A2H", ffmod, value = TRUE)
# process multiple files
lapply(1:length(fmod), getGPP, fmod, pols, vidir)

This workflow can be easily modified and applied to other MODIS products as well.