R statistics/rgrass: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(New, draft version)
 
(rewrite to rgrass7 and North Carolina sample dataset)
Line 16: Line 16:


Get the GRASS environment (mapset, region, map projection, etc.); you can display the metadata for your location by printing G:
Get the GRASS environment (mapset, region, map projection, etc.); you can display the metadata for your location by printing G:
  > gmeta()
  > G <- gmeta()
<pre>
<pre>
gisdbase    /home/neteler/grassdata  
gisdbase    /home/neteler/grassdata  
Line 37: Line 37:


List available vector maps:
List available vector maps:
  execGRASS("g.list", parameters = list(type = "vector"))
  > execGRASS("g.list", parameters = list(type = "vector"))


List selected vector maps (wildcard):
List selected vector maps (wildcard):
  execGRASS("g.list", parameters = list(type = "vector", pattern = "precip*"))
  > execGRASS("g.list", parameters = list(type = "vector", pattern = "precip*"))


Save selected vector maps into R vector:
Save selected vector maps into R vector:
  my_vmaps <- execGRASS("g.list", parameters = list(type = "vector", pattern = "precip*"))
  > my_vmaps <- execGRASS("g.list", parameters = list(type = "vector", 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 = "raster"))
  > execGRASS("g.list", parameters = list(type = "raster"))


List selected raster maps (wildcard):
List selected raster maps (wildcard):
  execGRASS("g.list", parameters = list(type = "raster", pattern = "lsat7_2000*"))
  > execGRASS("g.list", parameters = list(type = "raster", pattern = "lsat7_2002*"))


=== Reading in data ===
=== Reading in data ===
Read in two raster maps (North Carolina sample dataset):
Read in two raster maps (North Carolina sample dataset):
  > ncdata <- readRAST(c("geology", "elevation"),
  > ncdata <- readRAST(c("geology_30m", "elevation"), cat=c(TRUE, FALSE), ignore.stderr=TRUE, plugin=NULL)
            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, in this case a SpatialGridDataFrame object filled with data from two North Carolina layers. Here is a plot of the elevation data:
  > image(spear, attr = 2, col = terrain.colors(20))
  > image(ncdata, attr = 2, col = terrain.colors(20))


Add a title to the plot:
Add a title to the plot:
  > title("Spearfish elevation")
  > title("North Carolina elevation")


[[Image:R_stats_elev.png|center]]
[[Image:ncdata.png|center]]


In addition, we can show what is going on inside the objects read into R:
In addition, we can show what is going on inside the objects read into R:
  > '''str(G)'''
  > str(G)
<pre>
<pre>
List of 26
List of 26
$ GISDBASE    : chr "/home/rsb/topics/grassdata"
$ LOCATION_NAME: chr "spearfish57"
$ MAPSET      : chr "rsb"
  $ DEBUG        : chr "0"
  $ DEBUG        : chr "0"
  $ GRASS_GUI    : chr "text"
  $ LOCATION_NAME: chr "nc_spm_08_grass7"
  $ projection  : chr "1 (UTM)"
  $ GISDBASE    : chr "/home/veroandreo/grassdata"
  $ zone        : chr "13"
$ MAPSET      : chr "PERMANENT"
  $ datum        : chr "nad27"
  $ GUI          : chr "wxpython"
  $ ellipsoid    : chr "clark66"
  $ projection  : chr "99"
  $ north        : num 4928010
  $ zone        : chr "0"
  $ south        : num 4913700
  $ n            : num 228500
  $ west        : num 589980
  $ s            : num 215000
  $ east        : num 609000
  $ w            : num 630000
  $ top          : num 1
  $ e            : num 645000
  $ bottom      : num 0
  $ t            : num 1
  $ nsres        : num 30
  $ b            : num 0
  $ nsres3      : num 30
  $ nsres        : num 27.5
  $ ewres        : num 30
  $ nsres3      : num 10
  $ ewres3      : num 30
  $ ewres        : num 37.5
  $ ewres3      : num 10
  $ tbres        : num 1
  $ tbres        : num 1
  $ rows        : int 477
  $ rows        : int 491
  $ rows3        : int 477
  $ rows3        : int 1350
  $ cols        : int 634
  $ cols        : int 400
  $ cols3        : int 634
  $ cols3        : int 1500
  $ depths      : int 1
  $ 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"
$ cells        : chr "196400"
$ cells3      : chr "2025000"
  $ proj4        : chr "+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 +"| __truncated__
- attr(*, "class")= chr "gmeta"
</pre>
</pre>


 
  > summary(ncdata)
  > '''summary(spear)'''
<pre>
<pre>
Object of class SpatialGridDataFrame
Object of class SpatialGridDataFrame
Coordinates:
Coordinates:
              min     max
        min   max
coords.x1  589980  609000
[1,] 630000 645000
coords.x2 4913700 4928010
[2,] 215000 228500
Is projected: TRUE  
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]
proj4string :
Number of points: 2
[+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]
Grid attributes:
Grid attributes:
   cellcentre.offset cellsize cells.dim
   cellcentre.offset cellsize cells.dim
