Abgrenzung von Wassereinzugsgebieten mit Pysheds

Das Einzugsgebiet eines Flusses ist die Grenze, innerhalb derer ein Wassertropfen an der gleichen Mündung in den Ozean abfließt¹. Ganz allgemein kann ein Wassereinzugsgebiet durch einen weiter flussaufwärts gelegenen "Schüttungspunkt" (die Mündung im Falle eines ganzen Flusses) definiert werden, was bedeutet, dass man, wenn man Niederschlagsdaten innerhalb des Wassereinzugsgebiets abfragt, Dinge wie das Volumen/die Fließgeschwindigkeit des Wassers am Schüttungspunkt vorhersagen kann.

Entscheidungsträger in verschiedenen Bereichen könnten Anwendungen für die Abgrenzung eines Wassereinzugsgebiets finden. In der Stadtplanung kann es beispielsweise nützlich sein, das Einzugsgebiet eines Flusses zu ermitteln, um die Vorhersage von Überschwemmungen zu unterstützen, während die Betreiber von Wasserkraftwerken die Wasserscheide nutzen können, um die Energielieferung abzuschätzen. In diesem Beitrag zeige ich, wie man mit den hochaufgelösten Höhendaten von Meteomatics, die über unsere API verfügbar sind, Einzugsgebiete und Durchflussraten für beliebige Gewässer bestimmen kann.

PySheds

Die pysheds-Bibliothek enthält Funktionen/Methoden, die es ermöglichen, .tiff-Dateien zu lesen und Wassereinzugsgebiete abzugrenzen. Vielleicht haben Sie keinen Zugang zu einer .tiff-Datei, aber als Meteomatics-Abonnent können Sie unsere API verwenden, um ein digitales Höhenmodell (DEM) zu erhalten und die Grenzen einer Wasserscheide selbst zu bestimmen. Ich zeige Ihnen, wie das geht und wie Sie die Niederschlagsdaten für das Einzugsgebiet abrufen können.

Hier sehen Sie die Struktur meiner API-Anfrage für Höhendaten. Standardmäßig werden die Daten hier aus dem NASA SRTM/ASTER-Modell bezogen, das eine Bodenauflösung von etwa 90m hat. Es handelt sich um dasselbe Modell, das wir für die Ausgabe grob aufgelöster Wettermodelle verwenden, so dass Daten, auf die über die Meteomatics-API zugegriffen wird, entsprechend der lokalen Geografie einer Anfrage aktualisiert werden. "Bounds" ist ein Wörterbuch, das die Randkoordinaten der Region enthält, an der ich interessiert bin, und "res" ist eine Variable, die die Auflösung der zu erfassenden Daten bestimmt.

Ich würde eine möglichst hohe Auflösung empfehlen, denn je mehr Informationen über die lokale Neigung man aus der API extrahieren kann, desto besser wird das Bild der Wasserscheide aussehen. Abgesehen davon:

a) Einige Flüsse haben riesige Einzugsgebiete, die bei hoher Auflösung die maximale Anzahl von Datenpunkten in einer einzigen API-Anfrage überschreiten können²

b) Eine Auflösung von 0,001 Grad entspricht in etwa 100m, sodass es unnötig ist, Daten mit einer höheren Auflösung als dieser anzufordern, da sie kleiner sind als unsere herunterskalierten Daten

Das datetime-Objekt, das ich als erstes Argument der Anfrage übergebe, ist willkürlich, da wir nicht erwarten, dass sich die Höhe im Laufe der Zeit sehr stark ändert.

Als Nächstes ersetze ich alle Orte, an denen die Höhe Null ist, durch NaN und wandle die von der API abgerufenen Daten in ein netCDF um (ein Industriestandardformat, das für den nächsten Schritt erforderlich ist). Da 0m Höhe gleichbedeutend mit "auf Meereshöhe" ist, wird das Meer aus unseren Daten entfernt. Das ist wichtig, denn sonst würde das Wasser ewig über diese vollkommen ebene Fläche fließen, was eindeutig unphysikalisch ist.

