GDAL/OGR - Geodaten automatisiert verarbeiten

Ursprünglich wurde der Blogbeitrag auf https://jakobmiksch.eu/de/post/gdal_ogr/ veröffentlicht.

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

GDAL/OGR ist eine Programm-Bibliothek die verschiedene Vektor- und Rasterformate lesen und schreiben kann. Sie ist Open Source und in hauptsächlich 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

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 ).

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)

Die Funktion -al führt ogrinfo für jeden Layer der Datasource aus. Dies kann manchmal praktisch sein, kann aber auch schnell recht unübersichtlich werden. In Kombination mit -so wird die Übersicht über jeden Layer einer Datasource angezeigt.

ogrinfo -so -al natural_earth_vector.gpkg

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

Shapefile

Shapefiles müssen immer aus mindestens drei Dateien bestehen. Nämlich .shp, .dbf und .shx. Es genügt aber bei OGR nur die .shp Datei anzugeben, wobei die anderen Dateien natürlich trotzdem vorhanden sein müssen. Seit GDAL 3.1 kann man auch gezippte Shapfiles gelesen werden.

Wenn man ein Shapefile mittels ogrinfo abfrägt bekommt man nur den Namen des einen Layers zurück. Man muss also den Namen des einen Layers explizit angeben um die gewünschten Infos zu bekommen. Der Grund dafür ist, dass ein Shapefile Datasource und Layer in einem ist.

# gibt nur den Namen des einen Layers zurück
ogrinfo countries.shp

# gibt die gewünschten Information zurück
ogrinfo countries.shp countries

Oftmals liegen Shapefiles gesammelt in einem Ordner. Praktischerweise kann OGR das als Datasource erkennen. Siehe folgendes Bespiel:

my_shapefiles
├── countries.dbf
├── countries.prj
├── countries.shp
├── countries.shx
├── lakes.dbf
├── lakes.prj
├── lakes.shp
└── lakes.shx

Die Abfrage ogrinfo my_shapefiles ergibt:

INFO: Open of 'my_shapefiles'
      using driver 'ESRI Shapefile' successful.
1: countries (Polygon)
2: lakes (Polygon)

Die einzelnen Shapefiles werden also als Layer angezeigt. Damit kann ich nun Inforationen mittels ogrinfo abfragen:

# Informationen über einen bestimmten Layer
ogrinfo my_shapefiles lakes

# Zusammenfassung von allen Shapefiles des Ordners
ogrinfo -al -so my_shapefiles

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

Hier eine Übersicht über die meisten Parameter einer Postgres-Verbindung:

ogr2ogr \
-f GPKG output.gpkg \
PG:"
  dbname=my_database
  host=localhost
  port=5432
  user=user_name
  password=my_password
  schemas=schema1
  tables=table1,table2"

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:

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

Oftmals werden die SQL-Abfragen sehr lange. Man kann diese praktischerweise auch in eine separate Textdatei schreiben und dann im ogr2ogr Befehl referenzieren. Der vorherige Befehl kann beispielsweise aufgeteilt werden. Einmal in die separate Textdatei names query.sql:

SELECT *
  FROM staedte
 WHERE "Breitengrad" < 50

Und in den ogr2ogr Befehl, wobei der Pfad zu der Textdatei innerhalb von -sql mit einem @ referenziert wird:

ogr2ogr \
-f GPKG output.gpkg \
staedte.csv \
-sql @path/to/query.sql \
-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.

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

ogr2ogr: 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.

Es können auch Dateien im Internet mit GDAL/OGR verarbeitet werden. Dieser Befehl lädt bespielsweise eine GeoJSON herunter und verwandelt es in ein GeoPackage:

ogr2ogr \
-f GPKG chapters.gpkg \
https://raw.githubusercontent.com/maptime/maptime.github.io/master/_data/chapters.json

Virtuelles Vektor Format

Ähnlich wie bei Raster gibt es auch für Vektoren einen Treiber mit dem man virtuelle Vektordateien erstellen kann (siehe VRT). Im Prinzip werden mithilfe einer XML-Datei eine oder mehrere Vektordatenquellen referenziert. Die VRT-Datei ist aber nicht nur eine einfache Verknüpfung zu der referenzierten Datei, es können auch zahlreiche Umwandlungen druchgeführt werden. Beispielsweise können Felder umbenannt oder weggelassen werden oder das Koordinatensystem umprojeziert werden.

Folgendes VRT erstellt eine einfache Referenz auf ein im Internet verfügbares GeoJSON:

<OGRVRTDataSource>
  <OGRVRTLayer name="my_new_layer_name">
    <SrcDataSource relativeToVRT="0" shared="1">https://raw.githubusercontent.com/maptime/maptime.github.io/master/_data/chapters.json</SrcDataSource>
    <SrcLayer>chapters</SrcLayer>
  </OGRVRTLayer>
</OGRVRTDataSource>

