153ZODH / 7. cvičení

Z GeoWikiCZ

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ů

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

Obr. č.1: Hromadný import rastrových dat ve wxGUI

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ů:

  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 georeferencování (rastrová či vektorová mapa) uložené v lokaci s cílovým souřadnicovým systémem.
  3. 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 location modulu r.in.gdal
r.in.gdal input=mss1.tif out=mss1 location=zod-xy
  • Z uvítací obrazovky spustit "Location Wizard"
Obr. č.2: wxGUI 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:

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

Tab č.1: 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

V našem případě použijeme body mřížky označené čísly 1, 2, 3, 4 a 5.

Tab. č.2: 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

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

Obr č.3: Výběr zdrojové lokace
Obr č.4: Definice obrazové skupiny obsahující souřadnicově nepřipojená data
Obr č.5: Rastrová vrstva pro vizualizaci 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.

Obr č.6: Volba identických bodů a 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.

Obr. č.7: 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".

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.

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

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