Script a Habitat Suitability Map

Bereiten Sie sich durch Klonen des Übungs-Repository:

git clone https://github.com/Ecohydraulics/Exercise-geco.git

Sacramento suckers in the South Yuba River (source: Sebastian Schwindt @hydroinformatics on YouTube).

Was ist Habitat Suitability?¶

Fisch und andere Wasserarten ruhen, orientieren und reproduzieren sich in einer fluvialen Umgebung, die ihren physischen Lebensraum darstellt. Während ihrer verschiedenen Lebensstufen weisen verschiedene Fische spezifische physikalische Lebensraumpräferenzen auf, die beispielsweise in Abhängigkeit von Wassertiefe, Strömungsgeschwindigkeit und Korngröße des Flussbettes definiert sind. Der sogenannte Habitat Suitability Index $HSI$ kann für hydraulische (Wassertiefe oder Strömungsgeschwindigkeit) und morphologische (z.B. Korngröße oder Abdeckung in Form von großem Holz) Parameter individuell berechnet werden, um die Qualität des physikalischen Lebensraums für einen Fisch und in einem bestimmten Lebensstadium zu beschreiben. Die folgende Abbildung zeigt beispielhafte $HSI$ Kurven für die fry, juvenile und erwachsene Lebensstadien der Regenbogenforelle in Abhängigkeit von der Wassertiefe. Die $HSI$-Kurven sehen in jedem Fluss unterschiedlich aus und sollten von einem aquatischen Ökologen individuell aufgestellt werden.

Das $HSI$-Konzept enthält auch den sogenannten Deckungsraum in Form des cover Habitat Suitability Index $HSI_{cov}$. Der Lebensraum ist das Ergebnis lokaler Turbulenzen, die durch Rauhigkeitselemente wie Holz, Boulder oder Brückenpier verursacht werden. Bei dieser Übung werden wir uns jedoch nur mit hydraulischen Lebensraumeigenschaften befassen (nicht mit Lebensraum).

Erwachsene Forellen schwimmen in Deckung Lebensraum durch eine Brücke Pier im oberen Neckar River geschaffen.

Die Kombination mehrerer $HSI$-Werte (z.B. wassertiefenbezogen $HSI_{h}$, strömungsgeschwindigkeitsbezogen $HSI_{u}$, Korngrößenbezogen $HSI_{d}$ und/oder Cover $HSI_{cov}$) ergibt den *kombinierten Habitat Suitability Index$cHSI$. Es gibt verschiedene Berechnungsmethoden für die Kombination verschiedener $HSI_{par}$-Werte in einen $cHSI$-Wert, wobei der geometrische Mittelwert und das Produkt die am weitesten verbreiteten deterministischen Kombinationsmethoden sind:

Sind daher aus einem zweidimensionalen (2d) hydrodynamischen Modell die pixelbasierten $HSI$-Werte für Wassertiefe und Strömungsgeschwindigkeit bekannt, so kann für jeden Pixel der $cHSI$-Wert entweder als Produkt oder geometrisches Mittel des Einzelparameters $HSI_{par}$-Rasters berechnet werden.

Dieses Konzept für die Lebensraumbewertung wurde von Bovee (1986) und Stalnaker et al. (1995)(direct download) eingeführt. Diese Autoren bauten jedoch ihre nutzbare (physikalische) Lebensraumbewertung auf Basis von eindimensionalen (1d) numerischen Modellen, die im letzten Jahrtausend häufig verwendet wurden. Heutzutage sind 2d numerische Modelle die modernsten, um geospatiell explizit auf Pixelbasis $cHSI$-Werten physikalischen Lebensraum zu bestimmen. Es gibt zwei verschiedene Möglichkeiten zur Berechnung des nutzbaren Lebensraums ($UHA$) basierend auf Pixel-basierten $cHSI$-Werten (und noch mehr Möglichkeiten finden Sie in der wissenschaftlichen Literatur).

Eine Alternative zur deterministischen Berechnung der $HSI$ und $cHSI$-Werte eines Pixels ist ein Fuzzy-Logik-Ansatz Noack et al., 2013. Im Fuzzy-Logik-Ansatz werden Pixel z.B. als low, medium oder high Lebensraumqualität in Abhängigkeit von der zugehörigen Wassertiefe oder Strömungsgeschwindigkeit mit kategorischen (low, medium, oder high), expert assessment-based $HSI$Kurven klassifiziert. Der $cHSI$-Wert ergibt sich aus dem Schwerpunkt der überlagerten Mitgliedsfunktionen der betrachteten Parameter (z.B. Wassertiefe und Strömungsgeschwindigkeit).

