Knihovny pro práci s geoprostorovými daty (PROJ, GDAL, PDAL)

PROJ

Knihovna PROJ (dříve označovaná jako PROJ.4) zajišťuje převod souřadnic mezi jednotlivými souřadnicovými systémy. Jde o jednu ze základních knihoven využívaných v mnoha open source GIS projektech jako GRASS GIS, QGIS, PostGIS a dalších.

Knihovna je napsána v C++, existují ale rozhraní do dalších programovacích jazyků jako je Javascript, PHP nebo Python.

Konzolové nástroje

Knihovna PROJ nabízí řadu konzolových nástrojů, viz dokumentace.

cct - převod 4D souřadnic

echo 15 50 0 0 | cct +proj=utm +zone=32 +datum=WGS84

projinfo - dotazování se na souřadnicové systémy
```
projinfo EPSG:32632 -o proj
```
cs2cs - převod souřadnic mezi jednotlivými souřadnicovými systémy
```
echo 929865.50 5555901.55 | cs2cs EPSG:32632 EPSG:5514
```

Úkol

Porovnejte výsledek předchozího příkazu s EPSG:5514 code 5239:

# projinfo EPSG:5514 -o proj
    echo 929865.50 5555901.55 | cs2cs EPSG:32632 +to +proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 +alpha=30.2881397527778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=572.213,85.334,461.94,-4.9732,-1.529,-5.2484,3.5378 +units=m +no_defs +type=crs

geod - první geodetická úloha

# Kladno, 50d, 1km
# -> lat,log,azimut
echo 50.1473350N 14.1028503E 50 1000 | geod +ellps=bessel +units=m

invgeod - první geodetická úloha (inverzní)

# Kladno, Mladá Boleslav
# -> azimut, azimut zpět, vzdálenost
echo 50.1473350N 14.1028503E 50.4113514N 14.9031836E | invgeod +ellps=bessel +units=m

gie - testovací nástroje knihovny

proj - převod souřadnic mezi jednotlivými souřadnicovými systémy

echo 15 50 | proj +proj=utm +zone=32 +datum=WGS84

echo 15 50 | proj EPSG:32632

echo 15 50 | proj EPSG:5514

invproj - opačný převod souřadnic

echo 929865.50 5555901.55 | invproj EPSG:32632

projsync - umožňuje stahovat a aktualizovat podpůrná data (výpočetní gridy)

Python API

Ukázka transformace zeměpisných souřadnic (EPSG:4326) do S-JTSK (EPSG:5514):

from pyproj import CRS
crs_from = CRS.from_epsg(4326)
crs_to = CRS.from_epsg(5514)

from pyproj import Transformer
transformer = Transformer.from_crs(crs_from, crs_to)
transformer.transform(50, 15)

Úkol

Porovnej výsledek s příkazem cs2cs.

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é Python API. 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

Rastrová data

gdalinfo: Vypíše informace o rastrovém datasetu

gdalinfo /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF

gdalinfo -stats -json /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF

gdalsrsinfo: Vypíše informace o souřadnicových systémech

gdalsrsinfo /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF

gdal_translate: Převede rastrová data mezi zvolenými datovými formáty

# gdalinfo --formats
gdal_translate -of GPKG /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF eu_dem_v11_E40N10.gpkg

gdal_translate -of GTiff -projwin 4276219 1748621 5041965 1493878 /mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF eu_dem_v11_E40N10_REG.TIF

Poznámka

Volba -projwin předpokládá souřadnice ve formě UL LR. Pozor, v QGISu se ale rozsah uzemí zobrazuje ve formě UR LL.

gdalwarp: Transformace rastrových dat mezi souřadnicovými systémy

gdalwarp -t_srs EPSG:4326 /mnt/repository/155FGIS/03/eu_dem_v11_E40N10_REG.TIF eu_dem_v11_E40N10_WGS_84.TIF

Tip

Příkaz gdalwarp má parametr -te (souřadnice v cílovém souřadnicovém systému) odpovídající parametru gdal_translate -projwin. Prostorový rozsah lze určit i na základě vstupní vektorové vrstvy (-cutline).

