Raster (Grid) Handling - Hydro-Informatics

This section introduces geospatial analysis of raster (gridded) data with gdal and rasterstats. For interactive reading and executing code blocks and find geo02-raster.ipynb, or install Python and JupyterLab locally.

Watch this section and the Python tutorials in video formats

Watch this section as a video on the @Hydro-Morphodynamics channel on YouTube.

Lastraster¶

Offene vorhandene Rasterdaten¶

Rasterdaten können im Python-Code als Instanz von gdal.Open("FILENAME") geöffnet werden. Der folgende Code-Block bietet eine Funktion, um alle mit dem file_name-Eingabe-Argument angegebenen Raster zu öffnen. Eines der wichtigsten Elemente beim Umgang mit Rasterdaten ist das Rasterband, das eine ähnliche Datenträgerrolle wie GetLayer in Shapefile Handling übernimmt. Um auf die Raster-Band zuzugreifen, ist die unten gezeigte open_rasterfunktion:

Ermöglicht Fehler- und Warn-Feedback mit gdal.UseExceptions() (dieser Schritt ist entscheidend bei der Verwendung von gdal).
Eröffnet die bereitgestellten raster file_name, die von try -except-Anweisungen umfasst werden, um zu informieren, ob und warum ein potenzieller Fehler beim Öffnen des Rasters aufgetreten ist.
Öffnet die raster Bandnummer, die im optionalen band_numberkeyword-Argument mit raster_band = raster.GetRasterBand(band_number) angegeben ist (der Standardwert ist 1).
Gibt die Raster- und Rasterbandobjekte zurück.

from osgeo import gdal
import numpy as np


def open_raster(file_name, band_number=1):
    """
    Open a raster file and access its bands
    :param file_name: STR of a raster file directory and name
    :param band_number: INT of the raster band number to open (default: 1)
    :output: osgeo.gdal.Dataset, osgeo.gdal.Band objects
    """
    gdal.UseExceptions()
    # open raster file or return None if not accessible
    try:
        raster = gdal.Open(file_name)
    except RuntimeError as e:
        print("ERROR: Cannot open raster.")
        print(e)
        return None
    # open raster band or return None if corrupted
    try:
        raster_band = raster.GetRasterBand(band_number)
    except RuntimeError as e:
        print("ERROR: Cannot access raster band.")
        print(e)
        return None
    return raster, raster_band

Um die open_raster-Funktion zu nutzen, rufen Sie sie mit einem Dateinamen wie im folgenden Codeblock mit dem h001000.tif raster aus der River Architect Sample data. Das Skript schließt den Raster sofort wieder durch Überschreiben der Variablen mit einem None-Instance, um zu vermeiden, dass die GeoTIFF-Datei nachher gesperrt wird.

import os
file_name = r"" + os.getcwd() + "/geodata/rasters/h001000.tif"
src, depth = open_raster(file_name)
print(src)
print(depth)
depth = None

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x77015c0a7b70> >
<osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x77015c0f39b0> >

Rasterband Statistiken und Toolbox Scripts¶

Sobald der Raster und seine Band mit der obigen open_rasterfunktion geladen sind, können wir auf statistische Informationen (z.B. das Minimum oder das Maximum) zugreifen, den no-data-Wert (d.h. einen vordefinierten Wert, der Pixeln ohne Wert zugewiesen wird) oder den Typ der verwendeten Einheiten identifizieren.

Python-Skripte zur Verarbeitung von Geospatialdaten können auch als Plugins in GIS-Desktopanwendungen eingebettet werden (z.B. als Plugins in QGIS oder Toolbox in ArcGIS Pro). Um ein Python-Skript in einer GIS-Desktop-Anwendung auszuführen, sollte es als eigenständiges Skript geschrieben werden, das Eingabeargumente empfangen kann. Der interessierte Leser kann mehr über die Implementierung von Plugins in QGIS in der QGIS docs.

Hier schreiben wir nur den nächsten Codeblock, damit er als eigenständiges Skript in einer Konsole/Endanwendung ausgeführt werden kann (Recall the instructions for writing standalone script).

import sys
from osgeo import gdal

# make sure to use exceptions
gdal.UseExceptions()

def how2use():
    # provide usage instructions for the script
    print("""
    $ raster_band_info.py [ band number ] input-raster
    """)
    # exit program if wrong input arguments provided
    sys.exit(1)
    

def get_color_bands(raster_band):
    """
    :param raster_band: osgeo.gdal.Band object
    :output: list of color bands used in raster_band
    """ 
    
    # get ColorTable and return False if None
    color_table = raster_band.GetColorTable()
    if color_table is None:
        print("Band has no ColorTable.")
        return None
    else:
        print("Found %i color definitions." % int(color_table.GetCount()))

    # iterate through color_table and append objects found to colors_bands list
    color_bands = []
    for c in range(0, color_table.GetCount() ):
        entry = color_table.GetColorEntry(c)
        if not entry:
            continue
        color_bands.append(str(color_table.GetColorEntryAsRGB(c, entry)))
    return color_bands

def main(band_number, input_file):
    src, band = open_raster(input_file)
    print("Band minimum: ", band.GetMinimum())
    print("Band maximum: ", band.GetMaximum())
    print("No-data value: ", band.GetNoDataValue())
    print("Band unit type: ", band.GetUnitType())    

    try:
        print(", ".join(get_color_bands(band)))
    except TypeError:
        print("ColorTable: None")

if __name__ == '__main__':
    # make standalone
    if len( sys.argv ) < 3:
        print("""
        ERROR: Provide two arguments:
        1) the band number (int) and 2) input raster directory (str)
        """)
        how2use()

    main(int(sys.argv[1]), str(sys.argv[2]))

Um dieses Skript auszuführen, speichern Sie es als raster band info.py (z.B. in C:\temp in Windows oder ~/temp/ in Linux) und navigieren Sie mit dem Befehl cd. Dann führen Sie das Skript, um Statistiken über die Wassertiefe raster h001000.tif mit (in Windows):

