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

Z GeoWikiCZ
Přejít na: navigace, hledání
m (Výpočet NDVI)
m (Výpočet NDVI)
Řádek 67: Řádek 67:
 
  > library(spgrass6)
 
  > library(spgrass6)
 
  > library(RColorBrewer)
 
  > library(RColorBrewer)
 +
 +
<source lang="bash">
 +
> GDALinfo('lsat7_2002.tif')
 +
</source>
 +
 +
</pre>
 +
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
 +
</pre>
  
 
Vyexportovaná data načteme.
 
Vyexportovaná data načteme.

Verze z 1. 12. 2009, 00:49

< Stránky předmětuPředchozí cvičení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í R pro GRASS GIS jsou dostupné 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_1
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(spgrass6)
> library(RColorBrewer)
> GDALinfo('lsat7_2002.tif')

</pre> 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 </pre>

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.

> 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))
Obr č.1.: NDVI (GRASS/R)



< Stránky předmětuPředchozí cvičeníDalší cvičení