#----------------------------------------------------------------------------------------------
# CODE R POUR UTILISER LE PACKAGE BFAST POUR L'ANALYSE D'UNE SERIE TEMPORELLE D'IMAGES NDVI
# MODIS TELECHARGEES DEPUIS LE SITE "MODIS/VIIRS Land Product Subsets"
#----------------------------------------------------------------------------------------------
# SOURCE PRINCIPALE DU TUTORIEL: https://verbe039.github.io/BFASTforAEO/#MODIS_based_time_series_analysis_using_BFAST
# Quelques modifications/améliorations ont été apportées dans ce code
#----------------------------------------------------------------------------------------------
# Tentez de comprendre comment fonctionne ce code R en l'exécutant morceau par morceau et 
# visualisez les résultats intermédiaires et finaux. Faites éventuellement varier certains 
# paramètres afin de comprendre leur impact sur les résultats de l'analyse. 
#----------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------
rm(list=ls()) # Remove Objects from a Specified Environment (nettoyage des tests précédents éventuels)

# SETTING WORKING DIRECTORY RELATIVE TO THE R CODE TEXT FILE, 1 LEVEL HIGHER, AS RCODE IS IN A SUBFOLDER OF THE MAIN EXCHANGEABLE FOLDER: in order to exchange easily RCODE and relative data with others: the RCODE and relative data should all be included in the same directory folder that can then be sent to a colleague and reused easily
# - EXCHANGEABLE FOLDER = "8_SERIE_TEMPORELLE_BFAST" 
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))  # SETTING DIRECTORY RELATIVE TO THE CURRENT R CODE TEXT FILE https://stackoverflow.com/questions/13672720/r-command-for-setting-working-directory-to-source-file-location-in-rstudio?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa
getwd() # to visually check the current directory

#----------------------------------------------------------------------------------------------
# INSTALLATION DES PACKAGES NECESSAIRES POUR L'EXERCICE grâce à la fonction "library()"
# The necessary add-on packages need to be installed within R by using the library() function. 
# Below we define a helper function that does installing and loading of the packages for us.
# Cette section n'est pas très importante à comprendre pour l'analyse des des données MODIS. Elle peut être simplement exécutée.
#----------------------------------------------------------------------------------------------
# pkgTest is a helper function to load packages and install packages only when they are not installed yet.
pkgTest <- function(x)
{
  if (x %in% rownames(installed.packages()) == FALSE) {
    install.packages(x, dependencies= TRUE)    
  }
  library(x, character.only = TRUE)
}
neededPackages <- c("strucchange","zoo", "bfast", "raster", "leaflet", "MODISTools")
for (package in neededPackages){pkgTest(package)}


#----------------------------------------------------------------------------------------------
# CREATION DE 2 FONCTIONS NECESSAIRES POUR LA SUITE (PAS DETAILLE AU TP)
# Loading extra functions to create: 
# - via "mt_to_raster" a raster stack from the MODIS subset downloaded with the mt_subset function 
# - via "timeser" a time series object in R ts
# Cette section n'est pas très importante à comprendre pour l'analyse des des données MODIS. Elle peut être simplement exécutée.
#----------------------------------------------------------------------------------------------
mt_to_raster <- function (df = subset) 
{
  dates <- unique(df$calendar_date)
  df$scale[df$scale == "Not Available"] <- 1
  r <- do.call("stack", lapply(dates, function(date) {
    m <- matrix(df$value[df$calendar_date == date] * as.numeric(df$scale[df$calendar_date == 
                                                                           date]), df$nrows[1], df$ncols[1], byrow = TRUE)
    return(raster::raster(m))
  }))
  bb <- MODISTools::mt_bbox(xllcorner = df$xllcorner[1], yllcorner = df$yllcorner[1], 
                            cellsize = df$cellsize[1], nrows = df$nrows[1], ncols = df$ncols[1])
  bb <- as(bb, "Spatial")
  raster::extent(r) <- raster::extent(bb)
  raster::projection(r) <- raster::projection(bb)
  names(r) <- as.character(dates)
  return(r)
}

timeser <- function(index, dt) {
  z <- zoo(index, dt)
  yr <- as.numeric(format(time(z), "%Y"))
  jul <- as.numeric(format(time(z), "%j"))
  delta <- min(unlist(tapply(jul, yr, diff))) # 16
  zz <- aggregate(z, yr + (jul - 1) / delta / 23)
  (tso <- as.ts(zz))
  return(tso)
}