C:\temp\ python raster_band_info.py 1 "C:\temp\geodata\rasters\h001000.tif"

Auf Linux und mit der flussenv activated z.B.:

python raster_band_info.py 1 "~/temp/geodata/rasters/h001000.tif"

Band minimum:  0.0
Band maximum:  7.0613012313843
No-data value:  -3.4028234663852886e+38
Band unit type:
Band has no ColorTable.
ColorTable: None

Erstellen und Speichern eines Rasters (von Array)¶

Rastertreiber¶

Ebenso wie bei Formfile-Dateien muss der entsprechende gdal-Treiber (analog an ogr-Treiber) geladen werden, um einen Raster zu speichern. Um eine vollständige Liste von gdal raster Treibern laufen:

driver_list = [str(gdal.GetDriver(i).GetDescription()) for i in range(gdal.GetDriverCount())]
driver_list.sort()
print(", ".join(driver_list[:]))

AAIGrid, ACE2, ADRG, AIG, AIVector, AVCBin, AVCE00, AirSAR, AmigoCloud, BIGGIF, BMP, BSB, BT, BYN, CAD, CALS, CEOS, COASP, COG, COSAR, CPG, CSV, CSW, CTG, Carto, DAAS, DERIVED, DGN, DIMAP, DOQ1, DOQ2, DTED, DXF, ECRGTOC, EDIGEO, EEDA, EEDAI, EHdr, EIR, ENVI, ERS, ESAT, ESRI Shapefile, ESRIC, ESRIJSON, Elasticsearch, FAST, FlatGeobuf, GDALG, GFF, GIF, GML, GMLAS, GNMDatabase, GNMFile, GPKG, GPSBabel, GPX, GRASSASCIIGrid, GS7BG, GSAG, GSBG, GSC, GTFS, GTI, GTX, GTiff, GXF, GenBin, GeoJSON, GeoJSONSeq, GeoRSS, HF2, HFA, HTTP, ILWIS, IRIS, ISCE, ISG, ISIS2, ISIS3, Idrisi, Interlis 1, Interlis 2, JAXAPALSAR, JDEM, JML, JPEG, JPEGXL, JSONFG, KML, KMLSUPEROVERLAY, KRO, L1B, LAN, LCP, LIBERTIFF, LIBKML, LOSLAS, LVBAG, Leveller, MAP, MBTiles, MEM, MFF, MFF2, MRF, MSGN, MVT, MapInfo File, MapML, MiraMonRaster, MiraMonVector, NAS, NDF, NGSGEOID, NGW, NITF, NOAA_B, NSIDCbin, NTv2, NUMPY, NWT_GRC, NWT_GRD, OAPIF, ODS, OGCAPI, OGR_GMT, OGR_PDS, OGR_VRT, OSM, OpenFileGDB, PAux, PCIDSK, PCRaster, PDS, PDS4, PGDUMP, PLMOSAIC, PLSCENES, PMTiles, PNG, PNM, PRF, RCM, RIK, RMF, ROI_PAC, RPFTOC, RRASTER, RS2, RST, S57, SAFE, SAGA, SAR_CEOS, SENTINEL2, SIGDEM, SNAP_TIFF, SNODAS, SQLite, SRP, SRTMHGT, STACIT, STACTA, SXF, Selafin, TGA, TIL, TSX, Terragen, TopoJSON, USGSDEM, VDV, VFK, VICAR, VRT, WAsP, WCS, WEBP, WFS, WMS, WMTS, XLSX, XYZ, ZMap, Zarr

Rasterdatentypen¶

Die Ausgaberasterpixel können eines der folgenden Datentypen sein (Quelle: gdal.org/doxygen/):

GDT_Unknown Unbekannter oder nicht näher bezeichneter Typ
GDT_Byte 8 bit unsigned Ganzheit
GDT_UInt16 16 bit unsigned Ganzzahl
GDT_Int16 16 bit signiert ganze
GDT_UInt32 32 bit unsigned Ganzzahl
GDT_Int32 32 bit signiert ganze
GDT_Float32 32 bit schwimmender Punkt
GDT_Float64
GDT_CInt16
GDT_CInt32
GDT_CFloat32 Complex Float32
GDT_CFloat64 Complex Float64

Erstellen eines Rasters (Array to Raster)¶

Wenn wir die Grundlagen der Rasterbearbeitung, Datentypen und Python kennen, können wir einen Raster aus einem numerischen Array erstellen. Da ein Raster im Grunde ein Georeferenz-Array ist, ist es bequem, ein numpy array in ein Raster (Band) umzuwandeln. Die Funktionsblöcke enthalten die Umwandlung eines numpy-Arrays in einen GeoTIFF-Raster nach folgendem Workflow:

Sehen Sie sich den GeoTIFF Treiber an (driver = gdal.GetDriverByName('GTiff')).
Retrieve die Array-Größe und (Anzahl der Zeilen rows und Spalten cols).
Erstellen Sie einen neuen GeoTIFF-Raster (new_raster = driver.Create(file_name, cols, rows, 1, eType=rdtype)), wo
- file_name ist das Verzeichnis und der Name der neuen Raster-Datei, die an .tif (z.B. "C:\\temp\\rasters\\new.tif") enden.
- cols, rows repräsentiert die Array-Form und eType ist der Geospatial-Datentyp (siehe oben)
Legen Sie den geographischen Ursprung fest, der im Parameter origin (tuple) gespeichert ist und definieren Sie den pixel_width und pixel_height (Pixeleinheiten, die mit srs - siehe unten definiert sind).
Ersetzen Sie np.nan-Werte im numpy-Array mit nan_value.
Instantiate a band object, setze das NoDataValue an nan_value und schreibe das Array an die band.
Erstellen Sie ein räumliches Referenzsystemobjekt (srs) in Abhängigkeit vom epsgEingabeparameter und exportieren Sie es in WKT-Format.
Lassen Sie den Raster (Grippe aus Cache) frei.

