3  Spatial covariates

Ce document vise à charger les données géophysiques des ménages à l’aide du package R mapme.biodiversity (https://mapme-initiative.github.io/mapme.biodiversity/index.html) Görgen and Bhandari (2022). Ce package va automatiser la récupération de données volumineuses au format raster et leur traitement pour produire une série d’indicateurs appliqués à des polygones.

3.1 Données DHS

Avant le chargement des données des clusters GPS, nous définissons les systèmes de coordonnées. Les coordonnées GPS sont fournies dans le système de référence géographique WGS 84 (EPSG: 4326). Toutefois, pour effectuer des calculs précis de distances ou de surfaces, nous utilisons le système de projection officiel de Madagascar, qui est le Laborde (EPSG: 29702). On travaille en 4326 puis on le transforme en 29702 lorsque les opérations métriques sont nécessaires. Après le chargement des données, certains clusters n’ont pas de coordonnées valides (0,0). On les retire pour éviter des erreurs dans les calculs géospatiaux. On regroupe ensuite les clusters de l’année 1997, 2008, 2011, 2013, 2016 et 2021 et 2021 dans un seul objet spatial pour faciliter les opérations ultérieures. Puis, on applique un buffer zone de 10 km autour de chaque cluster, qui serviront à extraire les données géophysiques et à analyser les expositions à certaines zones.

Code 
# if (packageVersion("mapme.biodiversity") < "0.9.4") {
#   remotes::install_github("mapme-initiative/mapme.biodiversity")
# }


library(tidyverse) #Manipulation et visualisation des données
library(sf) #Analyse des données spatiales
library(tmap) #Analyse cartographique
library(labelled) # Manipulation des labels
library(mapme.biodiversity)
library(progressr) # Pour avoir des barres de progression
library(tictoc) # Pour minuter le temps d'exécution
library(future) # Pour permettre du calcul parallèle


# Systèmes de coordonnées de référence 
standard_crs <- 4326
mdg_crs <- 29702 

# On charge les données
gps_1997_initial <- st_read("data/raw/dhs/DHS_1997/MDGE32FL/MDGE32FL.shp", quiet = TRUE)
gps_2008_initial <- st_read("data/raw/dhs/DHS_2008/MDGE53FL/MDGE53FL.shp", quiet = TRUE) 
gps_2011_initial <- st_read("data/raw/dhs/DHS_2011/MDGE61FL/MDGE61FL.shp", quiet = TRUE) 
gps_2013_initial <- st_read("data/raw/dhs/DHS_2013/MDGE6AFL/MDGE6AFL.shp", quiet = TRUE) 
gps_2016_initial <- st_read("data/raw/dhs/DHS_2016/MDGE71FL/MDGE71FL.shp", quiet = TRUE)
gps_2021_initial <- st_read("data/raw/dhs/DHS_2021/MDGE81FL/MDGE81FL.shp", quiet = TRUE)

# Fonction qui vérifie que les coordonnées ne sont pas nulles
check_coordinates <- function(dhs_gps, country_polygon, negate = FALSE) {
  dhs_gps %>%
    filter(LONGNUM != 0 | LATNUM != 0)
}

gps_1997 <- check_coordinates(gps_1997_initial, contour_mada)
gps_2008 <- check_coordinates(gps_2008_initial, contour_mada)
gps_2011 <- check_coordinates(gps_2011_initial, contour_mada)
gps_2013 <- check_coordinates(gps_2013_initial, contour_mada)
gps_2016 <- check_coordinates(gps_2016_initial, contour_mada)
gps_2021 <- check_coordinates(gps_2021_initial, contour_mada)

gps_all <- bind_rows(gps_1997, gps_2008, gps_2011, gps_2013,  gps_2016, 
          gps_2021)

# On crée les buffers
buffer_all <- gps_all %>%
  st_transform(mdg_crs) %>% # projection dans le CRS malgache
  st_buffer(dist = 10000) %>% # 10 km
  st_transform(standard_crs) # Retour en WGS 84

# tmap_mode("view")
# tm_shape(buffer_all) +
#   tm_polygons() +
#   tm_facets("DHSYEAR")

4 Chargement des données géophysiques des ménages

Nous allons charger et extraire les données sur les données géophysiques des ménages qui serviront de variables d’appariement et de variables de contrôle pour l’analyse. Ces covariables sont des covariables succeptibles d’influencer le traitement (aires protégées créées à proximité du ménage) et le résultat (wealth index et inégalité des ménages):

  • Le couvert forestier en 2000 Hansen et al. (2013)

  • La pente et l’altitude NASA JPL (2020)

  • La densité de population en 2000 WorldPop and CIESIN (2018)

  • L’accessibilité en 2000 Uchida and Nelson (2011)

  • Les données worldclim Vicente-Serrano, Beguería, and López-Moreno (2010)

Afin d’automatiser le téléchargement et la gestion de ces données spatiales, nous configurons notre environnement de travail à l’aide de mapme.biodiversity. Nous utilisons la fonction get_resources () pour cela en démarrant un chronomètre avec le package tictoc pour suivre le temps d’exécution du bloc de code.

Code 
# Définir le chemin relatif pour ton répertoire local
outdir <- "data/raw/mapme"
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
mapme_options(outdir = outdir, verbose = TRUE)


# Couvert forestier
tic()
with_progress({
  buffer_all <- buffer_all  %>% 
    get_resources(get_gfw_treecover())
})
toc() # 620,11 sec elapsed

# Perte de couvert
tic()
with_progress({
  buffer_all <- buffer_all %>% 
  get_resources(get_gfw_lossyear()) 
})
toc() # 266,57 sec elapsed

# NASA SRTM
tic()
with_progress({
  buffer_all <- buffer_all %>% 
  get_resources(get_nasa_srtm()) 
})
toc() # 4578,35 sec elapsed

# Worldpop 2000
options(timeout = 7200)  

tic()
with_progress({
  buffer_all <- buffer_all %>%
    get_resources(get_worldpop(years = 2000))
})
toc()
  
# Accesibility
tic()
with_progress({
  buffer_all <- buffer_all %>% 
  get_resources(get_accessibility_2000()) 
})
toc()

# Maximum temperatures
tic()
with_progress({
  buffer_all <- buffer_all %>% get_resources(
    get_worldclim_max_temperature(years = 1980:2021, resolution = "2.5m")
  )
})
toc() # 17746 sec elapsed

# Minimum temperatures
tic()
with_progress({
  buffer_all <- buffer_all %>% get_resources(
    get_worldclim_min_temperature(years = 1980:2021, resolution = "2.5m")
  )
})
toc() 

# Precipitations
tic()
with_progress({
  buffer_all <- buffer_all %>% get_resources(
    get_worldclim_precipitation(years = 1980:2021, resolution = "2.5m")
  )
})
toc() # 17 855,95 sec elapsed

5 Calcul des indicateurs

Nous utilisons la fonction calc_indicators() du package pour calculer les indicateurs des variables environnementales dans un rayon de 10 km autour de chaque grappe d’enquête.

Code 
if (!file.exists("data/derived/spatial_covars_staggered.rds")) {
  cat("Fichier introuvable, début du traitement...\n")
  
  # Créer un plan pour paralléliser les calculs
  plan(sequential)
  plan(multisession, workers = 4)
  
  
  
  # Maximum temperatures
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_temperature_max_wc(engine = "extract", stats = "mean")
      )
  })
  toc() # 20 391,84 sec elapsed
  
  # Minimum temperatures
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
         calc_temperature_min_wc( engine = "extract", stats = "mean")
      )
  })
  toc() # 28 357,5 sec elapsed 
  
  
  # Precipitations
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_precipitation_wc(engine = "extract", stats = "mean")
      )
  })
  toc() # 19636.08 sec elapsed 
  
  
  # Forest cover rate in 2000
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_treecover_area(years = 2000, min_size = 1, min_cover = 10)
      )
  })
  toc()   # 2602.96 sec elapsed
  
  # slope 
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_slope(engine = "extract", stats = "mean")
      )
  })
  toc() # 2398.6 sec elapsed    
  
  # Elevation
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_elevation(engine = "extract", stats = "mean")
      )
  })
  toc() # 4396,93 sec elapsed  
  
  # Population density
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_population_count(engine = "extract", stats = "mean")
      )
  })
  toc() # 5233,48 sec elapsed
  
  # Accessibility
  tic()
  with_progress({
    buffer_all <- buffer_all %>%
      calc_indicators(
        calc_traveltime_2000(engine = "extract", stats = "mean")
      )
  })
  toc() # 87.63 sec | 2594.06 sec elapsed   
  
  # Enregistrement des données 
  write_rds(buffer_all, "data/derived/spatial_covars_staggered.rds")
  cat("Données enregistrées dans spatial_covars_staggered.rds")
}