Nachhaltiges Flussmanagement beinhaltet die Herausforderung, einen aquatischen Lebensraum für Zielfischarten in verschiedenen Lebensphasen zu gestalten. Das Konzept des nutzbaren physischen Lebensraums stellt ein leistungsfähiges Werkzeug dar, um die Bewertung der ökologischen Integrität von Flussmanagement- und Engineering-Maßnahmen zu nutzen. So können beispielsweise durch die Berechnung des nutzbaren Lebensraums vor und nach der Durchführung von Maßnahmen wertvolle Rückschlüsse auf die ökologische Integrität der Sanierungsbemühungen gezogen werden.

Diese Übung demonstriert die Verwendung von 2d hydrodynamischen Modellierungsergebnissen zur algorithmischen Bewertung des nutzbaren Lebensraums basierend auf der Berechnung von geospatially expliziten $cHSI$-Werten.

Verfügbare Daten und Codestruktur¶

Das folgende Flussdiagramm verdeutlicht den bereitgestellten Code und die Daten. Funktionen, Methoden und Dateien, die in dieser Übung erstellt werden, werden in fetter, italischer, YELLOW-ish Schriftart hervorgehoben.

Die bereitgestellte QGIS-Projektdatei visualize_with_QGIS.qgz hilft, Eingaberasterdatensätze und -ergebnisse zu überprüfen.

Zweidimensionale (2d) hydrodynamische Modellierung (Folder: BASEMENT)¶

Diese Übung verwendet (hydraulische) Strömungsgeschwindigkeits- und Wassertiefenraster (GeoTIFFs) mit der Software ETH Zürich*s BASEMENT. Lesen Sie mehr über hydrodynamische Modellierung mit BASEMENT im Kapitel ABSCHNITT. Die hydraulischen Raster wurden mit den Beispieldaten des Flaz River in der Schweiz) des BASEMENT-Entwicklers erstellt (weitere Informationen finden Sie auf der Website). Die Wassertiefe water_depth.tif und Fließgeschwindigkeit flow_velocity.tif rasters sind für diese Übung im Ordner /basement/ vorgesehen.

Habitat Suitability Index HSI Curves (folder: habitat)¶

Der /habitat/-Ordner im Übungsarchiv enthält $HSI$-Kurven in Form eines xlsx-Arbeitsbuchs (trout.xlsx) und in Form einer JSON-Datei (trout.json). Beide Dateien enthalten die gleichen Daten für Regenbogenforelle eines hypothetischen Cobble-Bett-Flusses und diese Übung verwendet nur die JSON-Datei (das Arbeitsbuch dient nur der visuellen Überprüfung).

Code¶

HSI Raster erstellen und kombinieren¶

Beenden Sie die __init__ Methode der RasterKlasse (raster.py)¶

The raster.py script imports the functions and libraries loaded in the fun.py script, and therefore, also the config.py script. For this reason, the NumP and Pandas libraries are already available (as np and pd, respectively), and the geo_utils package is already imported as geo (import geo_utils as geo in config.py).

Die Raster-Klasse lädt jeden GeoTIFF-Dateinamen als Geo-Referenzed Array-Objekt, das mit mathematischen Operatoren verwendet werden kann. Zunächst ergänzen wir die __init__ Methode durch eine Raster.name (Auszug aus dem file_name-Argument) sowie Georeferenzen und Array-Datensätze:

    # __init__(...) of Raster class in raster.py
        self.name = file_name.split("/")[-1].split("\\")[-1].split(".")[0]

