Rasterverarbeitung mit Python und ArcGIS

Seit dem Sommer ist ArcGIS 10 verfügbar und bei den meisten unserer Mitarbeitenden schon im Einsatz. Die neue Version bringt einige Änderungen mit sich – unser Eindruck ist bisher positiv. Ein Beispiel: Im Bereich des Geoprocessing mit Python wird arcgisscripting durch die ArcPy-Site Package abgelöst. Für mich als Spezialist im Bereich der Rasterprozessierung und Geländemodellierung ist dies sehr interessant: Erstmals kann man mit Python aus ArcGIS heraus einzelne Rasterzellen und deren Nachbarschaft ansprechen und Berechnungen anstellen. Damit ist das Scripting mit Python in ArcGIS im Bereich Raster deutlich verbessert worden; der Entwicklung komplexerer Rasterfunktionen steht nichts mehr im Weg.

Konkret erlaubt ArcPy die Umwandlung eines ESRI Grids in einen Array, der im Python-Package NumPy verwendet werden kann. In NumPy kann man einfach Operatoren erstellen, die auf einer gewissen Nachbarschaft um eine Rasterzelle wirken. Nach abgeschlossener Berechnung konvertiert ArcPy den NumPy-Array wieder zurück in ein ESRI Grid und speichert dieses ab. Als einfaches Beispiel habe ich einen Laplace-Filter zur Detektion von Kanten implementiert (angenommen ist ein Raster ohne NoData-Werte):

import arcpy
from arcpy.sa import *
import numpy

# Determine ESRI grid properties
ras = Raster("raster")
nd = ras.noDataValue
x_cellsize = ras.meanCellWidth
y_cellsize = ras.meanCellHeight
lowerleft = ras.extent.lowerLeft

# Convert ESRI grid to NumPy array
nparray = arcpy.RasterToNumPyArray(ras, "", "", "", numpy.nan)

# Iterate through array and apply filter
nparray[1:-1, 1:-1] = 4*nparray[1:-1, 1:-1] - nparray[0:-2, 1:-1] \
                      - nparray[2:, 1:-1] - nparray[1:-1, 0:-2] \
                      - nparray[1:-1, 2:]
nparray[0:,0] = nd
nparray[0:,-1] = nd
nparray[0,0:] = nd
nparray[-1,0:] = nd

# Save NumPy array to ESRI grid with appropriate properties
raster_edges = arcpy.NumPyArrayToRaster(nparray, lowerleft,
                                        x_cellsize, y_cellsize, numpy.nan)
raster_edges.save("C:/temp/raster_edges")

Dieser Code kann direkt ins Python-Fenster in ArcMap kopiert und ausgeführt werden. Das Resultat zeigen die folgenden Abbildungen:

Das Original-Raster

Kanten-Raster

Ralph Straumann

Ralph Straumann

Ralph Straumann (Dr. sc. nat.) hat an der Universität Zürich Geographie mit Vertiefung in GIS, Wirtschaftsgeographie und Politologie studiert.

Seit 2010 arbeitet er im Tätigkeitsfeld Systemberatung + Analytik von EBP Informatik als Senior Consultant.

Er berät Kunden bei strategischen Fragen, zu Geschäftsprozessen und Organisation sowie bezüglich Quellen, Modellierung, Workflows und Analyse mit verschiedenartigen Daten im Schnittbereich zwischen IT/GIS und Anwendungsfeldern wie Verkehr und Raumplanung.

Mail: ralph.straumann@ebp.ch

Ralph Straumann auf:

Das könnte dich auch interessieren...