from osgeo import osr


def create_raster(file_name, raster_array, origin=None, epsg=4326, pixel_width=10, pixel_height=10,
                  nan_value=-9999.0, rdtype=gdal.GDT_Float32, geo_info=False):
    """
    Convert a numpy.array to a GeoTIFF raster with the following parameters
    :param file_name: STR of target file name, including directory; must end on ".tif"
    :param raster_array: np.array of values to rasterize
    :param origin: TUPLE of (x, y) origin coordinates
    :param epsg: INT of EPSG:XXXX projection to use - default=4326
    :param pixel_height: INT of pixel height (multiple of unit defined with the EPSG number) - default=10m
    :param pixel_width: INT of pixel width (multiple of unit defined with the EPSG number) - default=10m
    :param nan_value: INT/FLOAT no-data value to be used in the raster (replaces non-numeric and np.nan in array)
                        default=-9999.0
    :param rdtype: gdal.GDALDataType raster data type - default=gdal.GDT_Float32 (32 bit floating point)
    :param geo_info: TUPLE defining a gdal.DataSet.GetGeoTransform object (supersedes origin, pixel_width, pixel_height)
                        default=False
    """
    # check out driver
    driver = gdal.GetDriverByName('GTiff')

    # create raster dataset with number of cols and rows of the input array
    cols = raster_array.shape[1]
    rows = raster_array.shape[0]
    new_raster = driver.Create(file_name, cols, rows, 1, eType=rdtype)    

    # apply geo-origin and pixel dimensions
    if not geo_info:
        origin_x = origin[0]
        origin_y = origin[1]
        new_raster.SetGeoTransform((origin_x, pixel_width, 0, origin_y, 0, pixel_height))
    else:
        new_raster.SetGeoTransform(geo_info)
    
    # replace np.nan values
    raster_array[np.isnan(raster_array)] = nan_value

    # retrieve band number 1
    band = new_raster.GetRasterBand(1)
    band.SetNoDataValue(nan_value)
    band.WriteArray(raster_array)
    band.SetScale(1.0)

    # create projection and assign to raster
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    new_raster.SetProjection(srs.ExportToWkt())

    # release raster band
    band.FlushCache()

Um die Funktion zum Schreiben eines zufälligen Numpy-Arrays zu rufen, können wir nun die create_raster()-Funktion nutzen (auch erhältlich unter geotools.create_raster()):

# set the name of the output GeoTIFF raster
raster_name = r"" + os.getcwd() + "/geodata/rasters/random_unis_dem.tif"
# create a random numpy array (DEM-like values) - can be replaced with any other numpy.array
unis_dem = np.random.rand(300, 300) + 455.0
# overwrite one pixel with np.nan
unis_dem[5, 7] = np.nan
# define a raster origin in EPSG:3857
raster_origin = (1013428.396233, 6231555.006177)
# call create_raster to create a 1-m-resolution raster in EPSG:4326 projection
create_raster(raster_name, unis_dem, raster_origin,  pixel_width=1,  pixel_height=1, epsg=3857)

python create raster file Geotiff — Figure 1:Der neue Raster enthält eine zufällige unis dem Punkthöhe.

Rasterkalkulus (Raster / Band nach Array)¶

Das mit der create_raster()-Funktion beschriebene Verfahren kann umgekehrt verwendet werden, um numpy array von raster bands zu erstellen. Der in ein numpy Array umgewandelte Raster ermöglicht es, algebraische oder andere logische Operationen auf vorhandene Rasterdaten anzuwenden.

