Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Vectoriser et rastériser

This section introduces geospatial dataset conversion with Python. In particular, the goal of this section is to guide to an understanding of conversions from raster to vector data formats and vice versa. For interactive reading and executing code blocks Binder and find geo03-conversion.ipynb, or install Python and JupyterLab locally.

Watch this section as a video

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

Importer les bibliothèques pertinentes

Assurez-vous d’importer les paquets pertinents pour la manipulation des rasters, des shapefiles et des références géospatiales:

from osgeo import gdal
from osgeo import osr
from osgeo import ogr

Vecteurs

Raster à Line

Dans cette section, nous convertissons l’ensemble de données least cost path raster (least cost.tif) en un fichier de forme (poly). À cette fin, nous écrivons d’abord une fonction appelée offset2coords(), qui représente l’inverse de la fonction coords2offset(), et convertit l’offset x/y (en numéros de pixel integer) en coordonnées de la géo-transformation d’un ensemble de données géospatial :

def offset2coords(geo_transform, offset_x, offset_y):
    # get origin and pixel dimensions from geo_transform (osgeo.gdal.Dataset.GetGeoTransform() object)
    origin_x = geo_transform[0]
    origin_y = geo_transform[3]
    pixel_width = geo_transform[1]
    pixel_height = geo_transform[5]
    
    # calculate x and y coordinates
    coord_x = origin_x + pixel_width * (offset_x + 0.5)
    coord_y = origin_y + pixel_height * (offset_y + 0.5)

    # return x and y coordinates
    return coord_x, coord_y

Ensuite, nous pouvons écrire une fonction de base pour convertir un jeu de données raster en un fichier de forme de ligne. Nous nommons cette fonction raster2line() et elle s’appuie sur le workflow suivant:

  • Ouvrir un raster, son groupe comme array, et geo_transform (géo-transformation) défini avec l’argument raster_file_name dans la fonction open_raster de la section raster.

  • Calculer la distance maximale (max_distance) entre deux pixels considérés comme pouvant être connectés*, en partant de l’hypothèse que la hauteur du pixel Δy et la largeur Δx sont les mêmes : img

  • Obtenez le trajectory des pixels qui ont un paramètre utilisateur défini pixel_value (par exemple, 1 pour tracer 1-pixels dans le binaire least cost.tif) et lancez une erreur si la trajectoire est vide (i.e., np.count_nonzero(trajectory) == 0).

  • Utilisez la fonction offset2coords ci-dessus pour ajouter les coordonnées des points à une liste points.

  • Créer un objet multi_line (installation de ogr.Geometry(ogr.wkbMultiLineString)), qui représente le chemin (éviter) final le moins cher.

  • Iterate through all possible combinations of points (excluding combinations of points with themselves) with itertools.combinations(iterable, r=number-of-combinations=2).

    • Les points sont stockés dans la liste points.

    • point1 et point2 sont nécessaires pour obtenir la distance entre les paires de points.

    • Si le distance entre les points est plus petit que max_distance, la fonction crée un objet en ligne à partir des deux points et l’ajoute à l’objet multi_line.

  • Créer un nouveau shapefile (nommé out_shp_fn) en utilisant la fonction create_shp() (avec vérification intégrée de la longueur du nom du shapefile à partir de flusstools geotools.create_shp()).

  • Ajouter l’objet multi_line comme nouvelle fonctionnalité au shapefile (selon les descriptions dans le shapefile section).

  • Créer un fichier de projection .prj (rappelez les descriptions dans le shapefile section) en utilisant le système de référence spatiale de l’entrée raster avec la fonction get_srs().

The raster2line function is also implemented in the flusstools.geotools.geotools script.

import os
import itertools
import numpy as np
from osgeo import ogr
from flusstools.geotools import (
    raster2array, open_raster, create_raster, create_shp, get_srs, make_prj,
)


