一、背景
隨着社會經濟發展需求，公路的重要性日益提高。在一些交通欠發達的地區，公路建設迫在眉睫。如何根據實際地形情況設計出比較合理的公路規划，是一個值得研究的問題。
二、實驗目的：
（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 os
path = 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”文件夾存放新生成柵格數據。
八、運行結果
實驗結束 byebye~~