Dieses VRT kann ganz normal mit ogrinfo oder ogr2ogr weiterverarbeitet werden oder auch mit QGIS geöffnet werden. Diese XML-Dateien werden üblicherweise per Hand erstellt. Allerdings gibt es auch ein Kommandozeilen-Werkzeug names ogr2vrt.py welches ein Grundgerüst der XML-Struktur erzeugt, welches wiederum manuell bearbeitet werden kann. Das Programm ist nich in der normalen GDAL-Installation vorhanden, sondern muss selber vom Repository heruntergeladen werden. Die Anwendung funkioniert wie folgt:

ogr2vrt.py \
  https://raw.githubusercontent.com/maptime/maptime.github.io/master/_data/chapters.json \
  chapters.vrt

erzeugt diese VRT-Datei:

<OGRVRTDataSource>
  <OGRVRTLayer name="chapters">
    <SrcDataSource relativeToVRT="0" shared="1">https://raw.githubusercontent.com/maptime/maptime.github.io/master/_data/chapters.json</SrcDataSource>
    <SrcLayer>chapters</SrcLayer>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>GEOGCS[&quot;WGS 84&quot;,DATUM[&quot;WGS_1984&quot;,SPHEROID[&quot;WGS 84&quot;,6378137,298.257223563,AUTHORITY[&quot;EPSG&quot;,&quot;7030&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6326&quot;]],PRIMEM[&quot;Greenwich&quot;,0,AUTHORITY[&quot;EPSG&quot;,&quot;8901&quot;]],UNIT[&quot;degree&quot;,0.0174532925199433,AUTHORITY[&quot;EPSG&quot;,&quot;9122&quot;]],AXIS[&quot;Latitude&quot;,NORTH],AXIS[&quot;Longitude&quot;,EAST],AUTHORITY[&quot;EPSG&quot;,&quot;4326&quot;]]</LayerSRS>
    <Field name="location" type="String" src="location"/>
    <Field name="title" type="String" src="title"/>
    <Field name="twitter" type="String" src="twitter"/>
    <Field name="website" type="String" src="website"/>
    <Field name="meetup" type="String" src="meetup"/>
    <Field name="comingSoon" type="String" src="comingSoon"/>
    <Field name="organizers" type="String" src="organizers"/>
    <Field name="moreInfo" type="String" src="moreInfo"/>
    <Field name="github" type="String" src="github"/>
    <Field name="facebook" type="String" src="facebook"/>
  </OGRVRTLayer>
</OGRVRTDataSource>

Diese Datei kopiert eigentlich nur die Eigenschaften die der Layer eh schon hat und eignet sich primär als Vorlage für eigene Veränderungen. Beispielsweise kann man einer Spalte einen neuen Namen zuweisen ohne die ursprüngliche Datenquelle zu ändern: <Field name="My new Name" type="String" src="title"/>

VRTs haben eine Vielzahl von nützlichen Anwendungsgebieten. So kann man eine Datenbank, oder WFS-Verbindung in eine VRT auslagern, und so tun als ob es eine einfache Datei wäre. Auch die Koordinaten in Excel-Dateien können referenziert werden, somit ist es möglich diese beliebig als räumliche Daten weiter zu verarbeiten oder beispielsweise in QGIS darzustellen.

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'
])

Docker

Neben der normalen Installation können wir die GDAL/OGR-Tools auch per Docker verwenden. Dies hat den Vorteil, dass wir die neueste Version oder eine beliebige ältere Version verwenden können. Es sind fünf Docker-Images verfügbar:

• alpine-normal
• alpine-small
• alpine-ultrasmall
• ubuntu-full
• ubuntu-small

Die detaillierte Beschreibung ist im offiziellen README zu finden. Diese enthalten eine unterschiedliche Anzahl von Treibern und haben eine unterschiedliche Downloadgröße. Daher müssen wir zuerst überlegen, welche Formate wir öffnen wollen, und dann entscheiden, welches Docker-Image wir wählen.

Die folgenden Beispiele verwenden alpine-small, da wir nur das GeoPackage und den Postgres-Treiber benötigen. Stelle sicher, dass docker auf dem Computer installiert ist. Navigiere zu einem Ordner, der Geodaten enthält. In den kommenden Beispielen nutzen wir ein GeoPackage namens sample.gpkg.

Der folgende Befehl öffnet eine sh-Shell innerhalb des Docker-Containers. Wir mounten das aktuelle Host-Verzeichnis in das Verzeichnis /data innerhalb des Docker-Containers.

docker run \
  -t \
  -i \
  -v $(pwd):/data \
  osgeo/gdal:alpine-small-latest \
  sh

Alle GDAL/OGR-Befehle sind in dieser Shell verfügbar. Damit können wir Informationen über das GeoPackage erhalten oder es umwandeln:

# die Version von GDAL/OGR
ogrinfo --version

