# M2 Sasha

## Réunion - 16/02/2024

- Objectifs du M2 :
- essayer de reproduire par simulation l'évolution entre 1950 et 2015 observée dans une parcelle de ST : est-on capable ? par quelle méthode ? pourquoi ?
- comparer les deux méthodes rangeshifter et SMA
- inconvénients et avantages des deux méthodes

- Choix de la parcelle
- Critères :
- monospécifique +
- dynamique++
- historique pastoralisme ?
- **Parcelle 716 (RNN Eyne)** ? 214 ? 370/524 (Orlu)

- Grandes tâches
- Etape 1 : sélection parcelle
- Etape 2 : identifier l'ensemble des informations disponibles sur le site (images, archives ST, etc) et sur l'espèce, se faire une idée de ce qu'il manque + données nécessaires
- Etape 3 : historique de la parcelle (évolution + déterminants -- pratiques RN) + gestion de la forêt ? replantation ?
- Etape 4 : Processus à simuler -> qu'est ce qu'on peut implementer ? [modèles / biblio]
- **Etape 5 : lancer les simulations ; comparer les résultats (rétrospective)**
- Etape 6 : choix des scénarios exploratoires

## Réunion - 02/04/2024*

/Présentation en avril devant labos/

Avancement surtout bibliographique (variables)

**Objectifs semaine 1:**
- informations sur la parcelle (topograhie, géologie, couverture + évolution)
- (+ explorer le sharedoc!)
- continuer bibliographie thématique (variables)
- faire l'exercice pratique présenté dans le manuel de Rangeshifter (prise en main)

**Variables** (liste en cours...)
- Climatique P,T, rayonnement -- min, max, moy, stress hydrique, vent + gradient orographique
- exposition
- Alt, pente, MOS, geol,
- Distance à l'océan
- Evenements ponctuels : feux ? mouvements de terrain ?
- Influences anthropiques : paturage, feux, entretien, plantation
- Forme de la treeline

[Notes :
Chevauchement des cycles phéno avec cdts météo
Non linéarié des phénomènes dans le temps
modification productivité
allogement durée saison croissance
avancement des événements phénologiques
Csq sequestration carbone CC Baisse seq Hausse CC]

*Recommandation biblio / tutoriels / manuels R ?*
- Carto R : https://rcarto.github.io/cartographie_avec_r/
- Géomatique avec R : https://rcarto.github.io/geomatique_avec_r/
- Initiation stat : https://tfeuillet.gitpages.huma-num.fr/cours/
- Régression : https://tfeuillet.gitpages.huma-num.fr/cours/cours_ed
- Geomorphons : https://rzine.fr/docs/20230425_geomorphon/index.html

## Notes utilisation du package Rangeshifter sur R

Notes accessibles ici :
- https://hackmd.io/m0w-hqeOQCy-QSDyumfURw
- https://docs.google.com/document/d/1G39S3ImcJNQedb6JhxmgZxd8Z78tTNTfR9fMFeaTGUw/edit?usp=sharing

https://rangeshifter.github.io/RangeshiftR-tutorials/overview_0.html
https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.05689

## Réunion 04/04/24
- Photo-interprétation
- Visio 9h30 vendredi 12

## Divers
Checker + savoir si peut être une approche pertinente ?
https://pro.arcgis.com/fr/pro-app/latest/tool-reference/spatial-statistics/how-presence-only-prediction-works.htm

## Réunion 12/04/2024
### Avancement et questions :
Données :
- check des données sur le sharedoc
- **Question** : quel MNT utilisé jusqu'ici ?
- RGE ALTI IGN 1m
- téléchargement :
- BD ORTHO : 2021, 2018, 2015,2012
- Photos historiques -- à géoréférencer
- 1953
- 1965
- 1966
- 1969
- 1985
- 1995
- 53, 65, 85, 95, 15, 18, 21
- **Questions** :
- Toutes les retenir ? (8 dates entre 1953 et 2015 (+2))
- Méthode pour treeline : digit ou classif type basée sur objet / SCP ?
- Quelle méthodo dans thèse de Déborah ?

