153YZOD Zpracování obrazových dat 2006 - 7. cvičení

Z GeoWikiCZ

[ Zpracování obrazových dat ]

Import, export a rektifikace dat

osnova

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-MSSLandSat7-ETM+) jsou volně šířitelná a pocházejí z datasetu Geografická data ČR pro GRASS.

seznam příkazů

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ýchvektorový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 lokacimapset. 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:

  1. založení lokace XY (tedy lokace s matematickým souřadnicovým systémem) a import vstupních dat do této lokace
  2. 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. První lokaci označíme jako zdrojovou,  druhou jako cílovou.
  3. 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.

Obr. č.1: Formulář pro nastavení výchozího regionu

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.greenmss1.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ím 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:

  1. obraz - obraz (matematický souřadnicový systém, souřadnicově nepřipojená data)
  2. obraz - rastrová mapa (pokladem pro georeferencování je rastrová mapa)
  3. obraz - vektorová mapa (pokladem pro georeferencování je vektorová mapa)
  4. 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.

Tab. č.1: Vlícovací body v S-JTSK
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.

Obr. č.2: Volba vlícovacích bodů v grafickém okně GRASSu

Práci s modulem ukončíme výběrem položky menu "QUIT" a potvrzením otázky "really quit?" - "yes".

Nutný počet vlícovacích bodů můžeme určit z následujícího vztahu

kde: p je stupeň polynomu transformace.

Tab č.2: Počet vlícovacích bodů při volbě stupně polynomu transformace p
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 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.

#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řírůstků (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.

Obr. č.3: RGB barevná syntéza 243: MSS, TM, ETM+

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


[ Zpracování obrazových dat ]