Pro některé z vás to nebude asi nic nového, ale nedávno jsem narazil na zajímavou vlastnost knihovny GDAL, resp. na zajímavý formát vstupních rastrových dat: GDAL Virtual Format.
Asi víte, že GDAL je knihovna pro práci s rastrovými (a díky části nazývané jako OGR i vektorovými) daty. GDAL se používá ve většině open source produktech, ale i např. v Esri ArcGIS. GDAL momentálně podporuje 142 rastrových a 84 vektorových formátů, mezi kterými můžete data konvertovat. GDAL, kromě toho že jde o programátorskou knihovnu, nabízí po instalaci i celou řadu nástrojů příkazové řádky, pomocí kterých můžete zjišťovat metadata dat, transformovat mezi souřadnicovými systémy nebo dotazovat se na prvky v datových sadách. GDAL je zkrátka skvělá věc jak do velkých projektů, tak do malých skriptů. Když máte na serveru hromadu dat a chcete zjistit, co jsou zač, GDAL tento úkol hravě vyřeší.
Jedním z podporovaných rastrových formátů je, jak jsem zmínil na začátku, GDAL Virtual Format. Je to vlastně hrozně jednoduché: Vstupním souborem je XML dokument, který obsahuje metadata a linky na další vstupní (rastrové) soubory, které mohou být z různých typů zdrojů. Kdy a jak se to dá použít? Pokud máte nástroj, který jako vstup vyžaduje jeden rastrový soubor, ale vy jste si z nějakého zdroje stáhli souborů víc – letecké snímky apod., potřebujete použít “nástroj typu patch“. Jenomže to by na vám disku vznikl jeden obří soubor – vedle několika původních. A další nástroje by musely pracovat s tímto velkým souborem.
Pomocí GDAL Virtual Format umístíte do XMLka odkaz na původní soubory, které se ale z pohledu GDAL budou tvářit jako jeden datový zdroj.
Pojďme si to zkusit na datech ze SRTM – jedná se digitální model terénu v rozlišení cca 30m, distribuovaných jako soubory GeoTIFF v rastrových dlaždicích o velikosti 5° zeměpisné šířky a délky. Data si můžete vybrat a stáhnout pomocí nástroje http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp.
Pro 4 dlaždice přibližně z oblasti střední Evropy stačí použít nástroj wget:
$ wget http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_40_03.zip $ wget http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_40_02.zip $ wget http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_39_02.zip $ wget http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_39_03.zip
Data samozřejmě musíme rozbalit
$ for i in *.zip; do unzip -o $i; done
Na disku nyní máme sadu souborů
srtm_39_02.tif srtm_39_03.tif srtm_40_02.tif srtm_40_03.tif srtm_39_02.hdr srtm_39_03.hdr srtm_40_02.hdr srtm_40_03.hdr srtm_39_02.tfw srtm_39_03.tfw srtm_40_02.tfw srtm_40_03.tfw
Můžeme použít gdalinfo pro zjištění základních metadat
$ gdalinfo srtm_39_03.tif
Driver: GTiff/GeoTIFF
Files: srtm_39_03.tif
Size is 6001, 6001
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (9.999583623875717,50.000417029852201)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 9.9995836, 50.0004170) ( 9d59'58.50"E, 50d 0' 1.50"N)
Lower Left ( 9.9995836, 44.9995837) ( 9d59'58.50"E, 44d59'58.50"N)
Upper Right ( 15.0004170, 50.0004170) ( 15d 0' 1.50"E, 50d 0' 1.50"N)
Lower Right ( 15.0004170, 44.9995837) ( 15d 0' 1.50"E, 44d59'58.50"N)
Center ( 12.5000003, 47.5000004) ( 12d30' 0.00"E, 47d30' 0.00"N)
Band 1 Block=6001x1 Type=Int16, ColorInterp=Gray
NoData Value=-32768
Nyní máme představu o souřadnicovém systému, hraničních souřadnicích a charakteru dat.
Takto 4 stažené a rozbalené soubory mají dohromady 276M. Kdybychom je “slepili” dohromady, výsledný soubor by nebyl o mnoho menší. Namísto toho použijeme GDAL virtual raster.
$ gdalbuildvrt srtm.vrt srtm*.tif
Výstupní soubor srtm.vrt je XML soubor, obsahující metadata o souřadnicovém systému, hraničních souřadnicích celé datové sady, základní statistiky pro případnou tvorbu histogramu a samozřejmě odkazy na původní soubory.
Můžeme ho použít pro gdalinfo
$ gdalinfo srtm.vrt
Size is 12001, 12001
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (9.999583623875717,55.000417150924136)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left ( 9.9995836, 55.0004172) ( 9d59'58.50"E, 55d 0' 1.50"N)
Lower Left ( 9.9995836, 44.9995838) ( 9d59'58.50"E, 44d59'58.50"N)
Upper Right ( 20.0004170, 55.0004172) ( 20d 0' 1.50"E, 55d 0' 1.50"N)
Lower Right ( 20.0004170, 44.9995838) ( 20d 0' 1.50"E, 44d59'58.50"N)
Center ( 15.0000003, 50.0000005) ( 15d 0' 0.00"E, 50d 0' 0.00"N)
Band 1 Block=128x128 Type=Int16, ColorInterp=Gray
Min=-6.000 Max=720.000
Minimum=-6.000, Maximum=720.000, Mean=180.106, StdDev=149.815
NoData Value=-32768
Metadata:
STATISTICS_MAXIMUM=720
STATISTICS_MEAN=180.10611254643
STATISTICS_MINIMUM=-6
STATISTICS_STDDEV=149.81495511855
Vidíte, že velikost je opravdu 2×2 dlaždice původního souboru. Na disku nám vznikl pouze malý soubor, původní data zůstala zachovaná.
Tímhle způsobem můžete kombinovat nejenom původní rastry, ale i další “virtuální formáty”, jako WMS nebo WMTS o tom až někdy příště.
Na závěr, pokud jste dočetli až sem, můžete se na výsledek za odměnu podívat třeba v QGISu:
$ qgis srtm.vrt
To samé můžete vyzkoušet, pokud máte např. satelitní nebo letecké snímky.
Za GISMentors Jáchym Čepický