利用水文分析方法提取山脊、山谷线（ArcPy实现）

本文主要是介绍利用水文分析方法提取山脊、山谷线（ArcPy实现），对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

山脊线和山谷线的提取实质上也是分水线和汇水线的提取。

对于山脊线来说，由于他同时也是分水线，而分水线的性质即为水流的起源点。所以，通过地表径流模拟计算后这些栅格的水流方向都应该只具有流出方向而不存在流入方向，即栅格的汇流积累量为0.因此，通过对零值得提取，就可以得到分水线，即山脊线。

对于山谷线，可以利用反地形计算，即利用一个较大的数值减去原始数值DEM数据，得到与原始数据相反的地形数据，使得原始的DEM山脊变成反地形的山谷，而原始数据DEM中的山谷在反地形中就变成了山脊。再利用山脊线的提取方法就可以实现山谷线的提取。但是此方法提取的山脊和山谷位置有些 偏差，可以利用正负地形加以纠正。

1、正负地形的提取

对原始的DEM数据进行焦点统计，计算像元平均值，dem数据变光滑。

用原始德玛数据减去焦点统计后的数据，再进行重分类，大于0的地形为正地形，小于0的地形为负地形。

2、山脊线的提取

填洼：操作对原始dem进行填洼，因为对所有的洼地填充所以Z limit为默认值。

流向：计算水流方向，

栅格计算器;提取汇流量为零值

做出原始dem 的等值线和山体阴影。

模型构建器

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")
# 用于提取dem所在的文件路径，例如 E:\data
File_Geodatabase_location = os.path.split(input_dem)[0]
Output_file_Geodatabase = File_Geodatabase_location
File_Geodatabase_name = raw_input("请输入文件地理数据库名称：").decode("utf-8")
# Contour_spacing = raw_input("请输入等值线间距：").decode("utf-8")
# 等值线间距，默认50，可自行修改
Contour_spacing = "50"
#开始计时
time_start = time.time()
print "以下的文件地理数据名称请不要重复，否则会删除再重新新建，或者自己修改一下代码也可"
# 新建的文件地理数据库的绝对路径
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"  # 取值value=0
Filter_Flo1 = "Filter_Flo1"  # 滤波器
Test_Filter_Tes1_tif = "Test_Filter_Flo"
Focus_statistics = "焦点统计"
D_M = "D_M"  # dem-meandem
zhengdixing = "正地形"
Times_Filter_Te1 = "Times_Filter_Te1"  # *
sgx = "山谷线"
sjx = "山脊线"
minuted = "3000"  # 用于求反地形，值需要比dem的最大值大即可
Negative_anti_terrain = "fu反地形"
Anti_terrain = "abs_反地形"  # 取绝对值
Output_descent_rate_raster_data2 = "Output_descent_rate_raster_data2"
FlowDir_fdx = "FlowDir_fdx"
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"  # *

# 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")

# 删除多余要素
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)

这篇关于利用水文分析方法提取山脊、山谷线（ArcPy实现）的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！