If you don't have this capability built into your GIS, but you can perform some basic grid operations ("map algebra"), there is still a solution.
The calculation comes down to finding the slope of the route at every point. If you knew this exactly, with no discretization error, you would integrate the secant of the slope. On a grid, the integral is estimated by obtaining the average of the secant for the cells intercepted by the route and multiplying the average by the route's length. (In map algebra-speak, that would be a "zonal average" multiplied by the route length.)
The slope of the route is not the same as the slope of the DEM! It depends on exactly how the route cuts across the surface. Thus, you need full information about the surface's "direction," which can be described in terms of strike and dip, slope and aspect, or by means of a unit normal vector (i.e., a 3D vector field perpendicular to the surface). The most reliable way is to reduce the problem to one where you know the normal vector field. This means you have a triple of numbers at every cell--represented as three separate grids, of course--which I will call (Nx, Ny, Nz). The direction of the route (in the plane) can be represented as a unit vector (x, y, t) where (x, y) gives its direction on the map. The value of t is the "rise" in the vertical direction: it tells us the rate at which the route must rise in order to remain on the surface. Thus, because the 2D speed of the route--its "run"--equals Sqrt(x^2 + y^2), the slope is given by
(1) tan(slope) = rise / run = t / Sqrt(x^2 + y^2).
In the calculations t will be a grid but the denominator, Sqrt(x^2+y^2), is just a number. If you calculate it using formula (4) below, it will equal 1, so you can forget about it: t will be the tangent of the route's slope grid and sec(slope) = sqrt(1 + t^2) will be the grid whose zonal average you compute.
It's easy to find t. By definition, the direction vector (x, y, t) is perpendicular to the normal vector. This means
0 = x*Nx + y*Ny + t*Nz, so
(2) t = -(x*Nx + y*Ny) / Nz.
In the calculation, Nx, Ny, and Nz are grids but x and y are numbers. Therefore t is a grid, as intended. (There won't be any trouble with the division, because it's not possible for Nz = 0: that would be a perfectly vertical cliff, which cannot be represented on a DEM.)
So: how do you find the normal vector (Nx, Ny, Nz) and the direction vector (x, y)? Typically a GIS will compute slope (s) and aspect (a) grids from the DEM. Express each as an angle. These are basically spherical coordinates for the unit normal vector. For aspects east of north, the unit normal is obtained by the usual spherical-to-cartesian coordinate conversion,
(3) (Nx, Ny, Nz) = (sin(s)*sin(a), sin(s)*cos(a), cos(s)).
In this calculation s and a are grids, so it describes three separate map algebra expressions to create three grids Nx, Ny, and Nz.
As a check, note that when the slope is zero (s=0), the normal vector is (0,0,1), pointing straight up, as it should. When the aspect is zero, the normal vector is (0, sin(s), cos(s)) which evidently points towards the north (y direction) and tilts from the vertical by an angle of s, which implies the surface tilts from the horizontal by an angle of s: that indeed is its slope.
Finally, let the bearing of the route be b (a constant angle, east of north). Its direction vector is
(4) bearing = (x, y) = (sin(b), cos(b)).
Note that the bearing is a pair of numbers, not a pair of grids, because it describes the direction of the route.
As the resolution of the DEM increases, you can observe more local variation in the slopes, causing the estimated slope to increase, as @johanvdw notes. I have studied this phenomenon by successively coarsening high-resolution DEMs and by comparing DEMs of one area obtained from different source. I found that in high-slope areas the differences in slope estimates can be substantial. These will translate to substantial differences in overland route length estimates. Otherwise, in uniformly low-slope areas the differences might be of no consequence.
One way you can assess the effect of resolution for your DEM is to perform a similar study. It takes little effort. For example, estimate the overland length of a route using the DEM, then re-estimate the length after aggregating that DEM into 2 x 2 blocks (coarsening by a factor of 2). If there is an inconsequential difference between the two estimates you should be fine; if the difference matters, then it might be worthwhile obtaining a finer-resolution DEM for your work. (There are more sophisticated methods to improve the slope and length estimates by exploiting the DEM you have, but it would take me too long to describe them here.)
Here is how to go in the other direction, computing azimuth, elevation based on known receiver and satellite locations (latitude, longitude, altitude). If this can be done, then given many azimuth and elevation data points (and given satellite location) one could solve for unknown receiver location using nonlinear optimization.
The best explanation on how to compute azimuth and elevation is given here. The implementation shown is in Pascal, but the pyorbital package wonderfully converted this code into Python. I pulled out the relevant code and am including it here:
import datetime
import numpy as np
F = 1 / 298.257223563 # Earth flattening WGS-84
MFACTOR = 7.292115E-5
EPS_COS = 1.5e-12
F = 1 / 298.257223563 # Earth flattening WGS-84
A = 6378.137 # WGS84 Equatorial radius
def jdays2000(utc_time):
"""Get the days since year 2000.
"""
return _days(utc_time - datetime(2000, 1, 1, 12, 0))
def jdays(utc_time):
"""Get the julian day of *utc_time*.
"""
return jdays2000(utc_time) + 2451545
def _fdays(dt):
"""Get the days (floating point) from *d_t*.
"""
return (dt.days +
(dt.seconds +
dt.microseconds / (1000000.0)) / (24 * 3600.0))
_vdays = np.vectorize(_fdays)
def _days(dt):
"""Get the days (floating point) from *d_t*.
"""
try:
return _fdays(dt)
except AttributeError:
return _vdays(dt)
def gmst(utc_time):
"""Greenwich mean sidereal utc_time, in radians.
As defined in the AIAA 2006 implementation:
http://www.celestrak.com/publications/AIAA/2006-6753/
"""
ut1 = jdays2000(utc_time) / 36525.0
theta = 67310.54841 + ut1 * (876600 * 3600 + 8640184.812866 + ut1 *
(0.093104 - ut1 * 6.2 * 10e-6))
return np.deg2rad(theta / 240.0) % (2 * np.pi)
def observer_position(time, lon, lat, alt):
"""Calculate observer ECI position.
http://celestrak.com/columns/v02n03/
"""
lon = np.deg2rad(lon)
lat = np.deg2rad(lat)
theta = (gmst(time) + lon) % (2 * np.pi)
c = 1 / np.sqrt(1 + F * (F - 2) * np.sin(lat)**2)
sq = c * (1 - F)**2
achcp = (A * c + alt) * np.cos(lat)
x = achcp * np.cos(theta) # kilometers
y = achcp * np.sin(theta)
z = (A * sq + alt) * np.sin(lat)
vx = -MFACTOR*y # kilometers/second
vy = MFACTOR*x
vz = 0
return (x, y, z), (vx, vy, vz)
def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt):
"""Calculate observers look angle to a satellite.
http://celestrak.com/columns/v02n02/
utc_time: Observation time (datetime object)
lon: Longitude of observer position on ground
lat: Latitude of observer position on ground
alt: Altitude above sea-level (geoid) of observer position on ground
Return: (Azimuth, Elevation)
"""
(pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = observer_position(
utc_time, sat_lon, sat_lat, sat_alt)
(opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \
observer_position(utc_time, lon, lat, alt)
lon = np.deg2rad(lon)
lat = np.deg2rad(lat)
theta = (gmst(utc_time) + lon) % (2 * np.pi)
rx = pos_x - opos_x
ry = pos_y - opos_y
rz = pos_z - opos_z
sin_lat = np.sin(lat)
cos_lat = np.cos(lat)
sin_theta = np.sin(theta)
cos_theta = np.cos(theta)
top_s = sin_lat * cos_theta * rx + \
sin_lat * sin_theta * ry - cos_lat * rz
top_e = -sin_theta * rx + cos_theta * ry
top_z = cos_lat * cos_theta * rx + \
cos_lat * sin_theta * ry + sin_lat * rz
az_ = np.arctan(-top_e / top_s)
az_ = np.where(top_s > 0, az_ + np.pi, az_)
az_ = np.where(az_ < 0, az_ + 2 * np.pi, az_)
rg_ = np.sqrt(rx * rx + ry * ry + rz * rz)
el_ = np.arcsin(top_z / rg_)
return np.rad2deg(az_), np.rad2deg(el_)
On Android using this package, I recorded some satellite values, one example is 22:0.0:331.0:7.0, the values are PRN:SNR:azimuth:elevation. For that PRN, I took the position of the satellite from here (updated regularly). I calculate,
millis = 1492964228447
from datetime import datetime
dt = datetime.fromtimestamp(millis/1000.0)
sat_lon = -109.861429
sat_lat = -6.725577
sat_alt = 20*1000*1000
lon = 13.442383333333332
lat = 52.483086666666665
alt = 0
az,el = get_observer_look(sat_lon, sat_lat, sat_alt, dt, lon, lat, alt)
print az, el
The result was
293.517066155 -25.1656865911
Not sure what to do with negative elevation but the azimuth seemed close to recorded values on Android. I checked few others, they were okay.
Going in the other direction would be an optimization problem. But of course cleaner solutions are probably possible.
Best Answer
Short answer, use the DEM.
Handheld GPS units of the type you mention (consumer grade), can potentially use two different methods to figure elevation - the calculated GPS position or an internal barometric altimeter. Because of the way GPS works (explained at the page mkennedy originally linked to in a comment, and elsewhere), without high-grade, corrected systems, the barometric method is typically more accurate than GPS. The next step up from Consumer is Differential, which according to that page has a vertical accuracy of 2-3x horizontal, and above that is Survey Grade which can get down to mm level depending on mode and range factors. Differential and Survey do use GPS calculation rather than barometric.
However the barometric method still has issues - readings can be affected by weather as the air pressure and temperature changes. You're unlikely to get consistently detailed or accurate elevation readings. I use a Garmin 60CSx and could provide multiple GPS track examples of this; I have tracks that are loops yet show a 200ft vertical difference between the start and end points even though they're the same place, as well as multiple versions of the same track with completely different elevations throughout. Some documentation on their Edge units puts vertical accuracy with calibrated barometric readings at +/-50 feet.
A couple of related links from Garmin's FAQ:
And there are many other places on the web you can find information/discussions about GPS calculated vs barometric altimeter accuracy for elevation readings, particularly related to handheld consumer units. Bottom line, every unit will return different values, and even the same unit will on different trips.
On the other hand, the DEM is basically a mapped value you're comparing the route against. You will have inaccuracies based on both the horizontal and vertical resolution of the DEM (so higher res is better), but it will be consistent from track to track (assuming said tracks are within a given horizontal tolerance of course). A DEM or map with contours or any other such data source is basically fixed - the values are what they were when measured/interpolated. But each GPS track is going to be whatever the unit recorded at whatever time, so it will differ from the map value at the same point. This is readily apparent on GPS units which are capable of displaying map data; you can see the elevation the map says you are at along with what the altimeter and potentially even GPS calculated elevation values are (and all three will probably differ).