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