Impact_PAs_Mada_07.Rmd

---
title: "Impact des aires protégées sur la déforestation à Madagascar"
author: "Florent Bédécarrats, Jeanne de Montalembert, Martin Ferry et Kenneth Houngbedji"
output:
  html_document:
    code_folding: hide
  pdf_document: 
    toc: yes
date: '2022-07-12'
editor_options:
  chunk_output_type: console
  markdown: 
    wrap: 80
bibliography: bibliography_tany_vao.bib
---

# Environnement et paramétrages

L'analyse est réalisée en R, qui est à la fois un logiciel et un langage open
sources dédiés à l'analyse de données. Les traitements sont réalisés en
Rmarkdown (à passer en Quarto). Le même code source peut générer un rendu en LaTeX/PDF, HTML ou Word.

```{r Evite affichage erreurs et messages pour tous les blocs}
# On supprime systématiquement les messages d'erreur ou d'info dans le document
# final généré (ils apparaissent dans la console) 
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
```

```{r Installation et chargement des librairies requises}
# # Le package est en cours de développement, toujours installer la version en cours
remotes::install_github("mapme-initiative/mapme.biodiversity", 
                        upgrade = "always")

librairies_requises <- c( # On liste les librairies dont on a besoin
  "dplyr", # Pour faciliter la manipulation de données tabulaires
  "tidyr", # Pour reformater les données (pivots...)
  "sf", # Pour faciliter la manipulation de données géographiques
  "wdpar", # Pour télécharger simplement la base d'aires protégées WDPA
  "tmap", # Pour produire de jolies carte
  "geodata", # Pour télécharger simplement les frontières administratives
  "tidygeocoder", # pour obtenir les coordo GPS d'un point à partir de son nom
  "maptiles", # Pour télécharger des fonds de carte 
  "purrr", # Pour utiliser des formes fonctionnelles de programmation (ex. map)
  "mapme.biodiversity", # Acquisition et traitement des données du projet
  
  "plm", # Linear Models for Panel Data and robust covariance matrices
  "stargazer", # Reformater de manière plus lisible les résumé des régressions
  "MatchIt", # Pour le matching
  # "glm", # Modèles linéaires généralisés (pour le PSM)
  "optmatch", # Fonctions d'optimisation du matching
  "cobalt") # Tables et graphs d'équilibre des groupes de matching
  
# On regarde parmi ces librairies lesquelles ne sont pas installées
manquantes <- !(librairies_requises %in% installed.packages())
# On installe celles qui manquent
if(any(manquantes)) install.packages(librairies_requises[manquantes])
# On charge toutes les librairies requises
invisible(lapply(librairies_requises, require, character.only= TRUE))
```

```{r Paramtères modifiables pour le document}
# Système de coordonnées géographiques utilisées pour le projet : EPSG:29739
mon_scr <- "EPSG:29739" # correspondant à Tananarive / UTM zone 39S
# Surface des hexagones en km2
taille_hex <- 5
# Taille des titres des cartes
taille_titres_cartes = 0.8
# on crée un dossier de données si pas déjà disponible
dir.create("data_s3")
# Désactiver les notations scientifiques
options(scipen =999)
```

Les analyses sont menées sur la plateforme SSP Cloud, mise à disposition par
l'INSEE pour les data scientist travaillant pour des administrations publiques.
Il s'agit d'un cluster Kubernetes disposant d'une interface simplifiée
permettant à l'utilisateur de configurer, lancer et administrer facilement des
pods. Le code est conçu pour s'exécuter de la même manière en local sur un PC,
mais la préparation des données sera certainement beaucoup plus longue à
exécuter.

# Idées pour la structuration du cours

Les ateliers dureront 5 jours, suivis d'1/2 journée de restitution, avec
25 participants par atelier (\~125 au total). On veillera à diviser l'atelier en
sous-groupes en essayant d'avoir dans chaque groupe au moins un participant
ayant quelques compétences économétrie et un participant capable de déchiffrer
un texte en anglais.

