CASA 0005 GIS and Science Seminar 1

Presenter: Nicolas Palominos
Duration: 60 minutes
Repository: https://github.com/npalomin/CASA_seminar_2020
URL: https://npalomin.github.io/CASA_seminar_2020/casa0005_seminar1.html

This seminar was presented on November 21, 2020 for the CASA GIS and Science Module at University College London. The topic presented was “Exploring the concept of 15-minute cities using service area analysis around train stations in London”

1 Learning outcomes

1. Understand and apply basic service area analysis in R
2. Conduct point-in-polygon queries
3. Use of the mapboxapi package
4. Understand basic concepts of accessibility analysis and diversity
5. Knowledge of the Points of Interest data set from Ordnance Survey

2 Introduction

The concept of 15-minutes cities has gained traction as a sustainable city planning strategy to create mixed-used and compact self-sufficient neighbourhoods. More recently, with the several urban mobility restrictions and lifestyle transformations originated from Covid-19, it has emerged as an urban planning framework for pandemic response and recovery.

The main idea behind the concept is that work, study, leisure, shopping and other basic human activities are available within a 15-minute trip mainly using active travel modes (walking and cycling).

This way long commutes are avoided originating important benefits to urban quality of life by both re-organising the use of dwellers’ time and reducing the externalities of urban transport among others.

Because of the predominant high dependency of private car transportation it is probable that many urban areas will not have the characteristics of good quality local environments across different urban domains (commercial and service activities, housing, public services, transport infrastructure, open and public spaces, employment).

Fig. Diagrammatic representation of a 15-minute city from ft.com

Fig. Multiple-dimensions of a 20-minute neighbourhood (Melbourne) from c40knowledgehub.org

Suburban areas often have less access to these features. Accessibility is a central concept in spatial analysis as it involves a cost often represented ‘spatially’ (distance, time), yet there is no agreed definition. For the purpose of this exercise, a useful definition is suggested by Batty (2009)

A well defined form of [accesibility] associates some measure of an opportunity at a place with the cost of actually realising that opportunity

Service area analysis finds all accessible streets (that is, streets that are within a specified impedance, cost expressed in time or distance). For example, the 10-minute service area for a given point on a network includes all the streets that can be reached within ten minutes from that point. Therefore, the geography of the street network is taken into account to calculate distance (see Gutierrez and Garcia-Palomares, 2008).

Fig. Service Area Isochrone vs ‘As-the-crow-Flies’ from Gutierrez and Garcia-Palomares, 2008. The isochrone shape (irregular polygon) is more accurate to represent equal cost of access (expressed in distance)

To conduct the analysis we will use Mapbox. The Mapbox platform provides mapping and navigation services with data mostly taken from open data sources such as OpenStreetMap. For example, the Isochrone API allows you to request polygon or line features that show areas that are reachable within a specified amount of time from a location.

#knitr::knit_exit()

3 Getting started

3.1 Install packages and load libraries

#install.packages("mapboxapi", dependencies = TRUE)
library(mapboxapi)
library(sf)
library(tidyverse)
library(tmap)
library(mapview) # to create fast interactive maps
library(ggplot2) # to create visualisations
library(ggthemes) # to style visualisations

3.2 Load and prepare data

For this exercise we will use the following data sets:

