GRASS/R

Z GeoWikiCZ

Tento článek není dokončen nebo potřebuje upravit


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
  • 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
  • RArcInfo - import dat ve formátu ArcInfo Coverage
  • RColorBrewer - tabulky barev optimalizované pro tématické mapy
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("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

Informaci o syntaxi GRASS modulů poskytuje funkce parseGRASS. 

> parseGRASS('g.list')

Command: g.list 
Description: Lists available GRASS data base files of the user-specified data type to standard output. 
Keywords: general, map management 
Parameters:
  name: type, type: string, required: yes, multiple: yes
  keydesc: datatype, keydesc_count: 1
[Data type]
  name: mapset, type: string, required: no, multiple: no
[Mapset to list (default: current search path)]
Flags:
  name: f [Verbose listing (also list map titles)]
  name: verbose [Verbose module output]
  name: quiet [Quiet module output]

Nápověda

> help.start()

Nápověda dané funkce 

> ?gmeta6

Nápověda je dostupná ve formě HTML stránek 

> options(htmlhelp = TRUE)

či textově orientovaných manuálových stránek 

> options(htmlhelp = FALSE)

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

Informace o atributech vytiskne také funkce table(). 

> table(elev_lu$landuse96_28m)

                      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 
                                1610

Přístup k vektorovým datům

Pro čtení vektorových vrstev se používá funkce readVECT6. 

> geology <- readVECT6('geology')
> summary(geology)

Object of class SpatialPolygonsDataFrame
Coordinates:
         min      max
r1 123971.19 930172.3
r2  10875.83 318117.4
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 +ellps=GRS80 +datum=NAD83
+units=m +no_defs +towgs84=0,0,0]
Data attributes:
      cat           onemap_pro          PERIMETER            GEOL250_     
 Min.   :   1.0   Min.   :1.133e+03   Min.   :    166.9   Min.   :   2.0  
 1st Qu.: 458.8   1st Qu.:6.330e+05   1st Qu.:   3700.8   1st Qu.: 459.8  
 Median : 916.5   Median :3.184e+06   Median :  10308.4   Median : 917.5  
 Mean   : 916.5   Mean   :6.966e+07   Mean   :  43054.7   Mean   : 917.5  
 3rd Qu.:1374.2   3rd Qu.:1.745e+07   3rd Qu.:  29277.1   3rd Qu.:1375.2  
 Max.   :1832.0   Max.   :1.401e+10   Max.   :2729482.2   Max.   :1833.0  
                                                                          
   GEOL250_ID        GEO_NAME      SHAPE_area          SHAPE_len        
 Min.   :   1.0   Qp     : 326   Min.   :1.133e+03   Min.   :    166.9  
 1st Qu.: 458.8   Tt     : 116   1st Qu.:6.330e+05   1st Qu.:   3700.8  
 Median : 916.5   PzZu   :  86   Median :3.184e+06   Median :  10308.4  
 Mean   : 916.5   CZms   :  59   Mean   :6.966e+07   Mean   :  43054.7  
 3rd Qu.:1374.2   CZbg   :  57   3rd Qu.:1.745e+07   3rd Qu.:  29277.1  
 Max.   :1832.0   CZfv   :  54   Max.   :1.401e+10   Max.   :2729482.1  
                  (Other):1134

Základní informace o typu vektorových prvků poskytuje vInfo. 

> vInfo('geology')
    nodes     points      lines boundaries  centroids      areas    islands 
     4556          0          0       3649       1832       1832        907 
    faces    kernels primitives      map3d 
        0          0       5481          0

Alternativně 

> system('v.info -t map=geology')

nodes=4556
points=0
lines=0
boundaries=3649
centroids=1832
areas=1832
islands=907
faces=0
kernels=0
primitives=5481
map3d=0

Externí odkazy

Citováno z „https://geo.fsv.cvut.cz/index.php?title=GRASS/R&oldid=15256“