2  Caractéristiques spatiales

Elle finit termine par enregistrer le jeu de données produit sur la machine qui exécute le code. Ce bloc ne s’exécute que si le jeu de données résultant n’est pas détecté sur la machine. Si le jeu de données résultant du script précédent est déjà disponible sur la machine, alors le bloc précédent ne s’exécute pas et celui qui s’exécute est le suivant

Code 
library(tidyverse)
library(mapme.biodiversity)
library(sf)

if (file.exists("data/Vahatra_poly.rds")) {
  load("data/Vahatra_poly.rds")
} else {

  Vahatra_poly <- AP_Vahatra %>%
    filter(st_geometry_type(.) == "MULTIPOLYGON") %>%
    st_cast("POLYGON")
  
  # Constitution d'un portefeuille (voir la documentation)
  Vahatra_poly <- init_portfolio(x = Vahatra_poly, 
                                 years = 2000:2020,
                                 outdir = "data/mapme_Vahatra",
                                 cores = 24,
                                 add_resources = TRUE,
                                 verbose = TRUE)
  
  # Données d'accessibilité de Nelson et al. (2018)
  Vahatra_poly <-  get_resources(x = Vahatra_poly, resource = "nelson_et_al",  
                                 range_traveltime = "5k_110mio")
  # Modèle numérique de terrain SRTM de la NASA
  Vahatra_poly <- get_resources(x = Vahatra_poly , resource = "nasa_srtm")
  
    # Indicateurs d'accessibilité
  Vahatra_poly <- calc_indicators(x = Vahatra_poly,
                                  "traveltime",  stats_accessibility = "mean",
                                  engine = "extract")
  # Indicateurs de relief de terrain
  Vahatra_poly <- calc_indicators(x = Vahatra_poly,
                                  indicators = c("tri", "elevation"),
                                  stats_tri = "mean", stats_elevation = "mean")
  
  #   # On récupère aussi les données de Global Forest Watch sur le couver forestier
  # Vahatra_poly <- get_resources(x = Vahatra_poly, 
  #                               resources = c("gfw_treecover", "gfw_lossyear"))
  #   # Indicateurs de couvert forestier
  # Vahatra_poly <- calc_indicators(x = Vahatra_poly,
  #                                 indicators = "treecover_area", 
  #                                 min_cover = 30, min_size = 1)
  
  save(Vahatra_poly, file = "data/Vahatra_poly.rds")
}

Mapme produit des colonnes imbriquées pour chaque observation, car dans bien des cas, on peut avoir plusieurs valeurs (par année) pour une même observation, voire plusieurs variables (par exemple, le calcul de l’indicateur traveltime produit des estimations de distance par rapport à une ville pour plusieurs tailles de ville possible. Lorsqu’on spécifie une taille, il produit deux variables : la distance estimée et la taille de la ville prise en compte pour l’estimation.

Cette imbrication n’est pas indispensable pour les trois variables calculées ici (indice de terrain accidenté, distance à une ville et altitude), car on ne cherche qu’une valeur par observation. On va donc dés-imbriquer les variables.

On va également procéder à une consolidation des données issues des traitements de {mapme.biodiversity}. Le jeu de données AP_Vahatra contenait 98 aires protégées, avec des géométries de types “multi-polygones”. Certaines de ces aires protégées étaient en effet composées de plusieurs polygones disjoints. Ces polygones disjoints ont été scindés pour être traités séparément dans le jeu de données AP_poly. Avant de repasser sur des analyses par aire protégée, on va agréger les statistiques d’AP_poly afin d’avoir pour chaque variable une valeur par aire protégée.

Code 
# Valeur agrégées par AP (moyennes pondérées par la surface)
Vahatra_vars_terrain <- Vahatra_poly %>%
  unnest(cols = c(tri, elevation, traveltime)) %>%
  st_drop_geometry() %>%
  select(nom, hectares, indice_accidente = tri_mean, dist_ville = minutes_mean, 
         altitude = elevation_mean) %>%
  group_by(nom) %>%
  summarise(indice_accidente = weighted.mean(indice_accidente, hectares,
                                             na.rm = TRUE),
            dist_ville = weighted.mean(dist_ville, hectares,
                                       na.rm = TRUE),
            altitude = weighted.mean(altitude, hectares,
                                     na.rm = TRUE))
# Valeurs qu'on insère dans le jeu de données de travail
load("data/ch1_AP_Vahatra.rds")
AP_Vahatra <- AP_Vahatra %>%
  left_join(Vahatra_vars_terrain, by = "nom")

save(AP_Vahatra, file = "data/ch2_AP_Vahatra.rds")

On doit aussi se rappeler que les aires protégées sont parfois composées de plusieurs polygones disjoints et que mapme.biodiversity a calculé chaque indicateur pour chaque polygone séparément. Pour chaque aire protégée, on va donc faire la moyenne de ces indicateurs, pondérée par la surface respective de chaque polygone.

Données d’accessibilité : attention car elles présentent un possible biais d’endogénéité. La construction de route au cours des dernières décennies peut être lié à l’établissement ou non d’aires protégées. L’inclusion d’une variable de contrôle qui peut être en partie affectée par notre variable de traitement (la conservation) est susceptible de problème. Il existe une carte de 2000 qui pourrait être mobilisée :

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