9 Changement de frontières

On charge les données de zones protégées historiques transmises par MNP et on les compare à la version de Novembre 2025 de WDPA.

Code

# We load the required libraries
library(tidyverse)
library(sf)
library(tmap)
library(wdpar)
library(gt)

# Chargement des données WDPA
wdpa_mdg <- wdpa_read("data/WDPA_WDOECM_Nov2025_Public_MDG.zip")

# The dataset seems to have some issues with missing .shx files
# We set an environment variable to restore them
Sys.setenv(SHAPE_RESTORE_SHX = "YES")

# Load the historical borders dataset
mnp_hist <- st_read(
  "data/no_id/allAP_ap2000_2010_Merge_Lab/allAP_ap2000_2010_Merge_Lab.shp",
  quiet = TRUE
) |>
  st_set_crs(29701) |>
  st_transform(st_crs(wdpa_mdg))

tm_shape(wdpa_mdg) +
  tm_polygons(fill = "blue", alpha = 0.5) +
  tm_shape(mnp_hist) +
  tm_polygons(fill = "red", alpha = 0.5)

On va maintenant retrouver, pour chaque polygone de

Code

# Extract names and WDPA IDs from the WDPA dataset
wdpa_ids <- wdpa_mdg |>
  st_drop_geometry() |>
  select(WDPAID, NAME)


# Match names
mnp_hist_ids <- mnp_hist |>
  st_drop_geometry() |>
  select(NAME_BASE) |>
  mutate(NAME = str_remove(NAME_BASE, " - 2000 2010")) |>
  left_join(wdpa_ids, by = "NAME")

Code

# Keep only what we need, make names consistent
wdpa <- wdpa_mdg %>%
  select(WDPAID, WDPA_NAME = NAME) %>%
  mutate(WDPA_NAME = str_squish(WDPA_NAME))

hist <- mnp_hist %>%
  mutate(
    HIST_ID = row_number(), # stable id for grouping
    HIST_NAME = str_remove(NAME_BASE, " - 2000 2010") %>% str_squish()
  ) %>%
  select(HIST_ID, HIST_NAME)

# Use a projected CRS
ea_crs <- 29702 # Equal area for Madagascar

wdpa_ea <- wdpa %>%
  st_make_valid() %>%
  st_transform(ea_crs) %>%
  mutate(wdpa_area = as.numeric(st_area(geometry)))

hist_ea <- hist %>%
  st_make_valid() %>%
  st_transform(ea_crs) %>%
  mutate(hist_area = as.numeric(st_area(geometry)))

# Build candidate pairs using spatial index, then compute intersections only for those
cand <- st_intersects(hist_ea, wdpa_ea)

pairs <- tibble(
  HIST_ID = rep(hist_ea$HIST_ID, lengths(cand)),
  wdpa_row = unlist(cand)
) %>%
  left_join(st_drop_geometry(hist_ea), by = "HIST_ID") %>%
  mutate(wdpa_row_id = wdpa_row) %>%
  left_join(
    wdpa_ea %>% st_drop_geometry() %>% mutate(wdpa_row_id = row_number()),
    by = "wdpa_row_id"
  )

# Compute overlap areas for each candidate pair
# (rowwise intersection on subsets keeps memory manageable)
overlaps <- pairs %>%
  rowwise() %>%
  mutate(
    inter_area = {
      g1 <- st_geometry(hist_ea)[HIST_ID]
      g2 <- st_geometry(wdpa_ea)[wdpa_row]
      as.numeric(st_area(st_intersection(g1, g2)))
    },
    pct_hist_within_wdpa = inter_area / hist_area,
    pct_wdpa_within_hist = inter_area / wdpa_area
  ) %>%
  ungroup() %>%
  filter(inter_area > 0)

