Métodos de análisis de datos espaciales

Práctica 6, 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

23 de noviembre de 2023

Fecha de modificación

4 de marzo de 2024

Introducción

En esta práctica aplicaremos diferentes técnicas de análisis espacial. Comenzaremos por métodos que permiten determinar el tipo de distribución espacial de puntos (función $K$) y continuaremos con la elaboración de correlogramas y la aplicación de técnicas de interpolación (kriging), que permiten abordar aspectos relacionados con la autocorrelación espacial. Finalmente aplicaremos procedimientos de análisis multivariante (ordenación y clasificación) a datos espaciales.

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( "bats" ) o info( "no2" ). Teclea info() para más información.

Para realizar esta práctica se necesita la instalación adicional de los paquetes de R terra, spatstat, spdep, pgirmess, gstat y geodata. Si no los tenemos instalados usaremos:

install.packages( c( "terra", "spatstat", "spdep", "pgirmess", "gstat", "geodata" ) )

Una vez instalados los cargaremos en memoria:

library( terra )
library( spatstat )
library( spdep )
library( pgirmess )
library( gstat )
library( geodata )

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

unwrap( murcia ) -> murcia
unwrap( bats ) -> bats
unwrap( cnv ) -> cnv
unwrap( mdt200 ) -> mdt200
project( murcia, crs( cnv ) ) -> murcia
ed50( bats ) -> bats

Análisis espacial de puntos: función K

Existen varias técnicas de análisis de la distribución espacial de localizaciones (mapas de puntos). Nos centraremos aquí en la aplicación de la denominada función $K$ de Ripley (Ripley’s $K$ function), utilizando funciones del paquete de R spatstat. Como datos de ejemplo utilizaremos las localizaciones de refugios de murciélagos en la Región de Murcia (objeto bats):

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

Analizaremos inicialmente la distribución del murciélago de cueva (Miniopterus schreibersii), para lo cual crearemos un nuevo mapa solo con las localizaciones de los refugios en los que está presente:

bats[ bats$msch > 0, ] -> msch
plot( murcia )
plot( msch, add = TRUE )

Necesitamos extraer la información de las coordenadas (con la función geom):

geom( msch )
xy <- geom( msch )[ , 3 : 4 ]
xy

A continuación las transformaremos al formato del paquete spatstat (con la función ppp), especificando una “ventana” espacial para el análisis. Finalmente aplicaremos la función envelope para realizar el análisis:

pmsch <- ppp( xy[ , 1 ], xy[ , 2 ], c( min( xy[ , 1 ] ), max( xy[ , 1 ] ) ),  c( min( xy[ , 2 ] ), max( xy[ , 2 ] ) ))
plot( envelope( pmsch ) )

Ejercicio 1

Realiza el ejercicio anterior con el resto de especies de murciélago (rfer, rhip y mcap).

Ver la solución

bats[ bats$rfer > 0, ] -> rfer
bats[ bats$rhip > 0, ] -> rhip
bats[ bats$mcap > 0, ] -> mcap

xy <- geom( rfer )[ , 3 : 4 ]
prfer <- ppp( xy[ , 1 ], xy[ , 2 ], c( min( xy[ , 1 ] ), max( xy[ , 1 ] ) ),  c( min( xy[ , 2 ] ), max( xy[ , 2 ] ) ))
xy <- geom( mcap )[ , 3 : 4 ]
prhip <- ppp( xy[ , 1 ], xy[ , 2 ], c( min( xy[ , 1 ] ), max( xy[ , 1 ] ) ),  c( min( xy[ , 2 ] ), max( xy[ , 2 ] ) ))
xy <- geom( msch )[ , 3 : 4 ]
pmcap <- ppp( xy[ , 1 ], xy[ , 2 ], c( min( xy[ , 1 ] ), max( xy[ , 1 ] ) ),  c( min( xy[ , 2 ] ), max( xy[ , 2 ] ) ))

