Extract routes between stations from a rail network shapefile

Question

Problem: I want to extract a route along the rail network between two stations from a shapefile and only plot this particular route, rather than the entire network.

This is what I have so far:

I have a shapefile with the entire rail network of the UK, plotted it looks like this:

library(maptools)
rail <- readShapeSpatial("railnetworkLine.shp")

I also have a list of stations with Eastings and Northings, for example:

 1) ABDARE 300400 202800
 2) DEIGHTN 416490  419140

I can add them to the map and it looks like this:

plot(rail)
plot(spdf.station, add=TRUE, col="red", pch=20)

So what I don't know, is how I can extract the route between them and just plot that route - the information is obviously in the shapefile and I have the coordinates of the station, but I don't understand how to extract it.

I managed to calculate the distance between them with this code:

SpacingInMetres <- 10000
require(secrlinear)
network <- read.linearmask(data=rail, spacing=SpacingInMetres)
distance <- (networkdistance (stations[1,], stations[2,], network))/1000

# Confirm distance:
distance
>311.7893

And I found that you can get the routes along roads with Google Maps with ggmaps (see here). But how can you do it when you have a shapefile as the network input rather than Google Maps?

I think maybe the packages 'shp2graph' + 'igraph' are useful, but I just can't figure it out. Any thoughts?

Answer 1

Shortest paths on route networks can be calculated using the stplanr package. I used a shapefile with the entire rail network for the Netherlands. This shapefile is available from:

https://mapcruzin.com/free-netherlands-arcgis-maps-shapefiles.htm

library(sf)
library(ggplot2)
library(stplanr)

# Read shapefile
nl_rails_sf <- sf::st_read("~/netherlands-railways-shape/railways.shp")

# Data frame with station locations
stations_df <- data.frame(station = c("Den Haag", "Den Bosch"), 
                          lat = c(52.080276, 51.690556), 
                          lon = c(4.325, 5.293611)) 

# Create sf object 
stations_sf <- sf::st_as_sf(stations_df, coords = c("lon", "lat"), crs = 4326)  

# Find shortest route
slnetwork <- SpatialLinesNetwork(nl_rails_sf)
find_nodes <- find_network_nodes(sln = slnetwork, 
                                 x = stations_df$lon, 
                                 y = stations_df$lat, 
                                 maxdist = 2e5)
route_dhdb_df <- data.frame(start = find_nodes[1], end = find_nodes[2])
route_dhdb_sf <- sum_network_links(sln = slnetwork, routedata = route_dhdb_df)

# Distance route in meters
distance_m <- sum(route_dhdb_sf$length) # 112189.5 [m]

# Plot results
ggplot(nl_rails_sf) +
    geom_sf() +
    theme_void() +
    geom_sf(data = stations_sf, color = "red") +
    geom_sf(data = route_dhdb_sf, color = "red")