R:计算两个地理点之间的距离

3

我对学习如何在R中处理道路网络文件很感兴趣。

例如,我想了解以下两个(加拿大)地址之间的驾车距离:

  • CN Tower: 290 Bremner Blvd, Toronto, ON M5V 3L9
  • Toronto Airport: 6301 Silver Dart Dr, Mississauga, ON L5P 1B2

过去,我会使用像OpenStreetMap(OSM)这样的API:

library(tmap)


library(tmaptools)
remotes::install_github("riatelab/osrm")

q1 = geocode_OSM("6301 Silver Dart Dr, Mississauga, ON L5P 1B2")
q2 = geocode_OSM("290 Bremner Blvd, Toronto, ON M5V 3L9")

q1 = as.numeric(q1$coords)
q2 = as.numeric(q2$coords)

q1_lat = q1[1]
q1_long = q1[2]
q2_lat = q2[1]
q2_long = q2[2]

route = osrmRoute(src = c(q1[1], q1[2]) ,  dst = c(q2[1], q2[2]), osrm.profile = "car")

> route$distance
[1] 26.2836

正如我们在这里所看到的,这两点之间的驾驶距离为26.2公里(与Google地图获取的距离相当接近)。

我的问题:我现在想尝试使用道路网络文件做类似的事情。

例如,我找到了以下包含有关道路网络信息的文件(https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/rnf-frr/index2021-eng.cfm?Year=21)。然后我将其以.shp格式(即shapefile)下载到了我的计算机上。

基于这样一个道路网络文件,是否可以找出任意两个地址(无论是语言地址还是地理坐标)之间的"行驶距离"?

谢谢!

注意: 这个文件似乎非常大,我不确定我的计算机能否完全加载它 - 是否可能命令计算机仅导入此文件的较小部分?(例如,导入省份=安大略省,导入城市=多伦多的部分)


2
你需要使用那个Statcan来源吗?或者基于另一个来源的答案,仅限于多伦多大都市地面交通网络,是否也可以? - Nicolás Velasquez PhD
@Nicolás Velásquez:非常感谢您的回复!我对了解更多关于这个主题和不同的解决问题的方法很感兴趣。我想也许Statcan来源可能更好,因为(理论上)它应该适用于加拿大任何地方...但我愿意学习不同的方法。非常感谢! - stats_noob
1个回答

5

是的,您可以下载地图并根据这些文件中的道路网络信息计算距离。您已经在使用osrm R包。这会向远程演示服务器发送请求,执行此操作。但是,docs文档中指出:

OSRM演示服务器不允许大型查询(超过10000个距离或持续时间)。

相反,您可以安装osrm-backend,这是一个用C++编写的高性能路由引擎。这将允许您基于提供的地图设置自己的路由服务器。然后,您可以从R中向本地服务器发出与上述相同的请求,而没有速率限制。

安装osrm-backend和构建地图

你可能是正确的,普通电脑很难轻松构建整个加拿大地图。如果设置了交换文件,可能可以实现,但可能需要很长时间(几小时到几天)。我使用过较小的地图,即由geofabrik托管的安大略省开放街道地图数据。您可以在这里下载其他地区。

最简单的方法是使用OSRM docker镜像。一旦安装了docker,以下内容将下载安大略省的数据和docker镜像,并启动服务器。以下内容应在终端中输入,而不是R。这应该适用于Windows(在PowerShell中),Mac或Linux。

# Download the map - may take about 15 mins - from https://download.geofabrik.de/north-america/canada/ontario.html
wget -O ontario-latest.osm.pbf https://download.geofabrik.de/north-america/canada/ontario-latest.osm.pbf

