153ZODH / 13. 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))
