4

I am trying to write a python script to calculate the mean aspect upslope from each cell in a DEM. I have the general workflow down, but my final output is restricted to angles between 0 and 90 degrees, which is not correct for the input data. The general procedure to average an angle is:

  1. Calculate the sine and cosine of each angle.
  2. Sum the sines and cosines.
  3. Use atan2 on the sums to find the mean angle.

I have successfully used this technique to find focal and zonal mean aspects, and am trying to adapt my code to calculate upslope aspect. I am summing up the sine and cosine rasters using a weighted flow accumulation, which is the only difference between this and my zonal/focal workflows. I have included the script below.

Does anyone have any idea what's different about flow accumulation that causes this, or is there something else in my code that I'm missing?

arcpy.CheckOutExtension("spatial")
import math
from arcpy.sa import *

#collect input parameters
inDEM = arcpy.GetParameter(0)
flowDir = arcpy.GetParameter(1)
outMean = arcpy.GetParameterAsText(2)
scratch = arcpy.GetParameter(3)

#set workspace and snap raster
arcpy.env.workspace=scratch
arcpy.env.snapRaster = inDEM

#Calculate aspect and set flat cells to NoData
rawAspect = Aspect(inDEM)
nullFlat = SetNull(rawAspect,rawAspect,"Value < 0")

#convert aspect to radians and calulate cos/Sin
Radians = Times(nullFlat,0.01745329)
cosAsp = Cos(Radians)
sinAsp = Sin(Radians)

#sum upslope Cos/Sin rasters, use ATan2 to average
cosAccum=FlowAccumulation(flowDir,cosAsp)
sinAccum=FlowAccumulation(flowDir,sinAsp)
ArcTan = ATan2(sinAccum,cosAccum)

#convert mean aspect back to degrees and save output
meanAspect = Mod(360+ArcTan*(180/math.pi),360)
meanAspect.save(outMean)
Improve this question
7
  • 1
    Since output is in [0,90], that suggests both cosAccum and sinAccum are positive. Perhaps FlowAccumulation is ignoring negative values? (If so, you could work around this by flow accumulating cosAsp+1, sinAsp+1, and 1 and subtracting the latter from the first two, and then proceed as before.)
    – whuber
    Commented Feb 20, 2013 at 20:16
  • 1
    It does look like FlowAccumulation does not support negative values. Your comment made me think of another thing to search, and I came up with this: forums.esri.com/Thread.asp?c=93&f=995&t=243378 It appears I should try to separately accumulate the the negative and positive portions. I'm going to modify my script, and if it works, I'll post it as the answer.
    – Jay Guarneri
    Commented Feb 20, 2013 at 20:21
  • 1
    The work-around I suggested uses fewer calculations (just three flow accumulations instead of four; the other calculations are all much faster). On the other hand, it is a little more sensitive to loss of floating point roundoff precision than your suggestion.
    – whuber
    Commented Feb 20, 2013 at 20:27
  • I see what you're saying. I don't think I quite grasped it before. Effectively, you are proprosing a +1 transform, then a -1 back transform after the flow accumulation. I've just tried it out, but it produces a different result than the pos/neg approach. For example, Cos accumulation ranges from -14804.2 to 64299.2 in the pos/neg, but ranges from -16002.2 to 58620.1 with the +1 transform. I'm guessing this has something to do with what values are considered to be 0 and how flow accumulation handles these. But which one is right?
    – Jay Guarneri
    Commented Feb 20, 2013 at 20:49
  • 1
    I would trust the positive/negative approach but am concerned that there's a difference between the two calculations--only the last one or two decimal places should differ. One check worth running is to replace the +1 offset with, say, +2 and verify that you get the same result (up to FP error), for then there couldn't be any zero values to accumulate. That at least would verify the transform-back transform calculations are consistently offsetting each other. Perhaps the best thing to do now is test on some tiny grids having integral values, which you can verify manually.
    – whuber
    Commented Feb 20, 2013 at 20:58

1 Answer 1

Reset to default
5

It looks like the problem is that weighted flow accumulation does not support negative values. The workaround is to split the sin and cos rasters into positive and negative portions, run a weighted accumulation on these, and subtract the negative accumulation from the positive. The negative raster should contain the absolute values of the negative portion.
To split, accumulate, and recombine the cos and sin rasters, I defined a function that replaced the single weighted accumulation step. I also included a final step to fill in any NoData values on the mean aspect raster with the values from the original raw aspect raster. The new code is below:

Edit: I wrote up a blog post that goes into script build in more detail 

arcpy.CheckOutExtension("spatial")
import math
from arcpy.sa import *

#define function to calculate flow accumulation
def angleAccum(flowdir,angle):
    pos_angle = Con(angle,angle,0,"Value > 0")
    neg_angle = Abs(Con(angle,angle,0,"Value < 0"))
    pos_accum = FlowAccumulation(flowdir, pos_angle,"FLOAT")
    neg_accum = FlowAccumulation(flowdir, neg_angle,"FLOAT")
    return Minus(pos_accum,neg_accum)


#collect input parameters
inDEM = arcpy.GetParameter(0)
flowDir = arcpy.GetParameter(1)
outMean = arcpy.GetParameterAsText(2)
scratch = arcpy.GetParameter(3)

#set workspace and snap raster
arcpy.env.workspace=scratch
arcpy.env.snapRaster = inDEM

#Calculate aspect and set flat cells to NoData
rawAspect = Aspect(inDEM)
nullFlat = SetNull(rawAspect,rawAspect,"Value < 0")

#convert aspect to radians and calulate cos/Sin
Radians = Times(nullFlat,0.01745329)
cosAsp = Cos(Radians)
sinAsp = Sin(Radians)

#sum upslope Cos/Sin rasters, use ATan2 to average
cosAccum=angleAccum(flowDir,cosAsp)
sinAccum=angleAccum(flowDir,sinAsp)
ArcTan = ATan2(sinAccum,cosAccum)

#convert mean aspect back to degrees and save output
meanAspect = Mod(360+ArcTan*(180/math.pi),360)
finalAspect = Con(IsNull(meanAspect),rawAspect,meanAspect)
finalAspect.save(outMean)
Improve this answer
1
  • 2
    As a side note: While I was researching this technique, I originally saw the average being calculated as ATan2(cos,sin) as opposed to ATan2(sin,cos) Apparently, the confusion is due to which parameter counts as the denominator, which will vary between implementations of ATan2. Looking at my results, it appears as though ATan2(sin,cos) is the correct input for ArcGIS, but this underscores the importance of checking your results to see if they make sense, and checking the documentation of the software you are using.
    – Jay Guarneri
    Commented Feb 21, 2013 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.