Collecting Highway Bearings

Given a point, what highway is it on and what direction does that highway go?

I’m kicking off this blog with a geospatial example. I did some work recently where I needed to know:

• the nearest section of highway to a given point; and,
• the highways’s approximate bearing.

I got this working pretty well and am repeating the crux of the solution here. As an aside, while I was doing this work I came across Geoff Boeing’s beautiful post visualising similar aggregate directional information for 25 cities - definitely worth a look.

The best approach I could think of was to use OpenStreetMap data. Firstly associate a point with it’s nearest highway ‘way’ feature. Then get the coordinates of the beginning and end of the feature and use the pair of points to define a bearing. Makes sense? OK, here we go:

Get the packages

library(geosphere)
library(sf)
library(leaflet)
library(osmdata)
library(data.table)
library(ggmap)

Organise the data

To keep things manageable I’m going to use 4 points dotted around the Northern Quarter in Manchester (UK). I’ll put them in a data.frame:

df <- data.frame(ID = 1:4,
 Longitude = c(-2.232086, -2.235675, -2.233628, -2.232453),
 Latitude = c(53.484510, 53.484318, 53.484103, 53.481176), stringsAsFactors = FALSE)

Then I’ll convert to an sf (spatial features) multi point object. I use this command over and over.

df.sf <- st_as_sf(df, coords = c("Longitude", "Latitude"), crs=4326)

The crs argument is telling st_as_sf that these coordinates are WGS84 by referencing the appropriate EPSG code 4326 (it’s a good idea to memorise that number!). This is needed later when I project the points out of latitude/longitude format and into distance-from-bottom-left format. If you’re not familiar with this concept here’s a good pdf tutorial.

Visualise

At this point I’ll take a look at the 4 points. I’m using leaflet which works well with Blogdown.

m <- leaflet(df) %>% addTiles() %>% addCircleMarkers(
  lng = ~Longitude,
  lat = ~Latitude,
    radius = 8,
    color = "green",
    stroke = FALSE, fillOpacity = 0.9
  )
m

Get information from OpenStreetMap

The opq query from the osmdata package forms an OSM overpass filter and sends it to OSM’s servers. The result is an osmdata object containing various items. The item I’m interested in is the osm_lines which is a set of sf linestrings. These are the way objects representing highway segments. The query has returned 44 highway segments that intersected the bounded box of our 4 points (at the time of posting).

this.bbox <- st_bbox(df.sf)
osm.road.types = c("motorway", "motorway_link", "trunk", "trunk_link", "primary", "primary_link", "secondary", 
                   "secondary_link", "tertiary", "unclassified", "residential")
out <- opq (this.bbox) %>%
  add_osm_feature (key = "highway", value = osm.road.types) %>%
  osmdata_sf(quiet = F)
## Issuing query to Overpass API ...
## Rate limit: 2
## Query complete!
## converting OSM data to sf format
out
## Object of class 'osmdata' with:
##                  $bbox : 53.481176,-2.235675,53.48451,-2.232086
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 199 points
##             $osm_lines : 'sf' Simple Features Collection with 44 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 0 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL
road.lines <- out$osm_lines
nrow(road.lines)
## [1] 44

Get candidate closest roads to points

Before finding the closest highway segment to a point, I’ll narrow down the search by first identifying all the highway segments within a fixed distance of each point. This won’t make much difference here but does provide a big speed up when there’s thousands of highways to choose from. Also you can, if needed, put this code in a loop and grow or shrink the buffer distance until a suitable number of candidate highways are found. First step is to project the points and the road segments to a suitable coordinate reference system (for Manchester, the UK Ordnance Survay EPSG: 27700 is a good choice), apply a 50m buffer to the points, then identify which road segments intersect with this buffer.

df.sf.tf <- st_transform(df.sf, 27700) # OSGB 1936 British National Grid CRS
road.lines.tf <- st_transform(road.lines, 27700)
buffer.size <- 50
df.sf.tf.bf <- st_buffer(df.sf.tf, buffer.size)
intersected.items <- st_intersects(df.sf.tf.bf, road.lines.tf, sparse = F)
dim(intersected.items)
## [1]  4 44

The dimensions show I have 4 rows, one for each of the 4 points, and 44 columns, one for each OSM road segment within the bounding box. The matrix contains TRUE/FALSE values where TRUE indicates a point of intersection.

Get the closest road segment

