Hi Spencer,I think that using numpy should be the way to go here. There is not much help on using numpy in ArcGIS, but I came across a nice blog post that explains some of the basics. I tried to calculate the E1 and E2 (see script and results below).import arcpy
from numpy import *
# based on
blog post:
# http://shreshai.blogspot.nl/2011/07/utilizing-numpy-to-perform-complex-gis.html
# input and output rasters
inputDEM = r'C:\Project\_Forums\numpy\grd\dem'
outRasE1 = r'C:\Project\_Forums\numpy\grd\dem01E1'
outRasE2 = r'C:\Project\_Forums\numpy\grd\dem01E2'
outRasAlt = r'C:\Project\_Forums\numpy\grd\dem01alt'
# determine dimensions of array
myArray = arcpy.RasterToNumPyArray(inputDEM)
[rows,cols]= myArray.shape
# cellsize = 10m, distance of 200m = 20, window size will be 41
# half the window size (20 = 20 pixels * 10m/pixel = 200 meters, moving window = 420 * 420)
winhalf = 20
# provide value in case of division by zero
infE1 = (winhalf * 2 + 1)^2
infE2 = (winhalf * 2 + 1)^2 * myArray.max()
IE1=zeros((cols,rows))
IE2=zeros((cols,rows))
Ialt=zeros((cols,rows))
cnt = 0
for r in range(0,rows):
if r % 10 == 0:
print "Processing row:{0}".format(r)
for c in range (1,cols):
cnt+=1
value = myArray[r,c]
r1 = r - winhalf
if r1 < 0:
r1 = 0
r2 = r + winhalf
if r2 > rows:
r2=rows
c1 = c - winhalf
if c1 < 1:
c1 = 1
c2 = c + winhalf
if c2 > cols:
c2 = cols
# size is number of pixels in moving window
size = (r2-r1+1) * (c2-c1+1)
# E1= (# cells in surrounding area with elev.< cell elev.)/(# cells in surrounding area with elev.> cell elev.)
data = myArray[r1:r2+1,c1:c2+1]
indicesGT = where(data>value,ones(data.shape),zeros(data.shape))
indicesLT = where(data<value,ones(data.shape),zeros(data.shape))
outGT = sum(indicesGT)
outLT = sum(indicesLT)
if outGT == 0:
outValE1 = infE1
else:
outValE1 = outLT / outGT
# E2= (sum[cell elev.-(elev.of surrounding cells<gage elev.)])/(sum[(elev.of surrounding cells>cell elev.)-cell elev.])
indicesGTE2 = where(data>value,data - (ones(data.shape)*value),zeros(data.shape))
indicesLTE2 = where(data<value,(ones(data.shape)*value) - data,zeros(data.shape))
outGTE2 = sum(indicesGTE2)
outLTE2 = sum(indicesLTE2)
if outGTE2 == 0:
outValE2 = infE2
else:
outValE2 = outLTE2 / outGTE2
# alternative output value could be:
# percentage of pixels in neighborhood with value higher than pixel value
indices = where(data>value,ones(data.shape),zeros(data.shape))
outVal = sum(indices)
outPerc = outVal * 100 / size
# assign values to output array
IE1[c,r] = outValE1 # E1
IE2[c,r] = outValE2 # E2
Ialt[c,r] = outPerc
# transpose
I2E1 = IE1.transpose()
I2E2 = IE2.transpose()
I2Alt = Ialt.transpose()
# retrieve meta from input raster
descData=arcpy.Describe(inputDEM)
cellSize=descData.meanCellHeight
extent=descData.Extent
spatialReference=descData.spatialReference
pnt=arcpy.Point(extent.XMin,extent.YMin)
#E1
newRasterE1 = arcpy.NumPyArrayToRaster(I2E1,pnt, cellSize,cellSize)
newRasterE1.save(outRasE1)
#E2
newRasterE2 = arcpy.NumPyArrayToRaster(I2E2,pnt, cellSize,cellSize)
newRasterE2.save(outRasE2)
#Alternative
newRasterAlt = arcpy.NumPyArrayToRaster(I2Alt,pnt, cellSize,cellSize)
newRasterAlt.save(outRasAlt)
print "ready..."
Result of analysis E1[ATTACH=CONFIG]29001[/ATTACH]In E2 there will be an nasty effect at the bottom and right side, since no compensation is applied for the reduction in the moving window size.Result of alternativeAs an alternative it might be interesting to look at the percentage of pixels in the neighborhood with a cell value higher than the value of the central pixel. This shows drainage with dark red colors and ridges in green. Don't know if this is useful...[ATTACH=CONFIG]29002[/ATTACH]Some useful links:http://shreshai.blogspot.nl/2011/07/utilizing-numpy-to-perform-complex-gis.htmlhttp://wiki.scipy.org/Tentative_NumPy_TutorialRasterToNumPyArray (arcpy)NumPyArrayToRaster (arcpy)Kind regards,Xander