Find the best path (ArcPy Implementation)

1, Background

With the demand of social and economic development, the importance of highway is increasing. In some areas with underdeveloped transportation, highway construction is imminent. How to design a more reasonable highway planning according to the actual terrain is a problem worthy of study.

2, Purpose of the experiment:

(1) Through practice, be familiar with ArcGIS grid data distance mapping, surface analysis, cost weight distance, data reclassification, shortest path and other spatial analysis functions;
(2) Master the basic process and operation process of using the above spatial analysis functions of ArcGIS to solve practical application problems.

3, Experimental data

(1) DEM (elevation data)

(2) startPot (path source point data)

(3) endPot (path end data)

(4) river (small watershed data)

4, Experimental software

Python

5, Experiment content:

The following points should be paid attention to in the site selection of the new highway:

1. The cost of new path is less;

2. The new path is a shorter path;

3. The new route should avoid the main river to reduce the cost;

4. When calculating the cost data of the new route, considering that the river cost (Reclass_river) is a key factor in the route cost, first merge the slope data (reclass_slope) and undulation data (reclass_QFD) according to the weight of 0.6:0.4, and then add and merge the same weight with the river cost. The formula is described as follows:
cost = reclassified watershed data + (reclassified slope data 0.6 + reclassified relief data 0.4)

5. The realization of finding the shortest path needs to be completed by using the functions of cost path and shortest path in distance mapping in ArcGIS Spatial Analyst, slope calculation and fluctuation calculation in surface analysis, reclassification and grid calculator

6. Finally, submit the shortest path and make a thematic map.

The main operation steps include:

(1) Create cost dataset

Slope cost dataset: the slope dataset is generated by Surface Analyst and divided into 10 levels at equal intervals. The minimum level of slope is assigned as 1 and the maximum level is assigned as 10.

Undulation cost data set: use 11 * 11 cell size as Neighborhood Statistics to obtain terrain undulation data, and reclassify it according to 10 levels of equal spacing. The more undulating the terrain is, the higher the level is assigned, that is, the minimum level is assigned as 1 and the maximum level is assigned as 10.

River cost data set: use the Reclassify command to classify according to the river level as follows: Level 4 is 10; The sequence is 8, 5, 2 and 1.

(2) Weighted consolidation of single factor cost data to generate the final cost data set.

(3) Use the Cost Distance function in Distance to calculate the cost weight Distance function.

(4) Select Cost Path in Distance to find the shortest path.

(5) Make thematic maps and output them.

Operation flow chart:

6, Model builder

7, ArcPy implementation

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-2 finding the best path.py
# Created on: 2021-10-10 10:02:59.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import os

path = raw_input("Please enter the absolute path where the data is located:").decode("utf-8")
paths = os.path.split(path)[0] + '\\result'
if not os.path.exists(paths):
    os.mkdir(paths)
# Local variables:
dem = path + "\\dem"

endPot = path + "\\endPot"
startPot = path + "\\startPot"
river = path + "\\river"
Reclass_rive = "Reclass_rive"
Slope_dem = "Slope_dem"
Reclass_Slop = "Reclass_Slop"
FocalSt_dem = "FocalSt_dem"
Reclass_Foca = "Reclass_Foca"
# rastercacu = "rastercacu"   # Look at the grid calculator, because it is defined there, so it is commented out here
CostDis_star = "CostDis_star"
huisuo = "huisuo"
result = "Cost path"

# Set Geoprocessing environments
print "Set Geoprocessing environments"
arcpy.env.scratchWorkspace = paths
arcpy.env.workspace = paths

# Process: slope
print "Process: slope"
arcpy.gp.Slope_sa(dem, Slope_dem, "DEGREE", "1", "PLANAR", "METER")

# Process: reclassification
print "Process: Reclassification"
arcpy.gp.Reclassify_sa(Slope_dem, "Value",
                       "0 5.730143 1;5.730143 11.460286 2;11.460286 17.190429 3;17.190429 22.920572 4;22.920572 28.650715 5;28.650715 34.380858 6;34.380858 40.111001 7;40.111001 45.841144 8;45.841144 51.571287 9;51.571287 57.301430 10",
                       Reclass_Slop, "DATA")

# Process: focus statistics
print "Process: Focus statistics"
arcpy.gp.FocalStatistics_sa(dem, FocalSt_dem, "Rectangle 3 3 CELL", "MEAN", "DATA", "90")

# Process: reclassification (2)
print "Process: Reclassification (2)"
arcpy.gp.Reclassify_sa(FocalSt_dem, "Value",
                       "179.462219 210.765332 1;210.765332 242.068445 2;242.068445 273.371558 3;273.371558 304.674670 4;304.674670 335.977783 5;335.977783 367.280896 6;367.280896 398.584009 7;398.584009 429.887122 8;429.887122 461.190234 9;461.190234 492.493347 10",
                       Reclass_Foca, "DATA")

# Process: reclassification (3)
print "Process: Reclassification (3)"
arcpy.gp.Reclassify_sa(river, "Value", "0 1;1 2;2 5;3 8;4 10", Reclass_rive, "DATA")

# Process: grid calculator
print "Process: Grid calculator"
# "\"%Reclass_rive%\" + (\"%Reclass_Slop%\" * 0.6 + \"%Reclass_Foca%\" * 0.4) "
rastercacu = arcpy.sa.Raster(Reclass_rive) + arcpy.sa.Raster(Reclass_Slop) * 0.6 + arcpy.sa.Raster(Reclass_Foca) * 0.4
rastercacu.save("rastercacu")

# Process: cost distance
print "Process: Cost distance"
arcpy.gp.CostDistance_sa(startPot, rastercacu, CostDis_star, "", huisuo, "", "", "", "", "")

# Process: cost path
print "Process: Cost path"
arcpy.gp.CostPath_sa(endPot, CostDis_star, huisuo, result, "EACH_CELL", "Id", "INPUT_RANGE")

save = [u"Cost path"]  # Layers to keep
features = arcpy.ListRasters()
for feature in features:
    if feature not in save:
        print u"Deleting layer{}".format(feature)
        arcpy.Delete_management(feature)
print "Operation completed~~~"

matters needing attention:
1. The grid calculator does not seem to be able to perform multiplication calculations directly.
arcpy.gp.RasterCalculator_sa("'%Reclass_rive%' + ('%Reclass_Slop%' * 0.6 + '%Reclass_Foca%' * 0.4) ", rastercacu)

When I change a floating-point number to an integer, I still report an error:

Let's just add. The results are as follows:

Check the following help documents:

Example expressions in the grid calculator tool dialog box
 In the grid calculator and directly in Python When using map algebra in, you should pay attention to some grammatical differences.

a,When using the grid calculator, the map algebra expression does not include the output name and the equal sign because there are output parameters specified in the grid calculator tool dialog box (=).  
b,Layer names can be used directly with operators only in the grid calculator tool dialog box. stay Python When processing in, layers must first be converted to grid objects. 
c,Similarly, the grid calculator variable can only be included in the percent sign in the tool dialog box (%) Or quotation marks (") Within.

The error is still reported. It may be that step b has not been done. If you are interested, you can continue to try again
Finally, I directly chose to use addition and multiplication instead of grid calculator.

2. Cost distance does not support output to PGDB workspace.
ERROR 010537: output parameter 'output distance raster data' is invalid: cost distance does not support output to PGDB workspace.
ERROR 010537: the output parameter 'output backtracking link grid data' is invalid: cost distance does not support output to PGDB workspace.
Execution (CostDistance) failed.
So instead, under the same path to the database, create a "result" folder to store the newly generated grid data.