Como utilizar el potencial de GRASS desde R

R tiene varios paquetes que pueden utilizarse para el análisis y visualización de datos espaciales (sp, maptools, rgdal, etc.). Sin embargo, R no es un SIG y, por tanto, tiene ciertas limitaciones a la hora de trabajar con información espacial. Por ejemplo, para algo tan simple como extraer la pendiente y la orientación a partir de un modelo digital de elevaciones, R necesita del paquete RSAGA, que a su vez necesita del software gratuito SAGA diseñado por Alexander Brenning, que funciona únicamente bajo Windows. Esto supone ciertas limitaciones al uso de R como SIG. Estas limitaciones no existen cuando se utiliza GRASS en combinación con R. GRASS ofrece un complemento perfecto a R para la implementación de operaciones típicamente realizadas en SIG y R permite potenciar las capacidades de GRASS para la implementación de análisis geoespaciales.

Para poder hacer uso de ambos programas en necesario tener instalados GRASS y una versión de R posterior a 2.1.0. Desde el 'shell' se puede abrir una sesión de R. Lo primero será instalar los paquetes necesarios para poder utilizar R en combinación con GRASS (sp, rgdal, maptools, spgrass6, spGDAL, spmaptools).

GRASS 6.3.0 (Filabres):~ > R

install.packages(c("sp", "rgdal", "maptools"), dependencies= TRUE)

rS <- "http://r-spatial.sourceforge.net/R" install.packages(c("spgrass6", "spGDAL", "spmaptools"), repos=rS, dependencies=TRUE)

Sin salir de la sesión de R, debemos de cargar el paquete spgrass6 y la 'location' de GRASS donde se encuentran todas las capas de información raster o vectorial.

library(spgrass6)
G <- gmeta6()
str(G)

Si ahora queremos acceder a alguna de las capas raster del 'location' en el que estamos trabajando dentro de GRASS tendremos que utilizar la función de R readCELL6sp() ò readRAST6(). Esto genera un objeto del tipo SpatialGridDataFrame.

mde <- readCELL6sp("mde")

Sin embargo, es importante tener en cuenta que si cargamos muchas capas simultáneamente, especialmente si son capas con un gran número de celdas, se puede agotar la memoria virtual. En este caso R dará el mensaje 'vector memory exhausted'. Cuando esto me ha ocurrido a mi, R se ha quedado colgado y he tenido que reiniciar la sesión. Así que lo que aconsejo es ir eliminando los objetos que ya no necesitemos una vez que extraigamos de ellos la información necesaria. Para conocer el uso que se va haciendo de la memoria se puede utilizar la función gc().

gc()

used (Mb) gc trigger (Mb) max used (Mb)
Ncells 279590 7.5 531268 14.2 361111 9.7
Vcells 3160949 24.2 15608917 119.1 30379099 231.8

Es decir que se están utilizando alrededor de 31 MB de memoria virtual (hay que mirar la primera columna). Ahora eliminamos la capa mde.

rm(mde)
gc()

used (Mb) gc trigger (Mb) max used (Mb)
Ncells 279513 7.5 531268 14.2 361111 9.7
Vcells 137524 1.1 12487133 95.3 30379099 231.8

Y comprobamos que disminuye considerablemente el uso de la memoria virtual.

Antes de acceder a las capas raster, vamos a leer una capa de puntos en formato *.txt (parcelas muestreadas en campo) para la cuál extraeremos la información de elevación, pendiente, altitud e insolación de las capas raster.

points <- read.table(url("http://archivos-para-subir.googlegroups.com/web/Analisis+correlacionales.txt?gda=0ywgME4AAACkOfwqHAVd4YqgfIB09GDRN3No94DJejqv6LWrffOcUtQTec6xGbbz9B_IWZmp37itXp5Ud9d8afGj09bAbTSQ47Cl1bPl-23V2XOW7kn5sQ"), header = T, sep = "\t", dec = ",")
str(points)

'data.frame': 76 obs. of 7 variables:
$ x : num 543726 543486 542903 535211 539086 ...
$ y : num 4124164 4124224 4122417 4118071 4121908 ...
$ Id : Factor w/ 76 levels "1_12","1_13",..: 1 2 4 6 7 8 9 10 11 12 ...
$ Especie : Factor w/ 2 levels "Pinus nigra",..: 1 1 1 1 1 1 1 2 1 1 ...
$ Defoliacion: num 55 55 30 55 30 55 55 60 45 55 ...
$ Area_basim : num 25.8 28.1 24.2 22.0 18.0 ...
$ Densidad : num 286 286 1592 923 1974 ...

El siguiente paso es convertir esta capa de puntos en un objeto espacial del tipo SpatialPointsDataFrame y asignarle una proyección (ha de ser la misma que la de las capas raster y sino habría que reproyectar esta capa).

coordinates(points) <- data.frame(points$x, points$y)
proj4string(points) <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")

Y ahora sí, importamos cada una de las capas raster de GRASS a R, extraemos los valores de las capas raster para los puntos incluidos en la capa vectorial points y, por último, borramos las capas raster para no agotar la memoria virtual.

mde <- readCELL6sp("mde")
extract.mde <- overlay(mde, points)
rm(mde)

slope <- readCELL6sp("slope")
extract.slope <- overlay(slope, points)
rm(slope)

aspect <- readCELL6sp("aspect")
extract.aspect <- overlay(aspect, points)
rm(aspect)

rad.summer <- readCELL6sp("beam_radiation_172")
extract.rad.sum <- overlay(rad.summer, points)
rm(rad.summer)

rad.winter <- readCELL6sp("beam_radiation_354")

extract.rad.win <- overlay(rad.winter, points)
rm(rad.winter)

wetness <- readCELL6sp("wetness_index")
extract.wetness <- overlay(wetness, points)
rm(wetness)

Y ahora juntamos toda la información extraida en un único data.frame.

points@data <- cbind(points@data, mde = extract.mde@data, slope = extract.slope@data, aspect = extract.aspect@data, rad.summer = extract.rad.sum@data, rad.winter = extract.rad.win@data, wetness = extract.wetness)

El resultado final, es por tanto, un data.frame con información sobre valores de decaimiento (incluidos en el archivo points) y toda una serie de variables físicas obtenidas a partir de capas raster (pendiente, altitud, etc) que podemos analizar estadísticamente en R, cosa que en principio no podríamos hacer en GRASS. En la siguiente entrada se muestra un análisis preliminar de estos datos en R.