Wie man eine schöne Schummerung herstellt (Teil 2)

SRTM Szene in QGIS
SRTM Szene in QGIS

Schöne Schummerung Teil 2

In dieser Serie von Blogposts in drei Teilen beschreibe ich, wie man mittels OpenSource-Werkzeugen optisch ansprechende Schummerungen herstellen kann. Ausgehend von einem Geländemodell wie einer SRTM-Szene, kann man mittels verschiedener Werkzeuge eine Schummerung herstellen, die den generischen Schummerungen von gdaldem oder von GRASS GIS r.shaded.relief überlegen ist. Die grundsätzlichen Schritte beinhalten Entrauschen des Geländemodells und das Verschneiden von Hangneigungs- und Reliefkarte.

Im ersten Teil habe ich die notwendige Werkzeuge und die Daten vorgestellt, die wir für das Erstellen einer schönen Schummerung benötigen. In diesem Teil werden wir das Geländemodell auf den Einsatz als Schummerung vorbereiten indem wir fehlerhafte Daten und eventuell vorhandene Muster aus dem Modell herausfiltern. Außerdem erstellen wir eine Standardschummerung mit GDAL um am Ende einen Vergleich ziehen zu können. Am Ende dieses Posts werden wir noch die notwendigen Einzelteile für eine schöne Schummerung erstellen.

Die Vorbereitung

Um an Ende einen Vergleich ziehen zu können zwischen einer normal erstellten Schummerung und unserer schönen Schummerung möchte ich am Anfang einmal eine normale Schummerung berechnen. Ich habe im ersten Post bereits ein Bild mit einer Schummerung gezeigt. Diese Schummerung war überlagert mit Höheninformationen (von grün – niedrig nach rot – hoch). Um aber einen neutraleren Vergleich ziehen zu können, rechnen wir nun eine Schummerung ohne die Überlagerung mit Höheninformationen. Wie im ersten Post bereits beschrieben öffnen wir dazu eine OSGeo4W-Shell im Arbeitsordner, in dem das Geländemodell liegt. Bevor wir die erste Schummerung berechnen, möchte ich das Geländemodell in ein kartesisches Koordinatensystem transformieren. Das hat mehrere Gründe, hauptsächlich aber weil das Programm zum Entrauschen der Satellitenszene ein kartesisches Koordinatensystem benötigt. Das Transformieren in ein anderes Koordinatensystem erledigt man mit dem Kommando gdal_translate. Bevor wir das aber starten, müssen wir noch ein geeignetes metrisches Koordinatensystem finden. Hier bietet sich UTM an, das Universal Transverse Mercator System, an.

UTM Umprojektion

Der Dateiname des Geländemodells, N42W001, gibt die Koordinaten der linken unteren Ecke der Szene an, hier 42° Nord und 1° West. Um jetzt die passende UTM-Szene zu finden, gibt es eine Vielzahl an Werkzeugen. In diesem Fall nutzen wir exemplarisch diesen Online-Konverter. In das Feld Latitude geben wir 42.5 an, in das Feld Longitude -0,5 um das Zentrum der Szene zu definieren. Wir erhalten die UTM-Zone 30 (T). Nun benötigen wir noch den EPSG-Code für dieses Koordinatensystem, das ist ein eindeutiger Schlüssel zur einfachen Verwendung des Koordinatensystems. Den finden wir bei Spatial Reference, einer Online-Bibliothek zum Nachschlagen von EPSG-Codes von Koordinatensystemen. Der UTM-Streifen 30N (N für nördlich des Äquators) hat den EPSG-Code 32630. Mittels gdalwarp nehmen wir nun die Umprojektion vor, wobei die einzelnen Parameter folgendes bedeuten:

  • -t_srs: Ziel-Koordinatensystem
  • -r: Bildumrechungsmethode
  • -dstnodata: Wert für Bereiche ohne Daten
> gdalwarp -t_srs "EPSG:32630" -r bilinear -dstnodata -9999 N42W001.hgt n42w001_epsg32630.tif

In QGIS sieht das umgerechnete Höhenmodell wie folgt aus:

Das Geländemodell in UTM Zone 30N umgerechnet
Das Geländemodell in UTM Zone 30N umgerechnet

Schummerungsberechnung

Nun können wir endlich die einfache Version einer Schummerung berechnen. Dazu nutzen wir das gerade eben umgerechnete Höhenmodell. Mittels gdaldem nehmen wir die Schummerungsberechnung folgender Maßen vor (der Parameter -z definiert eine Überhöhung des Geländes):

> gdaldem hillshade -z 1.5 n42w001_epsg32630.tif simple_hillshade.tif

In QGIS sieht die Schummerung wie folgt aus:

Eine einfache Schummerung
Eine einfache Schummerung

Vor allem im nördlichen Teil der Szene fallen weiße Stellen auf. Das sind Regionen, in denen vom Satelliten keine Daten aufgenommen werden konnten. Bevor wir weiter Fortschreiten mit der Verwendung der Satellitenszene müssen wir diese Stellen mit Daten füllen. Dazu gibt es mehrere Methoden: man kann andere Geländedaten benutzen um die Löcher zu füllen oder man kann die Löcher mittels Interpolation füllen. In diesem Fall werden wir wieder GDAL verwenden, genauer gdal_fillnodata.