Question autre : Comment prendre en compte directement le climat et l'impact anthropique dans RS ? Pas possible ou attente SMA ?

### Bilan réunion :
MNT : RGE 1m
Images : 53, 65, 85, 95, 15, 18, 21 (alternance 10 - 20 ans + précision après 2015)
- Faire géoréférencement
- Faire essai digit ou classification

RangeShifter : 2 étapes identifiées
- Traits spécifiques de l'espèce
- Telabotanica
- TRY Plant Trait Database (mondial)
- BROT (Bassin méditérannéen)
- (Données Matthieu Alpes du Sud)
- Mise en place de la carte habitat :
- Modèles aire de viabilité d'une espèce / quelles sont les cdts écologiques dans lesquelles l'espèce est présente / optimum écologique ==> une seule information (0 - 1 : 0 - 100)
- BioMod2 ou MaxEnt (ArcGis)
- Permet de prendre en compte plusieurs variables pour obtenir une carte de la qualité des habitats
- Limite : considère habitat vivable là où sont déjà implantés les individus (présence)
- Décisions
- StageStructure (3 stades : juvénils, jeune individu, mature)

test multicolinéarité ? (utilisé par https://environmentalsystemsresearch.springeropen.com/articles/10.1186/s40068-022-00248-6)

------

## Géoréférencement
Erreurs
| Date | Nb points | Méthode |Inverse |Inverse-avant|
|:-------- |:--------:| --------:|---------:|---------:|
|1962|94|Ajuster|0,019 | 0,007|
|1974|131|Ajuster|0,103| 0,082|
|1985|106|Ajuster|0,101|0,078|
|1995|36|Ajuster|0,008|0,014|

Ajustement sur base 1953

Garder :
- 1953
- 1962
- 1974
- 1985
- 1995
- 2012
- 2015
- 2021

-------------------
## Classification

Résultats variables...
2018 plûtot cohérent // 2021...compliqué !
Taux de validation pas très parlant ?



## à voir :
Deux types de données à intégrer :
- les données de présence de l'espèce (ici, les pins) : intégrées par une couche d'entités ponctuelles ? (Ce qui nécessite donc de créer la couche en "poitant" chaque pin à chaque date -- la classification d'occupation du sol étant faite en entités surfaciques)

Et des variables relatives à l'environnement :

1 - celles qui seraient intéresantes mais qui ne sont pas utiles à notre échelle de travail ou indisponibles (climat, ensoleillement, humidité du sol, utilisation du sol...)

**2 - celles qui varient dans le temps et dont on a accès à la donnée
-> occupation du sol - qui pourra être obtenue par la classification des images aériennes : sol nu / herbacé / forêt ouverte / fermée ?) + est ce que ces classes sont cohérentes ? **

3 - celles qui ne varient """pas""" dans le temps et dont on a l'accès ("pas" car ...pas d'informations historiques)
-> MNT => altitude, pente, exposition (+ géomorphons ?)

==> PAYSAGE DYNAMIQUE ?

## Réunion 26/04
- Validation du géoréférencement
- Discussion classification :
- Essai d'une classification non supervisée (résultats assez bons)
- Essayer classification supervisée emboitée
- Faire plus de classes puis les regrouper ensuite
- Toutes dates pour validation du modèle ultérieurement
- Discussion carte de viabilité :
- SOL : A faire sur 2021 (aide BRGM ?) : différents types de sols (éboulis, roches, sol nu, sol occ herbe, sol riche occ forêt)
- Choix 2021 car permet de tenir compte de l'emprise maximale connue
- Choix du sol et non de l'occupation du sol : type de sol (ou présence / absence) définit probabilité d'expansion des arbres
- Discussion autour de l'impossible intégration de données climatiques : postulat que le gradient orographique n'est pas un facteur limitant puisque des arbres sont présents dans les plus hautes altitudes de la parcelle
- RELIEF : Pentes, exposition, géomorphons
- Tester la multicolinéarité
- Pointer chaque arbre (2021)

