Recently, I have been helping Jessica E. Marsh to develop a method of calculating distances along linear features. Specifically, we were measuring the distances between salmonid redd locations and juvenile survey locations on the river Frome, Dorset, UK.
Together, we devise a generalisable method to calculate the least cost distance between points located on a complex river network.
Our method includes the following steps:
- convert the river polygon dataset to a high resolution raster;
- thin that high resolution raster to represent the river as a network of single pixel width linear features;
- shift the sample site coordinates to the closest river pixel in the thinned raster; and
- calculate the least cost distance between a set of coordinates.
Our analysis is quite Frome-specific and so rather than outline the whole procedure Jess and I decided to present a contrived example to illustrate the last and most important step of the analysis: step 4.
For the remainder of the steps, we used a combination of R and GRASS GIS
An example
We provide a small example script with reproducible data (see downloads below) to run step 4.
It calculates distances along the local roads from the FBA River Laboratory (our address) to the local shop and pub.
R script
# start clean
rm(list = ls())
# start the clock
timer <- proc.time()
# libraries
require(rgdal)
require(raster)
require(rasterVis)
require(gdistance)
# read in EastStoke to Wool road network; in your current path
v <- readOGR('EastStoke_to_Wool_roads.shp', 'EastStoke_to_Wool_roads')
# add "val" field as 1 everywhere
v@data$val <- 1
# get extent
ex <- extent(v)
# make empty raster at 1m resolution
r <- raster(ex, res = 0.0001, crs = proj4string(v))
# rasterize road network
rp <- rasterize(v, r, 'val')
# read in stops
stops <- readOGR('esw.shp')
# find least cost distances for stops
tr <- transition(rp, function(x) 1 / mean(x), 8)
tr1 <- geoCorrection(tr)
sl <- shortestPath(tr1, origin = stops[1, ], goal = stops[-1, ], output = 'SpatialLines')
# prepare plot
# x11(width = 16, height = 10)
par(mfrow = c(1, length(sl)))
# plot path loop
for(i in 1:length(sl)){
## plot raster
plot(rp, legend = FALSE)
## add road network
plot(v, add = TRUE)
## add stops to plot
points(stops[c(1, i + 1), ], pch = 16, col = 'red')
## add shortest path
plot(sl[i], col = 'blue', add = TRUE)
}
# get path lengths
l <- SpatialLinesLengths(sl)
print(l)Download the data here.
The final result is shown on a map at the start of this post.
You can check the calculations against those given by Google Maps here
We hope this is of interest to someone. If you have any questions, then contact Jess Marsh or I by email.
Reuse
Citation
@online{d. gregory2017,
author = {D. Gregory, Stephen},
title = {Calculate Distances Along a River in {R} and {GRASS}},
date = {2017-01-31},
url = {https://stephendavidgregory.github.io/posts/least-cost-dist-along-river},
langid = {en}
}