## ----aoi, message=FALSE-------------------------------------------------------
aoi <- agrodata::data_ibli("marsabit") 
aoi


## ----ken, message=FALSE-------------------------------------------------------
library(terra)
ken <- geodata::gadm("KEN", level=1, path=".")
plot(ken, border="gray")
plot(ken[ken$NAME_1 == "Marsabit", ], col="beige", border="gray", add=TRUE)
lines(aoi, col="red")


## ----downloaddir--------------------------------------------------------------
datadir <- "modis"
dir.create(datadir, FALSE)
# datadir


## ----dates--------------------------------------------------------------------
start <- "2010-01-01" 
end <- "2010-01-07"


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


## ----moddown, message=FALSE---------------------------------------------------
library(luna)
fmod <- luna::getModis(product="MOD09A1", start, end, aoi, download=TRUE, 
			path=datadir, username=up$user, password=up$pwd)
fmod


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


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


## ----Preprocess---------------------------------------------------------------
for(i in seq_along(fmod)) {
  # Load a MODIS file
  r <- rast(fmod[i])

  # Generate quality mask (layer 12)
  quality_mask <- modis_mask(r[[12]], 16, qmat)

  #Select only red and NIR
  r <- r[[1:2]]
  names(r) <- c("red", "NIR")

  # remove low quality pixels
  r <- mask(r, quality_mask)

  # remove areas outside the AOI rectangle
  r <- crop(r, pols)

  # ensure that all values are between 0 and 1 
  r <- clamp(r, 0, 1) 

  # Compute NDVI
  ndvi <- (r$NIR - r$red)/(r$NIR + r$red)

  filename <- gsub(".hdf$", "_ndvi.tif", fmod[i])

  writeRaster(ndvi, filename=filename, overwrite=TRUE)
}