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.
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 -> Import raster data), 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 -> Bulk import of raster data
Předzpracování dat
Po importu dat je vhodné nahradit hodnotu '0' za 'null' (žádná data, no-data) - r.null.
r.null map=B10 setnull=0
hromadně pro všechny mapsety (vyžaduje, vlastnictví adresáře mapsetu).
#!/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)
Nastavení vyhledávací cesty
Mapset s naimportovanými daty přidáme do vyhledávací cesty pomocí g.mapsets
g.mapsets addmapset=LM41890261983200FFF03
nebo pomocí wxGUI
Settings | GRASS working environment | Mapset access
Vizualizace dat
Zájmové území
Nejprve vytvoříme vektorovou vrstvu s polygony zájmových měst (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' \ nazev_eng = 'Znojmo'"
Anebo pro každé město vytvoříme vlastní vektorovou vrstvu.
#!/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)
Výpočetní region nastavíme pomocí 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
Interpolace chybějících dat
g.region rast=B10 r.fillnulls input=B10 output=B10_fill
Atmosferické korekce
Další poznámky najdete na GRASS Wiki.
Nejprve nastavíme cestu k mapsetu a výpočetní region (s odstupem 10km)
g.mapsets addmapset=LM41910251984105FFF03
g.region vect=Jablonec_nad_Nisou rast=B1 n=n+1e4 s=s-1e4 w=w-1e4 e=e+1e4
Provedeme konverzi digitálních hodnot (DH) na hodnoty záření na vrcholu atmosféry (TOA) pomocí modulu i.landsat.toar