Visualizar los registros de GBIF con R de una forma rápida

Hace poco ha llegado a mis manos (vía theBioBucket, de Kay Cichini) un aporte sobre la utilización de R en la gestión de registros de biodiversidad, que me parece sencillamente fantástico. Se trata de la consulta de los registros de biodiversidad de GBIF mediante R. 

Imaginad, que estamos trabajando con una especie y necesitamos consultar los datos existentes sobre esta especie en la plataforma GBIF (Global Biodiversity Information Facility). Pues, podemos utilizar, por ejemplo...R. Sí, existe un paquete llamado dismo, creado por Robert Hijmans, Steven Phillips, John Leathwick y jane Elith, que realiza diversas tareas relacionadas con los modelos de distribución de especies.


Pues, vamos a ver un ejemplo de alguna de las cosas que podemos hacer. Hace poco, leí mi Trabajo de Fin de Master, que versaba sobre los melojares de Quercus pyrenaica en Sierra Nevada (hablaré de ello en otra entrada) y voy a utilizar esa especie como consulta. 

En primer lugar, hemos de instalar y cargar el paquete dismo. Igualmente vamos a cargar otros paquetes necesarios para visualizar los resultados en mapas. 

Una vez cargado, utilizamos la función gbif() para que nos consulte el portal de GBIF y nos devuelva los resultados de la especie que deseemos. Para ello hemos de especificar el nombre del género (en este caso "Quercus") y el epiteto de la especie ("pyrenaica"). Además le indicamos que nos devuelva los datos geográficos, para lo cual fijamos el parámetro geo en TRUE. 

En este caso nos devuelve un dataframe con 667 resultados de ocurrencia de esta especie. En este objeto, observamos que nos devuelve información interesante como la institución responsable de la cita, el número de identificación de la cita, que tipo de registro es: si es observación o cita, etc. 

Esta sencilla consulta, en una sola línea de código, nos agiliza muchísimo la consulta en la base de datos de GBIF y nos permite llevar a cabo consultas de muchas especies en el caso en el que estemos llevando a cabo trabajos donde necesitamos recopilar información de citas de muchas especies. 

Pero, ademas de esto (que no es poco), podemos con dos o tres líneas mas de código, realizar mapas de los resultados de esta consulta. Para ello necesitamos otros packages. En este caso os dejo dos mapas de la distribución de las citas para Quercus pyrenaica, tanto en la Península Ibérica, como en Sierra Nevada. En este último caso, es para resaltar que también podemos utilizarlo incorporando shapefiles de ESRI. 

Figura 1. Ejemplo de consulta de los registros de Quercus pyrenaica en la Península Ibérica y su representación con R. 

Figura 2. Ejemplo de consulta de los registros de Quercus pyrenaica en Sierra Nevada. En este caso se ha importado un shapefile con los límites del Espacio Natural de Sierra Nevada. 

Pues nada más, solamente el código utilizado: 

######### Consulta y visualización de datos de GBIF 
 
# Cargamos la libreria dismo para la consulta a GBIF
library(dismo) 
 
# Consultamos los datos de gbif de la especie que deseemos
melojo <- gbif("Quercus", "pyrenaica", geo=TRUE)
 
# Cargamos la libreria para plotear en un mapa 
library(maps)
library(maptools)
library(mapdata)
 
# Ploteamos los resultados para la Península Ibérica
m <- map('worldHires', regions=c('Spain', 'Portugal'), exact=T, interior=FALSE, col="grey") 
area.map(m) 
# Aquí añadimos los resultados de nuestra consulta a GBIF en el mapa
points(melojo$lon, melojo$lat, cex=1.5, pch=19, col="blue") 
# Ahora la leyenda
text(-.001,36.3, "Resultados consulta GBIF\nQuercus pyrenaica")
 
 
 
# Ahora lo mismo pero con datos para Sierra Nevada
# Importamos un shape del shapefile del límite del Espacio Natural de Sierra Nevada
sn_lim <- readShapePoly("C:/consulta_datos_R_gbif/EENNPP/InfGeografica/InfVectorial/Shapes/EspacioNaturalSierraNevada.shp"
, proj4string=CRS("+proj=latlong +ellps=WGS84"))
 
# Ploteamos los límites de Sierra Nevada
plot(sn_lim, border="grey", axes=TRUE)
points(melojo$lon, melojo$lat, cex=1.3, pch=19, col="blue") 
text(-3.1, 37.4, "Consulta de datos de Quercus pyrenaica\nen GBIF mediante R", cex=1.5)
text(-3.1, 36.4, "Espacio Natural de Sierra Nevada", cex=1.5)

Created by Pretty R at inside-R.org