Brauchen Sie ein Beispiel? In der RiverArchitect SampleData befinden sich die Einheiten des Wassertiefenrasters h001000.tif in den üblichen Füßen der USA und die Einheiten des Strömungsgeschwindigkeitsrasters u001000.tif sind in Fuß pro Sekunde. Zur Berechnung der Froude number (um die Gravitationskonstanten) für jedes Pixel basierend auf den beiden Rastern (Wassertiefe und Strömungsgeschwindigkeit) ist es jedoch zweckmäßig, beide Raster in m bzw. m/s umzuwandeln. Zu diesem Zweck verfügt der folgende Codeblock über eine weitere wiederverwendbare, benutzerdefinierte Funktion, die einen Raster als Array lädt und NoDataValues mitnp.nan (raster und band mit der oben genannten open_raster-Funktion überschreibt:

def raster2array(file_name, band_number=1):
    """
    :param file_name: STR of target file name, including directory; must end on ".tif"
    :param band_number: INT of the raster band number to open (default: 1)
    :output: (1) ndarray() of the indicated raster band, where no-data values are replaced with np.nan
             (2) the GeoTransformation used in the original raster
    """
    # open the raster and band (see above)
    raster, band = open_raster(file_name, band_number=band_number)
    # read array data from band
    band_array = band.ReadAsArray()
    # overwrite NoDataValues with np.nan
    band_array = np.where(band_array == band.GetNoDataValue(), np.nan, band_array)
    # return the array and GeoTransformation used in the original raster
    return raster, band_array, raster.GetGeoTransform()

Die raster2array()-Funktion ist auch in flusstools enthalten: geotools.raster2array()

Um schließlich eine Froude-Nummer GeoTIFF-Raster zu erstellen, nutzt der folgende Codeblock die raster2array-Funktion zur Umwandlung der Wassertiefe und Strömungsgeschwindigkeit GeoTIFF-Raster in ein numpy-Array, wobei einfache algebraische Berechnungen durchgeführt werden, um die Raster in m bzw. m/s umzuwandeln und die resultierende GeoTIFF-Datei zu speichern. Im Detail beinhaltet der Workflow:

Legen Sie die Eingabe-Raster-Dateinamen mit Verzeichnissen fest (h_file und u_file),
Original-Raster als ndarray mit der raster2array()-Funktion laden und die ursprüngliche GeoTransform-Beschreibung erhalten,
Konvertieren Sie alle Werte von U.S. üblichen Füßen zu S.I. metrischen Einheiten (Recall the feet_to_meter() function from the Python basics), und
Speichern Sie eine neue Kopie des Rasters.

h_file = r"" + os.getcwd() + "/geodata/rasters/h001000.tif"
u_file = r"" + os.getcwd() + "/geodata/rasters/u001000.tif"

# load both rasters as arrays
h_ras, h, h_geo_info = raster2array(h_file)
u_ras, u, u_geo_info = raster2array(u_file)

#convert to metric system
h *= 0.3048
u *= 0.3048

# calculate the Froude number as array and avoid zero-division warning messages
with np.errstate(divide="ignore", invalid="ignore"):
    Froude = u / np.sqrt(h * 9.81)

# create Froude raster from array
create_raster(file_name= r"" + os.path.abspath("") + "/geodata/rasters/Fr1000cfs.tif",
              raster_array=Froude, epsg=6418, geo_info=h_geo_info)

python create raster froude number — Figure 2:Der Froude-Nummernraster berechnet mit Strömungsgeschwindigkeit und Wassertiefe Raster.

Reproject a Raster¶

Die Transformation (und Reprojektion) eines Rasters in ein anderes Koordinatensystem beinhaltet Dreh-, Schalt- und Scherpixel. Wenn einer dieser Operationen übersprungen ist, kann der reprojizierte Raster möglicherweise gedrückt, verdreht oder überall in der Welt platziert werden, aber nicht wo er platziert werden sollte. Der Ansatz zur Reprojektion eines Rasters in ein anderes Koordinatenreferenzsystem impliziert daher die folgenden Schritte:

Retrieve die Quell- und Zielraumreferenzsysteme (z.B. abgeleitet von einem gdal.Dataset oder einem EPSG Authority Code).
Lesen Sie die Geotransformation des Quelldatensatzes (gdal.Dataset.GetGeoTransform()).
Ableiten der Anzahl der Pixel und des Abstandes zwischen Pixeln in dem neuen (reprojektierten) Datensatz.
Instantiate den neuen (reprojektierten) Datensatz.
Projektieren Sie ein Bild des Quelldatensatzes auf den neuen (reprojektierten) Datensatz (gdal.ReprojectImage()).

Das räumliche Referenzsystem kann aus einem Datensatz mit den Erläuterungen im Abschnitt shapefile durch Schreiben einer get_srs()-Funktion abgeleitet werden. Der folgende Codeblock zeigt die get_srs()-Funktion (verwendet die osr-Bibliothek von osgeo /gdal ), die auch in flusstools (geotools.get_srs()) integriert ist.

def get_srs(dataset):
    """
    Get the spatial reference of any gdal.Dataset
    :param dataset: osgeo.gdal.Dataset (raster)
    :output: osr.SpatialReference
    """
    sr = osr.SpatialReference()
    sr.ImportFromWkt(dataset.GetProjection())
    # auto-detect epsg
    auto_detect = sr.AutoIdentifyEPSG()
    if auto_detect != 0:
        sr = sr.FindMatches()[0][0]  # Find matches returns list of tuple of SpatialReferences
        sr.AutoIdentifyEPSG()
    # assign input SpatialReference
    sr.ImportFromEPSG(int(sr.GetAuthorityCode(None)))
    return sr

Mit den open_raster() und get_srs()Funktionen haben wir alle notwendigen Zutaten, um den raster-Reprojektions-Workflow in einer weiteren Funkation unter dem Namen reproject_raster() zu erreichen. Eine zusätzliche Funktion ist, dass sie die korrekte Nutzung von osr.CoordinateTransformation gewährleistet, die sich unter gdal 3.0 im Vergleich zu älteren gdal-Versionen unterschiedlich verhält (weitere Informationen zu OSGeo’s GitHub page]).

def reproject_raster(source_dataset, source_srs, target_srs):
    """
    Reproject a raster dataset (preferably use through reproject function)
    :param source_dataset: osgeo.ogr.DataSource (instantiate with ogr.Open(SHP-FILE))
    :param source_srs: osgeo.osr.SpatialReference (instantiate with get_srs(source_dataset))
    :param target_srs: osgeo.osr.SpatialReference (instantiate with get_srs(DATASET-WITH-TARGET-PROJECTION))
    """
    # READ THE SOURCE GEO TRANSFORMATION (ORIGIN_X, PIXEL_WIDTH, 0, ORIGIN_Y, 0, PIXEL_HEIGHT)
    src_geo_transform = source_dataset.GetGeoTransform()
    
    # DERIVE PIXEL AND RASTER SIZE
    pixel_width = src_geo_transform[1]
    x_size = source_dataset.RasterXSize
    y_size = source_dataset.RasterYSize

    # ensure that TransformPoint (later) uses (x, y) instead of (y, x) with gdal version >= 3.0
    source_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    # get CoordinateTransformation
    coord_trans = osr.CoordinateTransformation(source_srs, target_srs)

    # get boundaries of reprojected (new) dataset
    (org_x, org_y, org_z) = coord_trans.TransformPoint(src_geo_transform[0], src_geo_transform[3])
    (max_x, min_y, new_z) = coord_trans.TransformPoint(src_geo_transform[0] + src_geo_transform[1] * x_size,
                                                       src_geo_transform[3] + src_geo_transform[5] * y_size,)

    # INSTANTIATE NEW (REPROJECTED) IN-MEMORY DATASET AS A FUNCTION OF THE RASTER SIZE
    mem_driver = gdal.GetDriverByName('MEM')
    tar_dataset = mem_driver.Create("",
                                    int((max_x - org_x) / pixel_width),
                                    int((org_y - min_y) / pixel_width),
                                    1, gdal.GDT_Float32)
    # create new GeoTransformation
    new_geo_transformation = (org_x, pixel_width, src_geo_transform[2],
                              org_y, src_geo_transform[4], -pixel_width)

    # assign the new GeoTransformation to the target dataset
    tar_dataset.SetGeoTransform(new_geo_transformation)
    tar_dataset.SetProjection(target_srs.ExportToWkt())

    # PROJECT THE SOURCE RASTER ONTO THE NEW REPROJECTED RASTER
    rep = gdal.ReprojectImage(source_dataset, tar_dataset,
                              source_srs.ExportToWkt(), target_srs.ExportToWkt(),
                              gdal.GRA_Bilinear)
    return tar_dataset

