## ----download_data------------------------------------------------------------
datadir <- file.path(dirname(tempdir()), "agrins")
vidir <- file.path(datadir, "vi")
dir.create(vidir, recursive = TRUE, showWarnings = FALSE)
ff <- agrodata::data_rice("MODIS", vidir)


## ----evi_file_details---------------------------------------------------------
library(luna)

# list VI files that we downloaded
fevi  <- list.files(path=vidir, pattern= "^MOD09A1.*_evi.tif$", full.names=TRUE)
fndfi <- list.files(path=vidir, pattern= "^MOD09A1.*_ndfi.tif$", full.names=TRUE)

# clean date of collections
fevi  <- luna::modisDate(fevi)
fndfi <- luna::modisDate(fndfi)

# check what the function returns
dim(fevi)
head(fevi)


## ----evi_subset_files---------------------------------------------------------
library(terra)
year <- 2012
startDate <- as.Date(paste0(year-1, "-07-01"))
endDate <- as.Date(paste0(year, "-06-30"))

evifiles  <- fevi[fevi$date >= startDate & fevi$date <= endDate, ]
ndfifiles <- fndfi[fndfi$date >= startDate & fndfi$date <= endDate, ]

# read evi and ndfi  
evi  <- rast(evifiles$filename)
ndfi <- rast(ndfifiles$filename)

# set appropriate names
# (these names should have been set in the data)
names(evi) <- paste0("evi_", evifiles$date)
names(ndfi) <- paste0("ndfi_", ndfifiles$date)



## ----rarap, fig.width=8-------------------------------------------------------
k <- seq(1, nlyr(evi), 4)
k
plot(evi[[k]])


## ----evi_smooth---------------------------------------------------------------
v <- unlist(evi[200])
plot(v, pch = "*", col = "red")
s <- filterVI(v)
lines(s, col = "darkgreen")
points(s, pch="+", col = "darkgreen")
legend("topleft", 
       legend = c("EVI (raw)", "EVI(fitted line)", "EVI(fitted)"),
       col = c("red","darkgreen", "darkgreen"),
       pch = c("*", NA, "+"), lty = c(NA, 1, NA), 
       bty = "n", cex = 0.75)


## ----evi_smoothing------------------------------------------------------------
# output direcotry
compdir <- file.path(datadir, "composite")
dir.create(compdir, recursive = TRUE, showWarnings = FALSE)

# run model and save output
fn <- file.path(compdir, paste0("/filter_evi_", year, ".tif"))
fevi <- app(evi, filterVI, filename = fn, overwrite=TRUE, wopt = list(names = names(evi)))

fn <- file.path(compdir, paste0("/filter_ndfi_", year, ".tif"))
fndfi <- app(ndfi, filterVI, filename=fn, overwrite =TRUE, wopt = list(names = names(ndfi)))

# also save the files used
saveRDS(evifiles, paste0(compdir, "/files_", year, "_evi.rds"))


## ----evi_plot, fig.width=8----------------------------------------------------
plot(fevi[[k]])


## ----phenoplot1---------------------------------------------------------------
# create random locations to test the rice detection method
set.seed(4)
cells <- sample(ncell(evi), 200)
xy <- xyFromCell(evi, cells)
s <- vect(xy, crs=crs(evi))

# extract time series 
ts_evi <- extract(evi, s, drop = TRUE)[,-1]
ts_fevi <- extract(fevi, s, drop = TRUE)[,-1]
ts_fndfi <- extract(fndfi, s, drop = TRUE)[,-1]

# Let's plot one sample that we already know is rice
i <- 121
dates <- evifiles$date
plot(dates, ts_evi[i,], pch = 16, col = "red", cex = 0.8, ylim = c(0,0.7), 
        ylab = "vegetation indices", xlab = "dates",
        main = paste0("location_",i))

# smoothed values
lines(dates, ts_fevi[i,], col = "green")
lines(dates, ts_fndfi[i,], col = "blue")