def raster2line(raster_file_name, out_shp_fn, pixel_value):
    """
    Convert a raster to a line shapefile, where pixel_value determines line start and end points
    :param raster_file_name: STR of input raster file name, including directory; must end on ".tif"
    :param out_shp_fn: STR of target shapefile name, including directory; must end on ".shp"
    :param pixel_value: INT/FLOAT of a pixel value
    :return: None (writes new shapefile).
    """

    # calculate max. distance between points
    # ensures correct neighbourhoods for start and end pts of lines
    raster, array, geo_transform = raster2array(raster_file_name)
    pixel_width = geo_transform[1]
    max_distance = np.ceil(np.sqrt(2 * pixel_width**2))

    # extract pixels with the user-defined pixel value from the raster array
    trajectory = np.where(array == pixel_value)
    if np.count_nonzero(trajectory) == 0:
        print("ERROR: The defined pixel_value (%s) does not occur in the raster band." % str(pixel_value))
        return None

    # convert pixel offset to coordinates and append to nested list of points
    points = []
    count = 0
    for offset_y in trajectory[0]:
        offset_x = trajectory[1][count]
        points.append(offset2coords(geo_transform, offset_x, offset_y))
        count += 1

    # create multiline (write points dictionary to line geometry (wkbMultiLineString)
    multi_line = ogr.Geometry(ogr.wkbMultiLineString)
    for i in itertools.combinations(points, 2):
        point1 = ogr.Geometry(ogr.wkbPoint)
        point1.AddPoint(i[0][0], i[0][1])
        point2 = ogr.Geometry(ogr.wkbPoint)
        point2.AddPoint(i[1][0], i[1][1])

        distance = point1.Distance(point2)
        if distance < max_distance:
            line = ogr.Geometry(ogr.wkbLineString)
            line.AddPoint(i[0][0], i[0][1])
            line.AddPoint(i[1][0], i[1][1])
            multi_line.AddGeometry(line)

    # write multiline (wkbMultiLineString2shp) to shapefile
    new_shp = create_shp(out_shp_fn, layer_name="raster_pts", layer_type="line")
    lyr = new_shp.GetLayer()
    feature_def = lyr.GetLayerDefn()
    new_line_feat = ogr.Feature(feature_def)
    new_line_feat.SetGeometry(multi_line)
    lyr.CreateFeature(new_line_feat)

    # create projection file
    srs = get_srs(raster)
    make_prj(out_shp_fn, int(srs.GetAuthorityCode(None)))
    print("Success: Wrote %s" % str(out_shp_fn))

La fonction raster2line() peut être appelée comme suit pour convertir le chemin le moins cher du pixel (raster) au format ligne (vecteur) :

source_raster_fn = r"" +  os.path.abspath("") + "/geodata/rasters/least-cost.tif"
target_shp_fn = r"" + os.path.abspath("") + "/geodata/shapefiles/least-cost.shp"
pixel_value = 1
raster2line(source_raster_fn, target_shp_fn, pixel_value)
Success: Wrote /home/schwindt/jupyter/geodata/shapefiles/least_cost.shp
python convert raster to line

Figure 1:Le raster du chemin le moins coûteux converti en un fichier de forme de ligne, montré dans QGIS.

Raster à Polygon

gdal est livré avec la puissante fonction Polygonize pour la conversion facile d’un jeu de données raster en polygone shapefile. Alors que gdal.Polygonize permet d’écrire un simple raster2polygon() Fonction Python, il a l’inconvénient qu’il ne peut gérer que des valeurs entières et il attribue simplement aléatoirement FID (identificateurs attribués automatiquement sans signification physique) valeurs par défaut. Puisque les valeurs FID ne sont pas significatives, nous créerons une fonction float2int() pour préserver la plage de valeurs originale (utilise les fonctions raster2array() et create_raster() de la section raster):

def float2int(raster_file_name, band_number=1):
    """
    :param raster_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: new_raster_file_name (STR)
    """
    # use raster2array function to get raster, np.array and the geo transformation
    raster, array, geo_transform = raster2array(raster_file_name, band_number=band_number)
    
    # convert np.array to integers
    try:
        array = array.astype(int)
    except ValueError:
        print("ERROR: Invalid raster pixel values.")
        return raster_file_name
    
    # get spatial reference system
    src_srs = get_srs(raster)
    
    # create integer raster    
    new_name = raster_file_name.split(".tif")[0] + "_int.tif"
    create_raster(new_name, array, epsg=int(src_srs.GetAuthorityCode(None)),
                  rdtype=gdal.GDT_Int32, geo_info=geo_transform)
    # return name of integer raster
    return new_name