1           589995       30      634
1         630018.8 37.50000       400
2           4913715      30       477
2         215013.7 27.49491       491
Data attributes:
Data attributes:
      geology     elevation.dem 
  geology_30m     elevation    
  sandstone:74959   Min.  : 1066  
  CZfg_217:70381   Min.  : 55.92  
  limestone:61355   1st Qu.: 1200  
  CZig_270:66861   1st Qu.: 94.78  
  shale    :46423   Median : 1316  
  CZbg_405:24561   Median :108.88  
  sand    :36561   Mean  : 1354  
  CZlg_262:19287   Mean  :110.38  
  igneous  :36534   3rd Qu.: 1488  
  CZam_862: 6017   3rd Qu.:126.78  
  (Other)  :37636   Max.  : 1840  
  CZbg_910: 4351   Max.  :156.25  
  NA's    : 8950  NA's  :10101 
  (Other) : 4942                 
</pre>
</pre>


Line 128: Line 128:


We can create a table of cell counts:
We can create a table of cell counts:
  > '''table(spear$geology)'''
  > table(ncdata$geology_30m)
{| class="wikitable" border="1"
{| class="wikitable" border="1"
!metamorphic
!CZfg_217
!transition
!CZlg_262
!igneous
!CZig_270
!sandstone
!CZbg_405
!limestone
!CZve_583
!shale
!CZam_720 
!sandy shale
!CZg_766
!claysand
!CZam_862
!sand
!CZbg_910 
!Km_921
!CZbg_945
!CZam_946
!CZam_948
|-
|-
|11693
|70381
|142
|19287
|36534
|66861   
|74959
|24561   
|61355
|2089     
|46423
|467     
|11266
|691   
|14535
|6017
|36561
|4351
|1211       
|1     
|398     
|85
|}
|}


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)'''
  > execGRASS("r.stats", flags=c("c", "l"), parameters=list(input="geology_30m"), ignore.stderr=TRUE)
<pre>
<pre>
1 metamorphic 11693
217 CZfg 70381
2 transition 142
262 CZlg 19287
3 igneous 36534
270 CZig 66861
4 sandstone 74959
405 CZbg 24561
5 limestone 61355
583 CZve 2089
6 shale 46423
720 CZam 467
7 sandy shale 11266
766 CZg 691
8 claysand 14535
862 CZam 6017
9 sand 36561
910 CZbg 4351
* no data 8950
921 Km 1211
945 CZbg 1
946 CZam 398
948 CZam 85
</pre>
</pre>


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)'''
  > boxplot(ncdata$elevation ~ ncdata$geology_30m, medlwd = 1)


[[Image:R_stats_boxplot_geo.png|center]]
[[Image:boxplot_geo.png|center]]


=== Exporting data back to GRASS ===
=== Exporting data back to GRASS ===
Line 177: Line 187:


First prepare some data:  (square root of elevation)
First prepare some data:  (square root of elevation)
  > spear$sqdem <- sqrt(spear$elevation.dem)
  > ncdata$sqdem <- sqrt(ncdata$elevation)


