Update: Für eine verbesserte Version siehe mein Jupyter Notebook.
Mit Python und PyGMT ist es ziemlich einfach, eine einfache Topokarte zu zeichnen, die auf den Höhendaten der NASA (Public Domain!) beruht. Allerdings beruht der einfachste Ansatz darauf, dass eine colormap alles unter 0 m NN in Blautönen, alles über 0 m NN in Grün-/Brauntönen einfärbt (z.B. mit cmap=“geo“). Bei einer Weltkarte oder Bergen sieht das Ergebnis gut aus, bei Küstenlinien gibt es aber Probleme. Eine Karte von Holland macht das Problem deutlich:
import pygmt
region=[3.2, 6.8, 51.4, 53.5]
grid = pygmt.datasets.load_earth_relief(resolution="03s",
region=region)
fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="M15c", cmap="geo")
fig.colorbar(frame=["a10", "x+lElevation", "y+lm"])
fig.coast(shorelines=True, region=region, resolution="f")
fig.show()
Land unter! In den Tutorials und Docs habe ich keine Anleitung gefunden, wie das Problem gelöst werden kann. Mithilfe von pygmt.grdlandmask war es dann ganz einfach. Der Ansatz beruht auf der Küstenlinie der Funktion fig.coast(). Die erzeugte Maske setzt für jedes Pixel den für „Land“ bzw. „Wasser“ festgelegten Wert, z.B. 0 oder 1. Der erste Trick ist, dass diese Maske einfach mit dem grid multipliziert werden kann, dann werden alle Wasser-Pixel auf 0 gesetzt. Jetzt können wir das Land mit einer geeigneten colormap zeichnen, z.B. cmap=“dem1″. Als spacing für die Maske verwende ich den Wert der resolution des verwendeten Datasets; resolution der Maske ist „f“ (full resolution), wie bei fig.coast().
(Update: Anstatt die Maske mit pygmt.grdlandmask für die Region zu berechnen, kann ab pygmt 0.9 mit pygmt.datasets.load_earth_mask ein fertiges Maskengrid heruntergeladen werden, mit 0 für Ozean, 1 für Land, 2 für See, 3 für Insel, 4 für „Pond“.)
Das Wasser plotten wir im zweiten Schritt darüber. Allerdings können wir jetzt nicht einfach „maskvalues=[1, 0]“ verwenden, weil dann das Land zu Flachwasser der Tiefe 0 wird. Es ist aber zum Glück möglich, bei der Maske „NaN“ (als string) als Wert anzugeben. Und in fig.grdimage können wir mit nan_transparent=True sicherstellen, dass das Land nun transparent ist.
land = grid * pygmt.grdlandmask(region=region,
spacing="03s",
maskvalues=[0, 1],
resolution="f")
wet = grid * pygmt.grdlandmask(region=region,
spacing="03s",
maskvalues=[1, "NaN"],
resolution="f")
fig = pygmt.Figure()
# Plotte Land (mit Hillshading)
fig.grdimage(grid=land, projection="M15c",
cmap="dem1", shading=True)
fig.colorbar(frame=["a10", "x+lElevation", "y+lm"])
# Plotte Wasser (NaN transparent)
fig.grdimage(grid=wet,
cmap="seafloor", nan_transparent=True)
fig.coast(shorelines=True, projection="M15c",
region=region, resolution="f")
fig.coast(rivers="1/0.5p,blue") # Flüsse
fig.show()
