Working with tbl_graphs using R

Question

I am working with the R programming language.

Based on a shapefile, I am trying to find out the driving distance between two coordinates (e.g. CN Tower and Toronto Pearson Airport).

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)

From here, I create a "tbl_graph":

path <- net %>%
    activate(edges) %>%
    mutate(weight = edge_length()) %>%
    as_tbl_graph()
> path
A tbl_graph: 430824 nodes and 570706 edges

A directed multigraph with 442 components

Edge Data: 570,706 x 25 (active)
from    to OBJECTID NGD_UID NAME      TYPE  DIR   AFL_VAL ATL_VAL AFR_VAL ATR_VAL CSDUID_L CSDNAME_L   CSDTYPE~ CSDUID~ CSDNAME_R  CSDTYP~ PRUID_L PRNAME~ PRUID_R PRNAME~ RANK  CLASS                  _ogr_geometry_ weight
  <int> <int>    <dbl> <chr>   <chr>     <chr> <chr> <chr>   <chr>   <chr>   <chr>   <chr>    <chr>       <chr>    <chr>   <chr>      <chr>   <chr>   <chr>   <chr>   <chr>   <chr> <chr>                  <LINESTRING [m]>    [m]
1     1     2        4 3434819 NA        NA    NA    NA      NA      NA      NA      3526003  Fort Erie   T        3526003 Fort Erie  T       35      Ontario 35      Ontario 5     23    (7269806 859183, 7269795 859166.~   25.9
2     3     4        5 1358252 South     LINE  NA    NA      NA      NA      NA      3551027  Gordon/Bar~ MU       3551027 Gordon/Ba~ MU      35      Ontario 35      Ontario 5     23    (6921247 1133452, 6921258 113345~ 1245. 
3     5     6        8 1778927 NA        NA    NA    NA      NA      NA      NA      3512054  Wollaston   TP       3512054 Wollaston  TP      35      Ontario 35      Ontario 5     23    (7320857 1089403, 7320903 108942~  254. 
4     7     8       11 731932  Albion    RD    NA    NA      NA      2010    2010    3520005  Toronto     C        3520005 Toronto    C       35      Ontario 35      Ontario 3     21    (7202555 935281.3, 7202653 93533~  111. 
5     9    10       18 3123617 County 41 RD    NA    640     708     635     709     3511015  Greater Na~ T        3511015 Greater N~ T       35      Ontario 35      Ontario 3     21    (7403627 1039328, 7403431 103963~  367. 
6    11    12       20 4814160 NA        NA    NA    NA      NA      NA      NA      3553005  Greater Su~ CV       3553005 Greater S~ CV      35      Ontario 35      Ontario 5     21    (7042806 1222708, 7042838 122273~  191.
... with 570,700 more rows

Node Data: 430,824 x 1
`_ogr_geometry_`
     &lt;POINT [m]&gt;

1   (7269806 859183)
2 (7269790 859162.7)
3  (6921247 1133452)
... with 430,821 more rows

From here, I am not sure how to use the tbl_graph to find the distance between the two points n1 and n2:

> n1
[1] 110393
> n2
[1] 319271

I think n1 and n2 correspond to the "to" and "from" columns in the "path" tbl_graph - but I am not sure how to use the tbl_graph to find out the distance between n1 and n2.

Can someone show me how to do this?

score 2 · Accepted Answer · answered Jun 20 '23 at 06:55

By default, as_sfnetwork creates a directed graph. This means all roads are one direction, defined (probably) by the order the vertices in the lines appear in the data.

To create an undirected graph, pass directed=FALSE.

I've cut down your road network to illustrate, so these numbers won't match the whole data set. I'm also working with the simplest igraph object before you mutate the weights etc.

net <- as_sfnetwork(a)

This is directed - all roads are one way - so no path between these two points:

net1  <- st_nearest_feature(pt1_transformed, net)
net2  <- st_nearest_feature(pt2_transformed, net)
sps = shortest_paths(net, net1, net2)
Warning message:
In shortest_paths(net, net1, net2) :
At core/paths/unweighted.c:368 : Couldn't reach some vertices.

Make it undirected:

net <- as_sfnetwork(a, directed=FALSE)

And now there's a way:

sps = shortest_paths(net, net1, net2)
sps
## $vpath
## $vpath[[1]]
## + > 
## 43/1040 vertices, from 8ad2ead:
##  [1]  547  501   74  402   99  107  108   82  350  164   33  177  229  113  114
## [etc]

If you have information about roads being one-way then you can probably add this to the graph structure at some point.

thank you so much for your answer! I am interested in learning how to find the distance of this shortest path - should I post a new question about this? Thank you so much! — stats_noob, Jun 20 '23 at 14:28
I asked it over here: https://gis.stackexchange.com/questions/462029/calculating-the-distance-of-a-path-on-a-road-network — stats_noob, Jun 20 '23 at 15:42

Working with tbl_graphs using R

Set the URL for the shapefile

Create a temporary folder to download and extract the shapefile

Download the shapefile to the temporary folder

Extract the shapefile from the downloaded zip file

Read the shapefile using the rgdal package

source: https://gis.stackexchange.com/questions/456748/strategies-for-working-with-large-shapefiles/456798#456798

Define the two points of interest in WGS84 (EPSG:4326)

Set the CRS of the points to WGS84 (EPSG:4326)

Transform the points to the CRS of the network

Find the nearest points on the network to the transformed points of interest

A tbl_graph: 430824 nodes and 570706 edges

A directed multigraph with 442 components

Edge Data: 570,706 x 25 (active)

... with 570,700 more rows

Node Data: 430,824 x 1

... with 430,821 more rows

1 Answers1

Warning message:

In shortest_paths(net, net1, net2) :

At core/paths/unweighted.c:368 : Couldn't reach some vertices.