Dos gráficos en uno. R

Retomo mi olvidado blog para escribir algunas notas y para que no se me olvide como procedí cuando quería realizar un gráfico como el que presento a continuación.

Una de las cosas que mas me gusta de R es que primero uno piensa el gráfico que necesita, y luego lo hace, es decir, que los gráficos no están tan cerrados como en otros softwares estadísticos. En este caso la idea era realizar un gráfico de densidad de dos conjuntos de datos y añadirle un boxplot que coincidiera en el eje x. Esta idea es una modificación de la que en su día tuvo el Dr. Blas Benito para representar requerimientos ecológicos de especies vegetales y que se puede encontrar aquí.

Para comprenderlo mejor vamos a entrar en faena. Supongamos que tenemos un dataset con una variable binaria (presencia/ausencia) de una especie vegetal y muchas otras variables ambientales de dichos puntos.

Como dijimos nuestro propósito es representar gráficos de densidad de los datos de presencia y ausencia respecto a una variable ambiental concreta, y añadirle después un boxplot debajo, que coincida en el eje x. Para ello en primer lugar creamos dos subsets de datos (uno para los datos de presencia y otro para los datos de ausencia) (existen otras formas de hacerlo, pero yo he procedido de este modo):

# Separo los datos de presencia de los no presencia
datos0 <- datos[datos$presencia==0,]
datos1 <- datos[datos$presencia==1,] 

A continuación calculamos el gráfico de densidad para la variable ambiental, en este caso, la Temperatura mínima del invierno para el año 2000.

# Calculamos la Densidad para la variable temperatura mínima 
d.datos0 <- density(datos0$tmniob2000)
d.datos0 <- density(datos1$tmniob2000)

Seguidamente para asegurarnos de que todos los gráficos de todas las variables ambientales que queremos hacer guarden la misma proporción entre el plot de densidad y el boxplot, vamos a obtener los límites para los ejes x e y.

# Obtener los valores máximos de los ejes y el mínimo del eje x
ymax <- max(c(d.datos0$y, d.datos1$y))
xmin <- min(c(d.datos0$x, d.datos1$x))
xmax <- max(c(d.datos0$x, d.datos1$x))

El valor ymax representa el valor positivo máximo del eje y. Esta parte del gráfico (es decir, la parte positiva) queremos que siempre suponga el 70 % (por ejemplo) del gráfico combinado (boxplot y density). La parte negativa del eje y, será donde ubiquemos los boxplots y queremos que siempre ocupen un 30 % de la proporción del gráfico combinado. Para ello lo que vamos a hacer es crear una variable llamada ymin que sea el 30 % de la región del gráfico, teniendo en cuenta para ello el valor de ymax (El fijar estas proporciones es importante porque cuando realizamos este gráfico para diferentes variables ambientales, puede ser que los boxplots no tengan el mismo tamaño).

# Establecer un mínimo en el eje y, que sea el 30 % de la región del gráfico
ymin <- (ymax*(10/7))*0.3

Una vez que tenemos fijados dichos valores, vamos a empezar a realizar el gráfico. En primer lugar dibujamos la densidad de los datos de ausencias. Seguidamente añadimos la densidad de datos de presencia.

plot(d.datos0, ylim=c(-ymin,ymax), main="", ylab="", xlab="Temperatura mínima ivierno (ºC*10)")
polygon(d.datos1, col="#669700", border="#669700")
lines(d.datos0,xlim=c(0, max(datos0$tmxvob2000)), ylim=c(0,ymax), col="#336699", lwd=2) 

Tras hacerlo nos queda una cosa similar a esta

Si os fijáis nos queda un espacio vacío en la zona negativa del eje y, que siempre será de la misma proporción para todos los gráficos de todas las variables ambientales. Allí es donde vamos a realizar el gráfico de boxplots. Para ello vamos a utilizar una función llamada subplots(), del package TeachingDemos. Su funcionamiento es muy sencillo y como argumentos (ver ayuda de la función) establecemos qué gráfico vamos a incluir y en qué lugar. Para eso último utilizaremos las coordenadas de la región gráfica, haciendo uso del argumento 'usr' de la función par (). Este argumento es un vector de 4 elementos dando las coordenadas de los extremos de la región gráfica, en el siguiente modo c(x1, x2, y1, y2). Para el caso de las coordenadas del eje x, utilizamos par('usr')[1:2], y para la región negativa del eje y, lo que hacemos es dividirla en dos para que cuadren los dos boxplot. Todo ello lo hacemos así:

