使用arcpy求取地形起伏度的最佳统计单元

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)

其次，使用均值变点法寻找拐点。

"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

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

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))
    

源代码在此