R statistics/rgrass
This page refers to the usage of R within a GRASS GIS 7 session. (see also R_statistics/spgrass6)
Startup
- First start a GRASS GIS session. Then, at the GRASS command line start R (for a 'rstudio' session, see below)
- In this example we will use the North Carolina sample dataset.
Reset the region settings to the defaults
GRASS 7.0.1svn (nc_spm_08_grass7):~ > g.region -d
Launch R from the GRASS prompt
GRASS 7.0.1svn (nc_spm_08_grass7):~ > R
Load the rgrass7 library:
> library(rgrass7)
Get the GRASS environment (mapset, region, map projection, etc.); you can display the metadata for your location by printing G:
> gmeta()
gisdbase /home/neteler/grassdata location nc_spm_08_grass7 mapset user1 rows 620 columns 1630 north 320000 south 10000 west 120000 east 935000 nsres 500 ewres 500 projection +proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +no_defs +a=6378137 +rf=298.257222101 +towgs84=0.000,0.000,0.000 +to_meter=1
Listing of existing maps
List available vector maps:
execGRASS("g.list", parameters = list(type = "vector"))
List selected vector maps (wildcard):
execGRASS("g.list", parameters = list(type = "vector", pattern = "precip*"))
Save selected vector maps into R vector:
my_vmaps <- execGRASS("g.list", parameters = list(type = "vector", pattern = "precip*"))
attributes(my_vmaps)
attributes(my_vmaps)$resOut
List available raster maps:
execGRASS("g.list", parameters = list(type = "raster"))
List selected raster maps (wildcard):
execGRASS("g.list", parameters = list(type = "raster", pattern = "lsat7_2000*"))
Reading in data
Read in two raster maps (North Carolina sample dataset):
> ncdata <- readRAST(c("geology", "elevation"),
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 | | | | | +----------------------------------------------------------------------------+