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.

Die kommerzielle arcpy Bibliothek

Geospatial-Analysen mit der kommerziellen ArcGIS-Software und dem Arcpy-Paket von ESRI. Dieses Notizbuch kann nicht in Jupyter ausgeführt werden und Codeblöcke erfordern die kommerzielle arcpy-Bibliothek von ESRI.

Einleitung

ESRIs ArcGIS Pro-Software verfügt über eine eigene Conda-Umgebung, die die Arbeit mit der kommerziellen arcpy-Bibliothek ermöglicht. arcpy ist entweder direkt unter ArcGIS Pro oder über die Conda-Umgebung von ESRI zugänglich, die über ArcGIS Pro verwaltet werden kann. Da arcpy an die Nutzung einer kommerziellen Umgebung gebunden ist, kann nur auf Jupyter-Notebooks zugegriffen werden, die direkt unter esri.com gehostet werden (diese funktionieren jedoch nicht mit den Beispielen in diesem Abschnitt). Diese Seite erklärt die grundlegende Verwendung von arcpy in externen IDEs (z. B. PyCharm, VS Code oder Spyder), die Integration von Lizenzen und die grundlegende Funktionalität von arcpy.

Hintergrund

***Warum kann die Arbeit mit arcpy in Python immer noch ein mächtiges Werkzeug sein, obwohl es viele Lizenz- und Plattformbeschränkungen gibt? ***

Wenn Sie für ein Unternehmen mit Engineering-Dienstleistungen arbeiten, verwendet das Unternehmen wahrscheinlich ArcGIS wegen seiner Popularität und kommerziellen Support-Service. In diesem Fall übernimmt ein Unternehmen (oder eine Forschungseinrichtung) die Lizenzgebühren und es ist ebenso wahrscheinlich, dass ein Windows-Betriebssystem aus ähnlichen Gründen verwendet wird. Zu den beliebten Funktionalitäten der ArcGIS-Software gehören Raster Calculator, Shapefile Feature Manager oder Tools zur statistischen Analyse von Geo-Datenbanken. Dank arcpy können solche beliebten ArcGIS-Tools in Python-Skripte eingebettet werden, wodurch Workflows automatisiert und die Effizienz deutlich gesteigert werden können.

This page shows the basics for manipulating raster and shapefile data in Python with arcpy. There are many more methods implemented in arcpy and only the fundamentals are explained here.

Use arcpy with External IDEs

Setup Interpreter

The Python installation section explains how to set up PyCharm IDE with a conda environment. To create a new conda environment in PyCharm with ESRI’s conda environment, the Location must be defined differently (see the original screenshot):

  • Wenn ArcGIS Pro wurde global vom Systemadministrator installiert, verwenden Sie: %PROGRAMFILES%\ArcGIS\Pro\bin\Python\Scripts\propy

  • Wenn ArcGIS Pro wurde für einzelne Benutzer installiert, verwenden Sie: %LOCALAPPDATA%\Programs\ArcGIS\Pro\bin\Python\Scripts\propy

Es kann vorkommen, dass es nicht notwendig ist, das propy Verzeichnis hinzuzufügen. Suchen Sie darüber hinaus das conda.exe oder python.exe in den obigen Verzeichnissen und definieren Sie es im Feld Conda executable. Beenden Sie die Erstellung der neuen Umgebung mit einem Klick auf Create.

Um weitere Pakete zu installieren, folgen Sie den Beschreibungen von Esri.]

Import arcpy und seine Module

arcpy verfügt über eigene Klassen für Geodatenobjekte (z. B. arcpy.Raster für gerasterte 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 Website des Entwicklers]. Diese Seite enthält Codeblöcke unter arcpy und arcpy.sa, wobei Spatial Analyst-Objekte unter * importiert werden, um beispielsweise direkten Zugriff auf Con() (anstelle von arcpy.sa.Con()) zu ermöglichen, was der bedingte Anweisung von arcpy für Raster entspricht.

import arcpy
from arcpy.sa import *

Arcpy Workspace einrichten

Geospatial calculations may produce many side products, which can be heavy in size. To better control which data is generated by arcpy and where it is a good idea to define a workspace for each arcpy script with:

arcpy.env.workspace = "C:\\workspace\\"  # or use os.path.dirname(__file__) to go to the script directory