Die Verwendung der reproject_raster()-Funktion in einem Python-Skript erfordert einen Quelldatensatz und einen weiteren (Ausrichtungs-)Datensatz mit dem neuen Koordinatensystem, in das der Quelldatensatz projiziert wird. Das folgende Beispiel zeigt, wie man den oben erstellten Froude-number-Raster in das EPSG=3857-Koordinatensystem neu projiziert, um es in QGIS auf einer Google Satellite-Basiskarte (recall basemaps in QGIS tutorial) anzuzeigen. Als Orientierungsdatensatz verwenden Sie web frame.tif, die auf der Basiskarte Google Satellite erstellt wurde.

Mit der get_srs()-Funktion, die die Rasterprojektion und das räumliche Referenzsystem automatisch erkennt, können wir die create_raster()-Funktion nutzen, um die oben erstellte Fr1000cfs.tif raster (z.B. an epsg=4326) neu zu projizieren.

# load original and orientation rasters
source_file_name = r"" + os.path.abspath("") + "/geodata/rasters/Fr1000cfs.tif"
orientation_file_name = r"" + os.path.abspath("") + "/geodata/rasters/web_frame.tif"

src_dataset, src_band = open_raster(source_file_name)
ort_dataset, ort_band = open_raster(orientation_file_name)

src_srs = get_srs(src_dataset)
new_srs = get_srs(ort_dataset)

print("Source EPSG: " + str(src_srs.GetAuthorityCode(None)))
print("Target EPSG: " + str(new_srs.GetAuthorityCode(None)))

# flush orientation dataset
ort_dataset, ort_band = None, None

# create re-projected raster and save as GeoTIFF
reproj_dataset = reproject_raster(src_dataset, src_srs, new_srs)
reproj_file_name = r"" + os.path.abspath("") + "/geodata/rasters/Fr1000cfs_reproj.tif"
array_data = reproj_dataset.ReadAsArray()
new_epsg = int(new_srs.GetAuthorityCode(None))
geo_transformation = reproj_dataset.GetGeoTransform()
create_raster(reproj_file_name, raster_array=array_data, epsg=new_epsg, geo_info=geo_transformation)
reproj_dataset = None

Source EPSG: 6418
Target EPSG: 3857

Der reprojizierte Froude-number-Raster sieht in QGIS so aus:

python reproject Froude number raster — Figure 3:Der projizierte Raster der Froude-Nummern.

Die reproject_raster-Funktion ist auch in flusstools (geotools.reproject_raster()) erhältlich, wo das Speichern des neuen reprojizierten Rasters in die Funktion eingebettet ist (automatisch die syllable "_epsg[NO]" an den ursprünglichen Dateinamen anhängen).

Note

To display multiple rasters with different coordinate systems on the same map in QGIS, the coordinate systems must be harmonized in most cases. While QGIS will automatically propose to convert a raster with a non-project-compliant CRS to the project CRS, it also has a dedicated function for adjusting raster coordinate systems: In QGIS, click on the Raster menu > Projections > Warp (Reproject).... Select the raster(s) to reproject (i.e., the raster(s) to harmonize with the project coordinate system). However, Warp may not perform all reprojection steps as desired and lead to wrong placements of the new raster. The Warp method is also available in Python through gdal.Warp (read the docs):

kwargs={'format': 'GTiff', 'geoloc': True}
gdal.Warp(TARGET_GEO_TIFF_FILE_NAME, SOURCE_GEO_TIFF_FILE_NAME, **kwargs)

Zonal Statistiken für morphologische Einheiten¶

Bei hydraulischen und georäumlichen Analysen stellt sich häufig die Frage nach statistischen Werten bestimmter Bereiche eines oder mehrerer Raster. Wir können beispielsweise an Mittelwerten und Standardabweichungen in bestimmten Zonen eines Wasserkörpers interessiert sein. Zu diesem Zweck ermöglichen Zonale Statistiken die Abgrenzung einer Fläche eines Rasters unter Verwendung einer Polygonformdatei.

Der River Architect-Datensatz umfasst eine Slackwater-Zone und zonale Statistiken helfen, die mittlere Wassertiefe und Strömungsgeschwindigkeit von Slackwaters zu identifizieren, die eine sogenannte morphologische Einheit sind (recall the create-shapefile section).

Um eine visuell sichtbare Slackwater-Einheit zu analysieren, können wir ein Polygon in einer neuen Formdatei zeichnen, die morphologische Einheiten abgrenzt. Die folgenden Figuren führen durch die Erstellung einer Polygon-Formdatei und die Abgrenzung der Riffle mit QGIS. Starten Sie mit der Eröffnung von QGIS und erstellen Sie ein neues Projekt. Importieren Sie die Wassertiefe und Strömungsgeschwindigkeitsraster mit der langsamen und flachen Wasserzone. Dann folgen Sie dem Workflow in den nachfolgenden Abbildungen.

qgis create new shapefile — Figure 4:QGIS: Erstellen von Ebenen > Neue Formdateiebene...

qgis define new shapefile — Figure 5:Neue Formdateiebene definieren

qgis shapefile toggle-editing — Figure 6:Formdatei bearbeiten aktivieren

qgis draw polygon shapefile — Figure 7:Zeichne Polygone in einer Formdatei in QGIS.

