Calculating mean upslope aspect from each cell in DEM using ArcPy?

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)

Answer

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)