Mapas ráster - aplicaciones ecológicas

Práctica 4, Datos Espaciales en Biodiversidad

Máster en Áreas Protegidas, Recursos Naturales y Biodiversidad

Profesores

José Francisco Calvo Sendín

María Victoria Jiménez Franco

Departamento

Ecología e Hidrología, Facultad de Biología

Fecha de publicación

28 de noviembre de 2023

Fecha de modificación

29 de febrero de 2024

Introducción

En esta práctica realizaremos diversos ejercicios con mapas ráster con los que aprenderemos a hacer recortes y ajustes de dimensiones y extensiones. Trabajaremos también con la estructura de capas de los objetos ráster, con diversos aspectos de la representación, y aprenderemos a modificar la resolución espacial y los valores de las celdas. Finalmente realizaremos operaciones algebraicas y calcularemos las predicciones de modelos estadísticos elaborados con la información almacenada en las capas de los mapas ráster.

Preparación

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

load( url( "http://webs.um.es/jfcalvo/deb.RData" ) )

Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "mdt200" ) o info( "taxiados" ). 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 geodata. Si no los tenemos instalados usaremos:

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

Una vez instalados los cargaremos en memoria:

library( terra )
library( geodata )

Descomprimiremos también los mapas que necesitaremos en la práctica, usando la función unwrap, y cambiaremos el sistema de referencia de coordenadas del mapa murcia:

unwrap( mdt200 ) -> mdt200
unwrap( murcia ) -> murcia
unwrap( strix_ih ) -> strix_ih
unwrap( taxiados ) -> taxiados
unwrap( cnv ) -> cnv
unwrap( zepa ) -> zepa
unwrap( zec ) -> zec
project( murcia, crs( mdt200 ) ) -> murcia

Recortes, dimensión y extensión

Trabajaremos incialmente con el mapa ráster de altitudes (modelo digital de elevaciones):

plot( mdt200 )
plot( murcia, add = TRUE )

El recorte, en función de un mapa vectorial (los límites territoriales de la Región de Murcia en este caso), se realiza con la función crop especificando el argumento mask = TRUE:

crop( mdt200, murcia, mask = TRUE ) -> mdt
mdt
plot( mdt )

No obstante, si queremos trabajar con varios mapas ráster conjuntamente, es necesario que todos tengan, además del mismo sistema de referencia de coordenadas, la misma resolución (tamaño de celda), dimensión (número de filas y columnas) y extensión (coordenadas máximas y mínimas). En este caso, si comparamos el resultado del recorte con el mapa ráster de idoneidad de hábitat del cárabo, comprobaremos que ni las dimensiones ni la extensión son las mismas:

strix_ih
plot( strix_ih )

Deberíamos, por tanto, utilizar un único mapa ráster de referencia, con el que dimensionar el resto. En este caso, la opción más cómoda sería realizar el recorte con el mapa ráster previamente ajustado a los límites de la Región de Murcia (strix_ih):

crop( mdt200, strix_ih, mask = TRUE ) -> mdt
plot( mdt )
plot( strix_ih, add = TRUE, alpha = 0.2, legend = "topleft" )

El recorte entre dos mapas ráster también puede realizarse con la función mask, pero solo cuando las extensiones de ambos mapas son las mismas.

Cuando los mapas presentan sistemas de referencia de coordenadas diferentes, el ajuste requiere, además del uso de la función project, aplicar la función resample. Para ver un ejemplo, descargaremos los mapas bioclimáticos de España desde la web de WorldClim:

worldclim_country( country = "Spain" , var = "bio", path = "./descargas/" ) -> spain_bio
spain_bio
plot( spain_bio )

A continuación cambiaremos el sistema de referencia de coordenadas (con project) y aplicaremos resample para igualar las dimensiones, la resolución y la extensión:

project( spain_bio, crs( strix_ih ) ) -> spain_bio
resample( spain_bio, strix_ih ) -> murcia_bio
murcia_bio
plot( murcia_bio )

Finalmente realizaremos el recorte:

crop( murcia_bio, strix_ih, mask = TRUE ) -> murcia_bio
plot( murcia_bio )

Recuerda que el mapa de variables bioclimáticas consta de 19 capas cuyos códigos pueden consultarse aquí.

Las coordenadas de un mapa ráster se obtienen con la función crds. Por ejemplo:

crds( strix_ih )
crds( mdt )[ data.frame( mdt ) > 1980, ]

Capas

Como acabamos de comprobar en la sección anterior, una característica importante de los mapas ráster es que pueden constar de diferentes capas que almacenan información sobre diferentes variables expresadas geográficamente sobre el mismo territorio. En R, estas capas se pueden manejar como los elementos de un objeto de tipo “lista”:

murcia_bio[[ 1 ]]
plot( murcia_bio[[ 1 ]] )

