GDAL/OGR - Geodaten automatisiert verarbeiten

GDAL/OGR - Geodaten automatisiert verarbeiten

Dieser Blogeintrag gibt eine Einführung in GDAL/OGR und erklärt anhand von Beispielen wie man die zahlreichen Kommandozeilenwerkzeuge nutzen kann.

/img/gdal_ogr/gdal_icon.png

GDAL/OGR ist eine Programm-Bibliothek die verschiedene Vektor- und Rasterformate lesen und schreiben kann. Sie ist komplett Open Source und in C/C++ geschrieben. Allerdings kann man auch mit anderen Programmiersprachen wie Python oder R darauf zugreifen. Es gibt zahlreiche Programme für die Kommandozeile. Diese werden in diesem Blogpost genauer vorgestellt.

Die meisten haben GDAL/OGR schon oft genutzt ohne es zu merken. Viele OpenSource Programme wie QGIS, PostGIS oder MapServer und auch einige properitäre Programm nutzen GDAL/OGR im Hintergrund um Geodaten zu lesen und zu schreiben.

Man kann mit GDAL/OGR durchaus auch einige räumliche Analysen machen. Jedoch wird es hauptsächlich für die Konvertierung von Datenformaten genutzt. Intern wandelt GDAL/OGR den eingelesen Datensatz in ein eigenes Datenmodell um schreibt daraus wieder das gewünschte Format.

Es gibt schon einige Blogeinträge, Videos und Cheatsheets (Spickzettel) zu diesem Thema. Diese sind unten bei den Links aufgeführt. Im folgenden stelle ich alle Befehle vor die schon oft gebraucht habe und zeige ein paar Features die ich spannend finde. Zum einen als Erinnerung für mich selbst und in der Hoffnung das andere auch davon profitieren können. Auf der AGIT 2019 habe ich diesen Vortrag dazu gehalten:

Die Dokumentation von GDAL ist sehr detailliert. Für jeden Treiber werden alle möglichen Optionen erklärt und am Schluss mit Beispielen für die Kommandozeile ergänzt.

Virtuelles Raster

Rasterdaten können schnell sehr groß werden. Deshalb werden häufig große Gebiete in mehreren einzelnen Rasterdateien unterteilt. Mit dem Programm gdalbuildvrt kann man sich eine “virtuelle Rasterdatei” erstellen lassen, die alle in ihr enthaltenen Rasterdateien referenziert. Das ist einfach eine Textdatei und kann in weiterer Folge wie eine ganz normale Rasterdatei verarbeitet werden. Dieser Befehl erstellt aus einem Ordner von GeoTIFF-Dateien ein virtuelles Raster:

gdalbuildvrt \
  output.vrt \
  my_folder/*.tif

Raster Pyramiden bauen

Das Darstellen von großen Rasterdateien dauert oft lange, weil im Prinzip jedes einzelne Pixel ausgelesen und gerendert werden muss. Um diesen Prozess zu beschleunigen kann mit gdaladdo (“GDAL add overview” ) aus den Pixeln Übersichten für niedrigere Zoomstufen erstellt werden (auch genannt Pyramiden). Diese können von den darstellenden Programmen direkt genutzt werden und das Anzeigen geht wesentlich schneller. Die ursprüngliche Datei wird dabei verändert bzw. sie wird größer, weil die zusätzlichen Übersichten in der Datei gespeichert werden.

gdaladdo \
  -r average \
  input.tif

bzw. für ältere GDAL Versionen

gdaladdo \
  -r average \
  input.tif \
  2 4 8 16
/img/gdal_ogr/build_pyramids.png

Raster umprojezieren

Mit gdalwarp kann man Raster auf verschiedenste Arten umprojezieren. In der Dokumenation sind sehr ausführliche Beispiele beschrieben. Ich habe bis jetzt aber nur diesen Befehl gebraucht:

gdalwarp \
  -t_srs "EPSG:4326" \
  input.tif \
  output.tif

MrSID Dateien

Oftmals bekommt man Rasterdateien im MrSID Format. Auf Linux ist es scheinbar möglich diese zu lesen. Einfacher geht es jedoch wenn man mit Windows auf MrSID zugreift. Dann kann man die Datei zum Beispiel in ein GeoTIFF umwandeln und auf Linux weiterarbeiten.

Kacheln

Das Programm gdal2tiles.py konvertiert eine Rasterdatei in einen Ordner voller Kacheln nach dem TMS (Tile Map Service) Standard. Dieser ist sehr ähnlich zu dem XYZ Standard. Bei TMS ist die Kachelnummer (0,0) unten links, beim XYZ hingegen oben links.

gdal2tiles.py \
  --zoom=2-5 \
  input.tif \
  output_folder

Die Ordnerstruktur sieht dann folgendermaßen aus. Der Übersichtlichkeit wegen ist nur die Zoomstufe 2 (und nicht 3,4 und 5) enthalten:

output_folder
├── 2
│   ├── 0
│   │   ├── 0.png
│   │   ├── 1.png
│   │   ├── 2.png
│   │   └── 3.png
│   ├── 1
│   │   ├── 0.png
│   │   ├── 1.png
│   │   ├── 2.png
│   │   └── 3.png
│   ├── 2
│   │   ├── 0.png
│   │   ├── 1.png
│   │   ├── 2.png
│   │   └── 3.png
│   └── 3
│       ├── 0.png
│       ├── 1.png
│       ├── 2.png
│       └── 3.png
├── googlemaps.html
├── leaflet.html
├── openlayers.html
└── tilemapresource.xml

5 directories, 20 files

Der Zielordner enthält nicht nur Kacheln sondern auch HTML-basierte Clients mit Leaflet und OpenLayers. Die Kacheln können dann in einen Webserver (z.B. Apache HTTP Server) gelegt werden und auch mit einem beliebigen Client aufgerufen werden. Ein typische URL sieht dann so aus: http://localhost/my_tiled_raster/{z}/{x}/{y}.png. Bei manchen Clients wie QGIS oder OpenLayers muss allerdings die Y-Ziffer invertiert werden, also: http://localhost/my_tiled_raster/{z}/{x}/{-y}.png (mehr Infos gibt es dazu in dieser QGIS issue ).

/img/gdal_ogr/gdal2tiles_leaflet_viewer.png

Auf Leaflet basierender Viewer für Kacheln

ogrinfo

ogrinfo gibt Informationen über eine räumliche Datenquelle zurück. In diesem Beispiel wird das GeoPackage von NaturalEarth abgefragt:

ogrinfo \
  natural_earth_vector.gpkg   

ergibt (Ausschnitt):

INFO: Open of `natural_earth_vector.gpkg'
      using driver `GPKG' successful.
1: ne_10m_admin_0_antarctic_claim_limit_lines (Line String)
2: ne_10m_admin_0_antarctic_claims (Polygon)
3: ne_10m_admin_0_boundary_lines_disputed_areas (Line String)
4: ne_10m_admin_0_boundary_lines_land (Line String)
5: ne_10m_admin_0_boundary_lines_map_units (Line String)
...

So bekommt man Informationen über einen Layer. Die Option -so steht für “summary only” und bewirkt dass nicht jedes einzelne Feature angezeigt wird:

ogrinfo \
  -so \
  natural_earth_vector.gpkg \
  ne_10m_admin_0_antarctic_claim_limit_lines

Dies gibt wichtige Informationen (Anzahl der Features, Ausdehnung, Projektion, …) über den Layer zurück.

INFO: Open of `natural_earth_vector.gpkg'
      using driver `GPKG' successful.

Layer name: ne_10m_admin_0_antarctic_claim_limit_lines
Geometry: Line String
Feature Count: 23
Extent: (-150.000000, -90.000000) - (160.100000, -60.000000)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
FID Column = fid
Geometry Column = geom
type: String (15.0)
scalerank: Integer (0.0)
featurecla: String (50.0)

ogr2ogr

Das Werkzeug ogr2ogr kann Vektordaten einlesen, verarbeiten und in andere Formate umwandeln. Dies wird im folgenden hauptsächlich am Beispiel vom GeoPackage erklärt. Dieses Format ist ein vielversprechender Kandidat um das Shapefile als Austauschformat für Geodaten abzulösen.

ogr2ogr Befehle können schnell sehr verwirrend aussehen, sind aber im Prinzip nach dem selben Muster aufgebaut.

Eine einfache Kopie eines Datensatzes sieht folgendermaßen aus:

ogr2ogr \
  -f GPKG output.gpkg \
  input.gpkg

Umwandlung der Projektion in EPSG:3857:

ogr2ogr \
  -t_srs EPSG:3857 \
  -f GPKG output.gpkg \
  input.gpkg

Dabei muss die Projektion des Input-Datensatzes allerdings bekannt sein. Ansonsten kann diese mit -s_srs angegeben werden:

ogr2ogr \
  -s_srs EPSG:4326 \
  -t_srs EPSG:3857 \
  -f GPKG output.gpkg \
  input.gpkg

GPX

GPS- bzw. GNSS-Geräte stellen ihre Daten zumeist als GPX-Datei bereit. Diese kann zwar in QGIS geöffnet werden, allerdings ist das immer mit etwas Aufwand verbunden. Die Umwandlung in ein GeoPackage funktioniert wie folgt:

ogr2ogr \
  -f GPKG output.gpkg \
  input.gpx

Im GeoPackage sind dann folgende Layer zu finden:

ogrinfo output.gpkg

ergibt:

INFO: Open of `output.gpkg'
      using driver `GPKG' successful.
1: tracks (Multi Line String)
2: track_points (Point)
3: waypoints (Point)
4: routes (Line String)
5: route_points (Point)

Wenn man diese Operation häufiger braucht, kann man sich daraus eine Funktion erstellen:

function gpx2gpkg {
  input_gpx="$1"
  filename_without_suffix=${input_gpx%.gpx}
  ogr2ogr -f GPKG "$filename_without_suffix".gpkg $input_gpx
}

CSV

Sehr häufig bekommt man Daten im CSV-Format. Diese können leicht in jedes beliebige Format konvertiert werden. Mittels der Option -oo (“output optiones”) können vielen Einstellungen vorgenommen werden. In unteren Beispiel werden die Spaltennamen der Koordinaten angegeben. Es ist wichtig das entsprechende Koordinatensystem des Datensatzes mittels -a_srs zuzuweisen.

ogr2ogr \
  -f GPKG output.gpkg \
  input.csv \
  -oo X_POSSIBLE_NAMES=laengengrad \
  -oo Y_POSSIBLE_NAMES=breitengrad \
  -a_srs 'EPSG:4326'

Der neu erstelle Layer wird durch eine Geometriespalte erweitert. Die ursprünglichen Spalten mit den Koordinaten bleiben erhalten. Wenn man dies allerdings nicht möchte kann man es mit der Option -oo KEEP_GEOM_COLUMNS=NO verhindern.

PostGIS

PostGIS ist die räumliche Erweiterung der PostgreSQL Datenbank und eignet sich hervorragend um Geodaten zu speichern (Link zum Treiber).

Eine Tabelle exportieren:

ogr2ogr \
  -f GPKG output.gpkg \
  PG:dbname="my_database" "my_table"

Mehrere Tabellen exportieren:

ogr2ogr \
  -f GPKG output.gpkg \
  PG:'dbname=my_database tables=table_1,table_3'

Die ganze Datenbank in ein GeoPackage exportieren:

ogr2ogr \
  -f GPKG output.gpkg \
  PG:dbname=my_database

Die ganze Datenbank in einen Ordner voller Shapefiles exportieren:

ogr2ogr \
  -f "ESRI Shapefile" output_folder \
  PG:dbname=my_database

Daten in PostGIS laden:

ogr2ogr \
  -f "PostgreSQL" PG:dbname="my_database" \
  input.gpkg \
  -nln "name_of_new_table"

Beim Import werden die Namen der Tabellen und Spalten zu Kleinbuchstaben konvertiert. Man kann allerdings mit der Option -lco LAUNDER=NO Großbuchstaben erzwingen. Allerdings empfiehlt es sich wenn möglich alles in Kleinbuchstaben zu schreiben (außer man hat gute Gründe die dagegen sprechen).

ogr2ogr \
  -f "PostgreSQL" PG:dbname="my_database" \
  input.gpkg \
  -lco LAUNDER=NO

Beim einlesen von Multigeometrien kann es zu folgendem Fehler kommen: ERROR: Geometry type (MultiPolygon) does not match column type (Polygon) Das kann mit der Option -nlt PROMOTE_TO_MULTI gelöst werden:

ogr2ogr \
  -f "PostgreSQL" PG:dbname="my_database" \
  input.gpkg \
  -nlt PROMOTE_TO_MULTI

Wenn man einen ganzen Ordner voller GeoPackages (oder Shapefiles, GeoJSONs, … ) hat. Kann man diese mit einer For-Schleife gesammelt in PostGIS laden.

for file in *.gpkg;
  do ogr2ogr \
  -f "PostgreSQL" PG:dbname="my_database" \
  $file;
done

Bestehende Tabellen kann man mit der Option -lco OVERWRITE=YES überschreiben:

ogr2ogr \
  -f "PostgreSQL" PG:dbname="my_database" \
  -lco OVERWRITE=YES \
  input.gpkg

Man kann sogar einen Postgres Dump erzeugen ohne überhaupt eine PostgreSQL Datenbank zu haben.

Abfragen und Filtern von Daten

ogr2ogr kann auch Geodaten filtern und bearbeiten. Eine häufige Anwendung ist die Eingrenzung eines Datensatzes auf ein bestimmtes Gebiet. Diese kann mit der Option -spat erzeugt werden:

ogr2ogr \
  -spat -13.931 34.886 46.23 74.12 \
  -f GPKG output.gpkg \
  natural_earth_vector.gpkg

Wenn man nur bestimmte Features eines Layers haben möchte, kann man diese mit -where definieren. Dieses Beispiel gibt nur Länder zurück deren Einwohnerzahl kleiner als 1 Million ist:

ogr2ogr \
  -where "\"POP_EST\" < 1000000" \
  -f GPKG output.gpkg \
  natural_earth_vector.gpkg \
  ne_10m_admin_0_countries

Die Bedingungen die nach dem -where steht muss in " stehen. Wenn darin allerdings weitere " vorkommen müssen diese mit einen \ gültig gemacht werden (“escaped”). So sieht das Ergebnis auf der Karte aus:

/img/gdal_ogr/countries_pop_less_one_million.png

Die beiden Optionen -spat und -where können natürlich auch kombiniert werden. Dies gibt dann die Länder in Europa unter 1 Million Einwohner zurück:

ogr2ogr \
  -where "\"POP_EST\" < 1000000" \
  -spat -13.931 34.886 46.23 74.12 \
  -f GPKG output.gpkg \
  natural_earth_vector.gpkg \
  ne_10m_admin_0_countries

Mit der Option -sql kann man Abfragen mittels SQL durchführen. Obwohl SQL eigentlich für Datenbanken gedacht ist, kann man mit ogr2ogr es auch für alle beliebigen Datensätze anwenden.

Das folgende Beispiel nimmt als Input diese CSV-Datei:

ID,Name,Breitengrad,Laengengrad
1,Salzburg, 47.8, 13.0
2,Heidelberg, 49.4, 8.7
3,Trondheim, 63.4, 10.4

Mit diesem Befehl werden nur Städte ausgegeben deren Breitengrad kleiner als 50° ist:

ogr2ogr \
-f "GPKG" output.gpkg \
staedte.csv \
-sql "SELECT * FROM staedte WHERE \"Breitengrad\" < 50" \
-oo X_POSSIBLE_NAMES=Laengengrad \
-oo Y_POSSIBLE_NAMES=Breitengrad \
-a_srs 'EPSG:4326'

Es gibt nicht das “eine” SQL, sondern es gibt verschiedenen Dialekte. ogr2ogr kann immer den Dialekt der entsprechenden Input-Datenqulle nehmen. Beispielsweise Postgres-SQL wenn es sich um eine Postgres/PostGIS Datenbank handelt. Alternativ auch ein eigenes OGR SQL oder SQLite SQL. Dies kann mit der Option -dialect eingestellt werden. Eine Besonderheit von SQLite SQL ist, dass man damit sogar Daten geocoden kann. Also für eine Adresse eine Koordinate oder für eine Koordinate eine Adresse bekommt.

Mit SQL kann man aus Koordinaten Geometrien erzeugen. Im folgenden Beispiel wird eine nicht-räumliche SQLite-Tabelle in eine räumliches GeoPackage umgewandelt. Dabei wird mit der Funktion MakePoint() aus den Koordinaten eine gültige Geometrie erzeugt.

ogr2ogr \
  -f "GPKG" output.gpkg \
  input.sqlite \
  -sql \
  "SELECT
     *,
     MakePoint(longitude, latitude, 4326) AS geometry
   FROM
     my_table" \
  -nln "location" \
  -s_srs "EPSG:4326"

Weitere Funktionen

Mit -append kann man verschiedenen Layer zusammenlegen. Diese müssen allerdings eine gleiche Struktur aufweisen. Das folgende Beispiel nimmt alle GeoPackages in einem Ordner und vereinigt diese:

for geopackage in *gpkg; do
  ogr2ogr \
  -f "GPKG" combined.gpkg \
  -append \
  $geopackage
done

Mit der Option -progress wird angezeigt wieviel Prozent der Umwandlung schon geschehen ist. Das mach vor allem Sinn für Operationen die sehr lange dauern.

Die Option –debug ON zeigt bei der Ausführung des Befehl sehr viele weitere Informationen an. Das ist besonders hilfreich, wenn der Befehl nicht funktioniert und man den Fehler finden möchte.

Python

Es gibt ein gdal Package für Python. Die Syntax der Funktionen ist allerdings sehr nah an der C++ API und deswegen eher anspruchsvoll. Für Rasterdaten gibt es den Wrapper rasterio. Für Vektordaten gibt es den Wrapper fiona. In dessen Dokumention wird beschrieben für welche Fälle fiona geeignet ist und in welchen Fällen das ogr2ogr direkt verwendet werden sollte.

Man kann auch innerhalb von Python auf ogr2ogr in der Kommandozeile zugreifen. Um das zu vereinfachen steht das Script ogr2ogr.py (Link) zu Verfügung. Dies wird wie folgt aufgerufen:

import ogr2ogr

ogr2ogr.main([
  'ogr2ogr',
  '-f', 'GPKG', 'output.gpkg' ,
  'input.gpkg'
])

Danke an die Entwickler

Wie schon erwähnt ist GDAL/OGR frei und OpenSource. An dieser Stelle möchte ich Frank Warmerdam danken, der das Projekt gestartet hat außerdem Even Rouault der mittlerweile der Hauptentwickler geworden ist. Und natürlich den vielen anderen Beitragenden die das Projekt weiterentwickeln und am laufen halten.

Deutsch

  • Im Vortrag “GeoPackages der freien Hamburger Geodaten” macht Johannes Kröger intensiven Gebrauch von ogr2ogr
  • Im Wiki der Hochschule Rapperswil gibt es einige Beispiele zu ogr2ogr
  • Der Vortrag von Jakob Miksch auf der AGIT 2019 gibt einen Überblick über die Kommandozeilenwerkzeuge von GDAL/OGR

English

CheatSheets