Cloudnative Geospatial: DIY-Höhenlinien aus der Cloud – Teil 2

Nachdem wir uns im letzten Blogpost durch verschiedene DEMs gewühlt haben, widmen wir uns jetzt der technischen Umsetzung, aus DEM abgeleitete Höhenlinien «cloudnative» in eine Webmap zu bringen. Konkret ist das Ziel, aus den einzelnen Rasterdaten des Copernicus-DEM einen Cloudnative Geospatial-Datensatz mit Protomaps zu erstellen, mit dem wir in einer Webapplikation mit MapLibre GL JS on-the-fly Höhenlinien und Reliefschattierung generieren können.

Wir verfolgen dabei einen ähnlichen Ansatz wie Syncpoint, wobei wir im Einzelnen aber davon abweichen.

Grober Ablauf und Tooling

Um ans Ziel zu kommen, müssen wir die folgenden Schritte durchführen:

1. Alle Rasterdatensätze zu einem einzelnen Datensatz zusammenfügen
2. Umwandlung der effektiven Höheninformationen in RGB-encodierte Höheninformationen
3. Umwandlung ins PMTiles-Format
4. Visualisieren in einer Webapplikation

Wie immer bei solchen Datentransformationen bietet es sich an, GDAL zu nutzen. Gleichzeitig verwenden wir die Tools Rasterio und Rasterio rgbify für die RGB-Encodierung sowie das pmtiles CLI für die Umwandlung in das PMTiles-Format für Protomaps.

Anmerkung: Da sich diese Toolchains am einfachsten in einer Linux Umgebung nutzen lassen und ich ein glühender Verfechter vom Windows Subsystem for Linux bin, sind die folgenden Anweisungen jeweils auf eine Linux-Umgebung gemünzt. Theoretisch sollte alles auch in Windows-Umgebungen laufen, allerdings gibt’s dort erfahrungsgemäss zum Teil Schwierigkeiten bei gewissen Tools.

Schritt 1: Datensätze zusammenfügen

Um verschiedene Rasterdaten zu einem einzelnen File zusammenzufügen, bietet sich gdal_merge an. Damit wir nicht alle Files einzeln spezifizieren müssen, können wir alle benötigten Files in ein Textfile schreiben und diese als Inputparameter angeben:

ls -1 *.tif > files.txt

gdal_merge.py -o merged_ch.tif --optfile files.txt

Mit unseren 16 Rasterdateien und auf meinem PC mit 32GB RAM und einer relativ guten CPU dauert die Operation knapp 6 Sekunden. Das Resultat sieht dann so aus:

Die schwarzen Flächen im Nordwesten bzw. Südwesten sind zwei Raster-Kacheln, die wir weggelassen haben, weil sie für die Schweiz nicht relevant sind. gdal_merge setzt ihren Wert auf 0; dies liesse sich aber anders parametrisieren.

Schritt 2: Umwandlung in RGB-encodierte Höheninformationen

Damit wir die Höhenlinien und Reliefschattierungen (Hillshades) generieren können, müssen die absoluten Höhenwerte des DEMs in ein RGB-encodiertes Format umgewandelt werden. rio rgbify ist ein Tool aus dem Rasterio Werkzeugkasten, welches uns diese Umrechnung abnimmt. Konkret werden alle Höheninformationen in einen RGB-Farbwert umgewandelt:

The Mapbox Terrain-DEM tileset contains Terrain-RGB tiles. Terrain-RGB tiles include elevation data that is encoded using each color channel as a position in a base-256 numbering system. This approach allows for 16,777,216 unique values which can be mapped to 0.1 meter height increments, enabling vertical precision necessary for cartographic and 3D applications.

Quelle: Mapbox

Der Aufruf ist wie folgt:

rio rgbify -b -10000 -i 0.1 merged_ch.tif merged_ch_rgb.tif