Wenn die bereitgestellte file_name nicht existiert, erstellt die __init__-Methode einen neuen Raster mit der file_name (dieses Verhalten wird bereits in der if not os.path.exists(file_name)-Anweisung umgesetzt. Anschließend laden Sie die osgeo.gdal.dataset, die np.array und die geo_transformation des Rasters ein. Zu diesem Zweck nutzen Sie das raster2array function aus diesem eBook, das auch im geo utils (geo)-Paket der Übung implementiert ist:

    # __init__(...) of Raster class in raster.py
        self.dataset, self.array, self.geo_transformation = geo.raster2array(file_name, band_number=band)

Um die EPSG number (Authority code) eines Rasters zu identifizieren, holen Sie das räumliche Referenzsystem (SRS) des Rasters ab. Auch zu diesem Zweck haben wir bereits eine Funktion im Vortrag mit der get_srs von der theory section on reprojection entwickelt. Laden Sie die SRS und die EPSG-Nummer mit der get srs-Funktion mit den folgenden zwei Codezeilen in der __init__ Methode:

    # __init__(...) of Raster class in raster.py
        self.srs = geo.get_srs(self.dataset)
        self.epsg = int(self.srs.GetAuthorityCode(None)

Die __init__-Methode der Raster-Klasse ist komplett.

Vollständige magische Methoden der Raster Klasse (raster.py)¶

Um mathematische Operationen zwischen mehreren Instanzen der Raster-Klasse zu ermöglichen, implementieren Sie Überladen und magische Methoden, die der Klasse sagen, was zu tun ist, wenn zwei RasterInstanzen zum Beispiel hinzugefügt werden (+-Zeichen), multipliziert (*-Zeichen), oder subtrahiert (--Zeichen). So können z.B. die magischen Methoden __truediv__ (für die Nutzung des /-Operators), __mul__ (für die Nutzung des *-Operators) und __pow__ (für die Nutzung des **-Operators) die Nutzung von Raster-Instanzen wie folgt ermöglichen:

# example for Raster instances, when operators are defined through magic methods

# load GeoTIFF rasters from file directory
velocity = Raster("/usr/geodata/u.tif")
depth = Raster("/usr/geodata/h.tif")

# calculate the Froude number using operators defined with magic methods
Froude = velocity / (depth * 9.81) ** 0.5

# save the new raster
Froude.save("/usr/geodata/froude.tif")

Die Raster-Klassenvorlage enthält bereits eine beispielhafte magische Methode, um Division zu ermöglichen (__truediv__):

    # Raster class in raster.py
    def __truediv__(self, constant_or_raster):
        try:
            self.array = np.divide(self.array, constant_or_raster.array)
        except AttributeError:
            self.array /= constant_or_raster
        return self._make_raster("div")

Hier ist, was die __truediv__ Methode tut:

• Das Eingabeargument constant_or_raster kann eine weitere RasterInstanz sein, die ein array Attribut oder eine numerische Konstante hat (z.B. 9.81).

• Die Methode versucht, das Array-Attribut von constant_or_raster anzurufen.

  • Wenn constant_or_raster ein Rasterobjekt ist, dann ist contant_or_raster.array erfolgreich. In diesem Fall wird self.array mit der elementaren Teilung des Arrays von contant_or_raster.array überschrieben. Die Element-weise Division baut auf NumPs integrierte Funktion np.divide, die eine rechnerisch effiziente Umschlingung von C/C++-Code ist (viel schneller als eine Python-Schleife über Array-Elemente).

  • Wenn constant_or_raster ein numerischer Wert ist, dann ruft contant_or_raster.array eine AttributeError an und die __truediv__-Methode fällt in die except AttributeError-Anweisung, wo self.array einfach durch constant_or_raster geteilt wird.

• Die Methode gibt das Ergebnis der Pseudo-Privatmethode self._make_raster("div") ([Recall PEP 8 Code Stil und Konventionen, die einer neuen RasterInstanz der eigentlichen RasterInstanz entspricht, die von constant_or_raster geteilt wird. Die neue RasterInstanz ist eine temporäre GeoTIFF-Datei im cache-Ordner (Recall the cache function). So sieht die pseudo-private Methode _make_raster(self, file_marker) aus:

Diese Funktion:

• Verwenden Sie das string-Argument file_marker, um es dem Dateinamen der Basis GeoTIFF zusammen mit einem zufälligen, vier Zeichen langen string hinzuzufügen (Recall the create_random_string). file_marker ist einzigartig für jeden implementierten Betreiber. Für die Methode __truediv__ verwenden Sie file_marker="div". So wird der temporäre GeoTIFF-Dateiname als cache_folder + self.name + f_ending (z.B. "C:\Excercise-geco\__cache__\velocity__divhjev__.tif") definiert.

• Setzt die Erstellen eines Rasters (Array to Raster)-Funktion von geo_utils um den temporären GeoTIFF an den __cache__-Ordner mit dem räumlichen Referenzsystem des ursprünglichen Rasters zu schreiben.

• returns eine neue RasterInstanz der temporären GeoTIFF-Datei.

from __future__ import annotations

    def __truediv__(self, other: Raster) -> Raster:
        f_ending = "__div%s__.tif" % create_random_string(4)
        return Raster(file_name=cache_folder + self.name + f_ending,
                      raster_array=np.divide(self.array, other.array),
                      epsg=self.epsg,
                      geo_info=self.geo_transformation)

Bei der Anwendung der _make_raster-Methode fügen Sie die folgenden magischen Methoden an die Raster-Klasse (Funktions-Platzhalter sind bereits in der raster.py-Vorlage vorhanden):

• __add__ (+ Operator):

        try:
            self.array += constant_or_raster.array
        except AttributeError:
            self.array += constant_or_raster
        return self._make_raster("add")

• __mul__ (* Operator):

        try:
            self.array = np.multiply(self.array, constant_or_raster.array)
        except AttributeError:
            self.array *= constant_or_raster
        return self._make_raster("mul")

• __pow__ (** Operator):

        try:
            self.array = np.power(self.array, constant_or_raster.array)
        except AttributeError:
            self.array **= constant_or_raster
        return self._make_raster("pow")

• __sub__ (- Operator):

        try:
            self.array -= constant_or_raster.array
        except AttributeError:
            self.array -= constant_or_raster
        return self._make_raster("sub")

Der letzte Punkt, der in der Raster-Klasse fertig gestellt wird, ist die integrierte save-Methode, die ein file_name(string)-Argument erhält, das das Verzeichnis definiert und den Namen der Raster-Instanz als Speicher bezeichnet:

        save_status = geo.create_raster(file_name, self.array, epsg=self.epsg, nan_val=0.0, geo_info=self.geo_transformation)
        return save_status

Warum brauchen wir die Variable save_status? Zuerst heißt es, wenn die Speicherung des Rasters erfolgreich war (save_status=0), und zweitens könnten diese Informationen verwendet werden, um den Raster aus dem __cache__-Ordner zu löschen und den Speicher zu spülen (für die Beschleunigung des Codes kostenlos).

Schreiben Sie HSI und cHSI Raster Creation Script¶

Das bereitgestellte create_hsi_rasters.pyScript enthält bereits die erforderlichen Paketimporte, eine if __name__ == '__main__'Stand-alone-Anweisung sowie die Leeremain,get_hsi_curve, get_hsi_raster und combine_hsi_rastersFunktionen:

Die if __name__ == '__main__'-Anweisung enthält einen Zeitzähler (perf_counter), der anfordert, wie lange das Skript läuft (typisch zwischen 3 bis 6 Sekunden). Stellen Sie sicher, dass

• Die parameters-Liste enthält "velocity" und "depth" (gemäß par_dict im config.pyScript),

• die Dateipfade korrekt definiert werden, und

• Eine Lebensphase ist definiert (d.h. entweder "fry", "juvenile", "adult", oder "spawning" gemäß /habitat/fish.xlsx Arbeitsmappe).

Die folgenden Absätze zeigen Schritt für Schritt, wie Sie die $HSI$Kurven aus der JSON-Datei (get_hsi_curve) laden können, wenden Sie sich an die flow_velocity und water_depth rasters (get_hsi_raster) und kombinieren Sie die resultierenden $HSI$ rasters in $cHSI$ rasters (combine_hsi_rasters).

The get_hsi_curve function will load the $HSI$ curve from the JSON file (/habitat/trout.json) in a dictionary for the two parameters "velocity" and "depth". Thus, the goal is to create a curve_data dictionary that contains one Pandas DataFrame object for all parameters (i.e., velocity and depth). For example, curve_data["velocity"]["u"] will be a Pandas Series of velocity entries (in m/s) that corresponds to curve_data["velocity"]["HSI"], which is a Pandas Series of $HSI$ values. Similarly, curve_data["depth"]["h"] is a Pandas Series of depth entries (in meters) that corresponds to curve_data["depth"]["HSI"], which is a Pandas Series of $HSI$ values (corresponds to the curves shown in the HSI graphs above). To extract the desired information from the JSON file, get_hsi_curve takes three arguments (json_file, life_stage, and parameters) in order to:

• Erhalten Sie die in der JSON-Datei gespeicherten Informationen mit der read_json-Funktion (see above).

• Instantiate a void curve_data Wörterbuch that will contain the Pandas DataFrames for "velocity" and "depth".

• Führen Sie eine Schleife über die (zwei) Parameter ("velocity" und "depth") aus, in denen es:

  • Creates a void par_pairs Liste for storing pairs of parameter (par) - $HSI$ values as nested lists.

  • Iterates through the length of provided curve data, where valid data pairs (e.g., [u_value, HSI_value]) are appended to the par_pairs Liste. This iteration is what actually creates the nested Liste.

  • Converts the final par_pairs list to a Pandas DataFrame that it adds to the curve_data Wörterbuch.

• return the curve_data Wörterbuch with its Pandas DataFrames.

# create_hsi_rasters.py
def get_hsi_curve(json_file, life_stage, parameters):
    # read the JSON file with fun.read_json
    file_info = read_json(json_file)
    # instantiate output dictionary
    curve_data = {}
    # iterate through parameter list (e.g., ["velocity", "depth"])
    for par in parameters:
        # create a void list to store pairs of parameter-HSI values as nested lists
        par_pairs = []
        # iterate through the length of parameter-HSI curves in the JSON file
        for i in range(0, file_info[par][life_stage].__len__():
            # if the parameter is not empty (i.e., __len__ > 0), append the parameter-HSI (e.g., [u_value, HSI_value]) pair as nested list
            if str(file_info[par][life_stage][i]["HSI"]).__len__() > 0:
                try:
                    # only append data pairs if both parameter and HSI are numeric (floats)
                    par_pairs.append([float(file_info[par][life_stage][i][par_dict[par]]),
                                      float(file_info[par][life_stage][i]["HSI"])])
                except ValueError:
                    logging.warning("Invalid HSI curve entry for {0} in parameter {1}.".format(life_stage, par)
        # add the nested parameter pair list as pandas DataFrame to the curve_data dictionary
        curve_data.update({par: pd.DataFrame(par_pairs, columns=[par_dict[par], "HSI"])})
    return curve_data

In der main-Funktion rufen Sie get_hsi_curves an, um die $HSI$-Kurven als Wörterbuch zu erhalten. Darüber hinaus setzen Sie die cache und die log_actions-Wrapper (siehe Beschreibungen von provided functions) für die main-Funktion ein:

# create_hsi_rasters.py

...

@log_actions
@cache
def main():
    # get HSI curves as pandas DataFrames nested in a dictionary
    hsi_curve = get_hsi_curve(fish_file, life_stage=life_stage, parameters=parameters)

...

With the provided HSIRaster (raster_hsi.py) class, the $HSI$ rasters can be conveniently created in the get_hsi_raster function. Before using the HSIRaster class, make sure to understand how it works. The HSIRaster class inherits from the Raster class and initiates its parent class in its __init__ method through Raster.__init__(self, file_name=file_name, band=band, raster_array=raster_array, geo_info=geo_info). Then, the class calls its make_hsi method, which takes an $HSI$ curve (nested Liste) of two equal Liste pairs (Liste of parameters and Liste of $HSI$ values) as argument. The make_hsi method:

• Extracts parameter values (e.g., depth or velocity) from the first element of the nested hsi_curves Liste, and $HSI$ values from the second element of the nested hsi_curves Liste.

• Uses NumPs built-in np.nditer function, which iterates through NumP arrays with high computational efficiency (read more about nditer).

  • The nditer loop passes the par_values as x_values Liste argument and the hsi_values as y_values Liste arguments to the interpolate_from_list function (recall the function descriptions above).

  • The array values (i.e., flow velocity or water depth) correspond to the xi_values Liste argument of the interpolate_from_list function.

  • The interpolate_from_list function then identifies for each element of the xi_values Liste the closest elements ($x_{i}$ values) in the x_values Liste and the corresponding positions in the y_values Liste.

  • Die interpolate_from_list-Funktion übergibt die ermittelten Werte an die interpolate_y-Funktion, die dann den entsprechenden yi-Wert linear interpoliert (d.h. einen $HSI$-Wert).

  • Die Strömungsgeschwindigkeit bzw. die Wassertiefe in self.array sind also reihenweise (row-by-row) ersetzt durch $HSI$-Werte.

• returns a RasterInstanz unter Verwendung der pseudo-privaten _make_raster-Methode (recall its contents).

# raster_hsi.py
from raster import *

class HSIRaster(Raster):
    def __init__(self, file_name, hsi_curve, band=1, raster_array=None, geo_info=False):
        Raster.__init__(self, file_name=file_name, band=band, raster_array=raster_array, geo_info=geo_info)
        self.make_hsi(hsi_curve)

    def make_hsi(self, hsi_curve):
        par_values = hsi_curve[0]
        hsi_values = hsi_curve[1]
        try:
            with np.nditer(self.array, flags=["external_loop"], op_flags=["readwrite"]) as it:
                for x in it:
                    x[...] = interpolate_from_list(par_values, hsi_values, x)
        except AttributeError:
            print("WARNING: np.array is one-dimensional.")
        return self._make_raster("hsi")

Ändern Sie die get_hsi_rasters-Funktion, um direkt ein HSIRasterObjekt zurückzugeben:

# create_hsi_rasters.py

...

def get_hsi_raster(tif_dir, hsi_curve):
    return HSIRaster(tif_dir, hsi_curve)
...

The get_hsi_raster function requires two arguments, which it must receive from the main function. For this reason, iterate over the parameters Liste in the main function and extract the corresponding raster directories from the tifs Wörterbuch (recall the variable definition in the standalone statement). In addition, save the Raster objects returned by the get_hsi_raster function in another Wörterbuch (eco_rasters) to combine them in the next step into a $cHSI$ raster.

# create_hsi_rasters.py

...

@log_actions
@cache
def main():
    # get HSI curves as pandas DataFrames nested in a dictionary
    hsi_curve = get_hsi_curve(fish_file, life_stage=life_stage, parameters=parameters)

        # create HSI rasters for all parameters considered and store the Raster objects in a dictionary
    eco_rasters = {}
    for par in parameters:
        hsi_par_curve = [list(hsi_curve[par][par_dict[par]]),
                         list(hsi_curve[par]["HSI"])]
        eco_rasters.update({par: get_hsi_raster(tif_dir=tifs[par], hsi_curve=hsi_par_curve)})
        eco_rasters[par].save(hsi_output_dir + "hsi_%s.tif" % par)
...

Selbstverständlich kann man auch die Parameter Liste direkt in der get_hsi_raster-Funktion einschleifen.

Als nächstes kommen wir zu dem Grund, warum wir magische Methoden für die Raster-Klasse definieren mussten: Kombinieren Sie die $HSI$ rasters mit beiden oben vorgestellten Kombinationsformeln (Recall the product and geometric meanFormel), wobei "geometric_mean"Standard verwendet werden sollte. Die combine_hsi_rasters-Funktion akzeptiert zwei Argumente (eine Liste von Raster Objekten entsprechend $HSI$ rasters und die method, um als string zu verwenden).

If the method corresponds to the default value "geometric_mean", then the power to be applied to the product of the Raster Liste is calculated from the nth root, where n corresponds to the number of Raster objects in the raster_list. Otherwise (e.g., method="product"), the power is exactly 1.0.

The combine_hsi_rasters function initially creates an empty $cHSI$ Raster in the cache_folder, with each cell having the value 1.0 (filled through np.ones). In a loop over the Raster elements of the raster_list, the function multiplies each $HSI$ raster with the $cHSI$ raster.

Schließlich gibt die Funktion das Produkt aller $HSI$ rasters an die Macht des zuvor ermittelten power-Wertes zurück.

# create_hsi_rasters.py
def combine_hsi_rasters(raster_list, method="geometric_mean"):
    if method is "geometric_mean":
        power = 1.0 / float(raster_list.__len__()
    else:
        # supposedly method is "product"
        power = 1.0

    chsi_raster = Raster(cache_folder + "chsi_start.tif",
                         raster_array=np.ones(raster_list[0].array.shape),
                         epsg=raster_list[0].epsg,
                         geo_info=raster_list[0].geo_transformation)
    for ras in raster_list:
        chsi_raster = chsi_raster * ras

    return chsi_raster ** power

Um das create_hsi_rasters.py-Skript zu beenden, setzen Sie den Anruf an die combine_hsi_rasters-Funktion in der main-Funktion um und speichern Sie das Ergebnis als $cHSI$ GeoTIFF raster im /habitat/Ordner:

# create_hsi_rasters.py

...

@log_actions
@cache
def main():
    ...

    for par in parameters:
        hsi_par_curve = [list(hsi_curve[par][par_dict[par]]),
                         list(hsi_curve[par]["HSI"])]
        eco_rasters.update({par: get_hsi_raster(tif_dir=tifs[par], hsi_curve=hsi_par_curve)})
        eco_rasters[par].save(hsi_output_dir + "hsi_%s.tif" % par)

    # get and save chsi raster
    chsi_raster = combine_hsi_rasters(raster_list=list(eco_rasters.values(),
                                      method="geometric_mean")
    chsi_raster.save(hsi_output_dir + "chsi.tif")
...

Führen Sie den HSI-Code aus¶

Ein erfolgreicher Lauf des Skripts create_hsi_rasters.py soll so aussehen (in PyCharm):

Plotted in QGIS, the $cHSI$ GeoTIFF raster should look like this:

Ergebnisse¶

Die Präsentation des $cHSI$ raster zeigt, dass bevorzugte Lebensräume für juvenile Forellen nur in der Nähe der Banken existieren. Auch sind numerische Artefakte des dreieckigen Netzes von ABSCHNITT sichtbar. Daher stellt sich die Frage, ob die berechneten Strömungsgeschwindigkeiten und Wassertiefen, und damit auch die $cHSI$-Werte, nahe den Banken, als repräsentativ angesehen werden können.

Berechnen Usable Habitat Area UHA¶

Schreiben Sie den Code¶

Die $cHSI$ rasters ermöglichen die Berechnung der verfügbaren nutzbaren Lebensraumfläche. Der vorherige Abschnitt umfasste Beispiele mit der Fischart trout und der juvenile-Lebensstufe, für die wir hier den nutzbaren Lebensraum $UHA$ (in m$^2$) unter Verwendung eines $cHSI$-Schwellwerts (anstelle des Pixelflächengewichtungsansatzes) bestimmen werden. So folgen wir dem threshold formula described above, mit einem Schwellenwert von $cHSI_{crit} = 0.4$. Jedes Pixel, das einen $cHSI$-Wert von 0,4 oder mehr als nutzbare Lebensraumfläche hat.

Aus technischer Sicht geht es in diesem Teil der Übung darum, einen Raster in eine Polygon-Formdatei umzuwandeln sowie auf die Beitragstabelle der Formdatei zuzugreifen und zu modifizieren.

Ähnlich wie bei der Erstellung des $cHSI$ raster gibt es ein Vorlagenskript für diesen Teil der Übung, genannt calculate_habitat_area.py, das Paket- und Modulimporte enthält, eine if __name__ == '__main__'Stand-alone-Anweisung sowie die Leere main und calculate_habitat_area Funktionen. (uha-template)= Das Template-Skript sieht so aus:

# this is calculate_habitat_area.py (template)
from fun import *
from raster import Raster


def calculate_habitat_area(layer, epsg):
    pass


def main():
    pass


if __name__ == '__main__':
    chsi_raster_name = os.path.abspath("") + "\\habitat\\chsi.tif"
    chsi_threshold = 0.4

    main()

In der if __name__ == '__main__'-Anweisung stellen Sie sicher, dass die globale Variable chsi_raster_name dem Verzeichnis des $cHSI$ raster entspricht, das im vorherigen Abschnitt erstellt wurde. Die andere globale Variable (chsi_threshold) entspricht dem $cHSI_{crit}$-Wert von 0,4, den wir mit der threshold formula verwenden werden.

In der Funktion main starten Sie mit dem Laden des $cHSI$ raster (chsi_raster) als Raster object. Dann nutzen Sie das NumP Array des $cHSI$ raster und vergleichen Sie es mit chsi_threshold unterNumPs integriertem greater equalfunktion. np.greater_equal nimmt ein Array als erstes Argument und ein zweites Argument an, das die Bedingung ist, die eine numerische Variable oder ein anderes NumPArray sein kann. Dann prüft np.greater_equal, ob die Elemente des ersten Arrays größer oder gleich dem zweiten Argument sind. Für den Fall, dass das zweite Argument ein Array ist, ist dies ein elementarer $\geq$-Vergleich. Das Ergebnis von np.greater_equal ist ein Boolen-Array (True, bei dem die mehr oder minderwertige Bedingung erfüllt ist und False andernfalls). Um jedoch aus dem Ergebnis von np.greater_equal ein osgeo.gdal.Dataset-Objekt zu erstellen, benötigen wir ein numerisches Array. Aus diesem Grund vervielfältigen Sie das Ergebnis von np.greater_equal by 1.0 und geben Sie es als neue NumP Array of Zeros (False) und eines (True) an eine Variable namens habitat_pixels an (siehe Codeblock unten).

With the habitat_pixels array and the georeference of chsi_raster, create a new integer GeoTIFF raster with the create_raster function (also available in flusstools.geotools); here, use geo.create_raster. In the following code block the new raster is saved in the /habitat/ folder of the exercise as habitat-pixels.tif.

# calculate_habitat_area.py
...

def main():
    # open the chsi raster
    chsi_raster = Raster(chsi_ras_name)
    # extract pixels where the physical habitat quality is higher than the user threshold value
    habitat_pixels = np.greater_equal(chsi_raster.array, chsi_threshold_value) * 1
    # write the habitat pixels to a binary array (0 -> no habitat, 1 -> usable habitat)
    geo.create_raster(os.path.abspath("") + "\\habitat\\habitat-pixels.tif",
                      raster_array=habitat_pixels,
                      epsg=chsi_raster.epsg,
                      geo_info=chsi_raster.geo_transformation)
...

Im nächsten Schritt wandeln Sie den Habitat-Pixel-Raster in eine Polygon-Formdatei um und speichern Sie ihn im /habitat/-Ordner als habitat-area.shp. Die Umwandlung eines Rasters in eine Polygon-Formdatei erfordert, dass der Raster nur integer-Werte enthält, was im Habitat-Pixel-Raster der Fall ist (nur Nullen und eins - Rückruf Raster nach Polygon). Verwenden Sie die raster2polygon-Funktion im Ordner geo utils (Paket), um die neue Polygon-Formdatei zu erstellen, habitat-pixels.tif als raster_file_name zu konvertieren, und /habitat/habitat-area.shp als Ausgabedateiname. Die geo.raster2polygon-Funktion gibt ein osgeo.ogr.DataSource-Objekt zurück und wir können ihre Schicht einschließlich der Informationen des EPSG-Behördecodes (vonchsi_raster) direkt an die nicht-yet-writtencalculate_habitat_area()-Funktion weitergeben:

# calculate_habitat_area.py
...

def main():
    ... (create habitat pixels raster)

    # convert the raster with usable pixels to polygon (must be an integer raster!)
    tar_shp_file_name = os.path.abspath("") + "\\habitat\\habitat-area.shp"
    habitat_polygons = geo.raster2polygon(os.path.abspath("") + "\\habitat\\habitat-pixels.tif",
                                          tar_shp_file_name)

    # calculate the habitat area (will be written to the attribute table)
    calculate_habitat_area(habitat_polygons.GetLayer(), chsi_raster.epsg)
...

For the calculate_habitat_area() function to produce what its name promises, we need to populate this function as well. For this purpose, use the epsg integer argument to identify the unit system of the shapefile.

# calculate_habitat_area.py
...

def calculate_habitat_area(layer, epsg):
    # retrieve units
    srs = geo.osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    area_unit = "square %s" % str(srs.GetLinearUnitsName()
...

To determine the habitat area, the area of each polygon must be calculated. For this purpose, add a new field to the layer in the Attribute Table, name it "area", and assign a geo.ogr.OFTReal (numeric) data type (recall how to create a field an data types). Then, create a void Liste called poly_size, in which we will write the area of all polygons that have a field value of 1. To access the individual polygons (features) of the layer, iterate through all features using a for loop, which:

• Extrahiert das Polygon jedes feature unter polygon = feature.GetGeometryRef()

• Appends the polygon’s area size to the poly_size Liste if the field "value" of the polygon (at position 0: feature.GetField(0)) is 1 (True).

• Schreibt die Flächengröße des Polygons an die Beitragstabelle mit feature.SetField("area", polygon.GetArea().

• Speichert die Änderungen (berechnete Fläche) an die formfilelayer mit layer.SetFeature(feature).

The last information needed after the for loop is the total area of the "value"=1 polygons, which we get by writing the sum of the poly_size Liste to the console. Therefore, the second and last part of the calculate_habitat_area function looks like this:

# calculate_habitat_area.py
...
def calculate_habitat_area(layer, epsg):

    ... (extract unit system information)

    # add area field
    layer.CreateField(geo.ogr.FieldDefn("area", geo.ogr.OFTReal)

    # create list to store polygon sizes
    poly_size = []

    # iterate through geometries (polygon features) of the layer
    for feature in layer:
        # retrieve polygon geometry
        polygon = feature.GetGeometryRef()
        # add polygon size if field "value" is one (determined by chsi_treshold)
        if int(feature.GetField(0):
            poly_size.append(polygon.GetArea()
        # write area to area field
        feature.SetField("area", polygon.GetArea()
        # add the feature modifications to the layer
        layer.SetFeature(feature)

    # calculate and print habitat area
    print("The total habitat area is {0} {1}.".format(str(sum(poly_size), area_unit)

...

Führen Sie den nutzbaren Habitat-Bereich Berechnungscode¶

Ein erfolgreicher Lauf des Skripts calculate_habitat_area.py soll so aussehen (in PyCharm):

Die in QGIS unterteilte habitat-area Formdatei sieht so aus (verwenden Sie Kategorisierte Symbologie):

Ergebnisse¶

Die $UHA$ des analysierten Flussabschnitts stellt einen sehr geringen Anteil der gesamten benetzten Fläche dar, die als ökologisch schlechten Zustand des Flusses interpretiert werden kann. Ein Blick auf eine Karte und die Simulationsdateien des Beispiels Flaz von BASEMENT deutet jedoch darauf hin, dass bei einer Entladung von 50 m$^3$/s eine Flutsituation angenommen werden kann. Wie bei Überschwemmungen gibt es in der Regel höhere Strömungsgeschwindigkeiten, die von jugendlichen Fischen abweichen, ist die kleine nutzbare Lebensraumfläche schließlich nicht überraschend.