menai-bridge-cycle.Rmd

---
title: "Cycling potential across the Menai Strait"
author: "Robin Lovelace on behalf of Beicio Bangor"
date: "`r Sys.Date()`"
output: rmarkdown::word_document
vignette: >
  %\VignetteIndexEntry{Cycling potential across the menai Strait}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: cycling-potential.bib
---

```{r, echo=FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
```

This report explores cycling potential across the Menai Strait, using the 'propensity to cycle tool' (PCT) method [@lovelace_propensity_2016]. The PCT is freely available for England thanks to a Department for Transport project at [www.pct.bike](http://pct.bike/) but has yet to be created for Wales.

To ensure reproducibility, and encourage the use of the methods for other areas, the code and data needed to reproduce this analysis have been made available [online, at github.com/Robinlovelace/pct-menai ](https://github.com/Robinlovelace/pct-menai/blob/master/vignettes/menai-bridge-cycle.Rmd). The remainder of this report discusses the input data, the current distribution of cycling in the study area, commuter cycling potential and cyclable commutes originating in Anglesey, which would have the greatest congestion benefits. Overall it was found that bridges crossing the Menai Strait have substantial commuter cycling potential. Trips for education, shopping and socialising are likely to have greater cycling potential due to economic and cultural exchange between the settlements of Porthaethwy and Llanfairpwllgwyngyll on Anglesey with Bangor and its surroundings on the mainland.

## Data

```{r preprocessing, echo=FALSE, include=FALSE}
# devtools::install_github("robinlovelace/pctWales")
library(pctWales)
library(stplanr)
library(tmap)
library(rgeos)
library(ggplot2)
data("wales")
data("od_wales")
data("od_wales_lines")
menai = wales[grep("Anglesey|Gwynedd", wales$name),]
angle = wales[grep("Anglesey", wales$name),]
proj4string(od_wales_lines) = proj4string(menai)
menai_outline = rgeos::gBuffer(menai, width = 0)
menai_outline = stplanr::buff_geo(menai, width = 1000)
# plot(menai_outline)
sel = gContains(menai_outline, od_wales_lines, byid = TRUE)
l = od_wales_lines[rowSums(sel) > 0,]
# mapview(l)

# subsetting by proximity to the menai bridge
# mgeo = tmap::geocode_OSM("menai bridge")
# data("mgeo")
mpt = SpatialPoints(matrix(mgeo$coords, ncol = 2), CRS(proj4string(wales)))
menai_5k = buff_geo(mpt, 5000)
menai_15k = buff_geo(mpt, 15000)
# plot(menai_10k)
# plot(l, add = T)
cents_5km = 
od_menai = l[rowSums(gContains(menai_15k, l, TRUE)) > 0,]
# plot(od_menai, add = T, col = "red")
# data(od_menai)
```

```{r, echo=FALSE}
l_twoways = od_menai[menai_straight,]
# ltf = line2route(l_twoways)
# devtools::use_data(ltf)
df = l_twoways@data
df$Distance = ltf$length / 1000
dbands = c(0, 5, 10, 15, 20, 40)
df$`Distance band` = cut(df$Distance, dbands)
library(dplyr)
dfdb =
  group_by(df, `Distance band`) %>% 
  summarise(All = sum(`All categories: Method of travel to work`),
            Cycle = sum(`Bicycle`),
            Drive = sum(`Driving a car or van`),
            Passenger = sum(`Passenger in a car or van`),
            Walk = sum(`On foot`)
              )
dfdb$Other = dfdb$All - rowSums(dfdb[3:ncol(dfdb)])
dfdb$`Distance band` = as.character(dfdb$`Distance band`)
dfdb = rbind(dfdb, c("All", colSums(dfdb[2:7])))
```


The input dataset for this report was origin-destination (OD) data from the 2011 Census. The number of people travelling by mode, cross tabulated by distance bands, is illustrated in Table 1.

```{r, echo=FALSE}
knitr::kable(dfdb, caption = "Numbers of people travelling by mode and route distance bands")
```


For the aggregate analysis, the OD dataset was preprocessed so that a single line represents travel in both directions (for the one way analysis one way flows were used). The dataset was converted to geographical lines with start and end points corresponding to population weighted middle layer super output areas (MSOAs). This work, and estimation of route distance, was calculated using the **stplanr** R package.

To highlight the fact that these are geographic entities, we will use the term 'flows' to refer to this data. 'Desire lines' is another commonly used term for these straight lines but we use the term flows here for brevity. The starting point of this report was to subset the relevant flows, in two stages:

1. Identify all flows with start *and* end points within a 15 km buffer of the Menai bridge.

2. Identify all flows that cross the Menai Strait, the commuters of which may potentially use the bridge.

These stages are represented in the figure below.

```{r selection, echo=FALSE, fig.cap="Flows selected for analysis of cycling potential across the Menai Strait. Grey lines represent flows that begin and end within 15 km of the Menai Bridge. ", warning=FALSE}
l = onewaygeo(l_twoways, attrib = 3:14)
l$`N. commuters (100s)` = l$`All categories: Method of travel to work` /
                           100
l$pcycle = l$Bicycle /  l$`All categories: Method of travel to work`
# osm_tiles = read_osm(bb(menai_15k, ext = 1.2))
# use_data(osm_tiles)
data("osm_tiles")
m = tm_shape(osm_tiles) + tm_raster() +
  tm_shape(menai_15k) + tm_borders() +
  tm_shape(od_menai) + tm_lines(col = "grey") +
  tm_shape(menai_straight) + tm_lines(col = "red", lwd = 3) +
  tm_shape(l) + tm_lines(col = "blue", lwd = "All categories: Method of travel to work", scale = 2) +
  tm_scale_bar()
m
save_tmap(m, "selection.pdf")
```

It is worth describing briefly this dataset. `r nrow(l)` flows were selected by the criteria above. The total number of commuters travelling to work along these flows by all modes, who could potentially use bridges crossing the Strait, was `r sum(l[[3]])`. Of these, `r sum(l$Bicycle)` reported using a pedal cycle as their main method of travel to work,
`r round(100 * sum(l$Bicycle) / sum(l[[3]]), 1)`%.

```{r, warning=FALSE, echo=FALSE}
lc1 = l[l$Bicycle >= 10,]
sel = l$Bicycle >= 5
lc4 = l[sel,]
lc3 = l[l$Bicycle >= 5 & l$Bicycle < 10,]
# lc1$Bicycle / sum(l$Bicycle)
# sum(lc4$Bicycle) / sum(l$Bicycle)
# sum(lc3$Bicycle) / sum(l$Bicycle)
```

\newpage

## Current levels of cycle commuting across the Menai Strait

Of those who report cycling to work in the flows plotted above, 47% (32 commuter cyclists) reported travelling between the single flow connecting the MSOAs of Bangor and Porthaethwy. 23% of reported travelling in one of the 3 flows connecting MSOAs centred in Bangor and Glanadda with MSOAs containing Beaumaris and Llangaffo. These flows (which together constitute 70% of commuter cycling among the subsetted flows) are represented in Figure 2.

```{r, echo=FALSE, fig.cap="Desire lines crossing the Menai Strait frequented by 5 or more cyclists in the 2011 census (blue) and all desire lines crossing the straight (green). Widths in both cases are proportional to the number of commutes by all modes."}
menai_5k = buff_geo(mpt, width = 5000)
l5 = l[menai_5k,]
l5$`N. Commuters (within 5 km)` = l5$`All categories: Method of travel to work`
lc4$`N. Commuters (more than 5 commuter cyclists)` = lc4$`All categories: Method of travel to work`
lc4$text = paste0(lc4$`All categories: Method of travel to work`, ", ", lc4$Bicycle, " cycle")
# osm_t4 = read_osm(bb(lc4, ext = 1.6))
# devtools::use_data(osm_t4)
m = tm_shape(osm_t4) + tm_raster() +
  tm_shape(l5) + tm_lines(col = "darkgreen", lwd = "N. Commuters (within 5 km)", scale = 2) +
  tm_shape(lc4) + tm_lines(col = "blue", lwd = "N. Commuters (more than 5 commuter cyclists)", scale = 4) +
  tm_shape(SpatialPointsDataFrame(coords = gCentroid(lc4, byid = T), data = lc4@data)) +
  tm_text(text = "text") +
  tm_scale_bar()
m
save_tmap(m, "current.pdf")
```

\newpage

## Potential cycle commuters

Of the `r sum(l[[3]])` commuters whose desire line crosses the Menai Strait, barely 1% of cycled. Is this due to the distance of average commuter trips and hills in the region? Or does inadequate infrastructure play a role? The graph below shows the trip distance-frequency distribution of the selected flows. This shows that there are 2 large flows under 5 km in distance and a number of flows between 10 and 15 km. Note that the numbers presented do not account for within-zone flows (people who live and work in the same zone) and people who have no fixed workplace.

<!-- , which accounted for  -->

```{r, echo=FALSE}
p = line2points(l)
p1 = p[seq(1, length.out = nrow(l), by = 2),]
p2 = p[seq(2, length.out = nrow(l), by = 2),]
lw = spTransform(l, CRS("+init=epsg:27700"))
# plot(gLength(lw, byid = T) / 1000, l$`All categories: Method of travel to work`)
```

```{r cycling-potential, echo=FALSE, fig.cap="Level of commuter cycling in the 2011 census (black) and under Government Target (blue) and Go Dutch (red) scenarios, omitting the impact of route gradient. Each point represents a desire line. The diameter of the points is proportional to the number of commuters by all modes for each desire line."}
# ??logit
distance = lr$length / 1000
gradient = lr$av_incline * 100

logit_pcycle = -3.894 + (-0.5872 * distance) + (1.832 * sqrt(distance) ) + (0.007956 * distance^2)
lr$govtarget = boot::inv.logit(logit_pcycle) + l$pcycle
logit_pcycle_dutch = logit_pcycle + 2.499 -0.07384 * distance
lr$godutch = boot::inv.logit(logit_pcycle_dutch)
lr$`Go Dutch` = lr$godutch * l$`All categories: Method of travel to work`
lr_nograd = lr

ggplot(lr@data) +
  geom_point(aes(distance, godutch * 100), color = "red",
             size = l$`All categories: Method of travel to work` / 100) +
  geom_point(aes(distance, govtarget * 100), color = "blue",
             size = l$`All categories: Method of travel to work` / 100) +
  geom_point(aes(distance, l$pcycle * 100), color = "black",
             size = l$`All categories: Method of travel to work` / 100) +
  xlab("Route distance (km)") + ylab("Proportion cycling to work (%)") +
  xlim(c(0,25))

# plot(distance, lr$govtarget)
```

```{r, echo=FALSE, fig.cap="Level of commuter cycling in the 2011 census (black) and under Government Target (blue) and Go Dutch (red) and E-bike scenarios, with the impact of hilliness included. Each point represents a desire line. The diameter of the points is proportional to the number of commuters by all modes for each desire line."}
# with hilliness
logit_pcycle = -3.894 + (-0.5872 * distance) + (1.832 * sqrt(distance) ) + (0.007956 * distance^2) +
(-0.2872 * gradient) + (0.01784 * distance * gradient) + (-0.09770 * sqrt(distance) * gradient)
lr$govtarget = boot::inv.logit(logit_pcycle) + l$pcycle
# plot(distance, lr$govtarget)
logit_pcycle_dutch = logit_pcycle + 2.499 -0.07384 * distance
lr$godutch = boot::inv.logit(logit_pcycle_dutch) + l$pcycle

logit_pcycle_ebike = logit_pcycle_dutch + (0.05710 * distance) + (-0.0001087 * distance^2) + (-0.67 * -0.2872 * gradient)
l$ebike = boot::inv.logit(logit_pcycle_ebike) 


ggplot(lr@data) +
  geom_point(aes(distance, godutch * 100), color = "red",
             size = l$`All categories: Method of travel to work` / 100) +
  geom_point(aes(distance, govtarget * 100), color = "blue",
             size = l$`All categories: Method of travel to work` / 100) +
  geom_point(aes(distance, l$pcycle * 100), color = "black",
             size = l$`All categories: Method of travel to work` / 100) +
  geom_point(aes(distance, l$ebike * 100), color = "green",
             size = l$`All categories: Method of travel to work` / 100) +
  xlab("Route distance (km)") + ylab("Proportion cycling to work (%)") +
  xlim(c(0,25))

# lr$length
  # (0.05710 * ebike * distance) + (-0.0001087 * ebike * distance sq ) + (-0.67 * -0.2872 *
# ebike * gradient).
lr$`Go Dutch` = lr$godutch * l$`All categories: Method of travel to work`
```

It is clear that the model, which was trained against English data, has a strong aversion to hills. Note that part of this aversion could be due to the fact that hilly places generally have a low level of cycling. For this reason results with and without hilliness included are presented.

To allocate these flows to the transport network, CycleStreets.net was used to find routes converting the straight lines in the 'fastest route' and the R package **stplanr** was used to aggregate overlapping routes. Cyclestreets.net does not recognize the Britannia Bridge as a cycleable route; all flows are therefore allocated to the Menai Bridge, and the cycleability of routes which would be optimal across the Britannia Bridge is systematically under-estimated. The results are illustrated in the map below.

```{r, echo=FALSE, message=FALSE, fig.cap="Cycle commutes allocated to the road network under the Go Dutch scenario, accounting for gradient."}
rnet = overline(lr_nograd, "Go Dutch")
# lr = line2route(l)
# devtools::use_data(lr)
# lb = line2route(l, plan = "balanced")
# plot(rnet, lwd = rnet$`Go Dutch` / 100)
# tmap_mode("view")
m = tm_shape(osm_t4) + tm_raster() +
  tm_shape(rnet) + tm_lines(col = "red", lwd = "Go Dutch", scale = 8) +
  tm_scale_bar()
m
# save_tmap(m, "rnet-dutch-nograd.pdf")
# save_tmap(m, "vignettes/rnet-dutch-nograd.png")
```

<!-- ```{r, echo=FALSE, message=FALSE, fig.cap="Cycle commutes allocated to the road network under the Go Dutch scenario, accounting for gradient."} -->
<!-- knitr::include_graphics("rnet-dutch-nograd.pdf") -->
<!-- ``` -->


```{r, echo=FALSE, message=FALSE, fig.cap="Cycle commutes allocated to the road network under the Go Dutch scenario, accounting for the impact of gradient in the English model."}
rnet = overline(lr, "Go Dutch")
# lr = line2route(l)
# devtools::use_data(lr)
# lb = line2route(l, plan = "balanced")
# plot(rnet, lwd = rnet$`Go Dutch` / 100)
m = tm_shape(osm_t4) + tm_raster() +
  tm_shape(rnet) + tm_lines(col = "red", lwd = "Go Dutch", scale = 8) +
  tm_scale_bar()
m
save_tmap(m, "rnet-dutch.pdf")
```

\newpage

## Flows originating in Anglesey

The rush hour traffic is worse on the mainland-bound routes over the Menai Strait. Thus cycling uptake among people travelling in this direction in the morning could be particularly beneficial. This is illustrated in the maps below, which show substantial potential for cycling update along these congested routes.

```{r, include=FALSE}
# plot(angle)
proj4string(angle) = proj4string(welsh_msoas)
angle_msoas = welsh_msoas[angle,]
# plot(angle_msoas)
# plot(l_twoways, add = T)
l_from_angle = l_twoways[l_twoways$`Area of residence` %in% angle_msoas$geo_code,]
# plot(l_from_angle, add = T, col = "red")
# nrow(l_from_angle)
# summary(l_from_angle$`All categories: Method of travel to work`)
# mapview::mapview(l_from_angle)
# summary(l_from_angle)
lfa = l_from_angle[l_from_angle$`All categories: Method of travel to work` > 50,]
# nrow(lfa)
# summary(lfa)
# rfa = line2route(lfa)
# devtools::use_data(rfa)
rfa@data = cbind(rfa@data, lfa@data)
rfa$Distance = rfa$length / 1000
rfa$pcycle = rfa$Bicycle / rfa$`All categories: Method of travel to work`
rfa_sub = rfa[rfa$Distance < 15,]
# nrow(rfa_sub)
# mapview::mapview(rfa_sub)
```

```{r uptake-rfa, echo=FALSE}
distance = rfa_sub$length / 1000
gradient = rfa_sub$av_incline * 100

logit_pcycle = -3.894 + (-0.5872 * distance) + (1.832 * sqrt(distance) ) + (0.007956 * distance^2)
rfa_sub$govtarget = boot::inv.logit(logit_pcycle) + rfa_sub$pcycle
logit_pcycle_dutch = logit_pcycle + 2.499 -0.07384 * distance
rfa_sub$godutch = boot::inv.logit(logit_pcycle_dutch)
rfa_sub$`Go Dutch` = rfa_sub$godutch * rfa_sub$`All categories: Method of travel to work`
rfa_sub_nograd = rfa_sub
```

```{r, echo=FALSE}
# with hilliness
logit_pcycle = -3.894 + (-0.5872 * distance) + (1.832 * sqrt(distance) ) + (0.007956 * distance^2) +
(-0.2872 * gradient) + (0.01784 * distance * gradient) + (-0.09770 * sqrt(distance) * gradient)
rfa_sub$govtarget = boot::inv.logit(logit_pcycle) + rfa_sub$pcycle
# plot(distance, rfa_sub$govtarget)
logit_pcycle_dutch = logit_pcycle + 2.499 -0.07384 * distance
rfa_sub$godutch = boot::inv.logit(logit_pcycle_dutch) + rfa_sub$pcycle

logit_pcycle_ebike = logit_pcycle_dutch + (0.05710 * distance) + (-0.0001087 * distance^2) + (-0.67 * -0.2872 * gradient)
rfa_sub$ebike = boot::inv.logit(logit_pcycle_ebike) 

rfa_sub$`Go Dutch` = rfa_sub$godutch * rfa_sub$`All categories: Method of travel to work`
```

```{r, echo=FALSE, message=FALSE, fig.cap="Cycle commutes from Anglesey under the Go Dutch scenario, not accounting for gradient."}
rnet = overline(rfa_sub_nograd, "Go Dutch")
# rfa_sub = line2route(l)
# devtools::use_data(rfa_sub)
# lb = line2route(l, plan = "balanced")
# plot(rnet, lwd = rnet$`Go Dutch` / 100)
m = tm_shape(osm_t4) + tm_raster() +
  tm_shape(rnet) + tm_lines(col = "red", lwd = "Go Dutch", scale = 8) +
  tm_scale_bar()
m
# save_tmap(m, "rnet-dutch-nograd.pdf")
# save_tmap(m, "vignettes/rnet-dutch-nograd.png")
```

<!-- ```{r, echo=FALSE, message=FALSE, fig.cap="Cycle commutes allocated to the road network under the Go Dutch scenario, accounting for gradient."} -->
<!-- knitr::include_graphics("rnet-dutch-nograd.pdf") -->
<!-- ``` -->


```{r, echo=FALSE, message=FALSE, fig.cap="Cycle commutes from Anglesey under the Go Dutch scenario, accounting for the impact of gradient in the English model."}
rnet = overline(rfa_sub, "Go Dutch")
# rfa_sub = line2route(l)
# devtools::use_data(rfa_sub)
# lb = line2route(l, plan = "balanced")
# plot(rnet, lwd = rnet$`Go Dutch` / 100)
m = tm_shape(osm_t4) + tm_raster() +
  tm_shape(rnet) + tm_lines(col = "red", lwd = "Go Dutch", scale = 8) +
  tm_scale_bar()
m
save_tmap(m, "rnet-dutch.pdf")
```

\newpage

## Discussion

```{r, eval=FALSE, echo=FALSE}
nrow(l)
nrow(lr)
l$Distance = lr$length / 1000
sum(l$`All categories: Method of travel to work`)
sel = l$Distance < 15
sum(l$`All categories: Method of travel to work`[sel])
```

The analysis presented in this report does not tell the full storey of cycling potential over the Menai Strait. True cycling numbers could be boosted substantially by trips for education, leisure and shopping, that are not represented in Census travel to work data. However, the analysis has shown that many (5757) people crossed the passage for work each day in 2011 with start and end points within 15 km of the Menai Bridge.
Of these, over half (3070) have an inferred trip distance of less than 15 km (9 miles). If only a small portion of these trips were replaced by cycling, as the model suggests is feasible, the benefits for congestion and health could be great.

Further work could involve the analysis of the potential impacts of electric cycles, which could further boost the uptake of cycling, and higher geographical resolution of OD data.


## References