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


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


## ----downloadMODIS, message=FALSE---------------------------------------------
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

length(fmod)
fmod[1:2]


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


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


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


## ----gpp_preprocess, eval=FALSE-----------------------------------------------
## 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)