# Selection rule:
# "max for both" implemented as lexicographic max on both percentages,
# tie-break by smallest absolute difference between the two
best_match <- overlaps %>%
  group_by(HIST_ID) %>%
  arrange(
    desc(pct_hist_within_wdpa),
    desc(pct_wdpa_within_hist),
    abs(pct_hist_within_wdpa - pct_wdpa_within_hist),
    desc(inter_area)
  ) %>%
  slice(1) %>%
  ungroup() %>%
  select(
    HIST_ID,
    HIST_NAME,
    WDPAID,
    WDPA_NAME,
    inter_area,
    pct_hist_within_wdpa,
    pct_wdpa_within_hist
  )

# Attach back to mnp_hist (keeps original geometry / CRS)
mnp_hist_matched <- mnp_hist %>%
  mutate(HIST_ID = row_number()) %>%
  left_join(best_match, by = "HIST_ID")

On voit qu’on a plusieurs polygones historiques transmis par MNP qui appartiennent à une même aire protégée dans la version WDPA de Novembre 2025. On va donc les regrouper dans une même entité, afin d’avoir, pour chaque aire protégée, d’un côté une entité rassemblant les polygones historiques et de l’autre côté le polygone WDPA actuel.

Code

mnp_hist_by_wdpaid <- mnp_hist_matched %>%
  filter(!is.na(WDPAID)) %>%
  st_make_valid() %>%
  mutate(geometry = st_buffer(geometry, 0)) %>% # often fixes duplicate vertices
  group_by(WDPAID, WDPA_NAME) %>%
  summarise(
    n_parts = n(),
    geometry = st_union(geometry),
    .groups = "drop"
  ) %>%
  st_collection_extract("POLYGON") %>%
  st_cast("MULTIPOLYGON")

Maintenant on va comparer leur similarité et repéré les 15 plus différents.

9.1 Revue des principales modifications

Code

# Equal area CRS for area ratios
ea_crs <- 29702 # Laborde

# WDPA as multipolygon by WDPAID (recommended)
wdpa_by_id <- wdpa_mdg %>%
  select(WDPAID, WDPA_NAME = NAME) %>%
  st_make_valid() %>%
  group_by(WDPAID, WDPA_NAME) %>%
  summarise(geometry = st_union(geometry), .groups = "drop") %>%
  st_collection_extract("POLYGON") %>%
  st_cast("MULTIPOLYGON")

# Transform + areas
hist_ea <- mnp_hist_by_wdpaid %>%
  st_make_valid() %>%
  st_transform(ea_crs) %>%
  mutate(hist_area = as.numeric(st_area(geometry)), hist_geom = geometry) %>%
  st_drop_geometry() %>%
  rename(WDPA_NAME_hist = WDPA_NAME)

wdpa_ea <- wdpa_by_id %>%
  st_make_valid() %>%
  st_transform(ea_crs) %>%
  mutate(wdpa_area = as.numeric(st_area(geometry)), wdpa_geom = geometry) %>%
  st_drop_geometry()

comp <- hist_ea %>%
  left_join(wdpa_ea, by = "WDPAID") %>%
  rowwise() %>%
  mutate(
    inter_area = if (is.na(wdpa_area)) {
      0
    } else {
      as.numeric(sum(st_area(st_intersection(hist_geom, wdpa_geom))))
    },
    pct_hist_within_wdpa = ifelse(
      hist_area > 0,
      inter_area / hist_area,
      NA_real_
    ),
    pct_wdpa_within_hist = ifelse(
      !is.na(wdpa_area) & wdpa_area > 0,
      inter_area / wdpa_area,
      NA_real_
    ),
    similarity = pmin(pct_hist_within_wdpa, pct_wdpa_within_hist, na.rm = TRUE),
    dissimilarity = 1 - similarity
  ) %>%
  ungroup() %>%
  select(
    WDPAID,
    WDPA_NAME_hist,
    WDPA_NAME,
    pct_hist_within_wdpa,
    pct_wdpa_within_hist,
    similarity,
    dissimilarity
  ) %>%
  arrange(desc(dissimilarity))

