山頂點提取（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~~