R statistics/spgrass6
This page refers to the usage of R within a GRASS GIS 6 session and the use of GRASS GIS 6 within an R session. (see also R_statistics/rgrass7)
R within GRASS
Startup
- First start a GRASS session. Then, at the GRASS command line start R.
- In this example we will use the sample Spearfish dataset.
Reset the region settings to the defaults
GRASS> g.region -d
Launch R from the GRASS prompt
GRASS> R
Load the spgrass6 library:
library(spgrass6)Get the GRASS environment (mapset, region, map projection, etc.); you can display the metadata for your location by printing G:
G <- gmeta6()Listing of existing maps
List available vector maps:
execGRASS("g.list", parameters = list(type = "vect"))List selected vector maps (wildcard):
execGRASS("g.list", parameters = list(type = "vect", pattern = "precip*"))Save selected vector maps into R vector:
my_vmaps <- execGRASS("g.list", parameters = list(type = "vect", pattern = "precip*"))
attributes(my_vmaps)
attributes(my_vmaps)$resOutList available raster maps:
execGRASS("g.list", parameters = list(type = "rast"))List selected raster maps (wildcard):
execGRASS("g.list", parameters = list(type = "rast", pattern = "lsat7_2000*"))Reading in data
Read in two raster maps (Spearfish sample dataset):
spear <- readRAST6(c("geology", "elevation.dem"), cat=c(TRUE, FALSE), ignore.stderr=TRUE, plugin=NULL)The metadata are accessed and available, but are not (yet) used to structure the sp class objects, here a SpatialGridDataFrame object filled with data from two Spearfish layers. Here is a plot of the elevation data:
image(spear, attr = 2, col = terrain.colors(20))Add a title to the plot:
title("Spearfish elevation")
In addition, we can show what is going on inside the objects read into R:
str(G)
List of 26
$ GISDBASE : chr "/home/rsb/topics/grassdata"
$ LOCATION_NAME: chr "spearfish57"
$ MAPSET : chr "rsb"
$ DEBUG : chr "0"
$ GRASS_GUI : chr "text"
$ projection : chr "1 (UTM)"
$ zone : chr "13"
$ datum : chr "nad27"
$ ellipsoid : chr "clark66"
$ north : num 4928010
$ south : num 4913700
$ west : num 589980
$ east : num 609000
$ top : num 1
$ bottom : num 0
$ nsres : num 30
$ nsres3 : num 30
$ ewres : num 30
$ ewres3 : num 30
$ tbres : num 1
$ rows : int 477
$ rows3 : int 477
$ cols : int 634
$ cols3 : int 634
$ depths : int 1
$ proj4 : chr "+proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +no_defs +nadgrids=/home/rsb/topics/grass61/grass-6.1.cvs/etc/nad/conus"
summary(spear)
Object of class SpatialGridDataFrame
Coordinates:
min max
coords.x1 589980 609000
coords.x2 4913700 4928010
Is projected: TRUE
proj4string : [+proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +no_defs +nadgrids=/home/rsb/topics/grass61/grass-6.1.cvs/etc/nad/conus]
Number of points: 2
Grid attributes:
cellcentre.offset cellsize cells.dim
1 589995 30 634
2 4913715 30 477
Data attributes:
geology elevation.dem
sandstone:74959 Min. : 1066
limestone:61355 1st Qu.: 1200
shale :46423 Median : 1316
sand :36561 Mean : 1354
igneous :36534 3rd Qu.: 1488
(Other) :37636 Max. : 1840
NA's : 8950 NA's :10101Summarizing data
We can create a table of cell counts:
table(spear$geology)| metamorphic | transition | igneous | sandstone | limestone | shale | sandy shale | claysand | sand |
|---|---|---|---|---|---|---|---|---|
| 11693 | 142 | 36534 | 74959 | 61355 | 46423 | 11266 | 14535 | 36561 |
And compare with the equivalent GRASS module:
execGRASS("r.stats", flags=c("c", "l"), parameters=list(input="geology"), ignore.stderr=TRUE)1 metamorphic 11693 2 transition 142 3 igneous 36534 4 sandstone 74959 5 limestone 61355 6 shale 46423 7 sandy shale 11266 8 claysand 14535 9 sand 36561 * no data 8950
Create a box plot of geologic types at different elevations:
boxplot(spear$elevation.dem ~ spear$geology, medlwd = 1)
Exporting data back to GRASS
Finally, a SpatialGridDataFrame object is written back to a GRASS raster map:
First prepare some data: (square root of elevation)
spear$sqdem <- sqrt(spear$elevation.dem)Export data from R back into a GRASS raster map:
writeRAST6(spear, "sqdemSP", zcol="sqdem", ignore.stderr=TRUE)Check that it imported into GRASS ok:
execGRASS("r.info", parameters=list(map="sqdemSP"))+----------------------------------------------------------------------------+ | Layer: sqdemSP Date: Sun May 14 21:59:26 2006 | | Mapset: rsb Login of Creator: rsb | | Location: spearfish57 | | DataBase: /home/rsb/topics/grassdata | | Title: ( sqdemSP ) | |----------------------------------------------------------------------------| | | | Type of Map: raster Number of Categories: 255 | | Data Type: FCELL | | Rows: 477 | | Columns: 634 | | Total Cells: 302418 | | Projection: UTM (zone 13) | | N: 4928010 S: 4913700 Res: 30 | | E: 609000 W: 589980 Res: 30 | | Range of data: min = 32.649654 max = 42.895222 | | | | Data Source: | | | | | | | | Data Description: | | generated by r.in.gdal | | | | | +----------------------------------------------------------------------------+
v.krige usage example
v.krige is a GRASS python script (Addon) which performs kriging operations in the GRASS environment, using R functions for the back-end interpolation. It is present in GRASS 6.5svn, and further developed in GRASS 7svn. It requires a number of dependencies: python-rpy2 (needs to be "Rpy2", "Rpy" will not do unless it is rpy 2.x), then the following R-CRAN packages:
* gstat, spgrass6 (as above)
install.packages(c("gstat","spgrass6"))
* maptools
install.packages("maptools")
* automap (optional), with gpclib (or rgeos)
install.packages("automap")
install.packages("rgeos")GRASS within R
To call GRASS GIS 6 functionality from R, use the initGRASS() function to define the GRASS settings:
> library(spgrass6)
# initialisation and the use of spearfish60 data
> initGRASS(gisBase = "/usr/local/grass-6.4.1", home = tempdir(),
gisDbase = "/home/neteler/grassdata/",
location = "spearfish60", mapset = "user1", SG="elevation.dem",
override = TRUE)
system("g.region -d")
# verify
gmeta6()
spear <- readRAST6(c("geology", "elevation.dem"),
cat=c(TRUE, FALSE), ignore.stderr=TRUE,
plugin=NULL)
summary(spear$geology)
Run this script with
R CMD BATCH batch.R
The result is (shorted here):
cat batch.Rout
R version 2.10.0 (2009-10-26)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
...
> library(spgrass6)
Loading required package: sp
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.2, released 2010/04/23
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: (GRASS not running)
>
> # initialisation and the use of spearfish60 data
> initGRASS(gisBase = "/usr/local/grass-6.4.1", home = tempdir(), gisDbase = "/home/neteler/grassdata/",
+ location = "spearfish60", mapset = "user1", SG="elevation.dem", override = TRUE)
gisdbase /home/neteler/grassdata/
location spearfish60
mapset user1
rows 477
columns 634
north 4928010
south 4913700
west 589980
east 609000
nsres 30
ewres 30
projection +proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +no_defs
+nadgrids=/usr/local/grass-6.4.1/etc/nad/conus +to_meter=1.0
Warning messages:
1: In dir.create(gisDbase) : '/home/neteler/grassdata' already exists
2: In dir.create(loc_path) :
'/home/neteler/grassdata//spearfish60' already exists
>
> system("g.region -d")
> # verify
> gmeta6()
gisdbase /home/neteler/grassdata/
location spearfish60
mapset user1
rows 477
columns 634
north 4928010
...
>
> spear <- readRAST6(c("geology", "elevation.dem"),
+ cat=c(TRUE, FALSE), ignore.stderr=TRUE,
+ plugin=NULL)
>
> summary(spear$geology)
metamorphic transition igneous sandstone limestone shale
11693 142 36534 74959 61355 46423
sandy shale claysand sand NA's
11266 14535 36561 8950
>
>
> proc.time()
user system elapsed
2.891 0.492 3.412