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

Z GeoWikiCZ

Import, export a rektifikace dat

< Zpracování obrazových 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-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".

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.

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 (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.

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