# grundlegende Infos über das GeoPackage
ogrinfo /data/sample.gpkg

# wir wandeln den 'layer_1' aus unserem GeoPackage in ein GeoJSON um
ogr2ogr \
  -f GeoJSON \
  /data/out.json \
  /data/sample.gpkg \
  layer_1

Der vorangegange Befehl erstellt eine neue Datei, die in dem Ordner erscheint, in dem wir die Shell gestartet haben. Wir können auch direkt Befehle ausführen, ohne die Docker-Shell öffnen zu müssen. Der folgende Befehl wandelt einen Layer unseres GeoPackage in ein Shapefile um.

docker run \
  -v $(pwd):/data osgeo/gdal:alpine-small-latest \
  ogr2ogr \
    -f "ESRI Shapefile" \
    /data/out.shp \
    /data/sample.gpkg \
    layer_1

Der Zugriff auf eine lokale Datenbank benötigt die Einstellung --network=host für unseren Docker-Befehl. Dies scheint jedoch nur für Linux-Computer zu funktionieren. Dieser Thread diskutiert Alternativen für Windows und Mac.

docker run \
  --network=host \
  -t \
  -i \
  -v $(pwd):/data \
  osgeo/gdal:alpine-small-latest \
  ogr2ogr \
    -f GeoJSON \
    /data/db_export.geojson \
    PG:"
    dbname=my_database
    host=localhost
    port=5432
    user=postgres
    password=postgres" \
    my_table

Bonus - Features

GDAL/OGR hat einige eher unbekannte Funktionen, die allerdings auch vermutlich nicht so häufig benötigt werden:

Viewshed Algorithmus: Ab GDAL Version 3.1 gibt es eine Funktion welche die sichbare Fläche auf einem digitalen Geländemodell von einem gewissen Betrachter aus berechnet

GeoCoding: Der SQLite Dialekt bietet eine Funktion zum Geocoden. Man kann also aus dem Namen eines Orts die Koordinate erhalten. Als Input nehmen wir diese CSV-Datei:

id,name,country
1,Salzburg,Austria
2,Heidelberg,Germany
3,Trondheim,Norway

Aus dieser Liste lassen wir uns die Koordinate berechnen und in ein GeoPackage speichern.

ogr2ogr \
  -f GPKG output.gpkg \
  geocode.csv \
  -dialect sqlite \
  -sql "
  SELECT
    id,
    name,
    country,
    ogr_geocode(name || ', ' || country) AS geometry
  FROM geocode"

Es gibt auch die Funktion ogr_geocode_reverse() mit der man aus einer Koordinate den Namen eines Ortes erhalten kann.

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.

Um die Entwicklung und Wartung von GDAL/OGR weiterhin zu ermöglichen, kann die Arbeit des Hauptenwicklers Even Rouault (Spatialys) mit einem Sponsoring unterstützt werden.

Es gibt auch ein GDAL T-Shirt zu kaufen. Der Erlös wird an OSGeo gespendet. Das ist die Dach-Organisation die GDAL und viele andere Projekte aus dem Geo-Umfeld unterstützt.

Links

Deutsch

• Eine 60min Live-Demo von Jakob Miksch auf der FOSSGIS 2020 zeigt wie man ogrinfo und ogr2ogr praktisch anwendet
• Der Vortrag von Jakob Miksch auf der AGIT 2019 gibt einen Überblick über die Kommandozeilenwerkzeuge von GDAL/OGR
• 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

Englisch

• Website von GDAL/OGR
• Sourcecode von GDAL/OGR
• Vortrag - Einführung in OGR mit Python von Nikolai Janakiev
• Vortrag - Einführung in GDAL mit Python von Hannah Augustin
• Vortrag - State of GDAL 2019 von Even Rouault
• Vortrag - An Introduction to GDAL von Robert Simmon
• ausführliches Intro zu GDAL mit einleitenden Geschichten
• Vortrag - GDAL Tips and Tricks
• Blogeinträge von Koko Alberti
• Blogeinträge von Hannes Kröger
• Buch GeoSpatial Power Tools (2012) von Tyler Mitchell
• Video Tutorials von “Making Sense Remotely”
• Twitter account @gdaltips mit vielen praktischen GDAL und OGR Befehlen
• Offener Kurs: Introduction to GDAL von GIS OpenCourseWare
• Präsentation zu GDAL/OGR Kommandozeilen Tools von Jakob Miksch

CheatSheets

• von Garret - glw
• von BostonGIS
• von Nathan Shaw
• von Derek Watkins - dwtkns
• von Lyzi Diamond

Überarbeitung dieses Artikels

• 2021-04-11: Anleitung für Docker
• 2020-05-10: Virtuelles Vektor Format (VRT)
• 2019-12-16: Zusatz zu Shapefile, SQL, GeoCoding, Viewshed, ogr2ogr Optionen, GeoJSON Download
• 2019-09-20: Erste Version des Artikels