In flow accumulation modeling, a common bit of data needed is the average upslope value of a parameter of interest, such as slope or curvature. In most cases, this is simply a matter of calculating a flow accumulation weighted by the parameter in question, then dividing this by the unweighted flow accumulation. But what if you want average upslope aspect? Since aspect is measured in degrees, and a value of 0° and 360° are the same, a simple arithmetic mean will be useless. So, what to do? In this post, I will walk you through the process of calculating the mean upslope aspect using ArcGIS, and leave you with a working python script that will automate the process.
The first step is to break the aspect angle into measures of “Northness” and “Eastness.” You calculate these as the cosine and sine of the angles, respectively. Cosine values range from -1 at 180° (South) to 1 at 0° (North) and sine values from -1 at 270° (West) to 1 at 90° (East). Once you’ve calculated the sine and cosine, you have meaningful numbers to average. However, mean angles are not calculated like ordinary arithmetic means. To calculate the mean, we need to sum up the cosine and sine values individually, then take the arctangent of their ratio.
Now that we have the basic idea down, we need to implement it. The first step is to preprocess the aspect raster in preparation for the sin and cos tools in Spatial Analyst. The Aspect tool assigns flat pixels a value of -1, which we want to ignore in our average. Also, we want to convert the aspect values from degrees to radians. We can accomplish this with a single statement, either in Raster Calculator or Python:
arcpy.sa.SetNull(rawAspect,Times(rawAspect,0.01745329),"Value < 0")
Finding these sums poses its own unique problem: the weighted flow accumulation tool in ArcMap does not support negative numbers! If you were to run the weighted flow accumulation for the cosine and sine of the angles, we would end up with all positive sums, which in turn would only yield us mean angles between 0° and 90°, where the cosine and sine are both positive. So how do we get around this? The trick is to separate the cosine and sine rasters one step further into their positive and negative components, then use the absolute values of each for a weighted flow accumulation, and subtract the “negative” raster from the “positive” raster. To minimize the amount of coding I would need to do, I created a function to split the rasters appropriately and return the 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)
Note that this function requires an input flow direction raster, which you should already have if you are working with flow accumulation. We can now call this function on the cosine and sine rasters to get the weighted flow accumulation. The result should look something like this:
Now, we need to calculate the arctangent of the ratio of the sine and cosine sums. We’ll do this using ATan2, a special version of arctangent that takes two parameters and remembers the sign of each input number so that the mean is assigned to the appropriate quadrant. The order of parameter input varies between implementations of ATan2; in ArcGIS it takes sine as the first parameter and cosine as the second. The code will look like this:
ArcTan = ATan2(sinAccum,cosAccum)
and then you need to convert your mean aspects into degrees:
meanAspect = Mod(360+ArcTan*(180/math.pi),360)
*Note: This snippet works in ArcGIS v 10.1. However, I’ve noticed that it produces only a limited set of values in 10.0 SP5. This can be remedied by breaking it into two steps: one to convert it to degrees, and one to set the range back to 0-360°.
Now, you have your mean upslope aspect! However, this process will leave some NoData gaps in your raster where there were flats and in some other areas. We can clean this up by filling in the missing data with the values from the original aspect raster:
finalAspect = Con(IsNull(meanAspect),rawAspect,meanAspect)
There you have it! A beautiful average upslope aspect raster, ready to feed into your favorite analysis method! Some variations we could do are averaging aspect over a block using zonal statistics or over a radius using focal statistics. You can also easily modify the script to average other tricky terrain measurements, such as curvature. I’ve included the completed script below, ready to be attached to an ArcGIS Python Script tool. If you want to copy this script for your own use, just expand the codeblock below and copy the code. This script is fully tested to work with ArcGIS v. 10.1, but it should work with v. 10.0 as well.
import arcpy 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) fillAspect = Con(IsNull(meanAspect),rawAspect,meanAspect) fillAspect.save(outMean)