Calculating Distances Between Points on a Shapefile

Question

I am working with the R programming language.

I am interested in learning about how 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 using the rgdal package
# source: https://gis.stackexchange.com/questions/456748/strategies-for-working-with-large-shapefiles/456798#456798
a = st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), query="select * from lrnf000r22a_e where PRUID_R ='35'")

When imported into R, the file looks something like this:

Simple feature collection with 570706 features and 21 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 5963148 ymin: 665490.8 xmax: 7581671 ymax: 2212179
Projected CRS: NAD83 / Statistics Canada Lambert
First 10 features:
   OBJECTID NGD_UID       NAME TYPE  DIR AFL_VAL ATL_VAL AFR_VAL ATR_VAL CSDUID_L                           CSDNAME_L CSDTYPE_L CSDUID_R                           CSDNAME_R CSDTYPE_R PRUID_L PRNAME_L PRUID_R PRNAME_R RANK
1         4 3434819       <NA> <NA> <NA>    <NA>    <NA>    <NA>    <NA>  3526003                           Fort Erie         T  3526003                           Fort Erie         T      35  Ontario      35  Ontario    5
2         5 1358252      South LINE <NA>    <NA>    <NA>    <NA>    <NA>  3551027                Gordon/Barrie Island        MU  3551027                Gordon/Barrie Island        MU      35  Ontario      35  Ontario    5
3         8 1778927       <NA> <NA> <NA>    <NA>    <NA>    <NA>    <NA>  3512054                           Wollaston        TP  3512054                           Wollaston        TP      35  Ontario      35  Ontario    5
4        11  731932     Albion   RD <NA>    <NA>    <NA>    2010    2010  3520005                             Toronto         C  3520005                             Toronto         C      35  Ontario      35  Ontario    3
5        18 3123617  County 41   RD <NA>     640     708     635     709  3511015                     Greater Napanee         T  3511015                     Greater Napanee         T      35  Ontario      35  Ontario    3
6        20 4814160       <NA> <NA> <NA>    <NA>    <NA>    <NA>    <NA>  3553005     Greater Sudbury / Grand Sudbury        CV  3553005     Greater Sudbury / Grand Sudbury        CV      35  Ontario      35  Ontario    5
7        21 1817031       <NA> <NA> <NA>    <NA>    <NA>    <NA>    <NA>  3537028                         Amherstburg         T  3537028                         Amherstburg         T      35  Ontario      35  Ontario    5
8        24 4825761       <NA> <NA> <NA>    <NA>    <NA>    <NA>    <NA>  3554094 Timiskaming, Unorganized, West Part        NO  3554094 Timiskaming, Unorganized, West Part        NO      35  Ontario      35  Ontario    5
9        25  544891     Dunelm   DR <NA>       1       9       2      10  3526053                      St. Catharines        CY  3526053                      St. Catharines        CY      35  Ontario      35  Ontario    5
10       28 1835384 Seven Oaks   DR <NA>     730     974     731     975  3515005             Otonabee-South Monaghan        TP  3515005             Otonabee-South Monaghan        TP      35  Ontario      35  Ontario    5
   CLASS                 _ogr_geometry_
1     23 LINESTRING (7269806 859183,...
2     23 LINESTRING (6921247 1133452...
3     23 LINESTRING (7320857 1089403..

Then, by consulting different references (e.g. https://cran.r-project.org/web/packages/sfnetworks/vignettes/sfn01_structure.html, Creating a route from a list of points and Overlaying the route (list of points) on the road segment/ road network in R) - 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!

Answer 1

You are very close to getting the desired result. The error is happening because st_network_paths function is expecting an sfnetwork object, but it is actually receiving a data.frame. The issue seems to be with the conversion of your "a" object to an sfnetwork object using as_sfnetwork() function.

Please try the following code to fix the issue:

library(sf)
library(rgdal)
library(sfnetworks)
library(units)

# Your previous code remains the same until this point

# convert the shapefile to an sfnetwork object
net <- as_sfnetwork(a, directed = FALSE)  # Set directed = FALSE if your road network is undirected

# Reproject the points to match the CRS of your shapefile
q1 <- st_transform(q1, st_crs(a))
q2 <- st_transform(q2, st_crs(a))

# Find the nearest node on the network for each point
from_node <- st_nearest_feature(q1, net)
to_node <- st_nearest_feature(q2, net)

# Find the shortest path between the two points
path <- st_network_paths(net, from_node, to_node)

# Calculate the distance of the path in meters
distance <- sum(path$length) %>% set_units("meters")

This should fix the error and give you the correct driving distance between the two points on the road network.

EDIT

Here's the revised code:

library(sf)
library(rgdal)
library(sfnetworks)
library(units)

# Your previous code remains the same until you read the shapefile into "a"

# Transform the data to WGS 84
a <- st_transform(a, 4326)

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

# Create the points with the same CRS as the shapefile
q1 <- st_point(c(-79.61203, 43.68312))
q2 <- st_point(c(-79.38709, 43.64256))
q1 <- st_sfc(q1, crs = 4326)
q2 <- st_sfc(q2, crs = 4326)

# Find the nearest node on the network for each point
from_node <- st_nearest_feature(q1, net)
to_node <- st_nearest_feature(q2, net)

# Find the shortest path between the two points
path <- st_network_paths(net, from_node, to_node)

# Calculate the distance of the path in meters
distance <- sum(path$paths[[1]]$length) %>% set_units("meters")

This should give you the correct distance between the two points. Note that the actual value will depend on the accuracy and completeness of the input road network data.