地形鞍部的提取（ArcPy實現）

1、背景

相鄰兩山頭之間呈馬鞍形的低凹部分稱為鞍部。鞍部點是重要的地形控制點，它和山頂點、山谷點及山脊線、山谷線等構成地形特征點線，對地形具有很強的控制作用。因此，因此，對這些地形特征點、線的分析研究在數字地形分析中具有很重要的意義。同時，由於鞍部點的特殊地貌形態，是的鞍部點的提取方法較山頂點低谷底點更難，目前都還存在一定的技術局限性。

2、目的

利用水文分析的方法提取地形鞍部點；

通過多種GIS空間分析方法的應用，提高對知識的綜合運用能力。

3、要求

利用水文分析模塊和空間分析模塊相應功能提取樣區地形鞍部點。

4、數據

25m分辨率的DEM數據，面積約為59平方千米。（數據來自湯國安《arcgis地理信息系統空間分析實驗教程（第二版））

5、算法思想

鞍部具有獨特的形態特征，可被認為是原始地形中的山脊和反地形中的山脊回合的地方，因此可通過提取正反地形的山脊線並求其交點，獲得鞍部點，鞍部點的提取流程如下圖所示。

6、模型構建器

注意：這里山脊和反山脊的提取只需要提取到流量等於0的地方。

7、ArcPy實現

這里的代碼和利用水文分析方法提取山脊、山谷線的部分重復。

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 11-2 地形鞍部提取.py
# Created on: 2021-10-11 13:11:26.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import os
import shutil
import time

print time.asctime()
path = raw_input("請輸入數據所在文件夾的絕對路徑：").decode("utf-8")
# 開始計時
time_start = time.time()
paths = path + "\\result"
if not os.path.exists(paths):
    os.mkdir(paths)
else:
    shutil.rmtree(paths)
    os.mkdir(paths)

# Local variables:
dem = path + "\\dem"
Fill_dem = "Fill_dem"
Output_descent_rate_raster_zdx = "Descent_zdx"
minuend = "30000"
fu_fdx = "fu_fdx"
Abs_fdx = "Abs_fdx"
Output_descent_rate_raster_fdx = "Descent_fdx"
FlowDir_fdx = "FlowDir_fdx"
FlowAcc_fdx = "FlowAcc_fdx"
FlowAcc0_fdx = "FlowAcc0_fdx"
FlowDir_zdx = "FlowDir_zdx"
FlowAcc_zdx = "FlowAcc_zdx"
FlowAcc0_zdx = "FlowAcc0_zdx"
Times_FlowAc00 = "Times_FAc00"
meanDem = "meanDem"
Dem_Mean = "Dem_Mean"
zhengdixing = "zhengdixing"
Times_saddle = "Times_saddle"
raster_saddle = "raster_saddle"
vector_saddle = "鞍部點"
Contour_dem = "Contour_dem"
HillSha_dem = "HillSha_dem"

# Set Geoprocessing environments
print "Set Geoprocessing environments"
arcpy.env.scratchWorkspace = paths  # 臨時工作空間
arcpy.env.workspace = paths  # 工作空間
arcpy.env.extent = dem  # 處理范圍
arcpy.env.cellSize = dem  # 像元大小
arcpy.env.mask = dem  # 掩膜

# Process: 焦點統計
print "Process: 焦點統計"
arcpy.gp.FocalStatistics_sa(dem, meanDem, "Rectangle 11 11 CELL", "MEAN", "DATA", "90")

# Process: 減
print "Process: 減"
arcpy.gp.Minus_sa(dem, meanDem, Dem_Mean)

# Process: 條件測試
print "Process: 條件測試"
arcpy.gp.Test_sa(Dem_Mean, "value>0", zhengdixing)

# Process: 填窪
print "Process: 填窪"
arcpy.gp.Fill_sa(dem, Fill_dem, "")

# Process: 流向
print "Process: 流向"
arcpy.gp.FlowDirection_sa(Fill_dem, FlowDir_zdx, "NORMAL", Output_descent_rate_raster_zdx, "D8")

# Process: 流量
print "Process: 流量"
arcpy.gp.FlowAccumulation_sa(FlowDir_zdx, FlowAcc_zdx, "", "FLOAT", "D8")

# Process: 條件測試 (2)
print "Process: 條件測試 (2)"
arcpy.gp.Test_sa(FlowAcc_zdx, "value=0", FlowAcc0_zdx)

# Process: 減 (2)
print "Process: 減 (2)"
arcpy.gp.Minus_sa(dem, minuend, fu_fdx)

# Process: Abs
print "Process: Abs"
arcpy.gp.Abs_sa(fu_fdx, Abs_fdx)

# Process: 流向 (2)
print "Process: 流向 (2)"
arcpy.gp.FlowDirection_sa(Abs_fdx, FlowDir_fdx, "NORMAL", Output_descent_rate_raster_fdx, "D8")

# Process: 流量 (2)
print "Process: 流量 (2)"
arcpy.gp.FlowAccumulation_sa(FlowDir_fdx, FlowAcc_fdx, "", "FLOAT", "D8")

# Process: 條件測試 (3)
print "Process: 條件測試 (3)"
arcpy.gp.Test_sa(FlowAcc_fdx, "value=0", FlowAcc0_fdx)


# Process: 乘
print "Process: 乘"
arcpy.gp.Times_sa(FlowAcc0_fdx, FlowAcc0_zdx, Times_FlowAc00)

# Process: 乘 (2)
print "Process: 乘 (2)"
arcpy.gp.Times_sa(Times_FlowAc00, zhengdixing, Times_saddle)

# Process: 重分類
print "Process: 重分類"
arcpy.gp.Reclassify_sa(Times_saddle, "VALUE", "0 NODATA;1 1", raster_saddle, "DATA")

# Process: 柵格轉點
print "Process: 柵格轉點"
tempEnvironment0 = arcpy.env.outputZFlag
arcpy.env.outputZFlag = "Disabled"
tempEnvironment1 = arcpy.env.outputMFlag
arcpy.env.outputMFlag = "Disabled"
arcpy.RasterToPoint_conversion(raster_saddle, vector_saddle, "VALUE")
arcpy.env.outputZFlag = tempEnvironment0
arcpy.env.outputMFlag = tempEnvironment1

# Process: 等值線
print "Process: 等值線"
arcpy.gp.Contour_sa(dem, Contour_dem, "40", "0", "1", "CONTOUR", "")

# Process: 山體陰影
print "Process: 山體陰影"
arcpy.gp.HillShade_sa(dem, HillSha_dem, "315", "45", "NO_SHADOWS", "1")

save = [u"鞍部點", "contour_dem", "hillsha_dem"]
rasters = arcpy.ListRasters()
for raster in rasters:
    if raster.lower() not in save:
        print u"正在刪除{}圖層".format(raster)
        arcpy.Delete_management(raster)
# 結束計時
time_end = time.time()
# 計算所用時間
time_all = time_end - time_start
print time.asctime()
print "執行完畢！>>><<< 共耗時{:.0f}分{:.2f}秒".format(time_all // 60, time_all % 60)

8、結果

注：生成的鞍部點有點多，得后續配合等高線數據和暈渲圖對矢量形式的鞍部點數據進行編輯，剔除那些處於樣區邊緣以及內部的偽鞍部點。最后得到的鞍部點數據如下圖所示：

實驗結束 byebye·~