Calculating Driving Distances using sfnetworks

[swaheera]

I am trying to calculate the driving distance (e.g. based on road networks) between two sets of coordinates.

For example:

• CN Tower: 290 Bremner Blvd, Toronto, ON M5V 3L9 (-79.61203, 43.68312)
• Toronto Airport: 6301 Silver Dart Dr, Mississauga, ON L5P 1B2 (-79.61203, 43.64256)

To solve this problem, I tried to download the shapefile for the Canadian Road Network and subset it for the Province of Ontario:

    library(sf)
    library(rgdal)
    library(sfnetworks)
    # Set the URL for the shapefile
    url <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"
    
    # Create a temporary folder to download and extract the shapefile
    temp_dir <- tempdir()
    temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")
    
    # Download the shapefile to the temporary folder
    download.file(url, temp_file)
    
    # Extract the shapefile from the downloaded zip file
    unzip(temp_file, exdir = temp_dir)
    
    # Read the shapefile  and subset to Ontario
    a = st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), query="select * from lrnf000r22a_e where PRUID_R ='35'")

Then, by consulting different references (e.g. https://cran.r-project.org/web/packages/sfnetworks/vignettes/sfn01_structure.html) - I tried to calculate the distance between these two points:

  # convert the shapefile to an sfnetwork object
   net <- as_sfnetwork(a)
   
   # define your start and end points
   q1 <- st_point(c(-79.61203, 43.68312))
   q2 <- st_point(c(-79.38709, 43.64256))
   
   # set the CRS of the points to match the CRS of your shapefile
   q1 <- st_sfc(q1, crs = st_crs(a))
   q2 <- st_sfc(q2, crs = st_crs(a))
   
   # find the shortest path between the two points
   path <- st_network_paths(net, q1, q2)
   
   # calculate the distance of the path in meters
   distance <- sum(st_length(path))

But I get the following error:
Error in UseMethod("st_geometry") : no applicable method for 'st_geometry' applied to an object of class "c('tbl_df', 'tbl', 'data.frame')"

Can someone please show me how to fix this problem?
Thanks!

[loreabad6]

Hi @swaheera,

After reviewing your code I see a couple of issues.

1. The way you define your start and end points is not entirely correct. When you set the CRS, you are basically telling sf to assign a EPSG:3347 to coordinates that actually are in EPSG:4326 as given above. To correct this you should do something like:

q1 <- st_transform(st_sfc(q1, crs = st_crs(4326)) , st_crs(a))
q2 <- st_transform(st_sfc(q2, crs = st_crs(4326)) , st_crs(a))

2. When you create the network as_sfnetwork() defaults to directed networks. This usually is a tricky setting knowing that many road networks are not correctly digitized. The algorithm doesnot find a path between your points because of the direction of the edges. To avoid this you can use:

net <- as_sfnetwork(a, directed = FALSE)

3. In your original code you were getting a result already stored on the paths object. This is where all the nodes and edges of you resulting path are located as list objects. If you are only looking for the distance however, I would recommend you to call st_network_cost(net, q1, q2). This will directly give you the shortest network distance in meters between your two points.

If hope you don't mind that I convert this issue to a discussion. If you have any question we can discuss there 😄

[swaheera]

@loreabad6 : thank you so much for your answer! I really appreciate it!

[swaheera]

@loreabad6 : Based on the comments and suggestions you kindly provided, would this be the correct way to do this?

library(sf)
library(rgdal)
library(sfnetworks)
library(igraph)
library(dplyr)
library(tidygraph)

# Set the URL for the shapefile
url <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"

# Create a temporary folder to download and extract the shapefile
temp_dir <- tempdir()
temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")

# Download the shapefile to the temporary folder
download.file(url, temp_file)

# Extract the shapefile from the downloaded zip file
unzip(temp_file, exdir = temp_dir)

# Read the shapefile using the rgdal package
a = st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), query="select * from lrnf000r22a_e where PRUID_R ='35'")

#############################################

  # convert the shapefile to an sfnetwork object
  net <- as_sfnetwork(a, directed = FALSE)

 # define your start and end points
   q1 <- st_point(c(-79.61203, 43.68312))
   q2 <- st_point(c(-79.38709, 43.64256))
   
q1 <- st_transform(st_sfc(q1, crs = st_crs(4326)) , st_crs(a))
q2 <- st_transform(st_sfc(q2, crs = st_crs(4326)) , st_crs(a))

st_network_cost(net, q1, q2)

Units: [m]
         [,1]
[1,] 22709.69

Is this correct?
Also - is there some way to visualize this analysis?

[loreabad6]

Yes, that should work. As for further info and visualizing, I suggest you take a look at the routing vignette https://luukvdmeer.github.io/sfnetworks/articles/sfn04_routing.html

[swaheera]

Some other comments:

[swaheera]

Here is the comparison with OSRM:

library(tmap)

library(tmaptools)
remotes::install_github("riatelab/osrm")

q1 = geocode_OSM("6301 Silver Dart Dr, Mississauga, ON L5P 1B2")
q2 = geocode_OSM("290 Bremner Blvd, Toronto, ON M5V 3L9")

q1 = as.numeric(q1$coords)
q2 = as.numeric(q2$coords)

q1_lat = q1[1]
q1_long = q1[2]
q2_lat = q2[1]
q2_long = q2[2]

route = osrmRoute(src = c(q1[1], q1[2]) ,  dst = c(q2[1], q2[2]), osrm.profile = "car")

> route$distance
[1] 26.2836

[swaheera]

There appears to be some differences in the path lengths from OSRM approach vs the shapefile approach:

st_network_cost(net, q1, q2)

Units: [m]
         [,1]
[1,] 22709.69

#############

route = osrmRoute(src = c(q1[1], q1[2]) ,  dst = c(q2[1], q2[2]), osrm.profile = "car")

> route$distance
[1] 26.2836

[loreabad6]

Hi @swaheera, yes probably the answers are different given taht osrm uses profiles for their computaiton. As you have set it to car, it must be using only roads that allow car usage. With sfnetworks you use the whole network. You can of course also set this profiles for yourself by subsetting the network and so on. A bit more fidling with the package will give you more possibilities.