Sun 18 March 2012
Ever try to figure out what the average aspect of an area is? i.e.
What direction does this hillside face?
Let's say we want to determine the average elevation of an area based on a raster DEM. Just take the arithmetic mean of all the elevation cells contained in the area - a simple zonal statistics problem.
Turns out that aspect is not quite as straightforward. True, we can easily use gdaldem to create an aspect map.
gdaldem aspect elevation.tif aspect.tif
This gives a raster with values in degrees: 0 is north, 90 is east, 180 is south, etc... but note that 360 is north as well. We're dealing with angular units, not linear units.
For example, take a nearly North facing hillside; the left edge is facing slightly NW (350 degrees) while the right edge faces slighty NE (10 degrees).
The arithmetic mean of the aspect values =
(350+350+10+10)/4 = 180°. Due south? That's entirely wrong! It doesn't take into account the angular units. For that we need to create grids representing the sin and cos of the aspect.
Luckily you can use the handy gdal_calc.py utility that comes with recent versions of gdal. This allows you to apply numpy's trigonometric functions to a raster...
gdal_calc.py -A aspect.tif --calc "cos(radians(A))" --format "GTiff" --outfile cos_aspect.tif
gdal_calc.py -A aspect.tif --calc "sin(radians(A))" --format "GTiff" --outfile sin_aspect.tif
Now we can look at the sum of the cos/sin grid cells for our area and take the arctangent according to this python code
avg_aspect_rad = math.atan2(sum(cos_cells), sum(sin_cells))
avg_aspect_deg = math.degrees(avg_aspect_rad)
In our example avg_aspect_deg comes out to an aspect of 0 degrees (due north) which is exactly what we'd expect.
Thanks to Dan Patterson for his forum post which clued me into this approach.