153ZODH / 13. cvičení: Porovnání verzí

Aktuální verze z 5. 9. 2014, 18:49

Tento text se vztahuje k předmětu, který již není součástí aktuálních studijních plánů.
Stránky navazujícího předmětu 155ZDDP jsou dostupné zde.

← předchozí cvičení • stránky předmětu • další cvičení →

Úvod do GRASS/R

Osnova

Toto cvičení je zaměřeno propojení GRASS GIS a prostředí pro statistickou analýzu dat R.

Seznam příkazů

Rozhraní GRASS/R

Základní informace o rozhraní zde.

Výpočet NDVI

Spustíme GRASS s lokací nc_spm_08. K dispozici máme snímky Landsat 7 z roku 2002.

g.mlist type=rast pattern='lsat7*'

lsat7_2002_10
lsat7_2002_20
lsat7_2002_30
lsat7_2002_40
lsat7_2002_50
lsat7_2002_61
lsat7_2002_62
lsat7_2002_70
lsat7_2002_80

Vytvoříme obrazovou skupinu, která obsahuje 1-3, 4-5 a 7 pásmo.

i.group group=lsat7_2002 input=lsat7_2002_10,lsat7_2002_20,lsat7_2002_30,lsat7_2002_40,lsat7_2002_50,lsat7_2002_70

Nastavíme region a data vyexportujeme do formátu GeoTIFF.

g.region rast=lsat7_2002_10
r.out.gdal input=lsat7_2002 format=GTiff type=Byte output=lsat7_2002.tif

Soubor by měl obsahovat 6 kanálů.

gdalinfo lsat7_2002.tif | grep Band

Band 1 Block=527x2 Type=Byte, ColorInterp=Gray
Band 2 Block=527x2 Type=Byte, ColorInterp=Undefined
Band 3 Block=527x2 Type=Byte, ColorInterp=Undefined
Band 4 Block=527x2 Type=Byte, ColorInterp=Undefined
Band 5 Block=527x2 Type=Byte, ColorInterp=Undefined
Band 6 Block=527x2 Type=Byte, ColorInterp=Undefined

Z příkazové řádky GRASSu spustíme prostředí R a nahrajeme balíček 'spgrass6' a 'RColorBrewer'.

GRASS (nc_spm_08):~ > R
> library(rgdal)
> library(spgrass6)
> library(RColorBrewer)

> GDALinfo('lsat7_2002.tif')

rows        475 
columns     527 
bands       6 
origin.x        629992.5 
origin.y        214975.5 
res.x       28.5 
res.y       28.5 
oblique.x   0 
oblique.y   0 
driver      GTiff 
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 +ellps=GRS80 +datum=NAD83
+units=m +no_defs 
file        lsat7_2002.tif 
apparent band summary:
  GDType Bmin Bmax
1   Byte    0  255
2   Byte    0  255
3   Byte    0  255
4   Byte    0  255
5   Byte    0  255
6   Byte    0  255

Vyexportovaná data načteme.

> lsat <- readGDAL('lsat7_2002.tif')
> summary(lsat)

Object of class SpatialGridDataFrame
Coordinates:
       min    max
x 629992.5 645012
y 214975.5 228513
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]
Number of points: 2
Grid attributes:
  cellcentre.offset cellsize cells.dim
x          630006.8     28.5       527
y          214989.8     28.5       475
Data attributes:
     band1            band2            band3           band4       
 Min.   : 42.00   Min.   : 28.00   Min.   :  1.0   Min.   :  1.00  
 1st Qu.: 72.00   1st Qu.: 55.00   1st Qu.: 41.0   1st Qu.: 81.00  
 Median : 80.00   Median : 64.00   Median : 57.0   Median : 93.00  
 Mean   : 85.47   Mean   : 70.92   Mean   : 66.7   Mean   : 93.16  
 3rd Qu.: 92.00   3rd Qu.: 79.00   3rd Qu.: 81.0   3rd Qu.:106.00  
 Max.   :255.00   Max.   :255.00   Max.   :255.0   Max.   :255.00  
     band5           band6       
 Min.   :  1.0   Min.   :  1.00  
 1st Qu.: 71.0   1st Qu.: 34.00  
 Median : 88.0   Median : 52.00  
 Mean   : 97.6   Mean   : 61.03  
 3rd Qu.:116.0   3rd Qu.: 79.00  
 Max.   :255.0   Max.   :255.00

Nahradíme názvy jednotlivých kanálů.

> names(lsat)

[1] "band1" "band2" "band3" "band4" "band5" "band6"

> names(lsat) <- c('blue1', 'green2', 'red3', 'nir4', 'mir5', 'mir7')

Vypočteme NDVI, viz 4. cvičení.

> lsat$ndvi <- (lsat$nir4 - lsat$red3) / (lsat$nir4 + lsat$red3)
> summary(lsat$ndvi)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.9565  0.0000  0.2152  0.1900  0.4140  0.9787

Nastavíme vhodnou tabulku barev.

> mypal <- brewer.pal(5, "Greens")
> greens <- colorRampPalette(mypal)

Data zobrazíme.

> image(lsat, "ndvi", col = greens(20))