There are multiple ways to use Python with ArcGIS. (Python window, Python script tool, Python toolbox, Python addin, ArcGIS for Python via Jupyter, etc.) Of these, not only can you combine multiple existing geoprocessing tools like the Model builder, but you can also combine multiple existing geoprocessing tools. The Script tool enables the automation of processing including conditional logic, and the module (or package) that provides classes and functions for Python to access the geoprocessing tool is ArcPy. is. You can create your own geoprocessing tool by using the script tool.
[Create Script Tool](https://pro.arcgis.com/en/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% 20provides% 20access% 20to% 20geoprocessing, complex% 20workflows% 20quickly% 20and% 20easily.) Setting script tool parameters Alphabetical list of ArcPy functions (ESRI)
Create a Python script and specify the script file from ** [Catalog Pane]> [Toolbox]> [Right-click the target toolbox]> [New]> [Script] **. When the new script setting screen is displayed, select ** [General]> [Script File] ** to specify the script file you want to import.
In the script, use the GetParameterAsText ()
function to store the data you want to get in a variable as shown below.
script1.py
# Get parameters
theBoundary = arcpy.GetParameterAsText(0)
theCont = arcpy.GetParameterAsText(1)
theSpot = arcpy.GetParameterAsText(2)
In the parameter settings, set the parameters according to the number (starting with 0 because it is Python). You can enter the parameter name in the label and set the data type, such as whether to enter it numerically or to acquire an existing layer. If you make the settings as shown above, the geoprocessing tool as shown below will be created.
I would like to show how to incorporate a series of geographic information analysis into multiple script tools with an example that actually refers to the following flood forecast. I am using the Spatial Analyst Extension. Learn ArcGIS Lesson: Flood Prediction Using Unit Flow Charts [Requirements]
It would be great if everything could be automated, but I have to manually operate it because I have to enter it manually and it is faster or more convenient to do it manually than to automate it with ArcPy. This time, I divided the series of work into three script tools (TASK1-3) and added manual work in between and creating the final result.
TASK1(Script1.py) --Parameter1: Enter the Boundary --Parameter2: Input contour data (theCont) --Parameter3: Enter the point elevation data (the Spot)
TASK2(Script2.py) --Manual: Set the outflow point (Pour Point) --Parameter1: Enter the DEM (DEM_input) created by TASK1 --Parameter2: Enter manually created outlet data (theOutlet) --Parameter3: Specify the maximum distance (the Distance) when snapping the runoff point to the river
TASK3(Script3.py) --Manual: Create a table (isochrones.txt) that instructs the simultaneous line map (isochrone map) to be output as a table at intervals of 30 minutes (1800 seconds). --Parameter1: Enter isochrones.txt --Parameter2: Specify the path of the table output destination
Create a unit flow graph at the outflow point using the CreateChart function based on the output table.
script.py
# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *
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")
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.")
Recommended Posts