Pour capter et retenir l'attention d'apprenants ayant des niveaux hétérogène, la
proposition serait d'alterner des séances très spécifiques/concrètes avec des
séances plus théoriques/méthodologiques. On pourrait ainsi avoir une approche
itérative, avec une durée des sections qui reste à déterminer (entre 1/2 et 1
journée ?). Le programme qui suit n'est qu'une proposition pour alimenter la
réflexion collective :

**Séance 1 - Définition des objectifs de l'évaluation**

Qu'est-ce qu'une évaluation d'impact au sens économétrique du terme ? Quels sont
les autres types d'évaluation et quelle est la différence ? Place centrale d'une
question évaluative précise pour les évaluations d'impact : impact de quoi
(intervention), sur qui (groupe de traitement), sur quoi (variable de résultat).
Les EI évaluent un mode d'intervention, pas un projet en particulier : enjeux de
validité interne et validité externe.

-   Exercice : Formuler des questions se prêtant à une évaluation d'impact sur
    des enjeux de conservation.

> Synthèse : introduction au formalisme des équation et aux DAG.

**Séance 2 - recherche bibliographique**

Pas pertinent de mener de nouvelles EI sur des sujets déjà bien étudiés :
important de commencer par une bonne revue de littérature. Présentation des
outils et méthodes de revue de littérature appliquées aux EI

Google Scholar : présentation de l'outil et recommandations pour trouver des
évaluations d'impact pertinentes.

-   Exercice : utiliser GS pour mener une revue de littérature sur le thème de
    l'impact des aires protégées, au niveau mondial, régional (ex. Afrique) ou
    national (Madagascar)

Présentation des bases d'évaluation d'impact (p. ex. Campbell Collaboration).
Focus sur la base d'évaluations d'impact du 3IE.

Exercice : explorer l'evidence gap map du sur la conservation des forêts :
<https://gapmaps.3ieimpact.org/evidence-maps/forest-conservation-gap-map> -
revue de littérature.

-   

> Synthèse théorique :

**Séance 3 - étude d'articles**

Les formateurs distribuent aux participant un article court sur l'évaluation de
l'impact des aires protégées sur la déforestation. Les apprenants le lisent. Les
formateurs en font une synthèse en : discipline des auteurs, périmètre de
l'étude, données mobilisées, unités comparées, taille d'échantillon, méthode,
spécifications du modèle (variable de traitement, de contrôle, et de résultat),
résultat.

-   Exercice : les étudiants se voient chacun remettre un article (en français)
    et remp

> Synthèse théorique : tableau complété de revue de la littérature.

**Séance 3 - Sources de données**

Exercice : recherche de source de données pertinentes pour étudier l'impact des
aires protégées sur la déforestation

**Traitement des données**

**Analyse des données** - Les logiciels et les ressources d'auto-formation. R
Stata Python Google Earth Engine

Concevoir un questionnaire

**Interprétation des résultats**

**Traitement des données**

```{r}
revue_litt <- tibble(
  `Référence` = character(),
  `Titre` = character(),
  `Discipline des auteurs` = character(),
  `Périmètre d'analyse` = character(),
  `Données mobilisées` = character(),
  `Unités comparées` = character(),
  `Taille d'échantillon` = character(),
  `Méthode d'attibution` = character(),
  `Variable de traitement` = character(),
  `Variables de contrôle` = character(),
  `Variable de résultat` = character(),
  `Résultats`= character())

revue_litt <- revue_litt %>% add_row(
  `Référence` = character(),
  `Titre` = character(),
  `Discipline des auteurs` = character(),
  `Périmètre d'analyse` = character(),
  `Données mobilisées` = character(),
  `Unités comparées` = character(),
  `Taille d'échantillon` = character(),
  `Méthode d'attibution` = character(),
  `Variable de traitement` = character(),
  `Variables de contrôle` = character(),
  `Variable de résultat` = character(),
  `Résultats`= character() 
  )

