ArcPy-Introduction to Geoprocessing with Python with Scripting Tools on ArcGIS Pro-

What is ArcPy

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.

How to make a 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 tried to automate with the flood prediction script tool

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]

• ArcGIS Pro
• Spatial Analyst Extension [Data used and survey area] The ESRI study course analyzes a tributary catchment in Stowe, Vermont, USA, but the process covers a river catchment in southern South Australia, Australia.

Distinguishing between work list and automation / manual operation

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)

1. Create a DEM of the survey area from the point elevation data and contour data (TopoToRaster ())
2. DEM completed by smoothing the depression (sink) that interferes with flow analysis (Fill ())
3. Create a flow direction raster from the created DEM (FlowDirection ())
4. Create a cumulative flow raster (FlowAccumulation ())

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

1. Snap the runoff point to the river (SnapPourPoint ())
2. Create a catchment raster that contributes to the runoff point (Watershed ())
3. After calculating the slope (Slope ()), calculate the outflow velocity from the cumulative flow rate raster with a raster calculator.
4. Exclude physically impossible (too late or too early) rasters using conditional evaluation using Con (Con ())
5. Evaluate the time it takes for water to flow through the channel and create an isochrone map

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

1. Convert the simultaneous line map to a table at the specified time interval (isochrones.txt) (TableToTable_conversion ())
2. Add an area field (Area_sqm) to the output table (AddField_management ()) and calculate the area with the raster calculator.
3. Add a Unit Hydrograph Ordinate field (UH_ordi) to the table (AddField_management ()) and calculate the flow rate per 1800 seconds with the raster calculator.
4. Output the table to the path specified in Parameter2

Final result creation

Create a unit flow graph at the outflow point using the CreateChart function based on the output table.

Tips for creating script tools

Import the required modules at the beginning

script.py

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

By displaying a message at the beginning and end of each process using the ʻAddMessage ()` function, you can see what process you are doing now.

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")

The script actually used

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.")