#----------------------------------------------------------------------------------------------
# TELECHARGEMENT DES DONNEES MODIS AVEC LE PACKAGE MODISTools 
# ET EN PARTICULIER LA FONCTION "mt_subset"
# Downloading MODIS data using the MODISTools package
# First we download the MODIS data via the mt_subset function and immediately create raster brick from it:
#---------------------------------------------------------------------------------------------- 
# !!! ATTENTION, CHOISISSEZ OU AJOUTEZ LE "Site ID" DU "MODIS FIXED SITE SUBSETS" A ANALYSER, par exemple:
# Confer le « Site ID » de l'onglet « Subset Summary » du produit "MOD13Q1" du site choisi sur le site https://modis.ornl.gov/sites/ 
MODIS_Fixed_Sites_Subsets_ID<-"nl_gelderland_loobos" # aux Pays-Bas
MODIS_Fixed_Sites_Subsets_ID<-"ke_laikipia_mpala"    # au Kenya

# mt_subset: Download MODIS Land Products subsets
# (Confer le document "MODISTools.pdf" ou la rubrique d'aide de RSTUDIO pour une description des paramètres de cette fonction)

# Télécharger le couche "250m_16_days_NDVI" CETTE ETAPE PEUT PRENDRE QUELQUES SECONDES: VOIR LA BARRE DE PROGRESSION
VI <- mt_subset(product = "MOD13Q1",
                site_id = MODIS_Fixed_Sites_Subsets_ID,
                band = "250m_16_days_NDVI",
                start = "2000-01-01",
                end   = "2019-10-10",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = TRUE)

str(VI)         # Accès à la structure de l'objet ('data.frame':	492 228 obs. of  21 variables); 492 228 observations = 33 pixels(lignes) * 33 pixels(colonnes) * 452 images; 452 images = (23 images/an * 18 ans (2001-2018) + 20 images en 2000 + 18 images en 2019)
names(VI)       # Accès aux noms de colonnes de l'objet
dim(VI)         # Accès aux dimensions de l'objet (492 228 lignes et 21 colonnes)
VI[1:5,]        # Affichage des 2 premières lignes de l'objet
head(VI, n=2L)  # Affichage des 2 premières lignes de l'objet
unique(VI$calendar_date) # Affichage des 452 dates des images, sans doublons: il y a 452 dates dans la série temporelle
VI[453,]        # Affichage de la lignes 453 de l'objet qui correspond au 2ème pixel (colonne pixel = 2).


# Téléchager la couche "250m_16_days_pixel_reliability" CETTE ETAPE PEUT PRENDRE QUELQUES SECONDES: VOIR LA BARRE DE PROGRESSION
# Plus d'information ici table 4 notamment:  https://lpdaac.usgs.gov/documents/103/MOD13_User_Guide_V6.pdf
QA <- mt_subset(product = "MOD13Q1",
                site_id = MODIS_Fixed_Sites_Subsets_ID,
                band = "250m_16_days_pixel_reliability",
                start = "2000-01-01",
                end   = "2019-10-10",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = TRUE)

str(QA)
names(QA)
dim(QA)
head(QA,n=2L)


#----------------------------------------------------------------------------------------------
# Creating a raster brick and cleaning the MODIS data using the reliability layer
# More information about the pixel reliability can be found here (see table 4): 
# https://lpdaac.usgs.gov/documents/103/MOD13_User_Guide_V6.pdf
#----------------------------------------------------------------------------------------------
# CONVERSION DU FORMAT "DATFRAME" AU FORMAT "RASTERSTACK"
# convert df to raster
VI_r <- mt_to_raster(df = VI)
VI_r
QA_r <- mt_to_raster(df = QA)
QA_r

## clean the data
# create mask on pixel reliability flag set all values <0 or >1 NA
m <- QA_r
m[(QA_r < 0 | QA_r > 1)] <- NA # obtain good data

# mask all values from VI raster NA
VI_m <- mask(VI_r, m,maskvalue=NA, updatevalue=NA)

dim(VI_m) # 452 raster of NDVI (452 dates) (365jours/16jours = 22.8 ; 22.8 images par année * 19 années = 437 images ; 22.8*20=456)

