Přeskočit obsah

09.2 GDAL

Knihovna GDAL umožňuje převádět data mezi různými rastrovými a vektorovými formáty. Seznam podporovaných rastrových a vektorových formátů je úctyhodný.

Knihovna je napsána v C++. Nabízí ale také API pro další programovací jazyky včetně Pythonu (dokumentace). Protože jsou však tato rozhraní vytvořena za využití generátoru SWIG, práce s těmito rozhraními občas bývá od úzu daného jazyka odlišná.

Součástí instalace je řada užitečných konzolových nástrojů.

Knihovna GDAL přistupuje k datovým zdrojům (soubor, adresář, databáze, protokol) abstraktně. Datový model pro rastrové zdroje dat:

Datový model pro vektorové zdroje dat:

Konzolové nástroje

Konzolové nástroje ve verzi GDAL 3.11 prošly kompletním přepsáním a konsolidací (návod na migraci). V našich materiálech se zaměříme na nové konzolové nástroje - viz dokumentace. Tradiční zápis používaný v dřívějších verzích knihovny uvádíme v poznámce. Je nutné ale poznamenat, že tradiční zápis příkazů funguje i v novějších verzích knihovny GDAL.

Note

Tradiční konzolové nástroje knihovny GDAL.

Rastrová data:

  • gdalinfo: Vypíše informace o rastrovém datasetu
  • gdalsrsinfo: Vypíše informace o souřadnicových systémech
  • gdal_translate: Převede rastrová data mezi zvolenými datovými formáty
  • gdalwarp: Transformace rastrových dat mezi souřadnicovými systémy
  • gdallocationinfo: Nástroj pro dotazování
  • gdaladdo: Vytváří nebo obnovuje přehledové obrazy
  • gdaldem: Nástroj pro analýzu a vizualizaci DEM
  • gdal_contour: Vytvoří vrstevnice z digitálního výškového modelu
  • gdal_rasterize: Rasterizace vektorových dat
  • gdalmanage: Nástroj pro správu rastrových souborů
  • gdal2tiles.py: Generuje adresář s dlaždicemi
  • gdal2xyz.py: Převádí rastrový soubor do xyz
  • gdal_calc.py: Rastorový kalkulátor
  • gdal_create: Vytvoří rastrový soubor
  • gdal_edit.py: Editace informací rastrového datasetu
  • gdal_fillnodata.py: Vyplní rastrové hodnoty
  • gdal_footprint: Vypočte zájmové území rastru
  • gdal_grid: Vytvoří pravidelnou mřížku hodnot z rozptýlených dat
  • gdal_merge.py: Vytvoří mozaiku z rastrových souborů
  • gdal_pansharpen.py: Provede operaci ostření
  • gdal_polygonize.py: Polygonizuje rastrová data
  • gdal_proximity.py: Vytvoří rastrovou mapu blízkosti
  • gdal_retile.py: Přepočte sada dlaždic a/nebo úrovní pyramidy z dlaždic
  • gdal_sieve.py: Odstraní malé rastrové plochy
  • gdal_viewshed: Vypočte masku viditelnosti
  • gdalattachpct.py: Připojí tabulku barev k rastrovému souboru
  • gdalbuildvrt: Vytvoří virtualní rastr ze vstupních souborů
  • gdalcompare.py: Porovná dva rastrové soubory
  • gdalmove.py: Provede změnu georeference rastrových souborů
  • gdaltindex: Vytvoří vektorový dataset z indexu dlaždic
  • gdaltransform: Transformace mezi souřadnicovými systémy
  • nearblack: Převede hodnoty bílé a černé barvy na hranách mřížky hodnot
  • pct2rgb.py: Převede 8bitová data na 24bitový RGB
  • rgb2pct.py: Převede 24bitový RGB na 8bitovou paletu
  • gdalmdiminfo: Vypíše strukturu vícerozměrného datasetu
  • gdalmdimtranslate: Převádí vícerozměrný dataset mezi různými formáty

Vektorová data:

  • ogrinfo: Vypíše informace o vektorovém datasetu
  • ogr2ogr: Převádí vektorové datasety mezi různými souřadnicovými systémy
  • ogrtindex: Vytvoří index dlaždic
  • ogrlineref: Vytvoří lineární referenční systém
  • ogrmerge.py: Spojí různé vektorové datasety do jednoho
  • ogr_layer_algebra.py: Provede operace mapové algebry na vektorových datech