Nous obtenons un dataframe contenant le couvert forestier en 2000, la pente et l’altitude, la densité de population, l’accessibilité, la température maximale et minimale, ainsi que la précipitation.

Görgen, Darius A., and Om Prakash Bhandari. 2022. “Efficient Monitoring of Global Biodiversity Portfolios,” April. https://doi.org/10.32614/CRAN.package.mapme.biodiversity.
Hansen, Matthew C., Peter V. Potapov, Rebecca Moore, Matt Hancher, S. A. Turubanova, Alexandra Tyukavina, David Thau, et al. 2013. “High-Resolution Global Maps of 21st-Century Forest Cover Change.” Science 342 (6160): 850853. http://www.sciencemag.org/content/342/6160/850.short.
NASA JPL. 2020. “NASADEM Merged DEM Global 1 Arc Second V001.” https://doi.org/10.5067/MEASURES/NASADEM/NASADEM_HGT.001.
Uchida, Hirotsugu, and Andrew Nelson. 2011. “Agglomeration Index: Towards a New Measure of Urban Concentration.” Urbanization and Development: Multidisciplinary Perspectives, January. https://doi.org/10.1093/acprof:oso/9780199590148.003.0003.
Vicente-Serrano, Sergio M., Santiago Beguería, and Juan I. López-Moreno. 2010. “A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index.” Journal of Climate 23 (7): 16961718. https://doi.org/10.1175/2009JCLI2909.1.
WorldPop, and CIESIN. 2018. “Global 1km Population.” https://doi.org/10.5258/SOTON/WP00647.