R: Calculating the Distance Between Two Geographical Points

Question

I am interested in learning how to work with Road Network Files in R.

For example, I would be interested in finding out the driving distance between the following two (Canadian) addresses:

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

In the past, I would have used an API such as the OpenStreetMap (OSM):

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

As we can see here, the driving distance between these two points is 26.2 KM (which is quite close to distance obtained from Google Maps)

My Question: I would now like to try and do something similar using a Road Network File.

For example, I found the following file which contains information about the Road Networks (https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/rnf-frr/index2021-eng.cfm?Year=21). I then downloaded this to my computer in .shp format (i.e. shapefile).

Based on such a file of Road Networks, is it possible to find out the "driving distance" between any two addresses (whether in "language address" or geographical coordinates)?

Thanks!

Note: This file appears to be quite large and I am not sure if my computer can fully load it - is it possible to command the computer to only import a smaller portion of this file? (e.g. import where province = ontario , import where city = toronto)

Answer 1

Yes, you can download maps and calculate distances based on the information about road networks in these files. You are already using the osrm R package. This sends requests to a remote, demo server, which does exactly that. However, the package docs state,

The OSRM demo server does not allow large queries (more than 10000 distances or durations).

Instead, you can install osrm-backend, a high performance routing engine written in C++. This will allow you to set up your own routing server, based on maps you provide. You can then make the same requests as above from within R to your local server, with no rate limits.

Installing osrm-backend and building the map

You are probably correct that you will not be able to easily build the entire Canada map on a standard PC. If you set up your swap file you may be able to but it could take a long time (hours to days). I have used a smaller map, the Ontario Open Street Map data hosted by geofabrik. You can download other regions here.

The easiest way is to use the OSRM docker image. Once docker is installed, the following will download the Ontario data and docker image, and start the server. The following is to be entered in a terminal, rather than R. This should work on Windows (in PowerShell), Mac or 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

This start a routing engine HTTP server on port 5000. You should see output similar to:

[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

You can then query this from R.

Querying the server from R

The important thing is to set your osrm.server to the local host:

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

Using the example coordinates from your question, we can do:

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...

For multiple coordinates, you can also use 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

Once you've stopped the osrm server Docker container, you only need the last line to run it again, i.e.

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

If you have a large number of coordinate pairs, you may hit the max-table-size parameter. You can increase this by passing it as an argument to docker run, e.g.:

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

There was a second part of your question about reverse geocoding addresses to lat/lon. It is a sufficiently different question that I'm not going to try to answer it here. However, the good news is there is a docker image for Nominatim, which is what tmaptools::geocode_OSM() is querying under the hood (again with rate limits). You can install that in a similar way - ask another question if you have trouble.