153ZODH / 7. cvičení
Import, export a rektifikace dat
Osnova
Import a export georeferencovaných (souřadnicově připojených) dat. Georeferencování souřadnicově nepřipojených rastrových dat.
Importovaná data (výřez družicové scény LandSat1-MSS a LandSat7-ETM+) jsou volně šiřitelná a pocházejí z datasetu FreeGeodataCZ.
Seznam příkazů
- r.in.gdal
- g.region
- r.colors
- r.composite
- i.group
- i.target
- i.points
- i.vpoints
- i.rectify
- r.out.gdal
- r.resamp.interp
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. Import dat si ukážeme na příkladu: Tarball etm-jtsk.tar.gz obsahuje družicová data LandSat7-ETM+. Tento soubor nejprve dekomprimujeme a 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 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í rastrových dat
Modul r.external (až od verze GRASS 6.4) umožňuje připojení rastrových dat přes knihovnu GDAL bez nutnosti konverze dat do nativního rastrového formátu GRASSu.
# připojení rastrového souboru etm1.tif r.external -o input=etm1.tif out=etm1_link
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
Georeferencová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 georeferencová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 -eo 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á o data 8-bitová, postačí importovat pouze jeden z kanálů.
r.in.gdal input=mss1.tif out=mss1 band=1
Volba vlícovacích bodů, rektifikace dat
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í. Nakonec určíme rastrovou vrstvu, která 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. č.6) a provede rektifikace dat.
CLI
Nejprve vytvoříme obrazovou skupinu (seskupuje obrazové vrstvy určené pro rektifikaci) 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. č.1) 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'. 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. č.7) 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. č.8.
Export rastrových dat
Exportovaná rastrová vrstva převezme prostorový rozsah a rozlišení výpočetního regionu. Jako výchozí metoda převzorkování je použita metoda nejbližšího souseda. Alternativně lze využít modul r.resamp.interp, který podporuje bilineární a bikubickou interpolaci.
# nastavení regionu na základě družicové scény (nedojde k převzorkování) # g.region rast=etm1 # # export dat, formát GeoTIFF # r.out.gdal input=etm1 out=etm1.tif type=Byte format=GTiff
# výřez # g.region n=-957500 s=-1007318 w=-830529 e=-763855 res=50 # # export výřezu, formát PNG (převzorkování metodou nejbližšího souseda) # r.out.gdal input=etm1 out=etm1.png type=Byte format=PNG