Shortest Path Calculations in R

[omario]

I am working with the R programming language.

Based on a geospatial shapefile, I am trying to find out the driving distance between two coordinates.

First I loaded the shapefile:

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
# 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'")

The shapefile 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 I tried to calculate the distance by treating the road network as a graph:

# Convert the shapefile to an sfnetwork object
net <- as_sfnetwork(a)

# Define the two points of interest in WGS84 (EPSG:4326)
pt1 <- st_point(c(-79.61203, 43.68312))
pt2 <- st_point(c(-79.61203, 43.64256))

# Set the CRS of the points to WGS84 (EPSG:4326)
pt1 <- st_sfc(pt1, crs = 4326)
pt2 <- st_sfc(pt2, crs = 4326)

# Transform the points to the CRS of the network
pt1_transformed <- st_transform(pt1, st_crs(net))
pt2_transformed <- st_transform(pt2, st_crs(net))

# Find the nearest points on the network to the transformed points of interest
n1 <- st_nearest_feature(pt1_transformed, net)
n2 <- st_nearest_feature(pt2_transformed, net)

Everything seems to work until the last part:

# Calculate the shortest path between n1 and n2 on the network
path <- net %>%
    activate(edges) %>%
    mutate(weight = edge_length()) %>%
    as_tbl_graph() %>%
    shortest_paths(n1, n2) %>%
    .$epath[[1]]

# Calculate the distance of the shortest path in meters
distance <- sum(path$weight)

# Print the distance in kilometers
distance_km <- distance / 1000
cat("The network distance between CN Tower and Toronto Airport is", distance_km, "km.")

But this gave me the following error:

Error in .[[.$epath, 1]] : incorrect number of subscripts
In addition: Warning message:
In shortest_paths(., n1, n2) :
  At structural_properties.c:4745 :Couldn't reach some vertices

Can someone please show me how to fix this?

Thanks!

Note: I am specifically trying to solve this problem without using OSRM

[mikecrobp]

I have never used these packages before.
But maybe work back through the steps of that pipeline until you get something that does not error. And then work forward again.
And look at the reason for the warnings.