Difference between revisions of "R statistics/spgrass6"

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)

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)

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