寻找起伏度最佳窗口、使用arcpy求取地形起伏度的最佳统计单元

1、引言

原文地址链接：https://blog.csdn.net/haichao062/article/details/38318525

地形起伏度相关的研究很多，而地形起伏度数据，需要从DEM中提取。查了下文献，有可靠依据的是“均值变点法”。均值变点法是一种对非线性数据进行处理的数量统计方法，该方法对恰有一个变点的检验最为有效。推荐详细读下2008年发表在《兰州大学学报》的《新​疆​地​势​起​伏​度​的​分​析​研​究》。顺便提一下，我首先找到的是发表在《地理科学》2012年01月，第32卷第1期的《基于变点分析法提取地势起伏度》这篇文章，也使用的这种方法。

2、技术流程

2.1 计算不同格网大小的地形起伏度

使用arcgis的邻域分析功能，依次计算2*2,3*3. ...... 45*45 格网大小的平均起伏度

import arcpy  
import os  
from arcpy import sa  
  
"""DEM路径"""  
DEM=ur'D:\BaihcWorkSpace\checkDEM\ASTGTM2_N36E114.tif'  
  
"""输出文件夹路径"""  
outputPath=ur'D:\BaihcWorkSpace\checkDEM'  
  
'''''邻域计算'''  
def BlockStatistics(DEM,nbr,statistics_type):  
    outBlockStat = sa.BlockStatistics(DEM,nbr,statistics_type)  
    return outBlockStat  
  
'''''计算制定格网的地形起伏度'''  
def calc_gridRA(gridLength,DEM,outputFolder=None):  
    nbr=sa.NbrRectangle(gridLength,gridLength,'CELL') #邻域分析的窗口大小  
    rasterMax=BlockStatistics(DEM,nbr,'MAXIMUM')  
    rasterMin=BlockStatistics(DEM,nbr,'MINIMUM')  
    RA=rasterMax-rasterMin  
      
    if outputFolder:  
        output=os.path.join(outputPath,ur'RA_%s'%(gridLength))  
        RA.save(output)  
    return RA.mean  
  
RAList=[]  
  
  
for i in range(2,45):  
    RAMean=calc_gridRA(i,DEM)  
    RAList.append(RAMean)  

2.2 使用均值变点法寻找拐点

"numrange"为2*2,3*3....格网对应的平均起伏度。 
gridLength 为 DEM数据的分辨率。 

import numpy  
import math  
  
class calcReliefAmplitude():  
    def __init__(self,numrange,gridLength ):  
        self.numrange=numrange  
        self.len=len(self.numrange)  
        self.gridLength=gridLength  
        self.bestWindow=self.calcBestWindow()  
    def _calcVar(self):  
        return numpy.var(self.numrange)*self.len  
      
    def calcBestWindow(self):  
        meanRAList=[math.log(a*10000/math.pow((self.griLength*(i+2)),2)) for i,a in enumerate(self.numrange)]  
        S=[]  
        for j in range(2,len(meanRAList)+1):  
            s1=numpy.var(meanRAList[:j-1])*(j-1)   
            s2=numpy.var(meanRAList[j-1:])*(len(meanRAList)-j+1)  
            S.append( s1+s2)  
        bestWindow= S.index( min(S) )+1+2  
        return bestWindow  

2.3 求得并输出最佳统计单元下的起伏度数据

DEMbestWindow= calcReliefAmplitude(rasterList,gridLength).bestWindow  
  
calc_gridRA(DEMbestWindow,DEM,outputPath) 

完整代码如下，安装过arcgis10.2的即可使用

import arcpy  
from arcpy import sa  
import os  
import numpy    
import math    
    
class calcReliefAmplitude():    
    def __init__(self,numrange,gridLength ):    
        self.numrange=numrange    
        self.len=len(self.numrange)    
        self.gridLength=gridLength    
    def _calcVar(self):    
        return numpy.var(self.numrange)*self.len    
        
    def calcBestWindow(self):    
        meanRAList=[math.log(a/math.pow((self.gridLength*(i+2)),2)) for i,a in enumerate(self.numrange)]    
        S=[]         
        for j in range(2,len(meanRAList)+1):    
            s1=numpy.var(meanRAList[:j-1])*(j-1)     
            s2=numpy.var(meanRAList[j-1:])*(len(meanRAList)-j+1)    
            S.append( s1+s2)  
              
        bestWindow= S.index( min(S) )+1+2    
          
        return [bestWindow,math.pow((bestWindow*self.gridLength),2)]    
  
  
def BlockStatistics(block,nbr,statistics_type):  
    outBlockStat = sa.BlockStatistics(block,nbr,statistics_type)  
    return outBlockStat  
  
  
def get_qfdRast(i):  
    nbr=sa.NbrRectangle(i,i,'CELL')          
    rasterMax=BlockStatistics(block,nbr,'MAXIMUM')  
    rasterMin=BlockStatistics(block,nbr,'MINIMUM')  
    qfd=rasterMax-rasterMin  
    return qfd  
  
  
  
  
def getRasterList():  
    arcpy.SetProgressor("step", "process...",0, 50, 1)  
    numrange=[]  
    rasterList=[]  
    for i in range(2,50):  
        arcpy.SetProgressorLabel("%s/%s"%((i,50)))   
        arcpy.SetProgressorPosition()  
        qfd=get_qfdRast(i)  
        if qfdparams=='mean':  
            qfdValue=qfd.mean  
        if qfdparams=='maximum':  
            qfdValue=qfd.maximum  
        numrange.append(qfdValue)    
        rasterList.append('dem_%s,%s'%(i, qfdValue))  
        arcpy.AddMessage('dem_%s,%s'%(i, qfdValue))  
    arcpy.ResetProgressor()  
    return numrange,rasterList  
  
  
if __name__ == '__main__':  
    """DEM路径"""  
    global block ,DEM_Length,qfdparams  
    block=arcpy.GetParameterAsText(0)  
#    block=ur'D:\work2014\起伏度\SX_DEM\dem'  
      
    """DEM路径"""  
    outputPath=arcpy.GetParameterAsText(1)  
#    outputPath=ur'D:\work2014\起伏度\New Folder'  
    DEM_Length=float(arcpy.GetParameterAsText(2))  
#    DEM_Length=30   #DEM的精度，即格网大小，单位m  
    qfdparams=arcpy.GetParameterAsText(3)  
    arcpy.CheckOutExtension("Spatial")  
    arcpy.env.overwriteOutput=True  
  
  
    numrange,rasterList=getRasterList()  
      
    bestWindow,bestArea= calcReliefAmplitude(numrange,DEM_Length).calcBestWindow()  
    output=os.path.join(outputPath,ur'dem_%s'%(bestWindow))  
    f=open(os.path.join(outputPath,'qfd.csv'),'a')  
    f.write('Raster,%s\n'%qfdparams)  
    for line in rasterList:  
        f.write(line+'\n')  
    f.close()  
    qfdResult=get_qfdRast(bestWindow)  
    qfdResult.save(output+'.img')  
    arcpy.AddMessage(u'The Best Window is %s*%s,Area is %s m²'%(bestWindow,bestWindow,bestArea))  

工具下载链接：https://pan.baidu.com/s/1zid7pYmZ2dOuSHV_xJqVcw 密码：mk65