地形鞍部的提取（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·~