R/GRASS connection: more than the sum of its parts

GRASS GIS is a very powerfull GIS offering an extensive set of tools for geospatial data management and analysis, image processing, and spatial modeling. The possibility to directly interact with R further improves the geospatial analysis capabilities of GRASS (see R Spatial View).

One of the advantages less often mentioned is that GRASS GIS can handle very large spatial data layers (especially raster data). It can thus be used to circumvent one of the drawbacks of R, its inability to deal with (very) large data sets (both analyses and visualization).

Suppose you want to create a species distribution model based on x number of sample points and 10 environmental variables. For building the model you can simple sample the raster maps at the locations of the vector points in GRASS using the v.what.rast function in GRASS. Next, you use the readVECT6() function in the R spgrass6 package to import the vector point layer in R (you need to run R from within GRASS).

Predictions from the fitted species distribution model can be easily mapped in the GRASS-R interface using the predict function and a data frame of raster variables that were used in the model. But this is where you might run into trouble if dealing with very large datasets. For example, 10 environmental layers of 500 MB each means you will need a lot of RAM to run a model in R.

One solution is to break up your data in smaller bits. Luckily this is very easy in GRASS using the r.tileset function. To access this function, look in the Raster menu –> Develop Raster map –> Tiling. Or type in r.tileset on the command line. The function produces an user-defined number of tiles (in the same or other projection then your current mapset).

It doesn’t actually produces the raster tiles themselves, but a text (file or console) with the coordinates of the bounding boxes of each tile and the number of columns and rows. By running the r.tileset through R (using system calls, you can run any GRASS command in R, as long as you start R on the GRASS GIS command line), you can capture those tile parameters and use this in a loop to import the environmental data layers and create the species suitability distribution maps tile by tile.

For example, if you have created a model gam1 based on 3 climatic variables, you can create a species distribution layers using something like:

require(spgrass6)
tilecoord <- system("r.tileset -g sourceproj=+init=epsg:4326 maxrows=1000 maxcols=1000 sf= ", intern=TRUE)
tilecoord <- gsub(";", " ", tilecoord)
for(i in 1:length(a)){
    system(paste("g.region", tilecoord[i]))
    r1 <- readRAST6("clim1@climate")
    r2 <- readRAST6("clim2@climate")
    r3 <- readRAST6("clim3@climate")
    predict = as.numeric(predict(gam1,newdata=list(clim1=r1@data, clim2=r2@data, clim3=r3@data),type="response"))
    r1@data <- predict
    writeRAST6(r1, "specdistr")
}