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..
I then converted the shapefile into a graph network, defined the coordinates and changed the coordinate systems of both points:
# Convert the shapefile to an sfnetwork object
net <- as_sfnetwork(a, directed = FALSE)
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)
Finally, I then identified the shortest path between these two points:
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]
My Question: Now, I am trying to calculate the distance of this path.
I know in general how to calculate the distance of a network path in R:
# example for path distance calculation
set.seed(123)
from <- sample(paste0("node_", 1:10), 20, replace = TRUE)
to <- sample(paste0("node_", 1:10), 20, replace = TRUE)
df <- data.frame(from, to)
df <- unique(df)
df$weight <- runif(nrow(df))
g <- graph_from_data_frame(df)
#plot(g, edge.width = df$weight * 5)
distance of the path : node_3 -> node_9 -> node_5
d1 <- distances(g,"node_3","node_9")
d2 <- distances(g,"node_9","node_5")
total_distance <- as.numeric(d1) + as.numeric(d2)
#0.3845944
But I am not sure how to perform to calculate for the shortest path for sps$vpath path.
Can someone please show me how to do this?
