Manipulation des grilles - 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.

Charger le raster¶

Ouvrir les données de grille existantes¶

Les données Raster peuvent être ouvertes en code Python en tant qu’exemple de gdal.Open("FILENAME"). Le bloc de code suivant fournit une fonction pour ouvrir tout raster spécifié avec l’argument d’entrée file_name. L’un des éléments les plus importants lors du traitement des données de raster est la bande de raster, qui assume un rôle de support de données similaire à GetLayer dans le traitement de shapefile. Pour accéder à la bande raster, la fonction open_raster:

Active l’erreur et la rétroaction d’avertissement avec gdal.UseExceptions() (cette étape est vitale lorsque vous utilisez gdal).
Ouvre les déclarations de raster file_name embrassées par try - except pour indiquer si et pourquoi une erreur potentielle s’est produite en ouvrant le raster.
Ouvre le numéro de bande de raster indiqué dans le paramètre optionnel band_number avec raster_band = raster.GetRasterBand(band_number) (la valeur par défaut est 1).
Retourne les objets raster et raster band.

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

Pour utiliser la fonction open_raster, appelez-la avec un nom de fichier comme indiqué dans le bloc de code suivant avec le h001000.tif raster de l’échantillon River Architect data. Le script ferme immédiatement le raster à nouveau en écraseant la variable avec un None-instance pour éviter que le fichier GeoTIFF ne soit verrouillé par la suite.

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> >

Statistiques de bande de raster et Scripts de boîte à outils¶

Une fois que le raster et sa bande sont chargés avec la fonction open_raster ci-dessus, nous pouvons accéder à des informations statistiques (p. ex. le minimum ou le maximum), identifier la valeur no-data (c.-à-d. une valeur prédéfinie qui est attribuée à des pixels sans valeur), ou le type d’unités utilisées.

Les scripts Python pour le traitement des données géospatiales peuvent également être intégrés en tant que plugins dans les applications de bureau SIG (par exemple, en tant que plugins dans QGIS ou Toolbox dans ArcGIS Pro). Pour exécuter un script Python dans une application de bureau SIG, il doit être écrit comme un script autonome qui peut recevoir des arguments d’entrée. Le lecteur intéressé peut en savoir plus sur la mise en oeuvre des plugins dans QGIS dans le QGIS docs.

Ici, nous n’écrirons que le prochain bloc de code afin qu’il puisse être exécuté dans une application console/terminale en tant que script autonome (appelez le 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]))

Pour exécuter ce script, enregistrez-le en tant que raster band info.py (p. ex., C:\temp dans Windows ou ~/temp/ dans Linux) et naviguez vers le répertoire du script dans une interface terminal (p. ex., dans le terminal de PyCharm, Anaconda Prompt ou le terminal Linux) en utilisant la commande cd. Puis exécutez le script pour obtenir des statistiques de la profondeur d’eau raster h001000.tif avec (dans Windows):

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

Sur Linux et avec le flussenv activated, par exemple:

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

Créer et enregistrer une grille (à partir de Array)¶

Chauffeurs¶

Tout comme pour les fichiers shapefile, le pilote gdal approprié (analogue à ogr drivers) doit être chargé pour enregistrer un raster. Pour obtenir une liste complète des pilotes gdal raster :

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

Types de données¶

Les pixels raster de sortie peuvent être l’un des types de données suivants (source: gdal.org/doxygen/):

GDT_Unknown Unknown or unspecified type
GDT_Byte 8 bits entier non signé
GDT_UInt16 16 bits entier non signé
GDT_Int16 16 bits signés entier
GDT_UInt32 32 bits entier non signé
GDT_Int32 32 bits signés entier
GDT_Float32 32 bits point flottant
GDT_Float64 64 bits point flottant
GDT_CInt16 Int16 complexe
GDT_CInt32 Complexe Int32
GDT_CFloat32 Complexe Float32
GDT_CFloat64 Complexe Float64

