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?
DEM Analysis – Calculating Surface Distance Between Two Coordinates
demdistancescale
Related Solutions
With regards to your question of the difference betwen SRTM and EU-DEM the ESA quotes:
The EU-DEM is a hybrid product based on SRTM and ASTER GDEM data fused by a weighted averaging approach and it has been generated as a contiguous dataset divided into 1 degree by 1 degree tiles, corresponding to the SRTM naming convention.
As they mention ASTER I and the term "hybrid" I would skip this data set. ASTER had major issues as it is an optical system. I would work with SRTM-1 data as it has exetnesive coverage in terms of quality analysis and the physics and methods behind the data are "defined".
I've written some comparisons for a non-urban area: ASTER vs. SRTM 3 and ALOS vs. SRTM1 I hope this helps somehow!
In R this is fairly straight forward for two points. There are functions to return the distance matrix and for extracting the elevation for the points based on a DEM.
First add the required packages and, here create a dummy raster for example purposes. You can use the "raster" function to read your elevation raster from disk.
library(raster)
library(sp)
e <- as(extent(-1942755,-1847241,2524840,2624901), "SpatialPolygons")
proj4string(e) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
r <- raster(e, res=1000)
r[] <- runif(ncell(r),500,2500)
Now, using raster::extract
we can extract the elevation values for the two points of your path and return the uncorrected straight line distance between the points using spDists
.
xy <- list( coordinates( sampleRandom(r, 2, sp=TRUE) ) )
elev <- extract(r, SpatialPoints(xy))
d <- spDists( SpatialPoints(xy) )[2]
Using the data generated above, it is a simple matter of deriving the corrected distance.
elev.delta <- max(elev) - min(elev)
sqrt(d^2 + elev.delta^2)
Where it gets more complicated is if you want to correct the distance along the entire path, which in this case would be defined by the cell size of the raster underlying the straight line or, even an irregular path. You would need to set up a looping structure to calculate the corrected distances for each line segment.
A starting point is to coerce your points to a sp line object, calculate the length of the line using sp::SpatialLinesLengths
(same result as sp::spDists
) and extract all of the underlying elevation values with raster::extract
.
At this point you would be set up to loop through each elevation value to return the corrected distance for all elevation deltas to then add to the uncorrected distance [d]. This would effectively be slope corrected distance.
l <- SpatialLines(list(Lines(list(Line(xy)),ID="1")))
( d <- SpatialLinesLengths(l) )
( elev <- unlist(extract(r, l)) )
Best Answer
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.)