Pro vypsání nápovědy lze použít přepínač --help.

gdal --help
gdal raster --help
gdal raster info --help

Ve výchozím nastavení nástroje vypisují informace ve datové struktuře JSON. Pokud prefererujete formátováný textový výstup tak použijte volbu --format text.

Rastrová data

  • gdal raster info: Vypíše informace o rastrovém datasetu 
    gdal raster info /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF

Tip

Lze zkratit na gdal info.

  • gdal raster convert: Převede rastrová data mezi zvolenými datovými formáty 
    # gdal raster info --formats
gdal raster convert --of GPKG /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF eu_dem_v11_E40N10.gpkg

Tip

Použijte přepínač --progress pro zobrazení postupu operace.

  • gdal raster clip: Umožňuje pracovat s částí rastrového datového zdroje 

    gdal raster clip --bbox 4445115,1510353,4896327,1734556 /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF eu_dem_v11_E40N10_REG.TIF

  • gdal raster reproject: Transformace rastrových dat mezi souřadnicovými systémy 

    gdal raster reproject --dst-crs EPSG:4326 eu_dem_v11_E40N10_REG.TIF eu_dem_v11_E40N10_WGS_84.TIF

  • gdal raster pixel-info: Nástroj pro dotazování 

    gdal raster pixel-info eu_dem_v11_E40N10_REG.TIF  11089 3843

    gdal raster pixel-info --position-crs dataset eu_dem_v11_E40N10_REG.TIF 4722347 1638498

  • gdal raster overview: Vytváří nebo obnovuje přehledové obrazy 

    # cp /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF .
# ls -l eu_dem_v11_E40N10.TIF
gdal raster overview add --levels=2,4,8,16 eu_dem_v11_E40N10.TIF

Úkol

Vyzkoušejte porovnat zrychlení načítání dat v QGIS.

  • gdal raster hillshade: Nástroj pro analýzu a vizualizaci DEM 

    gdal raster hillshade eu_dem_v11_E40N10_REG.TIF eu_dem_v11_E40N10_REG_HS.TIF

  • gdal raster contour: Vytvoří vrstevnice z digitálního výškového modelu 

    gdal raster contour --elevation-name elev --interval 10.0 eu_dem_v11_E40N10_REG.TIF vrstevnice.shp

  • gdal raster pipeline: Umožňuje sestavit vlastní procesní linku z jednotlivých výpočetních kroků 

    # hillshade + transformace
gdal raster pipeline "read eu_dem_v11_E40N10_REG.TIF ! hillshade --altitude 55 ! reproject --dst-crs EPSG:4326 ! write  eu_dem_v11_E40N10_REG_HS_WGS84.TIF"

Vektorová data

  • gdal vector info: Vypíše informace o vektorovém datasetu 
    gdal vector info vrstevnice.shp

    # datový zdroj s více vrstvami
gdal vector info /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg

    # seznam vrstev
# porovnej s `ogrinfo /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg`
# pripadne s `gdal vector info --format text /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg | grep 'Layer name'`
# tip: ve verzi GDAL 3.12+ existuje volba --summary
gdal vector info /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg | jq -r '.layers[].name'

    gdal vector info --layer U2018_CLC2018_V2020_20u1 /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg

    # přistup k prvkům
gdal vector info --layer U2018_CLC2018_V2020_20u1 --where "OBJECTID = 1" /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg 

    # počet polygonů spadající do skupiny jehličnatý les
gdal vector info --layer U2018_CLC2018_V2020_20u1 --where "Code_18 = 312" /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg | grep featureCount

Poznámka

Nomenklatura pro atribut Code_18

  • gdal vector convert: Převádí vektorové datasety mezi různými souřadnicovými systémy 

    gdal vector convert --format "ESRI Shapefile" --layer U2018_CLC2018_V2020_20u1_FR_REU /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg clc2018_fr_reu.shp

  • gdal vector filter: Umožňuje filtrovat data na základě atributů anebo prostorového rozsahu 

    # atributový filter
gdal vector filter --layer U2018_CLC2018_V2020_20u1 --where "Code_18 = 312" /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg clc2018_312.shp

    # kombinace atributového a prostorového filtru
gdal vector filter --layer U2018_CLC2018_V2020_20u1 --where "Code_18 = 312" --bbox 4445115,1510353,4896327,1734556 /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg clc2018_312_reg.shp

  • gdal vector pipeline: Umožňuje sestavit vlastní procesní linku z jednotlivých výpočetních kroků 

    # kombinace atributového a prostorového filtru + transformace