```

On pourrait partir par exemple d'un article paru dans une revue généraliste
évaluant cet impact à l'échelle internationale\
L'idée pourrait être de

Lire 1 ou 2 articles et étudier la manière dont ils ont abordé le problème. Tour
d'horizon théorique de s

On réutilise le code fourni par Johannes Schielein: Jochen Kluve, Johannes
Schielein, Melvin Wong, Yota Eilers, The KfW Protected Areas Portfolio: a
Rigorous Impact Evaluation, KfW, 2022-07-08.

On s'appuie sur le package R {mapme.biodiversity}, développé par la KfW dans le
cadre de l'initiative commune MAPME qui associe la KfW et l'AFD. Le package
{mapme.biodiversity} facilite l'acquisition et la préparation d'un grand nombre
de données (CHIRPS, Global Forest Watch, FIRMS, SRTM, Worldpop...) et calculer
un grand nombre d'indicateurs de manière harmonisée (active_fire_counts, biome
classification, land cover classification, population count, precipitation, soil
properties, tree cover loss, travel time...). Une documentation riche est
disponible [sur le portail dédié du package en
question](https://mapme-initiative.github.io/mapme.biodiversity/index.html).

On mobilise aussi les codes d'analyse d'impact développés par la même équipe et
mise à disposition dans le dépôt Github:
<https://github.com/openkfw/mapme.protectedareas> Le code développé par l'équipe
est assez complexe. A des fins pédagogiques et pour s'assurer qu'on l'a bien
compris, on propose ici une version simplifiée (en cours de développement)

# Prépration des données

Les données spatialisées à croiser son, pour certaines, des données vectorielles
(aires protégées, frontières administratives) et, pour d'autres, des données
matricielles ("raster data", en anglais).

## Aires protégées

Les données d'aires protégées sont issues de la base WDPA, consultable en ligne
sur protectedplanet.org.

```{r Préparation aires protégées}
# Ce qui suit jusqu'à la commande "save" ne s'execute que si le résultat n'a pas
# déjà été généré lors d'une exécution précédente.
if (file.exists("data_s3/aires_prot_mada.rds")) {
  load("data_s3/aires_prot_mada.rds")
} else {
  # Téléchargement et chargement dans R des données d'aires protégées malgaches
  aires_prot_mada <- wdpa_fetch("Madagascar", wait = TRUE,
                                download_dir = "data_s3/WDPA") %>%
    st_transform(crs = mon_scr) %>%
    filter(STATUS != "Proposed") %>%
    filter(DESIG != "Locally Managed Marine Area", DESIG != "Marine Park") 
  
  # Téléchargement du contour des zones émergées de Madagascar
  contour_mada <- gadm(country = "Madagascar", resolution = 1, level = 0,
                       path = "data_s3/GADM") %>%
    st_as_sf() %>% 
    st_transform(crs = mon_scr)
  # On sauve les objets créés pour ne pas avoir à refaire cette étape
  save(aires_prot_mada, contour_mada, file = "data_s3/aires_prot_mada.rds")
}

2006 - 2006 %% 5

tb <- aires_prot_mada %>%
  filter(STATUS != "Proposed", MARINE != 2) %>%
  mutate(decennie_creation = STATUS_YR -  STATUS_YR %% 10,
       strict = IUCN_CAT %in% c("I", "II", "III", "IV"),
       surface_terrestre = REP_AREA - REP_M_AREA) %>%
  group_by(decennie_creation, strict) %>%
  summarise(N = n(),
            aire_totale = sum(surface_terrestre, na.rm = TRUE))

# On génère un rendu cartographique
tm_shape(contour_mada) +
  tm_polygons() +
  tm_shape(filter(aires_prot_mada)) + 
  tm_polygons(col = "IUCN_CAT", alpha = 0.6, title = "Catégorie IUCN") +
  # NB : on note les positions en majuscules quand on veut coller aux marges
  tm_credits("Sources: WDPA et GADM", position = c("RIGHT", "BOTTOM"),
             size = 0.6) +
  tm_layout(main.title = "Aires protégées de Madagascar",
            # NB : position en minuscules pour laisser un espace avec la marge
            main.title.position = c("center", "top"),
            main.title.size = taille_titres_cartes,
            legend.position = c("left", "top"),
            legend.outside = TRUE)
