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.

Vectorize und Rasterize

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

Importrelevante Bibliotheken

Importieren Sie die relevanten Pakete für die Handhabung von Rastern, Shapefiles und Georeferenzen:

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

Vektor

Raster nach Line

In diesem Abschnitt konvertieren wir den least cost path Raster-Datensatz (least cost.tif) in eine (Poly-)Zeilen-Shapefile. Zu diesem Zweck schreiben wir zunächst eine Funktion namens offset2coords(), die die Inverse der coords2offset()-Funktion darstellt, und konvertiert den x/y-Offset (in * ganzzahligen* Pixelzahlen) in Koordinaten der Geotransformation eines Geodatensatzes:

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

Next, we can write a core function to convert a raster dataset to a line shapefile. We name this function raster2line() and it builds on the following workflow:

  • Öffnen Sie ein raster, dessen Band als array und geo_transform (Geo-Transformation) definiert mit dem raster_file_name-Argument in der open_raster-Funktion aus dem Rasterbereich.

  • Berechnen Sie den maximalen Abstand (max_distance) zwischen zwei Pixeln, die als * verbindungsfähig * betrachtet werden, basierend auf der Hypothese, dass die Pixelhöhe * & Delta; y * und die Breite * & Delta; x * gleich sind.

  • Get the trajectory of pixels that have a user parameter-defined pixel_value (e.g., 1 to trace 1-pixels in the binary least_cost.tif) and throw an error if the trajectory is empty (i.e., np.count_nonzero(trajectory) == 0).

  • Use the above-defined offset2coords function to append point coordinates to a points list.

  • Erstellen Sie ein multi_line-Objekt (Instanz von ogr.Geometry(ogr.wkbMultiLineString)), das den (void) letzten kostengünstigsten Pfad darstellt.

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

    • Points are stored in the points list.

    • point1 und point2 sind erforderlich, um den Abstand zwischen Punktenpaaren zu erhalten.

    • If the distance between the points is smaller than max_distance, the function creates a line object from the two points and appends it to the multi_line object.

  • Erstellen Sie eine neue Shapefile (genannt out_shp_fn) mit der Funktion create_shp() (mit integrierter Shapefile-Namenslängenüberprüfung von flusstools geotools.create_shp()).

  • Fügen Sie das multi_line-Objekt als neues Feature zur Shapedatei hinzu (gemäß den Beschreibungen im shapefile section).

  • Create a .prj projection file (recall descriptions in the shapefile section) using the spatial reference system of the input raster with the get_srs() function.

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

import itertools
import numpy as np
from flusstools.geotools import raster2array
from flusstools.geotools import create_shp
from flusstools.geotools import get_srs
from flusstools.geotools import 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 = abs(geo_transform[1])
    pixel_height = abs(geo_transform[5])
    max_distance = np.hypot(pixel_width, pixel_height) * 1.001

    # extract pixels with the user-defined pixel value from the raster array
    trajectory = np.where(array == pixel_value)
    if trajectory[0].size == 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)))

    # close the output shapefile and release the input raster
    new_line_feat = None
    lyr = None
    new_shp = None
    raster = None
    print("Success: Wrote %s" % str(out_shp_fn))

Die raster2line()-Funktion kann wie folgt aufgerufen werden, um den kostengünstigsten Pfad vom Pixel- (Raster-) in das Zeilen- (Vektor-)Format zu konvertieren:

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:Das Raster des kostengünstigsten Pfades, der in eine Linienformdatei umgewandelt wurde, wird in QGIS angezeigt.

Raster nach Polygon

gdal.Polygonize converts connected raster regions to polygons. It reads source values through an integer buffer, so floating-point values are truncated unless they are scaled and converted first; gdal.FPolygonize is the floating-point alternative. The float2int() helper below applies an optional scale factor before rounding so decimal precision can be retained in the stored integer classes. It uses the raster2array() and create_raster() functions from the raster section:

from flusstools.geotools import *
def float2int(raster_file_name, band_number=1, scale_factor=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)
    :param scale_factor: positive numeric factor applied before rounding (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)
    
    # scale, round, and convert valid pixels to integers
    try:
        if scale_factor <= 0:
            raise ValueError("scale_factor must be positive")
        array = np.where(np.isnan(array), -9999, np.rint(array * scale_factor)).astype(np.int32)
    except (TypeError, 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

Next, we will create the raster2polygon() function that implements the following workflow:

  1. Use the float2int() function to ensure that any raster file_name provided can be converted to purely integer values.

  2. Erstellen Sie eine neue Shapefile (genannt out_shp_fn) mit der Funktion create_shp() (auch bei flusstools verfügbar: geotools.create_shp()).

  3. Add a new ogr.OFTInteger field (recall how field creation works) in the shapefile section) named by the optional field_name input argument.

  4. Run gdal.Polygonize with:

    • hSrcBand=raster_band

    • hMaskBand=None (optionales Rasterband zur Definition von Polygonen)

    • hOutLayer=dst_layer

    • iPixValField=0 (wenn kein Feld hinzugefügt wurde, setzen Sie auf -1, um ein FID-Feld zu erstellen; wenn mehr Felder hinzugefügt wurden, setzen Sie auf 1, 2, ...)

    • papszOptions=[] (kein Effekt für ESRI Shapefile Treibertyp)

    • callback=None für die Nichtverwendung des Berichtsalgorithmus (GDALProgressFunc())

  5. Create a .prj projection file (recall descriptions in the shapefile section) using the spatial reference system of the input raster with the get_srs() function.

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

Die Funktion raster2polygon() kann zum Beispiel implementiert werden, um das Wassertiefenraster für 1000 CFS (h001000.tif aus den River Architect Beispieldatensätzen] in ein Polygon-Shapefile umzuwandeln:

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:Das Raster der Wassertiefen wurde in eine Polygonformdatei mit Zonen umgewandelt, dargestellt in QGIS.

Rasterize (Vector Shapefile zu 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. Öffnen Sie den vom Benutzer bereitgestellten Shapefile-Namen und -Layer.

  2. Lesen Sie die räumliche Ausdehnung der Schicht.

  3. Ableiten der x-y-Auflösung als Funktion der räumlichen Ausdehnung und eines benutzerdefinierten pixel_size (optionales Keyword-Argument mit Standardwert).

  4. Erstellen Sie ein neues GeoTIFF-Raster mit dem

    • benutzerdefinierte output_raster_file_name,

    • berechnete x- und y-Auflösung und

    • eType (Standard ist gdal.GDT_Float32 - erinnern Sie sich an alle Datentypoptionen, die unter raster section aufgeführt sind.)

  5. Apply the geotransformation defined by the source layer extents and the pixel_size.

  6. Create one raster band, fill the band with the user-defined no_data_value (default is -9999), and set the no_data_value.

  7. Stellen Sie das räumliche Bezugssystem des Rasters mit der Quellformdatei ein.

  8. Apply gdal.RasterizeLayer with

    • dataset=target_ds (Zielrasterdatensatz),

    • bands=[1] (list(integer) - erhöhen Sie auf definierte mehr Rasterbänder und weisen Sie andere Werte zu, zum Beispiel aus anderen Feldern der Quellformdatei),

    • layer=source_lyr (Schicht mit Funktionen zum Brennen im Raster),

    • pfnTransformer=None (lesen Sie mehr im gdal docs]),

    • pTransformArg=None (lesen Sie mehr im gdal docs]),

    • burn_values=[0] (a default value that is burned to the raster),

    • options=["ALL_TOUCHED=TRUE"] definiert, dass alle Pixel, die von einem Polygon berührt werden, den Feldwert des Polygons erhalten - wenn nicht gesetzt: Nur Pixel, die sich vollständig im Polygon befinden, erhalten einen zugewiesenen Wert,

    • options=["ATTRIBUTE=" + str(kwargs.get("field_name"))] defines the field name with values to burn.

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: None (writes the raster defined by out_raster_file_name)
    """

    # open data source
    source_ds = ogr.Open(in_shp_file_name)
    if source_ds is None:
        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(np.ceil((x_max - x_min) / pixel_size))
    y_res = int(np.ceil((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 = source_lyr.GetSpatialRef()
    if srs is None:
        raise ValueError("The input layer has no spatial reference.")
    target_ds.SetProjection(srs.ExportToWkt())

    field_name = kwargs.get("field_name")
    if not field_name:
        raise ValueError("field_name is required when burning attribute values.")

    gdal.RasterizeLayer(
        target_ds, [1], source_lyr,
        options=["ALL_TOUCHED=TRUE", "ATTRIBUTE=" + str(field_name)],
    )

    # flush and close datasets
    band.FlushCache()
    band = None
    target_ds = None
    source_lyr = None
    source_ds = None

Finally, the rasterize() function can be called to convert the polygonized water depth polygon shapefile /geodata/shapefiles/h_poly_cls.shp (download it as a zip file) back to a raster (this is practically useless but an illustrative exercise). Pay attention to the data type, which is gdal.GDT_Int32 in combination with the correctly defined field_name argument.

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:Das rekonvertierte Raster der Wassertiefen basierend auf dem Polygon-Shapefile mit Tiefenzonen, gezeigt in QGIS.

Geodatenanalysen
Raster (Grid) Handling
Geodatenanalysen
Die kommerzielle arcpy Bibliothek