GRASS GIS / Poznámky pro družicová data Landsat
(přesměrováno z GRASS GIS - Poznámky pro družicová data Landsat)
Tato stránka obsahuje poznámky pro stažení dat LandSat Thematic Mapper a následné zpracování těchto dat v prostřední GRASS GIS. Více v rámci předmětu Zpracování obrazových dat.
Postup pro manuální stažení dat
- Data jsou dostupná na adrese: http://eros.usgs.gov (bližší informace 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:
 |
 |
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ář).
#!/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
Archiv obsahuje několik souborů ve formátu 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 gdalinfo poslouží pro kontrolu metadat, např.
gdalinfo L5191025_02520090801_B10.TIF
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
Vytvoření lokace
Hromadný import dat
Data lze naimportovat pomocí modulu r.in.gdal (File -> Import raster data -> Common import formats), ukázka z příkazové řadky (přepínač -e 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 r.info.
r.info map=b10
+----------------------------------------------------------------------------+ | 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" | | | +----------------------------------------------------------------------------+
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).
#!/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
Příklad spuštění:
- ./import.sh projde všechny adresáře
- ./import.sh LM41890261983200FFF03 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.
#!/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()
wxGUI
Dále lze data naimportovat pomocí wxGUI.
File -> Import raster data -> Common formats import
Vytvoření přehledky
Následující skript v jazyce Python vytvoří vektorovou přehledku včetně atributové tabulky s názvy mapsetů.
#!/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')
# 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)
for m in maps:
grass.run_command('g.remove',
quiet = True,
vect = m)
Test lokace
Následující skript v jazyce Python vypíše pro danou lokaci (město) relevantní mapsety včetně roku, ve kterém byly snímky Landsat pořízeny.
#!/usr/bin/python
#%Module
#% description: Vypis mapsety pro danou lokaci
#%End
#%Option
#% key: mesto
#% description: Nazev mesta (bez diakritiky)
#% type: string
#% required: yes
#%End
import os
import sys
import atexit
import grass.script as grass
tmp_maps = []
def cleanup():
grass.run_command('g.remove',
vect = '%s' % ','.join(tmp_maps))
def main():
os.environ['GRASS_VERBOSE'] = '0'
os.environ['GRASS_OVERWRITE'] = '1'
grass.run_command('v.extract',
input = 'obce',
output = 'obce_tmp',
where = "nazev_eng = '%s'" % options['mesto'])
tmp_maps.append('obce_tmp')
mapsets = []
# 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')
tmp_maps.append(mapset)
grass.run_command('v.select',
flags = 't',
ainput = mapset,
atype = 'area',
binput = 'obce_tmp',
btype = 'area',
output = 'mapset_tmp')
if bool(grass.vector_info_topo('mapset_tmp')['areas']):
ret = grass.read_command('r.timestamp',
map = 'B1@%s' % mapset)
year = ''
if ret:
year = ret.split(' ')[-1].strip()
print mapset, year
grass.run_command('g.remove',
vect = 'mapset_tmp')
return 0
if __name__ == "__main__":
atexit.register(cleanup)
options, flags = grass.parser()
sys.exit(main())
Příklad:
test.py mesto='Ceske Budejovice'
LE71900262010185ASN00 2010 LM41910261984105FFF03 1984 LT41910261992191XXX02 1992 LT51900261991173XXX02 1991 LT51900262002171MTI00 2002 LT51900262010193MOR00 2010 LT51910262002178MTI00 2002 LT51910262010184MOR00 2010 LT51920262002137MTI00 2002
Předzpracování dat (no-data)
Po importu dat je vhodné nahradit hodnotu '0' za 'null' (žádná data, no-data) - r.null.
g.copy rast=B1@ceske-budejovice-2000,B1 g.region rast=B1 r.null map=B1 setnull=0
Tuto operaci lze provést hromadně pro všechny kanály družicové scény z vybraného mapsetu (více k tématu skriptování v jazyce Python na portálu FreeGIS).
#!/usr/bin/python
from grass.pygrass.gis import Mapset
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
mapset = 'ceske-budejovice-2000'
for band in sorted(Mapset(mapset).glist('rast')):
g.message("Kopiruji %s..." % band)
g.copy(rast = '%s@%s,%s' % (band,mapset,band))
g.region(rast = band)
g.message("Nastavuji no-data pro %s..." % band)
r.null(map = band, setnull='0')
Další možností je vytvořit model pomocí grafického modeleru 
File → Graphical Modeler
Příklad modelu - pro každý kanál
- vytvoří kopii kanálu v aktuálním mapsetu (nutné pro r.null)
- nastaví region
- nahradí hodnotu 0 hodnotou NULL
Nastavení vyhledávací cesty
Mapset s naimportovanými daty přidáme do vyhledávací cesty pomocí g.mapsets
g.mapsets mapset=LM41890261983200FFF03 operation=add
nebo pomocí wxGUI
Settings -> GRASS working environment -> Mapset access
Vizualizace dat
Zájmové území
Nejprve vytvoříme vektorovou vrstvu s polygonem zájmového města (v tomto případě 'Pardubice' - název města zadejte bez diakritiky).
v.extract input=obce output=par where="nazev_eng = 'Pardubice'
Výpočetní region nastavíme pomocí g.region, např. pro město "Pardubice" (rastrovou vrstvu zadáváme pro nastavení prostorového rozlišení) s offsetem 10km.
g.region rast=B1 vect=par n=n+1e4 s=s-1e4 w=w-1e4 e=e+1e4 save=par
Opětovné nastavení výpočetního regionu
g.region region=par
Interpolace chybějících dat
Pro interpolaci chybějících hodnot lze použít modul r.fillnulls.
g.region rast=B1 r.fillnulls input=B1 output=B1_fill
Atmosferické korekce
Další poznámky najdete na GRASS Wiki.
Provedeme konverzi digitálních hodnot (DH) na hodnoty záření na vrcholu atmosféry (TOA) pomocí modulu i.landsat.toar (soubory s metadaty jsou dostupné na adrese: http://geo102.fsv.cvut.cz/zodh/metadata/).
i.landsat.toar -r input_prefix=Brg output_prefix=B_rad \
metfile=/work/geodata/zodh2010/data/jablonec-nad-nisou-1980/L4191025_02519840414_MTL.txt
Připravíme konfigurační soubor pro 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
cd /work/geodata/zodh2010/data/jablonec-nad-nisou-1980/
grep -a 'ACQUISITION_DATE' L4191025_02519840414_MTL.txt
ACQUISITION_DATE = 1984-04-14
grep -a 'SCENE_CENTER_SCAN_TIME' L4191025_02519840414_MTL.txt
SCENE_CENTER_SCAN_TIME = 09:17:32.4240000Z
- Střed scény v zeměpisných souřadnicích
g.region -lg
... center_long=15.16068438 center_lat=50.72648033 ...
nebo
g.region -c
center easting: 511341.280354 center northing: 5619440.568602
echo '511341|5619440' | m.proj -ogd --q
15.16068045|50.72664675
Poznámka: V prvním případě se používá elipsoid lokace, v druhém WGS84.
- Průměrná výška
r.univar map=dem --q | grep 'mean:'
mean: 552.711
Příklad konfiguračního souboru:
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
Atmosferické korekce nakonec aplikujeme příkazem
i.atcorr iimg=B_rad1 icnd=atcorr.txt oimg=B_cor1
Detekce oblačnosti
Detekovat mraky (a jejich teplotu) umožňuje modul i.landsat.acca.
Nejprve konvertuje DH na hodnoty odrazivosti (soubory s metadaty (MTL) jsou dostupné zde).
i.landsat.toar input_prefix=Brg output_prefix=B_ref metfile=L4191025_02519840414_MTL.txt
Provedeme výpočet (pouze pro TM/ETM+)
i.landsat.acca -5 input_prefix=B_ref output=B_acca
r.report map=B_acca units=p
+-----------------------------------------------------------------------------+ | 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| +-----------------------------------------------------------------------------+
Nastavení masky
Příklad nastavení masky oblačnosti pomocí r.mask
r.mask -i input=B_acca maskcats="6 9"
Výpočet pozice slunce
Pozici slunce můžeme určit pomocí r.sunmask. Nejprve z metadat určíme vstupní parametry.
cd /work/geodata/zodh2010/data/jablonec-nad-nisou-1980/
grep -a 'ACQUISITION_DATE' L4191025_02519840414_MTL.txt
ACQUISITION_DATE = 1984-04-14
grep -a 'SCENE_CENTER_SCAN_TIME' L4191025_02519840414_MTL.txt
SCENE_CENTER_SCAN_TIME = 09:17:32.4240000Z
r.sunmask -sg elevation=dem year=1984 month=4 day=14 hour=9 minute=17 second=32 timezone=0
sunazimuth=143.930862 sunangleabovehorizon=44.188877
Vypočtena data můžeme porovnat s hodnotami uvedenými v metadatovém souboru.
grep -a 'SUN_' L4191025_02519840414_MTL.txt
SUN_AZIMUTH = 143.8553375
SUN_ELEVATION = 44.1660847
Klasifikace dat
Neřízená klasifikace
Více na stránkách 10. cvičení předmětu Zpracování obrazových dat.
Řízená klasifikace
Více na stránkách 11. cvičení předmětu Zpracování obrazových dat.
Pokud máme k dispozici výsledky analýzy oblačnosti, muže je použít pro vytvoření masky
r.mask -i input=B_acca maskcats="6 9"
Nejprve vytvoříme obrazovou skupinu (předpokládejme, že Brg jsou obrazové kanály po atmosferické korekci)
i.group group=Brg subgroup=Brg input=Brg1,Brg2,Brg3,Brg4,Brg5,Brg7
nebo pomocí dialogu wxGUI Imagery → Develop images and groups → Create/edit group
 |
 |
Barevnou kompozici lze vytvořit pomocí r.composite
g.region rast=Brg4
r.composite red=Brg4 green=Brg5 blue=Brg3 output=Brg453
Provedeme digitalizaci vybraných trénovacích ploch, viz cvičení předmětu Zpracování obrazových dat. Trénovací plochy definujeme pro následující třídy:
- Zástavba
- Nelesní zeleň
- Lesní porosty
- Zemědělská půda
- Vodní plochy
Digitalizovat trénovací plochy lze ve specializovaném nástroji g.gui.iclass (více zde). Anebo nepřímo pomocí digitalizačního nástroje - viz návod níže.
 |
 |
Po přidání atributového sloupce GRASSRGB (v.db.addcolumn) můžeme pro vektorovou vrstvu trénovacích ploch nastavit tabulku barev.
v.db.addcol map=tp column="GRASSRGB varchar(11)"
Vector -> Manage colors -> Color rules
Pro aktivaci tabulky barev použijeme přepínač -a modulu d.vect
d.vect -a map=tp
Dále provedeme rasterizaci trénovacích ploch
v.to.rast input=tp type=area output=tp use=attr column=cat rgbcolumn=GRASSRGB labelcolumn=trida
A na závěr klasifikaci dat
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
- pomocí klasifikátoru SMAP (i.gensigset, i.smap)
i.gensigset trainingmap=tp group=Brg subgroup=Brg signaturefile=sig_smap
i.smap group=Brg subgroup=Brg signaturefile=sig_smap output=Brg_smap
Můžeme převzít tabulku barev pomocí příkazu
r.colors map=Brg_mlc raster=tp
r.colors map=Brg_smap raster=tp
Využití dat z projektu OpenStreetMap
Pro určení trénovacích ploch lze také využít tématické vrstvy z mapsetu osm vytvořeného z dat OpenStreetMap. Příklad vytvoření pomocných tématických vrstev (viz wiki OSM):
v.extract input=landuse output=landuse_zastavba where="type = 'brownfield' or type = 'cemetery' or type = 'commercial' \
or type = 'construction' or type = 'farmyard' or type = 'garages' or type = 'industrial' or type = 'residential' \
or type = 'retail'" new=1
v.extract input=landuse output=landuse_zelen where="type = 'allotments' or type = 'grass' or type = 'greenfield' or \
type = 'greenhouse_horticulture' or type = 'recreation_ground' or type='village_green' or type = 'park'" new=2
v.extract input=natural output=landuse_les where="type = 'forest'" new=3
v.extract input=landuse output=landuse_zemedel where="type = 'farm' or type = 'farmland' or type = 'orchard' \
or type = 'vineyard'" new=4
v.extract input=landuse output=landuse_voda where="type = 'reservoir'" new=5
Alternativně můžeme vytvořit "reklasifikovanou" landuse rastrovou mapu:
v.proj location=osm_shp mapset=PERMANENT input=landuse_zastavba
v.proj location=osm_shp mapset=PERMANENT input=landuse_zelen
v.proj location=osm_shp mapset=PERMANENT input=landuse_les
v.proj location=osm_shp mapset=PERMANENT input=landuse_zemedel
v.proj location=osm_shp mapset=PERMANENT input=landuse_voda
g.region rast=B1@beroun-2000 vect=cr
v.to.rast input=landuse_zastavba output=landuse_zastavba type=area use=val value=1
v.to.rast input=landuse_zelen output=landuse_zelen type=area use=val value=2
v.to.rast input=landuse_les output=landuse_les type=area use=val value=3
v.to.rast input=landuse_zemedel output=landuse_zemedel type=area use=val value=4
v.to.rast input=landuse_voda output=landuse_voda type=area use=val value=5
r.patch input=landuse_zastavba,landuse_zelen,landuse_les,landuse_zemedel,landuse_voda out=landuse_tridy
r.colors map=landuse_tridy rules=- << EOF
1 139:105:20
2 59:227:59
3 7:121:7
4 191:191:191
5 0:0:255
nv 255:255:255
default 255:255:255
EOF
Využití dat z veřejně dostupných WMS
Jako podkladovou vrstvu můžeme také použít data z WMS serveru
http://www.bnhelp.cz/ows/crtopo
Seznam dostupných vrstev získáme příkazem r.in.wms v kombinaci s přepínačem -l
r.in.wms -l mapserver=http://www.bnhelp.cz/ows/crtopo
Postklasifikační úpravy
Viz 10.cvičení.
r.reclass.area input=Brg_mlc output=tmp1 greater=1
r.grow.distance input=tmp1 value=tmp2
r.neighbors input=tmp2 output=Brg_mlc_post method=mode
r.colors map=Brg_mlc_post raster=tp
g.remove rast=tmp1,tmp2
Skript na závěr
Tento skript pro vybraný mapset (scénu) provede výše uvedené analýzy.
#!/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())