153YZOD Zpracování obrazových dat 2006 - 7. cvičení
Import, export a rektifikace dat
osnona
V rámci tohoto kurzu bude objasněn princip importu a exportu rastrových dat v GRASSu a to jak georeferencovaných (tj. souřadnicově připojených), tak souřadnicově nepřipojených dat. Na modelovém příkladu bude ukázán proces georefencování dat v GRASSu.
Importovaná data (výřez družicové scény LandSat1-MSS a LandSat7-ETM+) jsou volně šířitelná a pocházejí z datasetu Geografická data ČR pro GRASS.
seznam příkazů
- r.in.gdal
- g.region
- r.colors
- d.rast
- r.composite
- g.remove
- i.group
- i.target
- i.points
- i.vpoints
- i.rectify
- r.out.gdal
rastrové datové 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 či BIL/BSQ. Obecně můžeme rozdělit datové formáty pro obrazová data na dva základní typy:
- datové formáty, které obsahují jedno pásmo (kanál) na jeden soubor (GeoTIFF, PNG, SUN raster formát, atd.)
- 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 vstupující do GISu - rastrová a vektorová data) mohou být již georeferencovaná (tj. souřadnicově přípojená) či souřadnicově nepřipojená. Tento fakt velmi významně ovlivňuje proces jejich importu do GRASSu (koneckonců do GISu jako takového). 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 rastrový soubor 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.
import georeferencovaných obrazových dat
Proces importu souřadnicově připojených dat je poměrně jednoduchý. V podstatě je nutné pouze určit v jakém datovém formátu jsou vstupní data uložena a podle toho vybrat korespondujícího modul pro import dat do GRASSu. Moduly pro import rastrových dat začínají předponou r.in.*. Tak například ve verzi 6.0 GRASS nabízí tyto moduly:
r.in.arc r.in.bin r.in.gdal r.in.mat r.in.srtm r.in.ascii r.info r.in.gridatb r.in.poly
V drtivé většině případů obstará import modul r.in.gdal, který pro svoji práci používá již výše zmíněnou knihovnu GDAL. Ta podporuje skutečně velký počet rastrových a vektorových formátů (a to jak v režimu čtení, tak zápisu) používaných ve většině GIS produktů.
Před vlastním importem je nutno založit novou GRASS lokaci a mapset. Vzhledem k náplni a časové dispozici tohoto kurzu tuto problematiku budeme ignorovat a pouze odkážeme na příslušný text v překladu GIS GRASS - Praktická rukověť, resp. GIS GRASS 6.0 - Praktická rukověť začínajících uživatelů. Pouze na okraj můžeme poznamenat, že modul r.in.gdal umožňuje automatické vytvoření lokace právě na základě importovaných souřadnicově připojených dat.
Spustíme tedy GRASS a vybereme lokaci/mapset sevcech/student. Stáhneme do pracovního adresáře (/home/student/work/) soubor etm-jtsk.tar.gz. Tento soubor dekomprimujeme, rozbalíme a posléze se přepneme do nově vzniklého adresáře etm-jtsk, který obsahuje jednotlivá pásma LandSat7-ETM+.
$ tar xvzf etm-jtsk.tar.gz $ cd etm-jtsk
Nejprve si ověříme, zda se jedná skutečně o data georefencovaná.
$ gdalinfo etm1.tif | head -n 20 Driver: GTiff/GeoTIFF Size is 2223, 1661 Coordinate System is `' Origin = (-830529.000000,-957500.000000) Pixel Size = (29.99277544,-29.99277544) 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 Metadata: COLOR_TABLE_RULES_COUNT=47 COLOR_TABLE_RULE_RGB_0=2.100000e+01 6.000000e+01 0 0 0 0 0 0 COLOR_TABLE_RULE_RGB_1=6.100000e+01 6.100000e+01 1 1 1 1 1 1 COLOR_TABLE_RULE_RGB_2=6.200000e+01 6.200000e+01 2 2 2 2 2 2 COLOR_TABLE_RULE_RGB_3=6.300000e+01 6.300000e+01 6 6 6 6 6 6 COLOR_TABLE_RULE_RGB_4=6.400000e+01 6.400000e+01 11 11 11 11 11 11
Poznámka: Pomocí systémové ulitky head vypíšeme pouze prvních 20 řádků výstupu z gdalinfo.
Z výše uvedeného je zřejmé, že jsou data georefencovaná a to v souřadnicovém systému S-JTSK.
Samotný import provedeme následujícím příkazem (zadáme název souboru obsahující vstupní data a název vzniknuvší vrstvy v databance GRASSu; za zmínku stojí přepínač -o, který ignoruje nastavení souřadnicového systému lokace GRASSu).
#import souboru ve formátu GeoTIFF # GRASS > r.in.gdal -o input=etm1.tif out=etm1 title="prvni pasmo LandSat7 ETM+"
Před zobrazením importované obrazové vrstvy nastavíme správný region a změníme tabulku barev na "grey.eq".
#nastavení regionu a tabulky barev # GRASS > g.region rast=etm1 GRASS > r.colors map=etm1 color=grey.eq GRASS > d.rast etm1
Stejným způsobem naimportujeme do aktuálního mapsetu zbývající soubory. Mapset student se nám tak rozroste o dalších 9 obrazových vrstev - 8 pásem LandSat7-ETM+ (2 vrstvy panchromatického pásma).
Hromadný import lze provést velmi jednoduše pomocí primitivního skriptu (pro BASH), např.:
#hromadný import dat ve formátu GeoTIFF
#
GRASS > for tiff in `ls *.tif`;do \
r.in.gdal -o input=$tiff out=${tiff%%.tif} title="LandSat7 ETM+";\
g.region rast=${tiff%%.tif};\
r.colors map=${tiff%%.tif} color=grey.eq;\
done
georeferencování obrazových dat
Nejprve stáhneme soubor mss-xy.tar.gz do našeho pracovního adresáře, soubor rozbalíme a následně se přepneme do vzniknuvšího adresáře mss-xy.
$ tar xvzf mss-xy.tar.gz $ cd mss-xy
Pomocí gdalinfo se přesvědčíme, že se skutečně jedná o souřadnicově nepřipojená data.
$ 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
Proces souřadnicového připojení (tzv. georeferencování) má tři základní fáze:
- 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 požadovaným souřadnicovým systémem. První lokaci označíme jako zdrojovou, druhou jako cílovou.
- transformace dat do cílové lokace
založení lokace XY a import dat
Spustíme GRASS a zadáme název zakládané lokace (např. sevcech-xy). Pokud spouštíme GRASS v grafickém módu, klikneme na "Create New Location". Mapset zvolíme základní - PERMANENT. První a druhou otázku potvrdíme - ENTER.
Would you like to create location <sevcech-xy> ? (y/n) [y] Do you have all this information? (y/n) [y]
Poté určíme souřadnicový systém naší lokace, v tomto případě zvolíme matematický - XY. Napíšeme "A" a potvrdíme.
Please specify the coordinate system for location <sevcech-xy> A x,y B Latitude-Longitude C UTM D Other Projection RETURN to cancel > A
Otázku
x,y coordinate system? (y/n) [y]
opět potvrdíme.
V dalším kroku zadáme krátký popis této lokace a následující otázku opět potvrdíme.
Please enter a one line description for location <sevcech-xy> > lokace xy ===================================================== lokace xy ===================================================== ok? (y/n) [y]
Následující formulář, viz obr. č.1 (nastavení výchozího regionu) můžeme přeskočit - ESC-ENTER. Region nastavéme později pomocí g.region.
Otázku
Do you accept this region? (y/n) [y] >
opět potvrdíme, podobně reagujeme na
Hit RETURN >
a poté povtrdíme již dobře známý formulář pro výběr lokace a mapsetu. GRASS se spustí s právě vytvořenou lokací.
Nyní provedeme import obrazových dat do této lokace. Použijeme již dobře známý modul r.in.gdal (přepínač -e má za následek rozšíření aktivního regionu podle importované vrstvy).
#import dat ve formátu TIFF # GRASS > r.in.gdal -e in=mss1.tif out=mss1
Takto se vytvoří tři obrazové vrstvy (mss1.red, mss1.green a mss1.blue), které složíme do výsledné vrstvy pomocí modulu r.composite. Nakonec můžeme již nepotřebné rastrové vrstvy odstranit.
#nutno nastavit nejdříve region podle aktuální vrstvy # GRASS > g.region rast=mss1.red # #složení RGB kanálů # GRASS > r.composite red=mss1.red green=mss1.green blue=mss1.blue out=mss1 # #odstranění nepotřebných vrstev # GRASS > g.remove rast=mss1.red,mss1.green,mss1.blue
nebo alternativně provést tyto operace hromadně pro všechny soubory ve formátu TIFF (dostupných v aktuální adresáři).
GRASS > for tiff in `ls *.tif`;do \
r.in.gdal -e in=$tiff out=${tiff%%.tif};\
g.region rast=${tiff%%.tif}.red;\
r.composite red=${tiff%%.tif}.red green=${tiff%%.tif}.green \
blue=${tiff%%.tif}.blue out=${tiff%%.tif};\
g.remove rast=${tiff%%.tif}.red,${tiff%%.tif}.green,${tiff%%.tif}.blue;\
done
volba vlícovacích bodů
Obecně můžeme rozlišit čtyři případy rektifikace obrazu:
- obraz - obraz (matematický souřadnicový systém, souřadnicově nepřipojená data)
- obraz - rastrová mapa (pokladem pro georeferencování je rastrová mapa)
- obraz - vektorová mapa (pokladem pro georeferencování je vektorová mapa)
- obraz - souřadnice (souřadnice v cílovém souřadnicovém systému jsou zadávány přímo z klávesnice)
Nejprve musíme vytvořit obrazovou skupinu, která bude obsahovat všechny obrazové vrstvy určené pro transformaci do lokace sevcech. K tomu nám poslouží modul i.group, skupinu nazveme mss a vložíme do ní jednotlivá pásma LandSat1-MSS.
#založení obrazové skupiny a vložení jednotlivých pásem MSS # GRASS > i.group group=mss input=mss1,mss2,mss3,mss4,mss1_mrizka
Před vlastní volbou vlícovacích bodů ještě nastavíme cílovou lokaci a mapset.
#nastavení cílové lokace a mapsetu # GRASS > i.target group=mss location=sevcech mapset=student
Vlícovací body (jejich nutný počet závisí na zvoleném stupni polynomu transformace, viz tab. č.2) se určí pomocí modulu i.points (je-li podkladem rastrová mapa) nebo i.vpoints (v případě vektorové podkladové mapy).
Z časových důvodů se uchýlíme k silně idealizovanému případu, kdy máme navíc k dispozici první pásmo LandSat1-MSS pokrytého geografickou mřížkou v souřadnicovém systému S-JTSK, 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 |
Mezitím ještě do lokace sevcech (mapset student) naimportujeme rastrovou vrstvu obsahující tuto mřížku a následně ji použijeme jako fiktivní podkladovou rastrovou vrstvu.
#import fiktivní podkladové rastrové vrstvy do sevcech/student # GRASS > r.in.gdal -o in=mrizka.tif out=mrizka
Poznámka: Současně můžeme spustit hned několik sezení GRASSu, nelze však spustit GRASS současně se stejnou lokací a mapsetem.
V lokaci sevcech-xy spustíme modul i.points a zadáme název obrazové skupiny.
GRASS > i.points
Z menu vybereme vrstvu mss1_mrizka (dvojklik myší). 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. Jako vstupní metody ("input method") se nám nabízí "keyboard" (vstup souřadnic z klavesnice) nebo "screen" (určení souřadnic z podkladové mapy - to je náš případ).
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). V tento okamžik můžeme označit bod nejprve ve zdrojové a posléze i korespondující bod v cílové lokaci. Vybraný vlícovací bod se nám označí zelenou značkou. 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". Takto deaktivovaný bod bude posléze označen červenou značkou. Příklad analýzy vlícovacích bodů:
LOCATION: sevcech-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".
transformace dat do cílové lokace
Samotný proces transformace je velmi jednoduchý, zvolíme název obrazové skupiny obsahující vrstvy určené pro transformaci, přívlastek pro název vrstev v cílové lokaci a stupeň polynomu transformace.
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 |
Pro obrazová data se obvykle používá neliniární transformace (stupeň polynomu dva nebo tři).
#transformace dat (stupeň polynomu 1) # GRASS > i.rectify group=mss ext= order=1
Spustíme-li GRASS s lokací sevcech a mapsetem student, mělo by se v seznamu rastrových vrstev objevit pět nových přírustků (vrstvu mss1_mrizka můžeme posléze odstranit).
Máme tak dispozici družicové snímky pořízené roku 1994 (LandSat5-TM), 2004 (LandSat7-ETM+) a konečně i roku 1975 (LandSat1-MSS) - viz obr. č.3.
export obrazových dat z GRASSu
Problematiky exportu rastrových dat z GRASSu se dotkneme pouze okrajově. Za upozornění stojí fakt, že je exportován pouze výřez vrstvy odpovídající nastavení aktuálního regionu.
Uvedeme si jeden ukázkový příklad - export prvního pásma LandSat7-ETM+ do souboru ve formátu GeoTIFF.
#nastavení regionu podle celkové družicové scény # GRASS > g.region rast=etm1 # #export dat do formátu GeoTIFF # GRASS > r.out.gdal input=etm1 out=etm1.tif type=Int16 format=GTiff