Vjylku

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?

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)





# 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)