```

Certaines améliorations doivent encore être apportées, pour préciser notamment
la date de création ou le statut de certaines aires =\> A travailler avec Jeanne
notamment.

Il faut aussi s'assurer qu'on filtre bien les entitées analysées selon un criète
pertinent. Actuellement, on ne garde que les aires qui ont encore un statut
"proposed" et on exclut les aires marines. Il pourrait toutefois sembler utile
d'écarter les aires dont le statut de protection est considéré comme trop
faible. En particulier, les aires de type "UNESCO MAB Biosphere reserves" sont 
exclues des analyses réalisées par Schielein et al., qui 
ont un statut faible de protection et des frontières floues. Les entités de niveau V 
Il pourrait aussi être pertinent de ne garder que les aires protégées
comportant un niveau minimum de couvert forestier : autrement, cela signifie que
la forêt n'est pas un habitat pertinent pour les écosystèmes que la démarche de
conservation cherche à protéger dans cette aire.



## Données satellitaires

Ici le package mapme.biodiversity développé par la KfW est particulièrement
utile pour l'analyse. Il automatise en large partie le processus d'acquisition
de données brutes issu de sources divers et le calcul d'indicateurs pour des
périmètres définies (ici, les 120 612 hexagones du maillage du territoire
malgache). Ce processus est toutefois très gourmand en ressources et on l'a
réalisé sur un environnement de calcul haute performance (la plateforme SSP
Cloud de l'INSEE). Les résultats de ces traitements ont été enregistrés et il ne
semble pas pertinent/utile de demander aux apprenants de le refaire, ce serait
beaucoup trop long.

```{r Données satellitaires}
# Ce qui suit jusqu'à la commande "save" ne s'execute que si le résultat n'a pas
# déjà été généré lors d'une exécution précédente.
if (file.exists("data_s3/grille_mada_donnees_raster.rds")) {
  load("data_s3/grille_mada_donnees_raster.rds")
} else {
  
  # Création d'un maillage du territoire émergé --------------------------------
  
  # On crée un cadre autour des aires protégées du pays
  cadre_autour_mada = st_as_sf(st_as_sfc(st_bbox(aires_prot_mada)))
  
  # Cellules de 5km de rayon
  surface_cellule <- taille_hex * (1e+6)
  taille_cellule <- 2 * sqrt(surface_cellule / ((3 * sqrt(3) / 2))) * sqrt(3) / 2
  grille_mada <- st_make_grid(x = cadre_autour_mada,
                              cellsize = taille_cellule,
                              square = FALSE)
  # On découpe la grille pour ne garder que les terres émergées
  cellules_emergees <- st_intersects(contour_mada, grille_mada) %>%
    unlist()
  grille_mada <- grille_mada[sort(cellules_emergees)] %>%
    st_sf()
  
  # Traitement des données satellitaires avec {mapme.bidiversity}---------------
  
  # Constitution d'un portefeuille (voir la documentation)
  grille_mada <- init_portfolio(x = grile_mada, 
                                years = 2000:2020,
                                outdir = "data_s3/mapme",
                                cores = 24,
                                add_resources = TRUE,
                                verbose = TRUE)
  
  # Acquisition des données satellitaires requises (rasters) ------------------- 
  # Données d'accessibilité de Nelson et al. (2018)
  grille_mada <-  get_resources(x = grille_mada, resource = "nelson_et_al",  
                                range_traveltime = "5k_110mio")
  # Données de qualité des sols (uniquement teneur )
  grille_mada <-  get_resources(x = grille_mada,
                                resources = "soilgrids",  layers = "clay", 
                                depths = "5-15cm", stats = "mean")
  # Données sur le couvert forestier de Global Forest Watch
  grille_mada <- get_resources(x = grille_mada, 
                               resources = c("gfw_treecover", "gfw_lossyear", 
                                             "gfw_emissions"))
  # Modèle numérique de terrain SRTM de la NASA
  grille_mada <- get_resources(x = grille_mada, resource = "nasa_srtm")
  # Données de feux
  grille_mada <- get_resources(x = grille_mada, resource = "nasa_firms",
                               instrument = "MODIS")
  
  # Calcul des indicateurs -----------------------------------------------------
  
  # Indicateurs d'accessibilité
  grille_mada <- calc_indicators(x = grille_mada,
                                 "traveltime",  stats_accessibility = "mean",
                                 engine = "extract")
  # Indicateurs de sols
  
  grille_mada <- calc_indicators(x = grille_mada,
                                 "soilproperties", stats_soil = "mean", 
                                 engine = "extract")
 
   # Indicateurs de couvert forestier
  grille_mada <- calc_indicators(x = grille_mada,
                                 indicators = "treecover_area_and_emissions", 
                                 min_cover = 10, min_size = 1)
  # Indicateurs de relief de terrain
  grille_mada <- calc_indicators(x = grille_mada,
                               indicators = c("tri", "elevation"),
                               stats_tri = "mean", stats_elevation = "mean")
  # Indicateurs d'incendies
  grille_mada <- calc_indicators(x = grille_mada,
                                 "active_fire_counts")
  grille_mada <- calc_indicators(x = grille_mada,
                                 "active_fire_properties")
  
  # Sauvegarde du résultat
  save(grille_mada, file = "data_s3/grille_mada_donnees_raster.rds")
}
```

Le maillage est trop fin pour être visible à l'échelle du pays, mais on peut
l'observer en zoomant sur une zone spécifique.

```{r Carte grille mada}
# On compte le nombre d'hexagones
n_hex <- length(grille_mada)
# Carte pour visualiser le résultat --------------------------------------------