# May take 5-10 mins to download the Docker image the first time
docker run -t -v "${PWD}:/data" ghcr.io/project-osrm/osrm-backend osrm-extract -p /opt/car.lua /data/ontario-latest.osm.pbf || "osrm-extract failed"
docker run -t -v "${PWD}:/data" ghcr.io/project-osrm/osrm-backend osrm-partition /data/ontario-latest.osm.pbf || "osrm-partition failed"
docker run -t -v "${PWD}:/data" ghcr.io/project-osrm/osrm-backend osrm-customize /data/ontario-latest.osm.pbf || "osrm-customize failed"

# Run the image
docker run -t -i -p 5000:5000 -v "${PWD}:/data" ghcr.io/project-osrm/osrm-backend osrm-routed --algorithm mld /data/ontario-latest.osrm

这将在5000端口启动一个路由引擎HTTP服务器。您应该看到类似以下的输出:

[2023-03-27T18:04:02.703986704] [info] starting up engines, v5.27.1
[2023-03-27T18:04:02.704179704] [info] Threads: 8
[2023-03-27T18:04:02.704216504] [info] IP address: 0.0.0.0
[2023-03-27T18:04:02.704250504] [info] IP port: 5000
[2023-03-27T18:04:07.973464309] [info] http 1.1 compression handled by zlib version 1.2.11
[2023-03-27T18:04:07.973828509] [info] Listening on: 0.0.0.0:5000
[2023-03-27T18:04:07.973957909] [info] running and waiting for requests

然后您可以从 R 中查询此内容。

从 R 查询服务器

重要的是将您的 osrm.server 设置为本地主机:

library(osrm)
options("osrm.server" = "http://127.0.0.1:5000/")
options("osrm.profile" = "car") # Easiest to set this here as well

使用您问题中的示例坐标,我们可以进行如下操作:

osrmRoute(src = src, dst = dst)
# Simple feature collection with 1 feature and 4 fields
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: -79.6122 ymin: 43.61352 xmax: -79.38643 ymax: 43.68883
# Geodetic CRS:  WGS 84
#     src dst duration distance                       geometry
# 1_1   1   1 23.25667  26.2836 LINESTRING (-79.61214 43.68...

对于多个坐标,您还可以使用osrmTable()

osrmTable(src = src, dst = dst)
# $durations
#      1
# 1 23.3

# $sources
#         lon      lat
# 1 -79.61214 43.68332

# $destinations
#         lon      lat
# 1 -79.38643 43.64182

一旦您停止了osrm服务器Docker容器,您只需要运行最后一行即可重新启动它,即:

docker run -t -i -p 5000:5000 -v "${PWD}:/data" ghcr.io/project-osrm/osrm-backend osrm-routed --algorithm mld /data/ontario-latest.osrm

如果您有大量的坐标对,可能会遇到max-table-size参数的限制。您可以通过将其作为参数传递给docker run来增加它,例如:
docker run -t -i -p 5000:5000 -v "${PWD}:/data" ghcr.io/project-osrm/osrm-backend osrm-routed --algorithm mld --max-table-size 100000 /data/ontario-latest.osrm

你的问题有第二部分,关于将地址反向地理编码为经纬度。这是一个足够不同的问题,我不会在这里尝试回答它。然而,好消息是有一个docker镜像可以用于Nominatim,这就是tmaptools::geocode_OSM()在幕后查询的内容(同样受到速率限制)。你可以类似地安装它-如果遇到问题,请提出另一个问题。


@SamR:非常感谢您的回答!使用我提供的道路网络文件是否可以实现这个功能? - stats_noob
@stats_noob 是的 - 你只需要将 shapefile 转换为 .pbf 格式。我不知道如何在 R 中实现这一点,但使用 JOSM 似乎很简单。我从未这样做过,因为我发现 geofabrik OSM 文件非常好用。 - SamR
@ SamR:非常感谢您的回答!我会研究如何进行这种转换。我一直在学习和尝试自学地理空间分析(例如,使用shapefiles)-我花了一些时间尝试解决这个问题:https://dev59.com/5WEmtIgBPY-HTNNjslxj ....如果您有时间,能否请您看一下?非常感谢! - stats_noob

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接