--------------
Données accessibles ici : https://drive.google.com/drive/folders/1adbtS_DqthEYI9ovEQ5yE0KUX57RhfZr?usp=sharing

## Bilan 22/05


Bilan : à régler :
- [x] Modifier exposition
- [x] Reclassifier les rasters d'occupation du sol :
- Sol nu, roches, éboulis, végétation
- [x] Tests VIF
- [ ] Relancer les modèles
- [x] Checker comment faire pour insérer toutes les dates dans le modèle

## Test colinéarité
Y = pente

library(raster)
library(usdm)
library(terra)
#Chargement raster NS
expoNS <- rast("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/Expo_NS.tif")
expoNS

# Garder uniquement valeurs relatives à l'expo NS
expoNS <- (expoNS$Expo_NS_1)
expoNS
plot(expoNS)

#Chargement raster EO
expoEO <- rast("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/Expo_EO.tif")

# Garder uniquement valeurs relatives à l'expo EO
expoEO <- (expoEO$Expo_EO_1)
plot(expoEO)

# Passer de SpatRaster à RasterLayer
expoNS <- raster(expoNS)
expoEO <- raster(expoEO)
summary(expoNS)
par(mfrow = c(1,2))
plot(expoNS, main="Exposition Nord-Sud")
plot(expoEO, main ="Exposition Est-Ouest")

##Charger rasters : Geomorphons et Pentes
#Geomorphons 
#Raster
geom <- rast("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/Geomorph716.tif")
#Conversion  vers RasterLayers
library(raster)
geom <- raster(geom)
#Pentes
#Rasters
pentes <- rast("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/Pentes716.tif")
#Conversion 
pentes <- raster(pentes)

### TEST MULTICOLINEARITE
library(raster)
library(usdm)

#1 : Empiler les rasters
#1.1 : Corriger étendue des rasters
# 1.1.a : définir étendue de référence
etendue_ref <- geom
# 1.1.b : Resampler les autres rasters
RexpoEO <- resample(expoEO, etendue_ref)
RexpoNS <- resample(expoNS, etendue_ref)
Rpentes <- resample(pentes, etendue_ref)
#1.2 : verification étendues
print(extent(RexpoEO))
print(extent(RexpoNS))
print(extent(geom))
print(extent(Rpentes))

#2 : Empiler les raters
Stack <- stack(RexpoEO, RexpoNS, Rpentes, geom)

#3 : Extraction des valeurs de pixels + conversion vers df
Valeurs <- as.data.frame(getValues(Stack))

#3 : Calcul de la matrice de corrélation 
Matrice_cor <- cor(Valeurs, use = 'complete.obs')
print(Matrice_cor)

### Test VIF 
# Supp NA
ValeursNA <- na.omit(Valeurs)


# Créer un modèle
library(car)
lmP <- lm(ValeursNA$Pentes716~ ValeursNA$Expo_EO_1 + ValeursNA$Expo_NS_1 + 
           ValeursNA$Geomorph716)
summary(lmP)
# Calcul VIF
vifP <- vif(lmP)
vifP

NOMP <- c("Expo E-O", "Expo N-S", "Géomorphons")

plot(vifP, pch = 16, main = "VIF", xlab = "Variables", ylab = "Valeurs", 
     xlim=c(1,3), ylim=c(0,3), xaxt="n")
axis(side=1, at=1:length(NOMP), labels=NOMP)
abline(h=2, lty=2, col="red")
text(1:length(vifP), vifP, labels = round(vifP, 2), pos = 3, cex=1.2)