• Rail and Tube stations in London: We will assume this are activity nodes or points from which we will calculate the service area (from the tubecreature project by Oliver O’Brien at CDRC)
• POI Points of Interest: Directory of all Businesses and Services in Great Britain created by Ordnance Survey and available through Digimap licence (log in to https://digimap.edina.ac.uk/, navigate to Ordnance Survey > Data Download > Boundary and Location Data).

# derived from OSM, stations are represented as one point. Lines are also simplified and loaded here to provide context
stations <- st_read("https://raw.githubusercontent.com/oobrien/vis/master/tube/data/tfl_stations.json")

## Reading layer `tfl_stations' from data source `https://raw.githubusercontent.com/oobrien/vis/master/tube/data/tfl_stations.json' using driver `GeoJSON'
## Simple feature collection with 529 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -0.9729138 ymin: 51.33461 xmax: 0.3306072 ymax: 51.7469
## CRS:            4326

lines <- st_read("https://raw.githubusercontent.com/oobrien/vis/master/tube/data/tfl_lines.json")

## Reading layer `tfl_lines' from data source `https://raw.githubusercontent.com/oobrien/vis/master/tube/data/tfl_lines.json' using driver `GeoJSON'
## Simple feature collection with 700 features and 2 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -0.9729138 ymin: 51.11713 xmax: 0.3305912 ymax: 51.90153
## CRS:            4326

# get from data folder in case you can't download it yourself
# warning: the whole of the UK data set is ~1.5GB and aroung 4 million obs. We will crop to the bounding box of the study area
# the poi.rds is croped to the study area and has less attributes and saved as R object(RDS) to reduce its size
# to load a RDS straight from a URL requires a built-in decompressor function called gzcon plus the readRDS function
# this file is 10MB ~480k obs.
poi <- readRDS(gzcon(url("https://npalomin.github.io/casa0005/seminar/poi.rds")))
# Data Copyright Note
# Ordnance Survey Digimap: Licensed Data: © Crown copyright and database rights 2020 Ordnance Survey (100025252). This material includes data licensed from PointX© Database Right/Copyright 2020.

Quick plot to check Stations and Lines

qtm(lines) + qtm(stations)

Clip lines to the station’s bounding box

lines <-st_crop(lines,stations)

Check POI data

sample_n(poi, 10)

3.2.1 Point density visualisation

The POI data set for the study area contains N = 479360 obs. Below is a point density visualisation to have an overview of the general spatial pattern

# create new sf with CRS 4326 and X and Y variables
poi_4326 <- poi %>%
  #sample_n(100000) %>%
  st_transform(4326) %>% # tranform to same CRS as stations and lines
  cbind(st_coordinates(.)) # get X and Y coordinates
# plot
ggplot() +
  geom_bin2d(data = poi_4326, aes(X, Y), binwidth = c(0.01, 0.01)) + # heatmap of 2d bin counts
  geom_sf(data = lines, col="white", size=0.1) +
  theme_tufte() +
  scale_fill_viridis_c(option = "plasma") +
  labs(x="", y="")

3.3 Setting mapbox access token

In case you haven’t yet created a Mapbox account you can create one here

Access your account and create and copy the “default public token” that appears on your screen

# paste token here
# my_token <- "YOUR TOKEN GOES HERE"

4 Service area analysis

4.1 Generate service area polygons

# create a unique id variable "sid"
stations <- stations %>%
  rowid_to_column(., "sid")

# sample 5 stations
#s5s <- sample_n(stations, 100)
s5s <- stations
# central stations
#c5s <- stations %>%
#  slice(c(150,49,147,23,215,500))

To calculate the isochrones we use mb_isochrone function with the following parameters

• location(s)
• profile = “walking”, “driving”, “cycling”
• time in minutes
• id = column identifier from location
• access_token = my_token

w_iso <- mb_isochrone(s5s, "walking", time = 10, id_column = "sid", access_token = my_token)

# get variables from station sf to join with isochrones
stations_kv <- stations %>%
  select("sid","name", "zone") %>%
  st_drop_geometry()
# join
w_iso <- w_iso %>%
  merge(., stations_kv, by.x="id", by.y="sid") %>%
  rename(s_name = name) # to avoid conflict with POI name

Plot a sample of isochrones

# select 9 random isochrones
r9s <- sample(w_iso$id, 9)

w_iso_9 <- w_iso %>%
  filter(id %in% r9s)

s9 <- stations %>%
  filter(sid %in% r9s)

tm_shape(w_iso_9) +
  tm_polygons(col = "tomato",
             border.lwd = 0.3,
             alpha = 0.7) +
  tm_facets(by="s_name", ncol = 3) +
tm_shape(s9) +
  tm_symbols(col = "navy",
             size = 0.2) +
tm_shape(lines) +
  tm_lines(lwd=0.2) +
  tm_layout(legend.show = F, outer.margins = 0, panel.label.bg.color = 'white')

4.2 Analyse service areas

4.2.1 What do we know about the POI attributes?

Plot a sample of POI

sample_POI <- sample_n(poi_4326, 500)
# handy package/function to plot interactive map
mapview(sample_POI)

How many variables do each of the columns have?

poi_4326 %>%
  st_drop_geometry() %>%
  select(groupname, categoryname, classname) %>%
  mutate(n_groupname = n_distinct(groupname),
         n_categoryname = n_distinct(categoryname),
         n_classname = n_distinct(classname)) %>%
  select(n_groupname, n_categoryname, n_classname) %>%
  unique()

As it is often the case with classification schemes, these are structured hierarchically. The ‘classname’ and ‘categoryname’ variables are nested into the 9 levels of the ‘groupname’ variable.

levels(poi_4326$groupname) 

## [1] "Accommodation, Eating and Drinking" "Attractions"                       
## [3] "Commercial Services"                "Education and Health"              
## [5] "Manufacturing and Production"       "Public Infrastructure"             
## [7] "Retail"                             "Sport and Entertainment"           
## [9] "Transport"

This summary will give us the count of data points in each level. A simple analysis of the service areas would be to count the number of observations within the area, but it would be interesting to examine to what group do those points belong. Is the accessibility to points of interest measured by counting the number of observations?, it worth having a deeper look into the kind of points are at reach within an area.

# how many observations do you have for 'groupname'
poi_4326 %>%
  st_drop_geometry() %>%
  count(groupname, sort = TRUE)

Now we can look at the composition of the groups (groupname variable). This table shows the amount of different values in each group

poi_4326 %>%
  st_drop_geometry() %>% 
  select(groupname, categoryname) %>% 
  unique() %>%
  count(groupname)

Here we are a step from building a dendogram, which is a netowork representation of hierarchical relationships such as the classification scheme of POI that we are studying.

# Libraries
library(ggraph)
library(igraph)
# create a hierarchy from data frame to examine classification scheme
poi_hier <- poi_4326 %>%
  st_drop_geometry() %>%
  select(groupname, categoryname, classname) %>%
  unique() 

rownames(poi_hier) <- NULL

# transform it to a edge list
edges_level0_1 <- poi_hier %>%
  select(groupname) %>% 
  mutate(head = "POI") %>%
  unique() %>% 
  rename(from=head, to=groupname) %>%
  select(2,1)

edges_level1_2 <- poi_hier %>% 
  select(groupname, categoryname) %>% 
  unique() %>% 
  rename(from=groupname, to=categoryname)

edges_level2_3 <- poi_hier %>% 
  select(categoryname, classname) %>% 
  unique() %>% rename(from=categoryname, to=classname)

edge_list  <- rbind(edges_level0_1, edges_level1_2)
edge_list  <- rbind(edge_list, edges_level2_3)

# Now we can plot the graph object
poi_graph <- graph_from_data_frame(edge_list)

ggraph(poi_graph, layout = 'dendrogram', circular = FALSE) + 
  geom_edge_diagonal(alpha = 0.1, size=0.5) +
  geom_node_point(color="tomato", size=0.01) +
  geom_node_text(aes(label=c("", poi_graph[[1]]$POI$name, rep("",662))), size=2.5) +
  theme_void() +
  coord_flip() +
  scale_y_reverse() 

Hierarchical struture of the POIs ‘group’, ‘categories’, ‘class’ parent-child relationships. This can help to better understand the richness of urban features available in the study area

4.2.2 POI counts per service area

# Join station identifier to POI and select only those intersecting with the isochrones (omit points with NA)
# this joins the isochrone data to the POIs
poi_sa <- poi_4326 %>%
  st_join(., w_iso) %>%
  na.omit() # omit POI that do not intersect with any service area

poi_count <- poi_sa %>%
  group_by(id, s_name) %>%
  summarise(n=n())
poi_count

# create a single point per service area to plot
poi_count_centroid <- st_centroid(poi_count)

Create a proportional symbol map weighted by count to examine the opportunities available within the study areas.

#tmap_mode("plot")
tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(poi_count_centroid) +
  tm_bubbles(col = "tomato",
             size = "n",
             alpha = 0.4,
             border.lwd = 0.3,
             style = "kmeans",
             title.size="Count") +
  tm_layout(title = "Points of Interest within a 10-minute walking service area from rail stations in London", title.size = 0.77) +
  tm_credits(size = 3, text = "Ordnance Survey Digimap: Licensed Data: © Crown copyright and database rights 2020 Ordnance Survey (100025252). This material includes data licensed from PointX© Database Right/Copyright 2020\nMap created by Nicolas Palominos", position = c("right", "BOTTOM"))

Let us investigate the richness of POIs in each area by looking at the type of class, category or group the points belong to. Richness is a measure of diversity which can help assess the variety of opportunities that are accessible within the study area which is more complex than the absolute count.

This table summarises the number of different ‘groups’ in each study area. The great majority of the service areas have POIs that belong to the 9 groups.

#
poi_sa %>%
  st_drop_geometry() %>%
  group_by(id, groupname, s_name) %>%
  summarise(n=n()) %>% # total number of observations per group
  group_by(s_name) %>%
  mutate(richness=n_distinct(groupname)) %>% # distinct number of groups
  select(3,5) %>% 
  unique()

What is the spatial pattern of the 9 different groups?

# create a sf object with the 'group' summary
poi_count_gr <- poi_sa %>%
  group_by(id, groupname) %>% # groupname count by station
  summarise(n=n())
poi_count_gr

# create a single point per service area to plot
poi_count_gr_cen <- st_centroid(poi_count_gr)

tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(poi_count_gr_cen) +
  tm_bubbles(col = "tomato",
             size = "n",
             alpha = 0.4,
             border.lwd = 0.3,
             style = "kmeans",
             title.size="Count") +
  tm_facets(by="groupname") +
  tm_layout(legend.show = F, outer.margins = 0, panel.label.bg.color = 'white', frame = F)

4.2.3 Diversity analysis: point types per service area

All the nine groups (Accommodation, Eating and Drinking; Attractions; Commercial Services; Education and Health; Manufacturing and Production; Public Infrastructure; Retail; Sport and Entertainment; Transport) seem to be relevant for a self-sufficient neighbourhood. Let’s quantify the richness of these groups by looking at the type of points that are located in each service area.

# create summary of obs. by station
pps <- poi_sa %>%
  st_drop_geometry() %>%
  count(id) %>%
  rename(station_count = n)

# create a 'Richness' metric by 'Group', "Category' and 'Class'
div_tab <- poi_sa %>%
  st_drop_geometry() %>%
  group_by(id, s_name) %>%
  summarise(rich_g = n_distinct(groupname),
            rich_cat = n_distinct(categoryname),
            rich_cla = n_distinct(classname))

# merge
div_tab <- div_tab %>%
  merge(., pps)

# merge to stations to get geometry back (the approach above was different, we created centroids)
sta_div <- stations %>%
  select("sid","name", "zone") %>%
  merge(., div_tab, by.x="sid", by.y="id")

Plot it

tmap_mode("plot")
tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(sta_div) +
  tm_symbols(col = "rich_cla",
             palette="RdPu",
             size = 0.1,
             shape = "station_count",
             alpha = 0.8,
             border.lwd = 0.3,
             style = "kmeans",
             n = 3,
             shapes.n = 3,
             shapes.legend.fill = "gray30",
             title.shape="Count",
             title.col="Diversity") +
  tm_layout(title = "Diversity of Points of Interest within a 10-minute walking service area from rail stations in London", main.title = "", title.size = 0.77) +
  tm_credits(size = 3, text = "Ordnance Survey Digimap: Licensed Data: © Crown copyright and database rights 2020 Ordnance Survey (100025252). This material includes data licensed from PointX© Database Right/Copyright 2020\nMap created by Nicolas Palominos", position = c("right", "BOTTOM"))

According to the pattern displayed in previous analysis, centrally located stations have a good provision of urban facilities according to both number and variety. Zone 1 in London more or less matches the planning definition of the Central Activity Zone which is effectively not only a self-sufficient area but also a provider of services for the whole of London. For a more meaningful analysis of potentially under-served areas we will filter out zone 1, and zoom-in to compare the top and tail areas according to richness quantified by the number of distinct ‘classes’ located within the service area.

# get stations with lowest diversity
tail <- sta_div %>%
  # keep all that are not zone = 1
  filter(zone != "1") %>%
  arrange(rich_g, rich_cat, rich_cla) %>% # arrange is ascendant by default (order of sorting factors affect result)
  slice(1:3)
tail

# get stations with highest diversity
head <- sta_div %>%
  # keep all that are not zone = 1
  filter(zone != "1") %>%
  arrange(desc(rich_g, rich_cat, rich_cla)) %>%
  slice(1:3)
head

# get POI and isochrones from sf by station id (sid)

# tail
tail_iso <- w_iso %>%
  filter(id %in% tail$sid)
tail_poi <- poi_sa %>%
  filter(id %in% tail$sid)

# head
head_iso <- w_iso %>%
  filter(id %in% head$sid)
head_poi <- poi_sa %>%
  filter(id %in% head$sid)

Plot head and tail

tmap_mode("plot")
tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(head) +
  tm_dots(col = "tomato",
             size = 0.3,
             alpha = 0.8,
             border.lwd = 0.3) +
  tm_text("rich_cla", size = 0.7, just=c("center", 2)) +
tm_shape(head) +
  tm_dots(col = "tomato",
             size = 0.1,
             alpha = 0.8,
             border.lwd = 0.3) +
  tm_text("name", size = 0.5, just=c("center", 4)) +
tm_shape(tail) +
  tm_dots(col = "pink",
             size = 0.3,
             alpha = 0.8,
             border.lwd = 0.3) +
  tm_text("rich_cla", size = 0.7, just=c("center", 2)) +
tm_shape(tail) +
  tm_dots(col = "pink",
             size = 0.1,
             alpha = 0.8,
             border.lwd = 0.3) +
  tm_text("name", size = 0.5, just=c("center", 4)) +
tm_shape(stations) +
  tm_dots(col = "gray20",
             size = 0.007,
             alpha = 0.3,
             border.lwd = 0.001) +
  tm_layout(title = "10-minute walking service areas with higher and lower diversity according to number of 'classes'", inner.margins = c(0.08,0,0,0), title.size = 0.77) 

Zooming-in

# create custom palette
g_col <- c('#e6194B', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#42d4f4', '#f032e6', '#fabed4', '#469990')

tmap_mode("plot")
head_m <- tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(head_iso) +
  tm_polygons(col = "grey90",
             border.lwd = 0.3,
             alpha = 0.7) +
tm_facets(by="s_name", nrow = 1) +
tm_shape(head_poi) +
  tm_dots(col = "groupname",
          palette = g_col,
             size = 0.2,
             alpha = 0.8,
             border.lwd = 0.3,
             title="Group") +
  tm_layout(title = "", frame=FALSE, panel.label.bg.color = 'white', main.title = "Higher diversity", main.title.size = 0.8) +
tm_legend(show = FALSE) 

tail_m <- tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(tail_iso) +
  tm_polygons(col = "grey90",
             border.lwd = 0.3,
             alpha = 0.7) +
tm_facets(by="s_name", nrow = 1) +
tm_shape(tail_poi) +
  tm_dots(col = "groupname",
          palette = g_col,
             size = 0.2,
             alpha = 0.8,
             border.lwd = 0.3,
             title="Group") +
  tm_layout(title = "", frame=FALSE, panel.label.bg.color = 'white', main.title = "Lower diversity", main.title.size = 0.8) +
tm_legend(show = FALSE)

legend <- tm_shape(lines) +
  tm_lines(lwd=0.1) +
tm_shape(tail_iso) +
  tm_polygons(col = "grey90",
             border.lwd = 0.3,
             alpha = 0.7) +
tm_shape(tail_poi) +
  tm_dots(col = "groupname",
          palette = g_col,
             size = 0.2,
             alpha = 0.8,
             border.lwd = 0.3,
             title="Group") +
  tm_layout(title = "", frame=FALSE, panel.label.bg.color = 'white', main.title = "Lower diversity", main.title.size = 0.8, legend.only =T)

tmap_arrange(head_m, legend, tail_m, nrow = 2, widths = c(0.75, 0.25))

What is the increase in diversity if we change the mode of transport to define the extent of the service area in the areas with lowest diversity?

# low diversity service area stations
ldsas <- stations %>%
  filter(sid %in% tail$sid)

# get isochrones for cycling mode
w_iso_low <- mb_isochrone(ldsas, "cycling", time = 10, id_column = "sid", access_token = my_token)

# join stations data
w_iso_low <- w_iso_low %>%
  merge(., stations_kv, by.x="id", by.y="sid") %>%
  rename(s_name = name) # to avoid conflict with POI name

#Join station identifier to POI and select only those intersecting with the isochrones (omit points with NA)
poi_low <- poi_sa %>%
  select(-id, -time, -s_name, -zone) %>% # already in poi_sa from previous join
  st_join(., w_iso_low) %>%
  na.omit() # omit POI that do not intersect with service area

# create a 'Richness' metric by 'Group', "Category' and 'Class'
div_low <- poi_low %>%
  st_drop_geometry() %>%
  group_by(id, s_name) %>%
  summarise(rich_g = n_distinct(groupname),
            rich_cat = n_distinct(categoryname),
            rich_cla = n_distinct(classname))

# merge to stations to get geometry back (the approach above was different, we created centroids)
sta_div <- stations %>%
  select("sid","name", "zone") %>%
  merge(., div_low, by.x="sid", by.y="id")

# print summary table: lowest diversity cycling mode
sta_div

tm_shape(lines) +
  tm_lines(lwd=0.3) +
tm_shape(tail_iso) +
  tm_polygons(col = "beige",
             border.lwd = 0.3,
             alpha = 0.8) +
tm_shape(poi_low) +
  tm_dots(col = "groupname",
          palette = g_col,
             size = 0.2,
             alpha = 0.8,
             border.lwd = 0.3,
             title="Group") +
tm_shape(w_iso_low) +
  tm_polygons(col = "beige",
             border.lwd = 0.3,
             border.col = "tomato",
             alpha = 0.1) +
tm_facets(by="s_name", nrow = 1) +
  tm_layout(title = "", frame=FALSE, panel.label.bg.color = 'white', main.title = "Lower diversity according to 'classes' type\nWalking and Cycling 10-minutes service area definitions", main.title.size = 0.8, legend.outside=TRUE, legend.outside.position = "bottom", legend.outside.size=0.25) +
  tm_legend(show =T)

The 10-minutes cycling service area of Coombe Lane has a POI diversity richness of 204 according to distinct classes which is in some way comparable to the 10-minutes walking service area of Archway, which has a richness of 175. However, further questions remain about other variables not considered in this analysis that characterise a 15-minute city, such as the quality of the street environment or the housing and demographic diversity, among others.

5 Final comments

• We used service area analysis to create a new geographic area definition to perform a comparative analysis of the number and diversity of urban facilities surrounding train stations.

• For what other applications could this research method could be useful for? For example, instead of train stations the centre points could be schools to investigate the kind of urban facilities students have access to.

• We used a licensed high quality data set of Point of Interest available for the the whole of the UK. However, a similar analysis could be replicated elsewhere using OpenStreetMap data (e.g. examining amenity and shop features).

• We only considered one variable to investigate the concept of 15-minute cities, however, it would have been interesting to incorporate other variables more directly related to quantifying the cost of trips using active travel modes.

• While we used the geography of the street network to define the study areas, there was no quantification of the characteristics of the street environment. There are plenty of data sets that could contribute to expand this analysis, such as street greenery, street hierarchy, cycle lanes and streetspace allocation.

6 Feedback

If you have any feedback or comments please feel free to contact me at @npalomin or n.palominos.16@ucl.ac.uk