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” (http://www.murcianatural.carm.es/geocatalogo) 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( "http://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" ) )

Una vez instalado lo cargaremos en memoria:

library( terra )
library( mapview )

Descomprimiremos inicialmente todos los mapas, 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" )
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:

plot( murcia )
plot( zec, add = TRUE, col = 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 ESC 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 inormación adicional, que podemos representar mediante:

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

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 850 filas y 800 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 ESC 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 2012), 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 (http://centrodedescargas.cnig.es/CentroDescargas). 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.

plot( corine, col = rainbow( 5 ) )

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 )

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:

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

Rasterización de archivos vectoriales

Para realizar determinadas operaciones entre diferentes mapas resulta conveniente trabajar con objetos ráster, por lo que en ocasiones necesitaremos “rasterizar” objetos vectoriales. La función rasterize convierte un objeto vectorial en uno ráster, utilizando las especificaciones de otro ráster ya existente. Por defecto la función asigna NA al entorno, salvo que se especifique un valor determinado (por ejemplo, background = 0). Con el argumento field se pueden asisgnar atributos contenidos en el mapa vectorial, como en este caso el nombre de la ZEPA (zepa$site_name):

rasterize( zepa, mdt, field = "site_name", touches = TRUE ) -> rzepa
plot( rzepa, col = rainbow( 26 ) )
plot( zepa, add = TRUE )

Ejercicio 2

Crea el mapa ráster de la red de LIC-ZEC de la Región de Murcia a partir del mapa vectorial zec.

Ver la solución

rasterize( zec, mdt, field = "site_name", touches = TRUE ) -> rzec
plot( rzec, col = rainbow( 50 ) )

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

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). El resultado puede agregarse en forma de variable al objeto:

extract( mdt, zepa, mean, na.rm = TRUE )
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 )

La función equivalente a extract, para usar entre dos archivos ráster es zonal:

zonal( mdt, rzepa, mean, na.rm = TRUE )

En el caso del mapa de usos del suelo, el atributo (el código CORINE) no tiene significado cuantitativo, por lo que la función zonal no resulta apropiada. En cambio, podemos usar:

crosstab( c( rzepa, corine ) ) -> zepacor
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 4, Estepas de Yecla), o en el conjunto, utilizaremos:

zepacor / rowSums( zepacor ) * 100
zepacor[ 4, ] / sum( zepacor[ 4, ] ) * 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[ 14, ] / sum( zepacor[ 14, ] )

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 )
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( rzepa, strix )
table( extract( rzepa, strix )[ , 2 ], 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 con crosstab:

crosstab( c( rzepa, strix_ih ) ) -> zepastrix
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" ).