elev[elev <= 0.0] = np.nan
da = elev.T.unstack().to_xarray()
nc = da.to_dataset(name='elevation')

Nun, da wir ein netCDF von elevation haben, müssen wir das rioxarray Paket installieren. Dieses Paket arbeitet mit xarray zusammen (ein Paket zur Manipulation von netCDF-Daten, das Sie ebenfalls installieren sollten), um das Schreiben von TIFF-Dateien zu erleichtern. Nach der Installation dieser Pakete mit dem

`pip install -c conda-forge <package>`

Befehl, kann man folgendes implementieren:

nc.rio.set_spatial_dims( 'lon', 'lat', inplace=True)
nc.rio.set_crs("epsg:4326", inplace=True)
nc.rio.to_raster('elevation.tif') 
grid = pysheds.grid.Grid.from_raster('elevation.tif')
dem = grid.read_raster('elevation.tif')

So erhalten Sie ein digitales Höhenmodell in dem von pysheds gewünschten Format. Werfen wir einen Blick auf unser DEM.

Wir sind nun fast bereit, eine Wasserscheide abzugrenzen. Zuvor müssen wir uns jedoch um alle kleinen Artefakte kümmern, die in unserem DEM vorhanden sein könnten. pysheds bezeichnet diese Artefakte als Gruben und Vertiefungen, d. h. Orte, an denen die gesamte Umgebung höher liegt als das Artefakt (daher kann sich dort Wasser ansammeln), und als Ebenen, d. h. Bereiche mit konstanter Höhe. Die schwierigste (und unrealistischste) Ebene - das Meer - haben wir bereits entfernt. Die anderen sowie die Gruben und Senken lassen sich mit den folgenden Methoden für das DEM-Objekt sehr leicht entfernen:

dem = grid.fill_pits(dem)
dem = grid.fill_depressions(dem)
dem = grid.resolve_flats(dem)

Wir sollten nun überprüfen, ob wir alle Artefakte wie oben beschrieben erfolgreich entfernt haben

assert not grid.detect_pits(dem).any()
assert not grid.detect_depressions(dem).any()
assert not grid.detect_flats(dem).any()

Wenn Sie hier Fehler bei der Behauptung erhalten, ist etwas schief gelaufen. Ich stellte fest, dass alle diese Fehler verschwanden, nachdem ich den Ozean korrekt aus meinen Daten maskiert hatte, aber wenn die Fehler weiterhin auftreten, müssen Sie das DEM im Detail untersuchen, um zu sehen, wo etwas schief läuft. Glücklicherweise hindert uns das Scheitern dieser Behauptungen nicht daran, die nächsten Schritte zu berechnen, und die Ergebnisse dieser Schritte können nützlich sein, wenn man versucht zu analysieren, wo der Fehler liegen könnte.

Wir können nun die Fließrichtung und die Akkumulation der DEM wie folgt berechnen:

fdir = grid.flowdir(dem)
acc = grid.accumulation(fdir)

Die Fließrichtung ist eine Zahl zwischen 0 und 360, die die Fließrichtung des Wassereinzugsgebiets an jedem gegebenen Punkt angibt. Die Kumulierung quantifiziert die Zellen im DEM, die der betreffenden Zelle vorgelagert sind. Diese Eigenschaft ist aus mehreren Gründen nützlich, wie wir sehen werden. Der folgende Code erzeugt eine Karte der Flüsse in unserem DEM, die ich in Abbildung 2 zeige.

f = plt.figure()
ax = f.add_subplot(111)
f.patch.set_alpha(0)
im = ax.pcolormesh(
 nc.lon, nc.lat, masked_acc, zorder=2, 
 cmap='cubehelix', norm=colors.LogNorm(1, acc.max())
)
plt.colorbar(im, ax=ax, label='Upstream Cells')
plt.title(
 f'Flow Accumulation in catchment\n of pour point at {x_snap},{y_snap}', 
 size=14
)
plt.tight_layout()
plt.show()

