在之前有很多小伙伴问我如何更准确的提取山顶点，在查阅了许多资料之后，我终于找出来了一篇相关算法（引用文章见文末），以下为此次分享的全部内容。
在GIS系统中，对于山顶点、洼地点、鞍部点的提取，都可以通过ArcGIS来完成，但是其实现原理我们却并不了解，因此，借此次机会，让我们一起来探索一下如何利用ArcPy模块完成对山顶点的提取吧。
首先，第一步，还是从需求分析开始，根据题意，我们知道需要从DEM数据中提取出山顶点，那么，需要我们对山顶点的概念有一定了解，而山顶点，也可以理解为是局部区域的最高点，因此，我们只需要提取局部区域最高点即可。
根据需求，我们可以将局部区域设置为窗口大小为3，5，7……,(2n+1)的矩形区域，接下来，只需要判断这个区域的中心点是否是这个区域中的最大值，最大值在这个区域内是否唯一就可以了，如果这个区域的中心点的值大于这个区域内的其他值，并且当前区域的最大值具有唯一性，那么我们便可以认为这个区域的中心点就是我们所要提取的山顶点。
以下，为此思路的核心代码
edge = int((Windows_length - 1) / 2)# Windows_length为窗口大小
dataX = list()#记录山顶点位置
dataY = list()#记录山顶点位置
newArr = np.zeros([rows,cols],dtype = int)#存储山顶点位置
for i in range(edge,rows-edge):#根据窗口大小，依次遍历数组
for j in range(edge,cols-edge):
ExtractData = data[i - edge:i + edge + 1, j - edge:j + edge + 1]#取得局部区域
if IsMiddleMax(ExtractData.flatten()) == True:
dataX.append(j)#存储山顶点位置
dataY.append(i)#存储山顶点位置
newArr[i,j] = 1#存储山顶点位置到数组
#判断局部区域的中心点是否是最大值，判断这个最大值在局部区域中是否唯一，两者都满足返回True，否则返回False
def IsMiddleMax(ArrMat):
arrLength = ArrMat.shape[0]
index = int(arrLength/2)
if ArrMat[index] == np.max(ArrMat) and np.count_nonzero(ArrMat == np.max(ArrMat)) == 1:
return True
else:
return False以下部分为将山顶点位置写入矢量数据的代码，其中MaxPointArr是记录山顶点位置的数组，lowerLeft为仿射矩阵左下角坐标位置，Widthcellsize,Heightcellsize分别为像元的宽度和高度，NoData为0
MaxPointRaster=arcpy.NumPyArrayToRaster(MaxPointArr,lowerLeft,Widthcellsize,Heightcellsize,value_to_nodata=NoData)#将数组转换为Raster
arcpy.ra.ConvertRasterToFeature(MaxPointRaster,'Value','POINT')#将Raster转换为Feature至此，所有核心步骤完成，接下来一部分，是利用matplotlib进行制图，展示一下我们的完成效果
plt.figure('ElevationMaxPoint')#图框标题
plt.subplots_adjust(left=0, right=1, top=0.99, bottom=0.01)#设置窗口位置及大小
plt.scatter(dataX,dataY,c='r',marker='.',s=120,alpha=1) #将局部最高点的坐标通过点叠加在图像上，dataX,dataY为存放山顶点位置的数据
plt.xticks([])#去除X轴刻度
plt.yticks([])#去除Y轴刻度
plt.show() #展示结果至此，所有计算步骤完成，以下为此次分享的代码的全部内容
import arcpy
import numpy as np
import matplotlib.pyplot as plt
import math
def ExtractMax(data,Windows_length,NoData):
rows,cols = data.shape#获取行数和列数
newArr = np.zeros([rows,cols],dtype = int)
edge = int((Windows_length - 1) / 2)
if rows - 1 - edge <= edge or cols - 1 - edge <= edge:
print("The parameter k is to large.")
return None,None,None
dataX = list()
dataY = list()
for i in range(edge,rows-edge):
for j in range(edge,cols-edge):
ExtractData = data[i - edge:i + edge + 1, j - edge:j + edge + 1]
if IsMiddleMax(ExtractData.flatten()) == True:
dataX.append(j)
dataY.append(i)
newArr[i,j] = 1
return dataX,dataY,newArr
#当窗口中间数值大于周围其他值时并且最大值为数组的唯一值时，返回True
def IsMiddleMax(ArrMat):
arrLength = ArrMat.shape[0]
index = int(arrLength/2)
if ArrMat[index] == np.max(ArrMat) and np.count_nonzero(ArrMat == np.max(ArrMat)) == 1:
return True else:
return False
######################计算坡向
#计算dx,dy
def calcFiniteSlopes(elevGrid, dx):
Zbc = assignBCs(elevGrid)
Sx = (Zbc[1:-1, :-2] - Zbc[1:-1, 2:]) / (2 * dx) # WE方向
Sy = (Zbc[2:, 1:-1] - Zbc[:-2, 1:-1]) / (2 * dx) # NS方向
return Sx, Sy
#栅格最外圈加一圈
def assignBCs(elevGrid):
ny, nx = elevGrid.shape
Zbc = np.zeros((ny + 2, nx + 2))
Zbc[1:-1, 1:-1] = elevGrid
Zbc[0, 1:-1] = elevGrid[0, :]
Zbc[-1, 1:-1] = elevGrid[-1, :]
Zbc[1:-1, 0] = elevGrid[:, 0]
Zbc[1:-1, -1] = elevGrid[:, -1]
Zbc[0, 0] = elevGrid[0, 0]
Zbc[0, -1] = elevGrid[0, -1]
Zbc[-1, 0] = elevGrid[-1, 0]
Zbc[-1, -1] = elevGrid[-1, 0]
return Zbc
#计算坡向
def AspectG(Sx,Sy,aspect):
for i in range(Sx.shape[0]):
for j in range(Sy.shape[1]):
sx = float(Sx[i, j])
sy = float(Sy[i, j])
if (sx == 0.0) & (sy == 0.0):
aspect[i, j] = -1
elif sx == 0.0:
if sy > 0.0:
aspect[i, j] = 0.0
else:
aspect[i, j] = 180.0
elif sy == 0.0:
if sx > 0.0:
aspect[i, j] = 90.0
else:
aspect[i, j] = 270.0
else:
aspect[i, j] = float(math.atan2(sy, sx) * 57.29578)
if aspect[i, j] < 0.0:
aspect[i, j] = 90.0 - aspect[i, j]
elif aspect[i, j] > 90.0:
aspect[i, j] = 360.0 - aspect[i, j] + 90.0
else:
aspect[i, j] = 90.0 - aspect[i, j]
return aspect
##################################
def Run(FileName,Windows_length,savepath=""):
tif = arcpy.Raster(FileName)#读取栅格
NoData = tif.noDataValue#读取Nodata
data = arcpy.RasterToNumPyArray(tif, nodata_to_value=NoData).astype(np.float)#读取栅格数据的数组
dataX,dataY,MaxPointArr = ExtractMax(data,Windows_length,NoData)
ExtentXmin = tif.extent.XMin#取x坐标最小值
ExtentYmin = tif.extent.YMin#取y坐标最小值
lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin)#取得数据起始点坐标(左下角坐标)
Widthcellsize=tif.meanCellWidth#像元宽度
Heightcellsize=tif.meanCellHeight#像元高度
###计算坡向坡度
Sx, Sy = calcFiniteSlopes(data, Widthcellsize)
slope = np.arctan(np.sqrt(Sx ** 2 + Sy ** 2)) * 57.29578#坡度
aspect = np.ones([Sx.shape[0], Sx.shape[1]]).astype(np.float32)
aspectA = AspectG(Sx,Sy,aspect)#坡向
######
###输出局部高程点位置
if savepath!="":
MaxPointRaster = arcpy.NumPyArrayToRaster(MaxPointArr,lowerLeft,Widthcellsize,
Heightcellsize,value_to_nodata=NoData)
arcpy.ra.ConvertRasterToFeature(MaxPointRaster,'Value','POINT')
return dataX,dataY,slope,aspectA
if __name__=='__main__':
dataX,dataY,slope,aspectA=Run(r"D:\作业\qq接的项目\Python提取局部高程点\DEM显示局部最高点\data\Dem2.tif",25)
plt.figure('ElevationMaxPoint')
plt.imshow(aspectA, cmap = 'gray')
plt.subplots_adjust(left=0, right=1, top=0.99, bottom=0.01)#设置窗口位置及大小
plt.scatter(dataX,dataY,c='r',marker='.',s=120,alpha=1) #将局部最高点的坐标通过点叠加在图像上
plt.xticks([])#去除X轴刻度
plt.yticks([])#去除Y轴刻度
plt.show() #展示结果
以下为此次代码效果展示
注：此算法仅能实现局部区域至高点提取，山顶点的提取需结合等值线共同完成，在理想情况下，此算法的精度会随着区域窗口的变大而提高，但在达到某临界点时，精度会骤然降低。
参考文献：[1]仲腾,汤国安,周毅,李若殷,张维.基于反地形DEM的山顶点自动提取[J].测绘通报,2009(04):35-37.
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
