From a Different Angle: Averaging Upslope Aspect

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.

Raw aspect raster, in ArcMap's default Technicolor scheme.

Raw aspect raster, in ArcMap’s default Technicolor scheme.

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:

Completed flow accumulation weighted by cosine.

Completed flow accumulation weighted by cosine.

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)
Still Technicolor, though.

Mean upslope aspect raster. Much nicer!

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)

About these ads
This entry was posted in How-To and tagged , , , , , . Bookmark the permalink.

4 Responses to From a Different Angle: Averaging Upslope Aspect

  1. Brian Gelder says:

    Many thanks for the idea of decomposing aspect into ‘eastness’ and ‘northness’ – I was calculating aspect range using radians and this will let me extend the concept to other statistics.

  2. Joe says:

    I attached the script to an ArcGIS Python tool, but when I tried to run it, I got the following error on line 4 (from arcpy.sa import *):
    “SyntaxError: invalid syntax (MeanUpslopeAspect.py, line 4) Failed to execute (MeanUpslopeAspect).”

    • jcguarneri says:

      Joe: I’m sorry it’s not working for you. I would recommend trying to run “from arcpy.sa import *” in a new Python window in ArcGIS to see if it works correctly (if it doesn’t, there’s something bigger wrong than the script). After that, I would double-check the syntax of the lines you had to add to make the Python tool.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s