par(mfrow = c(2,4))
#VIF Relief + Types Sol
S1953 <- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/TYPE_SOL/Sol1953.tif")
RS1953 <- resample(S1953,etendue_ref)
Stack1953 <- stack(RS1953, RexpoEO, RexpoNS, Rpentes, geom)
Val1953 <- as.data.frame(getValues(Stack1953))
Matrice_cor1953 <- cor(Val1953, use = 'complete.obs')
print(Matrice_cor1953)
ValeursNA53 <- na.omit(Val1953)
lmP1953 <- lm(ValeursNA53$Pentes716 ~ ValeursNA53$Expo_EO_1 + ValeursNA53$Expo_NS_1
              + ValeursNA53$Geomorph716 + ValeursNA53$TYPE)
summary(lmP1953)
vifP1953 <- vif(lmP1953)
vifP1953


NOM2 <- c("ExpoE-O", "Expo N-S", "Géomorphons", "Type sol")

plot(vifP1953, pch = 16, main = "VIF1953", xlab = "Variables", 
     ylab = "Valeurs", xlim=c(1,4), ylim=c(0,3), xaxt="n")
axis(side=1, at=1:length(NOM2), labels=NOM2)
abline(h=2, lty=2, col="red")
text(1:length(vifP1953), vifP1953, labels = round(vifP1953, 2), pos = 3, cex=0.7)


S1962 <- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/TYPE_SOL/Sol1962.tif")
RS1962 <- resample(S1962,etendue_ref)
Stack1962 <- stack(RS1962, RexpoEO, RexpoNS, Rpentes, geom)
Val1962 <- as.data.frame(getValues(Stack1962))
Matrice_cor1962 <- cor(Val1962, use = 'complete.obs')
print(Matrice_cor1962)
ValeursNA62 <- na.omit(Val1962)
lmP1962 <- lm(ValeursNA62$Pentes716 ~ ValeursNA62$Expo_EO_1 + ValeursNA62$Expo_NS_1
               + ValeursNA62$Geomorph716 + ValeursNA62$TYPE)
summary(lmP1962)
vifP1962 <- vif(lmP1962)
vifP1962
plot(vifP1962, pch = 16, main = "VIF1962", xlab = "Variables", 
     ylab = "Valeurs", xlim=c(1,4), ylim=c(0,3), xaxt="n")
axis(side=1, at=1:length(NOM2), labels=NOM2)
abline(h=2, lty=2, col="red")
text(1:length(vifP1962), vifP1962, labels = round(vifP1962, 2), pos = 3, cex=0.7)

etc.



## Biomod2

### BIOMOD
library(biomod2)
library(sf)
library(raster)
library(terra)
library(ggplot2)
library(gridExtra)

#Charger données
Arbres2021 <- st_read("D:/A_STAGE/DATA_SIG/A_IMAGES/5_CENTROIDES/Arbres_2021.shp")
EO<- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/A_Export/EO_Q.tif")
GeomF <- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/A_Export/Geo_Q.tif")
Pentes <- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/A_Export/Pentes_Q.tif")
NS <- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/EL_RELIEF/A_Export/NS_Q.tif")
s21 <- raster("D:/A_STAGE/DATA_SIG/B_Carto_viab/Entrees/TYPE_SOL/Sol2021.tif")

GeomF <- as.factor(GeomF)
s21 <- as.factor(s21)

str(GeomF)
str(s21)

# Corriger étendue des rasters
#définir étendue de référence
etendue_ref <- GeomF
#: Resampler les autres rasters
EO <- resample(EO, etendue_ref)
NS<- resample(NS, etendue_ref)
Pentes <- resample(Pentes, etendue_ref)
s21 <- resample(s21, etendue_ref)

#: verification étendues
print(extent(EO))
print(extent(NS))
print(extent(GeomF))
print(extent(Pentes))

#Stack variables environnementales 
Stackbio <- stack(Pentes, NS, EO, GeomF, s21)
 summary(Stackbio)
 str(Stackbio)
 #Coordonnées