Die Farbskala des Bildes ist die Variable acc. Man sieht, dass Zellen, die sich an den Hängen von Tälern befinden, sehr wenige stromaufwärts gelegene Zellen haben - ihnen fließt kaum etwas zu -, während die Akkumulation von Flüssen zunimmt, je weiter sie in Richtung Küste fließen - jede Zelle stromabwärts wird von mindestens einer stromaufwärts gelegenen Zelle gespeist. Das bedeutet, dass wir die Akkumulation als Proxy dafür verwenden können, wie weit wir flussabwärts sind, und wir können auch einen vernünftigen Wert für die Akkumulation festlegen, um fließende Flüsse von Hängen zu unterscheiden.

Sie können sehen, dass unsere Flüsse, wie wir bereits besprochen haben, nicht wirklich wissen, wie sie sich verhalten sollen, wenn sie auf den Ozean treffen. Das ist unwichtig, da die Strömung hier nicht quantifiziert wird, und wir können diese Teile des Bereichs einfach ausblenden.

Nun ist es ganz einfach, ein Einzugsgebiet für einen Fluss abzugrenzen. Wir brauchen lediglich die Koordinaten des "Einlaufpunktes". Das kann der Punkt sein, an dem der Fluss ins Meer mündet, oder ein weiter flussaufwärts gelegener Punkt, an dem man interessiert ist. Für mein Beispiel in Montpellier interessiere ich mich für die Stelle, an der der Lez und der Mosson zusammenfließen. Das ist ungefähr bei 3.81, 43.56. Wir geben die Koordinate in unseren Code ein, achten aber auch darauf, dass wir einen nahegelegenen Punkt mit ausreichender Akkumulation einrasten lassen: Es wäre peinlich, wenn das Ergebnis des nächsten Schritts uns ein winziges Gebiet zeigen würde, weil wir versehentlich einen nahegelegenen Hügel statt eines Punktes im Fluss gewählt haben...

# Bestimmen Sie den Gießpunkt
x, y = 3.8095, 43.5588
 
# Gießpunkt an hoher Akkumulationszelle fangen
x_snap, y_snap = grid.snap_to_mask(acc > 10000, (x, y)) 
catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, xytype='coordinate')
masked_acc = np.ma.masked_array(acc, catch==False)

Dann ermitteln wir das Einzugsgebiet mithilfe der catch-Methode des Grid-Objekts. Dies ist ein boolesches Feld, das beschreibt, ob eine Zelle im Einzugsgebiet liegt oder nicht. Ich maskiere meine Akkumulation entsprechend dem Einzugsgebiet und zeige das Ergebnis unten

Aus diesem Grund ist es einfach, das Gebiet nach der Entfernung zum Auslasspunkt zu kategorisieren

dist = grid.distance_to_outlet(
 x=x_snap, y=y_snap, fdir=fdir, xytype='coordinate'
)

Und es braucht nicht viel mehr, um die Region danach zu kategorisieren, wie viel Zeit das Wasser in jeder Zelle brauchen würde, um zum Schüttpunkt zu gelangen. Alles, was wir tun müssen, ist, eine einfache Annahme zu implementieren - sagen wir, dass Wasser in einem Fluss (d.h. eine Zelle mit vielen stromaufwärts gelegenen Zellen) etwa 10x so schnell fließt wie Wasser an einem Hang

>

Gewichte = acc.copy()

gewichte[acc >= 100] = 0.1

gewichte[(0 < acc) & (acc < 100)] = 1.

Das pysheds-Paket ist potenziell sehr nützlich für jede Anwendung, die mit Flusspegeln und -strömungen zu tun hat. Einige davon werden in diesem Artikel gezeigt oder beschrieben, aber es gibt noch viele mehr, die von Nutzern mit Zugang zu einem so hochwertigen DEM, wie es den API-Abonnenten von Meteomatics zur Verfügung steht, genutzt werden können.

Informationen über unsere zusätzlichen geografischen Daten finden Sie hier. Wie immer, wenn Sie Fragen zur Nutzung unserer Daten oder zu den in diesem Artikel vorgestellten Techniken haben, können Sie mich gerne über das untenstehende Formular kontaktieren!