subplot(boxplot(datos1$tmxvob2000, horizontal=T, ylim=c(xmin,xmax), axes=F, col="#669700"),par('usr')[1:2], c((par('usr')[3] - 0)/2,0))
subplot(boxplot(datos0$tmxvob2000, horizontal=T, ylim=c(xmin,xmax), axes=F, col="#336699"),par('usr')[1:2], c(par('usr')[3],(par('usr')[3] - 0)/2))

Y conseguimos el gráfico combinado que nos proponíamos al inicio. Como se puede apreciar los boxplots coincíden en el eje x con los gráficos de densidad.

Pues nada más por ahora, espero os sirva.

Pd: Gracias al Dr. Benito por los datos.

Pd2: Para presentar el código he utilizado Pretty-R

[GeoSTA] Semivariogramas teóricos con R

Empezando a aprender conceptos básicos de geoestadística, estoy lidiando con algunos de ellos duros inicialmente, pero interesantes e importantes para entender bien la modelización espacio-temporal. Así que en estas tareas, creo que he asimilado el concepto de variograma y semivariograma (pronto publicare unos pequeños apuntes sobre conceptos básicos); y he tenido que realizar un resumen de los principales modelos teóricos de semivariogramas. La cuestión es que me apetecía hacer unas gráficas básicas sobre los modelos teóricos y en primer lugar pense en hacerlos con power point, pero luego, pense que mejor los hacía con R, y así de camino aprendía algo más.

Para ello he utilizado el paquete geoR sobre geoestadística y he conseguido representar mediante una simulación tres modelos teóricos del semivariograma: exponencial, gausiano y esférico). A continuación os dejo la gráfica que he producido y el código que he utilizado para ello.

La figura muestra la semivarianza (eje y) frente a la distancia (eje x). En el siguiente código se explica como se ha realizado tal gráfico:

plot(0:2, 0:2, type="n",axes=F, xlab="|h|", ylab=expression(gamma(h)), font.lab=2)
axis(1, at=c(0,6), labels=F)
axis(2, at=c(-.5,6), labels=F, pos=0)
lines.variomodel(cov.m="sph", cov.p =c(1, .8), nug=.06, max.dist=3, lty=1, scaled=T, col="black", lwd=2)
lines.variomodel(cov.m="exp", cov.p =c(1, .8), nug=.06, max.dist=3, lty=1, scaled=T, col="red", lwd=2)
lines.variomodel(cov.m="gau", cov.p =c(1, .8), nug=.06, max.dist=3, lty=1, scaled=T, col="blue", lwd=2)
legend("top",
expression(list("Esférico"),list("Exponencial"), list("Gaussiano")),
lty=c(1,1,1), lwd=c(2,2,2), col=c("black","red","blue"))

Los parámetros mas importantes de lines.variomodel son:

• cov.m que nos indica el modelo teórico, como por ejemplo: esférico ("sph"), exponencial ("exp"), gaussiano ("gau"), etc.
• cov.p nos indica los parámetros meseta o sill y el rango
• nug nos permite especificar el efecto "pepita" (nugget) (vamos, la separación que observarmos en el eje de ordenadas respecto del origen)

Pues nada mas por ahora, pronto vendrán los apuntes de los conceptos básicos en geoestadística.

Escribir código R en tu blog... y que quede bonico.

Pues sí, mucha gente utiliza R para llevar a cabo sus análisis, y casi todos coinciden en la potencialidad de esta herramienta. Son muchos los investigadores que lo han incorporado como una herramienta diaria en su trabajo, y otros muchos, los que además de hacer increíbles cosas en R, lo publican en sus websites o en sus blogs (ej.: L. Cayuela, B. Benito, por citar gente cercana).

Como neo-neofito en esto de R, me ha sorprendido una simple aplicación que te permite poner bonico el código que publicas en tu blog, es decir, realizas tu script para explicar cualquier análisis y luego lo publicas en tu blog, y además, queda bonito. Pues si, es posible. Esta herramienta se llama Pretty R y transforma tu código R en código html que puedes fácilmente publicar en tu blog.

Así que ya tenemos otra excusa menos para seguir poniendo a disposición de la sociedad nuestros avances y las metodologías utilizadas para conseguirlos.

Baste un ejemplo tonto:

Imaginemos que tenemos que realizar un análisis cluster, lo hacemos con R, y queremos publicar nuestro código. La diferencia entre usar o no esta herramienta la podemos ver a continuación:

A. Código si usar Pretty R:

data <- read.table("datos_variables.txt", header=T)

d1 <- dis(data, method="euclidean")

cluster1 <- hclust(d, method = "complete")

B. Código usando Pretty R:

data <- read.table("datos_variables.txt", header=T)
d1 <- dis(data, method="euclidean")
cluster1 <- hclust(d, method = "complete") 

Pues nada, juzguen y decidan.