Ensuite, nous allons créer la fonction raster2polygon() qui implémente le workflow suivant:

  1. Utilisez la fonction float2int() pour vous assurer que tout raster file_name fourni peut être converti en valeurs purement entières.

  2. Créer un nouveau shapefile (nommé out_shp_fn) en utilisant la fonction create_shp() (également disponible à partir de flusstools: geotools.create_shp()).

  3. Ajouter un nouveau champ ogr.OFTInteger (rappeler how field creation works) dans la section shapefile) nommé par l’argument optionnel field_name en entrée.

  4. Run gdal.Polygonize with:

    • hSrcBand=raster_band

    • hMaskBand=None (bande de raster optionnelle pour définir les polygones)

    • hOutLayer=dst_layer

    • iPixValField=0 (si aucun champ n’a été ajouté, fixer à -1 afin de créer un champ FID; si d’autres champs ont été ajoutés, fixer à 1, 2, ... )

    • papszOptions=[] (aucun effet pour le type de pilote ESRI Shapefile)

    • callback=None pour ne pas utiliser l’algorithme de déclaration (GDALProgressFunc())

  5. Créer un fichier de projection .prj (rappelez les descriptions dans le shapefile section) en utilisant le système de référence spatiale de l’entrée raster avec la fonction get_srs().

def raster2polygon(file_name, out_shp_fn, band_number=1, field_name="values"):
    """
    Convert a raster to polygon
    :param file_name: STR of target file name, including directory; must end on ".tif"
    :param out_shp_fn: STR of a shapefile name (with directory e.g., "C:/temp/poly.shp")
    :param band_number: INT of the raster band number to open (default: 1)
    :param field_name: STR of the field where raster pixel values will be stored (default: "values")
    :return: None
    """
    # ensure that the input raster contains integer values only and open the input raster
    file_name = float2int(file_name)
    raster, raster_band = open_raster(file_name, band_number=band_number)

    # create new shapefile with the create_shp function
    new_shp = create_shp(out_shp_fn, layer_name="raster_data", layer_type="polygon")
    dst_layer = new_shp.GetLayer()

    # create new field to define values
    new_field = ogr.FieldDefn(field_name, ogr.OFTInteger)
    dst_layer.CreateField(new_field)

    # Polygonize(band, hMaskBand[optional]=None, destination lyr, field ID, papszOptions=[], callback=None)
    gdal.Polygonize(raster_band, None, dst_layer, 0, [], callback=None)

    # create projection file
    srs = get_srs(raster)
    make_prj(out_shp_fn, int(srs.GetAuthorityCode(None)))
    print("Success: Wrote %s" % str(out_shp_fn))

La fonction raster2polygon() peut être implémentée, par exemple, pour convertir le raster de profondeur d’eau pour 1000 CFS (h001000.tif de l’échantillon River Architect datasets) en polygone shapefile:

src_raster = r"" +  os.path.abspath("") + "/geodata/rasters/h001000.tif"
tar_shp = r"" + os.path.abspath("") + "/geodata/shapefiles/h_poly_cls.shp"
raster2polygon(src_raster, tar_shp)
Success: Wrote /home/schwindt/jupyter/geodata/shapefiles/h_poly_cls.shp
python convert raster to polygon shapefile

Figure 2:Le raster des profondeurs d’eau converti en des zones contenant un fichier polygone, montré dans QGIS.

Rasterize (Vector Shapefile to Raster)