Export data from ''R'' back into a GRASS raster map:
Export data from ''R'' back into a GRASS raster map:
  > writeRAST6(spear, "sqdemSP", zcol="sqdem", ignore.stderr=TRUE)
  > writeRAST(ncdata, "sqdemNC", zcol="sqdem", ignore.stderr=TRUE)


Check that it imported into GRASS ok:
Check that it imported into GRASS ok:
  > '''execGRASS("r.info", parameters=list(map="sqdemSP"))'''
  > execGRASS("r.info", parameters=list(map="sqdemNC"))


<pre>
<pre>
  +----------------------------------------------------------------------------+
  +----------------------------------------------------------------------------+
  | Layer:   sqdemSP                       Date: Sun May 14 21:59:26 2006   |
  | Map:     sqdemNC                       Date: Sun Jul 19 13:06:34 2015   |
  | Mapset:  rsb                            Login of Creator: rsb            |
  | Mapset:  PERMANENT                      Login of Creator: veroandreo      |
  | Location: spearfish57                                                      |
  | Location: nc_spm_08_grass7                                                |
  | DataBase: /home/rsb/topics/grassdata                                      |
  | DataBase: /home/veroandreo/grassdata                                      |
  | Title:    ( sqdemSP )                                                    |
  | Title:    ( sqdemNC )                                                    |
| Timestamp: none                                                            |
  |----------------------------------------------------------------------------|
  |----------------------------------------------------------------------------|
  |                                                                            |
  |                                                                            |
  |  Type of Map:  raster             Number of Categories: 255              |
  |  Type of Map:  raster               Number of Categories: 0              |
  |  Data Type:    FCELL                                                     |
  |  Data Type:    DCELL                                                     |
  |  Rows:        477                                                       |
  |  Rows:        491                                                       |
  |  Columns:      634                                                       |
  |  Columns:      400                                                       |
  |  Total Cells:  302418                                                     |
  |  Total Cells:  196400                                                     |
  |        Projection: UTM (zone 13)                                          |
  |        Projection: Lambert Conformal Conic                                |
  |            N:   4928010   S:   4913700   Res:   30                    |
  |            N:     228500   S: 215000.0002   Res: 27.49490794              |
  |            E:    609000   W:    589980   Res:   30                     |
  |            E:    645000   W:    630000   Res: 37.5                     |
  |  Range of data:    min =  32.649654 max = 42.895222                      |
  |  Range of data:    min = 7.47818253045085 max = 12.5000787351036        |
|                                                                            |
|  Data Source:                                                            |
|                                                                            |
|                                                                            |
  |                                                                            |
  |                                                                            |
  |  Data Description:                                                        |
  |  Data Description:                                                        |
  |    generated by r.in.gdal                                                  |
  |    generated by r.in.bin                                                  |
  |                                                                            |
  |                                                                            |
|  Comments:                                                                |
|    r.in.bin -d input="/home/veroandreo/grassdata/nc_spm_08_grass7/PERMA\  |
|    NENT/.tmp/localhost.localdomain/X666" output="sqdemNC" bytes=8 heade\  |
|    r=0 bands=1 order="native" north=228500 south=215000.0002 east=64500\  |
|    0 west=630000 rows=491 cols=400 anull=6                                |
  |                                                                            |
  |                                                                            |
  +----------------------------------------------------------------------------+
  +----------------------------------------------------------------------------+

Revision as of 16:15, 19 July 2015

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:

> 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_2002*"))

Reading in data

Read in two raster maps (North Carolina sample dataset):