Es kann auch nützlich sein, das Überschreiben bereits vorhandener Dateien mit dem gleichen Namen aufgrund der Dateigröße zu ermöglichen. Dieses Verhalten kann oder auch nicht gewünscht sein und kann wie folgt gesteuert werden.

arcpy.gp.overwriteOutput = True   # enable overwriting
arcpy.gp.overwriteOutput = False  # disable overwriting

Räumliche Ausdehnungen festlegen

Geodatensätze können große Ausdehnungen haben, ohne dass Daten in große Teile geschrieben werden. Die Verarbeitung von leeren Datenzellen kann zu einer unnötig langen Rechenzeit führen. Daher ist es ratsam, den Berechnungsumfang zu begrenzen 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 coordinates

Checkout-Lizenzen

Many arcpy methods need licenses such as Spatial Analyst or 3D. In stand-alone scripts, licenses can be activated (checked out) with arcpy.CheckOutExtension('NAME') and deactivated (checked in) with arcpy.CheckInExtension('NAME'). In light of the object orientation, arcpy operations should be embedded in functions or methods of classes. Therefore, it is recommended to wrap functions or methods using arcpy with decorators that activate necessary licenses. The following code block provides a wrapper to activate a Spatial Analyst license for a function.

def spatial_license(func):
    def wrapper(*args, **kwargs):
        arcpy.CheckOutExtension('Spatial')
        result = func(*args, **kwargs)
        arcpy.CheckInExtension('Spatial')
        return result
    return wrapper

Gleisfehler

The Python Errors, Logging, and Debugging section provides useful instructions for troubleshooting errors in code or code usage. To identify problems in object-oriented arcpy scripts, an additional wrapper function is recommended, which writes arcpy errors to a log file (logger). The logger.*(MESSAGE) expressions can also be replaced with print(MESSAGE).

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 wrapper

Spatial Analyst und Raster Operations

Grundlagen

arcpy bietet verschiedene Optionen zum Laden von Geodaten rasters (gridded data), die verschiedene Formate haben können, wie Esri Grid (keine Endung, das Raster ist ein Ordner mit anderen Dateien), GeoTIFF, DAT und viele mehr. Das folgende Skript lädt ein Flusstiefenraster h und ein Flussgeschwindigkeitsraster u:

h = arcpy.Raster("geodata/input/rasters/h")
u = arcpy.Raster("geodata/input/rasters/u")

The Froude-Zahl can be calculated pixel by pixel from the flow depth and velocity and the gravity constant gg=9.81 m/s2^2. The following script calculates the Froude number for all pixels where the flow depth is at least 0.1 m. The raster comparison is achieved through the Spatial Analyst’s Con(if_condition, then, else) conditional statement. The two rasters (u and h) are passed as arcpy.sa.Float() objects to ensure that the script uses the correct pixel data type.

froude = Con(Float(h) > 0.1, Float(u) / SquareRoot(Float(h) * Float(9.81)))

Better than the 0.1-meter criterion is to calculate the Froude number everywhere where the flow depth and the velocity have a numerical value. arcpy.sa.IsNull() evaluates where pixels are non-numeric. However, we are interested in the opposite (i.e., pixels that are not non-numeric), which we get through the ~ (not) sign. The below-functionalized calculation of the Froude number takes advantage of the possibility to nest multiple Con() expressions to check both rasters (u and h) for numeric pixels. In addition, we need a Spatial Analyst license to run this script. Therefore, it makes sense to rewrite the above code block into a function that uses the spatial_license wrapper function from the above checkout Licenses section. To make the code robust, we also add the above-defined *err_info* wrapper function.

@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 Website des Entwicklers (esri).

Zellstatistik

Bei der Auswertung numerischer Modelldaten werden häufig statistische Werte (z. B. Minimum oder Maximum) eines oder mehrerer Raster berechnet. Der Vergleich mehrerer ähnlicher Raster ist beispielsweise dann sinnvoll, wenn derselbe Parameter mit zwei verschiedenen Modellen oder zu unterschiedlichen Zeitpunkten berechnet wurde (z. B. um die morphodynamische Entwicklung von Flüssen zu beurteilen). arcpy.sa.CellStatistics([Raster1, Raster2, ... RasterN], TYPE, ...) ermöglicht solche statistischen Auswertungen. Der folgende Codeblock veranschaulicht den Vergleich von Strömungsgeschwindigkeiten, die mit zwei verschiedenen hydrodynamischen numerischen Modellen durch die Berechnung der MEAN (Durchschnitt) und Standardabweichung (STD) berechnet wurden.

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 Entwickler-Website (Esri).