-i 0.1 definiert dabei die Präzision der Codierung, womit wir (theoretisch) 10 cm genaue Rasterzellen erhalten. -b -10000 definiert die Basis, von der aus gerechnet wird; dies lässt sich dem Abschnitt «Decode Data» aus der Dokumentation entnehmen.

Das Resultat dieses Aufrufs sieht wie folgt aus:

Die dunkelblauen Artefakte in den Schweizer Alpen machen zunächst stutzig; allerdings sind dies Höhen über 3000m, die aufgrund der Codierung einen Bruch in der Farbgebung darstellen.

Schritt 3: PMTiles generieren

Als nächstes muss dieses encodierte Raster in MBTiles umgewandelt werden, bevor es ins PMTiles-Format überschrieben werden kann. In einem ersten Versuch habe ich das mit gdal_translate gemacht:

gdal_translate merged_ch_rgb.tif merged_ch_rgb.mbtiles

Allerdings gibt es dann im Endprodukt Artefakte, weil gdal_translate (bzw. das im Prozess verwendete gdal2tiles) ein Resampling vornimmt, was dann natürlich falsche Resultate liefert: Da unser Raster nun RGB-Werte und keine effektiven Höheninformationen mehr aufweist, produziert ein Resampling (im obigen Beispiel mit dem Default average) Artefakte, weil wir keine graduellen Werte mehr haben:

Mit einer anderen Resampling Methode (z.B. mit near, via -r near) liesse sich das umgehen, würde aber wiederum unschöne Darstellungen bei hohen Zoomstufen generieren (vgl. hier). Letztenendes habe ich dann den vorherigen Schritt mit der Umwandlung in MBTiles durch den Einsatz von rgbify kombiniert und einen sauberen, artefaktfreien Datensatz erhalten – so kann man sich auch den Zwischenschritt mit dem TIFF sparen. Der Aufruf ist dann wie folgt:

rio rgbify -b -10000 -i 0.1 merged_ch.tif merged_ch_rgb.mbtiles --max-z=13 --min-z=10

Die letzten beiden Parameter definieren dabei die Zoomstufen, die wir kacheln wollen. Je mehr und je detaillierter man dies benötigt, desto länger dauert die Umwandlung und desto grösser wird das File. Für unseren Testlauf reichen diese Zoomstufen.

Je nach gewählten Zoomstufen dauert die Operation eine ganze Weile. Anschliessend fehlt noch das Umschreiben in das PMTiles-Format:

pmtiles convert merged_ch_rgb.mbtiles merged_ch_rgb.pmtiles

Da es sich hier um ein simples Umschreiben handelt, geht dies in der Regel etwas schneller.

Schritt 5: Visualisierung

Damit sind wir eigentlich schon an dem Punkt angelangt, an dem wir die Daten visualisieren können. In unserer lokalen Umgebung bietet PMTiles einen einfachen Server an, welcher zur Visualisierung genutzt werden kann (mit ausgehebeltem CORS für lokale Zwecke):

pmtiles serve --cors=* .

Zur einfachen Visualisierung halten wir uns an das Add Contour Lines-Beispiel von MapLibre GL JS, bauen aber noch einen OpenStreetMap-Kartendienst ein, damit wir uns orientieren können. Das Endresultat lässt sich sehen: Es beinhaltet sowohl die Höhenlinien als auch Hillshading, und alles wird on-the-fly generiert!

In den folgenden Beispielen ist jeweils links nur der Kartenausschnitt angezeigt, auf der rechten Seite sind die on-the-fly generierten Hillshade und Höhenlinien überlagernd dargestellt.

Schritt 5: Bonus – 3D Visualisierung

Noch interessanter ist aber, dass wir nun auch 3D Terrains darstellen können und der Aufwand dafür praktisch bei Null liegt. MapLibre GL JS nimmt dies wiederum on-the-fly im Client vor.

Ausblick