Créer un Raster (Array to Raster)¶

Connaissant les bases de la manipulation de raster, les types de données et Python, nous pouvons créer un raster à partir d’un tableau numérique. Comme un raster est essentiellement un tableau géoréférencé, il est pratique de convertir un numpy array en raster (bande). Les blocs de fonctions disposent de la conversion d’un tableau numpy en un raster GeoTIF selon le workflow suivant:

Consultez le pilote GeoTIFF (driver = gdal.GetDriverByName('GTiff')).
Récupérer la taille du tableau et (nombre de lignes rows et colonnes cols).
Créer un nouveau master GeoTIFF (new_raster = driver.Create(file_name, cols, rows, 1, eType=rdtype)), où
- file_name est le répertoire et le nom du nouveau fichier raster se terminant par .tif (par exemple, "C:\\temp\\rasters\\new.tif").
- cols, rows représentent la forme du tableau, et eType est le type de données géospatiales (voir la liste ci-dessus)
Définir l’origine géographique, qui est stockée dans le paramètre origin (tuple) et définir le pixel_width et pixel_height (unités de pixels définies avec srs - voir ci-dessous).
Remplacer np.nan valeurs dans le tableau numpy par nan_value.
Activer un objet band, définir le NoDataValue à nan_value, et écrire le tableau à band.
Créez un objet système de référence spatiale (srs) en fonction du paramètre d’entrée epsg et exportez-le au format WKT.
Relâchez le raster (flush from cache).

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()