Vervollständigen Sie die Zeichnung, indem Sie auf die * Save Edits Disk Taste (zwischen Toggle Editing und Add Polygon*) klicken. Nur für den Fall, ist die Slackwater-Delineation Polygon Shapefile auch unter dem jupyter-python repository dieses eBook.

Zonal-Statistiken können mit den gdal und ogr Bibliotheken berechnet werden, aber das ist ein wenig umständlich. Die Bibliothek rasterio (conda install -c conda-forge rasterio) bietet eine viel bequemere Methode, um zonale Statistiken mit der Methode rasterstats.zonal_stats(SHP-FILE, RASTER, STATISTICS-TYPES) zu berechnen. Mit zonal_stats können wir problemlos Statistiken über die Wassertiefe und Strömungsgeschwindigkeitsraster in den Grenzen des neuen Slackwater-Polygons erhalten.

import rasterstats as rs
# make file names
h_file = r"" + os.getcwd() + "/geodata/rasters/h001000.tif"
u_file = r"" + os.getcwd() + "/geodata/rasters/u001000.tif"
zone = r"" + os.getcwd() + "/geodata/shapefiles/slackw-poly.shp"

# get water depth stats in zone
h_stats = rs.zonal_stats(zone, h_file, stats=["min", "max", "median", "majority", "sum"])
# get flow velocity stats in zone - note the different stats assignment
u_stats = rs.zonal_stats(zone, u_file, stats="min max median majority sum")

print(h_stats)
print(u_stats)

[{'min': 0.0, 'max': 5.423915386199951, 'sum': 1709.34521484375, 'median': 1.6403688192367554, 'majority': 0.0}]
[{'min': 0.0, 'max': 5.139162540435791, 'sum': 1609.26318359375, 'median': 1.879171371459961, 'majority': 0.0}]

Beachten Sie, dass beide Raster in dem US-üblichen Einheitssystem (d.h. Füße und Füße pro Sekunde) liegen. Weitere Statistiken können mit zonal_stats:
berechnet werden min,max,mean,count,sum,std,median,majority,minority,unique, range,nodata,percentile_<q> (wo <q> jede Schwimmernummer zwischen 0 und 100 sein kann).

Darüber hinaus können benutzerdefinierte Statistiken hinzugefügt werden, wobei das Modul numpy.ma mit seinen Array-Handling-Kapazitäten besonders nützlich ist, z.B. zur Übertragung oder Angabe von Statistiken entlang einer Achse. So können wir beispielsweise eine bestimmte Funktion definieren, um die Standardabweichung wie folgt zu berechnen:

def raster_std(raster_array):
    return np.ma.std(raster_array)

Um die raster_std-Funktion in zonal_stats zu nutzen, schreiben Sie etwas Ähnliches:

u_stats = rs.zonal_stats(
    zone, u_file,
    stats="min max",
    add_stats={"stdev": raster_std}
)
print(u_stats)

[{'min': 0.0, 'max': 5.139162540435791, 'stdev': np.float64(1.1065991101701524)}]

Clip eines Rasters¶

Die oben vorgestellte rasterstats.zonal_stats Methode arbeitet mit “Mini-Rasters”, die Clips des Eingaberasters zur benutzerdefinierten Polygonformdatei darstellen. Darüber hinaus kann ein Miniraster selbst durch die Definition des optionalen Keyword-Arguments raster_out=True erhalten werden. Für den Fall, dass wir den originalen Raster ohne und statistische Operation verclipsen möchten, können wir einen kleinen Trick verwenden, indem wir eine zusätzliche Statistikfunktion definieren, die das ursprüngliche Array zurückgibt:

def original(raster_array):
    return raster_array

Mit raster_out=True und der original()-Funktion können wir den Clip-Raster in den folgenden Array-Typen abrufen:

mini_raster_array - ein abgecliptes und maskiertes Numpy-Array,
mini_raster_affine - eine Transformation als Affine-Objekt (d.h. mit affine 6-Parameter transform) und
mini_raster_nodata - NoData

Der folgende Codeblock verdeutlicht die Verwendung zum Abrufen eines Numpy-Arrays:

import rasterstats as rs

h_file = r"" + os.getcwd() + "/geodata/rasters/h001000.tif"
h_stats = rs.zonal_stats(zone, h_file, stats="count",
                         add_stats={"original": original},
                         raster_out=True)
print(h_stats[0].keys())
print(h_stats[0]["mini_raster_array"])

dict_keys(['count', 'original', 'mini_raster_array', 'mini_raster_affine', 'mini_raster_nodata'])
[[-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 ...
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]]

Slope / Aspect Maps und integrierte Kommandozeilenskripte¶

Hillslope-Karten sind ein wichtiger Parameter in der Hydraulik, Hydrologie und Ökologie. Beispielsweise bestimmt die Steigung die Strömungsrichtung des Wassers und es ist auch ein Kriterium zur Abgrenzung des Lebensraums vieler Arten. Zur Berechnung von Steigungsgradienten und -richtungen hat gdal ein Kommandozeilen-Tool namens gdaldem, das einen Digitales Oberflächenmodell (DOM) (Digital Elevation Model)-Raster benötigt. Die allgemeine Nutzung von gdaldem in der Befehlszeile ist (Argumente in Klammern sind optional):

gdaldem slope input_dem output_slope_map  [-p use percent slope (default=degrees)] [-s scale* (default=1)] [-alg ZevenbergenThorne] [-compute_edges] [-b Band (default=1)] [-of format] [-co "NAME=VALUE"]* [-q]

To call the command-line tool, we can either use a terminal/command prompt or Python’s standard library subprocess. The following code block illustrates the usage of the gdaldem command line tool through subprocess.call() to create a slope raster (in percent) from the River Architect sample data’s dem.tif. Note that subprocess.call() returns 0 if the command execution was successful and any other return value indicates an error.

import subprocess, os

