ArcPy-Introduction au géotraitement avec Python avec Script Tool sur ArcGIS Pro-

Qu'est-ce que ArcPy

Il existe plusieurs façons d'utiliser Python avec ArcGIS. (Fenêtre Python, outil de script Python, boîte à outils Python, complément Python, ArcGIS for Python via Jupyter, etc.) Parmi ceux-ci, non seulement plusieurs outils de géotraitement existants (outils de géotraitement) peuvent être combinés comme le générateur de modèle, mais aussi L'outil Script permet l'automatisation du traitement, y compris la logique conditionnelle, et ArcPy est un module (ou package) qui fournit des classes et des fonctions permettant à Python d'accéder aux outils de géotraitement. est. Vous pouvez créer votre propre outil de géotraitement à l'aide de l'outil de script.

Comment créer un outil de script

[Créer un outil de script](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. Setting script tool parameters Alphabetical list of ArcPy functions (ESRI)

Créez un script Python et spécifiez le fichier de script dans ** Volet Catalogue> Boîte à outils> Cliquez avec le bouton droit sur Boîte à outils cible> Nouveau> Script **. Lorsque le nouvel écran de paramétrage du script s'affiche, sélectionnez ** [Général]> [Fichier de script] ** pour spécifier le fichier de script que vous souhaitez importer.

Dans le script, utilisez la fonction GetParameterAsText () pour stocker les données que vous souhaitez obtenir dans une variable comme indiqué ci-dessous.

script1.py

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

Dans les réglages des paramètres, définissez les paramètres en fonction du nombre (car c'est Python, il commence par 0). Vous pouvez saisir le nom du paramètre dans l'étiquette et définir le type de données, par exemple si vous souhaitez le saisir numériquement ou acquérir une couche existante. Si vous définissez les paramètres comme indiqué ci-dessus, l'outil de géotraitement comme indiqué ci-dessous sera créé.

J'ai essayé d'automatiser avec l'outil de script de prédiction d'inondation

Je voudrais montrer comment incorporer une série d'analyses d'informations géographiques dans plusieurs outils de script avec un exemple qui se réfère en fait à la prévision d'inondation suivante. J'utilise l'extension Spatial Analyst. Learn ArcGIS Lesson: Predicting floods using unit flow charts [Exigences]

• ArcGIS Pro
• Spatial Analyst Extension [Données utilisées et zone d'enquête] Le cours d'étude ESRI analyse dans un bassin affluent à Stowe, Vermont, États-Unis, mais le processus couvre un bassin hydrographique dans le sud de l'Australie.

Distinction entre liste de travail et automatisation / opération manuelle

J'aimerais pouvoir tout automatiser, mais je dois le faire fonctionner manuellement car je dois le saisir manuellement et il est plus rapide ou plus pratique de le faire manuellement que de l'automatiser avec ArcPy. Cette fois, la série de travaux a été divisée en trois outils de script (TASK1-3), et le travail manuel a été ajouté à l'intervalle et à la création du résultat final.

TASK1(Script1.py) --Paramètre1: Entrez la limite --Parameter2: saisissez les données de contour (theCont) --Paramètre3: Entrez les données d'élévation du point (le Spot)

1. Créez un DEM de la zone de levé à partir des données d'élévation du point et des données de contour (TopoToRaster ())
2. DEM complété par le lissage de la dépression (puits) qui interfère avec l'analyse du débit (Fill ())
3. Créez un raster de direction d'écoulement à partir du DEM créé (FlowDirection ())
4. Créez un raster de flux cumulatif (FlowAccumulation ())

TASK2(Script2.py) --Manuel: définir le point de sortie (point d'écoulement) --Parameter1: Entrez le DEM (DEM_input) créé par TASK1 --Paramètre2: Entrez les données de sortie créées manuellement (theOutlet) --Paramètre3: Spécifiez la distance lors de l'accrochage du point de ruissellement à la rivière

1. Accrochez le point de ruissellement à la rivière (SnapPourPoint ())
2. Créez un raster de bassin versant qui contribue au point de ruissellement (Watershed ())
3. Après avoir calculé la pente (Pente ()), calculez la vitesse de sortie à partir du raster de débit cumulé avec un ordinateur raster.
4. Excluez les rasters physiquement impossibles (trop tard ou trop tôt) à l'aide de l'évaluation conditionnelle à l'aide de Con (Con ())
5. Évaluez le temps nécessaire à l'eau pour s'écouler dans le canal et créez une carte isochrone

TASK3(Script3.py) --Manuel: Créez une table (isochrones.txt) qui demande de sortir la carte de ligne simultanée (carte isochrone) à des intervalles de 30 minutes (1800 secondes) --Paramètre1: entrez isochrones.txt --Paramètre2: Spécifiez le chemin de la destination de sortie de la table

1. Convertissez la carte de ligne simultanée en table à l'intervalle de temps spécifié (isochrones.txt) (TableToTable_conversion ())
2. Ajoutez un champ de zone (Area_sqm) à la table de sortie (AddField_management ()) et calculez la zone avec l'ordinateur raster.
3. Ajoutez un champ Unit Hydrograph Ordinate (UH_ordi) à la table (AddField_management ()) et calculez le débit toutes les 1800 secondes avec l'ordinateur raster.
4. Sortez la table vers le chemin spécifié dans le paramètre 2

Création du résultat final

Créez un graphique de flux unitaire au point de sortie à l'aide de la fonction Créer un graphique basée sur la table de sortie.

Conseils pour créer des outils de script

Importez les modules requis au début

script.py

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

En affichant un message au début et à la fin de chaque processus en utilisant la fonction ʻAddMessage () `, vous pouvez voir quel processus vous êtes en train de faire.

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

Le script réellement utilisé

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