# Plot the first image: 
# VISUALISEZ LES IMAGES DANS L'ONGLET "PLOTS" de RSTUDIO
plot(m,1) # plot mask

plot(VI_m,1) # plot masked NDVI raster
plot(VI_m,2)
plot(VI_m,3)
plot(VI_m,4)
plot(VI_m,450)

#----------------------------------------------------------------------------------------------
# OPTIONAL: You can extract data from the cleaned VI raster brick via the click function:
# A NE PAS FAIRE DANS LE CADRE DU TP
#----------------------------------------------------------------------------------------------
# ()
# click(VI_m, id=TRUE, xy=TRUE, cell=TRUE, n= 1)

#----------------------------------------------------------------------------------------------
# OPTIONAL: Creating a nice map with the leaflet package in R
#----------------------------------------------------------------------------------------------
library(leaflet)
r <- raster(VI_m,1)
pal <- colorNumeric(c("#ffffff", "#4dff88", "#004d1a"), values(r),
                    na.color = "transparent")

m <- leaflet() %>% addTiles() %>%
  addRasterImage(r, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(r),
            title = "NDVI")
m

#----------------------------------------------------------------------------------------------
# Apply BFASTmonitor
# Below we extract the data from the raster as a vector and create a time series using the timeser function:
#----------------------------------------------------------------------------------------------
## check VI data at a certain pixel e.g. 1 row, complete left hand site:
## the dimensions of the raster are: 33x33
## Confer the "Pixel Numbering Scheme" here for example https://modis.ornl.gov/cgi-bin/sites/site/?id=nl_gelderland_loobos&product=MOD13Q1
## Pour visualiser la position d'un numéro de pixel au sein d'une grille de 33*33
## Pixel du centre  = le numéro 545 

#-----------
# Quelques tests de consultation de VI_m (pour info)
dim(VI_m)   # 33 33 452 (452refers to number of dates of acquisition)
dim(VI_m[]) # 1089 (= 33*33) 452 (nombre de date)
str(VI_m)
VI_m[]
VI_m[1]     # 1er pixel pour les 452 images
VI_m[34]    # 34ème pixel pour les 452 images
VI_m@data
#-----------

#----------------------------------------------------------------------------------------------
# CHOIX D'UN PIXEL, CREATION D'UNE SERIE TEMPORELLE ET CREATION D'UN GRAPHIQUE
#----------------------------------------------------------------------------------------------
# CHOISISSEZ LE NUMERO DE PIXEL QUE VOUS VOULEZ ETUDIER
px <- 545 # N° 545 = Pixel du centre de la grille 33*33

# CREATION D'UNE SERIE TEMPORELLE POUR LE PIXEL SELECTIONNE
tspx                <- timeser(as.vector(VI_m[px]),as.Date(names(VI_m), "X%Y.%m.%d")) # convert selected pixel to a time series for the 452 dates
tspx
tspx_minus_12_dates <- timeser(as.vector(VI_m[px][1:440]),as.Date(names(VI_m)[1:440], "X%Y.%m.%d")) # convert selected pixel to a time series for the dates 1 to 440

# REALISATION D'UN GRAPHIQUE A PARTIR DE LA SERIE TEMPORELLE DU PIXEL SELCTIONNE
# VISUALISEZ LES IMAGES DANS L'ONGLET "PLOTS" de RSTUDIO
plot(tspx, main = 'NDVI') # NDVI time series cleaned using the "reliability informaiton"


#----------------------------------------------------------------------------------------------
# ANALYSE DE LA SERIE TEMPORELLE AVEC LA FONCTION BFAST
#----------------------------------------------------------------------------------------------
# --> IDENTIFICATION DES TENDANCES, SAISONS, CASSURES, BRUITS
# !! ## The data should be a regular ts() object without NA's --> confer BFAST package manuel, page 6

# REMPLACEMENT DES NA DE LA SERIE PAR INTERPOLATION AVEC "na.spline"
tspx              # Il y a des NA ! 
anyNA(tspx)       # Pour vérifier rapidement s'il y a au moins un NA
sum(is.na(tspx))  # Compte du nombre de Na

tspx_NO_NA <-timeser(na.spline(as.vector(VI_m[px]),na.rm = FALSE) ,as.Date(names(VI_m), "X%Y.%m.%d")) # # remplacement des NA par une valeur interpolée : confer l'aide de cette fonction & convert selected pixel to a time series
tspx_NO_NA[is.na(tspx)] # valeurs remplaçantes des NA

