地形鞍部的提取（ArcPy实现）

2023-10-23 20:53:58

1、背景

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

2、目的

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

​ 通过多种GIS空间分析方法的应用，提高对知识的综合运用能力。

3、要求

​ 利用水文分析模块和空间分析模块相应功能提取样区地形鞍部点。

4、数据

​ 25m分辨率的DEM数据，面积约为59平方千米。（数据来自汤国安《arcgis地理信息系统空间分析实验教程（第二版））

5、算法思想

​ 鞍部具有独特的形态特征，可被认为是原始地形中的山脊和反地形中的山脊回合的地方，因此可通过提取正反地形的山脊线并求其交点，获得鞍部点，鞍部点的提取流程如下图所示。
地形鞍部的提取（ArcPy实现）

6、模型构建器

地形鞍部的提取（ArcPy实现）
注意：这里山脊和反山脊的提取只需要提取到流量等于0的地方。
下面的只是构建了主要部分，另外的部分与利用水文分析方法提取山脊、山谷线的重复
地形鞍部的提取（ArcPy实现）

7、ArcPy实现

这里的代码和利用水文分析方法提取山脊、山谷线的部分重复。

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 山谷线和山脊线提取(new).py
# Created on: 2021-10-04 09:45:49.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

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

print time.asctime()
# Local variables:
print "绝对路径输入格式，例如：E:\\data\\dem"
input_dem = raw_input("请输入表面栅格数据路径(绝对路径)：").decode("utf-8")
File_Geodatabase_location = os.path.split(input_dem)[0]
Output_file_Geodatabase = File_Geodatabase_location
print "以下的文件地理数据名称请不要重复，否则会删除再重新新建，或者自己修改一下代码也可"
File_Geodatabase_name = raw_input("请输入文件地理数据库名称：").decode("utf-8")
# Contour_spacing = raw_input("请输入等值线间距：").decode("utf-8")
Contour_spacing = "50"  # 默认50

time_start = time.time()
path = File_Geodatabase_location + "\\" + File_Geodatabase_name + ".gdb"
if os.path.exists(path):
    print path + " 已存在，正在删除..."
    shutil.rmtree(path)

contour = "等高线"
Mountain_shadow = "山体阴影"
fill_dem = "dem填洼"
Output_descent_rate_raster_data1 = "Output_descent_rate_raster_data1"
FlowDir_dem = "FlowDir_dem"
FlowAcc_FlowDir1 = "FlowAcc_FlowDir1"
Test_FlowAcc_Fl1 = "Test_FlowAcc_Flow0"
Filter_Flo1 = "Filter_Flo1"
Test_Filter_Tes1_tif = "Test_Filter_Flo"
Focus_statistics = "焦点统计"
D_M = "D_M"
zhengdixing = "正地形"
Times_Filter_Te1 = "Times_Filter_Te1"
sgx = "山谷线"
sjx = "山脊线"
minuted = "3000"
Negative_anti_terrain = "fu反地形"
Anti_terrain = "abs_反地形"
Output_descent_rate_raster_data2 = "Output_descent_rate_raster_data2"
FlowDir_fdx = "FlowDir_反地形"
FlowAcc_fdx = "FlowAcc_fdx"
Test_FlowAcc_fdx = "Test_FlowAcc_fdx0"
Filter_fdx = "Filter_fdx"
Test_Filter_fdx = "Test_Filter_fdx"
fudixing = "负地形"
Times_Filt_fudixing = "Times_Filt_fudixing"

# 提取鞍部的参数
Test_FlowAcc_fdx0 = "Test_FlowAcc_fdx0"
Test_FlowAcc_Flow0 = "Test_FlowAcc_Flow0"
Times_00 = "Times_00"
zhengdixing = "正地形"
anbu = "anbudian"
Reclass_anbudian = "Reclass_anbudian"
anbudian = "鞍部点"

# Process: 创建文件地理数据库
print "Process: 创建文件地理数据库{}.gdb".format(File_Geodatabase_name)
arcpy.CreateFileGDB_management(File_Geodatabase_location, File_Geodatabase_name, "CURRENT")
arcpy.env.workspace = path

# Process: 等值线
print "Process: 等值线"
arcpy.gp.Contour_sa(input_dem, contour, Contour_spacing, "0", "1", "CONTOUR", "")

# Process: 山体阴影
print "Process: 山体阴影"
arcpy.gp.HillShade_sa(input_dem, Mountain_shadow, "315", "45", "NO_SHADOWS", "1")

# Process: 填洼
print "Process: 填洼"
arcpy.gp.Fill_sa(input_dem, fill_dem, "")

