Análisis de redes de áreas protegidas

Práctica 2, Gestión y Conservación de Ecosistemas

Grado en Biología

Introducción

En esta práctica analizaremos diversas características de las redes de áreas protegidas mediante un paquete específico de SIG (Sistema de Información Geográfica) en R. Utilizaremos mapas de redes de áreas protegidas en formato vectorial, mapas de puntos (distribución de especies) y mapas ráster de variables ambientales. Aprenderemos a representarlos y manejarlos, calcularemos áreas y perímetros y realizaremos conversiones de formatos de datos espaciales. Analizaremos también la representatividad de las redes en función de variables ambientales y de uso del suelo, y realizaremos un sencillo análisis gap para evaluar la eficacia de una red en la protección de los territorios de una determinada especie.

Trabajaremos principalmente con archivos de mapas oficiales de la Red Natura 2000 de la Región de Murcia, disponibles en el “Geocatálogo” de la Comunidad Autónoma, donde podemos usar el buscador de la página para encontrar los mapas correspondientes a las redes ZEPA y LIC-ZEC (en formato shp), que se proporcionan en sendos archivos comprimidos (formato zip).

Preparación

Los mapas necesarios para el desarrollo de la práctica se encuentran disponibles, en formato comprimido, en el archivo RData de la asignatura. Una vez iniciado R, cargaremos este archivo desde el servidor:

load( url( "https://webs.um.es/jfcalvo/gce.RData" ) )

Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "zepa" ) o info( "mdt" ). Teclea info() para más información.

Para realizar esta práctica se necesita la instalación adicional de los paquetes de R terra y mapview. Si no los tenemos instalados ejecutaremos:

install.packages( c( "terra", "mapview" ) )

Si estamos trabajando en el Aula Virtual, por problemas de compatibilidad entre versiones de paquetes de R, necesitaremos ejecutar lo siguiente:

Sys.unsetenv( "PROJ_LIB" ) 

A continuación, cargaremos los paquetes en memoria:

library( terra )
library( mapview )

Los mapas que usaremos en la práctica están en un formato comprimido, por lo que los descomprimiremos usando la función unwrap:

unwrap( zepa ) -> zepa
unwrap( zec ) -> zec
unwrap( murcia ) -> murcia
unwrap( strix ) -> strix
unwrap( mdt ) -> mdt
unwrap( corine ) -> corine
unwrap( strix_ih ) -> strix_ih

Mapas vectoriales

En primer lugar, trabajaremos con el objeto que contiene los límites de las Zonas de Especial Protección para las Aves (ZEPA) de la Región de Murcia. Visualizaremos esta cartografía con:

plot( zepa )

Los límites de la Región de Murcia están disponibles en el objeto murcia. Si queremos representarlos sobre el mapa de las ZEPA utilizaremos la función plot con el argumento add = TRUE, pudiendo especificar un color distinto (argumento border):

plot( murcia, add = TRUE, border = 4 )

El objeto zepa contiene la información geográfica sobre las áreas designadas ZEPA en la Región de Murcia, que podemos consultar mediante:

zepa
data.frame(zepa)

En el primer caso accedemos a la información sobre la clase del objeto (polígonos SpatVector), el número de elementos que lo componen (en este caso 26 áreas protegidas con cuatro variables), e información sobre coordenadas y proyección geográfica. Podemos comprobar que estamos usando la proyección UTM 30 y el datum ETRS89). Para trabajar juntamente con diversos mapas es necesario que todos tengan la misma proyección y datum.

En el segundo caso, obtenemos los datos asociados al objeto: código, nombre y comunidad autónoma del espacio. Es importante señalar que esta información no siempre es la misma, no siempre está disponible, o no siempre es correcta, en todos los archivos a los que podamos acceder desde distintos servidores cartográficos.

Dado que el objeto zepa tiene la estructura de una tabla de datos, podemos referirnos en particular a una “columna” (una determinada variable) o a una “fila” (una determinada ZEPA):

zepa$site_name
zepa[ 11, ]
plot( zepa [ 11, ], add = TRUE, col = "blue" )

También podemos seleccionar los espacios mediante operadores lógigos:

plot( zepa [ zepa$site_name == "Sierra de la Pila", ], add = TRUE, col = 6 )

Por su parte, la cartografía de los espacios LIC/ZEC (Lugares de Importancia Comunitaria/Zonas de Especial Conservación) de la Región de Murcia está disponible en el objeto zec:

zec
plot( murcia )
plot( zec, add = TRUE, col = map.pal( "rainbow",  50 ) )

Existen diversas funciones que nos permiten obtener información de un mapa vectorial en la ventana gráfica:

click()

Si queremos obtener más información, además de las coordenadas, especificaremos el nombre del objeto. En este caso es necesario presionar la tecla ESCESC para finalizar.

click( zec )

Otra función útil es zoom:

zoom( zec )

Ejercicio 1