Am Ende dieser Schritte wurde aus 16 einzelnen Rasterbildern mit gesamthaft 643 MB Grösse ein einziges PMTiles-File mit 3 GB Grösse generiert. Alles, was jetzt noch übrig bliebe, um wirklich in die Cloud zu gehen, wäre, dieses File auf einen Object Storage wie S3 zu laden. Dann lässt sich das File direkt einbinden, z.B. mit dem offiziellen MapLibre GL JS pmtiles Protokoll.

Noch spannender und truly cloudnative ist aber ein CDN-Deployment: Dabei wird das File ebenfalls in einen Bucket im Object Storage gespeichert, aber mit einer Serverless Function (aktuell für Cloudflare und AWS) in einen (scheinbar) echten Service umgewandelt. Damit lässt sich das File genauso einbinden, wie wir das in Schritt 4 und 5 gemacht haben, nämlich mit normalen https://tiles.myglobalcdn.com/myvectordata/{z}/{x}/{y}.png-URLs. Damit braucht man auch keine zusätzlichen Plugins oder Hacks mehr, sondern man kann ganz normale Layer einbinden. Durch die CDNs gibt es zudem vielfältige Caching-Methoden, um niedrige Latenzen zu erreichen.

Optimierungen und Einschränkungen

Das hier beschriebene Vorgehen klappt gut für begrenzte Datensätze. Werden mehr Daten benutzt, stossen Desktop-PCs schnell an ihre Grenzen. Ein Problem, in das man relativ schnell läuft, zeigt sich in Schritt 2: gdal_merge lädt alle Daten in den Arbeitsspeicher und nimmt dann den Merge vor. Das funktioniert natürlich nicht, sobald eine Vielzahl an TIFFs gemerged werden muss. Abhilfe verschafft hier aber ein weiterer Zwischenschritt, indem man zunächst ein VRT erstellt und dieses dann mit gdal_translate in ein einzelnes Raster schreibt:

gdalbuildvrt -input_file_list ./input_files.txt ./virtual_input.vrt

gdal_translate -of GTiff --config GDAL_CACHEMAX 8192 --config GDAL_NUM_THREADS ALL_CPUS --config GDAL_MAX_DATASET_POOL_RAM_USAGE 18000 --config GDAL_MAX_DATASET_POOL_SIZE 500  virtual_input.vrt ./merged.tif

Eine andere Optimierung gäbe es hinsichtlich der Daten an sich. Da wir dies in unserem Beispiel im Client machen, haben wir keine effektiven Vektoren, die wir wie bei MapTiler- oder Mapbox-Daten stylen und verwenden können. Eine einfache Möglichkeit, die Höhenlinien als Vektor zu speichern, wäre PostGIS: Dabei würden wir alle Rasterdaten nach PostGIS importieren und mit ST_Contour die Höhenlinien generieren, zum Beispiel als materialisierte View. Mit ogr2ogr liessen sich diese Daten direkt als MBTiles exportieren und dann weiterverarbeiten. Dann würde sich natürlich auch die Grösse des Endproduktes verringern, da wir dann keine Rasterdaten mehr kacheln müssen.

Fazit

Das hier gezeigte Beispiel hat überraschend gut und einfach funktioniert. Wenngleich dieser Case selten auftreten dürfte (weil es genügend gute Services gibt, die diese Daten liefern), ist das Vorgehen einfach adaptierbar.

Insbesondere in Kombination mit sensiblen Daten bieten sich hier natürlich viele Möglichkeiten: Beispielsweise wenn wir sensible Vektordaten haben, können wir diese relativ einfach über einen internen Object Storage (wie etwa Minio) in eine Applikation integrieren, ohne dass wir dafür Drittdienste benutzen müssen, wo unter Umständen Privacy-Bedenken bestehen.

Gerade für Daten, die sich selten ändern, ist Protomaps ein unglaublich gutes Framework, weil es die Kosten drastisch reduzieren kann und gleichzeitig die modernen Webtechnologien nutzbar und effizient einsetzbar macht.

Hast Du einen Use Case, für den das cloudnative-Paradigma vielleicht wertvoll sein könnte? Ich würde mich über eine Kontaktaufnahme freuen.