gdal vector pipeline ! read /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg --layer U2018_CLC2018_V2020_20u1 ! filter --where "Code_18 = 312" --bbox 4445115,1510353,4896327,1734556 ! reproject --dst-crs=EPSG:4326 ! write clc2018_312_reg_wgs84.shp

  • gdal vector rasterize: Rasterizace vektorových dat 

    gdal vector rasterize --attribute-name elev --format GTiff --datatype Byte --resolution 10,10 vrstevnice.shp vrstevnice.tiff

Různé

  • gdalmanage: Nástroj pro správu datových zdrojů 
    gdalmanage identify /mnt/repository/155FGIS/03/*

Tip

Ve verzi GDAL 3.12 je tento nástroj nahrazen gdal dataset identify

Síťové analýzy

  • gnmmanage: Správa síťového datasetu
  • gnmanalyse: Provede síťovou analýzu

Python API

Příklad práce s rastrovým datasetem:

from osgeo import gdal

with gdal.Open('/mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF') as ds:
    band = ds.GetRasterBand(1)
    stats = band.GetStatistics(True, True)
    print("minimum={:.3f}, maximum={:.3f}, mean={:.3f}, stddev={:.3f}".format(
          stats[0], stats[1], stats[2], stats[3]))

Příklad práce s vektorovým datasetem:

from osgeo import ogr

with ogr.Open('/mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg') as ds:
    for idx in range(ds.GetLayerCount()):
        lyr = ds.GetLayer(idx)
        print(lyr.GetName(), ":", lyr.GetFeatureCount())

Příklad konverze dat:

# gdal vector pipeline ! read /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg --layer U2018_CLC2018_V2020_20u1 ! filter --where "Code_18 = 312" ! write clc2018_312.shp

from osgeo import gdal
gdal.VectorTranslate(destNameOrDestDS='clc2018_312.shp',
                     srcDS='/mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg',
                     layers=['U2018_CLC2018_V2020_20u1'],
                     format='ESRI Shapefile',
                     where="Code_18 = 312")

Tip

Ve verzi GDAL 3.12 bylo rozšířeno Python API o možnost spouštět konzolové nástroje přímo.

from osgeo import gdal
with gdal.alg.pipeline(pipeline='read /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg --layer U2018_CLC2018_V2020_20u1 ! filter --where "Code_18 = 312" ! write clc2018_312.shp') as alg:
    ds = alg.Output()

Další materiály:

Příklady

Převod souborů SHP do GPKG

Nejprve stáhněme vstupní Data50:

wget https://openzu.cuzk.cz/opendata/Data50/epsg-5514/data50.zip

Vypišme základní udáje vstupního datového zdroje:

gdal vector info --format text /vsizip/data50.zip/shp | grep 'Layer name'

Tip

gdal vsi list -R /vsizip/data50.zip

# nebo `gdal vector info /vsizip/data50.zip/shp/Zamek.shp`
gdal vector info --layer Zamek /vsizip/data50.zip/shp/

Převod z formátu Esri Shapefile do OGC GeoPackage se provedeme jednoduše pomocí gdal vector convert:

gdal vector convert --format GPKG /vsizip/data50.zip/shp data50.gpkg

Úkol

Zkuste opravit chyby detekované při převodu z Esri Shapefile do OGC GeoPackage.

Převést můžeme také pouze vybrané vrstvy:

gdal vector pipeline ! read /vsizip/data50.zip/shp --layer Zamek --layer Hrad ! reproject --dst-crs EPSG:4326 ! write data50_zamek_hrad_wgs84.gpkg

Převod dlaždic do COG

Nejprve projdeme vstupní datatové soubory:

gdalmanage identify /mnt/repository/155FGIS/09/raster_tiles/*

gdalmanage identify /mnt/repository/155FGIS/09/raster_tiles/* | wc -l

Dále na základě vstupních dlaždic vytvoříme virtualní rastrový dataset:

gdal raster mosaic /mnt/repository/155FGIS/09/raster_tiles/*.tif osm_lc.vrt

gdal raster info --no-ct osm_lc.vrt

Tento virtualní rastrový dataset převedem do formátu COG:

gdal raster convert --format COG osm_lc.vrt osm_lc.tif

Poznámka

Vytvoření výstupního souboru o velikosti 1.2GB (255000x150000) zabere několik minut.

gdal raster info --no-ct /vsicurl/http://gislab.fsv.cvut.cz/barrel/osm_lc.tif

Úkol

Zkuste načíst rastrový soubor přes protokol v QGISu:

Úkol

Nahraďte hodnotu 110 za no-data.

Stahování dat pomocí ČÚZK Atom

V našem příkladě se pokusíme stáhnouti ortofoto data z ČÚZK prostřednictvím stahovací služby ATOM dostupné na adrese https://atom.cuzk.cz/ . Toho můžeme pomocí knihovny GDAL docíliti buď použitím konzolových nástrojů, nebo použitím rozhraní pro Python.

Konzolové nástroje

Na stránce stahovací služby ATOM najděme odkaz odpovídající ortofoto snímkům a dotažme se na jeho vrstvy za pomoci příkazu gdal vector info:

gdal vector info https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO

Tip

Pokud příkaz končí chybou ERROR 1: Failed to connect to atom.cuzk.cz port 443 after XXX ms: Connection timed out, upravte příkaz tak, aby specifikoval dobu dotazu: gdal vector info --config GDAL_HTTP_CONNECTTIMEOUT=666666 https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO

Podívejme se tedy na data obsažená v jediné vrácené vrstvě (georss):

gdal vector info --layer georss https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO

Tento příkaz už nám vrátí jednotlivé prvky vrstvy:

...

OGRFeature(georss):16299
  id (String) = https://atom.cuzk.cz/ORTOFOTO/datasetFeeds/CZ-00025712-CUZK_ORTOFOTO_WRTO24.2023.TREN42.xml
  title (String) = Ortofoto České republiky - mapový list: Třeboň 4-2
  updated (DateTime) = 2024/02/14 08:52:51+02
  author_name (String) = Český úřad zeměměřický a katastrální
  author_uri (String) = https://geoportal.cuzk.cz
  author_email (String) = cuzk.helpdesk@cuzk.cz
  rights (String) = žádné podmínky neplatí
  link_href (String) = https://geoportal.cuzk.cz/SDIProCSW/service.svc/get?REQUEST=GetRecordById&SERVICE=CSW&VERSION=2.0.2&OUTPUTSCHEMA=http://www.isotc211.org/2005/gmd&ELEMENTSETNAME=full&Id=CZ-00025712-CUZK_ATOM-MD_ORTOFOTO
  link_title (String) = metadata datasetu
  link_type (String) = application/xml
  link_rel (String) = describedby
  link2_href (String) = https://atom.cuzk.cz/ORTOFOTO/datasetFeeds/CZ-00025712-CUZK_ORTOFOTO_WRTO24.2023.TREN42.xml
  link2_rel (String) = alternate
  link2_title (String) = dataset feed
  link2_type (String) = application/atom+xml
  category_label (String) = S-JTSK
  category_term (String) = http://www.opengis.net/def/crs/EPSG/0/5514
  inspire_dls_spatial_dataset_identifier_code (String) = CZ-00025712-CUZK_ORTOFOTO_WRTO24.2023.TREN42
  inspire_dls_spatial_dataset_identifier_namespace (String) = ČÚZK
  POLYGON ((14.7186 48.998,14.7561 48.998,14.7561 49.0188,14.7186 49.0188,14.7186 48.998))

OGRFeature(georss):16300
...

V našem případě ovšem nepotřebujeme všechny prvky, ale vybereme si pouze jeden zájmový. Filtrovat můžeme např. pomocí atributu -where. Vyfiltrujme si pouze mapový list Praha 4-0:

gdal vector info --layer georss --where "title LIKE '%Praha 4-0'" https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO

Tip

Pokud chcete získat pouze odkaz link2_href, můžete jej z výstupu gdal vector info vyfiltrovat pomocí gdal vector info --format text --features --layer georss --where "title LIKE '%Praha 4-0'" https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO | grep -o 'https://atom.*xml'

Informace o produktu samotném se skrývají uvnitř pole link2_href. Bohužel, tento odkaz již nelze čísti přes nástroje GDAL, a musíme se na ně jetdy podívat ručně. Jeho obsah vypadá následovně:

<?xml version='1.0' encoding='UTF-8'?>
<feed xmlns:georss="http://www.georss.org/georss" xmlns="http://www.w3.org/2005/Atom" xml:lang="cs">
  <id>https://atom.cuzk.gov.cz/ORTOFOTO/datasetFeeds/CZ-00025712-CUZK_ORTOFOTO_WRTO24.2025.PRAH40.xml</id>
  <title>Ortofoto České republiky - mapový list: Praha 4-0</title>
  <updated>2026-03-11T08:58:38+02:00</updated>
  <author>
    <name>Český úřad zeměměřický a katastrální</name>
    <uri>https://geoportal.cuzk.gov.cz</uri>
    <email>cuzk.helpdesk@cuzk.gov.cz</email>
  </author>
  <rights>žádné podmínky neplatí</rights>
  <link href="https://atom.cuzk.gov.cz/ORTOFOTO/ORTOFOTO.xml" hreflang="cs" type="application/atom+xml" rel="up" />
  <link href="https://geoportal.cuzk.gov.cz" rel="describedby" hreflang="cs" type="text/html" />
  <entry>
    <id>https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip</id>
    <title>Ortofoto České republiky - mapový list: Praha 4-0</title>
    <updated>2026-03-11T08:58:38+02:00</updated>
    <link href="https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip" length="59466093" rel="alternate" type="image/jpeg" hreflang="cs" />
    <category label="S-JTSK" term="http://www.opengis.net/def/crs/EPSG/0/5514" />
    <georss:polygon>50.103 14.4891 50.103 14.5275 50.1238 14.5275 50.1238 14.4891 50.103 14.4891</georss:polygon>
  </entry>
</feed>

Úkol

Proč nelze tento odkaz čísti pomocí gdal vector info. Co v jeho hlavičce chybí?

Tip

Přestože se na xml nelze dotázati pomocí nástrojů GDAL, existuje možnost velice neintuitivného automatického zpracování pomocí xmllint: curl https://atom.cuzk.cz/ORTOFOTO/datasetFeeds/CZ-00025712-CUZK_ORTOFOTO_WRTO24.2023.TREN42.xml | xmllint --xpath "//*[name() = 'feed']/*[name() = 'entry']/*[name() = 'link']/@href" -

Zde nalezneme odkaz ke stažení zazipovaného souboru uvnitř tagu <id>. Dotažme se na tento produkt pomocí gdal raster info (rastr), abychom získali jeho metadata. Abychom jej nemuseli rozbalovati, je nutné užití virtuálních driverů vsicurl (soubor je online) a vsizip (dotazujeme se na jeden soubor uvnitř archivu .zip bez jeho rozbalování).

gdal raster info /vsizip//vsicurl/https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip/WRTO24.2025.PRAH40.jpg

Tip

Seznam podadresářů a souborů, a tedy i relativní cestu k nim získáme dotazem na samotný zip: gdalinfo /vsizip//vsicurl/https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip (bohužel, zde gdal raster info ve verzi GDAL 3.11 nefunguje).

Dotažme se tedy na statistiky daného produktu:

gdal raster info --stats /vsizip//vsicurl/https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip/WRTO24.2025.PRAH40.jpg

Zadáním do webového prohlížece si ortofoto stáhneme. 

gdal raster convert --format JP2OpenJPEG /vsizip//vsicurl/https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip/WRTO24.2025.PRAH40.jpg WRTO24.2025.PRAH40.jpg

Rozhraní Python

Načtěme si data virtuálně, přetypujme je do formátu Byte (hodnoty 0 - 255), přidejme jim souřadnicový systém a zapišme je jako GeoTIFF.

import numpy as np

from osgeo import gdal, osr

in_file = '/vsizip//vsicurl/https://openzu.cuzk.gov.cz/opendata/ORTOFOTO/WRTO24.2025.PRAH40.zip/WRTO24.2025.PRAH40.jpg'
out_file = '/tmp/test.tif'
nr_bands = 3

im = gdal.Open(in_file)

arr = im.ReadAsArray()
# porovnejme prumer s prumerem ziskanym pomoci gdalinfo - jsou totozne?
print(np.mean(arr[0]))

driver = gdal.GetDriverByName('GTiff')
srs = osr.SpatialReference()

out = driver.Create(out_file, arr.shape[2], arr.shape[1], nr_bands, gdal.GDT_Byte)

out.SetGeoTransform(im.GetGeoTransform())
srs.ImportFromEPSG(5514)
out.SetProjection(srs.ExportToWkt())

for i in range(nr_bands):
    band = im.GetRasterBand(i + 1).ReadAsArray()
    band_int_enhanced = band / band.max() * 255
    out.GetRasterBand(i + 1).WriteArray(band_int_enhanced)

# zavreme data
im = None
out = None

Další materiály