-
Notifications
You must be signed in to change notification settings - Fork 11
DEM Processing
Matthew Perry edited this page May 22, 2013
·
6 revisions
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!