1

I'm running into some issues with 250m 16-day MOD13Q1 NDVI data using the MODIS package in R.

Looking at the following raster plot, there are certain days where the NDVI does not match the true surface characteristics. Specifically ndvi.2017.11.01 and ndvi.2018.01.17. Even ndvi.2018.03.06 and ndvi.2018.02.18 are questionable. Is it possible to programmatically exclude these layers? Is there something to gain by using the following sds layer: MODIS_Grid_16DAY_250m_500m_VI:250m 16 days VI Quality?

enter image description here

Below is the code to produce the above plot:

library(raster)
library(MODIS)
library(parallel)
library(pbapply)
library(magrittr)
library(dplyr)
library(stringr)

# Final Study Area Boundary
study_area <- readRDS(gzcon(url('http://web.pdx.edu/~porlando/study_area.RDS')))  
# Study CRS
epsg_26910 <- "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "

# Download Data
runGdal(product = "MOD13Q1", collection = "006"
           ,tileH = 9, tileV = 4
           ,begin = str_replace_all("2017-09-01", "-", ".")
           ,end = str_replace_all("2018-09-30", "-", ".")
           ,overwrite = TRUE
           ,job="temporalComposite"
           #,SDSstring="" # not sure about which SDS string to use...
)


# Process data
ndvi_path <- "./MODIS/MOD13Q1.006" # 250 m

ndvi_files <- list.files(path = ndvi_path, pattern = "\\.hdf$"
                     ,all.files = FALSE, full.names = TRUE
                     ,recursive = TRUE
                     ,ignore.case = FALSE)


processNDVI <- function(file_path, study_area = study_area, proj = epsg_26910) {

  cat("\n")

  # extract date from file path
  date <- stringr::str_extract(file_path, '[0-9][0-9][0-9][0-9].[0-9][0-9].[0-9][0-9]')
  date_clean <- stringr::str_replace_all(date, "\\.", "-")

  # coerce 16 day average to daily composite?
  dates <- seq.Date(from = as.Date(date_clean)
                    #,to = as.Date(date_clean) + lubridate::days(15)
                    ,to = as.Date(date_clean)
                    ,by = "1 day")

  sds <- get_subdatasets(file_path)

  ndvi <- sds[grepl("16DAY_250m_500m_VI:250m 16 days NDVI", sds)] %>% readGDAL %>% raster
  evi <- sds[grepl("16DAY_250m_500m_VI:250m 16 days EVI", sds)] %>% readGDAL %>% raster
  # VI Quality may be useful?
  vi_quality <- sds[grepl("16DAY_250m_500m_VI:250m 16 days VI Quality", sds)] %>% readGDAL %>% raster

  # combine NDVI and EVI into a single stack
  r <- stack(ndvi, evi)
  names(r) <- c(paste0("ndvi.", date), paste0("evi.", date))

  r <- raster::projectRaster(from = r, res = 250, method = 'ngb', crs = proj)
  study_area <- spTransform(study_area, CRSobj = proj)

  if(identical(crs(study_area), crs(r))) {

    m <- mask(r, study_area)
    cr <- crop(m, study_area)

    pblapply(dates, function(x) {
      pblapply(1:nlayers(cr), function(y) {
        layer_name <- names(cr[[y]])
        var <- gsub("\\..*$", "", layer_name)
        date <- gsub("-", "\\.", x)
        writeRaster(cr[[y]]
                    ,filename = paste0("./output/", var, "/", var, ".", date, ".tif")
                    ,format = "GTiff"
                    ,overwrite = TRUE)
      })
    })
}

pblapply(ndvi_files, function(x) {
  processNDVI(file_path = x
              ,study_area = study_area
              )
  }
  ,cl = detectCores()-1
    )

ndvi <- stack(list.files(path = "./output/ndvi/", pattern = ".tif$", full.names = T))
plot(ndvi)
  • 1
    The VIQ is a 16-bit binary layer and can be useful but a good first check is the Pixel Reliability layer; values on that layer of -1, 2, and 3 are definitely problematic and the associated VI values are probably not reliable. 0 and 1 are good and marginal data respectively, and when PR = 1 then you should check out the VIQ layer for further information. Also note that the timestamps you show are for the layer but are not reflective of the date of composite for each pixel. – JepsonNomad Mar 29 at 3:27
  • Thank you for the info on Pixel Reliability. I will try to work this into my data processing. Can you elaborate a little more on the timestamp discrepancy? What do you mean by the "date of composite for each pixel"? – philiporlando Mar 29 at 20:46
  • 1
    MOD13Q1 data is reported in 16-day intervals at the tile level but is collected every day. The maximum NDVI value is reported for the 16-day window, and the date of NDVI composite is reported on the DOY band. See user guide for more info. – JepsonNomad Mar 29 at 20:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.