## Carte de droite : zoom sur une zone spécifique-------------------------------
# On part d'un dataframe contenant une adresse
nom_centre_zoom <- "Maroantsetra"
zoom_centre <- data.frame(address = nom_centre_zoom) %>%
  geocode(address, method = "osm") %>% # on retrouve sa localisation xy
  select(long, lat) %>% # on ne garde que le xy
  as.numeric() %>% # qu'on passe en format numérique attendu par st_point
  st_point() %>% # On le spécifie en point
  st_sfc(crs = "EPSG:4326") 

# On crée une boîte de 100km 
zoom_boite <- zoom_centre %>% # On repart du centre
  st_buffer(dist = 50000) %>% # On crée un cercle de 50km de rayon
  st_make_grid(n = 1) 

# On filtre les alvéoles pour ne garder que celles qui sont dans le zoom
grille_zoom <- st_intersection(grille_mada, zoom_boite)

# On télécharge un fond de carte pour la carte de droite
fond_carte_zoom <- get_tiles(zoom_boite, provider = "Stamen.Terrain", 
                             zoom = 10, crop = TRUE)
# On génère la carte de droite
carte_zoom <- tm_shape(fond_carte_zoom) + 
  tm_rgb() +
  tm_shape(grille_zoom) +
  tm_borders() +
  tm_shape(zoom_boite) + 
  tm_borders(col = "red") +
  tm_layout(frame = FALSE,
            main.title = paste("Zoom sur la zone de", nom_centre_zoom),
            main.title.size = taille_titres_cartes) +
  tm_credits(get_credit("Stamen.Toner"),
             bg.color = "white",
             align = "right",
             position = c("right", "BOTTOM"))

## Carte de gauche : simple à réaliser mais hexagones non visibles -------------
carte_grille <- tm_shape(grille_mada) +
  tm_polygons() + 
  tm_shape(zoom_boite) +
  tm_borders(col = "red") +
  tm_layout(frame = FALSE) +
  tm_layout(main.title = paste("Découpage en", n_hex,
                               "hexagones de", taille_hex*2, "km2"),
            main.title.size = taille_titres_cartes)

