GDAL Virtual driver

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ý