CoodArbres <- st_coordinates(Arbres2021)
#Df Présence 
DataPresence <- data.frame(longitude=CoodArbres[,1], latitude =CoodArbres[,2], 
                           species_presence = rep(1, nrow(CoodArbres)))

#Formatage des données pour Biomod1 
PresenceArbres <- BIOMOD_FormatingData(resp.var=DataPresence$species_presence,
                                       expl.var = Stackbio, 
                                       resp.xy = DataPresence[,c('longitude', 'latitude')], 
                                       resp.name = "ArbresBiomod", #Nom du projet
                                       PA.nb.rep = 1, #nb de répétition de points de pseudoabs
                                       PA.nb.absences = 530,#nb de points, ici ratio 3:1
                                       PA.strategy='disk',
                                       PA.dist.min = 10) #Aléatoire (hors points de présence)


PresenceArbres
summary(PresenceArbres)
library(tidyterra)
library(ggtext)


#Paramètres des modèles
OptDef <- BIOMOD_ModelingOptions(
  MAXENT= list(path_to_maxent_jar = "D:/A_STAGE/DATA_SIG/B_Carto_viab/A_TestMulticol
  /maxent/maxent/maxent.jar"))

  
# Lancement 
Models <- BIOMOD_Modeling(
  bm.format = PresenceArbres,
  modeling.id = "Modeles",
  models=c('GLM','GAM','GBM','CTA','ANN','SRE','FDA','MARS','RF','MAXENT','MAXNET','XGBOOST'),
  bm.options=OptDef,
  CV.strategy = 'random',
  CV.nb.rep = 10,
  CV.perc = 0.7,
  metric.eval = c('ROC', 'TSS', 'KAPPA'),
  var.import = 3,
  seed.val = 42)
  
 # Evaluation 
eval <- get_evaluations(Models)

eval_df <- as.data.frame(eval)

print(colnames(eval_df))
head(eval_df)

tss_scores <- eval_df %>%
  filter(metric.eval == "TSS") %>%
  select(full.name, PA, run, algo, cutoff, sensitivity, specificity, calibration, validation)

tss_scores$run <- factor(tss_scores$run, levels = c("RUN1", "RUN2", "RUN3", "RUN4", "RUN5", "RUN6", "RUN7", "RUN8", "RUN9", "RUN10"))
couleurs <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf")

point_tss <- ggplot(tss_scores, aes(x = run, y = calibration, color = algo)) +
  geom_point(size = 2, position = position_jitter(width = 0.2, height = 0)) +
  geom_text(aes(label = sprintf("%.3f", calibration), color = algo), size = 3, vjust = -1.5, position = position_jitter(width = 0.2, height = 0)) +
  labs(title = "Scores de calibration par Run et par Modèle",
       x = "Run",
       y = "Score de calibration",
       color = "Algorithme") +
  scale_color_manual(values = couleurs) +  # Utiliser les couleurs personnalisées
  scale_shape_manual(values = 1:length(unique(tss_scores$algo)), breaks = unique(tss_scores$algo)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "right")   # Positionner la légende à droite
print(point_tss)
bm_PlotEvalMean(Models, group.by="algo", metric.eval = c("ROC", "TSS"), xlim=c(0.8,1), ylim=c(0.8, 1))


# Courbes de réponse 
CourbesRF <- c("ArbresBiomod_PA1_RUN1_RF", 
                    "ArbresBiomod_PA1_RUN2_RF",
                    "ArbresBiomod_PA1_RUN3_RF", 
                    "ArbresBiomod_PA1_RUN4_RF", 
                    "ArbresBiomod_PA1_RUN5_RF", 
                    "ArbresBiomod_PA1_RUN6_RF", 
                    "ArbresBiomod_PA1_RUN7_RF",
                    "ArbresBiomod_PA1_RUN8_RF",
                    "ArbresBiomod_PA1_RUN9_RF",
                    "ArbresBiomod_PA1_RUN10_RF")

RF <- bm_PlotResponseCurves(bm.out = Models, models.chosen = CourbesRF)