So now, using the location of TRUE values in each row or the intersected.items matrix, I can subset the road.lines so that the nearest highway search is constrained to a much smaller number of highway segments (ones within 50m). The st_coordinates function is useful to return just the coordinate information from an sf object, and the dist2Line function from the geosphere package calculates the shortest distance between a point and a polyline (There are some good details about doing this type of calculation at the Movable Type website.

number.of.point <- nrow(intersected.items)
distances.vect <- matrix(NA_real_, nrow = nrow(df.sf.tf), ncol = 4)
colnames(distances.vect) <- c("distance", "lon", "lat", "ID")
for (i in seq_len(number.of.point)) {
  # Out of all 44 road segments, get the index of the ones that are within 50m of point i. 
  original.road.idx <- which(intersected.items[i, ] == TRUE)
  distances.vect.temp  <- dist2Line(p = st_coordinates(df.sf[i, ]), line = as_Spatial(road.lines[original.road.idx, ]))
  distances.vect.temp[1, "ID"] <- original.road.idx[distances.vect.temp[1, "ID"]]
  distances.vect[i, ] <- distances.vect.temp
}
distances.vect
##      distance       lon      lat ID
## [1,] 5.061288 -2.232137 53.48448 43
## [2,] 3.500525 -2.235712 53.48430 10
## [3,] 6.111912 -2.233558 53.48407  5
## [4,] 5.541079 -2.232491 53.48122 19

So now I’ve got the nearest road segments and the coordinates of the nearest point on the road segments to the original points. These are stored in the lon/lat columns of distances.vec. So the original vague points are now snapped to actual highway locations which is a bonus.

Get the road section bearings

Next I’ll extract the chosen 4 road segments from the full set of 44 and add the useful features I’ve just created so far. Getting the coordinates of the first and last points is a bit fiddly. I’m using data.table for this. When retrieving the coordinates via st_coordinates a third column labelled L1 is created with a group ID of each polyline object. This can be used as a grouping variable in data.table. The .SD[...] syntax holds the sub-data.table once the full data.table has been split by the grouping variable, finally 1 and .N are the first row and last row of each sub-data.table. There’s other ways to do this but the speed of data.table would pay off if I had thousands of points. Once we have the beginning and end coordinate pairs it’s straightforward to call the bearing function from the geosphere package.

out <- road.lines[as.integer(distances.vect[, "ID"]), ]
out$dist.to.nearest <- distances.vect[, "distance"]
out$snap.lat <- distances.vect[, "lat"]
out$snap.lon <- distances.vect[, "lon"]
matching.road.lines.coords <- data.table::as.data.table(st_coordinates(out))
matching.road.lines.coords.begin <- as.matrix(matching.road.lines.coords[, .SD[1, .(X, Y)], by = "L1"][, L1 := NULL])
matching.road.lines.coords.end <- as.matrix(matching.road.lines.coords[, .SD[.N, .(X, Y)], by = "L1"][, L1 := NULL])
matching.road.lines.bearing <- bearing(matching.road.lines.coords.begin, matching.road.lines.coords.end)
out$bearing <- (matching.road.lines.bearing + 360) %% 360
out$bearing
## [1] 133.2409 316.6075 217.5574 258.1974

Plot with estimated bearing

In order to visualise the estimated road bearing for my 4 points I need to add arrow endpoints information. Some trigonometry wil accomplish this (the 0.00025 value is determines the length of the direction arrow):

out$delta.x <- sin((pi*out$bearing) / 180) * 0.00025
out$delta.y <- cos((pi*out$bearing) / 180) * 0.00025

Now I can plot the snapped locations and show an arrow indicating the highway bearing over a basemap which is not only pretty but allows for a quick visual validation.

# A handcoded bounding box for tile retrieval. A function to make this would be a nice addition.
n.quarter.bb <- c(left = -2.236, bottom = 53.48, right = -2.23, top = 53.485)
n.quarter <- get_stamenmap(bbox = n.quarter.bb, zoom = 18, messaging = FALSE)
ggmap(n.quarter) + 
  geom_point(data = out, aes(x=snap.lon, y=snap.lat), colour = "darkgreen", size = 4) +
  geom_segment(data = out, aes(x=snap.lon, y=snap.lat, xend = snap.lon + delta.x, yend = snap.lat + delta.y), colour = "darkgreen", arrow = arrow(length = unit(0.2,"cm")))

Summary

The results are pretty good, definitely good enough to classify a point into one of the 8 compass bearing directions: N, NE, E… etc. In this post I’ve explored working with sf spatial points, generated and run an OSM overpass filter, and done some useful visualisation with Leaflet and ggmap. I’m sure there are things I could tidy and optimise here but I hope you found this a little useful.