legend("topleft", legend = c("EVI (raw)", "EVI", "NDFI"),
         pch = c(16, NA, NA), lty = c(NA, 1, 1),
         col = c("red", "green", "blue"), bty = "n", cex = 0.75)


## ----phenoplot_ask, eval = FALSE, include=FALSE-------------------------------
## par(ask=TRUE)
## for(i in 1:nrow(ts_evi)){
##   # raw EVI
##   plot(dates, ts_evi[i,], pch = 16, col = "red", cex = 0.8, ylim = c(0,0.7),
##        xlab = "vegetation indices", ylab = "dates",
##        main = paste0("location_",i))
##   # gap-filled and smoothed EVI
##   lines(dates, ts_fevi[i,], col = "green")
##   # gap-filled NDFI
##   lines(dates, ts_fndfi[i,], col = "blue")
##   # legend
##   legend("topleft", legend = c("raw_EVI", "fitted_EVI", "gap-filled_NDFI"),
##          pch = c(16, NA, NA), lty = c(NA, 1, 1),
##          col = c("red","green", "blue"), bty = "n", cex = 0.75)
## }
## 
## par(ask=FALSE)


## ----phenorice1---------------------------------------------------------------
library(phenorice)
p <- list()
p$evi_meanth <- 0.5 # threshold for annual mean EVI
p$evi_maxth <- 0.4  # threshold for maximum evi 
p$evi_minth <- 0.3  # threshold for minimum evi 
p$pos_start <- 20   # start of heading, after December
p$pos_end <- 40     # end of heading, not possible after April
p$vl1 <- 5          # the shortest vegetative growth length
p$vl2 <- 20         # the longest vegetative growth length
p$winfl <- 3        # period for flooding around min EVI
p$minndfi <- -0.1   # threshold for ndfi
p$windecr <- 15     # period after EVI maximum 
p$decr <- 0.4      # fraction decrease of EVI before and after EVI maximum for EVI change ratecheck
p$tl1 <- 8          # shortest total growing length
p$tl2 <- 30         # the longest total growing length
p$lst_th <- 15      # the minimum temperature for planting 


## ----phenorice2---------------------------------------------------------------
rice <- matrix(nrow=nrow(ts_fevi), ncol=5)
colnames(rice) <- c("start", "peak", "head", "flower", "mature")
ts_fevi <- as.matrix(ts_fevi)
ts_fndfi <- as.matrix(ts_fndfi)

for(i in 1:nrow(ts_fevi)){
  rice[i, ] <- phenorice(evi = ts_fevi[i,], ndfi = ts_fndfi[i,], p)
}

# output
head(rice)

#ricedates <- dates[rice]


## ----phenorice4---------------------------------------------------------------
# Only keep records that has been detected as rice
s <- rowSums(rice, na.rm=TRUE) > 0
peak <- rice[s,1]

# distribution of peak
pdates <- as.Date(dates[peak])
hist(pdates, main = "Distribution of peak dates", xlab = "dates", 
     breaks = "months", las = 2, cex.axis = 0.5) 

## save phenology and rice information for future use
ricepheno <- list(ts_evi, ts_fevi, ts_fndfi, dates, rice)
names(ricepheno) <- c("evi_raw", "evi_fit", "ndfi", "dates", "phenorice")
saveRDS(ricepheno, paste0(compdir, "/ricepheno_", year, "_samples.rds"))


## ----rice_map_raster, fig.width=10--------------------------------------------
# save result to a file
phenoraster <- file.path(compdir, paste0("/phenorice_", year, ".tif")) 
x <- sds(fevi, fndfi)
phenomatrix <- function(x, y, p) {
  sapply(1:nrow(x), \(i) phenorice(x[i,], y[i,], p=p)) |> t()
}

r <- lapp(x, phenomatrix, filename=phenoraster, overwrite = TRUE, p=p,
         wopt = list(names = paste0(c("start", "peak", "head", "flower", "mature"),"_",year)))
plot(r)