寻找最佳路径（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 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp8\Ex2"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
a = Reclassify(r'Road.mdb\river', "Value", "0 1;1 2;2 5;3 8;4 10")  # Process: river重分类
b = Slope('Road.mdb\dem')  # Process: dem坡度
c = Reclassify(b, "Value", '0.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.30143 10') # Process: Slope_dem重分类
d = FocalStatistics('Road.mdb\dem', "Rectangle 3 3 CELL", "MEAN")  # Process: dem焦点统计
e = Reclassify(d, "Value", '179.462219 210.7653318 1;210.7653318 242.0684446 2;242.0684446 273.3715574 3;273.3715574 304.6746702 4;304.6746702 335.977783 5;335.977783 367.2808958 6;367.2808958 398.5840086 7;398.5840086 429.8871214 8;429.8871214 461.1902342 9;461.1902342 492.493347 10')  # Process: FocalSt_dem重分类
f = a + c * 0.6 + e * 0.4  # Process: 栅格计算
g = CostDistance('Road.mdb\startPot.shp', f, "", 'backlink')  # Process: 成本距离
CostPath('Road.mdb\endPot.shp', g, 'backlink').save('成本路径')  # Process: 成本路径

注意事项：
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”文件夹存放新生成栅格数据。