Similar to gdal.Polygonize, gdal.RasterizeLayer represents a handy option to convert a shapefile into a raster. However, to be precise, a shapefile is not really converted into a raster but burned onto a raster. Thus, values stored in a field of a shapefile feature are used (burned) as pixel values for creating a new raster. Attention is required to ensure that the correct values and data types are used. To this end, the below shown rasterize() function implements the following workflow that avoids potential conversion headaches:

  1. Ouvrez le nom et la couche de shapefile fournis par l’utilisateur.

  2. Lisez l’étendue spatiale de la couche.

  3. Dériver la résolution x-y en fonction de l’étendue spatiale et d’une définition utilisateur pixel_size (mot-clé optionnel avec valeur par défaut).

  4. Créer un nouveau master GeoTIF en utilisant le

    • définition de l’utilisateur output_raster_file_name,

    • résolution x et y calculée, et

    • eType (par défaut est gdal.GDT_Float32 - rappeler toutes les options de type de données énumérées dans le raster section.

  5. Appliquer la géotransformation définie par l’étendue de la couche source et le pixel_size.

  6. Créez un raster band, remplissez le band avec le no_data_value (par défaut est -9999), et définissez le no_data_value.

  7. Définissez le système de référence spatiale du raster comme le fichier de forme source.

  8. Appliquer gdal.RasterizeLayer avec

    • dataset=target_ds (ensemble de données de raster ciblé),

    • bands=[1] (list(integer) - augmenter pour définir plus de bandes de raster et attribuer d’autres valeurs, par exemple, à partir d’autres champs du fichier de formes source),

    • layer=source_lyr (couche avec des fonctionnalités à graver au raster),

    • pfnTransformer=None (en savoir plus dans les documents gdal),

    • pTransformArg=None (en savoir plus dans les documents gdal),

    • burn_values=[0] (une valeur par défaut qui est brûlée dans le raster),

    • options=["ALL_TOUCHED=TRUE"] définit que tous les pixels touchés par un polygone obtiennent la valeur de champ du polygone - sinon définie: seuls les pixels qui sont entièrement dans le polygone obtiennent une valeur attribuée,

    • options=["ATTRIBUTE=" + str(kwargs.get("field_name"))] définit le nom du champ avec des valeurs à graver.

def rasterize(in_shp_file_name, out_raster_file_name, pixel_size=10, no_data_value=-9999,
              rdtype=gdal.GDT_Float32, **kwargs):
    """
    Converts any shapefile to a raster
    :param in_shp_file_name: STR of a shapefile name (with directory e.g., "C:/temp/poly.shp")
    :param out_raster_file_name: STR of target file name, including directory; must end on ".tif"
    :param pixel_size: INT of pixel size (default: 10)
    :param no_data_value: Numeric (INT/FLOAT) for no-data pixels (default: -9999)
    :param rdtype: gdal.GDALDataType raster data type - default=gdal.GDT_Float32 (32 bit floating point)
    :kwarg field_name: name of the shapefile's field with values to burn to the raster
    :return: produces the shapefile defined with in_shp_file_name
    """

    # open data source
    try:
        source_ds = ogr.Open(in_shp_file_name)
    except RuntimeError as e:
        print("Error: Could not open %s." % str(in_shp_file_name))
        return None
    source_lyr = source_ds.GetLayer()

    # read extent
    x_min, x_max, y_min, y_max = source_lyr.GetExtent()

    # get x and y resolution
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)

    # create destination data source (GeoTIff raster)
    target_ds = gdal.GetDriverByName('GTiff').Create(out_raster_file_name, x_res, y_res, 1, eType=rdtype)
    target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
    band = target_ds.GetRasterBand(1)
    band.Fill(no_data_value)
    band.SetNoDataValue(no_data_value)

    # get spatial reference system and assign to raster
    srs = get_srs(source_ds)
    try:
        srs.ImportFromEPSG(int(srs.GetAuthorityCode(None)))
    except RuntimeError as e:
        print(e)
        return None
    target_ds.SetProjection(srs.ExportToWkt())

    # RasterizeLayer(Dataset dataset, int bands, Layer layer, pfnTransformer=None, pTransformArg=None,
    # int burn_values=0, options=None, GDALProgressFunc callback=0, callback_data=None)
    gdal.RasterizeLayer(target_ds, [1], source_lyr, None, None, burn_values=[0],
                        options=["ALL_TOUCHED=TRUE", "ATTRIBUTE=" + str(kwargs.get("field_name"))])

    # release raster band
    band.FlushCache()

Enfin, la fonction rasterize() peut être appelée pour convertir le polygone polygone de profondeur d’eau shapefile /geodata/shapefiles/h poly cls.shp (téléchargez-le sous forme de fichier zip) en un raster (ceci est pratiquement inutile mais un exercice illustratif). Faites attention au type de données, qui est gdal.GDT_Int32 en combinaison avec l’argument field_name correctement défini.

src_shp = r"" + os.path.abspath("") + "/geodata/shapefiles/h_poly_cls.shp"
tar_ras = r"" +  os.path.abspath("") + "/geodata/rasters/h_re_rastered.tif"
rasterize(src_shp, tar_ras, pixel_size=5, rdtype=gdal.GDT_Int32, field_name="values")
python convert polygon to raster with rasterize

Figure 3:Le raster re-converti des profondeurs d’eau basé sur le polygone shapefile contenant des zones de profondeur, montré dans QGIS.

Géoanalyses
Manipulation des grilles
Géoanalyses
La bibliothèque commerciale arcpy