par( mfrow = c( 1, 3 ) )
plot( envelope ( prfer ) )
plot( envelope ( prhip ) )
plot( envelope ( pmcap ) )

Autocorrelación y geoestadística

Una de las principales características de los datos espaciales está relacionada con procesos de autocorrelación, es decir, con la circunstancia de que los elementos u objetos cercanos tienden a parecerse más entre sí, que los más alejados. Nos referiremos aquí, por tanto, a información espacial de naturaleza cuantitativa.

El análisis estadístico de los procesos de autocorrrelación es complejo. Aquí nos centraremos en la elaboración de correlogramas con el índice $I$ de Moran y en la aplicación de técnicas de interpolación basasadas en el uso de semivariogramas (kriging).

Correlogramas: índice de Moran

Trabajaremos con el ejemplo del mapa de carnívoros (objeto cnv), aplicando el análisis a los datos de riqueza de especies por cuadrícula.

plot( cnv )

Comenzaremos elaborando una nueva capa de riqueza, de forma que cada cuadrícula tenga un valor (entre 0 y 7):

subst( cnv, NA, 0 ) -> cnv
crop( cnv, murcia, mask = TRUE, snap = "out" ) -> cnv
sum( cnv ) -> cnv$riqueza
plot( cnv )

Para hacer el análisis hay que extraer las coordenadas de las cuadrículas por un lado (crds( cnv )), y los valores de riqueza por otro (data.frame(cnv$riqueza)[,1]). Aplicaremos la función correlog del paquete pgirmess (que recurre a la función moran.test del paquete spdep):

par( mfrow = c( 1, 2 ) )
plot( cnv$riqueza )
correlog( crds( cnv ), data.frame(cnv$riqueza)[,1] )
plot( correlog( crds( cnv ), data.frame(cnv$riqueza)[,1] ) )

Ejercicio 2

Realiza un correlograma con la riqueza de especies de murciélagos en los refugios de la Región de Murcia.

Ver la solución

rowSums( data.frame( bats[ , 2 : 5 ] ) ) -> bats$riqueza
par( mfrow = c( 1, 2 ) )
plot( murcia )
plot( bats, "riqueza", add = TRUE )
correlog( geom( bats )[ , 3 : 4 ], bats$riqueza )
plot( correlog( geom( bats )[ , 3 : 4 ], bats$riqueza ) )

Interpolación: kriging

Las técnicas de krigeado, o kriging son una herramienta fundamental en el campo de la geoestadística, desarrolladas para predecir los valores de una variable en el espacio. Su aplicación se fundamenta en la estimación de semivariogramas (o, simplificando, variogramas). En esta sección realizaremos dos ejercicios con variables de diferente naturaleza, utilizando funciones del paquete gstat.

Para el primer ejercicio usaremos los datos de concentraciones medias anuales de NO₂ atmosférico (expresadas en $\mu g /m^3$) medidas en 8 estaciones de la Región de Murcia.

no2

El primer paso es obtener el (semi)variograma. Utilizaremos los datos de NO₂ transformados logarítmicamente:

v <- variogram( log( NO2 ) ~ 1, ~ x + y, data = no2  )
plot( v )

El siguiente paso es ajustar el variograma obtenido a un modelo matemático. Existe una amplia variedad de modelos de variogramas para ajustarlo:

show.vgms()

En este caso elegiremos el modelo exponencial ("Exp"):

mv <- fit.variogram( v, vgm( "Exp" ) )
kmv <- gstat( NULL, "logNO2", log( NO2 ) ~ 1, ~ x + y, data = no2, model = mv )

Finalmente realizaremos el krigeado con la función interpolate del paquete terra (es necesario especificar un mapa raster de referencia):

kno2 <- interpolate( mdt200, kmv )
plot( kno2 )

Para una mejor interpretación recortaremos el mapa de predicción con los limites regionales y superpondremos las localizaciones de las estaciones de muestreo:

plot( mask( kno2[[ 1 ]], murcia ) )
points( no2[ , 2 : 3 ] )

El segundo ejercicio lo realizaremos con el mapa de riqueza de especies de murciélago en refugios de la Región de Murcia. Necesitamos extraer la información (coordenadas de los refugios y valores de riqueza) en una tabla de datos:

rb <- data.frame( riqueza = bats$riqueza, x = geom( bats )[ , 3 ], y = geom( bats )[ , 4 ] )
v <- variogram( riqueza ~ 1, ~ x + y, data = rb  )
plot( v )

A continuación realizaremos los siguientes pasos, ajustando en este caso el variograma a un modelo esférico ("Sph"):

mv <- fit.variogram( v, vgm( "Sph" ) )
kmv <- gstat(NULL, "riqueza", riqueza ~ 1, ~ x + y, data = rb, model = mv )
krb <- interpolate( mdt200, kmv )
plot( mask( krb[[ 1 ]], murcia ) )
plot( bats, "riqueza", add = TRUE, legend = "topleft" )

Ejercicio 3

Realiza un krigeado con la riqueza de mamíferos carnívoros.

Ver la solución

rc <- data.frame( riqueza = cnv$riqueza, x = crds( cnv )[ , 1 ], y = crds( cnv )[ , 2 ] )
v <- variogram( riqueza ~ 1, ~ x + y, data = rc  )
plot( v )

mv <- fit.variogram( v, vgm( "Sph" ) )
kmv <- gstat(NULL, "riqueza", riqueza ~ 1, ~ x + y, data = rc, model = mv )
krc <- interpolate( aggregate(mdt200, 25), kmv )
plot( mask( krc[[ 1 ]], murcia ) )

Análisis multivariante

Las características de las técnicas multivariantes (ordenación y clasificación) las hacen apropiadas para su aplicación a datos espaciales. Un procedimiento habitual es realizar una ordenación de un conjunto amplio de capas de información ambiental, y aplicar una clasificación sobre los resultados de la ordenación. Seguiremos este procedimiento en el ejercicio de este apartado, para el que usaremos como ejemplo el mapa de variables bioclimáticas de WorldClim, que tendremos que descargar y recortar previamente:

worldclim_country( country = "Spain" , var = "bio", path = "./descargas/" ) -> spain_bio
crop( project( spain_bio, crs( murcia ) ), murcia, mask = TRUE ) -> murcia_bio
plot( murcia_bio )

A continuación realizaremos una ordenación mediante un análisis de componentes principales (principal component analysis) con la función de R prcomp. Los resultados (las coordenadas de cada cuadrícula en cada componente o eje del análisis) se pueden aplicar sobre el mapa con la función predict del paquete terra:

prcomp( data.frame( murcia_bio ) ) -> pca
predict( murcia_bio, pca ) -> rpca
rpca
plot( rpca )

Finalmente usaremos el método de clasificación de K-medias (K-means) aplicado sobre un mapa ráster. Usaremos la función ckmeans, disponible en el archivo deb.RData, que recurre a la función kmeans de R. Especificaremos el número de grupos (clústers) deseados:

ckmeans( rpca, grupos = 5 ) -> rkm
plot( rkm )

Evaluación

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

Bibliografía

Burrough PA, McDonell RA, Lloyd CD (2015) Principles of Geographical Information Systems. 3ª ed. Oxford University Press.
Dale MRT, Fortin M-J (2014) Spatial Analysis. A Guide For Ecologists. 2ª ed. Cambridge University Press.
Hijmans R (2023) _terra: Spatial Data Analysis_. R package version 1.7-55, https://CRAN.R-project.org/package=terra
Moraga P (2023) Spatial Statistics for Data Science: Theory and Practice with R. Chapman & Hall/CRC Data Science Series. https://www.paulamoraga.com/book-spatial/

Descripción de los mapas y funciones

Utiliza la función info. Por ejemplo: info( "bats" ) o info( "no2" ).