相邻两山头之间呈马鞍形的低凹部分称为鞍部。鞍部点是重要的地形控制点,它和山顶点、山谷点及山脊线、山谷线等构成地形特征点线,对地形具有很强的控制作用。因此,因此,对这些地形特征点、线的分析研究在数字地形分析中具有很重要的意义。同时,由于鞍部点的特殊地貌形态,是的鞍部点的提取方法较山顶点低谷底点更难,目前都还存在一定的技术局限性。
利用水文分析的方法提取地形鞍部点;
通过多种GIS空间分析方法的应用,提高对知识的综合运用能力。
利用水文分析模块和空间分析模块相应功能提取样区地形鞍部点。
25m分辨率的DEM数据,面积约为59平方千米。(数据来自汤国安《arcgis地理信息系统空间分析实验教程(第二版))
鞍部具有独特的形态特征,可被认为是原始地形中的山脊和反地形中的山脊回合的地方,因此可通过提取正反地形的山脊线并求其交点,获得鞍部点,鞍部点的提取流程如下图所示。
注意:这里山脊和反山脊的提取只需要提取到流量等于0的地方。
下面的只是构建了主要部分,另外的部分与利用水文分析方法提取山脊、山谷线的重复
这里的代码和利用水文分析方法提取山脊、山谷线的部分重复。
# -*- 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)