Análisis de uso del hábitat

Práctica 7, Datos Espaciales en Biodiversidad

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

Profesores
Departamento

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

Fecha de publicación

5 de diciembre de 2023

Fecha de modificación

9 de abril de 2024

Introducción

En esta práctica realizaremos un ejercicio de análisis de uso del hábitat a partir de datos de localizaciones procedentes de estudios de seguimiento de animales por telemetría. Estimaremos áreas de campeo y realizaremos análisis de selección de hábitat, además de diversos ejercicios previos con mapas vectoriales y raster para obtener y gestionar la información cartográfica necesaria para la obtención de datos sobre las características del hábitat.

Preparación

Los datos, 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( "rl" ) o info( "corine90" ). Teclea info() para más información.

Para realizar esta práctica se necesita la instalación adicional de los paquetes de R terra, adehabitatHR, ctmm y adehabitatHS. Si no los tenemos instalados usaremos:

install.packages( c( "terra", "adehabitatHR", "ctmm", "adehabitatHS" ) )

Una vez instalados los cargaremos en memoria:

library( adehabitatHR )
library( ctmm )
library( adehabitatHS)
library( terra )

Es importante cargar el paquete terra después del paquete ctmm, porque ambos incluyen una función plot, y la versión de ctmm no permite representar los objetos con formato de terra. Ocurre lo mismo con la función buffer del paquete adehabitatMA que se carga automáticamente con el paquete adehabitatHR. No obstante, siempre se puede revertir la carga en memoria de cualquier paquete con, por ejemplo, detach( package : ctmm ), o utilizar cualquier función de cualquier paquete utilizando, por ejemplo: terra::plot().

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( murcia ) -> murcia
unwrap( zepa ) -> zepa
unwrap( corine90 ) -> corine90
project( murcia, crs( zepa ) ) -> murcia

Radiolocalizaciones

Utilizaremos el objeto rl que contiene los datos de radiolocalizaciones de 6 águilas calzadas (Hieraaetus pennatus), procedentes de un estudio de seguimiento realizado en la ZEPA Sierras de Burete, Lavia y Cambrón (ES0000267) entre 1999 y 2004 (Martínez et al. 2007). En total son 408 radiolocalizaciones de seis individuos (cuatro hembras y dos machos):

rl

La tabla incluye también la variable caza que indica si el animal estaba realizando o no actividades de caza. Para la estimación de las áreas de campeo utilizaremos todas las localizaciones, mientras que para el análisis de uso del hábitat utilizaremos solo las registradas como actividad de caza.

Comenzaremos elaborando un mapa de puntos con la función vect del paquete terra. Las coordenadas de las radiolocalizaciones corresponden al datum ED50, por lo que los transformaremos a ETRS89 con la función ed50:

vect( rl, geom = c( "x", "y" ) ) -> rl
ed50( rl ) -> rl
rl
plot( murcia )
plot( zepa, add = T )
plot( rl, "individuo", add = TRUE, cex = 0.5, legend = "topright" )

Calcularemos ahora, para cada individuo, el porcentaje de radiolocalizaciones registradas dentro y fuera de la ZEPA. Primero seleccionaremos la ZEPA:

zepa$site_name
zepa[ 11, ] -> burete

A continuación procederemos a extraer la información de las radiolocalizaciones:

extract( burete, rl )
extract( burete, rl )$site_code -> rl$zepa
rl
table(rl$individuo, rl$zepa, useNA = "always" ) -> inout
inout
inout / rowSums( inout ) * 100

Áreas de campeo: MCP y kernel

Para estimar las áreas de campeo (home range) usaremos funciones del paquete adehabitatHR. Este paquete requiere objetos de la clase SpatialPoints, por lo que es necesario convertir las radiolocalizaciones (mediante la función SpatialPointsDataFrame). En primer lugar, estimaremos las áreas de campeo mediante el método del polígono mínimo convexo (MCP):

hr <- SpatialPointsDataFrame( geom( rl )[  ,3 : 4 ], data = data.frame( rl[ , 1 ] ) )
hr
mcp.area( hr )
mcp( hr, 100 )

También estimaremos las áreas de campeo mediante el método de kernel:

kernelUD( hr ) -> ud
kernel.area( ud )
kernel.area( ud, 95 )

La representación conjunta de las áreas estimadas con ambos métodos se puede realizar con:

plot( mcp( hr, 100 ), border = 1 : 6 )
plot( getverticeshr( ud, 95 ), add = TRUE , border = 1 : 6, lwd = 2 )
plot( burete, add = TRUE )

Áreas de campeo: kernel autocorrelacionado