# Process: 流向
print "Process: 流向"
arcpy.gp.FlowDirection_sa(fill_dem, FlowDir_dem, "NORMAL", Output_descent_rate_raster_data1, "D8")

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

# Process: 条件测试 (3)
print "Process: 条件测试 (3)"
arcpy.gp.Test_sa(FlowAcc_FlowDir1, "value=0", Test_FlowAcc_Fl1)

# Process: 滤波器
print "Process: 滤波器"
arcpy.gp.Filter_sa(Test_FlowAcc_Fl1, Filter_Flo1, "LOW", "DATA")

# Process: 条件测试 (4)
print "Process: 条件测试 (4)"
arcpy.gp.Test_sa(Filter_Flo1, "value>0.5541", Test_Filter_Tes1_tif)

# Process: 焦点统计
print "Process: 焦点统计"
arcpy.gp.FocalStatistics_sa(input_dem, Focus_statistics, "Rectangle 11 11 CELL", "MEAN", "DATA", "90")

# Process: 减
print "Process: 减"
arcpy.gp.Minus_sa(input_dem, Focus_statistics, D_M)

# Process: 条件测试
print "Process: 条件测试"
arcpy.gp.Test_sa(D_M, "value>0", zhengdixing)

# Process: 乘
print "Process: 乘"
arcpy.gp.Times_sa(Test_Filter_Tes1_tif, zhengdixing, Times_Filter_Te1)

# Process: 重分类
print "Process: 重分类"
arcpy.gp.Reclassify_sa(Times_Filter_Te1, "VALUE", "0 NODATA;1 1", sjx, "DATA")

# Process: 减 (2)
print "Process: 减 (2)"
arcpy.gp.Minus_sa(input_dem, minuted, Negative_anti_terrain)

# Process: Abs
print "Process: Abs"
arcpy.gp.Abs_sa(Negative_anti_terrain, Anti_terrain)

# Process: 流向 (2)
print "Process: 流向 (2)"
arcpy.gp.FlowDirection_sa(Anti_terrain, FlowDir_fdx, "NORMAL", Output_descent_rate_raster_data2, "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", Test_FlowAcc_fdx)

# Process: 滤波器 (2)
print "Process: 滤波器 (2)"
arcpy.gp.Filter_sa(Test_FlowAcc_fdx, Filter_fdx, "LOW", "DATA")

# Process: 条件测试 (6)
print "Process: 条件测试 (6)"
arcpy.gp.Test_sa(Filter_fdx, "value>0.65677", Test_Filter_fdx)

# Process: 条件测试 (2)
print "Process: 条件测试 (2)"
arcpy.gp.Test_sa(D_M, "value<0", fudixing)

# Process: 乘 (2)
print "Process: 乘 (2)"
arcpy.gp.Times_sa(Test_Filter_fdx, fudixing, Times_Filt_fudixing)

# Process: 重分类 (2)
print "Process: 重分类 (2)"
arcpy.gp.Reclassify_sa(Times_Filt_fudixing, "VALUE", "0 NODATA;1 1", sgx, "DATA")

# 提取鞍部的函数
# Process: 乘
arcpy.gp.Times_sa(Test_FlowAcc_fdx0, Test_FlowAcc_Flow0, Times_00)

# Process: 乘 (2)
arcpy.gp.Times_sa(Times_00, zhengdixing, anbu)

# Process: 重分类
arcpy.gp.Reclassify_sa(anbu, "VALUE", "0 NODATA;1 1", Reclass_anbudian, "DATA")

# Process: 栅格转点
tempEnvironment0 = arcpy.env.outputZFlag
arcpy.env.outputZFlag = "Disabled"
tempEnvironment1 = arcpy.env.outputMFlag
arcpy.env.outputMFlag = "Disabled"
arcpy.RasterToPoint_conversion(Reclass_anbudian, anbudian, "VALUE")
arcpy.env.outputZFlag = tempEnvironment0
arcpy.env.outputMFlag = tempEnvironment1


# 删除多余要素
rasters = arcpy.ListRasters()

# 需要保留的要素图层
# save = [u"山谷线", u"山脊线", u"正地形", u"负地形", u"等值线", u"山体阴影"]
save = [u"山谷线", u"山脊线", u"正地形", u"等值线", u"山体阴影", u"Test_FlowAcc_fdx0", u"Test_FlowAcc_Flow0"]
for raster in rasters:
    if raster 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)
上一篇：Class path contains multiple SLF4J bindings.解决方案


下一篇：Diff-prime Pairs 牛客