Werkzeuge zur Rasterverarbeitung

Dieses Beispiel zeigt einige der Raster- und Bildverarbeitungstools, die KDRM zur Automatisierung von Raster-Workflows erstellt hat. Diese Tools sind nützlich für Aufgaben wie das Batch-Zuschneiden von Rastern auf einen interessierenden Bereich, das Resampling von Rastern auf eine erforderliche Ausgabeauflösung, das Vorbereiten von Rastereingaben für eine spätere Analyse und das Entfernen langsamer, sich wiederholender manueller Verarbeitungsschritte aus größeren Lieferabläufen.

Die folgenden Beispiele konzentrieren sich nur auf Open-Source-Python-Pakete. In der Praxis können Tools wie dieses als eigenständige Skripte, geplante Verarbeitungsaufträge oder eingebettet in umfassendere Desktop-Workflows in ArcGIS Pro oder QGIS ausgeführt werden, abhängig von der Client-Umgebung und der Art und Weise, wie der Workflow ausgeführt werden muss.

Wofür werden diese Tools verwendet?

• Batch-Zuschneiden von DEMs, DTMs und anderen Rastern auf einen interessierenden Bereich.
• Resampling von Rastern, um Zellgrößen für spätere Berechnungen und Modellierungen auszurichten.
• Filtern und Validieren von Rastereingaben vor der nachgelagerten Analyse.
• Automatisierung von Rasterverarbeitungsschritten, die sonst stundenlange manuelle Arbeit erfordern würden.
• Erstellen wiederholbarer Rasterverarbeitungspipelines, die in Desktop-, Server- oder Befehlszeilenumgebungen ausgeführt werden können.

Warum wir solche Tools bauen

Wir erstellen regelmäßig benutzerdefinierte Automatisierungen für rasterlastige Projekte, da Leistung, Wiederholbarkeit und Portabilität wichtig sind. Manchmal besteht die Anforderung darin, als reiner Open-Source-Workflow ohne ArcGIS-Abhängigkeit ausgeführt zu werden. In anderen Fällen muss dieselbe Logik in einem ArcGIS Pro- oder QGIS-Workflow integriert sein, damit Analysten den Prozess in einer vertrauten Desktop-Umgebung ausführen können. Der zugrunde liegende Ansatz ist derselbe: Automatisieren Sie die sich wiederholenden Teile, halten Sie die Logik transparent und erleichtern Sie die Wiederholung und Unterstützung des Workflows.

Das folgende Clipping-Beispiel wurde für einen Arbeitsablauf erstellt, bei dem das Batch-Clipping von Rastern mit einem herkömmlichen Desktop-Ansatz zu langsam war. Durch das Umschreiben des Prozesses mit Open-Source-Python-Paketen wurde der Arbeitsablauf erheblich schneller und praktischer für größere Rastersammlungen. Das Resampling-Beispiel wurde aus einem ähnlichen Grund erstellt, bei dem große Rasterdatensätze auf eine feinere Ausgabezellengröße ausgerichtet werden mussten, damit spätere Rasterberechnungen mit der erforderlichen Auflösung ausgeführt werden konnten.

Beispiel 1. Batch-Clipping von Rastern mit Open-Source-Python

import json
import numpy as np
import os
import rasterio
import shutil

from rasterio.mask import mask
from shapely.geometry import shape
# the input and output directories for the rasters
input_directory = "input"
output_directory = "output"

# delete and recreate the output directory
if(os.path.exists(output_directory)):
    shutil.rmtree(output_directory)

os.makedirs(output_directory)

feature_collection = dict()

with open('national_park.geojson', 'r') as f:
    feature_collection = json.load(f)
area_of_interest = feature_collection['features'][0]['geometry']
clip_area = shape(area_of_interest)

for file in os.listdir(input_directory):
    # filter by tif files and only files with 2019 in the name
    if file.endswith('.tif') and '2019' in file:
        print(file)

        input_raster = os.path.join(input_directory, file)
        output_raster=os.path.join(output_directory, file)
        with rasterio.open(input_raster) as src:
            try:
                out_image, out_transform = mask(src, [clip_area], crop=True)

                # Set all nodata values to -99999
                out_image[out_image == src.nodata] = -99999
                raster_intersected = not np.all(out_image == -99999)
                # Check if the clip area intersects the raster and only save the raster if the clip area intersects with it
                if raster_intersected:
                    out_meta = src.meta.copy()
                    out_meta.update(
                        {
                            "driver": "GTiff",
                            "height": out_image.shape[1],
                            "width": out_image.shape[2],
                            "transform": out_transform,
                        }
                    )
                    with rasterio.open(output_raster, "w", **out_meta) as dest:
                        dest.write(out_image)

            except Exception as err:
                if ('Input shapes do not overlap raster' not in str(err)):
                    print('### Error with {} ###'.format(input_raster))
                    print(err)
                pass

Beispiel 2. Resampling von Rastern mit Open-Source-Python

import os
from osgeo import gdal

import numpy as np
import os

# the input and output directories
input_directory = 'input/DTM'
output_directory = 'output/DTM'

if not os.path.exists(output_directory):
    os.makedirs(output_directory)
for root, dirs, files in os.walk(input_directory):
    for file in files:
        if(file.endswith('.tif')):
            input_file = os.path.join(root, file)
            output_file = os.path.join(output_directory, file)

            print(file)

            # Read the input raster file
            raster = gdal.Open(input_file)
            # Get the original cell size
            cell_size = raster.GetGeoTransform()[1]
            # print('Original cell size: {}'.format(cell_size))

            # Calculate the number of cells for the output
            num_cells = int(np.ceil(cell_size / 0.5))
            # Create the output raster file
            driver = gdal.GetDriverByName('GTiff')
            out_raster = driver.Create(output_file, raster.RasterXSize * num_cells,
                                    raster.RasterYSize * num_cells, 1,
                                    gdal.GDT_Float32)

            # Set the output raster's projection
            out_raster.SetProjection(raster.GetProjection())
            # Set the output raster's geotransform
            out_raster_geo = list(raster.GetGeoTransform())
            out_raster_geo[1] = 0.5
            out_raster_geo[5] = -0.5
            out_raster.SetGeoTransform(out_raster_geo)

            nodata = -99999

            for i in range(1, out_raster.RasterCount + 1):
                # set the nodata value of the band
                out_raster.GetRasterBand(i).SetNoDataValue(nodata)
            # Copy the input raster data to the output raster
            gdal.ReprojectImage(raster, out_raster,
                                raster.GetProjection(),
                                out_raster.GetProjection(),
                                gdal.GRA_NearestNeighbour)

            # Close the files
            out_raster = None
            raster = None

Lieferoptionen

Die obigen Beispiele sind bewusst einfach gehalten, aber der gleiche Ansatz kann auf größere Produktionsabläufe mit Protokollierung, Validierung, Parallelverarbeitung, Metadatenprüfungen, Paketierung und Integration in umfassendere GIS-Bereitstellungspipelines erweitert werden. Je nach Kundenanforderung können wir diese Tools als eigenständige Skripte, in den Desktop integrierte Tools für ArcGIS Pro oder QGIS oder als produktive Dienste als Teil einer umfassenderen Geodatenplattform bereitstellen.