Difference between revisions of "R statistics/spgrass6"
m |
Veroandreo (talk | contribs) (add GRASS within R section) |
||
Line 1: | Line 1: | ||
− | '''This page refers to the usage of R within a GRASS GIS 6 session.''' (see also [[R_statistics/rgrass7]]) | + | '''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 === | === Startup === | ||
Line 21: | Line 23: | ||
List available vector maps: | List available vector maps: | ||
− | execGRASS("g.list", parameters = list(type = "vect")) | + | > execGRASS("g.list", parameters = list(type = "vect")) |
List selected vector maps (wildcard): | List selected vector maps (wildcard): | ||
− | execGRASS("g.list", parameters = list(type = "vect", pattern = "precip*")) | + | > execGRASS("g.list", parameters = list(type = "vect", pattern = "precip*")) |
Save selected vector maps into R vector: | Save selected vector maps into R vector: | ||
− | my_vmaps <- execGRASS("g.list", parameters = list(type = "vect", pattern = "precip*")) | + | > my_vmaps <- execGRASS("g.list", parameters = list(type = "vect", pattern = "precip*")) |
− | attributes(my_vmaps) | + | > attributes(my_vmaps) |
− | attributes(my_vmaps)$resOut | + | > attributes(my_vmaps)$resOut |
List available raster maps: | List available raster maps: | ||
− | execGRASS("g.list", parameters = list(type = "rast")) | + | > execGRASS("g.list", parameters = list(type = "rast")) |
List selected raster maps (wildcard): | List selected raster maps (wildcard): | ||
− | execGRASS("g.list", parameters = list(type = "rast", pattern = "lsat7_2000*")) | + | > execGRASS("g.list", parameters = list(type = "rast", pattern = "lsat7_2000*")) |
=== Reading in data === | === Reading in data === | ||
Read in two raster maps (Spearfish sample dataset): | Read in two raster maps (Spearfish sample dataset): | ||
− | > spear <- readRAST6(c("geology", "elevation.dem"), | + | > 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: | 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: | ||
Line 83: | Line 83: | ||
</pre> | </pre> | ||
− | + | > summary(spear) | |
− | > | ||
<pre> | <pre> | ||
Object of class SpatialGridDataFrame | Object of class SpatialGridDataFrame | ||
Line 112: | Line 111: | ||
We can create a table of cell counts: | We can create a table of cell counts: | ||
− | > | + | > table(spear$geology) |
{| class="wikitable" border="1" | {| class="wikitable" border="1" | ||
!metamorphic | !metamorphic | ||
Line 136: | Line 135: | ||
And compare with the equivalent GRASS module: | And compare with the equivalent GRASS module: | ||
− | > | + | > execGRASS("r.stats", flags=c("c", "l"), parameters=list(input="geology"), ignore.stderr=TRUE) |
<pre> | <pre> | ||
1 metamorphic 11693 | 1 metamorphic 11693 | ||
Line 152: | Line 151: | ||
Create a box plot of geologic types at different elevations: | Create a box plot of geologic types at different elevations: | ||
− | > | + | > boxplot(spear$elevation.dem ~ spear$geology, medlwd = 1) |
[[Image:R_stats_boxplot_geo.png|center]] | [[Image:R_stats_boxplot_geo.png|center]] | ||
Line 167: | Line 166: | ||
Check that it imported into GRASS ok: | Check that it imported into GRASS ok: | ||
− | > | + | > execGRASS("r.info", parameters=list(map="sqdemSP")) |
<pre> | <pre> | ||
Line 209: | Line 208: | ||
install.packages("automap") | install.packages("automap") | ||
install.packages("rgeos") | 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 | ||
+ | |||
[[Category:FAQ]] | [[Category:FAQ]] |
Revision as of 09:16, 19 July 2015
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)
Contents
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