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)$resOut
List 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 :10101
Summarizing 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