cmd_create_slope = "gdaldem slope {0}/geodata/rasters/dem.tif {0}/geodata/rasters/slope-percent.tif -p".format(os.path.abspath(""))
subprocess.call(cmd_create_slope, shell=True)

0...10...20...30...40...50...60...70...80...90...100 - done.

0

python create slope raster Geotiff — Figure 8:Der neue Hangraster.

Neben der absoluten Steigung (d.h. Neigung in Grad) oder anstelle der Steigung kann es wichtig sein, die Neigungsrichtung (z.B. Neigung nach Süden, Westen, Osten oder Norden) zu kennen. Ein Raster, der die Neigungsrichtung anzeigt, wird als Aspektraster bezeichnet, wobei der Süden 0° (und 360°), westlich von 90°, nördlich von 180° und östlich von 270° entspricht. Ein Aspekt-Raster kann auch mit gdaldem erstellt werden:

gdaldem aspect input_dem output_aspect_map [-trigonometric] [-zero_for_flat] [-alg ZevenbergenThorne] [-compute_edges] [-b Band (default=1)] [-of format] [-co "NAME=VALUE"]* [-q]

Um einen Aspektraster der River Architect Musterdaten zu erstellen, führen Sie Digitales Oberflächenmodell (DOM) aus:

cmd_create_aspect = "gdaldem aspect {0}/geodata/rasters/dem.tif {0}/geodata/rasters/slope-aspect.tif".format(os.path.abspath(""))
subprocess.call(cmd_create_aspect, shell=True)

0...10...20...30...40...50...60...70...80...90...100 - done.

0

python create aspect raster Geotiff — Figure 9:Der neue Aspekt raster.

Least Cost Path zwischen Pixeln (und einer anderen Art von Reprojektion)¶

Ökohydraulik Hintergrund¶

Tiefstkostenpfade sind wichtig, um effiziente Routen für die Navigation (z.B. im Auto) zu planen und sie können auch bei Ökohydraulik hilfreich sein. Nehmen wir einen Moment die Position eines Fisches, der nach einer Flut mit abnehmender Entladung so schnell wie möglich vom Hochwasserboden zurück in den Hauptkanal schwimmen will, wo genügend Wasser vorhanden ist. In der nachfolgenden Abbildung zeigt Punkt 1 den Ausgangspunkt auf dem Hochwasser und Punkt 2 das Ziel im Hauptkanal. Der rötliche Hintergrund stellt den oben produzierten Hangraster (slope-percent.tif) dar und die Wassertiefe bei durchschnittlichen jährlichen Entladungen in blau gefärbt.

Natürlich entspricht der Weg der geringsten Kosten dem Weg der steilsten, monoton nach unten gerichteten Steigung und wir werden annehmen, dass ein Fisch ihn finden kann.

Funktionen und Bibliotheken beteiligt¶

The skimage (scikit-image) library (cf. Other packages in the Open source libraries section) provides with skimage.graph.route_through_array a smart method to calculate the least cost path by summing up pixel-wise connections from point 1 to point 2.

Hier ist, wie es funktioniert: Nehmen Sie ein Numpy-Array (z.B. mit zufälligen Steigungswerten) an, das so aussieht:

slope_image = np.random.randint(100, size=(3, 5))
slope_image

array([[43, 68, 67, 88, 60],
       [22, 95, 90, 97,  0],
       [38, 18, 83, 65, 83]])

Um den schnellsten Weg vom Array-Index [0][0] (oben point_1 = (0, 0)) bis zum Array-Index [2][4] (unten rechten Punktpoint_2 = (2, 4)) zu finden, können wir mit route_through_array() eine list (least_cost_path_indices) mit den Array-Koordinaten des Pfades zu erhalten und die Kosten (weight) involviert (sum aller Pixel des geringsten Kostenpfads):

from skimage.graph import route_through_array

point_1 = (0, 0)
point_2 = (2, 4)
least_cost_path_indices, weight = route_through_array(slope_image, point_1, point_2)
least_cost_path_indices, weight

([(0, 0), (1, 0), (2, 1), (2, 2), (2, 3), (2, 4)],
 np.float64(259.2842712474619))

Um die Mindestkostenpfadliste in ein Array zu integrieren, das wir rasterize (geotools.create_raster()) einbinden können, können wir least_cost_path_indices in ein numpy-zeros-Array des ursprünglichen Hangrasters (Bild) als transponierte Liste einfügen.

least_cost_path_indices = np.array(least_cost_path_indices).T
least_cost_path_array = np.zeros_like(slope_image)
least_cost_path_array[least_cost_path_indices[0], least_cost_path_indices[1]] = 1
least_cost_path_array

