In a recent #rspatial Twitter conversation I discussed analytical problem of finding establishments of a certain type along a walking route from point A to point B.

It is an interesting problem, one that deserves a bit more space than 280 characters allow.

The issue can be split into two analytical problems:

  • finding the walking route from point A to point B (as a line)
  • finding points of interest of a given type within distance from this line

Both analytical use cases are lend themselves easily to techniques from spatial R toolset.

In this demonstration I will be using two packages / APIs:

This blog post will focus more on technique used rather then the actual result; however a toy problem is necessary to demonstrate the tools. I will take inspiration from the “Tipperary” song and look for a place to have a drink along the route between two London landmarks - Piccadily Circus and Leicester square (it’s not a long walk by the way).

The first step is loading the packages and setting up the API key via hereR::set_key().

library(dplyr)   # ... because dplyr :)
library(osmdata) # data part
library(hereR)   # routing part / needs API
library(sf)      # spatial data processing

# I keep my key as an environment variable to prevent accidental sharing
secret_key <- Sys.getenv("HERE_API_KEY") # insert your API key here...

set_key(secret_key)

The second step is geocoding the locations and calculating the walking route (transport mode = pedestrian) between the two.

It will be returned in format of {sf} package, with vertical dimension (height) reflected. As the Z dimension can lead to issues in plotting and buffering it is best to remove it immediately.

# It's a long way, to Tipperary...
destinations <- c("Piccadilly Circus, London",
                  "Leicester Square, London") %>% 
  hereR::geocode() # geocode

# routing call - shortest walking path between the first and second destination
walking_route <- hereR::route(origin = destinations[1, ],
                      destination = destinations[2, ],
                      transport_mode = "pedestrian") %>% 
  st_zm() # vertical dimesion is neither required, nor helpful

walking_route
## Simple feature collection with 1 feature and 11 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -0.13451 ymin: 51.51008 xmax: -0.13022 ymax: 51.51057
## geographic CRS: WGS 84
##   id rank section           departure             arrival       type       mode
## 1  1    1       1 2021-02-22 16:37:09 2021-02-22 16:42:29 pedestrian pedestrian
##   distance duration duration_base consumption                       geometry
## 1      318      320           320          NA LINESTRING (-0.13451 51.510...

Third step is creating a buffer along the line object; a the distance will be in meters I am briefly reprojecting my line from WGS84 (in degrees) to Ordnance Survey CRS (any metric CRS would do).

As the OSM results come in EPSG:4326 and degrees are required for leaflet visualization I am reverting back to the safety of WGS84 immediately.

# a buffer of 20 meters from route; a transformation to metric 
# e.g. Ordnance Survey EPSG is required to use meters 
# (units of WGS84 are degrees, making a buffer somewhat complicated)
buffer <- walking_route %>% 
  # this could be another planar CRS relevant to the area of interest
  st_transform("EPSG:27700") %>% 
  st_buffer(20) %>% # size of the buffer / in units of the planar CRS
  st_transform(4326) # back to the safety of WGS84!

The fourth step involves creating an overpass query to OSM backend. This takes two arguments

  • area of interest, defined by the bounding box of my buffer object
  • type of OSM features, by key & value as listed on the OSM wiki.

Once the search result is obtained it is spatially joined to the buffer object via inner join technique (not left means inner join).

This will leave only those amenities that are within the buffer polygon.

# search result  
search_res <- opq(bbox = st_bbox(buffer), nodes_only = T) %>% 
  # 1) area of interest = the buffer
  add_osm_feature(key = "amenity",  
                  value = c("bar", "restaurant", "pub")) %>%
  # 2) search values - a watering hole, in OSM parlance
  osmdata_sf(quiet = F) # 3) result format

# osmdata returns multiple types; only points are relevant in our search
watering_holes <- search_res$osm_points 

# and finally - the intersection of buffer & search result
result <- watering_holes %>% 
  st_join(buffer, left = F) %>%  # not left = filtering join!
  select(osm_id, name)

result
## Simple feature collection with 4 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -0.1329428 ymin: 51.51005 xmax: -0.1324365 ymax: 51.51041
## geographic CRS: WGS 84
##                osm_id                 name                    geometry
## 348904266   348904266     Angus Steakhouse POINT (-0.1329228 51.51005)
## 348905523   348905523                HAPPY POINT (-0.1324716 51.51013)
## 4376186998 4376186998 Bubba Gump Shrimp Co POINT (-0.1324365 51.51041)
## 6255364830 6255364830     HaiDiLao Hot Pot POINT (-0.1329428 51.51033)

My work could in principle stop here, but in order to make certain we can perform a visual check.

I am also showing off a little by using the font awesome tankard icon, just because.

# these watering holes did not make the cut...
holes_oos <- watering_holes %>% # oos = out of scope
  anti_join(st_drop_geometry(result), by = "osm_id") # in scope are in scope

# a visual check / 

library(leaflet)

# a beer tankard icon
awesome <- makeAwesomeIcon(
  icon = "beer",
  iconColor = "orange",
  markerColor = "green",
  library = "fa"
)

leaflet(result) %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addAwesomeMarkers(icon = awesome,
                    label = ~ name) %>% 
  addPolygons(data = buffer,
              stroke = F,
              fillColor = "blue",
              fillOpacity = 1/5,
              label = "20 meter neigbourhood") %>% 
  addPolylines(data = walking_route,
               color = "red",
               label = "walking route") %>% 
  addCircleMarkers(data = holes_oos,
                   stroke = F,
                   fillColor = "gray",
                   fillOpacity = 1/2,
                   label = ~ name)


I believe my walkthrough demonstrated that both routing and finding points of interest are tasks that can be handled without leaving the comfort of your R session.

More advanced topics in the same vein include finding points of interest within walking distance from a coordinate (using hereR::isoline() approach) and querying historical OSM data (by utilizing the datetime and datetime2 arguments of your osmdata::opq() call).