一、背景
作為地形特征線的山脊線、山谷線對地形、地貌具有一定的控制作用。它們與山頂點、谷底點以及鞍部點等一起構成了地形起伏變化的骨架結構。同時由於山脊線具有分水性，山谷線具有合水性特征，使得它們在地形分析中具有特殊的意義。
二、目的
了解基於DEM水文分析方法提取山脊線和山谷線的原理;掌握水流方向、匯流累積量提取原理及方法。
三、要求
利用ArcGIS水文分析模塊提取出樣區的山脊線和山谷線。
四、數據
25m分辨率的DEM數據，區域面積約140km²（\ChP11 \Ex1目錄中)。
五、算法思想
山脊線和山谷線的提取實質上也是分水線與匯水線的提取。因此，可以利用水文分析的方法進行提取。
對於山脊線而言，由於它同時也是分水線，而分水線的性質即為水流的起源點。所以，通過地表徑流模擬計算之后，這些柵格的水流方向都應該只具有流出方向而不存在流入方向，即柵格的匯流累積量為零。因此，通過對零值的提取，就可得到分水線，即山脊線。
對於山谷線而言，可以利用反地形計算。即利用一個較大的數值減去原始DEM數據，得到與原始DEM地形相反的地形數據，使得原始的DEM中的山脊變成反地形的山谷，而原始DEM中的山谷在反地形中就變成了山脊。再利用山脊線的提取方法就可以實現山谷線的提取。但是此方法提取出的山脊和山谷位置有些偏差，可以利用正、負地形加以糾正。
流程圖
六、模型構建器
七、ArcPy實現
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 11-1 利用水文分析方法提取山脊線和山谷線.py
# Created on: 2021-10-11 10:09:40.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"
FlowDir_zdx = "FlowDir_zdx"
FlowAcc_zdx = "FlowAcc_zdx"
FlowAcc0_zdx = "FlowAcc0_zdx"
Filter_zdx = "Filter_zdx"
threshold_zdx0_5541 = "thresh_zdx0"
meanDem = "meanDem"
Dem_Mean = "Dem_Mean"
zhengdixing = "zhengdixing"
sjx = "sjx"
Ridge_line = "山脊線"
minuend = "5000"
fu_fdx = "fu_fdx"
Abs_fdx = "Abs_fdx"
Output_descent_rate_raster_fdx = "Descent_fdx"
fudixing = "fudixing"
FlowDir_fdx = "FlowDir_fdx"
FlowAcc_fdx = "FlowAcc_fdx"
FlowAcc0_fdx = "FlowAcc0_fdx"
Filter_fdx = "Filter_fdx"
threshold_fdx0_65677 = "thresh_fdx0"
sgx = "sgx"
Valley_line = "山谷線"
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.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: 條件測試 (3)
print "Process: 條件測試 (3)"
arcpy.gp.Test_sa(FlowAcc_zdx, "value=0", FlowAcc0_zdx)
# Process: 濾波器
print "Process: 濾波器"
arcpy.gp.Filter_sa(FlowAcc0_zdx, Filter_zdx, "LOW", "DATA")
# Process: 條件測試 (4)
print "Process: 條件測試 (4)"
arcpy.gp.Test_sa(Filter_zdx, "value>0.5541", threshold_zdx0_5541)
# 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.Times_sa(threshold_zdx0_5541, zhengdixing, sjx)
# Process: 重分類
print "Process: 重分類"
arcpy.gp.Reclassify_sa(sjx, "VALUE", "0 NODATA;1 1", Ridge_line, "DATA")
# 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.Test_sa(Dem_Mean, "value<0", fudixing)
# Process: 流量 (2)
print "Process: 流量 (2)"
arcpy.gp.FlowAccumulation_sa(FlowDir_fdx, FlowAcc_fdx, "", "FLOAT", "D8")
# Process: 條件測試 (5)
print "Process: 條件測試 (5)"
arcpy.gp.Test_sa(FlowAcc_fdx, "value=0", FlowAcc0_fdx)
# Process: 濾波器 (2)
print "Process: 濾波器 (2)"
arcpy.gp.Filter_sa(FlowAcc0_fdx, Filter_fdx, "LOW", "DATA")
# Process: 條件測試 (6)
print "Process: 條件測試 (6)"
arcpy.gp.Test_sa(Filter_fdx, "value>0.65677", threshold_fdx0_65677)
# Process: 乘 (2)
print "Process: 乘 (2)"
arcpy.gp.Times_sa(fudixing, threshold_fdx0_65677, sgx)
# Process: 重分類 (2)
print "Process: 重分類 (2)"
arcpy.gp.Reclassify_sa(sgx, "VALUE", "0 NODATA;1 1", Valley_line, "DATA")
# Process: 等值線
print "Process: 等值線"
arcpy.gp.Contour_sa(dem, Contour_dem, "50", "0", "1", "CONTOUR", "")
# Process: 山體陰影
print "Process: 山體陰影"
arcpy.gp.HillShade_sa(dem, HillSha_dem, "315", "45", "NO_SHADOWS", "1")
save = ["hillsha_dem", "contour_dem", u"山脊線", u"山谷線"]
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)
八、結果
九、其它
在上實驗課的時候，在老師那覓得相對上面，另一種更快捷的方式，而且也是比較通用，因為上面的那種方法，需要設置分界閾值，且的根據等值線和山體陰影來人工判斷，而下面這種方法，雖然也是用水文分析方法，但不需要設置分界閾值，且對所有dem較為通用。
下面讓我們來看看吧
模型構建器
跟上面的那種方法相比，主要省略了圈起來部分的步驟：
ArcPy實現
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 11-1 利用水文分析方法提取山脊線和山谷線2.py
# Created on: 2021-10-11 10:47:48.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"
meanDem = "meanDem"
Dem_Mean = "Dem_Mean"
zhengdixing = "zhengdixing"
FlowDir_zdx = "FlowDir_zdx"
FlowAcc_zdx = "FlowAcc_zdx"
FlowAcc0_zdx = "FlowAcc0_zdx"
sjx = "sjx"
Ridge_line = "山脊線"
minuend = "5000"
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"
fudixing = "fudixing"
sgx = "sgx"
Valley_line = "山谷線"
# 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.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.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.FlowAccumulation_sa(FlowDir_zdx, FlowAcc_zdx, "", "FLOAT", "D8")
# Process: 條件測試 (3)
print "Process: 條件測試 (3)"
arcpy.gp.Test_sa(FlowAcc_zdx, "value=0", FlowAcc0_zdx)
# Process: 乘
print "Process: 乘"
arcpy.gp.Times_sa(zhengdixing, FlowAcc0_zdx, sjx)
# Process: 重分類
print "Process: 重分類"
arcpy.gp.Reclassify_sa(sjx, "VALUE", "0 NODATA;1 1", Ridge_line, "DATA")
# 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: 條件測試 (5)
print "Process: 條件測試 (5)"
arcpy.gp.Test_sa(FlowAcc_fdx, "value=0", FlowAcc0_fdx)
# Process: 條件測試 (2)
print "Process: 條件測試 (2)"
arcpy.gp.Test_sa(Dem_Mean, "value<0", fudixing)
# Process: 乘 (2)
print "Process: 乘 (2)"
arcpy.gp.Times_sa(FlowAcc0_fdx, fudixing, sgx)
# Process: 重分類 (2)
print "Process: 重分類 (2)"
arcpy.gp.Reclassify_sa(sgx, "VALUE", "0 NODATA;1 1", Valley_line, "DATA")
save = [u"山脊線", u"山谷線"]
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)
結果
感覺上述二者差別不是很大，但是第一種的分界閾值對不同dem數據，需要再做判斷，而第二種不需要設置分界閾值，就很方便，所以，比較推薦第二種。
實驗結束 byebye~~~