# Assemblage des deux cartes ---------------------------------------------------
tmap_arrange(carte_grille, carte_zoom, ncol = 2) 
```

On peut également représenter les différentes valeurs des indicateurs générés à
partir des données satellitaires.

```{r Synthèse données satellitaires, fig.fullwidth = TRUE}
if (file.exists("data_s3/grille_mada_summary.rds")) {
  load("data_s3/grille_mada_summary.rds")
} else {
  grille_mada_summary <- grille_mada %>%
    # On met à plat les données de distance
    unnest(cols = c(traveltime, soilproperties, tri, elevation),
           names_repair = "universal") %>%
    select(-distance, -layer, -depth, -stat,  -active_fire_counts, 
           -active_fire_properties) %>%
    rename(distance_minutes_5k_110mio = minutes_mean, mean_clay_5_15cm = mean) 
  
  grille_mada_summary <- grille_mada_summary %>%
    unnest(cols = treecover_area_and_emissions) %>%
    pivot_wider(names_from = "years", values_from = c("treecover", "emissions")) %>%
    mutate(var_treecover = (treecover_2020 - treecover_2000)/treecover_2000,
           sum_emissions = rowSums(across(starts_with("emission")), na.rm = T)) %>%
    rename(init_treecover_2000 = treecover_2000) %>% # pour le garder
    select(-starts_with("treecover"), -starts_with("emission")) %>%
    rename(treecover_2000 = init_treecover_2000) %>%
    relocate(geometry, .after = last_col())
  
  save(grille_mada_summary, file = "data_s3/grille_mada_summary.rds")
}

carte_acces <- tm_shape(grille_mada_summary) +
  tm_fill("distance_minutes_5k_110mio",
          title = "Distance ville (>5K hab)",
          palette = "Oranges",
          style = "fisher",
          n = 8,
          legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE,
            # legend.title.size = 0.8,
            # legend.text.size = 0.6,
            legend.hist.width = 1,
            legend.hist.height = 1)

carte_sol <- tm_shape(grille_mada_summary) +
  tm_fill("mean_clay_5_15cm",
          title = "Sol argileux (5-15cm prof)",
          palette = "YlOrBr",
          n = 8,
          legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE,
            # legend.title.size = 0.8,
            # legend.text.size = 0.6
            legend.hist.width = 1,
            legend.hist.height = 1)

carte_TRI <- tm_shape(grille_mada_summary) +
  tm_fill("tri_mean",
          title = c("Terrain accidenté (TRI)"),
          palette = "Blues",
          n = 8,
          legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE,
            # legend.title.size = 0.8,
            # legend.text.size = 0.6,
            legend.hist.width = 1,
            legend.hist.height = 1)

carte_alt <- tm_shape(grille_mada_summary) +
  tm_fill("elevation_mean",
          title = "Altitude",
          palette = "Purples",
          n = 8,
          legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE,
            # legend.title.size = 0.8,
            # legend.text.size = 0.6,
            legend.hist.width = 1,
            legend.hist.height = 1)

carte_cover <- graph_alt <- tm_shape(grille_mada_summary) +
  tm_fill("treecover_2000",
          title = "Couvert arboré en 2000",
          palette = "Greens",
          n = 8,
          legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE,
            # legend.title.size = 0.8,
            # legend.text.size = 0.6,
            legend.hist.width = 1,
            legend.hist.height = 1)

carte_loss <- graph_alt <- tm_shape(grille_mada_summary) +
  tm_fill("var_treecover",
          title = "Perte couvert (2000-2020)",
          palette = "Reds",
          n = 8,
          legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE,
            # legend.title.size = 0.8,
            # legend.text.size = 0.6,
            legend.hist.width = 1,
            legend.hist.height = 1)

tmap_arrange(carte_acces, carte_sol, 
             carte_alt, carte_TRI, 
             carte_cover, carte_loss,
             ncol = 2, nrow = 3) 

