英文:
Calculating Distances Between Points on a Shapefile
问题
我正在使用R编程语言。
我对如何计算两组坐标之间的驾驶距离（例如，基于道路网络）感兴趣。
例如：
- CN塔：290 Bremner Blvd，Toronto，ON M5V 3L9（-79.61203，43.68312）
- 多伦多机场：6301 Silver Dart Dr，Mississauga，ON L5P 1B2（-79.61203，43.64256）
为了解决这个问题，我尝试下载加拿大道路网络的形状文件并将其子集化为安大略省：
library(sf)
library(rgdal)
library(sfnetworks)
# 设置形状文件的URL
url <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"
# 创建一个临时文件夹来下载和提取形状文件
temp_dir <- tempdir()
temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")
# 下载形状文件到临时文件夹
download.file(url, temp_file)
# 从下载的zip文件中提取形状文件
unzip(temp_file, exdir = temp_dir)
# 使用rgdal包读取形状文件
a = st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), query="select * from lrnf000r22a_e where PRUID_R ='35'")
导入R后，文件看起来像这样：
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> <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> <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> <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> <NA> 3537028 Amherstburg T 3537028 Amherstburg T 35 Ontario 35 Ontario 5
8 24 4825761 <NA> <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..
然后，通过参考不同的参考资料（例如 https://cran.r-project.org/web/packages/sfnetworks/vignettes/sfn01_structure.html，https://stackoverflow.com/questions/66513615/creating-a-route-from-a-list-of-points-and-overlaying-the-route-list-of-points） - 我尝试计算这两个点之间的距离：
# 将形状文件转换为sfnetwork对象
net <- as_sfnetwork(a)
# 定义起点和终点
q1 <- st_point(c(-79.61203, 43.68312))
q2 <- st_point(c(-79.38709, 43.64256))
# 将点的CRS设置为与形状文件的CRS匹配
q1 <- st_sfc(q1, crs = st_crs(a))
q2 <- st_sfc(q2, crs = st_crs(a))
# 找到两点之间的最短路径
path <- st_network_paths(net, q1, q2)
#
<details>
<summary>英文:</summary>
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, https://stackoverflow.com/questions/66513615/creating-a-route-from-a-list-of-points-and-overlaying-the-route-list-of-points) - 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!
</details>
# 答案1
**得分**: 1
你很接近得到期望的结果。错误发生是因为`st_network_paths`函数期望一个`sfnetwork`对象，但实际上它收到了一个`data.frame`。问题似乎出在将你的 "a" 对象转换为`sfnetwork`对象时使用`as_sfnetwork()`函数。
请尝试以下代码以修复这个问题：
```R
library(sf)
library(rgdal)
library(sfnetworks)
library(units)
# 你之前的代码一直保持不变，直到这一点
# 将 shapefile 转换为 sfnetwork 对象
net <- as_sfnetwork(a, directed = FALSE) # 如果你的道路网络是无向的，设置 directed = FALSE
# 将点投影到与 shapefile 的 CRS 匹配
q1 <- st_transform(q1, st_crs(a))
q2 <- st_transform(q2, st_crs(a))
# 找到每个点在网络上的最近节点
from_node <- st_nearest_feature(q1, net)
to_node <- st_nearest_feature(q2, net)
# 找到两点之间的最短路径
path <- st_network_paths(net, from_node, to_node)
# 计算路径的距离（以米为单位）
distance <- sum(path$length) %>% set_units("meters")
这应该修复错误，并为你提供道路网络上两个点之间的正确驾驶距离。
编辑
这是修订后的代码：
library(sf)
library(rgdal)
library(sfnetworks)
library(units)
# 你之前的代码一直保持不变，直到你将 shapefile 读入 "a" 中
# 将数据转换为 WGS 84
a <- st_transform(a, 4326)
# 将 shapefile 转换为 sfnetwork 对象
net <- as_sfnetwork(a, directed = FALSE)
# 创建具有与 shapefile 相同 CRS 的点
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)
# 找到每个点在网络上的最近节点
from_node <- st_nearest_feature(q1, net)
to_node <- st_nearest_feature(q2, net)
# 找到两点之间的最短路径
path <- st_network_paths(net, from_node, to_node)
# 计算路径的距离（以米为单位）
distance <- sum(path$paths[[1]]$length) %>% set_units("meters")
这应该为你提供两个点之间的正确距离。请注意，实际值将取决于输入道路网络数据的准确性和完整性。
英文:
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.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论