Pour appeler la fonction pour écrire un tableau aléatoire, nous pouvons maintenant utiliser la fonction create_raster() (également disponible à partir de 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:Le nouveau raster contenant une hauteur de point unis dem aléatoire.

Raster Calculus (Raster / Band to Array)¶

La procédure décrite avec la fonction create_raster() peut être utilisée vice-versa, pour créer numpy array à partir des bandes raster. Le raster converti en un tableau numpy permet d’appliquer des opérations algébriques ou d’autres opérations logiques aux données de raster existantes.

Besoin d’un exemple ? Dans le RiverArchitect SampleData, les unités de la profondeur d’eau raster h001000.tif sont en pieds habituels aux États-Unis et les unités de la vitesse d’écoulement raster u001000.tif sont en pieds par seconde. Cependant, pour calculer la Froude number (impliquant la constante de gravité) pour chaque pixel en fonction des deux rasters (profondeur d’eau et vitesse de débit), il est pratique de convertir les deux rasters en m et m/s, respectivement. À cette fin, le bloc de code suivant comporte une autre fonction réutilisable et personnalisée qui charge un raster comme tableau et écrase NoDataValues avec np.nan (raster and band peut être instantanée avec la fonction open_raster ci-dessus):

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()

La fonction raster2array() est également incluse dans Flusstools: geotools.raster2array()

Enfin, pour créer un numéro de Froude GeoTIFF raster, le bloc de code suivant utilise la fonction raster2array pour convertir la profondeur d’eau et la vitesse d’écoulement des rasters GeoTIFF en un tableau numpy, effectuant des calculs algébriques simples pour convertir les rasters en m et m/s, respectivement, et enregistrer le fichier GeoTIFF résultant. En détail, le workflow consiste à :

Définir les noms de fichiers master avec les répertoires (h_file et u_file),
Charger les rasters originaux comme ndarray avec la fonction raster2array() et obtenir la description GeoTransform originale,
Convertir toutes les valeurs des pieds habituels des États-Unis en unités métriques S.I. (appeler la fonction feet_to_meter() à partir des bases de Python), et
Enregistrez une nouvelle copie du raster.

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:Le nombre Froude raster calculé avec la vitesse d’écoulement et la profondeur de l’eau rasters.

Reprojeter une grille¶

La transformation (et la reprojection) d’un raster en un autre système de coordonnées implique la rotation, le déplacement et le cisaillement des pixels. Si l’une de ces opérations est ignorée, le raster reprojecté peut être serré, tordu ou placé n’importe où dans le monde, mais pas où il devrait être placé. Par conséquent, l’approche pour la reprojection d’un raster dans un autre système de référence de coordonnées implique les étapes suivantes:

Récupérer les systèmes de référence spatiale source et cible (p. ex., dériver d’un code d’autorité gdal.Dataset ou EPSG).
Lire la transformation géographique de l’ensemble de données source ( gdal.Dataset.GetGeoTransform()).
Calculer le nombre de pixels et l’espacement entre les pixels dans le nouvel ensemble de données (réaménagé).
Activez le nouvel ensemble de données (réaménagé).
Projeter une image de l’ensemble de données source sur le nouvel ensemble de données (réaménagé) (gdal.ReprojectImage()).

Le système de référence spatiale peut être dérivé d’un ensemble de données avec les explications dans la section shapefile en écrivant une fonction get_srs(). Le bloc de code suivant indique la fonction get_srs() (utilise la bibliothèque osr de osgeo / gdal ), qui est également intégrée à flusstools (geotools.get_srs()).

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

Avec les fonctions open_raster() et get_srs(), nous avons tous les ingrédients nécessaires pour accomplir le travail de reprojection de raster dans une autre funcation appelée reproject_raster(). Une autre caractéristique de la fonction est qu’elle assure l’utilisation correcte de osr.CoordinateTransformation, qui se comporte différemment sous gdal 3.0 par rapport aux anciennes versions gdal (en savoir plus sur la page GitHub d’OSGeo).

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

L’utilisation de la fonction reproject_raster() dans un script Python nécessite un jeu de données source et un autre ensemble de données (orientation) avec le nouveau système de coordonnées dans lequel le jeu de données source sera projeté. L’exemple suivant montre comment re-projecter le raster Froude-numéro créé ci-dessus dans le système de coordonnées EPSG=3857 pour le visionner en QGIS sur une carte de base Google Satellite (recall basemaps in QGIS tutorial). Comme ensemble de données d’orientation, utilisez web frame.tif, qui a été créé sur la carte de base Google Satellite.

Avec la fonction get_srs() qui détecte automatiquement la projection de raster et le système de référence spatiale, nous pouvons utiliser la fonction create_raster() pour reprojecter la fonction Fr1000cfs.tif raster (par exemple, epsg=4326).

# 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

Plotté dans QGIS, le reprojected Froude-number raster ressemble à ceci:

python reproject Froude number raster — Figure 3:Le raster prévu des numéros Froude.

La fonction reproject_raster est également disponible (légèrement modifiée) dans flusstools (geotools.reproject_raster()), où la sauvegarde du nouveau raster reprojecté est intégrée dans la fonction (Appendise automatiquement la syllabe "_epsg[NO]" au nom du fichier original).

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)

Statistiques zonales pour les unités morphologiques¶

Dans les analyses hydrauliques et géospatiales, la question des valeurs statistiques de certaines zones d’un ou de plusieurs rasters se pose souvent. Par exemple, nous pourrions nous intéresser aux valeurs moyennes et aux écarts types dans des zones spécifiques d’un plan d’eau. À cette fin, les statistiques zonales* permettent de délimiter une zone de raster en utilisant un fichier de forme polygone.

L’ensemble de données River Architect comprend une zone d’eau douce et des statistiques zonales permettant d’identifier la profondeur moyenne d’eau et la vitesse d’écoulement des eaux douces, qui sont une unité morphologique (recall the create-shapefile section).

Background

