Packages

R is very capable as a code-led tool for geospatial manipulation and mapping. Along with {tidyverse} for data wrangling, there’s a few geospatial packages we need: {geojsonio} lets us read GeoJSON files, {sf} is for handling ‘special features’ geometry in a ‘tidy’ way, and {leaflet} lets us create interactive maps. {osmextract} was the main motivation for this post; it fetches OpenStreetMap features pretty painlessly.

library(geojsonio)
library(leaflet)
library(leaflet.extras)
library(osmextract)
library(sf)
library(tidyverse)

All these packages can all be downloaded from CRAN with install.packages().

The boundary

For purposes of this post, I wanted to look at a small, contained, ‘regularly-shaped’ geographic area to keep things simple. It didn’t have to be a country.

I settled on the Isle of Wight (IOW)¹, a small island off the south-coast of England. It’s mostly rural, with farms, hedges and waterways to cross, but there are certainly more built-up areas. It also helps that the IOW is featured in the {osmextract} documentation!

The Isle of Wight: a vexillologist’s delight.

First thing is to fetch a polygon that represents the extent of the island. Fortunately, Local Authority District (LAD) boundaries for the UK are available to download from the Official for National Statistics (ONS) in GeoJSON form. We can filter for the IOW, specifically.

# LAD boundaries (December 2020) UK BGC from the ONS Open Geography Portal
# Generalised (20m) - clipped to the coastline (Mean High Water mark)
# https://geoportal.statistics.gov.uk/datasets/local-authority-districts-december-2020-uk-bgc
geojson_path <- 
  "https://opendata.arcgis.com/datasets/db23041df155451b9a703494854c18c4_0.geojson"

# Read LAD boundaries, filter to IOW
iow_extent <- geojson_read(geojson_path, what = "sp") %>% 
  st_as_sf(crs = 4326) %>%
  filter(LAD20NM == "Isle of Wight")

I’ve used boundaries that are ‘clipped to the coastline’ because I don’t think you should have to swim out to the low-water mark to complete such a challenge.

OpenStreetMap features

We want to identify features like hedgerows, buildings and waterways that will become obstructions for our imaginary walk across the island. The oe_get() function from {osmextra} is an easy way to pull features from OpenStreetMap en masse. You can supply a location and receive features within that area.

First, the polygonal features, which you can get with argument layer = "multipolygons". You can see that a geometry column is returned, which contains the coordinates for the polygons.

# Fetch polygonal features for IoW
osm_polys <- oe_get(
  "Isle of Wight", # geographic area of interest
  layer = "multipolygons",  # fetch polygons
  stringsAsFactors = FALSE,  # return character-class
  quiet = TRUE  # don't print info
) %>%
  st_transform(crs = 4326)  # latlong coord reference system

# Limited preview
glimpse(select(osm_polys, osm_id, name, type, geometry))

