Protected Areas data sources

We have four different data sources showing various boundaries and categories of protected areas. The question is: which one is correct?

Data sources :

MINEF-ICF Map of Forêts classées 2019

# LOADING DATA

# CIV departments ----------------------------
civ_departments <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/DEPARTEMENT/OUT/CIV_DEPARTEMENTS.geojson",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) %>%
  dplyr::select(LVL_4_CODE, LVL_4_NAME, LVL_3_CODE, LVL_3_NAME, geometry) %>% 
  st_transform(4326) %>%  # the cooperative dots do not appear on the map when using the same crs as civ departments (EPSG 32630) - I don't know why. So I rather change the crs of civ departments into EPSG 4326 (WGS84)
  # Correcting an invalid polygon (self-intersection)
  as("Spatial") %>% 
  rgeos::gBuffer(byid=TRUE, width=0) %>% 
  st_as_sf()


# All categories WDPA ---------------------------------
civ_protected_areas_wdpa <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/WDPA/OUT/WDPA_PROTECTED_AREAS_CIV_AUG21.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) %>% 
  st_transform(4326) # making sure of the CRS 

civ_protected_areas_wdpa$IUCN_CAT[civ_protected_areas_wdpa$IUCN_CAT == "Not Reported"] <- NA
civ_protected_areas_wdpa$IUCN_CAT[civ_protected_areas_wdpa$IUCN_CAT == "Not Applicable"] <- NA


# Forêts classées SST GEOPORTAIL ---------------------------------
civ_forets_classees_sst <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/SST_GEOPORTAIL/OUT/SST_FORETS_CLASSEES_CIV_AUG21.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T))


# Forêts classées, Parcs and Réserves BNETD ---------------------------------
civ_forets_classees_bnetd <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/BNETD/OUT/BNETD_FORETS_CLASSEES_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) 

civ_parcs_reserves_bnetd <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/BNETD/OUT/BNETD_PARCS_RESERVES_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) 

civ_protected_areas_bnetd <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/BNETD/OUT/BNETD_PROTECTED_AREAS_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) 


# Parcs & Reserves : Source OIPR through EFI ---------------------------------
civ_protected_areas_oipr <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/OIPR/OUT/OIPR_PROTECTED_AREAS_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) 
# PLOTTING MAPS

# Plot protected areas WDPA and Forêts classées SST
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_wdpa) + 
  tm_polygons(border.col = NULL, col = "gray20") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_forets_classees_sst) +
  tm_polygons(border.col = NULL, col = "coral1", alpha = 0.6) +
  tm_add_legend(labels = c("All protected areas (WDPA)", "Forets classees (SST)"), col = c("gray20", "coral1"), border.col = NULL)

# Plot protected areas BNETD
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_wdpa) + 
  tm_polygons(border.col = NULL, col = "gray20") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_forets_classees_bnetd) +
  tm_polygons(border.col = NULL, col = "steelblue", alpha = 0.6) +
tm_shape(civ_parcs_reserves_bnetd) +
  tm_polygons(border.col = NULL, col = "forestgreen", alpha = 0.6) +
  tm_add_legend(labels = c("All protected areas (WDPA)", "Forets classees (BNETD)", "Parcs & Reserves (BNETD)"), col = c("gray20", "steelblue", "forestgreen"), border.col = NULL) 

# Plot Forêts classées SST and BNETD
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_forets_classees_bnetd) + 
  tm_polygons(border.col = NULL, col = "steelblue") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_forets_classees_sst) +
  tm_polygons(border.col = NULL, col = "coral1", alpha = 0.6) +
  tm_add_legend(labels = c("Forets classees (BNETD)", "Forets classees (SST)"), col = c("steelblue", "coral1"), border.col = NULL)

# Plot protected areas WDPA & reported IUCN Categories
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_wdpa) + 
  tm_polygons(border.col = NULL, col = "IUCN_CAT", palette = c("green", "forestgreen", "gold2"), title = "IUCN Category (WDPA)") +
  tm_text(text = "IUCN_CAT", size = 0.5) + 
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)

# Plot protected areas BNETD & reported Categories
tm_shape(civ_departments) + 
  tm_polygons(col = "gray95") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_bnetd) + 
  tm_polygons(border.col = NULL, col = "CATEGORY", title = "Protected area category (BNETD)", palette = c("steelblue", "yellow", "forestgreen", "gold", "gold2", "gold3", "goldenrod4")) +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)

# Plot protected areas OIPR & reported Categories
tm_shape(civ_departments) + 
  tm_polygons(col = "gray95") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_oipr) + 
  tm_polygons(border.col = NULL, col = "CATEGORY", title = "Protected areas (OIPR)", palette = c("forestgreen", "green", "gold2")) +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)

# CIV WDPA Protected areas IUCN categories
table(civ_protected_areas_wdpa$IUCN_CAT)
## 
## Ia II IV 
##  3  9  3
#Number of forêts classées (WDPA)
sum(is.na(civ_protected_areas_wdpa$IUCN_CAT))
## [1] 228
# CIV BNETD Protected areas categories
table(civ_protected_areas_bnetd$CATEGORY)
## 
##                Forêt classée               Parc Animalier 
##                          222                            1 
##                Parc National            Réserve Botanique 
##                            8                            1 
##             Réserve de Faune Réserve de Faune et de Flore 
##                            1                            1 
##            Réserve Naturelle 
##                            2
# CIV OIPR Protected areas categories
table(civ_protected_areas_oipr$CATEGORY)
## 
##     Parc National Reserve integrale Reserve naturelle 
##                 9                 1                 7