comp %>%
  slice_head(n = 15) %>%
  mutate(
    pct_hist_within_wdpa = scales::percent(
      pct_hist_within_wdpa,
      accuracy = 0.1
    ),
    pct_wdpa_within_hist = scales::percent(
      pct_wdpa_within_hist,
      accuracy = 0.1
    ),
    similarity = round(similarity, 3)
  ) %>%
  select(
    WDPA_NAME,
    pct_hist_within_wdpa,
    pct_wdpa_within_hist,
    similarity
  ) %>%
  gt() %>%
  tab_header(
    title = "Comparaison des limites historiques et WDPA",
    subtitle = "15 aires protégées les plus dissemblables"
  ) %>%
  cols_label(
    WDPA_NAME = "Aire protégée",
    pct_hist_within_wdpa = "Part du périmètre historique incluse dans WDPA",
    pct_wdpa_within_hist = "Part de WDPA incluse dans le périmètre historique",
    similarity = "Indice de similarité"
  ) %>%
  fmt_number(columns = similarity, decimals = 3) %>%
  data_color(
    columns = similarity,
    domain = c(0, 1),
    palette = c("#d73027", "#fee08b", "#1a9850")
  ) %>%
  opt_table_outline()

Aire protégée	Part du périmètre historique incluse dans WDPA	Part de WDPA incluse dans le périmètre historique	Indice de similarité
Comparaison des limites historiques et WDPA
15 aires protégées les plus dissemblables
Baie d'Antogil	100.0%	0.3%	0.003
Complexe des AP Ambohimirahavavy Marivorahona	100.0%	8.9%	0.089
Beza Mahafaly	99.5%	11.4%	0.114
Tsimanampesotse	99.4%	23.2%	0.232
Zones Humides de Sahamalaza	95.4%	46.9%	0.469
Kirindy Mite	96.3%	48.9%	0.489
Cap Sainte Marie	98.4%	61.4%	0.614
Anjanaharibe_sud	97.9%	63.5%	0.635
Manongarivo	99.1%	80.5%	0.805
Ambatovaky	82.6%	80.9%	0.809
Befotaka Midongy	88.7%	100.0%	0.887
Mananara Nord	90.7%	94.8%	0.907
Montagne d'Ambre	97.1%	90.8%	0.908
Masoala	92.9%	95.2%	0.929
Ranomafana	93.4%	100.0%	0.934

Maintenant, on visualise sous forme de carte simple.

Code

plot_comp_facets <- function(comp, rows, hist_sf, wdpa_sf, ncol = 1) {
  sub <- dplyr::slice(comp, rows) %>%
    dplyr::transmute(WDPAID, panel_name = WDPA_NAME) %>%
    dplyr::filter(!is.na(WDPAID))

  hist <- hist_sf %>%
    dplyr::semi_join(sub, by = "WDPAID") %>%
    dplyr::select(WDPAID, geometry) %>%
    dplyr::mutate(layer = "Historical")

  wdpa <- wdpa_sf %>%
    dplyr::semi_join(sub, by = "WDPAID") %>%
    dplyr::select(WDPAID, geometry) %>%
    dplyr::mutate(layer = "WDPA")

  dat <- dplyr::bind_rows(wdpa, hist) %>%
    dplyr::left_join(sub, by = "WDPAID") %>%
    dplyr::mutate(
      WDPAID = factor(WDPAID, levels = sub$WDPAID),
      layer = factor(layer, levels = c("WDPA", "Historical"))
    )

  tmap::tmap_mode("plot")

  tmap::tm_shape(dat) +
    tmap::tm_polygons(
      fill = "layer",
      fill.scale = tmap::tm_scale(
        values = c(WDPA = "blue", Historical = "red")
      ),
      fill_alpha = 0.5,
      col = NA
    ) +
    tmap::tm_facets(by = "WDPAID", ncol = ncol) +
    tmap::tm_layout(
      legend.outside = TRUE,
      panel.labels = sub$panel_name
    )
}

# example
plot_comp_facets(
  comp,
  rows = 1:15,
  hist_sf = mnp_hist_by_wdpaid,
  wdpa_sf = wdpa_by_id,
  ncol = 3
)

Code

legal_texts.rds <- read_rds("data/id/legal_texts.rds")

9.2 Questions en suspens

Validation de ces différences
Quelles dates d’entrée en vigueur des modifications de limites ?