array([[1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [0, 1, 1, 1, 1]])

In der Praxis ist der Hangraster georeferiert, weshalb wir Pixelkoordinaten relativ zum Koordinatensystemursprung verwenden müssen. Dazu benötigen wir zwei weitere Funktionen:

Eine Funktion zur Berechnung des Pixel-Index-verwandten Offsets, den wir coords2offset: Die coords2offset()-Funktion wird die x-y-Schicht in Form von “Zahl der Pixel” zurückgeben (zwei Integers, eine für x und eine für y-Schicht).
Die above-defined get_srs()-Funktion (d.h. geotools.get_srs()).

Die coords2offset()-Funktion sieht so aus:

def coords2offset(geo_transform, x_coord, y_coord):
    """
    Returns x-y pixel offset
    :param geo_transform: osgeo.gdal.Dataset.GetGeoTransform() object
    :param x_coord: FLOAT of x-coordinate
    :param y_coord: FLOAT of y-coordinate
    :return: offset_x, offset_y (both integer of pixel numbers)
    """
    origin_x = geo_transform[0]
    origin_y = geo_transform[3]
    pixel_width = geo_transform[1]
    pixel_height = geo_transform[5]
    offset_x = int((x_coord - origin_x) / pixel_width)
    offset_y = int((y_coord - origin_y) / pixel_height)
    return offset_x, offset_y

Die coords2offset()-Funktion konvertiert ein Rasterfeld (z.B. mit der oben definierten geotools.raster2array()-Funktion) in ein Array, das mit route_through_array() mit folgendem Workflow verwendet werden kann:

Verwenden Sie die raster’s geo_transform (gdal.Dataset.GetGeoTransform = (origin_x, pixel_width, 0, origin_y, 0, pixel_height)) und die Start- und Endpunktkoordinaten (d.h.start_coord von Punkt 1 und stop_coord von Punkt 2) in coords2offset(), um ihre Pixel-Indizes (start_index_x, start_index_y, stop_index_x und stop_index_y) im Rasterfeld zu erhalten.
Ersetzen Sie np.nan-Werte im Rasterfeld mit Werten, die höher sind als der maximale Wert des Arrays. Verwenden Sie keine Nullen, denn wir möchten die np.nanpixels später durch Überschreiben np.nan mit sehr hohen Pixelkosten aus dem kostengünstigsten Pfad ausschließen.
Verwenden Sie route_through_array(), wie oben mit den optionalen Argumenten geometric=True erklärt (verwenden Sie die MCP Geometric class anstatt MCP base, um Kosten zu berechnen) und fully_connected=True (verwenden Sie Diagonalpixel als direkte Nachbarn).
Integrieren Sie die Mindestkostenpfadliste (index_path) in ein numpy-zeros-Array (Kind von raster_array, wie oben erläutert) und geben Sie die path_array zurück.

def create_path_array(raster_array, geo_transform, start_coord, stop_coord):
    # transform coordinates to array index
    start_index_x, start_index_y = coords2offset(geo_transform, start_coord[0], start_coord[1])
    stop_index_x, stop_index_y = coords2offset(geo_transform, stop_coord[0], stop_coord[1])

    # replace np.nan with max raised by an order of magnitude to exclude pixels from least cost
    raster_array[np.isnan(raster_array)] = np.nanmax(raster_array) * 10

    # create path and costs
    index_path, cost = route_through_array(raster_array, (start_index_y, start_index_x),
                                               (stop_index_y, stop_index_x),
                                               geometric=True, fully_connected=True)


    index_path = np.array(index_path).T
    path_array = np.zeros_like(raster_array)
    path_array[index_path[0], index_path[1]] = 1
    return path_array

Anwendung¶

Recall, wir haben die folgenden Funktionen definiert (alle sind in flusstools durch from flusstools import geotools) verfügbar, die wir für die Berechnung des am wenigsten kostenpfads nutzen können, um von Punkt 1 bis Punkt 2 im slope-percent.tif raster zu erhalten:

geotools.raster2array()
geotools.create_path_array()
geotools.get_srs()
geotools.create_raster()

Der folgende Codeblock verwendet diese Funktionen wie folgt:

Define Input (slope-percent.tif) und Output (least cost.tif) Rasternamen (mit Verzeichnissen).
Definieren Sie die Koordinaten der Punkte 1 und 2 als tuples (x, y) in der EPSG:6418 Projektion.
Laden Sie die Eingabe raster(src_raster), ihre Band als Array (raster_array) und Geotransformation (geo_transform) mit der Funktion raster2array() ein.
Holen Sie sich den mit einem in einer Reihe von Nullen angezeigten Mindestkostenpfad (d.h. ein on-off path_array) mit der create_path_array() Funktion.
Erhalten Sie die osgeo.osr.SpatialReference des Eingaberasters (src_raster = osgeo.gdal.Dataset("slope-percent.tif")).
Erstellen Sie den Mindestkostenpfad GeoTIFF raster mit der create_raster()-Funktion als gdal.GDT_Byte-Band.

from skimage.graph import route_through_array

# define raster input and out names
in_raster_name = r"" + os.path.abspath("") + "/geodata/rasters/slope-percent.tif"
out_raster_name = r"" + os.path.abspath("") + "/geodata/rasters/least_cost.tif"
# define coordinates of points 1 and 2 (in EPSG:6418)
point_1_coord = (6749261.94092826917767525, 2206970.35179582564160228)  
point_2_coord = (6749016.82820663042366505, 2207050.61491037486121058)

# get source raster (osgeo.gdal.Dataset), the raster as nd.array, and the geotransformation tuple
src_raster, raster_array, geo_transform = raster2array(in_raster_name)
# get the zeros-like array with least cost pixels = 1
path_array = create_path_array(raster_array, geo_transform, point_1_coord, point_2_coord)
# get the spatial reference system of the input raster (slope-percent.tif)
src_srs = get_srs(src_raster)
# project the least cost path_array into a Byte (only zeros and ones) raster
create_raster(out_raster_name, path_array, epsg=int(src_srs.GetAuthorityCode(None)),
              rdtype=gdal.GDT_Byte, geo_info=geo_transform)

python create least cost path — Figure 12:Der am wenigsten in QGIS dargestellte Kostenpfad.

Legitimativ, Sie können sich fragen, ob es besser war, den geringsten Kostenpfad als Linie zu repräsentieren. Natürlich ist das richtig. Diese Operation ist jedoch eine Umwandlung eines Rasters in eine Linienformdatei, die im nächsten Abschnitt unter geodata conversion erläutert wird. Kuriose Leser können auch direkt die raster2line() Funktion flusstools(oder im flusstools.geotools/dataset mgmt.pyskript) nutzen.

References¶

Wyrick, J. R., & Pasternack, G. B. (2014). Geospatial organization of fluvial landforms in a gravel-cobble river: Beyond the riffle-pool couplet. Geomorphology, 213(Supplement C), 48–65. 10.1016/j.geomorph.2013.12.040
Schwindt, S., Larrieu, K., Pasternack, G. B., & Rabone, G. (2020). River Architect. SoftwareX, 11, 100438. 10.1016/j.softx.2020.100438