Zobrazení zdroje stránky GRASS GIS / Poznámky pro družicová data Landsat

Tato stránka obsahuje poznámky pro stažení dat [http://en.wikipedia.org/wiki/Landsat_program LandSat] [http://en.wikipedia.org/wiki/Thematic_Mapper Thematic Mapper] a následné zpracování těchto dat v prostřední [[GRASS GIS]]. Více v rámci předmětu [[153ZODH|Zpracování obrazových dat]].
__TOC__
== Postup pro manuální stažení dat ==

* Data jsou dostupná na adrese: http://eros.usgs.gov (bližší informace [http://landsat.usgs.gov/products_data_at_no_charge.php zde])
** Data mohou být stažena také z adres: http://glovis.usgs.gov nebo http://edcsns17.cr.usgs.gov/EarthExplorer

Návod pro http://glovis.usgs.gov:

<center>
{|
| [[Image:glovis-0.png|thumb|400px]] || [[Image:glovis-1.png|thumb|400px]]
|}
[[Image:glovis-2.png|thumb|center|400px]]
</center>

== Import dat ==

Data je nejprve nutno dekomprimovat, např.

 tar xvzf LT51910252009213MOR00.tar.gz

Příklad skriptu pro [[Bash]], který rozbalí všechny dostupné archivy v aktuálním adresáři (pro každý rozbalený archiv vytvoří samostatný adresář).

<source lang="bash">
#!/bin/sh

for file in *.tar.gz ; do
    echo "Zpracovavam $file..."
    dir=${file%.tar.gz}
    if [ -d $dir ] ; then
        rm -rf $dir
    fi
    mkdir $dir
    tar -xzf $file -C $dir
done
</source>

Archiv obsahuje několik souborů ve formátu [http://en.wikipedia.org/wiki/Geotiff geotiff] (jeden soubor na kanál), seznam identických bodů (GCP, ground control points), metadata v textové podobě (MTL) a soubor README.GTF.

 ls -w1
 L5191025_02520090801_B10.TIF
 L5191025_02520090801_B20.TIF
 L5191025_02520090801_B30.TIF
 L5191025_02520090801_B40.TIF
 L5191025_02520090801_B50.TIF
 L5191025_02520090801_B60.TIF
 L5191025_02520090801_B70.TIF
 L5191025_02520090801_GCP.txt
 L5191025_02520090801_MTL.txt
 README.GTF

Nástroj knihovny [[GDAL/OGR|GDAL]] <tt>gdalinfo</tt> poslouží pro kontrolu metadat, např.

 gdalinfo L5191025_02520090801_B10.TIF

<pre>
Driver: GTiff/GeoTIFF
Files: L5191025_02520090801_B10.TIF
Size is 8441, 7441
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["unknown",1],
    AUTHORITY["EPSG","32633"]]
Origin = (398699.999999999941792,5679000.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  398700.000, 5679000.000) ( 13d32'54.25"E, 51d15'12.06"N)
Lower Left  (  398700.000, 5455770.000) ( 13d36'29.29"E, 49d14'46.73"N)
Upper Right (  651930.000, 5679000.000) ( 17d10'35.91"E, 51d14'31.60"N)
Lower Right (  651930.000, 5455770.000) ( 17d 5'13.62"E, 49d14'9.03"N)
Center      (  525315.000, 5567385.000) ( 15d21'18.49"E, 50d15'29.05"N)
Band 1 Block=8441x1 Type=Byte, ColorInterp=Gray
</pre>

=== Vytvoření lokace ===

<center>
<gallery widths=400px heights=300px perrow=2>
Image:loc-tm-0.png
Image:loc-tm-1.png
Image:loc-tm-2.png
Image:loc-tm-3.png
Image:loc-tm-4.png
Image:loc-tm-5.png
Image:loc-tm-6.png
Image:loc-tm-7.png
</gallery>
</center>

=== Hromadný import dat ===

Data lze naimportovat pomocí modulu {{GrassPrikaz|r.in.gdal}} (<tt>File -> Import raster data -> Common import formats</tt>), ukázka z příkazové řadky (přepínač <tt>-e</tt> nastaví po importu region podle vstupní dat).

 r.in.gdal -e input=L5191025_02520090801_B10.TIF output=b10 title="kanal 1"

Základní metadata vypíše {{GrassPrikaz|r.info}}.

 r.info map=b10

<pre>
 +----------------------------------------------------------------------------+
 | Layer:    b10                            Date: Thu Sep 16 17:48:23 2010    |
 | Mapset:   PERMANENT                      Login of Creator: martin          |
 | Location: tm                                                               |
 | DataBase: /home/martin/grassdata                                           |
 | Title:    kanal 1 ( b10 )                                                  |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    CELL                                                       |
 |   Rows:         7441                                                       |
 |   Columns:      8441                                                       |
 |   Total Cells:  62809481                                                   |
 |        Projection: UTM (zone 33)                                           |
 |            N:    5679000    S:    5455770   Res:    30                     |
 |            E:     651930    W:     398700   Res:    30                     |
 |   Range of data:    min = 0  max = 255                                     |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal -e input="L5191025_02520090801_B10.TIF" output="b10"          |
 |                                                                            |
 +----------------------------------------------------------------------------+
</pre>

Data lze hromadně naimportovat pomocí jednoduchého skriptu. Následuje ukázka pro [[Bash]] a [[Python]].

==== Bash ====

Ukázka skriptu pro [[Bash]] (pro každou družicovou scénu vytvoří vlastní mapset).

<source lang="bash">
#!/bin/sh
 
function import()
{
    idx=1
    for file in `ls *.TIF 2>/dev/null`; do
        oname=`echo $file | cut -d'_' -f3 | cut -d'.' -f1`
	oname=${oname:0:2}
        g.message "Importuji $file -> $oname@$mapset..."
        g.mapset -c mapset=$mapset --quiet 2>/dev/null
        r.in.gdal -e input=$file output=$oname title="kanal $idx" --overwrite --quiet
        idx=$(($idx+1))
    done
}
 
eval `g.gisenv`
 
if test -n "$1" ; then
    mapset=$1
    cd $1
    import
else
    for dir in `find . -maxdepth 1 -mindepth 1 -type d`; do
	mapset=`echo $dir | cut -d'/' -f2`
	cd $dir
	import
	cd ..
    done
fi
 
g.mapset mapset=$MAPSET --quiet
</source>

Příklad spuštění:

* <tt>./import.sh</tt> projde všechny adresáře
* <tt>./import.sh LM41890261983200FFF03</tt> naimportuje data pouze z vybraného adresáře

==== Python ====

Ukázka skriptu pro [[Python]] (pro každou družicovou scénu vytvoří vlastní mapset). Tento skript navíc nastaví timestamp rastrových dat.

<source lang="python">
#!/usr/bin/python
 
import os
import sys
import glob
import grass.script as grass
 
def get_timestamp(mapset):
    try:
        metafile = glob.glob(mapset + '/*MTL.txt')[0]
    except IndexError:
        return
 
    result = dict()
    try:
        fd = open(metafile)
        for line in fd.readlines():
            line = line.rstrip('\n')
            if len(line) == 0:
                continue
            if 'ACQUISITION_DATE' in line:
                result['date'] = line.strip().split('=')[1].strip()
    finally:
        fd.close()
 
    return result
 
def import_tifs(mapset):
    meta = get_timestamp(mapset)
    for file in os.listdir(mapset):
        if os.path.splitext(file)[-1] != '.TIF':
            continue
        ffile = os.path.join(mapset, file)
        name = os.path.splitext(file)[0].split('_')[-1] #B10, B61
        if name[-1] == '0':
            name = name[:2]
        kanal = int(name[1]) # B
        grass.message('Importuji %s -> %s@%s...' % (file, name, mapset))
        grass.run_command('g.mapset',
                          flags = 'c',
                          mapset = mapset,
                          quiet = True,
                          stderr = open(os.devnull, 'w'))
        grass.run_command('r.in.gdal',
                          input = ffile,
                          output = name,
                          quiet = True,
                          overwrite = True,
                          title = 'kanal %d' % kanal)
        if meta:
            year, month, day = meta['date'].split('-')
            if month == '01':
                month = 'jan'
            elif month == '02':
                month = 'feb'
            elif month == '03':
                month = 'mar'
            elif month == '04':
                month = 'apr'
            elif month == '05':
                month = 'may'
            elif month == '06':
                month = 'jun'
            elif month == '07':
                month = 'jul'
            elif month == '08':
                month = 'aug'
            elif month == '09':
                month = 'sep'
            elif month == '10':
                month = 'oct'
            elif month == '11':
                month = 'nov'
            elif month == '12':
                month = 'dec'
 
            grass.run_command('r.timestamp',
                              map = name,
                              date = ' '.join((day, month, year)))
 
def main():
    if len(sys.argv) == 1:
        for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
            import_tifs(directory)
    else:
        import_tifs(sys.argv[1])
 
if __name__ == "__main__":
    main()
</source>

==== wxGUI ====
Dále lze data naimportovat pomocí wxGUI.

 File -> Import raster data -> Bulk import of raster data

<center>
<gallery widths=400px heights=300x perrow=3>
Image:import-tm-0.png
Image:import-tm-1.png
</gallery>
</center>

=== Vytvoření přehledky ===

Následující skript v jazyce [[Python]] vytvoří vektorovou přehledku včetně atributové tabulky s názvy mapsetů.

<source lang=python>
#!/usr/bin/python

import os
import sys
import grass.script as grass

maps = [] # tmp maps
cats = {} # attribute table
cat = 1

os.environ['GRASS_VERBOSE'] = '0'
os.environ['GRASS_OVERWRITE'] = '1'

# list mapsets
for mapset in grass.mapsets():
    if mapset[:2] not in ('LE', 'LM', 'LT'):
        continue
    grass.run_command('g.region',
                      rast = "B1@%s" % mapset)
    grass.run_command('v.in.region',
                      output = mapset,
                      type = 'area',
                      cat = cat)
    cats[cat] = mapset
    maps.append(mapset)
    
    cat += 1

# patch maps
grass.run_command('v.patch',
                  input = ','.join(maps),
                  output = 'prehledka_tmp')
grass.run_command('v.type',
                  input = 'prehledka_tmp',
                  output = 'prehledka',
                  type = 'boundary,line')

# create attribute table
grass.run_command('v.db.addtable',
                  map = 'prehledka',
                  columns = 'cat integer, mapset varchar(255)')
for cat, mapset in cats.iteritems():
    grass.run_command('v.db.update',
                      map = 'prehledka',
                      column = 'mapset',
                      value = mapset,
                      where = "cat = %d" % cat)

maps.append('prehledka_tmp')
for m in maps:
    grass.run_command('g.remove',
                      quiet = True,
                      vect = m)
</source>

[[Image:landsat-prehledka.png|center|640px|thumb|Vektorová přehledka mapsetů s importovanými daty Landsat]]

== Předzpracování dat (no-data) ==

Po importu dat je vhodné nahradit hodnotu '0' za 'null' (žádná data, no-data) - {{GrassPrikaz|r.null}}.

 r.null map=B10 setnull=0

hromadně pro všechny mapsety (vyžaduje, vlastnictví adresáře mapsetu).

<source lang="python">
#!/usr/bin/python

import os
import grass.script as grass

cur_mapset = grass.gisenv()['MAPSET']

for mapset in grass.mapsets(accessible = False):
    if mapset != cur_mapset:
        grass.run_command('g.mapset',
                          mapset = mapset,
                          quiet = True)
    for rast in grass.read_command('g.mlist', type = 'rast', mapset = mapset).splitlines():
        grass.message("Zpracovavam %s@%s..." % (rast, mapset))
        grass.run_command('r.null', map = rast, setnull = 0)
    
grass.run_command('g.mapset',
                  mapset = cur_mapset,
                  quiet = True)
</source>

Další možností je vytvořit ''model'' pomocí grafického modeleru [[Image:wxgui-modeler.png]] - {{GrassPrikaz|wxGUI.Modeler}}

 File -> Graphical Modeler

Příklad modelu - pro každý kanál

# vytvoří kopii kanálu v aktuálním mapsetu (nutné pro {{GrassPrikaz|r.null}})
# nastaví region
# nahradí hodnotu 0 hodnotou NULL

{{YouTube|LFdn1VozMfo|desc=Model pro předzpracování družicových dat (hodnoty NULL)}}

== Nastavení vyhledávací cesty ==

Mapset s naimportovanými daty přidáme do vyhledávací cesty pomocí {{GrassPrikaz|g.mapsets}}

 g.mapsets addmapset=LM41890261983200FFF03

nebo pomocí {{GrassPrikaz|wxGUI}}

 Settings -> GRASS working environment -> Mapset access

[[Image:wxgui-mapset-access.png|center|thumb|300px|Nastavení přístupu k mapsetům]]

== Vizualizace dat ==

[[Image:wxgui-tm453.png|center|640px|thumb|Příklad vizualizace kompozitního snímku (kanály 453)]]

== Zájmové území ==

Nejprve vytvoříme vektorovou vrstvu s polygony zájmových měst ({{GrassPrikaz|v.extract}}).

 v.extract input=obce output=mesta where="nazev_eng = 'Ceske Budejovice' or \
  nazev_eng = 'Vsetin' or nazev_eng = 'Teplice' or nazev_eng = 'Jablonec nad Nisou' or \
  nazev_eng = 'Znojmo'"

Anebo pro každé město vytvoříme vlastní vektorovou vrstvu.

<source lang=python>
#!/usr/bin/python
 
import grass.script as grass

for mesto in ('Ceske Budejovice', 'Vsetin', 'Teplice', 'Jablonec nad Nisou', 'Znojmo'):
    grass.run_command('v.extract',
                      input = 'obce',
                      output = mesto.replace(' ', '_'),
                      where="nazev_eng = '%s'" % mesto)
</source>


[[Image:grass-region-vect.png|center|640px|thumb|Nastavení výpočetního regionu podle vektorové vrstvy]]

Výpočetní region nastavíme pomocí {{GrassPrikaz|g.region}}, např. pro město "Vsetín" (rastrovou vrstvu zadáváme pro nastavení prostorového rozlišení).

 g.region rast=B10 vect=Vsetin

alternativně lze nastavit výpočetní region širší, např. o 1 kilometr

 g.region rast=B10 vect=Vsetin n=n+1e3 s=s-1e3 w=w-1e3 e=e+1e3 

Dále vytvoříme masku zájmového území vybraného města.

 v.to.rast input=Vsetin output=MASK use=val value=1

[[Image:grass-region-mask.png|center|640px|thumb|Nastavení masky zájmové oblasti]]

== Interpolace chybějících dat ==

Pro interpolaci chybějících hodnot lze použít modul {{GrassPrikaz|r.fillnulls}}.

 g.region rast=B10
 r.fillnulls input=B10 output=B10_fill

== Atmosferické korekce ==

Další poznámky najdete na [http://grass.osgeo.org/wiki/Atmospheric_correction GRASS Wiki].

Nejprve nastavíme cestu k mapsetu a výpočetní region (s odstupem 10km).

<source lang="bash">
g.mapsets addmapset=jablonec-nad-nisou-1980
g.region vect=Jablonec_nad_Nisou rast=B1 n=n+1e4 s=s-1e4 w=w-1e4 e=e+1e4 save=region
</source>

''Poznámka:'' Následující moduly ignorují výpočetní region a operuje na celém gridu! Pokud si přejete pracovat pouze v aktuálním výpočetním regionu, musíte vytvořit kopie rastrových dat, např.

<source lang="python">
#!/usr/bin/python

import grass.script as grass

for band in grass.read_command('g.mlist',
                               type = 'rast',
                               pattern = 'B*').splitlines():
    oname = 'Brg' + band[1:]
    grass.message("Kopiruji %s -> %s..." % (band, oname))
    grass.mapcalc(oname + ' = ' + band)
    grass.run_command('r.null',
                      map = oname,
                      setnull = 0)
</source>

Provedeme konverzi digitálních hodnot (DH) na hodnoty záření na vrcholu atmosféry ([http://www.google.com/url?sa=t&source=web&cd=2&ved=0CCQQFjAB&url=http%3A%2F%2Fwww.delmar.edu%2Figett%2FResources%2FRemote%2520Sensing%2520Technology%2520Training%2FCalculation-DN_to_Reflectance_Irish_20June08.pdf&rct=j&q=remote%20sensing%20toa&ei=Jy_PTP6ML8OZOomAseAB&usg=AFQjCNEWFo9vnvWm20N7FQI4XPhmXkVw-g&cad=rja TOA]) pomocí modulu {{GrassPrikaz|i.landsat.toar}}.

<source lang="bash">
i.landsat.toar -tr input_prefix=Brg output_prefix=B_rad \
 metfile=/work/geodata/zodh2010/data/jablonec-nad-nisou-1980/L4191025_02519840414_MTL.txt
</source>

Připravíme konfigurační soubor pro {{GrassPrikaz|i.atcorr}}. Tento modul aplikuje atmosferické korekce, přepočte hodnoty záření na vrcholu atmosféry na hodnoty absolutní odrazivosti objektů.

Některé údaje získáme z metadatového souboru (MTL), jiné spočteme.

* Dat a čas pořízení dat

<source lang=python>
cd /work/geodata/zodh2010/data/jablonec-nad-nisou-1980/
grep -a 'ACQUISITION_DATE' L4191025_02519840414_MTL.txt
</source>
<pre>
    ACQUISITION_DATE = 1984-04-14
</pre>
<source lang=python>
grep -a 'SCENE_CENTER_SCAN_TIME' L4191025_02519840414_MTL.txt
</source>
<pre>
    SCENE_CENTER_SCAN_TIME = 09:17:32.4240000Z
</pre>

* Střed scény v zeměpisných souřadnicích

<source lang="bash">
 g.region -lg
</source>

<pre>
...
center_long=15.16068438
center_lat=50.72648033
...
</pre>

nebo 

<source lang="bash">
 g.region -c
</source>
<pre>
center easting:  511341.280354
center northing: 5619440.568602
</pre>
<source lang="bash">
echo '511341|5619440' | m.proj -ogd --q
</source>
<pre>
15.16068045|50.72664675
</pre>

Poznámka: V prvním případě se používá elipsoid lokace, v druhém WGS84.

* Průměrná výška

<source lang="bash">
r.univar map=dem --q | grep 'mean:'
</source>
<pre>
mean: 552.711
</pre>

Příklad konfiguračního souboru:

<pre>
7                                    # geometricke podminky Landsat TM
04 14 9.28 15.16068045 50.72664675   # mesic den h.mmm (GMT) zem. sirka delka
2                                    # atmosfericky model midlatitude summer
3                                    # aerosol model urban
50                                   # viditelnost [km]
-0.552                               # prumerna vyska [km]
-1000                                # sensor na urovni nosice
31                                   # 1. kanal Landsat 5 TM
</pre>

Atmosferické korekce nakonec aplikujeme příkazem

<source lang="bash">
i.atcorr iimg=B_rad1 icnd=atcorr.txt oimg=B_cor1
</source>

== Detekce oblačnosti ==

Detekovat mraky (a jejich teplotu) umožňuje modul [http://grass.osgeo.org/wiki/GRASS_AddOns#i.landsat.acca i.landsat.acca].

Nejprve konvertuje DH na hodnoty odrazivosti

<source lang="bash">
i.landsat.toar -t input_prefix=Brg output_prefix=B_ref metfile=L4191025_02519840414_MTL.txt
</source>

Provedeme výpočet (pouze pro TM/ETM+)

<source lang="bash">
i.landsat.acca -5 input_prefix=B_ref output=B_acca
</source>

<source lang="bash">
r.report map=B_acca units=p
</source>

<pre>                                              
+-----------------------------------------------------------------------------+
|                         RASTER MAP CATEGORY REPORT                          |
|LOCATION: zod2010                                    Mon Nov  8 23:05:41 2010|
|-----------------------------------------------------------------------------|
|          north: 5633223.50868776    east: 525403.19584718                   |
|REGION    south: 5605657.62851621    west: 497279.36486102                   |
|          res:        29.99551705    res:      30.01476092                   |
|-----------------------------------------------------------------------------|
|MASK:none                                                                    |
|-----------------------------------------------------------------------------|
|MAP: LANDSAT-5 TM Automatic Cloud Cover Assessment (B_acca@landa in landa)   |
|-----------------------------------------------------------------------------|
|                        Category Information                          |   %  |
|#|description                                                         | cover|
|-----------------------------------------------------------------------------|
|6|Cold cloud . . . . . . . . . . . . . . . . . . . . . . . . . . . . .|  0.01|
|9|Warm cloud . . . . . . . . . . . . . . . . . . . . . . . . . . . . .|  0.10|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 99.89|
|-----------------------------------------------------------------------------|
|TOTAL                                                                 |100.00|
+-----------------------------------------------------------------------------+
</pre>

=== Nastavení masky ===

Příklad nastavení masky oblačnosti pomocí {{GrassPrikaz|r.mask}}

<source lang=bash>
r.mask -i input=B_acca maskcats=6,9
</source>

[[Image:B321.png|thumb|center|640px|Barevná syntéza 321]]
[[Image:B_acca.png|thumb|center|640px|Výsledek detekce oblačnosti]]
[[Image:B321_cloud.png|thumb|center|640px|Barevná syntéza 321 s maskou oblačnosti]]

== Výpočet pozice slunce ==

Pozici slunce můžeme určit pomocí {{GrassPrikaz|r.sunmask}}. Nejprve z metadat určíme vstupní parametry.

<source lang=python>
cd /work/geodata/zodh2010/data/jablonec-nad-nisou-1980/
grep -a 'ACQUISITION_DATE' L4191025_02519840414_MTL.txt
</source>
<pre>
    ACQUISITION_DATE = 1984-04-14
</pre>
<source lang=python>
grep -a 'SCENE_CENTER_SCAN_TIME' L4191025_02519840414_MTL.txt
</source>
<pre>
    SCENE_CENTER_SCAN_TIME = 09:17:32.4240000Z
</pre>

<source lang=python>
r.sunmask -sg elevation=dem year=1984 month=4 day=14 hour=9 minute=17 second=32 timezone=0
</source>

<pre>
sunazimuth=143.930862
sunangleabovehorizon=44.188877
</pre>

Vypočtena data můžeme porovnat s hodnotami uvedenými v metadatovém souboru.

<source lang=python>
grep -a 'SUN_' L4191025_02519840414_MTL.txt
</source>

<pre>
    SUN_AZIMUTH = 143.8553375
    SUN_ELEVATION = 44.1660847
</pre>

== Klasifikace dat ==

Více na stránkách cvičení předmětu [[153YZOD Zpracování obrazových dat - cvičení 11|Zpracování obrazových dat]].

Pokud máme k dispozici výsledky [[GRASS GIS - Poznámky pro družicová data Landsat#Detekce oblačnosti|analýzy oblačnosti]], muže je použít pro vytvoření masky

<source lang="bash">
r.mask -i input=B_acca maskcats="6 9"
</source>

Vytvoříme obrazovou skupinu

<source lang="bash">
i.group group=Brg subgroup=Brg input=Brg1,Brg2,Brg3,Brg4,Brg5,Brg7
</source>

Barevnou kompozici lze vytvořit pomocí {{GrassPrikaz|r.composite}}

<source lang="bash">
g.region rast=Brg4
r.composite red=Brg4 green=Brg5 blue=Brg3 output=Brg453
</source>

Provedeme digitalizaci vybraných trénovacích ploch, viz cvičení předmětu [[153YZOD Zpracování obrazových dat - cvičení 11#Neinteraktivní|Zpracování obrazových dat]]. Trénovací plochy definujeme pro následující třídy:

* Zástavba
* Nelesní zeleň
* Lesy
* Zemědělské území
* Vodní plochy

[[Image:grass-attrb.png|thumb|640px|center|Atributová tabulka vektorové vrstvy trénovacích ploch]]
[[Image:grass-digit-settings.png|thumb|640px|center|Nastavení digitizéru při vektorizaci trénovacích ploch]]
[[Image:grass-digit-tp.png|thumb|640px|center|Atributy vektorové vrstvy trénovacích ploch]]

Jako podkladovou vrstvu můžeme použít data z WMS serveru

 http://www.bnhelp.cz/ows/crtopo

Seznam dostupných vrstev získáme příkazem {{GrassPrikaz|r.in.wms}} v kombinaci s přepínačem <tt>-l</tt>

<source lang="bash">
r.in.wms -l mapserver=http://www.bnhelp.cz/ows/crtopo
</source>

Po přidání atributového sloupce GRASSRGB ({{GrassPrikaz|v.addcolumn}}) můžeme pro vektorovou vrstvu trénovacích ploch nastavit tabulku barev.

<source lang="bash">
 v.db.addcol map=tp column="GRASSRGB varchar(11)"
</source>

[[Image:grass-vcolors.png|thumb|640px|center|Nastavení tabulky barev vektorové vrstvy trénovacích ploch]]

Pro aktivaci tabulky barev použijeme přepínač <tt>-a</tt> modulu {{GrassPrikaz|d.vect}}

<source lang="bash">
 d.vect -a map=tp
</source>

Dále provedeme rasterizaci trénovacích ploch

<source lang="bash">
 v.to.rast input=tp type=area output=tp use=attr column=cat rgbcolumn=GRASSRGB labelcolumn=trida
</source>

Na závěr provedeme klasifikaci 

* pomocí klasifikátoru MLC ({{GrassPrikaz|i.gensig}}, {{GrassPrikaz|i.maxlik}})

<source lang="bash">
i.gensig trainingmap=tp group=Brg subgroup=Brg signaturefile=sig
i.maxlik group=Brg subgroup=Brg sigfile=sig class=Brg_mlc reject=Brg_mlc_ref
</source>

* pomocí klasifikátoru SMAP ({{GrassPrikaz|i.gensigset}}, {{GrassPrikaz|i.smap}})

<source lang="bash">
i.gensigset trainingmap=tp group=Brg subgroup=Brg signaturefile=sig_smap
i.smap group=Brg subgroup=Brg signaturefile=sig_smap output=Brg_smap
</source>

Můžeme převzít tabulku barev pomocí příkazu

<source lang="bash">
r.colors map=Brg_mlc raster=tp
</source>

[[Image:grass-mlc.png|thumb|640px|center|Výsledek klasifikace (MLC)]]
[[Image:grass-smap.png|thumb|640px|center|Výsledek klasifikace (SMAP)]]

=== Postklasifikační úpravy ===

Viz {{153YZODCv|10|Postklasifikační úpravy|10.cvičení}}.

<source lang="bash">
r.reclass.area input=Brg_mlc output=Brg_mlc1 greater=1
r.surf.idw input=Brg_mlc1 output=Brg_mlc2
r.neighbors input=Brg_mlc2 output=Brg_mlc_post method=mode
r.colors map=Brg_mlc_post raster=tp
</source>

[[Image:grass-postkl.png|thumb|640px|center|Výsledek postklasifikačních úprav]]

== Skript na závěr ==

Tento skript pro vybraný mapset (scénu) provede výše uvedené analýzy.

<source lang=python>
#!/usr/bin/python

#%module
#% description: Davkovy vypocet pro ZODH.
#%end
#%option
#% key: mapset
#% description: Nazev mapsetu
#% gisprompt: old,mapset,mapset
#% key_desc: name
#% required: yes
##% answer: ceske-budejovice-2000
#%end
#%option
#% key: city
#% description: Nazev mesta
#% required: yes
#% options: Ceske Budejovice,Jablonec nad Nisou,Teplice,Vsetin,Znojmo
##% answer: Ceske Budejovice
#%end
#%option
#% key: mtl
#% description: Cesta k MTL souboru
#% gisprompt: old,file,file
#% key_desc: path
#% required: yes
##% answer: /work/geodata/zodh2010/data/ceske-budejovice-2000/L5191026_02620020627_MTL.txt
#%end
#%option
#% key: output_prefix
#% description: Prefix pro vystupni datove vrstvy
#% answer: B
#%end
#%flag
#% key: r
#% description: Pred vypoctem vycistit mapset
#%end

import os
import sys
import tempfile
 
import grass.script as grass
 
def create_atcorr_file(mtlfile):
    atcorrfile = list()
    text = ''
    bandid = None
    try:
        fd = open(mtlfile)
        for line in fd.readlines():
            if 'SENSOR_ID' in line:
                sensor = line.split('=')[1].strip().replace('"', '')
                if sensor != 'ETM+':
                    atcorrfile.append('7')
                    if sensor == 'MSS':
                        bandid = 31
                    else:
                        bandid = 25
                else:
                    atcorrfile.append('8')
                    bandid = 61
            elif 'ACQUISITION_DATE' in line:
                year, month, day = line.split('=')[1].strip().split('-')
                text = '%s %s ' % (month, day)
            elif 'SCENE_CENTER_SCAN_TIME' in line:
                time = line.split('=')[1].strip().replace('"', '')
                h, m, s = time.split(':')
                m = int(m) * 100 / 60
                text += '%d.%.2d ' % (int(h), m)
                region = grass.parse_command('g.region',
                                             flags = 'lg')
                text += '%f %f' % (float(region['center_long']),
                                   float(region['center_lat']))
                atcorrfile.append(text)
 
        atcorrfile.append('2')      # midlatitude summer
        atcorrfile.append('3')      # urban model
        atcorrfile.append('15')     # visibility
        stats = grass.parse_command('r.univar',
                                    flags = 'g',
                                    map = 'dem')
        atcorrfile.append('-%f' % (float(stats['mean']) / 1000))
        atcorrfile.append('-1000')
    finally:
        fd.close()
    
    return atcorrfile, bandid
 
def main():
    mapset  = options['mapset']
    mesto   = options['city']
    mtlfile = options['mtl']
    prefix  = options['output_prefix']
    
    if flags['r']:
        grass.info("Odstranuji vsechny rastrove a vektorove mapy...")
        grass.run_command('g.mremove',
                          flags = 'f',
                          rast = '*',
                          vect = '*')
    
    path = grass.mapsets()
    current_mapset = grass.gisenv()['MAPSET']
    
    # globalni promenne
    os.environ['GRASS_OVERWRITE'] = "1"
    os.environ['GRASS_VERBOSE']   = "1"
    
    # pridej mapset do cesty
    grass.run_command('g.mapsets',
                      mapset = '%s,PERMANENT,%s' % \
                          (current_mapset, mapset))
    
    # vytvor polygon mesta
    grass.info("Vytvarim polygon mesta...")
    grass.run_command('v.extract',
                      input = 'obce',
                      output = mesto.replace(' ', '_'),
                      where="nazev_eng = '%s'" % mesto)
    mesto = mesto.replace(' ', '_')
    
    # nastaveni regionu
    grass.info("Nastavuji vypocetni region...")
    grass.run_command('g.region',
                      vect = mesto,
                      rast = 'B1',
                      n = 'n+1e4',
                      s = 's-1e4',
                      w = 'w-1e4',
                      e = 'e+1e4')
    
    # vytvor kopii rastrovych dat v ramci vypocetniho regionu
    grass.info("Vytvarim kopie rastrovych dat...")
    for band in grass.read_command('g.mlist',
                                   type = 'rast',
                                   pattern = 'B*',
                                   mapset = mapset).splitlines():
        oname = '%srg%s' % (prefix, band[1:])
        grass.info(" %s -> %s" % (band, oname))
        grass.mapcalc(oname + ' = ' + band)
        grass.run_command('r.null',
                          map = oname,
                          setnull = 0)
        grass.run_command('r.colors',
                          map = oname,
                          color = 'grey.eq')
    
    # vypocet zareni
    grass.info("Konvertuji rastrova data DH -> hodnoty zareni...")
    grass.run_command('i.landsat.toar',
                      flags = 'tr',
                      input_prefix = '%srg' % prefix,
                      output_prefix = '%s_rad' % prefix,
                      metfile = mtlfile)
    
    # atmosfericke korekce
    grass.info("Aplikuji atmosfericke korekce...")
    atcorr, bandid = create_atcorr_file(mtlfile)
    isMss = bandid == 31
    i = 1
    for band in grass.read_command('g.mlist',
                                   type = 'rast',
                                   pattern = '%s_rad*' % prefix).splitlines():
        oband = band.split('_', 1)[0] + '_cor' + str(i)
        grass.info("Zpracovavam <%s> -> <%s>..." % (band, oband))
        atcorr.append('%d' % bandid)
        atcorrfile = tempfile.NamedTemporaryFile()
        atcorrfile.write('\n'.join(atcorr))
        atcorrfile.seek(0)
        grass.run_command('i.atcorr',
                          iimg = band,
                          icnd = atcorrfile.name,
                          oimg = oband)
        bandid += 1 
        i += 1
        atcorr.pop()
    
    # hodnoty odrazivost
    grass.info("Konvertuji rastrova data DH -> hodnoty odrazivosti...")
    grass.run_command('i.landsat.toar',
                      flags = 't',
                      input_prefix = '%srg' % prefix,
                      output_prefix = '%s_ref' % prefix,
                      metfile = mtlfile)
    
    # detekce mraku
    grass.info("Analyzuji oblacnost...")
    if isMss:
        grass.info("...nelze pro MSS")
    else:
        grass.run_command('i.landsat.acca',
                          flags = '5',
                          input_prefix = '%s_ref' % prefix,
                          output = '%s_acca' % prefix)
    
    # obnov vyhledavaci cestu
    grass.run_command('g.mapsets',
                      mapset = ','.join(path))
    
    return 0
 
if __name__ == "__main__":
    options, flags = grass.parser()
    sys.exit(main())
</source>

== Externí odkazy ==

* [http://geo.fsv.cvut.cz/~landa/vyuka/153ZODH/load%20USGS.pdf Poznámky z cvičení YZD1]
* [http://grass.osgeo.org/wiki/LANDSAT GRASS User Wiki - some notes about LandSat]
* [http://grass.osgeo.org/wiki/MODIS#Removing_holes Removing holes]
* [http://toroid.org/ams/landsat-and-grass An introduction to using Landsat data with GRASS]
{{GRASS}}