山顶点提取（ArcPy实现）

一、背景

山顶点指哪些在特定邻域分析范围内，该点都比周围点高的区域。山顶点是地形的重要特征点，它的分布与密度反映了地貌的发育特征，同时也制约着地貌发育。因此，如何基于DEM数据正确有效的提取山顶点，在数字地形分析中具有重要意义。

二、目的

通过等高线、山顶点、洼地点的提取和配置、引导读者熟练掌握利用ArcGIS栅格数据空间分析中等高线的提取、栅格数据邻域分析和窗口计算功能，完成栅格数据表面分析。

三、数据:

黄土丘陵地区 1:10000 DEM 数据（Chp8\Ex5\）。

四、要求

(1)应用栅格数据空间分析模块中的等高线提取功能，分别提取等高距为15米和75米的等高线图，并按标准地形图绘制等高线方法绘制等高线，作为山顶点提取的地形背景;

(2)通过邻域分析和栅格计算器提取山顶点。

五、流程
加载DEM数据->加载Spatial Analyst模块，邻域分析->地图代数->重分类->栅格转点

用DEM生成间隔为15m和75m的等高线，生成山体阴影结果图，二者构成地形晕渲图以辅助判断山顶点位置。对DEM数据进行焦点统计分析，以21*21的窗口进行处理，将生成的结果与DEM数据做差重分类后可得到栅格形式的山顶点数据。将栅格数据转为矢量后结合地形晕渲图删除不合理的山顶点，即得到山顶点的分布。

流程图如下所示：

六、模型构建器

七、ArcPy实现

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-5 山顶点提取.py
# Created on: 2021-10-10 22:29:58.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
from arcpy import Raster, ListRasters, Delete_management
from arcpy.sa import Test
import os
import shutil

path = raw_input("请输入数据所在文件夹的绝对路径：").decode("utf-8")
paths = path + "\\result"
if not os.path.exists(paths):
    os.mkdir(paths)
else:
    shutil.rmtree(paths)
    os.mkdir(paths)

# Local variables:
dem = path + "\\dem"
Contour_dem15 = "Contour_dem15"
Contour_dem75 = "Contour_dem75"
HillSha_dem = "hillSha_dem"
FocalSt_dem = "FocalSt_dem"
F_dem = "F_dem"
raster_peak = "raster_peak"
vector_peak = "山顶点"

# 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: 等值线 15m"
arcpy.gp.Contour_sa(dem, Contour_dem15, "15", "0", "1", "CONTOUR", "")

# Process: 等值线 (2)
print "Process: 等值线 75m"
arcpy.gp.Contour_sa(dem, Contour_dem75, "75", "0", "1", "CONTOUR", "")

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

# Process: 焦点统计
print "Process: 焦点统计"
arcpy.gp.FocalStatistics_sa(dem, FocalSt_dem, "Rectangle 21 21 CELL", "MAXIMUM", "DATA", "90")

# Process: 栅格计算器
# arcpy.gp.RasterCalculator_sa("\"%FocalSt_dem%\" - \"%dem%\" == 0", F_dem)

# Process: FocalSt_dem - dem == 0
print "Process: FocalSt_dem - dem == 0"
Test(Raster(FocalSt_dem) - Raster(dem), "VALUE=0").save(F_dem)

# Process: 重分类
print "Process: 重分类"
arcpy.gp.Reclassify_sa(F_dem, "Value", "0 NODATA;1 1", raster_peak, "DATA")

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

save = ["contour_dem15", "contour_dem75", "hillsha_dem", u"山顶点"]
rasters = ListRasters()
for raster in rasters:
    if raster.lower() not in save:
        print u"正在删除{}图层".format(raster)
        Delete_management(raster)
print "运行完毕~~~"

八、结果

注：提取的山顶点的点多不多，是受块统计中分析窗口大小的影响，窗口越大提取的点越少，但是窗口过大的话，将会漏掉一些重要的山顶点。对提取的结果，可人工判断删除局部点。

实验结束 byebye~~