Luis Verde Arregoitia bio photo

Luis Verde Arregoitia

Personal research page - ecology . evolution . conservation . biogeography

Twitter Publons ResearchGate Google Scholar

Extraer valores de uno o más raster a partir de coordenadas xy usando R.

En temas de biodiversidad, es común trabajar con registros puntuales de la distribución de una o más especies, generalmente representadas en coordenadas geográficas xy (longitud/latitud). Para ciertos análisis en macroecología, conservación, y biogeografía, necesitamos datos sobre las condiciones ambientales en cada registro. Estos datos de clima, vegetación, uso de suelo, etc. generalmente existen como capas raster (una representación del área de estudio en formato de matriz dividida en celdas con valores únicos).

Extraer el valor de las celdas en donde ocurren nuestros registros puntuales es fácil usando sistemas de información geográfica (SIG) con interfaces gráficas, que casi siempre tienen algun botón o herramienta dedicada a esta tarea. Para uno de mis proyectos, tuve que extraer los valores de varios raster para los registros puntuales de algunas especies de murciélagos en México. En algún momento decidí que mi regreso a los SIG iba a ser exclusivamente en R, y por eso investigué cómo logar esta extracción de datos sin programas adicionales.

En este caso usé los raster de temperatura y precipitación del proyecto CHELSA Climate. Estos datos se encuentran disponibles desde hace un par de semanas, y parece que tienen mejor exactitud en zonas montañosas que otras opciones existentes como WorldClim. Aquí se pueden descargar los datos de clima en formato geotiff. Son archivos grandes, pero yo no tuve problema para manipularlos en una laptop que ya tiene 5 años.

En este ejemplo, obtuve simultáneamente los valores de dos raster para 150 localidades: la temperatura promedio en enero y la precipitacion promedio en julio. Para reproducirlo, hay que tener los archivos de clima que son Mean january temperature y Mean july precipitation en esta página de descargas, y las coordenadas se pueden leer directamente de este archivo csv.

# cargar paquetes, instalarlos si hace falta usando install.packages()
library(raster)
library(dplyr)
library(sp)
library(rgdal)

# leer los raster de temperatura, en este caso los archivos estan en el directorio de trabajo, que se puede 
   # averiguar usando getwd()
janTemp <- raster("CHELSA_temp_1_1979-2013.tif") #temperaturas enero
julPrec <- raster("CHELSA_prec_7_1979-2013.tif") # precipitaciones julio
# apilarlos como ingredientes en una hamburguesa
climStack <- stack(janTemp,julPrec)
# leer el archivo con las localidades 
localidades <- read.csv("https://raw.githubusercontent.com/luisDVA/codeluis/master/localidades.csv")%>% select(long=2,lat=3)

Podemos visualizar las diferentes capas (llamándolas con índices), en este caso usando una extensión geográfica definida con extent y encima de esto se pueden graficar los puntos.

# extension geografica para graficar
	# definida por límites máximos y mínimos de long y lat 
boxextent <- extent(-115,-90,10,25)

# dibujar el raster y los puntos (en este caso el primer elemento de la pila)
plot(climStack$CHELSA_temp_1_1979.2013,ext=boxextent)
points(localidades) #añade los puntos sobre la gráfica existente
utilizando puros valores gráficos predeterminados

Nota: los conjuntos de capas se pueden recortar con rasterVis::crop por cuestiones de memoria.

La extracción es fácil, y se hace para todos los elementos de la pila al mismo tiempo.

# extrer valores de todas las capas en la pila
extractedClim <- extract(climStack,localidades)
# juntar valores con las coordenadas de las localidades
bioclimat<-cbind(localidades,extractedClim)

# escrbir la tabla en la memoria
write.csv(bioclimat,"biolcimatic.csv")

Al final nos queda algo que se ve más o menos así:

longitude lat CHELSA_temp_1_1979-2013 CHELSA_prec_7_1979-2013
-96.385 19.60055556 -99.7333 17.33361
-96.385 19.60055556 25.20189667 1548.687744
-99.5689 21.19917 19.83035 940.9037
-96.385 19.60055556 25.20189667 1548.687744
-98.7333 21.1 17.33361 545.4098

Como tenemos dos variables climáticas, las podemos graficar en 2d como un perfil bioclimático. Este ejemplo es con ggplot pero hay otras opciones para generar gráficas.

library(ggplot2)


ggplot(bioclimat,aes(x=CHELSA_temp_1_1979.2013,y=CHELSA_prec_7_1979.2013))+
  geom_point()+geom_density2d()+
  xlab("temperatura")+ylab("precipitación")+
  theme_minimal()
envoltura bioclimática

!Listo! Sabiendo extraer datos, los más tardado será la descarga de los raster.

No duden en contactarme si tienen cualquier duda o sugerencia.