Conclusion from the comparison of the maps

For Parcs and reserves: We can assume that OIPR - office in charge of Protected areas and Reserves in CIV - has the most up to date map. The shapefile we have has the same boundaries and categories as the image from their website below (WDPA has similar results).

16 Parcs (green) and Réserves (orange), image from the OIPR website

For Forêts Classées: From visual comparison, the SST map of forêts classées is the same as the MINEF pdf. However, the pdf is dated 2019 and does not include the latest change - i.e. does not yet list Mabi-yaya as a Reserve - Decret October 2019 - but the OIPR data does. WDPA does not have Mabi-Yaya as a reserve, it seems to be the map differing the most from the MINEF pdf map. BNETD is in between, having Mabi-Yaya as a Reserve but other Forêts classées are mapped in BNETD map and not in the MINEF pdf.

With the creation of the Mabi-Yaya reserve, the previous Mabi and Yaya Forêts classées did not completely disappeared but have been reshaped - Decret October 2019 on the new boundaries of the Mabi and Yaya Forêts classées. It seems that the BNETD map includes these new boundaries for those two Forêts classées, but we don’t know their category.

mabi_yaya_fc <- filter(civ_protected_areas_bnetd, NAME %in% "Mabi" | NAME %in% "Yaya") %>% 
  dplyr::select(-IUCNCAT, -YEAR) %>% 
  #Remove dimension Z (creates error otherwise)
  st_zm(drop = TRUE, what = "ZM") %>% 
  # MISTAKE in the BNETD dataset, Mabi and Yaya are reversed 
  # (check for comparison: https://www.researchgate.net/figure/Situation-geographique-de-la-zone-detude-et-de-trois-regions-du-Sud-Est-de-la-Cote_fig1_350592790)
  mutate(NAME = case_when(
    NAME == "Yaya" ~ "Forêt classée de MABI", 
    NAME == "Mabi" ~ "Forêt classée de YAYA")) %>% 
  mutate(SUP_HA = drop_units(units::set_units(st_area(geom), value = hectare)))

# Plot Mabi & Yaya Forêts classées according to BNETD
tmap_mode("view")
tm_shape(mabi_yaya_fc ) + 
  tm_polygons(border.col = NULL, title = "Mabi & Yaya Forêts classées (BNETD)", col = "forestgreen") +
  tm_text(text = "NAME", col = "black", size = 0.8) + 
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_departments) + 
  tm_polygons(col = "gray97", alpha = 0.2) + 
  tmap_options(check.and.fix = TRUE) 

Conclusion: we decide to work with the merger of forêts classées from the SST layer, and parcs and reserves from the OIPR layer to get a single file with all protected areas in CIV. For Mabi and Yaya Forêts classées : we remove them from SST (old boundaries) and replace them with the new versions from the BNETD map.

Confirmation should be asked to officials in CIV, specially to confirm if SST forêts classées layer is better than BNETD.

# Merging SST forets classees and OIPR parcs and reserves
civ_forets_classees_sst_CLEAN <- civ_forets_classees_sst %>% 
  transmute(NAME = nom, 
            CATEGORY = case_when(
              cate %in% c(1:4) ~ paste("Forêt classée de catégorie", civ_forets_classees_sst$cate, sep=" "),
              cate == "Forêt classée de savane" ~ "Forêt classée de savane",
              TRUE ~ as.character(NA)),
            SUP_HA = aire) %>% 
  # Remove dimension Z (creates error otherwise)
  st_zm(drop = TRUE, what = "ZM") %>% 
  # Remove old versions of Forêts classées Mabi and Yaya 
  dplyr::filter(NAME != "Forêt classée de MABI" & NAME != "Forêt classée de YAYA") %>% 
  # Add new versions of Mabi and Yaya (from BNETD)
  rbind(mabi_yaya_fc) %>% 
  # Add the area for all 
  mutate(SUP_HA = drop_units(units::set_units(st_area(geom), value = hectare)))

civ_protected_areas_oipr_sst <- rbind(dplyr::select(civ_protected_areas_oipr, -GID), civ_forets_classees_sst_CLEAN) #SUP_HA are coming from the root oipr and sst files and were not computed here

# Plot merged file: parc & reserves OIPR + forêts classées SST 
tmap_mode("view")
tm_shape(civ_departments) + 
  tm_polygons(col = "gray95") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_oipr_sst) + 
  tm_polygons(border.col = NULL, col = "CATEGORY", title = "Protected areas OIPR & SST", palette = c("coral", "coral1", "coral2","coral3", "coral4","forestgreen", "green", "gold2")) +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)
# Plot merged file: parc & reserves OIPR + forêts classées SST 
tmap_mode("plot")
tm_shape(civ_departments) + 
  tm_polygons(col = "gray95") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_oipr_sst) + 
  tm_polygons(border.col = NULL, col = "CATEGORY", title = "Protected areas OIPR & SST", palette = c("coral", "coral1", "coral2","coral3", "coral4","forestgreen", "green", "gold2")) +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)

Save merged OIPR and SST layer in S3

s3write_using(civ_protected_areas_oipr_sst, FUN = write_sf,
               bucket = "trase-storage",
               object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/OIPR_SST/OUT/OIPR_SST_PROTECTED_AREAS_CIV.gpkg")

# Load
civ_protected_areas_oipr_sst <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/OIPR_SST/OUT/OIPR_SST_PROTECTED_AREAS_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T))