```

On notera que plusieurs autres indicateurs peuvent être calculés à partir du
pabkage mapme.biodiversity:

-   active_fire_counts: Calculate active fire counts based on NASA FIRMS
    polygons\
    active_fire_properties: Calculate active fire properties based on NASA FIRMS
    polygons\
-   biome: Calculate biomes statistics (TEOW) based on WWF\
-   drought_indicator: Calculate drought indicator statistics\
-   ecoregion: Calculate terrestrial ecoregions statistics (TEOW) based on WWF\
-   landcover: Calculate area of different landcover classes\
-   mangroves_area: Calculate mangrove extent based on Global Mangrove Watch
    (GMW)\
-   population_count: Calculate population count statistics (Worldpop)
-   precipitation_chirps: Calculate precipitation statistics based on CHIRPS\
-   precipitation_wc: Calculate precipitation statistics\
-   soilproperties: Calculate Zonal Soil Properties\
-   temperature_max_wc: Calculate maximum temperature statistics\
-   temperature_min_wc: Calculate minimum temperature statistics based on
    WorldClim\
-   traveltime: Calculate accessibility statistics\
-   treecover_area: Calculate treecover statistics\
-   treecover_area_and_emissions: Calculate treeloss statistics\
-   treecoverloss_emissions: Calculate emission statistics\
-   tri: Calculate Terrain Ruggedness Index (TRI) statistics

## Croisement des données d'aires protégées et satellitaires

On peut maintenant associer les données d'aires protégées aux hexagones afin de
les croiser avec les indicateurs issus des données satellitaries déjà calculés
pour ces hexagones.

```{r Jointure aires protégées et données satellitaires}
if (file.exists("data_s3/grille_mada_summary_AP.rds")) {
  load("data_s3/grille_mada_summary_AP.rds")
} else {
  # Le code suivant va asocier les hexagones aux aires protégées en se référant
  # aux AP par leur rang dans la table des AP. On voudra plutôt leur identifiant, 
  # alors on crée une table d'équivalence rang/identifiant 
  aires_prot_mada_rang_id <- aires_prot_mada %>%
    st_drop_geometry() %>% # Enlève l'information spatiale
    mutate(AP_ligne = row_number()) %>% # Intègre le numéro de ligne dans un champ
    select(AP_ligne, WDPAID) # On ne garde que le numéro de ligne et l'identifiant
  
  # Pour chaque hexagone, on va maintenant identifier s'ils touchent ("intersect")
  # ou s'ils sont strictiement inclus dans ("within") une aire protégé
  grille_mada_summary_AP <- grille_mada_summary %>%
    st_transform(crs = mon_scr) %>%
    mutate(AP_ligne = st_intersects(., aires_prot_mada), # liste des n° de lignes d'AP qui recoupent
           AP_ligne = map(AP_ligne, 1), # On extrait le 1° élément de la liste (toutes n'ont qu'1 élément)
           AP_ligne = as.integer(as.character(AP_ligne))) %>%  # formattage en numérique
    left_join(aires_prot_mada_rang_id, by = "AP_ligne") %>% # récupère l'id de l'AP
    rename(WDPAID_touche = WDPAID) %>% # on renomme pour différentier
    mutate(AP_ligne = st_within(., aires_prot_mada),
           AP_ligne = map(AP_ligne, 1),
           AP_ligne = as.integer(as.character(AP_ligne))) %>%
    left_join(aires_prot_mada_rang_id, by = "AP_ligne") %>%
    rename(WDPAID_inclus = WDPAID) %>%
    select(-AP_ligne) 
  
  grille_mada_summary_AP <- grille_mada_summary_AP %>%
    st_sf() %>%
    mutate(position_ap = ifelse(is.na(WDPAID_touche), "Extérieur",
                                ifelse(!is.na(WDPAID_inclus), "Intérieur",
                                       "Frontière"))) %>%
    relocate(geometry, .after = last_col()) 
  save(grille_mada_summary_AP, file = "data_s3/grille_mada_summary_AP.rds")
  haven::write_dta(st_drop_geometry(grille_mada_summary_AP), 
                   path = "data_s3/grille_mada_summary_AP.dta")
}

