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. コメント 2020-08-15 1800000.png

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. コメント 2020-08-15 180300.png If you make the settings as shown above, the geoprocessing tool as shown below will be created. コメント 2020-08-15 175900.png

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]

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.

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

Recommended Posts

ArcPy-Introduction to Geoprocessing with Python with Scripting Tools on ArcGIS Pro-
Strategy on how to monetize with Python Java
Introduction to Python with Atom (on the way)
Connect to MySQL with Python on Raspberry Pi
I tried to implement Minesweeper on terminal with python
Yum command to access MySQL with Python 3 on Linux
I want to AWS Lambda with Python on Mac!
Something to enjoy with Prim Pro (X-Play) and Python
Connect to BigQuery with Python
How to install Python2.7 python3.5 with pyenv (on RHEL5 CentOS5) (2016 Nov)
Use Python 3 introduced with command line tools on macOS Catalina
A collection of competitive pro techniques to solve with Python
Connect to Wikipedia with Python
Post to slack with Python 3
[Ev3dev] How to display bmp image on LCD with python
Update python on Mac to 3.7-> 3.8
Switch python to 2.7 with alternatives
Write to csv with Python
Save images on the web to Drive with Python (Colab)
A note on what you did to use Flycheck with Python
I tried Python! ] Can I post to Kaggle on my iPad Pro?
Linking Python and Arduino to display IME On / Off with LED
How to use python put in pyenv on macOS with PyCall
I tried to summarize everyone's remarks on slack with wordcloud (Python)
Python: How to use async with
Link to get started with python
[Python] Write to csv file with Python
Create folders from '01' to '12' with python
Nice to meet you with python
Try to operate Facebook with Python
Introduction to Python Hands On Part 1
Output to csv file with Python
Convert list to DataFrame with python
MP3 to WAV conversion with Python
To do tail recursion with Python2
How to get started with Python
What to do with PYTHON release?
Unable to install Python with pyenv
How to use FTP with Python
Notes on using rstrip with python.
How to calculate date with python
Easily post to twitter with Python 3
Steps to install python3 on mac
Getting started with Python 3.8 on Windows
I want to debug with Python
[Memo] Tweet on twitter with python
Update Python on Mac from 2 to 3
[Python] I tried to visualize the night on the Galactic Railroad with WordCloud!
Put Cabocha 0.68 on Windows and try to analyze the dependency with Python
I want to tweet on Twitter with Python, but I'm addicted to it
Data integration from Python app on Linux to Amazon Redshift with ODBC
Add 95% confidence intervals on both sides to the diagram with Python / Matplotlib
Use python on Raspberry Pi 3 to light the LED with switch control!
Steps to create a Python virtual environment with VS Code on Windows
How to draw a vertical line on a heatmap drawn with Python seaborn
I tried with the top 100 PyPI packages> I tried to graph the packages installed on Python
Data integration from Python app on Windows to Amazon Redshift with ODBC
Try to poke DB on IBM i with python + JDBC using JayDeBeApi
Migrate Django applications running on Python 2.7 to Python 3.5
Try logging in to qiita with Python
Change Python 64bit environment to 32bit environment with Anaconda