ArcPy-Einführung in die Geoverarbeitung mit Python mit dem Skript-Tool in ArcGIS Pro-

Was ist ArcPy?

Es gibt mehrere Möglichkeiten, Python mit ArcGIS zu verwenden. (Python-Fenster, Python-Skript-Tool, Python-Toolbox, Python-Add-In, ArcGIS für Python über Jupyter usw.) Von diesen können nicht nur mehrere vorhandene Geoverarbeitungswerkzeuge (Geoverarbeitungswerkzeuge) wie der Modellbauer kombiniert werden, sondern auch Das Skript-Tool ermöglicht die Automatisierung der Verarbeitung einschließlich der bedingten Logik. Das Modul (oder Paket), das Klassen und Funktionen für Python bereitstellt, um auf das Geoverarbeitungs-Tool zuzugreifen, ist ArcPy. ist. Mit dem Skript-Tool können Sie Ihr eigenes Geoverarbeitungswerkzeug erstellen.

So erstellen Sie ein Skript-Tool

[Skript-Tool erstellen](https://pro.arcgis.com/de/pro-app/help/analysis/geoprocessing/basics/create-a-python-script-tool.htm#:~:text=Script% 20tools% 20are% 20geoprocessing% 20tools% 20that% 20execute% 20a% 20script% 20or% 20executable% 20file. & Text = ArcPy% 20 bietet% 20access% 20to% 20geoprocessing, komplexe% 20workflows% 20quickly% 20und% 20easily.) Setting script tool parameters Alphabetical list of ArcPy functions (ESRI)

Erstellen Sie ein Python-Skript und geben Sie die Skriptdatei unter ** [Katalogfenster]> [Toolbox]> [Klicken Sie mit der rechten Maustaste auf die Ziel-Toolbox]> [Neu]> [Skript] ** an. Wenn der neue Bildschirm für die Skripteinstellungen angezeigt wird, wählen Sie ** [Allgemein]> [Skriptdatei] **, um die zu importierende Skriptdatei anzugeben.

Verwenden Sie im Skript die Funktion "GetParameterAsText ()", um die Daten, die Sie abrufen möchten, in einer Variablen zu speichern, wie unten gezeigt.

script1.py

# Get parameters
theBoundary = arcpy.GetParameterAsText(0)
theCont = arcpy.GetParameterAsText(1)
theSpot = arcpy.GetParameterAsText(2)

Stellen Sie in den Parametereinstellungen die Parameter entsprechend der Nummer ein (beginnend mit 0, da es sich um Python handelt). Sie können den Parameternamen in die Beschriftung eingeben und den Datentyp festlegen, z. B. ob Sie ihn numerisch eingeben oder eine vorhandene Ebene abrufen möchten. Wenn Sie die Einstellungen wie oben gezeigt vornehmen, wird das unten gezeigte Geoverarbeitungswerkzeug erstellt.

Ich habe versucht, mit dem Skript-Tool für die Hochwasservorhersage zu automatisieren

Ich möchte anhand eines Beispiels, das sich tatsächlich auf die folgende Hochwasservorhersage bezieht, zeigen, wie eine Reihe von Analysen geografischer Informationen in mehrere Skriptwerkzeuge integriert werden kann. Ich verwende die Spatial Analyst Extension. ArcGIS-Lektion lernen: Hochwasser mithilfe von Einheitsflussdiagrammen vorhersagen [Bedarf]

• ArcGIS Pro
• Spatial Analyst Extension [Verwendete Daten und Erhebungsgebiet] Der ESRI-Studiengang analysiert in einem Nebeneinzugsgebiet in Stowe, Vermont, USA, aber der Prozess zielt auf ein Flusseinzugsgebiet in Südaustralien ab.

Unterscheidung zwischen Arbeitsliste und Automatisierung / manueller Bedienung

Ich möchte in der Lage sein, alles zu automatisieren, aber ich muss es manuell bedienen, weil ich es manuell eingeben muss und es schneller oder bequemer ist, es manuell zu tun, als es mit ArcPy zu automatisieren. Dieses Mal wurde die Arbeitsserie in drei Skriptwerkzeuge (TASK1-3) unterteilt, und manuelle Arbeit wurde zur Intervall- und Endergebniserstellung hinzugefügt.

TASK1(Script1.py) --Parameter1: Geben Sie die Grenze ein --Parameter2: Geben Sie die Konturdaten ein (theCont) --Parameter3: Geben Sie die Punkthöhendaten (den Punkt) ein.

1. Erstellen Sie ein DEM des Vermessungsbereichs aus den Punkthöhen- und Konturdaten (TopoToRaster ()).
2. DEM wird durch Glätten der Vertiefung (Senke) abgeschlossen, die die Durchflussanalyse stört (Fill ()).
3. Erstellen Sie ein Flussrichtungsraster aus dem erstellten DEM (FlowDirection ()).
4. Erstellen Sie ein kumulatives Flow-Raster (FlowAccumulation ())

TASK2(Script2.py) --Handbuch: Stellen Sie den Abflusspunkt (Pour Point) ein. --Parameter1: Geben Sie das von TASK1 erstellte DEM (DEM_input) ein --Parameter2: Geben Sie manuell erstellte Outlet-Daten ein (theOutlet) --Parameter3: Geben Sie die maximale Entfernung (die Entfernung) an, wenn Sie den Abflusspunkt zum Fluss einrasten lassen

1. Schnappen Sie den Abflusspunkt zum Fluss (SnapPourPoint ())
2. Erstellen Sie ein Einzugsgebietsraster, das zum Abflusspunkt beiträgt (Watershed ()).
3. Berechnen Sie nach der Berechnung der Steigung (Slope ()) die Abflussgeschwindigkeit aus dem Raster für die kumulative Flussrate mit einem Rastercomputer.
4. Schließen Sie physikalisch unmögliche (zu späte oder zu frühe) Raster mit der bedingten Bewertung mit Con (Con ()) aus.
5. Bewerten Sie die Zeit, die Wasser benötigt, um durch den Kanal zu fließen, und erstellen Sie eine Isochronen-Karte

TASK3(Script3.py) --Handbuch: Erstellen Sie eine Tabelle (isochrones.txt), die anweist, die simultane Linienkarte (Isochronenkarte) in Intervallen von 30 Minuten (1800 Sekunden) an die Tabelle auszugeben. --Parameter1: Geben Sie isochrones.txt ein --Parameter2: Geben Sie den Pfad des Tabellenausgabeziels an

1. Konvertieren Sie die simultane Linienzuordnung im angegebenen Zeitintervall (isochrones.txt) in eine Tabelle (TableToTable_conversion ())
2. Fügen Sie der Ausgabetabelle (AddField_management ()) ein Flächenfeld (Area_sqm) hinzu und berechnen Sie die Fläche mit dem Rastercomputer.
3. Fügen Sie der Tabelle (AddField_management ()) ein Unit Hydrograph Ordinate-Feld (UH_ordi) hinzu und berechnen Sie die Durchflussrate pro 1800 Sekunden mit dem Rastercomputer.
4. Geben Sie die Tabelle unter dem in Parameter2 angegebenen Pfad aus

Erstellung des Endergebnisses

Erstellen Sie am Abflusspunkt ein Einheitsflussdiagramm mit der Funktion Diagramm erstellen basierend auf der Ausgabetabelle.

Tipps zum Erstellen von Skriptwerkzeugen

Importieren Sie zu Beginn die erforderlichen Module

script.py

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

Verwenden Sie die Funktion "AddMessage ()", um zu Beginn und am Ende jedes Prozesses eine Nachricht anzuzeigen, damit Sie sehen können, welchen Prozess Sie ausführen.

script3.py

### Calculate the area of Raster according to each attribute
arcpy.AddMessage("If the resolution of the raster is not 30m, modify the python script")
fieldname = "Area_sqm"
expression = '(!Count! * 30 * 30 )'
arcpy.CalculateField_management(outReclassTable, fieldname, expression)
arcpy.AddMessage("The area of each classes calculated")

Das tatsächlich verwendete Skript

TASK1(Script1.py)

script1.py

# -*- coding: utf-8 -*-
"""
@author: hiro-ishi
TASK1:Script 1
"""

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("Spatial")

# Get parameters
theBoundary = arcpy.GetParameterAsText(0)
theCont = arcpy.GetParameterAsText(1)
theSpot = arcpy.GetParameterAsText(2)

# Set env variables
env.overwriteOutput = True

# Set local variables
theBnd = "Boundary"
arcpy.CopyFeatures_management(theBoundary, theBnd)

#Create output variables for names
theClipCont = "Contours_clip"
theClipSpot = "Spotheight_clip"
theDEM = "DEM_fill"
theFDir = "DEM_fill_FlowDir"
theFAccum = "DEM_fill_FlowAccum"

# =============================================================================
# Create and fill a DEM by point data
# =============================================================================
###Clip out contours
arcpy.Clip_analysis(theCont,theBnd,theClipCont)
arcpy.Clip_analysis(theSpot,theBnd,theClipSpot)
    
#Set topo variables
inContours = TopoContour([[theClipCont, "ELEVATION"]])
inSpotHeight = TopoPointElevation ([[theClipSpot, "ELEVATION"]])
inBoundary = TopoBoundary([theBnd])
inFeatures = [inContours, inSpotHeight, inBoundary]
    
###Create DEM
arcpy.AddMessage("Creating DEM")
outTTR = TopoToRaster(inFeatures, cell_size = 30)
 
###Fill DEM and save
arcpy.AddMessage("Running fill")
outFill = Fill(outTTR)
outFill.save(theDEM)

###Add a filled DEM to the current map
aprx = arcpy.mp.ArcGISProject("CURRENT")
aprxMap = aprx.listMaps("Map")[0] 
aprxMap.addDataFromPath(outFill)
arcpy.AddMessage("The filled DEM named " + theDEM + " was added to the map")
arcpy.AddMessage("Then, let's delineate the watershed")

TASK2(Script2.py)

script1.py

# -*- coding: utf-8 -*-
"""
@author: hiro-ishi
TASK2: Script 2
"""

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("Spatial")

# Get parameters
DEM_input = arcpy.GetParameterAsText(0)
theOutlet = arcpy.GetParameterAsText(1)
theDistance = arcpy.GetParameterAsText(2)

# Set env variables
env.overwriteOutput = True

#Create output variables for names
theDEM = "DEM_fill"
theFDir = "DEM_fill_FlowDir"
theFAccum = "DEM_fill_FlowAccum"
thePPoint = "PPoint_raster"
theWShed = "DEM_fill_WShed"
theSlope = "DEM_fill_Slope"
sbAc = "DEM_fill_sbAc"
velUnlim = "DEM_VelField_unlim"
velField = "DEM_VelField"
theFDirClip = "DEM_fill_WShed_FlowDir"
theLength = "DEM_time"
isochrones = "DEM_time_isochrones"
isochrones_table = "DEM_time_isochrones_table"

# =============================================================================
# Delineate the watershed
# =============================================================================
### Snap the pour point to raster
arcpy.AddMessage("Let's delineate the watershed")
arcpy.AddMessage("Snapping Poor Point")
outSnapPour = SnapPourPoint(theOutlet, theFAccum, theDistance)
outSnapPour.save(thePPoint)
arcpy.AddMessage("Snapped Poor Point Raster generated")

### Add the Pour point raster to the current map
aprx = arcpy.mp.ArcGISProject("CURRENT")
aprxMap = aprx.listMaps("Map")[0] 
aprxMap.addDataFromPath(outSnapPour)
arcpy.AddMessage("The snapped pour point (raster data) was added to the map")

### Deliineate the Watershed
arcpy.AddMessage("Delineating the watershed")
outWatershed = Watershed(theFDir, outSnapPour)
outWatershed.save(theWShed)

### Add the watershed to the current map
aprxMap.addDataFromPath(outWatershed)
arcpy.AddMessage("The watershed area was added to the map")
arcpy.AddMessage("You finished to delineate the watershed")

# =============================================================================
# Create a velocity field	
# =============================================================================
### Create a Slope field
arcpy.AddMessage("Calculating slope")
outSlope = Slope(DEM_input, "Percent rise")
outSlope.save(theSlope)

### Raster calculation to acquire sbAc
arcpy.AddMessage("Creating a velocity field")
out_term = SquareRoot(theSlope) * SquareRoot(theFAccum)

### Extract by mask
out_term_clip = ExtractByMask(out_term, theWShed)
out_term_clip.save(sbAc)

### To acquire statistics of the mean slop-area
arcpy.AddMessage("Acquire a Mean Value")
sbAc_mean = arcpy.GetRasterProperties_management(out_term_clip, "MEAN")
arcpy.AddMessage("Mean Value: " + sbAc_mean[0])

### Create a velocity field
# If there is no "float()", the value 15.97.... comes as a string
velField = 0.1 * out_term_clip / float(sbAc_mean[0])
velField.save(velUnlim)

velLower = Con(velUnlim, velUnlim, 0.02, "Value >= 0.02")
velocity = Con(velLower, velLower, 2, "Value <= 2")

arcpy.AddMessage("Velocity Field was acquired!!")

# =============================================================================
# Create an isochrone map
# =============================================================================
###Raster calculator
arcpy.AddMessage("Acquiring weight raster")
outWeight = 1 / velocity

theFDir_clip = ExtractByMask(theFDir, theWShed)
theFDir_clip.save(theFDirClip)

### Create a Flow length field
arcpy.AddMessage("Acquiring flow length")
outFlowLength = FlowLength(theFDirClip, "downstream", outWeight)
outFlowLength.save(theLength)

### Add the flow length field to the current map
aprxMap.addDataFromPath(outFlowLength)
arcpy.AddMessage("The flowtime map named " + theLength + " was added to the map")

arcpy.AddMessage("Finish")
DEM_time_max = arcpy.GetRasterProperties_management(outFlowLength, "MAXIMUM")
arcpy.AddMessage("The time it takes water to flow to the outlet ranges from 0 seconds (rain that falls on the outlet itself) to " + DEM_time_max[0] + "seconds!!")
arcpy.AddMessage("Finally, let's reclassify the table to create a graph of a Hydrograph!! Let's go to the script3")

TASK3(Script3.py)

script3.py

# -*- coding: utf-8 -*-
"""
@author: hiro-ishi
TASK3:Script3 
"""

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("Spatial")

# Set env variables
env.overwriteOutput = True

# Get parameters
in_raster = arcpy.GetParameterAsText(0)
tableoutpath = arcpy.GetParameterAsText(1)

#Create output variables for names
isochrones = "DEM_time_reclass"
isochrones_table = "DEM_time_reclass_table"

# =============================================================================
# Create a unit hydrograph
# =============================================================================
### Convert a raster attribute into a table
outReclassTable = arcpy.TableToTable_conversion(in_raster, tableoutpath, isochrones_table)
arcpy.AddMessage( isochrones_table + ".dbf was saved to the folder")

### Add a Field
arcpy.AddField_management(outReclassTable, "Area_sqm", "DOUBLE", field_alias = "Area(Sq.Meters)")
arcpy.AddMessage("A field added")

### Calculate the area of Raster according to each attribute
arcpy.AddMessage("If the resolution of the raster is not 30m, modify the python script")
fieldname = "Area_sqm"
expression = '(!Count! * 30 * 30 )'
arcpy.CalculateField_management(outReclassTable, fieldname, expression)
arcpy.AddMessage("The area of each classes calculated")

### Add Field
arcpy.AddField_management(outReclassTable, "UH_ordi", "DOUBLE", field_alias = "UnitHydrographOrdinate")
arcpy.AddMessage("A field added")

### Calculate the area of Raster according to each attribute
fieldname = "UH_ordi"
expression = '(!Area_sqm! / 1800 )'
arcpy.CalculateField_management(outReclassTable, fieldname, expression)
arcpy.AddMessage("Unit Hydrograph Ordinate was calculated")

### Add the table to the current map
aprx = arcpy.mp.ArcGISProject("CURRENT")
aprxMap = aprx.listMaps("Map")[0] 
aprxMap.addDataFromPath(outReclassTable)

arcpy.AddMessage("Then, please create a graph manually.")