Lücken füllen

Wir beginnen wir noch einmal mit der originalen Szene, um das Füllen von Nullwerten außerhalb des Datenbereiches zu vermeiden. Anschließend projizieren wir die Szene wieder in das UTM-System und überschreiben die weiter oben erstellte.

> gdal_fillnodata N42W001.hgt n42w001.tif
> gdalwarp -t_srs "EPSG:32630" -r bilinear -dstnodata -9999 n42w001.tif n42w001_epsg32630.tif

Berechnung der benötigten Karten

Nun gehen wir hin und berechnen endlich die benötigten Karten um die schöne Schummerung zu erstellen. Der erste Schritt ist das Entrauschen der Satellitenszene um Unregelmäßigkeiten, kleinere Strukturen und eventuelle Muster des Satelliten zu entfernen.

Entrauschen mit Mdenoise

Das Entrauschen ist ein dreiteiliger Prozess, der aus der Konvertierung des Geländemodells in ein für das Entrauschen geeignetes Format, dem Entrauschen selber und dem Zurückumwandeln in das Urpsrungsformat besteht. Das für Mdenoise benötigte Format nennt sich AAIGrid, ArcInfo ASCII Grid. Zu beachten ist bei der Konvertierung mit GDAL, dass man ein Ziel- und Startkoordinatensystem angibt, da es bei älteren Versionen von GDAL zu einem Versatz des Geländemodells von einem halben Pixel kommt. Dieser Fehler ist in Version 2.02 behoben. Beim Entrauschen halte ich mich an die Parameter der Anleitung, aber Experimentieren ist erlaubt! Zu beachten ist der Parameter -z, der die Glättung auf die vertikale Richtung beschränkt. Der dreiteilige Prozess sieht dann wie folgt aus:

> gdal_translate -a_srs "EPSG:32630" -of AAIGrid n42w001_epsg32630.tif n42w001_epsg32630.asc
> mdenoise -i n42w001_epsg32630.asc -n 5 -t 0.99 -z -o n42w001_denoised.asc
> gdal_translate -a_srs "EPSG:32630" n42w001_denoised.asc n42w001_denoised.tif

Zwei Schummerungskarten

Wie im ersten Teil der Serie beschrieben wollen wir zwei Schummerungskarten miteinander verschneiden, um das beste aus beiden Versionen herauszuholen und eine plastische Schummerung zu erstellen. Die zwei Schummerungskarten werden mit unterschiedlichen Parametern gerechnet. Für die erste Karte sollen uns die Standardwerte für Überhöhung, Azimuth und Sonnenhöhe reichen, Azimuth 315, Sonnenhöhe 45 und Überhöhung 1. Für die zweite Karte bleiben Azimuth und Sonnenhöhe unverändert, die Überhöhung beträgt aber lediglich 0.25. Der Parameter -co TFW=YES sagt GDAL, für das GeoTIFF ein Worldfile zu erzeugen. Da Imagemagick die Geodaten im TIFF nicht versteht, müssen wir die Georeferenzierung aus der Datei in eine kleine Textdatei auslagern um die Lage und Ausdehnung des Geländemodells zu erhalten.

> gdaldem hillshade -z 1 -az 315 -alt 45 -co TFW=YES n42w001_denoised.tif nice_hillshade_map1.tif
> gdaldem hillshade -z 0.25 -az 315 -alt 45 -co TFW=YES n42w001_denoised.tif nice_hillshade_map2.tif

Die beiden erzeugten Schummerungskarten sollten dann wie im Bildschirmfoto gezeigt aussehen. Gemacht habe ich die Screenshots aus QGIS heraus mit der Ansicht 100%, auf eigene Auflösung zoomen. Diese Option erreicht man, indem per Rechtsklick auf den Layer im Layerfenster klickt.

Wenn man zwischen den Screenshots hin- und herwechselt, kann man in kleinen Details bereits Unterschiede erkennen, die in dem jeweils anderen nicht sichtbar sind.

Abschluss

In diesem Post haben wir gelernt, wie man Lücken im Geländemodell schließen kann, wie man Störungen und Artefakte herausrechnen kann und wie man die Parameter zur Schummerungsberechnung anpassen muss, um zwei verschiedene Schummerungskarten zu erhalten, die unterschiedliche Details betonen. Außerdem haben wir gelernt, wie man Karten in ein anderes Koordinatensystem umwandeln kann, wie man Bereiche ohne Daten klassifizieren kann und wie man die richtige UTM Zone bestimmt und die passende EPSG-Id herausbekommt.

Die verwendeten Kommandos noch einmal im Überblick:

  • gdal_translate
  • gdalwarp
  • gdaldem
  • mdenoise

Im nächsten Post werden wir die beiden Schummerungskarten zu einer einzigen verschneiden und einen Vergleich mit der herkömmlich erzeugten Schummerungskarte ziehen. Außerdem werde ich noch andere Methoden kurz anreißen, mit denen man auch verbesserte Schummerungskarten herstellen kann.