寻找最佳路径（ArcPy实现）

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

一、背景

随着社会经济发展需求，公路的重要性日益提高。在一些交通欠发达的地区，公路建设迫在眉睫。如何根据实际地形情况设计出比较合理的公路规划，是一个值得研究的问题。

二、实验目的：

（1）通过练习，熟悉 ArcGIS 栅格数据距离制图、表面分析、成本权重距离、数据重分类、最短路径等空间分析功能；

（2）熟练掌握利用 ArcGIS 上述空间分析功能解决实际应用问题的基本流程和操作过程。

三、实验数据

（1）dem(高程数据)

（2）startPot (路径源点数据)

（3）endPot (路径终点数据)

（4）river (小流域数据)

四、实验软件

Python

五、实验内容：

新公路选址需注意如下几点：

1、新建路径成本较少；

2、新建路径为较短路径；

3、新建路径的选择应该避开主干河流，以减少成本；

4、新建路径的成本数据计算时，考虑到河流成本（Reclass_river）是路径成本中较关键因素，先将坡度数据（reclass_slope）和起伏度数据（reclass_QFD）按照 0.6：0.4 权重合并，然后与河流成本作等权重的加和合并，公式描述如下：

cost =重分类流域数据+ (重分类坡度数据0.6+重分类起伏度数据0.4)

5、寻找最短路径的实现需要运用 ArcGIS 的空间分析（Spatial Analyst）中距离制图中的成本路径及最短路径、表面分析中的坡度计算及起伏度计算、重分类及栅格计算器等功能完成

6、最后提交寻找到的最短路径路线，并制作专题图。

主要操作步骤包括：

（1）创建成本数据集

坡度成本数据集：使用Surface Analyst生成坡度数据集，并采用等间距分为 10 级，坡度最小一级赋值为 1，最大一级赋值为 10。

起伏度成本数据集:使用11*11单元格大小做Neighborhood Statistics，得到地形起伏度数据，并按 10 级等间距实施重分类，地形越起伏，级数赋值越高，即最小一级赋值为 1，最大一级赋值为 10。

河流成本数据集：使用Reclassify 命令，按照河流等级如下进行分类：4 级为 10；如此依次为8，5，2，1。

（2）加权合并单因素成本数据，生成最终成本数据集。

（3）使用 Distance 中的 Cost Distance功能计算成本权重距离函数。

（4）选择 Distance 中的 Cost Path求取最短路径。

（5）制作专题地图并输出。

操作流程图：

六、模型构建器

七、ArcPy实现

# -*- coding: utf-8 -*- # --------------------------------------------------------------------------- # 8-2寻找最佳路径.py # Created on: 2021-10-10 10:02:59.00000 # (generated by ArcGIS/ModelBuilder) # Description: # ---------------------------------------------------------------------------# Import arcpy module import arcpy import ospath = raw_input("请输入数据所在绝对路径：").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" # 看栅格计算器那里，因为那里定义了，所以这里注释掉 CostDis_star = "CostDis_star" huisuo = "huisuo" result = "成本路径"# Set Geoprocessing environments print "Set Geoprocessing environments" arcpy.env.scratchWorkspace = paths arcpy.env.workspace = paths# Process: 坡度 print "Process: 坡度" arcpy.gp.Slope_sa(dem, Slope_dem, "DEGREE", "1", "PLANAR", "METER")# Process: 重分类 print "Process: 重分类" 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: 焦点统计 print "Process: 焦点统计" arcpy.gp.FocalStatistics_sa(dem, FocalSt_dem, "Rectangle 3 3 CELL", "MEAN", "DATA", "90")# Process: 重分类 (2) print "Process: 重分类 (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: 重分类 (3) print "Process: 重分类 (3)" arcpy.gp.Reclassify_sa(river, "Value", "0 1;1 2;2 5;3 8;4 10", Reclass_rive, "DATA")# Process: 栅格计算器 print "Process: 栅格计算器" # "\"%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: 成本距离 print "Process: 成本距离" arcpy.gp.CostDistance_sa(startPot, rastercacu, CostDis_star, "", huisuo, "", "", "", "", "")# Process: 成本路径 print "Process: 成本路径" arcpy.gp.CostPath_sa(endPot, CostDis_star, huisuo, result, "EACH_CELL", "Id", "INPUT_RANGE")save = [u"成本路径"] # 需要保留的图层 features = arcpy.ListRasters() for feature in features: if feature not in save: print u"正在删除图层{}".format(feature) arcpy.Delete_management(feature) print "运行完毕~~~"

注意事项：

1、栅格计算器好像不可以直接执行乘计算。

arcpy.gp.RasterCalculator_sa("'%Reclass_rive%' + ('%Reclass_Slop%' * 0.6 + '%Reclass_Foca%' * 0.4) ", rastercacu)

当我把浮点数改为整数，仍报错：

那就只运算加法呢，结果如下：

查阅了一下帮助文档：

栅格计算器工具对话框中的示例表达式在“栅格计算器”和直接在 Python 中使用“地图代数”时，您应注意语法上存在一些差异。a、当使用栅格计算器时，由于在栅格计算器工具对话框中有指定的输出参数，“地图代数”表达式不包括输出名称和等号 (=)。 b、只有在栅格计算器工具对话框中，图层名称才可以直接与运算符一起使用。在 Python 中进行处理时，图层必须首先转换为栅格对象。 c、同样，“栅格计算器”变量仅在工具对话框中才能包含在百分号 (%) 或引号 (") 之内。

还是报错，可能是b步骤没有做，不弄了，有兴趣的可以再继续尝试哈

最后我直接选择使用加法，乘法运算来代替栅格计算器。

2、成本距离不支持输出到 PGDB 工作空间。

ERROR 010537: 输出参数 '输出距离栅格数据' 无效：成本距离不支持输出到 PGDB 工作空间。

ERROR 010537: 输出参数 '输出回溯链接栅格数据' 无效：成本距离不支持输出到 PGDB 工作空间。

执行(CostDistance)失败。

所以改成在给该数据库同路径下，新建一个“result”文件夹存放新生成栅格数据。