英文:
How to Use "tbl_graphs" in R?
问题
我正在使用R编程语言工作。
基于一个shapefile，我试图可视化/找出两个坐标之间的驾驶距离（例如，CN塔和多伦多皮尔逊机场）。
首先我加载了shapefile：
library(sf)library(rgdal)library(sfnetworks)library(igraph)library(dplyr)library(tidygraph)# 设置shapefile的URLurl <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"# 创建一个临时文件夹来下载和解压缩shapefiletemp_dir <- tempdir()temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")# 下载shapefile到临时文件夹download.file(url, temp_file)# 从下载的zip文件中解压shapefileunzip(temp_file, exdir = temp_dir)# 使用rgdal包读取shapefile# 来源：https://gis.stackexchange.com/questions/456748/strategies-for-working-with-large-shapefiles/456798#456798a = st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), query="select * from lrnf000r22a_e where PRUID_R ='35'")
shapefile看起来像这样：
包含570,706个要素和21个字段的简单要素集几何类型：LINESTRING维度：XY边界框：xmin：5963148 ymin：665490.8 xmax：7581671 ymax：2212179投影的CRS：NAD83 / 加拿大统计局兰伯特前10个要素：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 RANK1 4 3434819 <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 3526003 Fort Erie T 3526003 Fort Erie T 35 Ontario 35 Ontario 52 5 1358252 South LINE <NA> <NA> <NA> <NA> <NA> 3551027 Gordon/Barrie Island MU 3551027 Gordon/Barrie Island MU 35 Ontario 35 Ontario 53 8 1778927 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 3512054 Wollaston TP 3512054 Wollaston TP 35 Ontario 35 Ontario 54 11 731932 Albion RD <NA> <NA> <NA> 2010 2010 3520005 Toronto C 3520005 Toronto C 35 Ontario 35 Ontario 35 18 3123617 County 41 RD <NA> 640 708 635 709 3511015 Greater Napanee T 3511015 Greater Napanee T 35 Ontario 35 Ontario 36 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 57 21 1817031 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 3537028 Amherstburg T 3537028 Amherstburg T 35 Ontario 35 Ontario 58 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 59 25 544891 Dunelm DR <NA> 1 9 2 10 3526053 St. Catharines CY 3526053 St. Catharines CY 35 Ontario 35 Ontario 510 28 1835384 Seven Oaks DR <NA> 730 974 731 975 3515005 Otonabee-South Monaghan TP 3515005 Otonabee-South Monaghan TP 35 Ontario 35 Ontario 5CLASS _ogr_geometry_1 23 LINESTRING (7269806 859183,...2 23 LINESTRING (6921247 1133452...3 23 LINESTRING (7320857 1089403..
然后我尝试通过将道路网络视为图来计算距离：
# 将shapefile转换为sfnetwork对象net <- as_sfnetwork(a)# 在WGS84（EPSG:4326）中定义两个兴趣点pt1 <- st_point(c(-79.61203, 43.68312))pt2 <- st_point(c(-79.61203, 43.64256))# 将点的CRS设置为WGS84（EPSG:4326）pt1 <- st_sfc(pt1, crs = 4326)pt2 <- st_sfc(pt2, crs = 4326)# 将点转换为网络的CRSpt1_transformed <- st_transform(pt1, st_crs(net))pt2_transformed <- st_transform(pt2, st_crs(net))# 找到变换后的兴趣点在网络上的最近点n1 <- st_nearest_feature(pt1_transformed, net)n2 <- st_nearest_feature(pt2_transformed, net)
从这里，我创建了一个“tbl_graph”：
path <- net %>%activate(edges) %>%mutate(weight = edge_length()) %>%as_tbl_graph()> path<details><summary>英文:</summary>I am working with the R programming language.**Based on a shapefile, I am trying to visualize/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 shapefileurl <- "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 shapefiletemp_dir <- tempdir()temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")# Download the shapefile to the temporary folderdownload.file(url, temp_file)# Extract the shapefile from the downloaded zip fileunzip(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#456798a = 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 fieldsGeometry type: LINESTRINGDimension: XYBounding box: xmin: 5963148 ymin: 665490.8 xmax: 7581671 ymax: 2212179Projected CRS: NAD83 / Statistics Canada LambertFirst 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 RANK1 4 3434819 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 3526003 Fort Erie T 3526003 Fort Erie T 35 Ontario 35 Ontario 52 5 1358252 South LINE <NA> <NA> <NA> <NA> <NA> 3551027 Gordon/Barrie Island MU 3551027 Gordon/Barrie Island MU 35 Ontario 35 Ontario 53 8 1778927 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 3512054 Wollaston TP 3512054 Wollaston TP 35 Ontario 35 Ontario 54 11 731932 Albion RD <NA> <NA> <NA> 2010 2010 3520005 Toronto C 3520005 Toronto C 35 Ontario 35 Ontario 35 18 3123617 County 41 RD <NA> 640 708 635 709 3511015 Greater Napanee T 3511015 Greater Napanee T 35 Ontario 35 Ontario 36 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 57 21 1817031 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 3537028 Amherstburg T 3537028 Amherstburg T 35 Ontario 35 Ontario 58 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 59 25 544891 Dunelm DR <NA> 1 9 2 10 3526053 St. Catharines CY 3526053 St. Catharines CY 35 Ontario 35 Ontario 510 28 1835384 Seven Oaks DR <NA> 730 974 731 975 3515005 Otonabee-South Monaghan TP 3515005 Otonabee-South Monaghan TP 35 Ontario 35 Ontario 5CLASS _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 objectnet <- 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 networkpt1_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 interestn1 <- 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.92 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_`<POINT [m]>1 (7269806 859183)2 (7269790 859162.7)3 (6921247 1133452)# ... with 430,821 more rows**My Question:** 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] 319271I 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 please show me how to do this?Thanks!</details># 答案1**得分**: 4`n1` 和 `n2` 是一个名为 `net` 的 `tibble` 的行号索引。这些值可以用于子集化图列表，或者作为 `igraph::shortest_paths()` 函数的 `from` 和 `to` 参数。在这种情况下，如果你使用 `n1` 和 `n2` 来查找最短路径，你可能会看到一些相当奇怪的结果，或者可能只是一个空列表。此外，`sfnetworks` 对象 `net` 继承自 `tbl_graph` 和 `igraph`，因此 `tiynetworks` 和 `igraph` 的方法应该在不需要转换的情况下正常工作。将 `net` 转换为 `path`，并添加权重，不会更接近所需的路径。它仍然是几乎与之前相同的 `tbl_graph`，只是没有了空间图层。`sfnetworks` 还实现了满足你的主要目标所需的所有内容：`st_network_paths()` 和/或 `to_spatial_shortest_paths()` 转换器。`sfnetwork` 也被配置为有向（这是默认设置），这很可能会破坏所有的路由尝试，至少在这种情况下。`st_network_paths()` 方法示例：``` rlibrary(sf)library(sfnetworks)library(mapview)library(dplyr)library(ggplot2)# ...（代码的其余部分，已经在上面提供了）
igraph 方法示例：
library(igraph)# ...（代码的其余部分，已经在上面提供了）
tidygraph 和 sfnetworks 的 morphers 示例：
library(tidygraph)# ...（代码的其余部分，已经在上面提供了）
请注意，以上示例代码中包含了其他内容，只提供了与翻译请求相关的部分。
英文:
n1 & n2 are indeces, row numbers of one of the net tibbles. Those values can be used to subset graph lists or used as from and to parameters for igraph::shortest_paths(). In this case, n1 & n2 probably targeted edges tibble as nodes were not activated first. Means that if you used those n1 & n2 for finding shortest path, you probably witnessed some rather funky results. Or perhaps just an empty list.
Also, the sfnetworks object, net, inherits from tbl_graph and also igraph, so tiynetworks and igraph methods should work just fine on netwithout converting it first. net -> path conversion, while adding weights, leads no closer to desired route. It's still (almost) the same tbl_graph it was before, now without a spatial layer. sfnetworks also implements all what's needed for your main objective: st_network_paths() and/or to_spatial_shortest_paths() morpher.
sfnetwork was also configured to be directed (well, it's a default setting), this most likely will ruin all routing attempts, at least in this case.
sfnetworks::st_network_paths() approach
library(sf)library(sfnetworks)library(mapview)library(dplyr)library(ggplot2)shp_url <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"if (!file.exists(basename(shp_url))) curl::curl_download(shp_url, basename(shp_url), quiet = FALSE)# keep archive zipped;# predefined bbox, small patch in Toronto,# covers cover both points and leaves plenty of routing optionsa <- read_sf(paste0("/vsizip/", basename(shp_url)),wkt_filter = "POLYGON ((7200615 918857, 7211499 918857, 7211499 933387, 7200615 933387, 7200615 918857))")a#> Simple feature collection with 7584 features and 21 fields#> Geometry type: LINESTRING#> Dimension: XY#> Bounding box: xmin: 7199959 ymin: 918285.7 xmax: 7211930 ymax: 934030.6#> Projected CRS: NAD83 / Statistics Canada Lambert#> # A tibble: 7,584 × 22#> OBJECTID NGD_UID NAME TYPE DIR AFL_VAL ATL_VAL AFR_VAL ATR_VAL CSDUID_L#> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>#> 1 668 4668376 Dixie RD <NA> 2520 2520 2531 2543 3521005#> 2 874 4617171 401 HWY <NA> <NA> <NA> <NA> <NA> 3521005#> 3 2252 4092570 Breton AVE <NA> 178 300 175 295 3521005#> 4 3320 485941 Kipling AVE <NA> <NA> <NA> <NA> <NA> 3520005#> 5 3698 4775771 Eringa… DR <NA> 99 105 108 120 3520005#> 6 3700 507598 Bristol RD W <NA> <NA> <NA> <NA> 3521005#> 7 4472 4683515 Airport RD <NA> <NA> <NA> <NA> <NA> 3521005#> 8 4759 5779033 401 HWY <NA> <NA> <NA> <NA> <NA> 3520005#> 9 5051 505570 Grand … RD <NA> 3496 3530 3511 3525 3521005#> 10 5533 5876974 Yorkda… CRES <NA> 31 67 30 58 3520005#> # ℹ 7,574 more rows#> # ℹ 12 more variables: CSDNAME_L <chr>, CSDTYPE_L <chr>, CSDUID_R <chr>,#> # CSDNAME_R <chr>, CSDTYPE_R <chr>, PRUID_L <chr>, PRNAME_L <chr>,#> # PRUID_R <chr>, PRNAME_R <chr>, RANK <chr>, CLASS <chr>,#> # geometry <LINESTRING [m]># from / to points for routing:pts <- tribble(~point, ~lon, ~lat,"pt1", -79.61203, 43.68312,"pt2", -79.61203, 43.64256) %>%st_as_sf(coords =c("lon", "lat"), crs = st_crs(4326)) %>%st_transform(st_crs(a))# create un-directed network (default is directed),# store edge lenghts as weightsnet <- as_sfnetwork(a, directed = FALSE) %>%activate("edges") %>%mutate(weight = edge_length())# get paths between pt1 & pt2 with sfnetworkpaths_sfn <- st_network_paths(net,from = pts[pts$point == "pt1", ],to = pts[pts$point == "pt2", ],weights = "weight")paths_sfn#> # A tibble: 1 × 2#> node_paths edge_paths#> <list> <list>#> 1 <int [54]> <int [53]># includes list of edge indeces that form the path:paths_sfn$edge_paths[[1]]#> [1] 1027 1032 1295 697 121 2704 2702 2778 1125 1280 4364 1360 235 1686 2048#> [16] 6115 2266 865 2508 7511 7179 2956 1823 5139 1055 5929 3684 3417 5960 1196#> [31] 3383 1603 1975 1159 6250 5814 58 1964 3593 4423 3477 4940 383 5784 2362#> [46] 4120 7251 1970 4834 5155 624 1935 5134# get edges with matching indeces:route <- net %>%activate("edges") %>%slice(paths_sfn$edge_paths[[1]]) %>%st_as_sf()# distance:sum(route$weight)#> 7505.43 [m]# visualizemapview(list(pts, route))
<!-- -->
Path-finding and distance with igraph, using the same net object
library(igraph)# to find node indeces, keep net nodes activatednet <- net %>% activate("nodes")n1 <- st_nearest_feature(pts[1,], net)n2 <- st_nearest_feature(pts[2,], net)node1; node2#> [1] 1560#> [1] 3938# locate those on a map:net %>%st_as_sf() %>%slice(c(n1, n2)) %>%mapview(layer.name = "net nodes") +mapview(pts)
<!-- -->
# as sfnetwork inherits from tbl_graph and igraph:class(net)#> [1] "sfnetwork" "tbl_graph" "igraph"# we can use igraph::shortest_paths() with our sfnetwork object:paths_ig <- shortest_paths(net, n1, n2, output = "epath")# paths_ig$epath is a list of edges:str(paths_ig)#> List of 4#> $ vpath : NULL#> $ epath :List of 1#> ..$ : 'igraph.es' int [1:53] 1027 1032 1295 697 121 2704 2702 2778 1125 1280 ...#> .. ..- attr(*, "env")=<weakref>#> .. ..- attr(*, "graph")= chr "aac7eae8-1b4f-11ee-a8b6-53eec7eb8cd4"#> $ predecessors : NULL#> $ inbound_edges: NULL# when converted to integer, we'll have the same index list as# from sfnetwork::st_network_paths() :e <- paths_ig$epath[[1]] %>% as.integer()# check if those are indeed identical:identical(e, paths_sfn$edge_paths[[1]])#> [1] TRUE# sum edge wights to get path distance:sum(E(net)[e]$weight)#> 7505.43 [m]# we have not altered our original sf object (a) nor# sfnetwork object (net), so indeces should still match and we can use# e for subsetting: a[e,]ggplot() +geom_sf(data = a, color = "grey80") +geom_sf(data = a[e,], linewidth = 1.5, alpha = .6) +geom_sf(data = pts, aes(color = point), size = 3) +coord_sf(datum = st_crs(a), ylim = c(922000, 930000)) +theme_light()
<!-- -->
Morphers by tidygraph and sfnetworks
There are also morphers, i.e. tidygraph::to_shortest_path (would work with previous n1 and n2 objects for from/ to) and sfnetworks::to_spatial_shortest_paths (takes sf points for from / to, just like st_network_paths()) , so there's no need to worry about picking relevant nodes / edges from the network, convert() output is already filtered:
library(tidygraph)converted_net <- net %>% convert(to_spatial_shortest_paths, pts[1,], pts[2,])route_edges <- st_as_sf(converted_net, "edges")route_nodes <- st_as_sf(converted_net, "nodes")# route length:sum(route_edges$weight)#> 7505.43 [m]ggplot() +geom_sf(data = a, color = "grey80") +geom_sf(data = route_edges) +geom_sf(data = route_nodes, alpha = .2) +coord_sf(datum = st_crs(a), ylim = c(922000, 930000)) +theme_light()
<!-- -->
<sup>Created on 2023-07-05 with reprex v2.0.2</sup>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论