Como ejercicio crearemos un mapa ráster con varias capas de variables ambientales: temperarura media anual (capa 1 de murcia_bio), precipitación anual (capa 12 de murcia_bio) y altitud (mapa mdt):

amb <- c( murcia_bio[[ 1 ]], murcia_bio[[ 12 ]], mdt )
names( amb ) <- c( "temp", "prec", "alt" )
amb
plot( amb )

Para añadir más capas se puede usar la función add. Por ejemplo, añadiremos una capa con las localizaciones de los taxiados, que rasterizaremos previamente usando como atributo los valores de riqueza (número de especies de aves registradas en cada taxiado):

add( amb ) <- rasterize( taxiados, mdt, field = "riqueza", touches = TRUE )
amb
plot( amb )

El argumento touches = TRUE asegura que la rasterización afecte a todas las celdas que son “tocadas” por las líneas o polígonos, no solo aquellas en la ruta de representación de la línea o cuyo punto central está dentro del polígono.

La extracción de información ambiental del mapa se realiza para todas las capas simultáneamente. Por ejemplo:

extract( amb, centroids( taxiados ) )
global( amb, mean, na.rm = TRUE )
Ejercicio 1

Crea un mapa ráster con las siguientes capas: rango anual de temperatura, idoneidad de hábitat del cárabo y localización de territorios del cárabo.

mapa <- c( murcia_bio[[ 7 ]], strix_ih, strix = rasterize( vect( strix ), strix_ih ) )
plot( mapa )

Modelos

En esta sección comprobaremos la utilidad de los mapas ráster para analizar información biológica en función de información ambiental. Utilizaremos el mapa amb elaborado en la sección anterior para realizar un sencillo modelo de regresión de Poisson que estime la riqueza de especies de aves en función de las variables ambientales:

glm( riqueza ~ temp + prec + alt, data = amb, family = poisson ) -> modelo
summary(modelo)
plot( predict( amb, modelo, type = "response" ) )

De forma alternativa, se puede extraer la información de las capas del mapa (con extract), almacenarla en una tabla de datos, ejecutar el modelo con las variables de la tabla, y utilizar finalmente predict con el mapa. Para que predict funcione hay que asegurarse de que los nombres de las variables de la tabla de datos coincidan con los de las capas del mapa ráster.

Resolución espacial

La resolución espacial de un mapa ráster hace referencia al tamaño de las celdas o píxeles que lo conforman. Cuando la resolución es excesivamente detallada, o cuando debemos ajustar la resolución de dos mapas para trabajar con ellos conjuntamente, necesitaremos realizar cambios de resolución. Por ejemplo, el objeto carnívoros es un mapa ráster de tipo átlas, que contiene datos de presencia de varias especies de mamíferos carnívoros en la Región de Murcia en cuadrículas UTM de 5 x 5 kilómetros:

cnv
plot( cnv )
plot( cnv[[ 4 ]], main = "Gato montés" )
plot( murcia, add = TRUE )

Si queremos analizar la distrubución de las especies en función de variables ambientales, tendremos que crear un mapa con menor resolución agregando las celdas del mapa original. Realizaremos un ejemplo con el mapa de altitudes. El proceso requiere varios pasos.

En primer lugar debemos asegurar que las extensiones coincidan. Usaremos crop sobre mdt200 con el argumento mask = FALSE, de modo que apliquemos el recorte antes de agregar las celdas:

crop( mdt200, cnv, mask = FALSE ) -> mdt5
mdt5
plot( mdt5 )

A continuación usaremosla función aggregate, especificando el número de celdas a agregar y la función para calcular el valor de las nuevas celdas (por defecto, la media).

aggregate( mdt5, 25, na.rm = TRUE ) -> mdt5
mdt5
plot( mdt5 )

Finalmente haremos el recorte (ahora sí con el argumento mask = TRUE), ajustaremos la extensión de mdt5 a la de cnv (función crop y argumento ext = TRUE), y añadiremos la capa de elevaciones al objeto cnv. Usamos el mapa vectorial murcia para aplicar la máscara porque ninguna de las capas del mapa cnv contiene todas las cuadrículas 5 x 5 de la Región. Además, añadimos el argumento snap = "out" para que ninguna cuadrícula limítrofe con escasa superficie regional se quede fuera del recorte.

crop( mdt5, murcia, mask = TRUE, snap = "out" ) -> mdt5
mdt5
plot( mdt5 )
crop( mdt5, cnv, ext = TRUE ) -> mdt5
mdt5
plot( mdt5 )
names( mdt5 ) <- "alt"
add( cnv ) <- mdt5
plot( cnv )

Como alternativa a la función aggregate puede usarse la función resample con el argumento method = "average". En este caso nos ahorraríamos el primer paso del procedimiento anterior:

resample( mdt200, cnv, "average" ) -> mdt5b
crop( mdt5b, murcia, mask = TRUE, snap = "out" ) -> mdt5b
crop( mdt5b, cnv, ext = TRUE ) -> mdt5b
mdt5b
plot(mdt5b)
Ejercicio 2

Repite el ejercicio realizado en esta sección con las variables 1 (temperatura media anual) y 12 (precipitación anual) del mapa spain_bio. En este caso es necesario utilizar resample en lugar de aggregate.

resample( spain_bio[[ c( 1, 12 ) ]], cnv ) -> bio5
crop( bio5, murcia, mask = TRUE, snap = "out" ) -> bio5
crop( bio5, cnv, ext = TRUE )-> bio5
names( bio5 )<- c( "temp", "prec" )
add( cnv ) <- bio5

Manejo de mapas ráster

Continuaremos trabajando en esta sección con el mapa de los mamíferos carnívoros para realizar unos modelos de regresión logística y estimar la probabilidad de presencia en función de variables ambientales.

Dado que las capas de distribución de las especies solo tienen valores 1, y que la regresión logística requiere la existencia de valores 0, tendremos que asignar este valor de base a las cuadrículas regionales 5 x 5. Procederemos sustituyendo los valores NA por 0 y aplicando un recorte:

subst( cnv[[ 1 : 7 ]], NA, 0 ) -> cnv[[ 1 : 7 ]]
plot( cnv )
mask( cnv[[ 1 : 7 ]], cnv$alt ) -> cnv[[ 1 : 7 ]]
plot( cnv )

A continuación aplicaremos el modelo de regresión logística con los datos de una de las especies y realizaremos la predicción de las probabilidades estimadas:

glm( vulpes ~ alt, data = cnv, family = binomial ) -> mzorro
summary( mzorro )
plot( predict( cnv, mzorro, type = "response" ) )
Ejercicio 3

Elabora un mapa de probabilidad de presencia estimada de la garduña (Martes foina), en función de la altitud, en cuadrículas UTM 5 x 5 de la Región de Murcia.

glm( martes ~ alt, data = cnv, family = binomial ) -> mmartes
plot( predict( cnv, mmartes, type = "response" ) )

El cambio de valores de un ráster puede realizarse también con la función classify. Por ejemplo, podemos crear un nuevo mapa de elevaciones de acuerdo con unos puntos de corte:

classify( mdt200, c( 1, 250, 500, 1000, 2000 ) ) -> mdt200c
plot( mdt200c, legend = "bottomright" )

Operaciones

Realizaremos en esta sección diversos ejercicios con los mapas rasterizados de la red Natura 2000 en la Región de Murcia (mapas zepa y zec):

rasterize( zepa, mdt, touches = TRUE ) -> rzepa
rasterize( zec, mdt, touches = TRUE ) -> rzec
plot( c( rzepa, rzec ) )

Al no especificar ningún valor para el argumento field, la rasterización proporciona mapas con valor 1 en las cuadrículas que pertenecen a algún espacio de la red.

Podemos crear un nuevo objeto uniendo las dos capas:

rn2k <- c( rzepa, rzec )
names( rn2k ) <- c( "zepa", "lic_zec" )
plot( rn2k )

Con las capas de este objeto (o con los mapas rasterizados) podemos realizar operaciones algebraicas. Por ejemplo:

plot( sum( rn2k, na.rm = TRUE ) )
freq( sum( rn2k, na.rm = TRUE ) )
plot( sum( rn2k ) )
freq( sum( rn2k ) )
plot( max( rn2k, na.rm = TRUE ) )
freq( max( rn2k, na.rm = TRUE ) )

Multiplicando por 4 los números de cuadrículas obtenemos la superficie total (en hectáreas) protegida por la red (incluidas las áreas marinas) y la superficie solapada por ambas redes (ZEPA y LIC/ZEC).

freq( max( rn2k, na.rm = TRUE ) ) * 4
freq( sum( rn2k ) ) * 4

Alternativamente podemos usar la función expanse:

expanse( sum( rn2k, na.rm = TRUE ), transform = FALSE )
expanse( sum( rn2k ), transform = FALSE )

La función cover produce el mismo resultado que max, aunque en este caso hay que especificar el nombre de ambos mapas:

plot( cover( rn2k$zepa, rn2k$lic_zec ) )
Ejercicio 4

Elabora un mapa de riqueza (número de especies) de carnívoros a partir de las capas 1-7 del mapa cnv y realiza un modelo de regresión de Poisson usando la altitud como variable explicativa.

sum( cnv[[ 1 : 7 ]] ) -> cnv$riqueza
plot( cnv )
glm( riqueza ~ alt, data = cnv, family = poisson ) -> mriqueza
summary( mriqueza )
plot( predict( cnv, mriqueza, type = "response" ) )

Evaluación

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

Bibliografía

Descripción de los mapas

Utiliza la función info. Por ejemplo: info( "mdt200" ) o info( "taxiados" ).