# Une vue après classification
tm_shape(grille_mada_summary_AP) +
  tm_fill(col = "position_ap", title = "par rapport aux aires protégées") +
  tm_layout(main.title = "Localisation des hexagones",
            # NB : position en minuscules pour laisser un espace avec la marge
            main.title.position = c("center", "top"),
            main.title.size = taille_titres_cartes,
            legend.position = c("left", "top"),
            legend.outside = FALSE)
```

```{r eval=FALSE, include=FALSE}
# Sauvegarde sur le serveur distant pour éviter de télécharger à chaque fois
aws.s3::s3sync(path = "data_s3",
               bucket = "fbedecarrats",
               prefix = "diffusion/deforestation_madagascar/data_s3/",
               create = FALSE,
               region = "",
               verbose = FALSE)
```

En plus d'un format natif R (rds), on a aussi enregistré l'export au format
Stata (.dta)

# Première approche d'appariement "naïve"

Une procédure détaillée est proposée dans
<https://github.com/openkfw/mapme.protectedareas>

On commence ici par une approche naïve, dans le sens où on apparie simplement
les zones dans les aires protégées avec les zones hors aires protégées pour
expliquer le principe du matching ("appariement", en français). On verra ensuite
que cette approche est trop simpliste pour être valide et qu'il faut réfléchir à
la population cible, aux variables d'appariement et au recouvrement entre les
groupes de traitement et de contrôle.

```{r}
# On référence le nom des variables qui vont servir à l'analyse
variables_analyse <- c("assetid","treatment","distance_minutes_5k_110mio",
                       "tri_mean", "elevation_mean", "mean_clay_5_15cm",
                       "treecover_2000", "var_treecover")
# On renomme le ficher 'df' (dataframe) : plus concis dans les commandes ensuite
df <- grille_mada_summary_AP %>%
  # On supprime toutes les lignes pour lesquelles au moins 1 valeur variable 
  # est manquante parmi les variables d'analyse
  drop_na(any_of(variables_analyse)) %>%
  mutate(treatment = position_ap == "Intérieur")
  

# Get propensity scores
glm_out <- glm(treatment ~ 
                 distance_minutes_5k_110mio + 
                 mean_clay_5_15cm + 
                 tri_mean +
                 elevation_mean + 
                 treecover_2000,  # Très étrange
               family = binomial(link = "probit"),
               data = df)

stargazer(glm_out,
          summary = TRUE,
          type = "text",
          title = "Probit regression for matching frame ")

m_out <- matchit(treatment ~ 
                   distance_minutes_5k_110mio + 
                   mean_clay_5_15cm + 
                   tri_mean +
                   elevation_mean + 
                   treecover_2000,
                 data = df,
                 method = "nearest",
                 replace = TRUE,
                 # exact = ~ as.factor(NAME_0),
                 distance = "glm", 
                 discard = "both", # common support: drop units from both groups 
                 link = "probit")

print(m_out)
print(summary(m_out, un = FALSE))
bal_table <- bal.tab(m_out, un = TRUE)
print(bal_table)
m_data <- match.data(m_out)

```

Première approche : critiquer la méthode employée ici.

# Pistes d'amélioration

-   Exclure les Aires protégées avant 2000, voire 2003.
-   On pourrait éventuelement prendre comme variables de contrôle le couvert
    forestier en 2000 et le taux de couverture entre 2000 et 2003.
-   Gros problème : toutes les aires protégées créées à partir de 2003 n'ont pas
    de
-   Schielein et al. excluent UNESCO MAB Biosphere reserves: pourquoi ?

# A faire

-   Passer le doc en Quarto\
-   Préciser les types d'aires protégées à conserver\
-   Corriger/enrichir les métadonnées des aires protégées avec les informations
    recueillies par Jeanne\
-   Ajouter une biblio\
-   Vérifier la préséance lorsque un hexa recoupe plusieurs AP (rang IUCN en
    premier)
-   Expliciter les unités
-   Expliciter le choix de la teneur en argile du sol pour l'analyse