Les unités morphologiques instream aident à décrire l’organisation géospatiale des formes de terres fluviales, qui jouent un rôle important dans les analyses écohydrauliques et la restauration des rivières. Par exemple, les unités pool décrivent les zones d’eau profonde à faible vitesse d’écoulement, les unités riffle sont typiquement caractérisées par des profondeurs d’eau peu profondes et une vitesse élevée, et les unités slackwater sont des zones d’écoulement peu profondes à faible vitesse d’écoulement (beaucoup de juvéniles aiment les eaux de relâche). Wyrick & Pasternack (2014) présentez la délimitation des unités morphologiques et un résumé de l’accès libre se trouve dans l’annexe 5](https://ars.els-cdn.com/content/image/1-s2.0-S235271101930281X-mmc1.pdf) à Schwindt et al. (2020).

Pour analyser une unité d’eau douce apparente visuellement, nous pouvons dessiner un polygone dans un nouveau shapefile qui délimite les unités morphologiques. Les figures suivantes guident à travers la création d’un fichier de forme polygone et la délimitation de la trémie avec QGIS. Commencez par ouvrir QGIS et créez un nouveau projet. Importez la profondeur d’eau et la vitesse d’écoulement des rasters montrant la zone d’eau lente et peu profonde. Suivez ensuite le flux de travail dans les illustrations ci-dessous.

qgis create new shapefile — Figure 4:QGIS : Créer un calque > Nouveau calque...

qgis define new shapefile — Figure 5:Définir le nouveau calque Shapefile

qgis shapefile toggle-editing — Figure 6:Activer l’édition de shapefile

qgis draw polygon shapefile — Figure 7:Dessiner les polygones dans un fichier de forme en QGIS.

Finalisez le dessin en cliquant sur le bouton de disque Enregistrer Éditions (entre Toggle Editing et Ajouter Polygon). Juste au cas où, le polygone de délimitation de l’eau shapefile est également disponible à le dépôt jupyter-python de ce livre électronique.

Les statistiques zonales peuvent être calculées en utilisant les bibliothèques gdal et ogr, mais c’est un peu lourd. La bibliothèque rasterio (conda install -c conda-forge rasterio) offre une méthode beaucoup plus pratique pour calculer les statistiques zonales avec sa méthode rasterstats.zonal_stats(SHP-FILE, RASTER, STATISTICS-TYPES). Avec zonal_stats, nous pouvons facilement obtenir des statistiques sur la profondeur de l’eau et la vitesse d’écoulement rasters dans les limites du nouveau polygone d’eau de relâche.

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}]

Rappelons que les deux rasters sont dans le système d’unité habituel américain (c.-à-d. pieds et pieds par seconde). D’autres statistiques peuvent être calculées avec zonal_stats:
min, max, mean, count, sum, std, median, majority, minority, unique, range, nodata, percentile_<q> (où <q> peut être n’importe quel numéro de flotteur entre 0 et 100).

En outre, des statistiques définies par l’utilisateur peuvent être ajoutées, où le module numpy.ma est particulièrement utile avec ses capacités de gestion de tableaux, par exemple pour transposer ou spécifier des statistiques le long d’un axe. Par exemple, nous pouvons définir une fonction spécifique pour calculer l’écart type comme suit:

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

Pour utiliser la fonction raster_std dans zonal_stats, écrivez quelque chose de similaire à ceci :

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)}]

Couper une grille¶

La méthode ci-dessus introduite rasterstats.zonal_stats fonctionne avec “Mini-Rasters”, qui représentent des clips du raster d’entrée au fichier de forme polygone défini par l’utilisateur. De plus, un mini-raster lui-même peut être obtenu en définissant le paramètre optionnel raster_out=True. Dans le cas où nous voulons obtenir le raster original clippé sans et opération statistique, nous pouvons utiliser un petit truc en définissant une fonction statistique supplémentaire qui retourne le tableau original:

def original(raster_array):
    return raster_array

Avec raster_out=True et la fonction original(), nous pouvons récupérer le raster clippé dans les types de tableau suivants:

mini_raster_array - un tableau de numpy coupé et masqué,
mini_raster_affine - une transformation en tant qu’objet Affine (c.-à-d. avec la transformation de affine 6 paramètres), et
mini_raster_nodata - NoData valeurs.

Le bloc de code suivant illustre l’utilisation pour récupérer un tableau numpy :

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 Cartes et scripts de ligne de commande intégrés¶

Les cartes des pentes sont un paramètre important en hydraulique, en hydrologie et en écologie. Par exemple, la pente détermine la direction d’écoulement de l’eau et c’est aussi un critère de délimitation de l’habitat de nombreuses espèces. Pour calculer les pentes et les directions, gdal a un outil en ligne de commande appelé gdaldem , qui nécessite un raster Modèle numérique de terrain (MNT) (Modèle d’élévation numérique). L’utilisation générale de gdaldem dans la ligne de commande est (les arguments entre parenthèses sont facultatifs):

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:La nouvelle pente raster.

En plus de la pente absolue (c.-à-d., inclinaison en degrés), ou au lieu de la pente, il peut être important de connaître la direction de la pente (p. ex., inclinaison vers le sud, l’ouest, l’est ou le nord). Un raster indiquant la direction de la pente s’appelle aspect raster, où le sud correspond à 0° (et 360°), à l’ouest à 90°, au nord à 180° et à l’est à 270°. Un aspect raster peut également être créé avec gdaldem:

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]

Pour créer un aspect raster de l’échantillon de données River Architect exécuter Modèle numérique de terrain (MNT):

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:Le nouvel aspect raster.

Chemin le moins coûteux entre Pixels (et une autre voie de reprojection)¶

Contexte écologique¶

Les itinéraires les moins coûteux sont importants pour planifier des itinéraires efficaces pour la navigation (par exemple, dans une voiture) et ils peuvent également être utiles en écohydraulique. Prenons un instant la position d’un poisson qui, après une inondation avec une décharge décroissante, veut nager aussi vite que possible de la plaine inondable de nouveau dans le chenal principal où il y a assez d’eau. Dans la figure ci-dessous, le point 1 indique le point de départ de la plaine inondable et le point 2 la destination du chenal principal. Le fond rougeâtre représente le raster de pente (p. cent.tif) et la profondeur de l’eau au débit annuel moyen est colorée en bleu.

Naturellement, la trajectoire des moindres coûts correspond à celle de la pente la plus raide et la plus monotone vers le bas et nous supposerons qu’un poisson est capable de la trouver.

Note

Dans de nombreuses rivières, les poissons font face au défi quotidien de s’échapper du piège mortel des dépressions latérales. La raison en est que de nombreuses centrales hydroélectriques provoquent des fluctuations brutales du débit en raison des fluctuations de la demande et de la production d’énergie (appelées hydro-peaking). Par conséquent, le risque d’échouement est élevé dans de nombreux cours d’eau réglementés qui sont soumis à l’hydropeaking. La figure suivante illustre les zones de risque d’échouement en fonction du rejet (en pieds cubes par seconde) dans le cours inférieur du fleuve Yuba (Californie, États-Unis).

Figure 11:Carte d’échouement de poissons créée avec River Architect (source : le Wiki Architecte de rivière / Kenneth Larrieu](https://riverarchitect.github.io/RA_wiki/StrandingRisk)).

Fonctions et bibliothèques impliquées¶

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.

Voici comment ça marche : Supposez un tableau numpy (par exemple, avec des valeurs aléatoires de pente) qui ressemble à ceci :

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]])

Pour trouver le moyen le plus rapide de l’index de tableau [0][0] (en haut à gauche point_1 = (0, 0)) à l’index de tableau [2][4] (en bas à droite point_2 = (2, 4)), nous pouvons utiliser route_through_array() pour obtenir une list (least_cost_path_indices) avec les coordonnées de tableau du chemin à suivre, et les coûts (weight) impliqués (somme de tous les pixels du chemin le moins coûteux):

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))