gdallocationinfo: Nástroj pro dotazování

gdallocationinfo -geoloc eu_dem_v11_E40N10_REG.TIF 4659075 1621250

gdaladdo: Vytváří nebo obnovuje přehledové obrazy

# ls -l eu_dem_v11_E40N10.TIF
gdaladdo eu_dem_v11_E40N10.TIF 2 4 8 16

Úkol

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

gdaldem: Nástroj pro analýzu a vizualizaci DEM

gdaldem hillshade eu_dem_v11_E40N10_REG.TIF eu_dem_v11_E40N10_REG_HS.TIF

gdal_contour: Vytvoří vrstevnice z digitálního výškového modelu

gdal_contour -a elev eu_dem_v11_E40N10_REG.TIF vrstevnice.shp -i 10.0

gdal_rasterize: Rasterizace vektorových dat

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

gdalmanage: Nástroj pro správu rastrových souborů
```
gdalmanage identify /mnt/repository/155FGIS/03/*
```

Další příkazy:

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

ogrinfo vrstevnice.shp

ogrinfo vrstevnice.shp -so -al

# datatový zdroj s více vrstvami
ogrinfo /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg

ogrinfo /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -so

# přistup k prvkům
ogrinfo /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -fid 1

ogrinfo /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -fid 1 -geom=SUMMARY

# počet polygonů spadající do skupiny jehličnatý les
ogrinfo U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -where "Code_18 = 312" -so | grep Count

Poznámka

Nomenklatura pro atribut Code_18

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

# atributový filtr
ogr2ogr -f "ESRI Shapefile" clc2018_312.shp /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -where "Code_18 = 312"

# kombinace atributového a prostorového filtru
# gdalinfo eu_dem_v11_E40N10_REG.TIF
ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:4326 -spat 4276200 1493875 5041950 1748625 clc2018_312_reg_wgs84.shp /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -where "Code_18 = 312"

Další příkazy:

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

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

ds = gdal.Open('/mnt/repository/155FGIS/03/eu_dem_v11_E40N10.TIF')
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]))
ds = None

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

from osgeo import ogr

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

Příklad konverze dat:

# ogr2ogr -f "ESRI Shapefile" clc2018_312.shp /mnt/repository/155FGIS/03/U2018_CLC2018_V2020_20u1.gpkg U2018_CLC2018_V2020_20u1 -where "Code_18 = 312"

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

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:

ogrinfo /vsizip/data50.zip

ogrinfo /vsizip/data50.zip Zamek -so

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

ogr2ogr -f GPKG data50.gpkg /vsizip/data50.zip

Úkol

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

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

ogr2ogr -f GPKG -t_srs EPSG:4326 zamek_hrad.gpkg /vsizip/data50.zip Zamek Hrad

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:

gdalbuildvrt osm_lc.vrt /mnt/repository/155FGIS/09/raster_tiles/*.tif

gdalinfo -noct osm_lc.vrt

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

gdal_translate -of COG osm_lc.vrt osm_lc.tif

Poznámka

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

gdalinfo -noct /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 ogrinfo:

ogrinfo -ro 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: ogrinfo --config GDAL_HTTP_CONNECTTIMEOUT=666666 -ro https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO

ogrinfo nám bez specifikování dotazované vrstvy vrací seznam vrstev. Podívejme se tedy na data obsažená v jediné vrácené vrstvě (georss):

ogrinfo -ro https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO georss

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:

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

Tip

Pokud chcete získat pouze odkaz link2_href, můžete jej z výstupu ogrinfo vyfiltrovat pomocí ogrinfo -ro -where "title LIKE '%Praha 4-0'" https://atom.cuzk.cz/getservicefeed.ashx?service=ORTOFOTO georss | 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ě:

<feed xml:lang="cs">
<id>
https://atom.cuzk.cz/ORTOFOTO/datasetFeeds/CZ-00025712-CUZK_ORTOFOTO_WRTO24.2023.TREN42.xml
</id>
<title>Ortofoto České republiky - mapový list: Třeboň 4-2</title>
<updated>2024-02-14T08:52:51+02:00</updated>
<author>
<name>Český úřad zeměměřický a katastrální</name>
<uri>https://geoportal.cuzk.cz</uri>
<email>cuzk.helpdesk@cuzk.cz</email>
</author>
<rights>žádné podmínky neplatí</rights>
<link href="https://atom.cuzk.cz/ORTOFOTO/ORTOFOTO.xml" hreflang="cs" type="application/atom+xml" rel="up"/>
<link href="https://geoportal.cuzk.cz" rel="describedby" hreflang="cs" type="text/html"/>
<entry>
<id>
https://openzu.cuzk.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip
</id>
<title>Ortofoto České republiky - mapový list: Třeboň 4-2</title>
<updated>2024-02-14T08:52:51+02:00</updated>
<link href="https://openzu.cuzk.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip" length="59045908" rel="alternate" type="image/jpeg" hreflang="cs"/>
<category label="S-JTSK" term="http://www.opengis.net/def/crs/EPSG/0/5514"/>
<georss:polygon>
48.998 14.7186 48.998 14.7561 49.0188 14.7561 49.0188 14.7186 48.998 14.7186
</georss:polygon>
</entry>
</feed>

Úkol

Proč nelze tento odkaz čísti pomocí ogrinfo. 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í gdalinfo (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í).

gdalinfo /vsizip//vsicurl/https://openzu.cuzk.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip/WRTO24.2023.tren42.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.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip

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

gdalinfo -stats /vsizip//vsicurl/https://openzu.cuzk.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip/WRTO24.2023.tren42.jpg

Zadáním do webového prohlížece si ortofoto stáhneme. Chceme-li zůstati věrni terminálu, můžeme jej stáhnouti též pomoc příkazů jako wget nebo curl:

wget https://openzu.cuzk.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip

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.cz/opendata/ORTOFOTO/WRTO24.2023.TREN42.zip/WRTO24.2023.tren42.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

PDAL

PDAL je knihovna pro práci s mračny bodů. Podobně jako knihovna GDAL je napsaná v C++. Nicméně umožňuje dopisovat výpočetní filtry i v jiných jazycích jako je např. Python. Kromě toho obsahuje i užitečné konzolové nástroje.

Konzolové nástroje

pdal --help

pdal info /mnt/repository/155FGIS/06/dmr5g/KRAV59.laz

# pdal --drivers
# pdal --options writers.las
pdal translate /mnt/repository/155FGIS/06/dmr5g/KRAV59.laz KRAV59.las

pdal translate /mnt/repository/155FGIS/06/dmr5g/KRAV59.laz KRAV59_wgs84.las -f filters.reprojection \
  --filters.reprojection.out_srs="EPSG:4326" --filters.reprojection.in_srs="EPSG:5514"

pdal merge /mnt/repository/155FGIS/06/dmr5g/*.laz dmr5g.laz

Python API

json = """
[
    "/mnt/repository/155FGIS/06/dmr5g/KRAV59.laz",
    {
        "type": "filters.sort",
        "dimension": "Z"
    }
]
"""

import pdal
pipeline = pdal.Pipeline(json)
count = pipeline.execute()
arrays = pipeline.arrays
print(arrays)

Implementace filtru v Pythonu

import numpy as np

def multiply_z(ins, outs):
    Z = ins['Z']
    Z *= 10.0
    outs['Z'] = Z

    return True

[
    "/mnt/repository/155FGIS/06/dmr5g/KRAV59.laz",
    {
        "type":"filters.smrf"
    },
    {
        "type":"filters.python",
        "script":"multiply_z.py",
        "function":"multiply_z",
        "module":"anything"
    },
    {
        "type":"writers.las",
        "filename":"KRAV59_multiplied.las"
    }
]

pdal pipeline ./pdal-pipeline.json

Další materiály

GISMentors: GeoPython pro začátečníky