Shapefil-Operationen

The geospatial shapefile vector format is an Esri invention. No wonder, arcpy is good at handling this vector data format. In hydraulic engineering, however, we usually create (draw) shapefiles manually either directly with ArcGIS or its open-source competitor QGIS to delineate, for example, particular flow regions. Examples can be found in the BASEMENT tutorial (explore the creation of elevation point, boundary polygon, and breakline polyline shapefiles). In codes, the processing of shapefiles only becomes important in the analysis of the output of numerical models (e.g., to classify morphological unit features, exactly calculate patch areas, or automatically place reinforcements in construction plans). At this stage, raster data (output of numerical models) must first be converted into shapefiles. This is why this tutorial starts with the conversion of raster data to shapefiles along with the illustration of other functions such as calculating patch area and accessing shapefile attribute tables.

Raster können in Polygon- und andere Shapefile-Typen (z. B. Punkt) umgewandelt werden. Das folgende Beispiel zeigt die Umwandlung eines Rasters in ein Polygon-Shapefil. Es verwendet ein ganzzahliges Raster aller Pixel, bei denen die Strömungstiefe und -geschwindigkeit kleiner als 1,4 m bzw. 0,15 m/s sind. Solche flachen und langsam fließenden Regionen werden Slackwater genannt (laut Wyrick & Pasternack (2014)). Da Slackwater als bevorzugter Lebensraum einiger aquatischer Arten bezeichnet wird, fragen wir uns jetzt, wie viel Slackwater-Gebiet das numerische Modell im simulierten Flussabschnitt vorhersagt. Dazu konvertieren wir das Slackwater-Raster in eine Shapefile und berechnen die Oberfläche der Shapefile mit den folgenden arcpy Methoden:

  • Convert the raster to a shapefile with arcpy.RasterToPolygon_conversion() with arguments:

    • in_raster ist ein arcpy.Raster()-Objekt von integer-Werten (die Verwendung eines Float-Rasters führt zu einem Fehler!).

    • out_polygon_features is a string of the output file name and directory.

    • simplify is an optional string that can be either "NO_SIMPLIFY" to force exact drawing of polygon boundaries along pixel border, or "SIMPLIFY" to enable polygon boundaries crossing pixels.

  • Add a new field to the new polygon shapefile with arcpy.AddField_management() with arguments:

    • field_name can be any string without blanks.

    • field_type ist ein String, der definiert, ob das Feld numerisch ist (z. B. "FLOAT" oder "LONG" für integer), Datum/Uhrzeit ("DATE"), "TEXT" oder "RASTER".

    • field_precision ist eine (optionale) * Integer* (Long), die die Anzahl der Ziffern definiert, die im neuen Feld gespeichert werden können.

    • More optional arguments can be set to define the number of decimals or characters, an alternative field name, enable NULL, a field domain, or if a field is required.

  • Calculate patch area with arcpy.CalculateGeometryAttributes_management() with arguments:

    • in_features ist ein String, der das Verzeichnis und den Namen eines Feature-Layers definiert.

    • geometry_property ist eine verschachtelte Liste von [[Target-field-name, Property], [Another-target-field-name, Another-property], ...] zur Berechnung geometrischer Eigenschaften wie "AREA", "HOLE_COUNT" oder "PART_COUNT" (und viele mehr).

    • area_unit kann "SQUARE_METERS" oder "SQUARE_KILOMETERS" sein (und viele andere Optionen für US-amerikanische Einheiten).

    • Weitere optionale Argumente können festgelegt werden, um Längeneinheiten (für Perimeterbewertungen) oder ein Koordinatensystem und -format zu definieren.

The below code block additionally features the application of these methods and also illustrates how the area value can be read row-by-row from the attribute table of a shapefile with arcpy.da.UpdateCursor(shapefile-name, field-name).

# 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))
References
  1. 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
Geodatenanalysen
Vectorize und Rasterize
Geodatenanalysen
Lernerfolgstest