Given two coordinates on a digital elevation model (DEM), how do I compute the actual distance traveled between the two locations, assuming a straight line route? How will changing the resolution of my DEM affect the results?
Add a comment
|
2 Answers
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.)
SAGA GIS has a module for this: Interactive profile
http://www.saga-gis.org/saga_modules_doc/ta_profiles/index.html
The resulting points will contain the distance and the overland distance. If the DEM has a coarser resolution your overland distance will always be a bit lower (unless you have strange border conditions), but in reality this difference is most likely not important. If the area is rather flat, even the overland distance and normal distance will be almost the same:if the slope between two points along you line is eg 20%, the overland distance will only be 2% higher than the normal distance (sqrt(1^2+0.2^2)=1.019 ).
