[NEW] 利用水文分析方法提取山脊线和山谷线（ArcPy实现） | flowacc – NATAVIGUIDES

作为地形特征线的山脊线、山谷线对地形、地貌具有一定的控制作用。它们与山顶点、谷底点以及鞍部点等一起构成了地形起伏变化的骨架结构。同时由于山脊线具有分水性，山谷线具有合水性特征，使得它们在地形分析中具有特殊的意义。

一、背景

作为地形特征线的山脊线、山谷线对地形、地貌具有一定的控制作用。它们与山顶点、谷底点以及鞍部点等一起构成了地形起伏变化的骨架结构。同时由于山脊线具有分水性，山谷线具有合水性特征，使得它们在地形分析中具有特殊的意义。

二、目的

了解基于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~~~

相邻两山头之间呈马鞍形的低凹部分称为鞍部。鞍部点是重要的地形控制点，它和山顶点、山谷点及山脊线、山谷线等构成地形特征点线，对地形具有很强的控制作用。因此，因此，对这些地形特征点、线的分析研究在数字地形分析中具有很重要的意义。同时，由于鞍部点的特殊地貌形态，是的鞍部点的提取方法较山顶点低谷底点更难，目前都还存在一定的技术局限性。

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