DEM Processing

See scripts/process_dems.sh

Handling aspects

$ gdaldem aspect elevation.tif aspect.tif
$ 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

http://forums.esri.com/Thread.asp?c=93&f=995&t=211300

Python code to calculate true average aspect (arctan2(sum(cos), sum(sin))

import math

# 9 cells, 3 facing NNW, the others NNE 
# average should be roughly north
aspects_deg = [357,357,358,45,21,13,55,28,24]
aspects_deg = [float(x) for x in aspects_deg]
print sum(aspects_deg)/9.0
# mean aspect of 139.7 (SE) is DEFINITELY NOT RIGHT

# convert to radians 
aspects_rad = [math.radians(x) for x in aspects_deg]
cos = [math.sin(x) for x in aspects_rad]
sin = [math.cos(x) for x in aspects_rad]

avg_aspect_rad = math.atan2(sum(cos), sum(sin))
avg_aspect_deg = math.degrees(avg_aspect_rad)
print avg_aspect_deg
#avg_aspect_deg == 19.625 .. CORRECT!