dem analysis

3314
13
11-06-2013 11:43 AM
spencersawaske
New Contributor
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
0 Kudos
13 Replies
XanderBakker
Esri Esteemed Contributor
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
0 Kudos
spencersawaske
New Contributor
Xander,

Wow, thank you so much for putting so much time into this.  I am incredibly impressed.

-spencer
0 Kudos
XanderBakker
Esri Esteemed Contributor
Xander,

Wow, thank you so much for putting so much time into this.  I am incredibly impressed.

-spencer


Hi Spencer,

I'm glad I could be of any assistance. If you think it was helpful you   can use the "arrow" button in order to help other  members find useful  information:



More info here: http://resources.arcgis.com/en/help/forums-mvp/

Kind regards,

Xander
0 Kudos
TimBarnes
Occasional Contributor III
Xander,

Wow, thank you so much for putting so much time into this.  I am incredibly impressed.

-spencer


+1
I find your posts very helpful Xander 🙂 Kudos!
0 Kudos
808707
by
New Contributor III
Incidently E2 looks like it might come in very handy.

I've been wanting to learn python for ages but just never get the time.

How would I go about running this code?
0 Kudos
XanderBakker
Esri Esteemed Contributor
The code is a standalone script and you can run it from any Python IDE. You can also load this into the Python window in ArcMap. This is done with the Load command (context sensitive menu of the Python window).

Make sure you adjust the names and paths to the datasets before you load the script into the Python window. I'm referring to the following lines of code:

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'


Kind regards,

Xander
0 Kudos
808707
by
New Contributor III
Great thanks Xander. I'll give it a go and let you know how it goes.
0 Kudos
808707
by
New Contributor III
Hello. I was wondering if I could get a bit of help with this code?

I've been trying to run it in PyScripter.

I'm getting the following error message:
 
TypeError: unsupported operand type(s) for ^: 'int' and 'numpy.float64'
0 Kudos
XanderBakker
Esri Esteemed Contributor
Hello. I was wondering if I could get a bit of help with this code?

I've been trying to run it in PyScripter.

I'm getting the following error message:
       
TypeError: unsupported operand type(s) for ^: 'int' and 'numpy.float64'


Not sure if this helps, but in this thread:
http://forums.arcgis.com/threads/54169-PyScripter-Question

... they say that you need to use the 32 bits version of PyScripter.

Kind regards,

Xander
0 Kudos