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