Pour intégrer la liste des chemins les moins coûteux dans un tableau que nous pouvons * rasterize* (geotools.create_raster()), nous pouvons coller least_cost_path_indices dans un tableau numpy-zeros du raster de pente original (image) comme liste transposée.

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]])

Dans la pratique, le raster de pente est géoréférencé, et par conséquent, nous devons utiliser les coordonnées de pixel par rapport à l’origine du système de coordonnées. À cette fin, nous avons besoin de deux autres fonctions :

Une fonction pour calculer le décalage lié à l’index pixel que nous nommerons coords2offset: La fonction coords2offset() retournera le décalage x-y sous la forme du « nombre de pixels » (deux integers, un pour x et un pour y shift).
La fonction above-defined get_srs() (i.e., geotools.get_srs()).

La fonction coords2offset() ressemble à ceci :

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

La fonction coords2offset() convertit un tableau raster (par exemple, produit avec la fonction geotools.raster2array() ci-dessus) en un tableau qui peut être utilisé avec route_through_array() avec le workflow suivant :

Utilisez les coordonnées geo_transform (gdal.Dataset.GetGeoTransform = (origin_x, pixel_width, 0, origin_y, 0, pixel_height)) et les coordonnées de début et de fin (i.e., start_coord du point 1 et stop_coord du point 2) à coords2offset() pour obtenir leurs indices de pixel (start_index_x, start_index_y, stop_index_x et stop_index_y) dans le tableau de raster.
Remplacer les valeurs np.nan dans le tableau raster par des valeurs supérieures à la valeur maximale du tableau. N’utilisez pas de zéros, parce que nous voulons exclure le np.nan pixels du chemin le moins cher plus tard en écraseant np.nan avec des coûts de pixel très élevés.
Utiliser route_through_array() comme expliqué ci-dessus avec les arguments optionnels geometric=True (utiliser la MCP Geometric class plutôt que MCP base pour calculer les coûts) et fully_connected=True (enable en utilisant des pixels diagonaux comme voisins directs).
Intégrer la liste des chemins les moins coûteux (index_path) dans un tableau numpy-zeros (enfant de raster_array, comme expliqué ci-dessus) et retourner le path_array.

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

Demande¶

Rappelons, nous avons défini les fonctions suivantes (toutes sont disponibles dans flusstools via from flusstools import geotools) que nous pouvons utiliser pour le calcul du chemin le moins coûteux pour obtenir du point 1 au point 2 dans le raster slope-percent.tif:

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

Le bloc de code ci-dessous utilise ces fonctions comme suit:

Définir les entrées (slope-percent.tif) et les sorties (least cost.tif) les noms de raster (avec répertoires).
Définir les coordonnées des points 1 et 2 comme tuples (x, y) dans la projection EPSG:6418.
Charger l’entrée raster(src_raster), sa bande en tableau (raster_array) et la géotransformation (geo_transform) avec la fonction raster2array().
Obtenez le chemin le moins cher indiqué avec ceux dans un tableau de zéros (c.-à-d. un path_array) avec la fonction create_path_array().
Obtenez le osgeo.osr.SpatialReference de l’entrée raster (src_raster = osgeo.gdal.Dataset("slope-percent.tif")).
Créer le chemin le moins coûteux GeoTIFF raster avec la fonction create_raster() comme une bande gdal.GDT_Byte.

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:Le chemin le moins coûteux illustré dans QGIS.

On peut légitimement se demander s’il valait mieux représenter le chemin le moins coûteux en tant que ligne. Bien sûr, c’est exact. Cependant, cette opération est une conversion d’un raster en un fichier de forme de ligne, qui est expliqué dans la section suivante sur geodata conversion. Les lecteurs curieux peuvent aussi utiliser directement la fonction raster2line() de flusstools (ou jeter un coup d’oeil dans le script flusstools.geotools/dataset mgmt.py).

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