We have four different data sources showing various boundaries and categories of protected areas. The question is: which one is correct?
Data sources :
WDPA (World Database of Protected Areas): all categories (parcs, reserves, forêts classées): downloaded on the wdpa website in August 2021. This map identifies more areas as protected compared to BNETD (below). Most protected areas do not have a reported IUCN category - cf. map and table below.
SST (Surveillance Spatiale des Terres): downloaded on the Geoportail in August 2021. Includes Forêts classées only (I could not find a shp for parcs & reserves on this website). It does not match the boundaries of the WDPA protected areas.
BNETD (Bureau national d’études techniques et de développement): sent by Nitidae in August 2021, but we don’t have the date of creation of this layer. All categories are gathered in two shapefiles, one with forêts classées and one with Parcs & Réserves. Variations can be observed, with areas indicated as protected in WDPA shp that are not present in BNETD. As compared to SST, some areas are indicated as forêts classées by SST and not by BNETD and vice-versa.
OIPR (Office Ivoirien des Parcs et Réserves): sent by EFI in August 2021. It is supposed to be the latest version of national parks and reserves boundaries. EFI did not have the shp for forêts classées, only a pdf map from the MINEF (Ministère des Eaux et Forêts), cf. below.
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
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)
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))