Las tecnologías de seguimiento modernas han incrementado sustancialmente el volumen de información sobre los procesos relacionados con el movimiento de los animales. Este importante avance tecnológico de la telemetría viene acompañado de nuevas técnicas de análisis de datos, que tienen en consideración especialmente las cuestiones relacionadas con la autocorrelación de las localizaciones obtenidas. El paquete ctmm y sus funciones proporcionan técnicas de análisis kernel autocorrelacionado que se adecúan a las características de este tipo de datos.

Como ejemplo utilizaremos un conjunto más reducido de radiolocalizacioes para las que disponemos la fecha y hora de cada registro. En concreto usaremos el objeto rlt que contiene los datos del individuo número 6 correspondientes al seguimiento en 2003:

rlt

Convertiremos la tabla de datos al formato requerido por el paquete ctmm. Ajustaremos la información temporal a un modelo autocorrelacionado, y aplicaremos la función akde (autocorrelated kernel density) para estimar áreas de campeo:

as.telemetry( rlt ) -> rlt
fit <- ctmm.select( rlt, ctmm.guess( rlt ) )
ud6 <- akde( rlt, fit )
plot( rlt , UD = ud6 )
summary( ud6 )

El objeto resultante puede convertirse al formato de mapas vectorial de terra aplicacando las funciones SpatialPolygonsDataFrame.UD y vect:

SpatialPolygonsDataFrame.UD( ud6 ) -> ud6
vect( ud6 ) -> ud6
project( ud6, crs( murcia ) ) -> ud6
plot( murcia )
plot( ud6, add =TRUE )

Análisis de selección de hábitat: preparación del mapa CORINE

Realizaremos un anális de uso del hábitat por parte de las águilas considerando por un lado las radiolocalizaciones dentro de las áreas de campeo (usaremos las áreas estimadas por MCP), y por otro las áreas de campeo dentro del área de estudio. Para trabajar con comodidad seleccionaremos este área de estudio en función de la extensión del mapa de áreas de campeo:

ext( mcp( hr, 100 ) )

Recortaremos el mapa corine90 a partir de estas coordenadas, pero ampliando los márgenes en, por ejemplo, 1000 metros cada lado:

crop( corine90, ext( mcp( hr, 100 ) ) + c( 1000, 1000, 1000, 1000 ) ) -> corine
plot(corine, "CODE_90" )
plot( mcp ( hr ), add = TRUE, lwd = 2, border = "red" )
plot( burete, add = TRUE, lwd = 3 )

Para nuestro análisis crearemos una nueva categoría de uso no incluida en el CORINE que represente los ecotonos, es decir, los bordes o fronteras entre polígonos de diferente tipo de uso. Por otra parte, es conveniente reducir el número total de usos con objeto de simplificar el análisis y facilitar la intepretación de los resultados. Esto podemos hacerlo agrupando al nivel 1, de forma similar a como lo hicimos en un ejercicio de la práctica 3. En este caso mantendremos los bordes con el argumento dissolve = FALSE, ya que los necesitaremos para crear la capa de ecotonos.

round( as.numeric( corine$CODE_90 ) / 100, 0 ) -> corine$n1
aggregate(corine, "n1", dissolve = FALSE ) -> corine
corine
plot( corine, "n1" )

Finalmente eliminaremos las variables innecesarias del objeto corine con objeto de facilitar la visualización de los siguientes análisis:

corine[ , 1 ] -> corine
corine

Análisis selección de hábitat: ecotonos

Ahora crearemos la capa de ecotonos, considerando un ancho de banda de 50 metros, con la función buffer. No existe ninguna función en terra que realice directamente el proceso que deseamos, por lo que debemos proceder en varios pasos. Con un valor positivo del argumento width creamos una capa que amplia la geometría del objeto por el exterior, y con un valor negativo por el interior. A continuación “restamos” ambas capas para dejar solo la banda deseada alrededor de las líneas de cada polígono y crear así nuestro mapa ecotonos, al que asignaremos un nuevo código para la variable n1. Agregaremos los polígonos de este mapa en uno solo y eliminaremos las variables innecesarias.

buffer( corine, width = 25 ) -> bout
buffer( corine, width = -25 ) -> bin
bout - bin -> ecotonos
ecotonos$n1 <- 4
aggregate( ecotonos, "n1" ) -> ecotonos
ecotonos
ecotonos[ , 1 ] -> ecotonos
plot( ecotonos )

El siguiente paso es superponer los ecotonos sobre nuestro mapa corine para componer el mapa de usos definitivo. Usaremos la función cover:

cover( corine, ecotonos ) -> corine
corine
plot( corine, "n1", border = 0 )
plot( mcp( hr, 100 ), add = TRUE, lwd = 2, border = 1 : 6 )
plot( rl, add = TRUE )

Finalmente eliminaremos el borde externo del área de estudio, ya que se trata de un ecotono “artificial”:

crop( corine, ext( corine ) - c( 500, 500, 500, 500 ) ) -> corine
plot( corine, "n1", border = 0 )
plot( mcp( hr, 100 ), add = TRUE )
plot( rl, add = TRUE )

Selección de hábitat de 3er orden: radiolocalizaciones

Realizaremos en esta sección un análsis de selección de hábitat de 3er orden (radiolocalizaciones dentro del área de campeo, Jones 2001).

Inicialmete seleccionaremos las radiolocalizaciones que corresponden con actividades de caza (279 en total):

rl[ rl$caza == 1, ] -> rlc
rlc
plot( rlc, add = TRUE, col = "yellow" )

A continuación extraeremos la información del hábitat en cada radiolocalización y la añadiremos como nueva variable:

extract( corine, rlc ) -> cor_rlc
cor_rlc
cor_rlc[ , 2 ] -> rlc$habitat
rlc

Para obtener la tabla de uso para cada individuo:

table( rlc$individuo, rlc$habitat ) -> uso
uso
uso / rowSums( uso ) * 100 -> uso
colnames(uso) <- c( "urbano", "agrícola", "natural", "ecotono" )
uso

Selección de hábitat de 3er orden: áreas de campeo

Continuaremos elaborando una tabla similar a la de la sección anterior, pero ahora para las áreas de campeo. En primer lugar debemos convertir las áreas MCP al formato de mapa vectorial de terra:

vect( mcp( hr, 100 ) ) -> ac
crs( ac ) <- crs( corine )

Ahora extraeremos la información de usos dentro de cada área de campeo. Al tratarse de dos mapas de polígonos es necesario usar intersect y crear un nuevo mapa:

intersect( ac, corine ) -> ac_cor
data.frame( ac_cor )
plot( ac_cor )

La variable area corresponde a la superficie total de cada área de campeo, por lo que la sustituiremos por la superficie correspondiente a cada nuevo polígono:

expanse ( ac_cor, transform = FALSE ) -> ac_cor$area
data.frame( ac_cor )

Comprobemos que la suma de las áreas de estos polígonos para cada área de campeo es igual a la superficie de cada MCP:

aggregate( ac_cor$area, list(ac_cor$id), sum )

La elaboración de la tabla final con estos datos requiere de los siguientes pasos:

xtabs( ac_cor$area ~ ac_cor$id + ac_cor$n1) -> campeo
campeo
campeo / rowSums( campeo ) * 100 -> campeo
colnames(campeo) <- c( "urbano", "agrícola", "natural", "ecotono" )
campeo

Selección de hábitat de 3er orden: análisis composicional

Ya tenemos las dos tablas de porcentajes de tipos de uso por individuo para comparar, lo que haremos mediante un análisis composicional, usando la función compana del paquete adehabitatHS:

compana( uso, campeo )
compana( uso, campeo )$rm

Selección de hábitat de 2º orden

Realizaremos finalmente un análsis de selección de hábitat de 2º orden (áreas de campeo dentro del área de estudio, Jones 2001). Necesitamos estimar los porcentajes de usos del mapa corine en el área de estudio:

expanse( corine )

Dado que esta disponibilidad es la misma para los 6 individuos (las 6 áreas de campeo), construiremos la tabla mediante:

matrix( rep( expanse( corine ), 6 ), 6, 4, byrow = TRUE ) -> total
colnames(total) <- c( "urbano", "agrícola", "natural", "ecotono" )
total

En porcentaje:

total / rowSums( total ) * 100 -> total
total

Por último, realizamos el correspondiente análisis composicional:

compana( campeo, total )
compana( campeo, total )$rm

Evaluación

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

Bibliografía

  • Calenge C (2006) The package adehabitat for the R software: a tool for the analysis of space and habitat use by animals. Ecological Modelling 197: 516-519.

  • Conroy MJ, Carroll JP (2009) Quantitative conservation of vertebrates. Wiley-Blackwell, Oxford.

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

  • Jones J (2001) Habitat selection studies in avian ecology: a critical review. Auk 118: 557-562.

  • Martínez JE, Pagán I, Palazón JA, Calvo JF (2007) Habitat use of booted eagles (Hieraaetus pennatus) in a Special Protection Area: implications for conservation. Biodiversity and Conservation 16: 3481-3488.

Descripción de los mapas y funciones

Utiliza la función info. Por ejemplo: info( "rl" ) o info( "corine90" ).