> ncdata <- readRAST(c("geology_30m", "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, in this case a SpatialGridDataFrame object filled with data from two North Carolina layers. Here is a plot of the elevation data:

> image(ncdata, attr = 2, col = terrain.colors(20))

Add a title to the plot:

> title("North Carolina elevation")

In addition, we can show what is going on inside the objects read into R:

> str(G)
List of 26
 $ DEBUG        : chr "0"
 $ LOCATION_NAME: chr "nc_spm_08_grass7"
 $ GISDBASE     : chr "/home/veroandreo/grassdata"
 $ MAPSET       : chr "PERMANENT"
 $ GUI          : chr "wxpython"
 $ projection   : chr "99"
 $ zone         : chr "0"
 $ n            : num 228500
 $ s            : num 215000
 $ w            : num 630000
 $ e            : num 645000
 $ t            : num 1
 $ b            : num 0
 $ nsres        : num 27.5
 $ nsres3       : num 10
 $ ewres        : num 37.5
 $ ewres3       : num 10
 $ tbres        : num 1
 $ rows         : int 491
 $ rows3        : int 1350
 $ cols         : int 400
 $ cols3        : int 1500
 $ depths       : int 1
 $ cells        : chr "196400"
 $ cells3       : chr "2025000"
 $ proj4        : chr "+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 +"| __truncated__
 - attr(*, "class")= chr "gmeta"
> summary(ncdata)
Object of class SpatialGridDataFrame
Coordinates:
        min    max
[1,] 630000 645000
[2,] 215000 228500
Is projected: TRUE 
proj4string :
[+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]
Grid attributes:
  cellcentre.offset cellsize cells.dim
1          630018.8 37.50000       400
2          215013.7 27.49491       491
Data attributes:
   geology_30m      elevation     
 CZfg_217:70381   Min.   : 55.92  
 CZig_270:66861   1st Qu.: 94.78  
 CZbg_405:24561   Median :108.88  
 CZlg_262:19287   Mean   :110.38  
 CZam_862: 6017   3rd Qu.:126.78  
 CZbg_910: 4351   Max.   :156.25  
 (Other) : 4942                   

Summarizing data

We can create a table of cell counts:

> table(ncdata$geology_30m)
CZfg_217 CZlg_262 CZig_270 CZbg_405 CZve_583 CZam_720 CZg_766 CZam_862 CZbg_910 Km_921 CZbg_945 CZam_946 CZam_948
70381 19287 66861 24561 2089 467 691 6017 4351 1211 1 398 85

And compare with the equivalent GRASS module:

> execGRASS("r.stats", flags=c("c", "l"), parameters=list(input="geology_30m"), ignore.stderr=TRUE)
217 CZfg 70381
262 CZlg 19287
270 CZig 66861
405 CZbg 24561
583 CZve 2089
720 CZam 467
766 CZg 691
862 CZam 6017
910 CZbg 4351
921 Km 1211
945 CZbg 1
946 CZam 398
948 CZam 85

Create a box plot of geologic types at different elevations:

> boxplot(ncdata$elevation ~ ncdata$geology_30m, 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)

> ncdata$sqdem <- sqrt(ncdata$elevation)

Export data from R back into a GRASS raster map:

> writeRAST(ncdata, "sqdemNC", zcol="sqdem", ignore.stderr=TRUE)

Check that it imported into GRASS ok:

> execGRASS("r.info", parameters=list(map="sqdemNC"))
 +----------------------------------------------------------------------------+
 | Map:      sqdemNC                        Date: Sun Jul 19 13:06:34 2015    |
 | Mapset:   PERMANENT                      Login of Creator: veroandreo      |
 | Location: nc_spm_08_grass7                                                 |
 | DataBase: /home/veroandreo/grassdata                                       |
 | Title:     ( sqdemNC )                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    DCELL                                                      |
 |   Rows:         491                                                        |
 |   Columns:      400                                                        |
 |   Total Cells:  196400                                                     |
 |        Projection: Lambert Conformal Conic                                 |
 |            N:     228500    S: 215000.0002   Res: 27.49490794              |
 |            E:     645000    W:     630000   Res:  37.5                     |
 |   Range of data:    min = 7.47818253045085  max = 12.5000787351036         |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.bin                                                   |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.bin -d input="/home/veroandreo/grassdata/nc_spm_08_grass7/PERMA\   |
 |    NENT/.tmp/localhost.localdomain/X666" output="sqdemNC" bytes=8 heade\   |
 |    r=0 bands=1 order="native" north=228500 south=215000.0002 east=64500\   |
 |    0 west=630000 rows=491 cols=400 anull=6                                 |
 |                                                                            |
 +----------------------------------------------------------------------------+