Knihovna GDAL pro převod mezi formáty a souřadnicovými systémy

Programátorská knihovna GDAL (Geospatial Data Abstraction Library) je určena pro práci s množstvím rastrových i vektorových formátů používaných v GIS. GDAL je využíván celou řadou dalších programů jako základní knihovna (GRASS GIS, QGIS, …), dokonce i v proprietárním produktu ArcGIS.

Poznámka

V dřívějšich verzích byla tato knihovna rozdělena na dvě části. GDAL pracující s rastrovými daty a OGR pro vektorová data. Ve verzi 2.0 byly tyto dvě větve sloučeny. Stále však můžete narazit na označení části pro práci s vektory jako OGR.

Knihovna je šířena s několika konzolovými programy, které můžeme použít na celou řadu operací. Detailnější práci s knihovnou z pohledu programátora rozebíráme v části věnované programovacímu jazyku Python z pohledu GIS. Knihovna GDAL kromě toho nabízí mnoho užitečných příkazů pro práci s rastrovými (gdalinfo, gdalsrsinfo, gdalwarp, gdaltransform, gdal_translate, gdaldem, gdallocationinfo, gdalmanage, gdaladdo, gdal_contour, gdaltindex) nebo vektorovými (ogrinfo, ogrtindex, ogrlineref, ogr2ogr) daty.

Zde si představíme pouze některé příkazy, které jsou distribuovány spolu s knihovnou GDAL. Úplný seznam naleznete na https://gdal.org/programs/index.html.

Příkazy pro práci s rastrovými daty

gdalinfo

Příkaz gdalinfo umožňuje zobrazit některá metadat rastrových dat.

Zobrazení metadat z rastrového souboru z příkazové řádky

gdalinfo lsat7_2002_nir.tiff

