## ----NDVIdownload, message=FALSE----------------------------------------------
library(terra)
datadir <- "modis"
files <- agrodata::data_ibli("marsabit_modis_ndvi", datadir)
head(basename(files))


## ----aoi----------------------------------------------------------------------
aoi <- agrodata::data_ibli("marsabit")
sinprj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs "
pols <- project(aoi, sinprj)


## ----NDVIs--------------------------------------------------------------------
# Define time period
startyear <- 2008
endyear   <- 2015

# Define a season (LRLD in this case)
sSeason <- "LRLD"
startSeason <- "03-01"
endSeason   <- "09-30"


## ----selectModFs, message=FALSE-----------------------------------------------
library(luna)
library(meteor)

selectModisFiles <- function(files, startdate, enddate) {
  base_names <- basename(files)
  dates <- luna::dateFromYearDoy(substr(base_names, 10, 16)) 
  i <- (dates >= as.Date(startdate)) & (dates <= as.Date(enddate) )
  files[i]
}


## ----NDVIs_Plots, fig.width=10, fig.height=10---------------------------------
season <- selectModisFiles(files, paste0(startyear, "-", startSeason), paste0(startyear, "-", endSeason))
sndvi <- rast(season)
plot(sndvi[[1:4]])


## ----NDVIs_Plots2-------------------------------------------------------------
ndvi_val <- terra::extract(sndvi, pols, fun=mean, na.rm=TRUE)
# the first region
plot(unlist(ndvi_val[1, -1]), ylab= "NDVI", pch=20, cex=2)


## ----UseselectModFs-----------------------------------------------------------
workdir <- "work"
dir.create(workdir, recursive=TRUE, showWarnings=FALSE)

for(year in startyear:endyear) {
  season <- selectModisFiles(files, paste0(year, "-", startSeason), paste0(year, "-", endSeason))
  
  sndvi <- rast(season)
  ndvimean <- mean(sndvi, na.rm = TRUE)
  
  # Shorten layer names 
  names(sndvi) <- paste0("ndvi_", luna::dateFromYearDoy(substr(names(sndvi), 10, 16)))
  filename=file.path(workdir, paste0( 'MOD09A1.h21v08.',year,"_",sSeason,'_ndvit.tif'))
  writeRaster(ndvimean, filename = filename , overwrite=TRUE)
}


## ----MorePlots----------------------------------------------------------------
par(mar = c(2.2, 2.2, 1.2, 0.5)) #c(bottom, left, top, right)
stitle <- paste(year, " NDVI", sep="")
plot(ndvimean, main = stitle)


## ----NDVISpatial--------------------------------------------------------------
season <- "LRLD"
files <- list.files(workdir, pattern=paste0(season,"_ndvit.tif$"), full.names=TRUE)
ndvi <- rast(files)

# Compute spatial aggregates by sub-location
output <- extract(ndvi, pols, fun=mean, na.rm=TRUE)
# remove the first column with IDs
output <- output[,-1]
colnames(output) <- substr(basename(files),16,19)
output[1:3, 1:6]


## ----comb---------------------------------------------------------------------
# Create a data.frame
res <- data.frame(SLNAME=pols$SUB_LOCATION, IBLI_Zone=pols$IBLI_UNIT, output, stringsAsFactors = FALSE, check.names=FALSE)
saveRDS(res, file.path(workdir, paste0(season,"_ndvi_st_mat.rds")))


## ----NDVIls_Plots-------------------------------------------------------------
output <- as.matrix(output)
plot(output[1,], ylab="NDVI", xlab="Year", type="l", ylim=c(0,.6))
for (i in 2:nrow(output)) {
  lines(output[i,])
}


## ----zscore-------------------------------------------------------------------
zscore <- function(y){
  (y - mean(y, na.rm=TRUE) ) / (sd(y, na.rm=TRUE))
}


## ----zScores------------------------------------------------------------------
scores <- apply(res[,-c(1:2)], 1, zscore)
scores <- t(scores)
scores[1:5, 1:6]

ndvi_st <- data.frame(SLNAME=pols$SUB_LOCATION, IBLI_Zone=pols$IBLI_UNIT, scores, stringsAsFactors = FALSE, check.names=FALSE)
saveRDS(ndvi_st, file.path(workdir, paste0(sSeason,"_zndvi_st_mat.rds")))



## ----PlotzScores, fig.width=12, fig.height=12---------------------------------
years <- 2008:2015

par(mfrow=c(2,2))
for(i in 1:4){
  
	plot(years, scores[i,], ylab="zNDVI", xlab="Year", main=pols$SUB_LOCATION[i])
	lines(years, scores[i,])

# Line showing where NDVi is 1 sd below the mean
	abline(h=-1, col="red",lty=2)
}