## ----load_rice_phenology, message=FALSE---------------------------------------
library(luna)

# set directories
datadir <- file.path(dirname(tempdir()), "agrins")
compdir <- file.path(datadir, "composite")
dir.create(compdir, recursive = TRUE, showWarnings = FALSE)

# rice phenology information
year <- 2012
ricepheno <- readRDS(paste0(compdir, "/ricepheno_",year,"_samples.rds"))

# we need the following three for this exercise
vi <- ricepheno[["evi_fit"]]
dates <- ricepheno[["dates"]]
rice <- ricepheno[["phenorice"]]


## ----estimate_indices1--------------------------------------------------------

# valid rice records
s <- which(rowSums(rice, na.rm=TRUE) > 0)

# let's consider one of them
p1 <- rice[s[1], ]

# peak
print(p1)
peak <- p1[2]
start <- p1[1]
mature <- p1[5]

# corresponding fitted time series
v1 <- unlist(vi[s[1],])

# remove background signal
v1 <- v1 - min(v1)

# plots
plot(v1, type = "l", col = "darkgreen", ylim = c(min(v1), max(v1)*1.25), 
      ylab = "EVI", xlab = "date", xaxt = "n", main = "EVI based index")

at <- seq(1, length(v1), 4) 
labs <- months(dates[at])
axis(1, at = at, labels = substr(labs, 1,3), las = 2)

# vegetative
polygon(x=c(start,start:peak,peak), 
        y=c(0, v1[start:peak], 0), 
        col="green")

# grain-filling
polygon(x=c(peak, peak:mature, mature), 
        y=c(0, v1[peak:mature], 0), col="yellow")

# add full growing season
polygon(x=c(start,start:mature, mature), 
        y=c(0, v1[start:mature], 0), 
        col="red", density=10)

# add legend
legend("topleft", legend=c("vegetative","grain","full"), fill=c("green", "yellow", "red"), density=c(NA, NA, 20), bty="n", border="black")


## ----auc1---------------------------------------------------------------------
# vegetative growth
i1 <- sum(v1[start:peak], na.rm = TRUE)
  
# reproductive
i2 <- sum(v1[peak:mature], na.rm = TRUE)
  
# full season
i3 <- sum(v1[start:mature], na.rm = TRUE)

ind <- c(i1, i2, i3)  
names(ind) <- paste0("idx_evi_", c("v","g","f"))


## ----auc_raster---------------------------------------------------------------
getAUC <- function(v, peak, start, mature){
  
  ind <- c(NA, NA, NA)
  peak <- as.vector(peak)
  start <- as.vector(start)
  mature <- as.vector(mature)
  
  # check both NA and NaN
  if (sum(is.na(v)) > 5) return(ind)
  if(is.na(sum(peak, start, mature))) return(ind)
  
  # vegetative growth
  i1 <- sum(v[start:peak], na.rm = TRUE)
  
  # reproductive
  i2 <- sum(v[peak:mature], na.rm = TRUE)
    
  # full season
  i3 <- sum(v[start:mature], na.rm = TRUE)

  ind <- c(i1,i2,i3)  
  return(ind)
}

# apply to single pixel
a1 <- getAUC(v1, peak = 37, start = 32, mature = 42)
names(a1) <- paste0("idx_evi_", c("v","g","f"))

# compare the value with the earlier index calculated based phenorice derived dates
ind
a1


## ----metrics_raster-----------------------------------------------------------
library(terra)
year <- 2012
fevi <- rast(file.path(compdir, paste0("/filter_evi_", year, ".tif")))
phenoraster <- rast(file.path(compdir, paste0("/phenorice_", year, ".tif")))

# for simplicity we assume average phenology dates for the region
ar <- app(fevi, getAUC, peak = 37, start = 30, mature = 43)
names(ar) <- paste0("idx_evi_", c("v","g","f"))

# but we could replace this with something like this
#ff <- c(s+3, p+3, m+3, ar)
#ar = app(ff, function(x) c(sum(x[x[1]:x[2]], na.rm=T),
#                            sum(x[x[2]:x[3]], na.rm=T),
#                            sum(x[x[1]:x[3]], na.rm=T)) )
			 

# mask metrics by rice growing pixels
maskfun <- function(x){
  x <- sum(x, na.rm = TRUE) 
  x[is.na(x)] <- 0
  x[x==0] <- NA
  x[!is.na(x)] <- 1
  return(x)
}
# apply rice mask
ricemask <- app(phenoraster, fun = maskfun)
arm <- mask(ar, ricemask)


## ----spatial_agg_metrics------------------------------------------------------
# get insurance zones
z <- agrodata::data_rice("zones")
plot(z, main = "zones")

# transform z to the crs of arm
z <- project(z, crs(arm))

# extract zone-level metrics and add the id
mz <- extract(arm, z, fun=mean, na.rm=TRUE)
# merge with zones
values(z) <- cbind(as.data.frame(z), mz)


## ----zone_metrics_plot--------------------------------------------------------

# remove zones with few pixels
n <- extract(ricemask, z, fun=sum, na.rm=TRUE)
n <- which(n[,2] > 10)
mzv <- z[n,]

# We only plot the first metrics ~ AUC under vegetative growth
par(mfrow = c(1,2))
plot(arm[["idx_evi_v"]], main = "idx_evi_v", col=topo.colors(25))
plot(z, add = T)
plot(mzv, "idx_evi_v", main = "zone metrics", col=topo.colors(25))