Cosas diferentes con R.

Buscando herramientas de la llamada web 2.0 aplicadas a la ciencia (Ciencia 2.0) hemos encontrado una manera elegante de acceder a los contenidos de cualquier sitio web mediante la creación de un archivo html con una nubes de etiquetas en Flash. Hasta ahí todo normal, pero lo interesante es que se ha hecho con el programa R, lo cual nos convence mas de la gran capacidad y versatilidad del mismo. Gracias al trabajo de Y. Xie, T. Wei y Y. Qui, que han desarrollado un paquete (fun) para R que permite realizar esta y otras acciones.

Para hacer esta nube de etiquetas podemos seguir las instrucciones que muestra Yihui Xie en su blog, que se resumen en unos pasos muy sencillos:

1. Cargar en R el paquete fun

2. Tener preparado un archivo de datos con las etiquetas de los enlaces que queremos. Este archivo ha de tener al menos 3 columnas: tag (nombre de la etiqueta), link (enlace a la etiqueta) y count (numero de veces que aparece la etiqueta). Además se pueden añadir dos columnas mas: color y hicolor, que son el color que mostraran las etiquetas y el color que tomarán cuando las señalemos, respectivamente.

3. Ejecutar la función tagCloud ( ) con los parámetros que nosotros definamos: nombre del archivo de salida, dimensiones, color de fondo, etc. En este sentido es importante atender a la opción target nos permite decidir si cuando pinchamos sobre una etiqueta se abre en una pestaña nueva del explorador (para ello especificar target="_blank"). No obstante para conocer todas las opciones de la nube de etiquetas a crear, podemos consultar la ayuda de la función: ?tagCloud

4. Una vez creado el archivo html hemos he abrirlo con un editor (ej.: Notepad++) y añadirle la siguiente línea, en cualquier parte de la etiqueta script, es decir entre <*script> y <*/script>.
so.addParam("allowScriptAccess", "always");

Nosotros lo hemos utilizado para crear una puerta de acceso mas amena a los diferentes ámbitos temáticos del programa de seguimiento del Observatorio de Cambio Global en Sierra Nevada, y aunque no es definitvo, puede ser una forma elegante y diferente de acceder a los contenidos. Para ver el ejemplo seguir este enlace.

Os dejamos el código utilizado y el archivo de datos. Para mas información sobre el funcionamiento del paquete fun, consultar:

Yihui Xie, Taiyun Wei and Yixuan Qiu (2009). fun: Use R for Fun. R package version 0.1-0/r14. http://R-Forge.R-project.org/projects/fun/

Riego tras la plantación, grafica con R.

Preparándo un poster sobre los trabajos de restauración de una planta amenazada (Arenaria delaguardiae), que vamos a presentar al próximo congreso de Biología de la Conservación de Plantas, organizado por la Sociedad Española de Biología de la Conservación de Plantas y que se celebrará en Septiembre en Almería; se nos ha ocurrido una forma gráfica de mostrar los días en los que se realizaron los riegos tras la plantación de ejemplares de esta especie en campo.

En las restauraciones de la cubierta vegetal, después de una plantación se aconseja realizar riegos de apoyo que minimicen la pérdida de ejemplares y se disminuyan las marras, pero encontrando un equilibrio que evite que las plantas introducidas dependan de éstos riegos, y sobreviva sin necesidad de apoyo continuado (que es el objetivo). La cantidad de riegos a realizar depende de muchos factores: el microclima de la zona, la orientación, la pendiente, la comunidad vegetal, el trascurso del año hidrológico, etc.

En este caso concreto se atendieron a algunos de esos factores, especialmente la climatología del año (2008) y se realizaron riegos en los cuales se disminuyó progresivamente la cantidad de agua aportada por planta. En la siguiente gráfica se muestran los días en los que se realizaron los riegos, además se acompañan de las medias mensuales del periodo 2001-2008 (con la desviación estandar) y la precipitación del año 2008.


Hemos querido trabajar con el programa R. Gracias especialmente al empeño y ánimo de L. Cayuela, vamos poco a poco manejándonos con el programa y realizando los análisis estadísticos y las gráficas que mas nos interesan. En este sentido vale la pena visitar su blog en el que se aportan unos manuales bastante sencillos y claros para trabajar con R.

Esta gráfica es una primera aproximación. Queremos añadir datos históricos con mayor rango (mas años) y de mas estaciones (se ha realizado con datos de la estación mas cercana: estación agroclimática de Padul). No obstante nos parece una forma diferente de expresar los criterios seguidos a la hora de la realización de los riegos.

Os adjunto el código en R que he utilizado para crear la gráfica asi como los datos utilizados.