## Rows: 77,337
## Columns: 4
## $ osm_id   <chr> "4763", "5922", "6022", "7141", "7316", "29744", "70978", "71…
## $ name     <chr> NA, NA, NA, NA, NA, "Ryde Canoe Lake", "Quarr Abbey", NA, "Is…
## $ type     <chr> "multipolygon", "multipolygon", "multipolygon", "multipolygon…
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-1.251425 5..., MULTIPOLYGON (((…

Of course, we can grab line features too. The default for layer is lines, so we don’t need to supply this argument.

# Fetch line features for IoW
osm_lines <-  oe_get(
  "Isle of Wight", 
  stringsAsFactors = FALSE,
  quiet = TRUE
) %>%
  st_transform(crs = 4326)

For convenience, I’m simplifying the features down to the ones we care about. For Tom, hedgerows were a constant nuisance, waterways were cold and perilous, and buildings could contain angry landowners. We want to avoid them all, ideally.

We can create a function to extract named features (we want rows containing ‘barriers’, ‘buildings’, ‘natural’ and ‘waterway’ features), and then iterate over our lines and polygons to isolate them. I’ve put them into a single object so it’s easier to reference them later.

# Filter for a single feature type from oe_get() output
# 'sf_in' is output from oe_get(); 'feature' is the feature we want
isolate_feature <- function(sf_in, feature) {
  sf_in %>% 
    filter(!is.na(sf_in[[feature]])) %>%  # filter for feature
    select(osm_id, type = all_of(feature), geometry)  # simplify object
}

# Get all the features as a list object with one element per feature
features <- map2(
  list(osm_lines, osm_polys, osm_polys, osm_lines),
  list("barrier", "building", "natural", "waterway"),
  isolate_feature
) %>% 
  set_names("barrs", "bldgs", "natur", "wways")

# Limited preview of the waterways data
glimpse(features$wways)

## Rows: 2,689
## Columns: 3
## $ osm_id   <chr> "3027797", "3127808", "4680059", "5171926", "5171930", "51722…
## $ type     <chr> "river", "river", "river", "stream", "ditch", "river", "strea…
## $ geometry <LINESTRING [°]> LINESTRING (-1.272943 50.62..., LINESTRING (-1.290…

The line

Of course, we need to specify a straight line route.² For this demonstration, and in interests fo speed, I’ve just chosen one that looks alright by eye in terms of obstructions.

Of course, you can use the approach outlined in this post to try other lines and discover quantitatively which ones have the fewest obstructions. That’s the subject of an upcoming Shiny app, which will allow the user to provide a line and get feedback on the number and count of hazards.

Crucially, the line is clipped to the boundary of the island, so it goes from coast to coast.

# Create a straight line
# 'x1', etc, are start/end latlongs; 'boundary_poly' is the GeoJSON
make_line <- function(x1, y1, x2, y2, boundary_poly) {
  x <- st_linestring(matrix(c(y1, y2, x1, x2), ncol = 2))  # to matrix
  y <- st_sfc(x, crs = 4326)  # set coord reference system
  st_intersection(y, boundary_poly)  # clip line to island boundary
}

# Hard-coded start/end latlongs
start_x <- 50.658; start_y <- -1.472
end_x <- 50.707; end_y <- -1.101

# Build line between the points, clip to IOW boundary
line <- make_line(start_x, start_y, end_x, end_y, iow_extent)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

# Preview
glimpse(line)

## sfc_LINESTRING of length 1; first list element:  'XY' num [1:2, 1:2] -1.47 -1.1 50.66 50.71

You’ll notice we get a warning about ‘planar’ coordinates. As far as I know (as a non-geospatial analyst), this is related to working with data that assume the earth is flat (!) versus curved. Of course, calculations will be slightly off if you assume no curvature. I’m not bothered with that on the scale of the IOW, so I’m ignoring it (and suppressing subsequent warnings in this post).

Here’s a quick preview of where our line is:

ggplot() +
  geom_sf(data = iow_extent) +
  geom_sf(data = line) +
  ggthemes::theme_map()

The platinum zone

In practice it’s very difficult to keep exactly to the line and some deviation may be required at a landowner’s request, for example. What’s an acceptable amount of wiggle room?

Tom spoke in one of his videos about assigning scores to a mission based on minimal deviation from the line. For example, staying within 25 m of the line either side would be a ‘platinum’ standard.

We can build a buffer around our line to create a platinum zone, which will be a useful visual aid in the final map.

# Create a buffer around the straight line
# 'line' as created with make_line(); 25m is 'platinum' standard
make_buffer <- function(straight_line, buffer_size = 25) { 
  x <- st_transform(straight_line, crs = 27700)  # line to draw buffer around
  y <- st_buffer(x, dist = buffer_size)
  st_transform(y, crs = 4326)  # reset coord reference system
}

# Generate a 25m buffer ('platinum standard') around the line
buffer <- make_buffer(line)

# Preview
glimpse(buffer)

## sfc_POLYGON of length 1; first list element: List of 1
##  $ : num [1:123, 1:2] -1.1 -1.1 -1.1 -1.1 -1.1 ...
##  - attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"

Obstructions

Our features object currently has all the lines and polygons within the IOW boundary, but we only want the ones that come in contact with (i.e. intersect) with our proposed route. The st_intersects() fcuntion from {sf} does exactly that.

# Extract features that intersect with the line or buffer geometry
# 'features_sf' contains the features; 'path_sf' for what to intersect
find_obs <- function(features_sf, path_sf) { 
  row_nums <- st_intersects(path_sf, features_sf)[[1]]  # rows with obstruction
  slice(features_sf, row_nums)  # extract matching rows
}

# Find intersection between features and line
obstructions <- map(features, ~find_obs(.x, line)) %>% 
  set_names("barrs", "bldgs", "natur", "wways")

# Limited preview
glimpse(obstructions$wways)

## Rows: 11
## Columns: 3
## $ osm_id   <chr> "471224575", "552565448", "552565555", "552565565", "55650679…
## $ type     <chr> "stream", "stream", "stream", "stream", "ditch", "ditch", "di…
## $ geometry <LINESTRING [°]> LINESTRING (-1.289448 50.67..., LINESTRING (-1.216…

We can create a quick table to see what the obstructions are for this route.

tribble(
  ~Type, ~Count,
  "Barriers", nrow(obstructions$barrs),
  "Buildings", nrow(obstructions$bldgs),
  "Waterways", nrow(obstructions$wways),
  "Water bodies", nrow(filter(obstructions$natur, type == "water")
  )
) %>% knitr::kable()

You can see how this information is useful if you wanted to try other straight-line paths and see how they stack up against each other. Maybe you want to reduce the number of barriers (typically hedgerows); maybe you don’t mind swimming across a large water body.

Map

So we have the objects we need: our straight-line route, the platinum-zone buffer and all the features that cross our path. Now to map it. We can use the {leaflet} package to layer these up on an interactive map, allowing the user to inspect the route and the hazards along the way.

# Set multi-use variables
marker_fill <- "darkblue"
icon_fill <- "white"
col_line <- "#000"
col_artifical <- "#F00"
col_water <- "#00F"
weight_line <- 1
weight_obstruction <- 2
alpha_all <- 0.5

# Interactive map
leaflet() %>% 
  # Base groups: map tiles
  addProviderTiles("CartoDB.PositronNoLabels", group = "Simple") %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Terrain") %>%
  addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
  # Overlay groups: line and buffer
  addPolygons(
    data = buffer, group = "Line/buffer", 
    color = col_line, weight = weight_line, dashArray = 4, 
    fill = TRUE, opacity = alpha_all
  ) %>% 
  addPolylines(
    data = line, group = "Line/buffer",
    color = col_line, weight = weight_line, 
    fill = FALSE
  ) %>%
  # Overlay groups: start and end points
  addAwesomeMarkers(
    lng = start_y, lat = start_x, group = "Start/end",
    icon = awesomeIcons(
      markerColor = "blue",
      icon = "play", library = "fa", iconColor = "#FFF"
    ),
    popup = paste0("<center>Start<br>", start_x, ", ", start_y, "<center>")
  ) %>% 
  addAwesomeMarkers(
    lng = end_y, lat = end_x, group = "Start/end",
    icon = awesomeIcons(
      markerColor = "blue",
      icon = "stop", library = "fa", iconColor = "#FFF"
    ),
    popup = paste0("<center>End<br>", end_x, ", ", end_y, "<center>")
  ) %>% 
  # Overlay groups: features in buffer
  addPolylines(
    data = obstructions$barrs, group = "Barriers",
    color = col_artifical, weight = weight_obstruction,
    label = paste("Barrier:", obstructions$barrs$type)
  ) %>% 
  addPolygons(
    data = obstructions$bldgs, group = "Buildings",
    color = col_artifical, weight = weight_obstruction, 
    fillColor = col_artifical, fillOpacity = alpha_all,
    label = paste("Building:", obstructions$bldgs$type)
  ) %>%
  addPolygons(
    data = filter(obstructions$natur, type == "water"), group = "Water",
    color = col_water, weight = weight_obstruction, 
    fillColor = col_water, fillOpacity = alpha_all,
    label = "Water body"
  ) %>% 
  addPolylines(
    data = obstructions$wways, group = "Water",
    color = col_water, weight = weight_obstruction,
    label = paste("Waterway:", obstructions$wways$type)
  ) %>% 
  # Control which layers are shown
  addLayersControl(
    baseGroups = c("Simple", "Terrain", "Satellite"),  # radio button
    overlayGroups = c(  # checkboxes
      "Line/buffer", "Start/end",  # line-related
      "Water", "Barriers", "Buildings"  # obstructions
    ),
    position = "topright",
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
  # Other map features
  addMeasure(  # tool for users to measure distances
    position = "topleft",
    primaryLengthUnit = "meters",
    primaryAreaUnit = "sqmeters"
  ) %>% 
  addFullscreenControl()  # clickable full-screen button

So, you can inspect the route interactively by zoom and dragging, by hovering over highlighted features to see what they are, and by turning on and off the different map and feature layers for more or less context.