dem analysis

spencersawaske

Hello,

I'm having trouble with what I think should be a relatively straightforward analysis of 10 meter dem data.  If you have a minute, could you offer any suggestions. 

I am trying to calculate a sort of "exposure" score for every cell in a raster.  The basic idea is to compare the value of each raster cell with those nearby(200-400m radius).   

E1= (# cells in surrounding area with elev.< cell  elev.)/(# cells in surrounding area with elev.> cell  elev.)
E2= (sum[cell elev.-(elev.of surrounding cells<gage  elev.)])/(sum[(elev.of surrounding cells>cell  elev.)-cell elev.])

For both indices, larger values represent locations with greater exposure and smaller values represent more protected, less exposed locations.  For instance, the E1 value for a location on top of a cone would be infinity, while that for the low point of an inverted cone would be 0.  I've been able to compute these indices for a small number of point locations using a combination of buffer and zonal tools.  Unfortunately, I'm now trying to perform this analysis for every cell in a much larger area(2-3 thousand cells) and my original method is no longer feasible.

I would greatly appreciate any suggestion you could offer.

Thanks,
spencer

XanderBakker

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 alternative
As 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.html
http://wiki.scipy.org/Tentative_NumPy_Tutorial
RasterToNumPyArray (arcpy)
NumPyArrayToRaster (arcpy)


Kind regards,

Xander

numpyDEM_alt.jpg

numpyDEM_E1.jpg