Driver: GTiff/GeoTIFF
Files: lsat7_2002_nir.tiff
Size is 1287, 831
Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",36.16666666666666],
    PARAMETER["standard_parallel_2",34.33333333333334],
    PARAMETER["latitude_of_origin",33.75],
    PARAMETER["central_meridian",-79],
    PARAMETER["false_easting",609601.22],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (630540.000000000000000,226980.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  630540.000,  226980.000) ( 78d46' 6.04"W, 35d47'45.23"N)
Lower Left  (  630540.000,  218670.000) ( 78d46' 6.81"W, 35d43'15.59"N)
Upper Right (  643410.000,  226980.000) ( 78d37'33.46"W, 35d47'43.96"N)
Lower Right (  643410.000,  218670.000) ( 78d37'34.70"W, 35d43'14.31"N)
Center      (  636975.000,  222825.000) ( 78d41'50.25"W, 35d45'29.85"N)
Band 1 Block=1287x1 Type=Float32, ColorInterp=Gray
Band 2 Block=1287x1 Type=Float32, ColorInterp=Undefined
Band 3 Block=1287x1 Type=Float32, ColorInterp=Undefined

gdalsrsinfo

Pokud vám stačí pouze informace o použitém souřadnicovém systému, tak stačí použít příkaz gdalsrsinfo, který vrátí definici souřadnicového systému rastru ve formátu knihovny Proj.4 a v tzv. Well Known Text (WKT) notaci:

Zobrazení informace o souřadnicovém systému z příkazové řádky

gdalsrsinfo lsat7_2002_nir.tiff

PROJ.4 : '+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79
 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs '

 OGC WKT :
 PROJCS["Lambert Conformal Conic",
     GEOGCS["NAD83",
         DATUM["North_American_Datum_1983",
             SPHEROID["GRS 1980",6378137,298.2572221010002,
                 AUTHORITY["EPSG","7019"]],
             AUTHORITY["EPSG","6269"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4269"]],
     PROJECTION["Lambert_Conformal_Conic_2SP"],
     PARAMETER["standard_parallel_1",36.16666666666666],
     PARAMETER["standard_parallel_2",34.33333333333334],
     PARAMETER["latitude_of_origin",33.75],
     PARAMETER["central_meridian",-79],
     PARAMETER["false_easting",609601.22],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

gdalwarp

Asi nejpoužívanější příkaz je gdalwarp. Tento příkaz má dvě funkce: práce se souřadnicovými systémy rastrových dat a jejich transformace mezi jednotlivými formáty.

Podporované formáty zjistíte pomocí parametru –formats:

Podporované formáty knihovny gdal z příkazové řádky

gdalwarp --formats

Supported Formats:
  VRT (rw+v): Virtual Raster
  GTiff (rw+vs): GeoTIFF
  NITF (rw+vs): National Imagery Transmission Format
  RPFTOC (rovs): Raster Product Format TOC format
  ECRGTOC (rovs): ECRG TOC format
  HFA (rw+v): Erdas Imagine Images (.img)
  SAR_CEOS (rov): CEOS SAR Image
  CEOS (rov): CEOS Image
  JAXAPALSAR (rov): JAXA PALSAR Product Reader (Level 1.1/1.5)
  GFF (rov): Ground-based SAR Applications Testbed File Format (.gff)
  ELAS (rw+v): ELAS
  AIG (rov): Arc/Info Binary Grid
  AAIGrid (rwv): Arc/Info ASCII Grid
  GRASSASCIIGrid (rov): GRASS ASCII Grid
  SDTS (rov): SDTS Raster
  ...

Syntaxe programu gdalwarp (i u tohoto programu funguje parametr --help a určitě se podívejte na manuálovou stránku programu man gdalarp) je následující:

gdalwarp [PŘEPÍNAČE A VOLBY] zdrojový_soubor výstupní_soubor

Transformace rastru ve formátu GeoTIFF do formátu Windows Bitmap při zachování souřadnicového systému vypadá následovně:

Transformace GDAL z GeoTIFF do BMP z příkazové řádky

gdalwarp -of BMP lsat7_2002_nir.tiff lsat7_2002_nir.bmp

Creating output file that is 1287P x 831L.
ERROR 1: Attempt to create BMP dataset with an illegal
data type (Float32), only Byte supported by the format.

Vidíme, že formát BMP nepodporuje zdrojová data - číslo s plovoucí desetinnou čárkou. Datový typ nastavíme pomocí parametru -type (samozřejmě tak přijdeme o hodnoty mimo rozsah tohoto datového typu).

gdalwarp -of BMP -ot Byte lsat7_2002_nir.tiff lsat7_2002_nir.bmp

Creating output file that is 1287P x 831L.
Processing input file lsat7_2002_nir.tiff.
0...10...20...30...40...50...60...70...80...90...100 - done.
../_images/lsat7_2002_nir.png

Obr. 41 Výsledný obrázek převodu rastrové mapy na formát BMP.

Poznámka

Vedle souboru lsat7_2002_nir.bmp vytvořil GDAL také souboru lsat7_2002_nir.bmp.aux.xml obsahující metadata, mimo jiné i informace o souřadnicovém systému. Pokud tento soubor smažete nebo změníte jeho jméno, dostanete následující výstup, tj. bez informace o souřadnicovém systému.

Ověření výsledného souboru pomocí gdalinfo z příkazové řádky

gdalinfo lsat7_2002_nir.bmp

Driver: BMP/MS Windows Device Independent Bitmap
Files: lsat7_2002_nir.bmp
Size is 1287, 831
Coordinate System is `'
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  831.0)
Upper Right ( 1287.0,    0.0)
Lower Right ( 1287.0,  831.0)
Center      (  643.5,  415.5)
Band 1 Block=1287x1 Type=Byte, ColorInterp=Red
Band 2 Block=1287x1 Type=Byte, ColorInterp=Green
Band 3 Block=1287x1 Type=Byte, ColorInterp=Blue

Dalším obvyklým krokem je transformace při změně souřadnicového systému (v našem případě zůstane vstupní formát GeoTIFF zachován i na výstupu). Při transformacích můžeme použít 2 parametry pro popis souřadnicových systémů ve vztahu ke vstupní resp. výstupní rastrové mapě:

-s_srs

definice souř. systému vstupní dat (source)

-t_srs

definice souř. systému výstupní dat (target)

Tyto parametry mají větší prioritu při zpracování vstupních dat, než případná metadata v těchto datech přítomná.

Transformace rastrových dat do jiného souřadnicového systému z příkazové řádky

Souřadnicový systém vstupních dat je známý, v našem příkladě nastavíme pouze souřadnicový systém pro výstupní data. Zápis souřadnicového systému je totožný se zápisem pro knihovnu Proj.4. My použijeme kód EPSG:4326, což je souřadnicový systém WGS84.

gdalwarp -t_srs +init=epsg:4326 lsat7_2002_nir.tiff lsat7_2002_nir-wgs84.bmp

Creating output file that is 1359P x 717L.
Processing input file lsat7_2002_nir.tiff.
0...10...20...30...40...50...60...70...80...90...100 - done.
../_images/lsat7_2002_nir-wgs84.png

Obr. 42 Výsledek převodu rastrových dat do souřadnicového systému WGS84.

gdaltransform

Funguje podobně jako program cs2cs knihovny Proj4, tj. transformuje souřadnice mezi souřadnicovými systémy.

gdal_translate

Převádí rastrová data mezi různými formáty. Na rozdíl od gdalwarp neumožňuje data transformovat do jiného souřadnicového systému. Lze ale nastavit souřadnicový systém výstupních dat pomocí parametru -a_srs (kdy nechodází k transformaci dat, ale pouze nastavení souřadnicového systému do metadat výstupního souboru).

gdaldem

Nástroj gdaldem vám pomůže zanalyzovat a vizualizovat digitální modely reliéfu (DMR). Ze vstupního DMR lze vygenerovat

  • stínovaný reliéf,

  • mapu sklonu svahu,

  • mapu expozice,

  • barevný reliéf,

  • a další …

Vytvoření mapy stínového reliéfu ze vstupního rastrového souboru z příkazové řádky

gdaldem hillshade dem.tiff hillshade.tiff
../_images/hillshade.png

Obr. 43 Mapa stínovaného reliéfu vytvořená pomocí utility gdaldem.

gdallocationinfo

Nástroj gdallocationinfo se umožňuje ptát se na hodnoty rastrových dat o daných rastrových souřadnicích.

Dotaz na hodnotu rastru podle souřadnic z příkazové řádky

gdallocationinfo lsat7_2002_nir-wgs84.tiff 15 50

Report:
  Location: (15P,50L)
  Band 1:
    Value: 110
  Band 2:
    Value: 221
  Band 3:
    Value: 189

gdalmanage

Program gdalmanage umožňuje práci s rastrovými soubory na úrovni operačního systému, jejich identifikaci, přejmenování, mazání a kopírování.

Použití z příkazové řádky

Obsah pracovního adresáře může vypadat z pohledu GDAL následovně:

gdalmanage identify *

dem_srtm.tiff: GTiff
hillshade.bmp: BMP
hillshade.png: PNG
hillshade.tiff: GTiff
lsat7_2002_nir.bmp: BMP
lsat7_2002_nir.png: PNG
lsat7_2002_nir.tiff: GTiff
lsat7_2002_nir-wgs84.bmp: BMP
lsat7_2002_nir-wgs84.png: PNG
lsat7_2002_nir-wgs84.tiff: GTiff

Poznámka

gdalmanage lze použít pro případné změny a mazání více souborových formátů (např. *.tfw soubory).

gdaladdo

Nástroj gdaladdo umožňuje pracovat s tzv. pyramidami - zmenšenými kopiemi rastrových dat uložených přímo uvnitř anebo externě rastrového souboru. Ve výsledku bude práce s rastrem u malých měřítek výrazně rychlejší - vznikne v podstatě prostorový index rastrových dat (používá např. QGIS pro zobrazování rastrů).

Vytvoření přehledových pyramid rastrového souboru z příkazové řádky

# ověření velikosti původního souboru
ls -lh lsat7_2002_nir.tiff

-rw-rw-r-- 1 user user 13M apr 18 00:00 lsat7_2002_nir.tiff

# vytvoření pyramid
gdaladdo lsat7_2002_nir.tiff 2 4 8 16

# opětovné ověření velikosti změněného souboru
ls -lh lsat7_2002_nir.tiff

-rw-rw-r-- 1 user user 19M apr 18 00:00 lsat7_2002_nir.tiff

gdal_contour

Nástroj gdal_contour vytvoří vektorové vrstevnice ze vstupního digitálního modelu reliéfu.

Vytvoření vrstevnic z příkazové řádky

gdal_contour -a elev dem.tiff vrstevnice.shp -i 10.0
../_images/vrstevnice.png

Obr. 44 Získané (a obarvené) vrstevnice z DMT.

gdal_rasterize

Nástroj gdal_rasterize provede rasterizaci vektorových dat (tj. převede data z vektorové reprezentace do rastru).

Převod vektorových vrstevnic na rastrová data z příkazové řádky

Výstupní formát BMP, prostorové rozlišení 10m

gdal_rasterize -a elev -of GTiff -ot Byte -tr 10 10 vrstevnice.shp vrstevnice.tiff

gdaltindex

Vytvoří tzv. tile-index vektorový soubor obsahující obalový polygon (obdélník) okolo každého rastrového souboru. Tento prostorový index lze pak použít do dalších operací v prostředí GDAL, stejně tak jako vrstvu v programu mapserver.

Příkazy pro práci s vektorovými daty

ogrinfo

Sesterským programem ke gdalinfo je ogrinfo - vypíše dostupné informace o vektorových datech.

Poznámka

OGR pracuje na abstraktním datovém modelu

  • Zdroj (data source)
    • Vrstva (layer)

      • Vektorový objekt (feature)

kde

  • Zdrojem může být soubor, adresář nebo prostorová databáze

  • Vrstvou může být tabulka v databázi nebo vlastní data v souboru

Jsou-li data uložena v souboru, bývá název souboru a název vrstvy totožný.

Toto na první pohled možná lehce matoucí uspořádání je způsobeno tím, že GDAL (resp. vektorová část OGR) se snaží přistupovat ke všem možným datovým zdrojům, z nichž některé umožňují do zdroje (souboru, databáze, …) uložit více dat (vrstev, tabulek) a jiné ne.

Podrobnější informace o datovém modulu knihovny GDAL najdete ve školení GeoPython pro pokročilé.

Dotaz na metadata vektorového souboru z příkazové řádky

Necháme si vypsat informace o souboru vrstevnice.shp (pokud vynecháme parametr -so (summary only), vypíší se informace o každém vektorovém prvku):

ogrinfo vrstevnice.shp vrstevnice -so

INFO: Open of `vrstevnice.shp'
      using driver `ESRI Shapefile' successful.

Layer name: vrstevnice
Geometry: Line String
Feature Count: 175150
Extent: (-904049.056059, -1227170.827189) - (-431499.549460, -935327.979496)
Layer SRS WKT:
PROJCS["Krovak",
    GEOGCS["GCS_bessel",
        DATUM["Militar_Geographische_Institut",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.8],
    PARAMETER["azimuth",0],
    PARAMETER["pseudo_standard_parallel_1",0],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
ID: Integer (8.0)
elev: Real (12.3)

Vidíme, že vektorová data jsou v souřadnicovém systému S-JTSK, hraniční souřadnice jsou (-904049.056059, -1227170.827189) - (-431499.549460, -935327.979496) a atributová tabulka má 2 atributy: ID a elev (obsahující výšku nad mořem každé vrstevnice). Jedná se o soubor s liniovou geometrií.

ogrtindex

ogrtindex je sesterským programem k programu gdaltindex. Máte-li adresář plný vektorových dlaždic a chcete-li s nimy rychle pracovat, vytvoříte vektrový soubor s hranicemi těchto souborů a odkazem do adresářové struktury.

ogrlineref

ogrlineref slouží k tvorbě souboru obsahujícím segmenty o daných délek. Umožňuje získávat jejich souřadnice, vzdálenosti, staničení atd., to vše v lineární referenční síti.

ogr2ogr

Nástroj ogr2ogr je obdobou rastrového gdalwarp, který umožňuje transformaci vektorových dat.

Obecná syntaxe je:

ogr2ogr [VOLBY] výstupní_soubor vstupní_soubor

Stejně jako u gdalwarp, můžete podporované formáty vypsat pomocí parametru –formats:

ogr2ogr --formats

Supported Formats:
  -> "ESRI Shapefile" (read/write)
  -> "MapInfo File" (read/write)
  -> "UK .NTF" (readonly)
  -> "SDTS" (readonly)
  -> "TIGER" (read/write)
  -> "S57" (read/write)
  -> "DGN" (read/write)
  ...

Pro práci se souřadnicovými systémy opět můžeme použít některý z následujících parametrů:

  • -a_srs - přiřadí informaci o souřadnicovém systému do metadat výstupnímu souboru

  • -t_srs - provode transformaci dat do souřadnicového systému výstupních dat

  • -s_srs - nastaví souřadnicový systém vstupních dat

Tyto parametry jsou kompatibilní se zápisem pro knihovnu Proj4.

Převod souboru vrstevnic ve formátu Esri Shapefile na formát KML z příkazové řádky

ogr2ogr -f KML -t_srs epsg:4326 vrstevnice.kml vrstevnice.shp

Výsledný soubor můžeme zkontrolovat pomocí ogrinfo:

ogrinfo vrstevnice.kml vrstevnice -so

INFO: Open of `vrstevnice.kml'
      using driver `LIBKML' successful.

Layer name: vrstevnice
Geometry: Unknown (any)
Feature Count: 175150
Extent: (12.060792, 48.554130) - (18.825375, 51.055295)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4326"]]
 ...