Introducción
Para quien trabaja con indicadores de seguimiento de los ODS, el ODS 15 representa un reto cuando hay de bajar al detalle municipal. Incluso la cobertura de bosque puede presentar dudas y dificultades, como ya vimos en una nota anterior.
Un indicador muy significativo recomendado por los informes de REDS es el Porcentaje de espacio natural con algún tipo de protección frente al total de la superficie municipal. Según sus informes (ver por ejemplo, Los Objetivos de Desarrollo Sostenible en 100 ciudades españolas), la fuente de datos son el SIOSE y el catastro.
Nuestra propuesta alternativa es obtener el mismo indicador utilizando una fuente de datos global sobre areas protegidas: la base de datos cartográficos de protectedplanet.net. Protected Planet es la fuente de datos más actualizada y completa sobre áreas protegidas, actualizada mensualmente con aportaciones de gobiernos, organizaciones no gubernamentales, propietarios de tierras y comunidades.
Si hay que elaborar indicadores para diversas zonas de del mundo o incluso de España, nos podemos encontrar falta de datos o con datos contradictorios entre diferentes fuentes. Para ello hay que encontrar fuentes homogéneas que cubran de manera igual todas las zonas que queremos analizar.
Protected Planet es un proyecto conjunto entre el Programa de las Naciones Unidas para el Medio Ambiente y la Unión Internacional para la Conservación de la Naturaleza (UICN), y está gestionado por el Centro de Monitoreo de la Conservación Mundial del Programa de las Naciones Unidas para el Medio Ambiente (PNUMA-WCMC), en colaboración con gobiernos, organizaciones no gubernamentales, instituciones académicas y industria. La WDPA se actualiza mensualmente.
Tal y como se podrá comprobar más adelante, la representación generada por los datos de dicha fuente coincide con la del MITECO.
Los datos
Los datos se obtienen en la página protectedplanet.net/en/thematic-areas/wdpa?tab=WDPA.
También se dispone de una biblioteca de R, wdpar, para
poderlos descargar directamente desde nuestro codigo, que se obtiene en
CRAN:
Aquí hay un breve manual de uso de la biblioteca.
libs <- c(
"wdpar", # Para obtener datos cartográficos de áreas protegidas
"tidyverse",
"sf",
"stars",
"terra",
"maptiles",
"tmap"
)
invisible(
lapply(
libs, library,
character.only = TRUE
)
)
terraOptions(progress=0) # Per a no veure la barra de progrés#====================================
# DESCARGA Y PREPARACIÓN DE LOS DATOS
#====================================
#esp_raw_data <- wdpa_fetch("ESP", wait = TRUE)
#esp_clean_dat <- wdpa_clean(esp_raw_data) # Opcional puede ser largo
#esp_raw_data <- sf::st_as_sf(esp_raw_data)
#esp_raw_data <- st_transform(esp_raw_data, crs = 4326)
# Proceso opcional: Se realizan transformaciones y simplificaciones y estandardizaciones según lo descriro en:
# https://prioritizr.github.io/wdpar/reference/wdpa_clean.html
#=============================================
comarques <- st_read("../repo/divisions-administratives/divisions-administratives-v2r1-comarques-1000000-20240118.shp")## Reading layer `divisions-administratives-v2r1-comarques-1000000-20240118' from data source `C:\DataPortatil\Laboratori_R\Mapes\repo\divisions-administratives\divisions-administratives-v2r1-comarques-1000000-20240118.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 260160.2 ymin: 4488767 xmax: 527397.6 ymax: 4747976
## Projected CRS: ETRS89 / UTM zone 31Ncomarques <- st_transform(comarques, crs = 4326)
# ENCAJAR LOS DATOS EN EL CUADRO DE CATALUNYA
#============================================
#ymax <- 42.916667
#ymin <- 40.522919
#xmin <- 0.15925
#xmax <- 3.332488
#ylims <- c(ymin, ymax)
#xlims <- c(xmin, xmax)
#cat_box <- tibble(x = xlims, y = ylims) %>%
# st_as_sf(coords = c("x", "y"), crs= 4326) %>%
# st_bbox()
#cat_bbox <- st_bbox(cat_box) %>% st_as_sfc()
#cat_subset <- st_intersection(esp_raw_data, cat_bbox)
#plot(st_geometry(cat_subset))
#cat_subset_clean <- wdpa_clean(cat_subset) # Opcional
#plot(st_geometry(cat_subset_clean))
# Mejor guardar los datos obtenidos para las ocasiones siguientes
# ya que el roceso es bastante lento
#st_write(cat_subset, "../repo/wdpa-catalunya/wdpa-catalunya2.gpkg", append = FALSE)
#===================================================
# PROCESO QUE SE REALIZA A PARTIR DE LA SEGUNDA VEZ
#===================================================
cat_subset_clean <-
st_read("../repo/wdpa-catalunya/wdpa-catalunya2.gpkg")## Reading layer `wdpa-catalunya2' from data source
## `C:\DataPortatil\Laboratori_R\Mapes\repo\wdpa-catalunya\wdpa-catalunya2.gpkg'
## using driver `GPKG'
## Simple feature collection with 532 features and 28 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 0.15925 ymin: 40.52292 xmax: 3.332488 ymax: 42.8615
## Geodetic CRS: WGS 84cat_subset_clean <- st_transform(cat_subset_clean, crs = 4326)
#plot(st_geometry(ap_cat))
# diferentes manipulaciones para eliminar errores y mejorar la usabilidad:
# https://prioritizr.github.io/wdpar/reference/wdpa_clean.html
#cat_clean_dat <- wdpa_clean(ap_cat)
#==================================
# LIMPIEZA Y SIMPLIFICACIÓN
#==================================
# 1. Apartar las áreas que son 100% marinas ya que no entran en el análisis.
# 2. Agupar un poco las tipologías de parque para mejorar la comprensibilidad.
ap_cat_clean <- cat_subset_clean %>%
filter(DESIG != "Parque Natural") %>% # IMPORTANT!! Apartem els polígons dels "Parque Natural" perquè ja coberts els PEIN
filter(MARINE != "2") %>% # Se apartan las áreas 100% marinas
mutate(TERR_AREA= REP_AREA-REP_M_AREA) %>%
mutate(TIPO= case_match(DESIG,
"Ramsar Site, Wetland of International Importance"~
"Humedales (Ramsar)",
"Specially Protected Areas of Mediterranean Importance (Barcelona Convention)"~
"Zona SPAMI",
"Special Protection Area (Birds Directive)"~
"Espais Interès Natual (PEIN)",
"Special Areas of Conservation (Habitats Directive)"~
"Espais Interès Natual (PEIN)",
"Site of Community Importance (Habitats Directive)"~
"Espais Interès Natual (PEIN)",
"World Heritage Site (natural or mixed)"~ "UNESCO",
"Plan Especial de Protección (PEIN)"~ "Espais Interès Natual (PEIN)",
"UNESCO-MAB Biosphere Reserve"~ "UNESCO",
.default= as.character(DESIG)
))Para Catalunya, se obtiene un dataframe en formato
sf de 499 líneas y 31
columnas. Particularmente importantes, además de los polígonos,
las variables DESIG, que indica la tipología del área, y
MARINE que indica si se trata de una área marina, terrestre
o ambas cosas. Sobre estas dos variables se realizan unas pocas
operaciones de limpieza que se pueden ver en el código.
En realidad, estos datos presentan algun pequeño problema: el más importante que nos puede afecar para calcular las extensiones de las áreas protegidas es que en el caso del municipio de Olot -si se describe la cobertura utilizando todas la posibles capas- el centro de Olot (una ciudad de más de 30.000 hab.) queda incluido en el polígono de Parque Natural. Es un problema que, como se puede comprobar, también presenta el mapa del MITECO. Es decir que, si se dejara tal com está, el municipio de Olot quedaría ocupado por casi el 100% de su superficie por superficie protegida.
En este caso, obviamos la incidencia eliminando los polígonos Parque Natural, ya que todos quedan incluídos en los Espacios de Interés Natural (PEIN).
De todos modos, trabajando con datos de otras zonas de España, habrá que vigilar para prevenir que este posible problema no se repita.
Visualización estática de las áreas protegidas
En una primera visualización se puede apreciar la cobertura de Catalunya por las áreas protegidas y comprobar su coincidencia general con las visualizaciones del MITECO.
#===================================================
cat_bbox <- sf::st_bbox(
comarques
)
bg_map <- maptiles::get_tiles(
x = cat_bbox,
provider = "Esri.NatGeoWorldMap",
zoom = 10,
crop = TRUE,
project = FALSE
)
p <- ggplot() +
tidyterra::geom_spatraster_rgb( # tidyterra passes raster files to images
data = bg_map,
maxcell = terra::ncell( # set the cell number.
bg_map
),
alpha = 1
)+
geom_sf(data = ap_cat_clean,
aes(fill= TIPO),
alpha= 0.85)+
geom_sf(data=comarques, colour= "grey40", alpha=0, lwd=1) +
theme_minimal()+
#scale_fill_discrete(labels= function(x) str_wrap(x, 15))+
guides(fill=guide_legend(
title = NULL,
direction = "horizontal",
label.position= "top",
nrow = 1
)) +
theme(
axis.line = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
plot.background = element_blank(),
legend.background = element_rect(
fill = "white", color = NA),
legend.position = "bottom",
legend.text = element_text(size = 12),
)
library(grid)
grob1 <- grobTree(textGrob(
'Áreas terrestres protegidas de Catalunya',
x=.45, y=0.92, hjust = 0.1,
gp=gpar(col="darkgreen", fontsize=28, fontfamily= "serif",
fontface= "bold"))
)
grob2 <- grobTree(textGrob(
'Data from: "Protected Planet" - The World Database on Protected Areas (WDPA)',
x=.65, y=0.07, hjust=0.3,
gp=gpar(col="black", fontsize=13, fontface="italic"))
)
p + annotation_custom(grob1) + annotation_custom(grob2)
Unificación y “rasterización” de los polígonos
Sin embargo, nuestro interés no es clasificar y visualizar áreas protegidas, sino encontrar la manera de medir la superficie de cada municipio ocupada por alguna área protegida, independentemente de su tipología.
Para ello, la metodología que hemos encontrado como factible (hemos intentado otras que no nos han dado resultado) ha sido:
- “Rasterizar” los poligonos con
terra; - Intersecar las áreas generadas con los polígonos municipales;
- Calcular la extensión de las áreas al interior de los polígonos con
las librerías
terrayexactextractr.
# Rasterización
#======================================================
ap <- ap_cat_clean %>%
mutate(value= 1) %>% select(32) %>%
st_make_valid()
v1 <- vect(st_geometry(ap_cat_clean))
r1 <- rast(v1, resolution= 1/5000)
ap_rast <- rasterize(v1, r1,
#field = value,
touches = TRUE
)
#plot(ap_rast)Visualización de los datos rasterizados
Mapa rasterizado de la distribución de las áreas protegida, sin clasificación
tm_shape(municipis)+
tm_polygons(id= "NOMMUNI",
border.col = "steelblue",
alpha=0,
lwd= 0.8,
popup.vars=c("Superficie municipal: "= "AREAM5000"),
popup.format = list(AREAM5000= list(digits = 2))
) +
tm_shape(ap_rast) +
tm_raster(
palette = c("#287233"),
#colorNA = NULL,
alpha = 0.85,
legend.show = FALSE
) +
tm_shape(comarques) +
tm_borders(col="#073b7b", lwd=1.3)+
tmap_mode('view') +
tm_view(set.view = 8.2) +
tmap_options(basemaps = "Esri.NatGeoWorldMap",
basemaps.alpha = 1)