Utiliza la cartografía de la red de LIC-ZEC de la Región de Murcia (objeto zec) para obtener las coordenadas UTM de un punto central del LIC Cabezo de Roldán.

Ver la solución

plot( zec[ zec$site_name == "Cabezo de Roldán", ] )
click()

Mapas de puntos

Frecuentemente la información sobre la distribución de una especie está referida a puntos o cuadrículas del territorio, y se recoge en archivos de coordenadas x, y. Para la práctica disponemos del archivo de presencias (territorios) del cárabo cómun (Strix aluco) en la Región de Murcia:

strix

Comprobamos que se trata también de un objeto de clase SpatVector, con 77 territorios, sin información adicional, que podemos representar mediante:

plot( strix )
plot( murcia, add = TRUE )
plot( zepa, add = TRUE, border = 4 )

Visualización interactiva

El paquete mapview permite la representación online interactiva de mapas:

mapview( zepa, col.regions = "red" ) -> mapa1
mapview( zec, col.regions = "green" ) -> mapa2 
mapview( strix ) -> mapa3
mapa1 + mapa2 + mapa3

Mapas ráster

Los mapas ráster se presentan en múltiples formatos de archivo. Para la práctica disponemos del objeto mdt que contiene el mapa de elevaciones de la Región de Murcia, elaborado a partir del modelo digital del terreno (o mapa digital de elevaciones con paso de malla de 200 m del Instituto Geográfico Nacional (descargable en http://centrodedescargas.cnig.es/CentroDescargas).

mdt

Podemos comprobar que la resolución de este mapa es de 200 m (cuadrículas o celdas de 200 x 200 m); en total 769 filas y 750 columnas.

Para representar un mapa ráster se usa también la función plot, con la posibilidad de superponerles mapas vectoriales:

plot( mdt )
plot( murcia, add = TRUE )
plot( zepa, add = TRUE, border = "red" )
plot( strix, add = TRUE, col = "blue" )

Con los mapas ráster podemos interaccionar de igual manera que con los vectoriales (recuerda presionar ESCESC para salir):

click( mdt, xy = TRUE )

El análisis cuantitativo de los valores del mapa puede realizarse con la función global, o alternativamente con la función values. Por ejemplo:

global( mdt, mean, na.rm = TRUE )
mean( values( mdt ), na.rm = TRUE )

Utilizaremos ahora otro mapa ráster, pero en este caso con información cualitativa. Se trata de un mapa regional de usos del suelo previamente “rasterizado” a partir del mapa de usos CORINE Land Cover (año 1990), la base de datos europea de ocupación del suelo con escala de referencia 1:100.000, basada en una nomenclatura de 44 clases (3 niveles), y al que podemos acceder desde el centro de descargas del Instituto Geográfico Nacional . El mapa (objeto corine) representa los usos al nivel 1 de detalle para la Región de Murcia. La correspondencia entre los códigos y las categorías del mapa es la siguiente: (1) superficies artificiales; (2) zonas agrícolas; (3) zonas forestales con vegetación natural y espacios abiertos; (4) zonas húmedas; y (5) superficies de agua.

corine
plot( corine )

En este caso no tiene interés calcular estadísticos como la media, pero sí conocer la frecuencia de cada categoría de uso. Para ello recurriremos a:

freq( corine )

Operaciones con mapas: cálculo de áreas y perímetros

El cálculo de áreas y perímetros de mapas vectoriales se puede realizar con las funciones expanse y perim, respecitvamente, con la posibilidad de añadirlos como nuevas variables del objeto. Por ejemplo, para la red de ZEPA:

expanse( zepa ) -> zepa$area
perim( zepa ) -> zepa$perimetro
data.frame( zepa )

Si queremos estimar la proporción de superficie regional protegida por la red de LIC/ZEC utilizaremos:

sum( expanse( zec ) ) / expanse( murcia )

El cálculo, no obstante, está sobreestimado debido a que hemos incluido también los espacios marinos (26: Valles submarinos del escarpe de Mazarrón; 37: Franja litoral sumergida de la Región de Murcia; 45: Mar Menor). Para excluirlos recurrimos a eliminar las filas 26, 37 y 45 del objeto zec:

zec$site_name
zec[ -c( 26, 37, 45 ), ] -> zec_terr
plot( murcia )
plot( zec_terr, add = TRUE )
sum( expanse( zec_terr ) ) / expanse( murcia )

Ejercicio 3

Calcula, para la red de ZEPA, el porcentaje de superficie regional protegida. Para ello deberás crear un nuevo objeto zepa_terr, a partir del objeto zepa, excluyendo las ZEPA “Espacio marino de Tabarca-Cabo de Palos”, “Espacio marino de los islotes litorales de Murcia y Almería” y “Mar Menor”.

Ver la solución

zepa[ -c( 14, 15, 19 ), ] -> zepa_terr
plot( murcia )
plot( zepa_terr, add = TRUE )
sum( expanse( zepa_terr ) ) / expanse( murcia ) * 100

Para calcular superficies con mapas ráster es muy importante tener en cuenta la resolución. En nuestro caso, cada cuadrícula en el mapa corresponde a una cuadrícula de 200 × 200 m en el terreno (4 hectáreas). Así, por ejemplo, podremos calcular la superficie de la Región de Murcia, en hectáreas mediante:

freq( corine )
freq( corine )[ , 3 ] * 4
sum( freq( corine )[ , 3 ] ) * 4

Alternativamente, podemos usar también la función expanse:

expanse( corine )

Operaciones con mapas: representatividad ambiental

El análisis de las características de un mapa vectorial (o de puntos) en función de los atributos de un mapa ráster se realiza con extract, especificando una función (media, máximo, mínimo…). También es necesario utilizar el argumento na.rm = TRUE, para eliminar de los cálculos las cuadrículas sin valor (NA).

extract( mdt, zepa, mean, na.rm = TRUE )

El resultado puede agregarse en forma de variable al objeto:

extract( mdt, zepa, mean, na.rm = TRUE )[ , 2 ] -> zepa$altitud
data.frame( zepa )

Cuando se trata de mapas de puntos, como en el caso del cárabo, no hace falta especificar una función:

extract( mdt, strix )[ , 2 ] -> strix$altitud
data.frame( strix )

En el caso del mapa de usos del suelo, el atributo (el código CORINE) no tiene significado cuantitativo, por lo que el cálculo de valores medios no resulta apropiado. En cambio, podemos usar:

extract( corine, zepa, table, touches = TRUE )

El argumento touches = TRUE asegura que la inclusión de todas las celdas con superficie dentro de los límites de los espacios, no solo aquellas cuyo punto central esté incluido.

Con los resultados de la función anterior construiremos una tabla con la que podremos realizar diversas operaciones:

extract( corine, zepa, table, touches = TRUE )[ , -1 ] -> zepacor
rownames( zepacor ) <- zepa$site_name
zepacor

Con esta tabla podemos operar como con cualquier otra en R. Así, para estimar el porcentaje de los diferentes usos en cada uno de los espacios de la red de ZEPA, en uno de ellos en particular (por ejemplo, el número 2, Estepas de Yecla), o en el conjunto, utilizaremos:

zepacor / rowSums( zepacor ) * 100
zepacor[ 2, ] / sum( zepacor[ 2, ] ) * 100
colSums( zepacor ) / sum( zepacor ) * 100

Ejercicio 4

• Calcula la media de los valores medios de altitud en los LIC-ZEC terrestres de la Región de Murcia usando la función extract.

• ¿Cuál es la altitud media de los territorios del cárabo común?

• Estima la proporción de usos del suelo en la ZEPA Saladares del Guadalentín.

Ver la solución

• Calcula la media de los valores medios de altitud en los LIC-ZEC terrestres de la Región de Murcia usando la función extract.

   mean( extract( mdt, zec, mean, na.rm = TRUE )[ , 2 ], na.rm = TRUE )

• ¿Cuál es la altitud media de los territorios del cárabo común?

   mean( strix$altitud )

• Estima la proporción de usos del suelo en la ZEPA Saladares del Guadalentín.

   zepacor[ 13, ] / sum( zepacor[ 13, ] )

Análisis gap: el cárabo común y la red de ZEPA

Además del mapa de puntos que representa la distribución de territorios de cárabos, disponemos de un mapa ráster de “idoneidad de hábitat” (objeto strix_ih) con cuatro clases que representan el hábitat inadecuado (1), marginal (2), adecuado (3) y óptimo (4).

plot( strix_ih )
plot( strix, add = TRUE, col = 2 )
plot( murcia, add = TRUE )

Analizaremos ambos mapas en relación con la red de espacios ZEPA. En este caso, como se trata de puntos, al usar extract no es necesario especificar una función:

extract( zepa, strix )
table( extract( zepa, strix )[ , 3 ], useNA = "always" )

Por su parte, la relación entre el mapa ráster de idoneidad de hábitat y la red ZEPA se puede analizar elaborando una tabla similar a la elaborada con el mapa corine:

extract( strix_ih, zepa, table, touches = TRUE )[ , -1 ] -> zepastrix
rownames( zepastrix ) <- zepa$site_name
zepastrix

Ejercicio 5

• ¿Qué ZEPA alberga el mayor número de territorios? ¿Cuántos?

• ¿Cuál es el porcentaje de territorios situados fuera de la red de ZEPA?

• ¿Qué ZEPA presenta la mayor proporción de hábitat óptimo? ¿Qué proporción?

Evaluación

Realiza la prueba de evaluación disponible en la herramienta Exámenes del Aula Virtual.

Bibliografía

• Hijmans R (2023) _terra: Spatial Data Analysis_. R package version 1.7-55, https://CRAN.R-project.org/package=terra

Descripción de los mapas

Utiliza la función info. Por ejemplo: info( "zepa" ) o info( "mdt" ).