GRASS/R

Tato stránka je věnována propojení GRASS GIS a R.

R je programovací jazyk a prostředí určené pro statistickou analýzu dat a jejich grafické zobrazení. Jde o implementaci programovacího jazyka S pod svobodnou licencí. R již předstihlo počtem uživatelů komerční S a stalo se faktickým standardem v řadě oblastí statistiky.

Používá se z příkazové řádky, existuje však několik frontendů s grafickým rozhraním jako RKWard, R Commander nebo rozšíření do OpenOffice.org Calcu R4Calc.

R bývá také propojováno či využíváno v komerčních softwarech, např. v prostředí PASW mohou uživatelé přímo psát a spouštět programy v jazyce R nad otevřenými daty.

Další informace na wikipedii.

Balíčky pro analýzu geoprostorových dat

R nabízí až více než 2000 různých rozšíření (tzv. balíčku) specializovaných na různé analýzy dat. My se zaměříme na balíčky určené pro práci s geoprostorovými (většinou obrazovými) daty.

• sp - základní balíček definující třídy a metody pro práci s prostorovými daty

• maptools - manipulace s prostorovými objekty

• maps - vykreslování map

• spatial - Kriging

• spatstat - analýza bodových výskytů

• splancs - analýza bodových výskytů v čase

• spdep - autokorelace prostorových objektů

• gstat - geostatistické modelování

• geoR - analýza geostatistických dat

• fields - analýza geoprostorových dat

• spatialCovariance - výpočet kovariační matice prostorových dat

• RArcInfo - import dat ve formátu ArcInfo Coverage

• shapefiles - čtení a zápis ESRI Shapefile

• RColorBrewer - tabulky barev optimalizované pro tématické mapy

• spgrass6 - rozhraní pro GRASS GIS 6+

• rgdal - rozhraní pro knihovnu GDAL

Tab. č. 1: Datové typy definované balíčkem sp

Datový typ	Třída	Rodičovská třída
body	SpatialPoints	Spatial
pixely	SpatialPixels	SpatialPoints
mřížka	SpatialGrid	SpatialPixels
linie	SpatialLines	Spatial, Line
hranice	SpatialRings	SpatialLines
polygon	SpatialPolygons	Spatial, Polygon

spgrass6

Spgrass6 je rozšířením balíčku sp a je koncipován jako rozhraní pro moduly GRASS GIS 6. Přístup k datům v nativním formátu GRASS je zajištěn balíčkem rgdal.

Instalace

$ R
> install.packages(c("sp", "rgdal", "spgrass6"), dependencies = TRUE)

Spuštění

Z příkazové řádky GRASS GIS spustíme interpret R.

GRASS 6.5.svn (nc_spm_08):~ > R

Dále nahrajeme balíček spgrass6 a všechny jeho závislosti.

> 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.0dev, released 2008/11/26
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.6.1, 21 August 2008
Path to PROJ.4 shared files: (autodetected)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: 6.5.svn
and location: nc_spm_08

Metainformace:

> gmeta6()

gisdbase    /home/martin/grassdata 
location    nc_spm_08 
mapset      landa 
rows        1350 
columns     1500 
north       228500 
south       215000 
west        630000 
east        645000 
nsres       10 
ewres       10 
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 

Pro spuštění GRASS modulů se používá funkce system, např.

> system('g.region -p')

projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      227412.33194791
south:      223694.5675084
west:       635197.26282579
east:       639950.9754717
nsres:      9.99399043
ewres:      10.0078161
rows:       372
cols:       475
cells:      176700

Přístup k rastrovým datům

Nejprve nastavíme region.

> system('g.region rast=elevation')

Na základě rastrové vrstvy bude vytvořen objekt 'SpatialGridDataFrame'.

> elev <- readRAST6('elevation')

Základní informace o objektu poskytuje funkce summary().

> summary(elev)

Object of class SpatialGridDataFrame
Coordinates:
     min    max
x 630000 645000
y 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]
Number of points: 2
Grid attributes:
  cellcentre.offset cellsize cells.dim
x            630005       10      1500
y            215005       10      1350
Data attributes:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  55.58   94.79  108.90  110.40  126.80  156.30 

V případě klasifikovaných rastrových dat použijeme parametr 'CAT', např.

> lu <- readRAST6('landuse96_28m', cat = TRUE)
> summary(lu)

Object of class SpatialGridDataFrame
Coordinates:
     min    max
x 630000 645000
y 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]
Number of points: 2
Grid attributes:
  cellcentre.offset cellsize cells.dim
x            630005       10      1500
y            215005       10      1350
Data attributes:
                      not classified             High Intensity Developed 
                                   9                               283200 
             Low Intensity Developed                           Cultivated 
                              309088                                17347 
            Managed Herbaceous Cover        Riverine/Estuarine Herbaceous 
                              205883                                  246 
                 Evergreen Shrubland                  Deciduous Shrubland 
                              133042                                 2596 
                     Mixed Shrubland                      Mixed Hardwoods 
                                 356                                64058 
Bottomland Hardwoods/Hardwood Swamps                 Southern Yellow Pine 
                              159491                               524056 
            Mixed Hardwoods/Conifers                         Water Bodies 
                              274120                                42854 
             Unconsolidated Sediment                                 NA's 
                                1610                                 7044 

Poznámka: Statistika kategorií odpovídá příkazu r.stats -lc.

> system('r.stats -lc input=landuse96_28m --q')

0 not classified 9
1 High Intensity Developed 283200
2 Low Intensity Developed 309088
3 Cultivated 17347
4 Managed Herbaceous Cover 205883
6 Riverine/Estuarine Herbaceous 246
7 Evergreen Shrubland 133042
8 Deciduous Shrubland 2596
9 Mixed Shrubland 356
10 Mixed Hardwoods 64058
11 Bottomland Hardwoods/Hardwood Swamps 159491
15 Southern Yellow Pine 524056
18 Mixed Hardwoods/Conifers 274120
20 Water Bodies 42854
21 Unconsolidated Sediment 1610
* no data 7044

Rastrové vrstvy lze kombinovat, např.

> elev_lu <- readRAST6(c('elevation', 'landuse96_28m'), cat = c(FALSE, TRUE))
> summary(elev_lu)

Object of class SpatialGridDataFrame
Coordinates:
     min    max
x 630000 645000
y 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]
Number of points: 2
Grid attributes:
  cellcentre.offset cellsize cells.dim
x            630005       10      1500
y            215005       10      1350
Data attributes:
   elevation                       landuse96_28m   
 Min.   : 55.58   Southern Yellow Pine    :524056  
 1st Qu.: 94.79   Low Intensity Developed :309088  
 Median :108.88   High Intensity Developed:283200  
 Mean   :110.38   Mixed Hardwoods/Conifers:274120  
 3rd Qu.:126.79   Managed Herbaceous Cover:205883  
 Max.   :156.33   (Other)                 :421609  
                  NA's                    :  7044  

Externí odkazy

• R, GRASS, and Spatial Analysis
• GRASS Newsletter vol. 3
  • Roger Bivand: Interfacing GRASS 6 and R
• OSGeo Journal Volume 1
  • Roger Bivand: Using R - GRASS interface