[GIS] Measuring distance between two points considering elevation raster

coordinatesdistancepyqgisqgisr

enter image description here

Is it possible to compute the distance between 2 coordinates (x1,y1,z1) – (x2,y2,z2) considering an elevation raster (cells containing elevation)? On the image, there are 2 coordinates (red crosses). I would like the green distance. I know how to get the direct distance taking the elevation into consideration (red dashed line on the picture), but this does not take that there can be a mountain or a valley between the 2 points into account. I am using R but if there is a nice way to do it in QGIS that would be great too.

Best Answer

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)) )
Related Question