153ZODH / 7. cvičení: Porovnání verzí
m →wxGUI: center |
m →wxGUI |
||
| Řádek 270: | Řádek 270: | ||
|} | |} | ||
</center> | </center> | ||
V mapovém okně aktuální (tj. cílové) a zdrojové lokace se určí identické body (viz obr. č.x). | |||
[[Image:wxGUI-georect-3.png|center|thumb|450px|Volba identických bodů a rektifikace dat]] | |||
==== CLI ==== | ==== CLI ==== | ||
Verze z 11. 11. 2008, 21:23
< Stránky předmětu • Předchozí cvičení • Další cvičení
Osnova
Rastrové formáty
Satelitní (či obecně obrazová) data mohou být poskytována v rozličných datových formátech. Mezi nejběžnější patří GeoTIFF, ECW či BIL/BSQ. Obecně můžeme rozdělit datové formáty na dva základní typy:
- Datové formáty, které obsahují jedno pásmo (kanál) na jeden soubor (TIFF, PNG, SUN raster formát, a další)
- Datové formáty, které mohou obsahovat více kanálů na jeden soubor (BIL/BSQ, CEOS, ERDAS/LAN, HDF a další)
Obrazová data (či obecně geodata) mohou být georeferencovaná (tj. souřadnicově přípojena) či souřadnicově nepřipojená. Tento fakt velmi významně ovlivňuje proces jejich importu. Informace o souřadnicovém připojení by měla být poskytnuta dodavatelem dat (jako součást metadat - "dat o datech"), v opačném případě lze použít nástroj gdalinfo, který je součástí knihovny GDAL/OGR. Jako příklad uvedeme výpis metadat rastrového souboru '94T1.tif':
gdalinfo 94T1.tif
Driver: GTiff/GeoTIFF
Size is 2223, 1661
Coordinate System is:
PROJCS["unnamed",
GEOGCS["unnamed",
DATUM["unknown",
SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT[,0.0174532925199433],
AUTHORITY["EPSG","-32768"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","-32768"]]
Origin = (-830529.00000,-957500.00000)
Pixel Size = (29.9927754,-29.9927754)
Metadata:
TIFFTAG_DATETIME=2003:03:20 17:05:41
TIFFTAG_XRESOLUTION=0
TIFFTAG_YRESOLUTION=0
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Corner Coordinates:
Upper Left ( -830529.00, -957500.00)
Lower Left ( -830529.00, -1007318.00)
Upper Right ( -763855.06, -957500.00)
Lower Right ( -763855.06, -1007318.00)
Center ( -797192.03, -982409.00)
Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
V tomto případě se jedná o datový formát GeoTIFF (georeferencovaný TIFF). Mezi důležité informace patří geometrické rozlišení dat a souřadnice rohů obrazové scény.
Seznam podporovaných rastrových formátů knihovnou GDAL/OGR naleznete zde. Pro import/export dat GRASS totiž až na několik málo výjimek využívá právě tuto knihovnu.
Více informací naleznete v přednáškách k předmětu Free Software GIS.
Georeferencovaná data
Import
Moduly pro import rastrových dat začínají v GRASSu předponou r.in.*. Tak například ve verzi GRASS 6.4:
r.in.arc r.in.ascii r.in.aster r.in.bin r.in.gdal r.in.gridatb r.in.mat r.in.poly r.in.srtm r.in.wms r.in.xyz
Ve většině případů obstará import dat modul r.in.gdal. Knihovna GDAL podporuje skutečně velký počet rastrových a vektorových formátů.
Import dat si ukážeme na příkladu: Tarball etm-jtsk.tar.gz obsahuje družicová data LandSat7-ETM+. Tento soubor nejprve dekomprimujeme, rozbalíme.
wget http://gama.fsv.cvut.cz/~landa/vyuka/153YZOD/dataset/etm-jtsk.tar.gz tar xvzf etm-jtsk.tar.gz cd etm-jtsk
Vytiskneme si metadata prvního kanálu družicové scény.
# parametr '-noct' potlačí tisk tabulky barev # parametr '-nomd' potlačí tisk metadat gdalinfo -noct -nomd etm1.tif Driver: GTiff/GeoTIFF Files: etm1.tif Size is 2223, 1661 Coordinate System is `' Origin = (-830529.000000000000000,-957500.000000000000000) Pixel Size = (29.992775436482237,-29.992775436484045) Corner Coordinates: Upper Left ( -830529.000, -957500.000) Lower Left ( -830529.000,-1007318.000) Upper Right ( -763855.060, -957500.000) Lower Right ( -763855.060,-1007318.000) Center ( -797192.030, -982409.000) Band 1 Block=2223x1 Type=Int16, ColorInterp=Gray NoData Value=0
# import souboru ve formátu GeoTIFF # přepínač '-o' ignoruje aktuální nastavení souřadnicového systému lokace # r.in.gdal -o input=etm1.tif out=etm1 title="Prvni pasmo LandSat7 ETM+"
# nastavení regionu a tabulky barev # g.region rast=etm1 r.colors map=etm1 color=grey.eq d.rast etm1
Hromadný import lze provést velmi jednoduše pomocí primitivního skriptu (pro BASH), např.:
# hromadný import dat ve formátu GeoTIFF
#
for file in *.tif; do \
map=${file%%.tif}
r.in.gdal -o input=$file out=$map title="LandSat7 ETM+"; \
g.region rast=$map; \
r.colors map=$map color=grey.eq; \
done
wxGUI nabízí dialog pro hromadný import (File->Import raster map->Multiple raster data import using GDAL).
Připojení
Od verze 6.4 nabízí GRASS možnost rastrová data připojit přímo přes knihovnu GDAL bez nutnosti data importovat do nativního formátu GRASS. K tomu slouží modul r.external.
Georeferencování dat
Nejprve si stáhneme tarball mss-xy.tar.gz s pracovními daty.
wget http://gama.fsv.cvut.cz/~landa/vyuka/153YZOD/dataset/mss-xy.tar.gz tar xvzf mss-xy.tar.gz cd mss-xy
Jedná se skutečně o data souřadnicově nepřipojená:
gdalinfo mss1.tif Driver: GTiff/GeoTIFF Size is 1111, 623 Coordinate System is `' Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 623.0) Upper Right ( 1111.0, 0.0) Lower Right ( 1111.0, 623.0) Center ( 555.5, 311.5) Band 1 Block=1111x2 Type=Byte, ColorInterp=Red Band 2 Block=1111x2 Type=Byte, ColorInterp=Green Band 3 Block=1111x2 Type=Byte, ColorInterp=Blue
Georeferncování dat lze v GRASSu rozdělit do tří základních kroků:
- Založení lokace XY (tedy lokace s matematickým souřadnicovým systémem) a import vstupních dat do této lokace
- Určení vlícovacích bodů a jejich identifikace v souřadnicově nepřipojené vrstvě (lokace XY) a současně v podkladové vrstvě pro georefencování (rastrová či vektorová mapa) uložené v lokaci s cílovým souřadnicovým systémem.
- Transformace dat do cílové lokace
Založení lokace XY a import dat
Založit lokaci lze dvěma způsoby:
- Importem dat z existující lokace, parametr
locationmodulu r.in.gdal
r.in.gdal input=mss1.tif out=mss1 location=zod-xy
- Z uvítací obrazovky spustit "Location Wizard"
Pro import dat použijeme modul r.in.gdal
# import dat ve formátu TIFF # přepínač '-e' nastaví výpočetní region podle importované vrstvy # r.in.gdal -e in=mss1.tif out=mss1
Po importu se vytvoří tři obrazové vrstvy (mss1.red, mss1.green, mss1.blue). RGB syntézu lze vytvořit pomocí modulu r.composite.
g.region rast=mss1.red # # složení RGB kanálů # r.composite red=mss1.red green=mss1.green blue=mss1.blue out=mss1 # # odstranění nepotřebných vrstev # g.remove rast=mss1.red,mss1.green,mss1.blue
Pro hromadné zpracování:
for file in *.tif;do \
map=${file%%.tif}
r.in.gdal -e in=$file out=$map; \
g.region rast=${map}.red; \
r.composite red=${map}.red green=${map}.green blue=${map}.blue out=$map; \
g.remove rast=${map}.red,${map}.green,${map}.blue; \
done
Vzhledem k tomu, že se jedná vstupní data obsahují 8-bitová data, postačí importovat pouze jeden z kanálu.
r.in.gdal input=mss1.tif out=mss1 band=1
Volba vlícovacích bodů
Rozlišujeme čtyři případy rektifikace obrazu:
- obraz - obraz (matematický souřadnicový systém, souřadnicově nepřipojená data)
- obraz - rastr (pokladem pro georeferencování je rastrová vrstva)
- obraz - vektor (pokladem pro georeferencování je vektorová vrstva)
- obraz - souřadnice (souřadnice v cílovém souřadnicovém systému jsou zadávány číselně)
Pro transformaci obrazových dat se obvykle používá neliniární transformace (stupeň polynomu dva nebo tři).
Nutný počet vlícovacích bodů můžeme určit z následujícího vztahu
kde: p je stupeň polynomu transformace.
| p | nutný počet bodů | doporučený počet bodů |
|---|---|---|
| 1 | 3 | 4 |
| 2 | 6 | 10 |
| 3 | 10 | 15 |
V našem případě použijeme body mřížky označené čísly 1, 2, 3, 4 a 5.
| bod | Y | X |
|---|---|---|
| 1 | 830 000 | 1 000 000 |
| 2 | 770 000 | 1 000 000 |
| 3 | 960 000 | 960 000 |
| 4 | 960 000 | 960 000 |
| 5 | 800 000 | 980 000 |
Do lokace zod naimportujeme podkladovou rastrovou vrstvu obsahující mřížku.
wget http://gama.fsv.cvut.cz/~landa/vyuka/153YZOD/dataset/mrizka.tif # import pomocné podkladové rastrové vrstvy (lokace zod) # r.in.gdal -o in=mrizka.tif out=mrizka
wxGUI
wxGUI obsahuje nástroj pro georeferencování dostupný z menu File->Georectify. Na první obrazovce zvolíme lokaci obsahující souřadnicově nepřipojená data, vytvoříme či modifikujeme obrazovou skupinu, tak aby obsahovala všechny rastrové vrstvy určené ke georeferencování a nakonec rastrovou vrstvu, kterou bude vizualizována v mapovém okně.
 |
 |
 |
V mapovém okně aktuální (tj. cílové) a zdrojové lokace se určí identické body (viz obr. č.x).
CLI
Nejprve vytvoříme obrazovou skupinu (seskupuje obrazové vrstvy určené pro rektifikace) pomocí modulu i.group.
# založení obrazové skupiny a vložení jednotlivých pásem MSS # i.group group=mss input=mss1,mss2,mss3,mss4,mss1_mrizka
Dále nastavíme cílovou lokaci a mapset modulem i.target.
# nastavení cílové lokace a mapsetu # i.target group=mss location=zod mapset=student
Vlícovací body (jejich nutný počet závisí na zvoleném stupni polynomu transformace, viz tab. č.2) se definují pomocí modulu i.points (je-li podkladem rastrová vrstva) nebo i.vpoints (vektorová vrstva).
V lokaci zod-xy spustíme modul i.points
i.points
Z menu vybereme vrstvu 'mss1_mrizka' a grafické okno se rozdělí na čtyři rámce (viz obr. č.2), ze spodního menu vybereme "Plot Raster" a kliknutím do pravé části okna (tj. části odpovídající cílové lokaci) vybereme dvojklikem vrstvu 'mrizka'.
Pomocí funkce "Zoom" (v podmenu zvolíme "Box") zvětšíme detail prvního bodu jak ve zdrojové (levá horní část), tak v cílové lokaci (pravá horní část, viz obr. č.2) a označíme daný vlícovací bod nejprve ve zdrojové a posléze v cílové lokaci. Podobně pokračujeme u dalších vlícovacích bodů. Po označení všech požadovaných bodů provedeme statistickou analýzu - položka menu "Analyze". Nevhodně či chybně určený vlícovací bod lze vypnout kliknutím na tento bod v tabulce "Analyze". Příklad analýzy vlícovacích bodů:
LOCATION: zod-xy GROUP: mss MAPSET: PERMANENT
Analysis of control point registration
error image target
# col row target east north east north
1 0.0 -0.0 2.1 10.0 129.1 -829979.1 -999962.0
2 -0.0 -0.0 1.2 1062.1 129.0 -770004.0 -999964.1
3 0.0 -0.0 2.1 1062.0 831.0 -770004.2 -959950.0
4 -0.0 -0.0 1.2 10.0 831.0 -829977.8 -959953.4
5 -0.1 0.0 3.7 536.1 480.0 -799991.5 -979957.6
Overall rms error: 2.27
Kromě souřadnic bodů v obou souřadnicových systémech je důležitá jejich směrodatná odchylka, která by neměla přesáhnout polovinu požadovaného rozlišení v cílové lokaci.
Práci s modulem ukončíme výběrem položky menu "QUIT" a potvrzením otázky "really quit?" - "yes".
Po zvolení vlícovacích bodů provedeme rektifikaci dat modulem i.rectify.
# transformace dat (stupeň polynomu 1) # i.rectify group=mss ext= order=1
Po úspěšné rektifikaci se lokaci zod objeví georeferencované obrazové kanály LandSat-MSS. Ve výsledku máme k dispozici družicové snímky pořízené roku 1994 (LandSat5-TM), 2004 (LandSat7-ETM+) a z roku 1975 (LandSat1-MSS) - viz obr. č.3.
