python实现均值变点法计算地形起伏度

# -*- coding:utf8 -*- #------------------------------------------------------------------------------- # Name: D:/Working2013/TempCoding/calc_qifudu.py # Abstract: abstract # Description: description # # Created: 2014-7-14 # Author: Baihc # Contact: baihc@esrichina.com.cn # # Copyright: (c) 2013 Gistech All Rights Reserved. # License: license #------------------------------------------------------------------------------- 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))