topic dem analysis in ArcGIS Spatial Analyst Questions

dem analysis

spencersawaske — Wed, 06 Nov 2013 19:43:12 GMT

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

Re: dem analysis

XanderBakker — Fri, 10 Dec 2021 20:18:16 GMT

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

Re: dem analysis

spencersawaske — Mon, 11 Nov 2013 13:42:24 GMT

Xander,

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

-spencer

Re: dem analysis

XanderBakker — Mon, 11 Nov 2013 13:45:21 GMT

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

Re: dem analysis

TimBarnes — Thu, 14 Nov 2013 01:29:16 GMT

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!

Re: dem analysis

808707 — Wed, 20 Nov 2013 10:57:57 GMT

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?

Re: dem analysis

XanderBakker — Fri, 10 Dec 2021 20:18:20 GMT

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

Re: dem analysis

808707 — Thu, 21 Nov 2013 14:05:30 GMT

Great thanks Xander. I'll give it a go and let you know how it goes.

Re: dem analysis

808707 — Fri, 29 Nov 2013 10:11:50 GMT

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'

Re: dem analysis

XanderBakker — Fri, 29 Nov 2013 11:42:39 GMT

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

Re: dem analysis

808707 — Fri, 29 Nov 2013 13:47:53 GMT

Ahh great thanks for the help Xander. I'll download the 32 bit version and see if that works.

Re: dem analysis

808707 — Fri, 10 Dec 2021 20:18:22 GMT

Hello again. It looks like I got it working. I installed the 32 bit PyScripter and got the same error again. However I realised my DEM had a cell type of float and converted it to an integer and it seems to work okay with a cell size of 10m.

However I would like to see if this will work with a cell size of 0.2m. I tried this and a memory allocation error occurs. The same happened with 1m. Is there a limit to the resolution of the raster?

This got me thinking to this part of the code

# 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

I was wondering how you calculated that the window size will be 41? and the winhalf will equal 20?

Re: dem analysis

XanderBakker — Mon, 02 Dec 2013 13:25:20 GMT

Hi Mark,

Basically the initial idea was a search distance of 200m (actually 200-400m, but I wanted to keep the array small). If you have a pixels size of 10m, a distance of 200m will be 20 pixels. So twenty pixels on each side means 20 * 2 = 40. When using filters and similar calculations, one uses an impair number of pixels to have a central pixel (the one where the value will be assigned). This is why the size is 41 pixels (2 * 20 + 1) and the winhalf is 20 pixels.

Kind regards,

Xander

Re: dem analysis

808707 — Mon, 02 Dec 2013 13:48:32 GMT

Hello Xander. Thanks very much for the information. I'll play around with this some more and see what I can do.