Geospatialanalysen mit ESRIs kommerzieller ArcGIS-Software und Arcpy-Paket. Dieses Notebook kann nicht in Jupyter ausgeführt werden und Codeblöcke erfordern die kommerzielle Arcpy-Bibliothek von ESRI.
Einleitung¶
Die ArcGIS Pro Software von ESRI verfügt über eine eigene Conda-Umgebung, die die Zusammenarbeit mit der kommerziellen arcpyBibliothek ermöglicht. arcpy ist entweder direkt unter ArcGIS Pro oder unter Verwendung von ESRIs conda-Umgebung zugänglich, die über ArcGIS Pro verwaltet werden kann. Da arcpy an die Nutzung eines kommerziellen Umfelds gebunden ist, ist es nur in jupyter Notebooks, die direkt unter esri.com (diese werden aber nicht mit den Beispielen in diesem Abschnitt funktionieren) zu erreichen. Diese Seite erklärt die grundsätzliche Nutzung von arcpy in extern IDEs (z.B. PyCharm, VS Code oder Spyder), die Integration von Lizenzen und die grundlegende Funktionalität von arcpy.
Hintergrund¶
Warum kann mit arcpy in Python arbeiten, obwohl es viele Lizenz- und Plattformbeschränkungen gibt?
Wenn Sie für ein Unternehmen mit Engineering-Services arbeiten, nutzt das Unternehmen wahrscheinlich ArcGIS wegen seiner Popularität und kommerziellen Support-Service. In diesem Fall umfasst ein Unternehmen (oder Forschungseinrichtung) die Lizenzgebühren und es ist ebenso wahrscheinlich, dass ein Windows-Betriebssystem aus ähnlichen Gründen verwendet wird. Unter den beliebten Funktionalitäten der ArcGIS-Software sind Raster Calculator, Shapefile Feature Manager oder Werkzeuge für die statistische Analyse von Geo-Datenbanken. Dank arcpy können solche beliebten ArcGIS-Tools in Python-Skripte eingebettet werden, wodurch Workflows automatisiert und Effizienz deutlich gesteigert werden können.
Diese Seite zeigt die Grundlagen für die Verarbeitung von Raster- und Formdateidaten in Python mit arcpy. Es gibt viele weitere Methoden, die in arcpy umgesetzt werden, und hier werden nur die Grundlagen erläutert.
Verwenden Sie arcpy mit externen IDEs¶
Setup Interpreter¶
Der Python-Installationsabschnitt erklärt how to set up PyCharm IDE mit einer conda-Umgebung. Um eine neue conda-Umgebung in PyCharm mit der conda-Umgebung von ESRI zu schaffen, muss die Location anders definiert werden (siehe die original screenshot):
ArcGIS Pro wurde weltweit vom Systemadministrator installiert, verwenden:
%PROGRAMFILES%\ArcGIS\Pro\bin\Python\Scripts\propyArcGIS Pro wurde für einzelne Benutzer installiert, verwenden:
%LOCALAPPDATA%\Programs\ArcGIS\Pro\bin\Python\Scripts\propy
Es kann vorkommen, dass es nicht notwendig ist, das propy-Verzeichnis hinzuzufügen. Darüber hinaus finden Sie die conda.exe oder python.exe in den obigen Verzeichnissen und definieren sie im Feld Conda ausführbar. Vervollständigen Sie die Schaffung der neuen Umgebung mit einem Klick auf Kreate.
Um weitere Pakete zu installieren, folgen Sie den Beschreibungen von Esri.
Import Arcpy und seine Module¶
arcpy kommt mit eigenen Klassen für geospatiale Datenobjekte (z.B. arcpy.Raster für netzgebundene Daten) und Module für Mapping (arcpy.mp), Emulation des Spatial Analyst (arcpy.sa) oder Zugriff auf Daten (arcpy.da). Eine vollständige Liste finden Sie auf der developer’s website. Diese Seite enthält Code-Blöcke mit arcpy und arcpy.sa, bei denen Spatial Analyst-Objekte mit * importiert werden, um einen direkten Zugriff auf z.B. Con() (anstelle von arcpy.sa.Con()) zu ermöglichen, was der bedingten Aussage von arcpy für Raster entspricht.
import arcpy
from arcpy.sa import *Setup arcpy Workspace¶
Geospatiale Berechnungen können viele Nebenprodukte produzieren, die schwer sein können. Um die Daten besser zu kontrollieren, die von arcpy generiert werden, und wo es eine gute Idee ist, einen Workspace für jedes Arcpy-Skript zu definieren:
arcpy.env.workspace = "C:\\workspace\\" # or use os.path.dirname(__file__) to go to the script directoryEs kann auch nützlich sein, das Überschreiben bereits bestehender Dateien mit dem gleichen Namen aufgrund der Dateigröße zu aktivieren. Dieses Verhalten kann oder darf nicht erwünscht sein und kann wie folgt gesteuert werden.
arcpy.gp.overwriteOutput = True # enable overwriting
arcpy.gp.overwriteOutput = False # disable overwritingSet Spatial Extents¶
Geospatial-Datensätze können große Ausmaße haben, ohne dass Daten in großen Teilen geschrieben wurden. Die Verarbeitung von Datenzellen kann zu unnötig langer Rechenzeit führen. Daher ist es ratsam, den Berechnungsumfang mit:
arcpy.env.extent = "MAXOF" # uses the combined extent of all input datasets
arcpy.env.extent = "MINOF" # uses only the overlap of all input datasets
arcpy.env.extent = arcpy.Extent(arcpy.Raster("base.tif")) # uses the controlled extent of a raster
arcpy.env.extent = "Xmin YMin XMax Ymax" # imposes user-defined minimum and maximum coordinatesCheckout-Lizenzen¶
Viele arcpy Methoden benötigen Lizenzen wie Spatial Analyst oder 3D. In eigenständigen Skripten können Lizenzen (checked out) mit arcpy.CheckOutExtension('NAME') und deaktiviert (checked in) mit arcpy.CheckInExtension('NAME') aktiviert werden. Angesichts der Objektorientierung sollten arcpy-Operationen in Funktionen oder Klassenmethoden eingebettet werden. Daher wird empfohlen, Funktionen oder Methoden mit arcpy mit decorators zu wickeln, die notwendige Lizenzen aktivieren. Der folgende Codeblock bietet eine Wrapper, um eine Spatial Analyst Lizenz für eine Funktion zu aktivieren.
def spatial_license(func):
def wrapper(*args, **kwargs):
arcpy.CheckOutExtension('Spatial')
result = func(*args, **kwargs)
arcpy.CheckInExtension('Spatial')
return result
return wrapperFehler verfolgen¶
Der Abschnitt Python Errors, Logging, and Debugging bietet nützliche Anweisungen für Fehlerbehebungen in Code- oder Codenutzung. Um Probleme in objektorientierten arcpyScripts zu identifizieren, wird eine zusätzliche Wrapper-Funktion empfohlen, die arcpyFehler an eine Logdatei (logger) schreibt. Die logger.*(MESSAGE)-Ausdrücke können auch durch print(MESSAGE) ersetzt werden.
import logging
def err_info(func):
def wrapper(*args, **kwargs):
arcpy.gp.overwriteOutput = True
logger = logging.getLogger("logfile")
try:
return func(*args, **kwargs)
except arcpy.ExecuteError:
logger.info(arcpy.GetMessages(2))
arcpy.AddError(arcpy.GetMessages(2))
except Exception as e:
logger.info(e.args[0])
arcpy.AddError(e.args[0])
except:
logger.info(arcpy.GetMessages())
return wrapperRaumanalyst und Rasteroperationen¶
Grundlagen¶
arcpy bietet verschiedene Optionen, um Geospatial zu laden rasters (gridded data), die verschiedene Formate wie Esri Grid haben kann (kein Ende ist der Raster ein Ordner mit anderen Dateien), GeoTIFF, DAT und vieles mehr. Das folgende Skript lädt eine Fließtiefe Grid rasterh und eine Strömungsgeschwindigkeit Grid raster u:
h = arcpy.Raster("geodata/input/rasters/h")
u = arcpy.Raster("geodata/input/rasters/u")Der Froude number kann aus der Strömungstiefe und Geschwindigkeit und der Gravitationskonstanten =9.81 m/s pixel berechnet werden. Das folgende Skript berechnet die Froude-Nummer für alle Pixel, bei denen die Durchflusstiefe mindestens 0,1 m beträgt. Der Rastervergleich wird durch die Spatial Analyst’s Con(if_condition, then, else)Zustandserklärung erreicht. Die beiden Raster (u und h) werden als arcpy.sa.Float()Objekte weitergegeben, um sicherzustellen, dass das Skript den richtigen Pixeldatentyp verwendet.
froude = Con(Float(h) > 0.1, Float(u) / SquareRoot(Float(h) * Float(9.81)))Besser als das 0,1-Meter-Kriterium ist die Berechnung der Froude-Zahl überall dort, wo die Strömungstiefe und die Geschwindigkeit einen Zahlenwert haben. arcpy.sa.IsNull() wertet aus, wo Pixel nicht numerisch sind. Wir interessieren uns jedoch für das Gegenteil (d.h. Pixel, die nicht numerisch sind), das wir durch das ~ (nicht)-Zeichen erhalten. Die unterfunktionalisierte Berechnung der Froude-Nummer nutzt die Möglichkeit, mehrere Con()-Ausdrücke zu verschachteln, um beide Raster (u und h) für numerische Pixel zu überprüfen.
Darüber hinaus benötigen wir eine Spatial Analyst Lizenz, um dieses Skript auszuführen. Daher ist es sinnvoll, den obigen Codeblock in eine Funktion umzuschreiben, die die spatial_licenseWrapper-Funktion aus dem oben genannten checkout Licenses-Bereich verwendet. Um den Code robust zu machen, fügen wir auch die above-defined *err_info*Wrapper-Funktion hinzu.
@err_info
@spatial_license
def calculate_froude(h, u):
return Con(~IsNull(h), Con(~IsNull(u), Float(u) / SquareRoot(Float(h) * Float(9.81))))
froude = calculate_froude(h, u)Erfahren Sie mehr über Raster Calculator und Map Algebra auf der developer’s website (esri).
Zellstatistik¶
Bei der Auswertung von numerischen Modelldaten werden oft statistische Werte (z.B. minimal oder maximal) eines oder mehrerer Raster berechnet. Der Vergleich mehrerer ähnlicher Raster ist z.B. dann sinnvoll, wenn derselbe Parameter mit zwei verschiedenen Modellen oder zu unterschiedlichen Zeitpunkten berechnet wurde (z.B. zur Beurteilung der morphodynamischen Flussentwicklung). arcpy.sa.CellStatistics([Raster1, Raster2, ... RasterN], TYPE, ...) ermöglicht solche statistischen Auswertungen. Der folgende Codeblock verdeutlicht den Vergleich der mit zwei verschiedenen hydrodynamischen numerischen Modellen berechneten Strömungsgeschwindigkeiten durch die Berechnung der MEAN (Durchschnitt) und Standardabweichung (STD).
u_basement = arcpy.Raster("geodata/bm/velocity.tif")
u_tuflow = arcpy.Raster("geodata/tf/velocity.dat")
u_mean = CellStatistics([u_basement, u_tuflow], "MEAN")
u_stdv = CellStatistics([u_basement, u_tuflow], "STD")Lesen Sie mehr Optionen Statistiktypen und Umgang mit nicht-numerischen Daten auf der developer’s website (Esri).
Shapefile Operationen¶
Das geospatial shapefile Vektorformat ist eine Esri-Erfindung. Kein Wunder, arcpy ist gut bei der Verarbeitung dieses Vektordatenformats. In der Hydrauliktechnik erstellen wir jedoch in der Regel (Zeichnen) Shapefiles manuell entweder direkt mit ArcGIS oder dessen Open-Source-Konkurrenten QGIS um z.B. bestimmte Strömungsregionen abzugrenzen. Beispiele sind in the BASEMENT tutorial (erweitern Sie die Erstellung von Elevationspunkt, Randpolygon und Breakline Polyline Shapefiles). In Codes wird die Verarbeitung von Shapefiles nur bei der Analyse der Ausgabe von numerischen Modellen wichtig (z.B. zur Klassifizierung von morphologischen Einheitsmerkmalen, zur exakten Berechnung von Patch-Bereichen oder zur automatischen Versteifung in Bauplänen). In diesem Stadium müssen Rasterdaten (Ausgang von numerischen Modellen) zunächst in Formdateien umgewandelt werden. Aus diesem Grund beginnt dieses Tutorial mit der Umwandlung von Rasterdaten in Formdateien zusammen mit der Darstellung anderer Funktionen wie der Berechnung von Patch-Bereich und dem Zugriff auf formfile Attributtabellen.
Raster können in Polygon und andere Shapefile-Typen (z.B. Punkt) umgewandelt werden. Das folgende Beispiel zeigt die Umwandlung eines Rasters in eine Polygonformdatei. Es verwendet einen ganzzahligen Raster aller Pixel, wobei die Strömungstiefe und -geschwindigkeit kleiner als 1,4 m bzw. 0,15 m/s sind. Diese flachen und langsam fließenden Regionen werden als slackwater (nach Wyrick & Pasternack (2014)) bezeichnet. Da Slackwater als bevorzugter Lebensraum einiger Wasserarten bezeichnet wird, fragen wir uns jetzt, wie viel Slackwater-Bereich das numerische Modell im simulierten Flussabschnitt vorhersagt. Dazu wandeln wir den Slackwater-Raster in eine Formdatei um und berechnen die Oberfläche der Formdatei mit den folgenden arcpy Methoden:
Convert the raster to a shapefile with
arcpy.RasterToPolygon_conversion()with arguments:in_rasterist einarcpy.Raster()Objekt von integer*-Werten (mitFloatraster ergibt sich ein Fehler!).out_polygon_featuresist ein String des Ausgabedateinamens und -verzeichnisses.simplifyist ein optionaler string, der entweder"NO_SIMPLIFY"sein kann, um die exakte Zeichnung von Polygongrenzen entlang der Pixelgrenze zu zwingen, oder"SIMPLIFY", um Polygongrenzen über Pixel zu ermöglichen.
Add a new field to the new polygon shapefile with
arcpy.AddField_management()with arguments:field_namekann alle string ohne Leerzeichen sein.field_typeist ein string, der definiert, ob das Feld numerisch ist (z.B."FLOAT"oder"LONG"für integer), Datum/Zeit ("DATE"),"TEXT"oder"RASTER".field_precisionist ein (optional) integer (Long), der die Anzahl der im neuen Feld gespeicherten Ziffern festlegt.Weitere optionale Argumente können festgelegt werden, um die Anzahl der Dezimals oder Zeichen, einen alternativen Feldnamen, aktivieren Sie
NULL, eine Felddomäne oder wenn ein Feld benötigt wird.
Calculate patch area with
arcpy.CalculateGeometryAttributes_management()with arguments:in_featuresist ein string, das das Verzeichnis und den Namen einer Feature-Schicht definiert.geometry_propertyist eine geschachtelte Liste von[[Target-field-name, Property], [Another-target-field-name, Another-property], ...]zur Berechnung geometrischer Eigenschaften wie"AREA","HOLE_COUNT"oder"PART_COUNT"(und vieles mehr).area_unitkann"SQUARE_METERS"oder"SQUARE_KILOMETERS"sein (und viele andere Optionen für die üblichen US-Einheiten).Weitere optionale Argumente können festgelegt werden, um Längeneinheiten (für Perimeterbeurteilungen) oder ein Koordinatensystem und Format zu definieren.
Der untenstehende Codeblock enthält zusätzlich die Anwendung dieser Methoden und zeigt auch, wie der Bereichswert aus der Attributtabelle einer Shapefile mit arcpy.da.UpdateCursor(shapefile-name, field-name) reihenweise gelesen werden kann.
# create a slackwater raster that is arcpy.sa.Int(1) where h and u criteria are true and NULL elsewhere
slackwater = Con((Float(h) <= 1.4) & (Float(u) <= 0.15), Int(1))
# define directory and name of the new shapefile
new_shp_file = "geodata/shapefiles/slackwater.shp"
# convert slackwater raster to polygon shapefile
arcpy.RasterToPolygon_conversion(in_raster=slackwater, out_polygon_features=new_shp_file, simplify="NO_SIMPLIFY")
# add a new field to the new shapefile's attribute table (name)
arcpy.AddField_management(new_shp_file, field_name="F_AREA", field_type="FLOAT", field_precision=9)
# calculate area of all polygons in attribute table
area_unit="SQUARE_METERS"
arcpy.CalculateGeometryAttributes_management(in_features=new_shp_file, geometry_property=[["F_AREA", "AREA"]], area_unit=area_unit)
area = 0.0
with arcpy.da.UpdateCursor(new_shp_file, "F_AREA") as cursor:
for row in cursor:
try:
area += float(row[0])
except ValueError:
print("WARNING: Patch with invalid area value (%s)." % str(row))
print("Sum of all patches = {0} {1}".format(str(area), area_unit))- Wyrick, J. R., & Pasternack, G. B. (2014). Geospatial organization of fluvial landforms in a gravel-cobble river: Beyond the riffle-pool couplet. Geomorphology, 213(Supplement C), 48–65. 10.1016/j.geomorph.2013.12.040