Une approche courante consiste à diviser le territoires en mailles, carrées ou en forme d’alvéoles d’abeilles (hexagones), et à calculer des indicateurs pour chacune de ces mailles.
4.1 Constitution d’un maillage
On montre ci-dessous comment cette approche fonctionne. La première étape consiste à dessiner un carré autour des aires protégées malgaches, puis à subdiviser ce grand carré en un damier de formes hexagonales. Enfin, on ne garde que les hexagones qui se trouvent dans les frontières terrestres de Madagascar.
Code
library(tidyverse)library(tmap)library(sf)library(mapme.biodiversity)library(tidygeocoder) # pour obtenir les coordo GPS d'un point à partir de son nomlibrary(maptiles) # Pour télécharger des fonds de carte# Surface des hexagones en km2taille_hex <-5# 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/grille_mada_donnees_raster.rds")) {load("data/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()}
Le maillage produit est trop fin pour être visible à l’échelle du pays, mais on peut l’observer en zoomant sur une zone spécifique.
Code
sf_use_s2(TRUE) # On compte le nombre d'hexagonesn_hex <-nrow(grille_mada)# Carte pour visualiser le résultat ------------------------------------------# A supprimer après validation, l'échappement ne fonctionne pas. # if (file.exists("data/carte_zoom.rds")) {# load("data/carte_zoom.rds")# load("data/elements_carte_zoom.rds")# } else {## Carte de droite : zoom sur une zone spécifique-----------------------------# On part d'un dataframe contenant une adressenom_centre_zoom <-"Maroantsetra"zoom_centre <-data.frame(address = nom_centre_zoom) %>%geocode(address, method ="osm") %>%# on retrouve sa localisation xyselect(long, lat) %>%# on ne garde que le xyas.numeric() %>%# qu'on passe en format numérique attendu par st_pointst_point() %>%# On le spécifie en pointst_sfc(crs ="EPSG:4326") # On crée une boîte de 100km zoom_boite <- zoom_centre %>%# On repart du centrest_buffer(dist =50000) %>%# On crée un cercle de 50km de rayonst_make_grid(n =1) # On filtre les alvéoles pour ne garder que celles qui sont dans le zoomgrille_zoom <-st_intersection(grille_mada, zoom_boite)# On télécharge un fond de carte pour la carte de droitefond_carte_zoom <-get_tiles(zoom_boite, provider ="Stamen.Terrain", zoom =10, crop =TRUE) save(n_hex, nom_centre_zoom, zoom_centre, zoom_boite, grille_zoom, fond_carte_zoom, file ="data/elements_carte_zoom.rds")# On était restés en mode interactif, on repasse en mode statique pour les # cartestmap_mode("plot")# On génère la carte de droitecarte_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 =1) +tm_credits(get_credit("Stamen.Toner"),bg.color ="white",align ="right",position =c("right", "BOTTOM"))save(carte_zoom, zoom_boite, file ="data/carte_zoom.rds")# }## Carte de gauche : simple à réaliser mais hexagones non visibles -------------load("data/contour_mada.rds")carte_grille <-tm_shape(contour_mada) +tm_polygons(col ="grey") +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, "km2"),main.title.size =1)# Assemblage des deux cartes ---------------------------------------------------tmap_arrange(carte_grille, carte_zoom, ncol =2)
4.2 Récupération des données pour le maillage
On va ensuite utiliser le package mapme.biodiversity pour calculer, pour chaque hexagones, une série d’indicateurs : temps de parcours jusqu’à la ville (définie comme toute localité de 5000 habitants) la plus proche en 2015, teneur du sol en argile et couvert forestier par année).
Code
if (file.exists("data/grille_mada_donnees_raster.rds")) {load("data/grille_mada_donnees_raster.rds")} else {# Traitement des données satellitaires avec {mapme.bidiversity}---------------# Constitution d'un portefeuille (voir la documentation) grille_mada <-init_portfolio(x = grille_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 =30, 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ésultatsave(grille_mada, file ="data_s3/grille_mada_donnees_raster.rds")}# Assemblage des deux cartes ---------------------------------------------------tmap_arrange(carte_grille, carte_zoom, ncol =2)
On peut représenter sous forme de cartes et d’histogrammes les différentes valeurs des indicateurs générés à partir des données satellitaires.
Les cartes et histogrammes ci-dessus illustrent la distribution des variables spatiales calculées par hexagones.
4.3 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 satellitaires déjà calculés pour ces hexagones.
Code
if (file.exists("data/grille_mada_AP.rds")) {load("data/grille_mada_AP.rds")} else {# On charge les données d'aires protégées élaborées au chapitre précédentload("data/ch3_AP_Vahatra.rds")# On prépare ces données pour les joindre avec celles en mailles aires_prot_mada <- AP_Vahatra %>%st_make_valid() %>%# Corrige les erreurs topo dans certains polygonesmutate(AP_ligne =row_number()) %>%# Intègre le numéro de ligne dans un champmutate(an_creation =year(date_creation)) # passe les dates en années# 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 spatialeselect(AP_ligne, nom)# On sélectionne des infos additionnelles qu'on va inclure dans les données info_vahatra_a_inclure <- aires_prot_mada %>%st_drop_geometry() %>%# Enlève l'information spatialeselect(AP_ligne, nom, an_creation, cat_iucn = cat__iucn, gestionnaire = gest_2) # On ne garde que les variables d'intérêt# 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_AP <- grille_mada %>%mutate(AP_ligne =st_intersects(., aires_prot_mada), # liste des n° de lignes d'AP qui recoupentAP_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ériqueleft_join(aires_prot_mada_rang_id, by ="AP_ligne") %>%# récupère l'id de l'APrename(AP_touche = nom) %>%# on renomme pour différentiermutate(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(AP_inclus = nom) %>%select(-AP_ligne) grille_mada_AP <- grille_mada_AP %>%st_sf() %>%mutate(position_ap =ifelse(is.na(AP_touche), "Extérieur",ifelse(!is.na(AP_inclus), "Intérieur","Frontière")),ref_AP =ifelse(position_ap =="Intérieur", AP_inclus, ifelse(position_ap =="Frontière", AP_touche, NA))) %>%left_join(info_vahatra_a_inclure, by =c("ref_AP"="nom")) %>%relocate(geometry, .after =last_col()) grille_mada_AP <- grille_mada_AP %>%# On met à plat les données de distanceunnest(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)save(grille_mada_AP, file ="data/grille_mada_AP.rds")}# Une vue après classificationtm_shape(grille_mada_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 margemain.title.position =c("center", "top"),main.title.size =1,legend.position =c("left", "top"),legend.outside =FALSE)
