一、背景

隨著社會經濟發展需求，公路的重要性日益提高。在一些交通欠發達的地區，公路建設迫在眉睫。如何根據實際地形情況設計出比較合理的公路規劃，是一個值得研究的問題。

二、實驗目的：

（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~~

Collection：因為這些東西是非常簡單的。不要抱怨自己學不會，那是因為你沒有足夠用心。