# APPLICATION DE LA FONCTION BFAST ET GRAPHIQUE DU RESULTAT
# VISUALISEZ LES IMAGES DANS L'ONGLET "PLOTS" de RSTUDIO
# CETTE ETAPE PEUT PRENDRE QUELQUES SECONDES
bfast_results<-bfast(tspx_NO_NA, h = 0.15, season ="harmonic", max.iter = 10, breaks = 10, hpc = "none", level = 0.05, type= "OLS-MOSUM")
plot(bfast_results, type = c("components", "all", "data", "seasonal","trend", "noise"))

#----------------------------------------------------------------------------------------------
# ANALYSE DE LA SERIE TEMPORELLE AVEC LA FONCTION BFASTMONITOR
#----------------------------------------------------------------------------------------------
# Now we apply the bfastmonitor function using a trend + harmon model with order 3 for the harmonics (i.e. seasonality modelling):
# --> FAIRE VARIER LE PARAMETRE "ORDER" et voir les différences 
# --> TESTER AVEC/SANS "trend" et "harmon"
bfastm <- bfastmonitor(tspx, response ~ trend + harmon, order = 3, start = c(2010,1)) # trend + harmon
bfastm <- bfastmonitor(tspx, response ~ trend         , order = 3, start = c(2010,1)) # trend only
bfastm <- bfastmonitor(tspx, response ~         harmon, order = 3, start = c(2010,1)) # harmon only

# GRAPHIQUE ISSU DE LA FONCTION BFASTMONITOR
plot(bfastm)

# SAUVEGARDE DU GRAPHIQUE
outdir<-paste0(getwd(),"/PLOT/",MODIS_Fixed_Sites_Subsets_ID,"/")
dir.create(outdir) #  Attention, le nom de répertoire correspondant au "Site ID" doit être créé avant l'exportation du graphique dans ce répertoire. Ce que fait cette ligne de commande.

# ATTENTION! CHANGER LE NOM DU GRAPHIQUE SELON LE PARAMETRAGE CHOISI!!
png(paste0(outdir,"pix_545_trend_harmo_ORDER_3_start2010_1.png"),width = 1000, height = 1000, units = "px")
plot(bfastm)
mtext(MODIS_Fixed_Sites_Subsets_ID, line=0.5)
dev.off()


#----------------------------------------------------------------------------------------------
# Critical questions for the BFAST AEO exercise:
#----------------------------------------------------------------------------------------------
# Confer le tutoriel d'origine: https://verbe039.github.io/BFASTforAEO/


#----------------------------------------------------------------------------------------------
# CI-DESSOUS: SECTION DU TUTORIEL EN LIGNE NON VUE AU COURS DU TP
#----------------------------------------------------------------------------------------------
# More information
# More information can be found on the BFAST website and in the BFAST papers mentioned on the website.
# seasonlity monitoring using harmonics
#----------------------------------------------------------------------------------------------
library(bfast)
## a demo ndvi time series:
ndvi <- ts(rowSums(simts$time.series))
tsp(ndvi) <- tsp(simts$time.series)
## input variable for the sinus and cosinus functions
f <- 23
w <- 1/f
tl <- 1:length(ndvi)
## 3th order harmonic model
co <- cos(2 * pi * tl * w)
si <- sin(2 * pi * tl * w)
co2 <- cos(2 * pi * tl * w * 2)
si2 <- sin(2 * pi * tl * w * 2)
co3 <- cos(2 * pi * tl * w * 3)
si3 <- sin(2 * pi * tl * w * 3)
# fit the seasonal model using linear regression
fitm<- lm(ndvi~co+si+co2+si2+co3+si3) 
predm <- fitted(fitm) ## predict based on the modelfit
plot(co, type = "l", ylab = "cos and sin")
lines(si, type = "l", lty = 2)
#----------------------------------------------------------------------------------------------
#create time series bfast on the 3th order harmonic function
predm <- ts(as.numeric(predm), start=c(2000,4), frequency=23) 
plot(ndvi, lwd = 3, col = "grey", ylab = "NDVI")
lines(predm, type = "l", col = "red") # fitted
#----------------------------------------------------------------------------------------------
