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:
- Calculate the sine and cosine of each angle.
- Sum the sines and cosines.
- 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)