this document descibes the process of calcualting walking and cycling exposure for a number of typical routes within the London Borough of Waltham Forest.

Preparation

Load libraries

rm(list = ls())

## Load libraries
library('knitr')
library("googleway")
library("gepaf")
library("sp")
library('raster')
library('ggplot2')
library('geosphere')
library('sf')
## Linking to GEOS 3.5.0, GDAL 2.2.2, proj.4 4.9.2
library('reshape2')
library('leaflet')
library('mapview')
library('kableExtra')

## rmarkdown::render('make_exposure.Rmd', output_dir = "docs", output_file ="exposure.html")

Data import

Import routes

exposure_routes <- read.csv('exposure_journeys.csv',
                            stringsAsFactors = F)

exposure_routes[exposure_routes$via_array == '','via_array'] <- NA

Import colours

source('https://raw.githubusercontent.com/KCL-ERG/colour_schemes/master/pm25_laei2013_colours_breaks.R')
source('https://raw.githubusercontent.com/KCL-ERG/colour_schemes/master/no2_laei2013_colours_breaks.R')
source('https://raw.githubusercontent.com/KCL-ERG/colour_schemes/master/pm10_laei2013_colours_breaks.R')

Route calculation, offest, and exposure calculation

Make a loop to calculate the routes and extract the concentrations. If seggregated==true, then the route is offset by the variable ‘offset’

start_time        <- Sys.time() + 1000
via_array         <- NA


for (i in 1:nrow(exposure_routes)) {

# Journey parameters
start_lat         <- exposure_routes[i,]$start_lat
start_lon         <- exposure_routes[i,]$start_lon

if (!is.na(exposure_routes[i,]$via_array)) {via_array <- as.list(strsplit(as.character(exposure_routes[i,]$via_array), ";")[[1]]) } else {via_array <- NA}

end_lat           <- exposure_routes[i,]$end_lat
end_lon           <- exposure_routes[i,]$end_lon
mode              <- as.character(exposure_routes[i,]$mode)
seggregated       <- exposure_routes[i,]$seggregated

## take account of where in the road the journey is taking place
if (mode == 'walk'    & seggregated == 'no')  {offset <- 6}
if (mode == 'walk'    & seggregated == 'yes') {offset <- 8}
if (mode == 'bicycle' & seggregated == 'no')  {offset <- 4}
if (mode == 'bicycle' & seggregated == 'yes') {offset <- 6}


# EPSG strings we might need
latlong                 <- "+init=epsg:4326"
ukgrid                  <- "+init=epsg:27700"
google_projected        <- "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

# Take the start/end data and get a route from Google
if ('walk' %in% mode)     { source("routing/walking.R") }
if ('bicycle' %in% mode)  { source("routing/cycling.R") }

rm(start_lat, start_lon, end_lat, end_lon)

if ('walk' %in% mode) {
  result <- walk_result
  rm(walk_result)
  duration <- walk_duration
  rm(walk_duration)
} else {
  result <- cycle_result
  rm(cycle_result)
  duration <- cycle_duration
  rm(cycle_duration)
}

##### In here is where we calculate the offset for where the journey is taking place compared to the centre of the road
##### https://stackoverflow.com/questions/50275195/draw-a-parallel-line-in-r-offset-from-a-line
#### Only want to do this if seggregated == yes

coordinates(result) <- ~lon + lat
proj4string(result) =  CRS(latlong)
result <- spTransform(result, ukgrid)

segment.shift <- function(x, y, d){
#  
  # calculate vector
  v <- c(x[2] - x[1],y[2] - y[1])
#  
  # normalize vector
  v <- v/sqrt((v[1]**2 + v[2]**2))
#  
  # perpendicular unit vector
  vnp <- c( -v[2], v[1] )
  
  return(list(x =  c( x[1] + d*vnp[1], x[2] + d*vnp[1]), 
              y =  c( y[1] + d*vnp[2], y[2] + d*vnp[2])))
  
}

x <- result@coords[,1]
y <- result@coords[,2]

xn <- numeric( (length(x) - 1) * 2 )
yn <- numeric( (length(y) - 1) * 2 )

for ( p in 1:(length(x) - 1) ) {
  xs <- c(x[p], x[p+1])
  ys <- c(y[p], y[p+1])
  new.s <- segment.shift( xs, ys, offset)
  xn[(p-1)*2+1] <- new.s$x[1] ; xn[(p-1)*2+2] <- new.s$x[2]
  yn[(p-1)*2+1] <- new.s$y[1] ; yn[(p-1)*2+2] <- new.s$y[2]
}

new_result   <- data.frame(x = xn,
                           y = yn)
coordinates(new_result) <- ~x+y
proj4string(new_result) = CRS(ukgrid)
new_result              <- spTransform(new_result, latlong)
new_result              <- data.frame(new_result)
new_result$id           <- unique(result$id)
new_result$line         <- unique(result$line)
new_result$mode         <- unique(result$mode)

result <- new_result
rm(new_result, new.s, offset, x, xn, xs, y, yn, ys, segment.shift)
names(result)[1:2] <- c('lon', 'lat')

##### End of the help from stackoverflow

coordinates(result) <- ~lon + lat
proj4string(result) = CRS(latlong)
result <- spTransform(result, ukgrid) 

## This makes a spatialline from the points for the result
x <- lapply(split(result, result$id), function(x) Lines(list(Line(coordinates(x))), result$id[1L]))
lines <- SpatialLines(x)
data <- data.frame(id = unique(result$id))
rownames(data) <- data$id
result <- SpatialLinesDataFrame(lines, data)
proj4string(result) = CRS(ukgrid)
rm(x, lines, data) 


## How many minutes long was the bicycle and walk journeys

timeslots                <- data.frame(id = unique(result$id), time = seq.POSIXt(start_time, start_time + duration, "min"))

## Now want to split the line into a point per minute

result                    <- st_as_sf(result)
result$id                 <- exposure_routes[i,]$journey_id

if (i == 1) {line_result <- result} else {line_result <- rbind(line_result, result)}

result                   <- st_line_sample(result, n = nrow(timeslots), type = 'regular')
result                   <- suppressWarnings(st_set_crs(result, 27700))
result                   <- data.frame(as(result, 'Spatial'))

## Now some harmonisation stuff to get the car, walk and bus all the same format so can join them together
result$id                 <- exposure_routes[i,]$journey_id
result$mode               <- mode
result$line               <- NA
result                    <- cbind(result, timeslots)
names(result)             <- c('easting', 'northing', 'id', 'mode', 'line', 'id_2', 'time')
result                    <- result[,c('id', 'mode', 'line', 'time', 'easting', 'northing')]
result$year               <- exposure_routes[i,]$year
coordinates(result)       <- ~easting + northing
proj4string(result)       <- CRS(ukgrid)

# Now get our pollutant files

## 2013 concentration maps
if (exposure_routes[i,]$year == 2013) {
wf_2013_no2              <- raster('air_quality/no2_2013.asc')
crs(wf_2013_no2)         <- CRS(ukgrid)
wf_2013_nox              <- raster('air_quality/nox_2013.asc')
crs(wf_2013_nox)         <- CRS(ukgrid)
wf_2013_pm10             <- raster('air_quality/pm10_2013.asc')
crs(wf_2013_pm10)        <- CRS(ukgrid)
wf_2013_pm25             <- raster('air_quality/pm25_2013.asc')
crs(wf_2013_pm25)        <- CRS(ukgrid)

# Extract concentrations
no2                                   <- result
no2@data$pollutant                    <- 'no2'
no2@data$concentration                <- extract(wf_2013_no2, no2@coords, method = 'bilinear')

pm25                                  <- result
pm25@data$pollutant                   <- 'pm25'
pm25@data$concentration               <- extract(wf_2013_pm25, pm25@coords, method = 'bilinear')

pm10                                  <- result
pm10@data$pollutant                   <- 'pm10'
pm10@data$concentration               <- extract(wf_2013_pm10, pm10@coords, method = 'bilinear')

nox                                   <- result
nox@data$pollutant                    <- 'nox'
nox@data$concentration                <- extract(wf_2013_nox, nox@coords, method = 'bilinear')

result                                <- rbind.SpatialPointsDataFrame(no2, pm25, nox, pm10)

rm(wf_2013_no2, wf_2013_nox, wf_2013_pm10, wf_2013_pm25,nox, no2, pm25, pm10)

}

## 2020 concentration maps 
if (exposure_routes[i,]$year == 2020) {
wf_2020_no2              <- raster('air_quality/no2_2020.asc')
crs(wf_2020_no2)         <- CRS(ukgrid)
wf_2020_nox              <- raster('air_quality/nox_2020.asc')
crs(wf_2020_nox)         <- CRS(ukgrid)
wf_2020_pm10             <- raster('air_quality/pm10_2020.asc')
crs(wf_2020_pm10)        <- CRS(ukgrid)
wf_2020_pm25             <- raster('air_quality/pm25_2020.asc')
crs(wf_2020_pm25)        <- CRS(ukgrid)

# Extract concentrations
no2                                   <- result
no2@data$pollutant                    <- 'no2'
no2@data$concentration                <- extract(wf_2020_no2, no2@coords, method = 'bilinear')

pm25                                  <- result
pm25@data$pollutant                   <- 'pm25'
pm25@data$concentration               <- extract(wf_2020_pm25, pm25@coords, method = 'bilinear')

pm10                                  <- result
pm10@data$pollutant                   <- 'pm10'
pm10@data$concentration               <- extract(wf_2020_pm10, pm10@coords, method = 'bilinear')

nox                                   <- result
nox@data$pollutant                    <- 'nox'
nox@data$concentration                <- extract(wf_2020_nox, nox@coords, method = 'bilinear')

result                                <- rbind.SpatialPointsDataFrame(no2, pm25, nox, pm10)

rm(wf_2020_no2, wf_2020_nox, wf_2020_pm10, wf_2020_pm25,nox, no2, pm25, pm10)

}

if (i == 1) {final_result <- result} else {final_result <- rbind.SpatialPointsDataFrame(final_result, result)}

rm(timeslots, duration, result, via_array)

}

Bind the results we’ve made together with the original parameter data, and delete some variables that aren’t needed anymore

final_result$time <- NULL
final_result$line <- NULL
final_result      <- merge(final_result, exposure_routes[,c('journey_id', 'name', 'seggregated')], by.x = 'id', by.y = 'journey_id')
rm(i, google_projected, latlong, mode, seggregated, ukgrid)

Results

The summary of the routes and years are shown.

summary_results <- aggregate(concentration ~ name + mode + year + pollutant + seggregated, data=final_result, FUN=mean)

names(summary_results) <- c('Route', 'Mode', 'Year', 'Pollutant', 'Seggregated', 'Mean Exposure (ug/m3)')

summary_results[order(summary_results$Route, summary_results$Pollutant, summary_results$Year),] %>% kable(row.names = F) %>% kable_styling()
Route Mode Year Pollutant Seggregated Mean Exposure (ug/m3)
Chingford to Chingford Police walk 2013 no2 no 40.53042
Chingford to Chingford Police walk 2020 no2 no 33.97841
Chingford to Chingford Police walk 2013 nox no 76.93478
Chingford to Chingford Police walk 2020 nox no 57.44137
Chingford to Chingford Police walk 2013 pm10 no 25.62571
Chingford to Chingford Police walk 2020 pm10 no 24.09539
Chingford to Chingford Police walk 2013 pm25 no 15.77246
Chingford to Chingford Police walk 2020 pm25 no 14.31240
Chingford to Leyton bicycle 2013 no2 no 53.72370
Chingford to Leyton bicycle 2020 no2 no 43.11749
Chingford to Leyton bicycle 2020 no2 yes 42.14660
Chingford to Leyton bicycle 2013 nox no 113.64550
Chingford to Leyton bicycle 2020 nox no 76.74637
Chingford to Leyton bicycle 2020 nox yes 74.32698
Chingford to Leyton bicycle 2013 pm10 no 28.43883
Chingford to Leyton bicycle 2020 pm10 no 26.76959
Chingford to Leyton bicycle 2020 pm10 yes 26.51731
Chingford to Leyton bicycle 2013 pm25 no 17.07426
Chingford to Leyton bicycle 2020 pm25 no 15.42368
Chingford to Leyton bicycle 2020 pm25 yes 15.34176
Coppermill to Wood Street bicycle 2013 no2 no 35.87851
Coppermill to Wood Street bicycle 2020 no2 no 29.82928
Coppermill to Wood Street bicycle 2013 nox no 59.58672
Coppermill to Wood Street bicycle 2020 nox no 43.94401
Coppermill to Wood Street bicycle 2013 pm10 no 25.04390
Coppermill to Wood Street bicycle 2020 pm10 no 23.46135
Coppermill to Wood Street bicycle 2013 pm25 no 15.89599
Coppermill to Wood Street bicycle 2020 pm25 no 14.41685
Green Man Roundabout to Leytonstone High Road Station walk 2013 no2 no 54.59607
Green Man Roundabout to Leytonstone High Road Station walk 2020 no2 no 43.18944
Green Man Roundabout to Leytonstone High Road Station walk 2013 nox no 114.81455
Green Man Roundabout to Leytonstone High Road Station walk 2020 nox no 75.16567
Green Man Roundabout to Leytonstone High Road Station walk 2013 pm10 no 28.28507
Green Man Roundabout to Leytonstone High Road Station walk 2020 pm10 no 26.62365
Green Man Roundabout to Leytonstone High Road Station walk 2013 pm25 no 17.15719
Green Man Roundabout to Leytonstone High Road Station walk 2020 pm25 no 15.49845
Lea Bridge to Whipps Cross bicycle 2013 no2 no 58.64818
Lea Bridge to Whipps Cross bicycle 2020 no2 no 44.59887
Lea Bridge to Whipps Cross bicycle 2020 no2 yes 43.86417
Lea Bridge to Whipps Cross bicycle 2013 nox no 132.07544
Lea Bridge to Whipps Cross bicycle 2020 nox no 79.32771
Lea Bridge to Whipps Cross bicycle 2020 nox yes 77.42049
Lea Bridge to Whipps Cross bicycle 2013 pm10 no 29.30678
Lea Bridge to Whipps Cross bicycle 2020 pm10 no 27.49929
Lea Bridge to Whipps Cross bicycle 2020 pm10 yes 27.30874
Lea Bridge to Whipps Cross bicycle 2013 pm25 no 17.52780
Lea Bridge to Whipps Cross bicycle 2020 pm25 no 15.76842
Lea Bridge to Whipps Cross bicycle 2020 pm25 yes 15.70525
Leytonstone to Lea Bridge via Ruckolt bicycle 2013 no2 no 49.29758
Leytonstone to Lea Bridge via Ruckolt bicycle 2020 no2 no 40.39847
Leytonstone to Lea Bridge via Ruckolt bicycle 2013 nox no 94.53851
Leytonstone to Lea Bridge via Ruckolt bicycle 2020 nox no 68.04068
Leytonstone to Lea Bridge via Ruckolt bicycle 2013 pm10 no 27.91970
Leytonstone to Lea Bridge via Ruckolt bicycle 2020 pm10 no 26.62305
Leytonstone to Lea Bridge via Ruckolt bicycle 2013 pm25 no 17.03712
Leytonstone to Lea Bridge via Ruckolt bicycle 2020 pm25 no 15.53458
Leytonstone to Stratford Drapers bicycle 2013 no2 no 58.02299
Leytonstone to Stratford Drapers bicycle 2020 no2 no 47.56987
Leytonstone to Stratford Drapers bicycle 2013 nox no 119.82710
Leytonstone to Stratford Drapers bicycle 2020 nox no 85.59571
Leytonstone to Stratford Drapers bicycle 2013 pm10 no 29.52559
Leytonstone to Stratford Drapers bicycle 2020 pm10 no 28.43967
Leytonstone to Stratford Drapers bicycle 2013 pm25 no 17.64129
Leytonstone to Stratford Drapers bicycle 2020 pm25 no 16.14284
Leyton to Blackhorse bicycle 2013 no2 no 56.37080
Leyton to Blackhorse bicycle 2020 no2 no 44.37406
Leyton to Blackhorse bicycle 2013 nox no 121.52751
Leyton to Blackhorse bicycle 2020 nox no 78.48987
Leyton to Blackhorse bicycle 2013 pm10 no 28.74930
Leyton to Blackhorse bicycle 2020 pm10 no 26.98430
Leyton to Blackhorse bicycle 2013 pm25 no 17.32166
Leyton to Blackhorse bicycle 2020 pm25 no 15.60624
Leyton to Drapers Fields walk 2013 no2 no 54.47035
Leyton to Drapers Fields walk 2020 no2 no 44.00757
Leyton to Drapers Fields walk 2013 nox no 111.94394
Leyton to Drapers Fields walk 2020 nox no 77.46930
Leyton to Drapers Fields walk 2013 pm10 no 28.57642
Leyton to Drapers Fields walk 2020 pm10 no 26.99013
Leyton to Drapers Fields walk 2013 pm25 no 17.24248
Leyton to Drapers Fields walk 2020 pm25 no 15.59055
Walthamstow Central to W altham Town Hall walk 2013 no2 no 50.31994
Walthamstow Central to W altham Town Hall walk 2020 no2 no 40.78683
Walthamstow Central to W altham Town Hall walk 2013 nox no 102.75337
Walthamstow Central to W altham Town Hall walk 2020 nox no 71.51593
Walthamstow Central to W altham Town Hall walk 2013 pm10 no 27.50929
Walthamstow Central to W altham Town Hall walk 2020 pm10 no 25.82083
Walthamstow Central to W altham Town Hall walk 2013 pm25 no 16.79701
Walthamstow Central to W altham Town Hall walk 2020 pm25 no 15.17257
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2013 no2 no 54.74540
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2020 no2 no 43.17686
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2013 nox no 117.34390
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2020 nox no 76.28655
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2013 pm10 no 28.54720
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2020 pm10 no 26.80878
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2013 pm25 no 17.21467
Walthamstow to Lea Bridge via Selbourne and Markhouse walk 2020 pm25 no 15.53184
Wood Street to Blackhorse bicycle 2013 no2 no 41.73367
Wood Street to Blackhorse bicycle 2020 no2 no 35.00155
Wood Street to Blackhorse bicycle 2013 nox no 77.24893
Wood Street to Blackhorse bicycle 2020 nox no 57.72396
Wood Street to Blackhorse bicycle 2013 pm10 no 25.98101
Wood Street to Blackhorse bicycle 2020 pm10 no 24.38233
Wood Street to Blackhorse bicycle 2013 pm25 no 16.24164
Wood Street to Blackhorse bicycle 2020 pm25 no 14.72489
Wood Street to Waltham Town Hall walk 2013 no2 no 46.98477
Wood Street to Waltham Town Hall walk 2020 no2 no 39.23297
Wood Street to Waltham Town Hall walk 2013 nox no 89.84098
Wood Street to Waltham Town Hall walk 2020 nox no 66.92266
Wood Street to Waltham Town Hall walk 2013 pm10 no 27.20727
Wood Street to Waltham Town Hall walk 2020 pm10 no 25.54458
Wood Street to Waltham Town Hall walk 2013 pm25 no 16.74899
Wood Street to Waltham Town Hall walk 2020 pm25 no 15.15275 
final_result    <- st_as_sf(final_result)
summary_results <- merge(summary_results, exposure_routes[,c('journey_id', 'name')], by.x = 'Route', by.y = 'name')

summary_results <- merge(summary_results, line_result, by.x = 'journey_id', by.y ='id')

summary_results <- st_sf(summary_results[,-c(8)],st_sfc(summary_results$geometry))

Now make some plots of the results

plot <- list()

for (i in 1:length(unique(final_result$name))) {
  
  route_name <- unique(final_result$name)[i]
  
  plot[[i]] <- ggplot(final_result[final_result$name == route_name & final_result$pollutant %in% c('no2', 'pm25'),], 
                      aes(x = as.factor(year), y = concentration, fill = seggregated)) +
    geom_boxplot() +
    facet_wrap(name~pollutant, ncol = 4, scales = 'free_y') +
    xlab('Year') +
    theme(plot.subtitle = element_text(vjust = 1), 
    plot.caption = element_text(vjust = 1), 
    panel.grid.major = element_line(colour = "gray92"), 
    panel.grid.minor = element_line(colour = "gray92"), 
    axis.title = element_text(size = 14), 
    axis.text = element_text(size = 14, colour = "black"), 
    axis.text.y = element_text(size = 14), 
    plot.title = element_text(size = 14), 
    legend.title = element_text(size = 13), 
    panel.background = element_rect(fill = NA, 
        colour = "black", linetype = "solid"), 
    plot.background = element_rect(colour = NA),
    strip.text.x = element_text(size = 12),
    strip.background = element_blank())

}

Example 1 - Chingford to Leyton (Concentration plots)

plot[[1]]

Example 1 - Chingford to Leyton (Spatial plots)

final_result <- st_set_crs(final_result, 27700)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
final_result <- st_transform(final_result, 4326)

example_2013 <- final_result[final_result$pollutant == 'no2' & final_result$name == 'Chingford to Leyton' & final_result$year == 2013,]

example_2020 <- final_result[final_result$pollutant == 'no2' & final_result$name == 'Chingford to Leyton' & final_result$year == 2020 & final_result$seggregated == 'yes',]

qpal <- colorQuantile(no2_laei2013_colours, c(example_2013$concentration, example_2020$concentration), n = 16)

qpal <- colorNumeric(
  palette = no2_laei2013_colours,
  domain = c(example_2013$concentration, example_2020$concentration))

leaflet(options = leafletOptions(attributionControl=F))%>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2013,
              group = '2013',
              fillColor = qpal(example_2013$concentration),
              color="transparent",
              fillOpacity = 0.7,
             popup = paste0(as.character(round(example_2013$concentration,0)),' ?g/m<sup>3</sup>')) %>%
  addFeatures(example_2020,
              group = '2020',
              fillColor = qpal(example_2020$concentration),
              color="transparent",
              fillOpacity = 0.7,
              popup = paste0(as.character(round(example_2020$concentration,0)), ' ?g/m<sup>3</sup>')) %>%
  addLayersControl(baseGroups = c('2013', '2020')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="NO2 (?g/m<sup>3</sup>)") %>%
  addControl('<strong>Exposure to NO<sub>2</sub> while cycling between Chingford <br> & Leyton has dropped by 22% between 2013<br>and 2020 due to improvements in air<br> quality, and the introduction of seggregated <br> cycle lanes.</strong>', position = "topleft")

Now make a static output of the same.

no2_2013_chingford_leyton <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2013,
              group = '2013',
              fillColor = qpal(example_2013$concentration),
              color="transparent",
              fillOpacity = 0.9,
             popup = paste0(as.character(round(example_2013$concentration,0)),' ?g/m<sup>3</sup>')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="NO2 (µg/m<sup>3</sup>)")

mapshot(no2_2013_chingford_leyton, file = "no2_2013_chingford_leyton.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))

no2_2020_chingford_leyton <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2020,
              group = '2020',
              fillColor = qpal(example_2020$concentration),
              color="transparent",
              fillOpacity = 0.9,
             popup = paste0(as.character(round(example_2020$concentration,0)),' ?g/m<sup>3</sup>')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="NO2 (µg/m<sup>3</sup>)")

mapshot(no2_2020_chingford_leyton, file = "no2_2020_chingford_leyton.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))

Example 2 - Leyton to Blackhorse (Spatial plots)

example_2013 <- final_result[final_result$pollutant == 'no2' & final_result$name == 'Leyton to Blackhorse' & final_result$year == 2013,]

example_2020 <- final_result[final_result$pollutant == 'no2' & final_result$name == 'Leyton to Blackhorse' & final_result$year == 2020,]

qpal <- colorQuantile(no2_laei2013_colours, c(example_2013$concentration, example_2020$concentration), n = 16)

qpal <- colorNumeric(
  palette = no2_laei2013_colours,
  domain = c(example_2013$concentration, example_2020$concentration))

leaflet(options = leafletOptions(attributionControl=F))%>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2013,
              group = '2013',
              fillColor = qpal(example_2013$concentration),
              color="transparent",
              fillOpacity = 0.7,
             popup = paste0(as.character(round(example_2013$concentration,0)),' µg/m<sup>3</sup>')) %>%
  addFeatures(example_2020,
              group = '2020',
              fillColor = qpal(example_2020$concentration),
              color="transparent",
              fillOpacity = 0.7,
              popup = paste0(as.character(round(example_2020$concentration,0)), ' µg/m<sup>3</sup>')) %>%
  addLayersControl(baseGroups = c('2013', '2020')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="NO2 (µg/m<sup>3</sup>)") %>%
  addControl('<strong>Exposure to NO<sub>2</sub> while cycling between Leyton <br> & Blackhorse has dropped by 21% between 2013<br>and 2020 due to improvements in air<br> quality, and the introduction of seggregated <br> cycle lanes.</strong>', position = "bottomleft")

Now make a static output of the same.

no2_2013_leyton_blackhorse <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2013,
              group = '2013',
              fillColor = qpal(example_2013$concentration),
              color="transparent",
              fillOpacity = 0.9,
             popup = paste0(as.character(round(example_2013$concentration,0)),' ?g/m<sup>3</sup>')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="NO2 (µg/m<sup>3</sup>)")

mapshot(no2_2013_leyton_blackhorse, file = "no2_2013_leyton_blackhorse.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))

no2_2020_leyton_blackhorse <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2020,
              group = '2020',
              fillColor = qpal(example_2020$concentration),
              color="transparent",
              fillOpacity = 0.9,
             popup = paste0(as.character(round(example_2020$concentration,0)),' ?g/m<sup>3</sup>')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="NO2 (µg/m<sup>3</sup>)")

mapshot(no2_2020_leyton_blackhorse, file = "no2_2020_leyton_blackhorse.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))

PM2.5 from Chingford to Leyton

example_2013 <- final_result[final_result$pollutant == 'pm25' & final_result$name == 'Leyton to Blackhorse' & final_result$year == 2013,]

example_2020 <- final_result[final_result$pollutant == 'pm25' & final_result$name == 'Leyton to Blackhorse' & final_result$year == 2020,]

qpal <- colorQuantile(pm25_laei2013_colours, c(example_2013$concentration, example_2020$concentration), n = 11)

qpal <- colorNumeric(
  palette = pm25_laei2013_colours,
  domain = c(example_2013$concentration, example_2020$concentration))

pm25_2013_chingford_leyton <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2013,
              group = '2013',
              fillColor = qpal(example_2013$concentration),
              color="transparent",
              fillOpacity = 0.9,
             popup = paste0(as.character(round(example_2013$concentration,0)),' µg/m<sup>3</sup>')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="PM2.5 (µg/m<sup>3</sup>)")

mapshot(pm25_2013_chingford_leyton, file = "pm25_2013_leyton_blackhorse.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))

pm25_2020_chingford_leyton <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(example_2020,
              group = '2020',
              fillColor = qpal(example_2020$concentration),
              color="transparent",
              fillOpacity = 0.9,
             popup = paste0(as.character(round(example_2020$concentration,0)),' µg/m<sup>3</sup>')) %>%
  addLegend(pal=qpal,values=c(example_2013$concentration, example_2020$concentration),title="PM2.5 (µg/m<sup>3</sup>)")

mapshot(pm25_2020_chingford_leyton, file = "pm25_2020_leyton_blackhorse.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))

All routes

line_result <- st_set_crs(line_result, 27700)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
line_result <- st_transform(line_result, 4326)

leaflet(options = leafletOptions(attributionControl=F))%>%
  #addProviderTiles(providers$Esri.WorldTopoMap)%>%
  #addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
  addTiles() %>%
  addFeatures(line_result[line_result$id %in% exposure_routes[exposure_routes$mode == 'bicycle',]$journey_id,],
              group = 'Bike Routes',
              color = 'red',
              popup = as.character(exposure_routes[exposure_routes$journey_id %in% line_result[line_result$id %in% exposure_routes[exposure_routes$mode == 'bicycle',]$journey_id,]$id,]$name),
              opacity = 0.3) %>%
    addFeatures(line_result[line_result$id %in% exposure_routes[exposure_routes$mode == 'walk',]$journey_id,],
                group = 'Walk Routes',
                color = 'blue',
                popup = as.character(exposure_routes[exposure_routes$journey_id %in% line_result[line_result$id %in% exposure_routes[exposure_routes$mode == 'walk',]$journey_id,]$id,]$name),
                 opacity = 0.3) %>%
  addLayersControl(baseGroups = c('Bike Routes', 'Walk Routes'))

Chingford Station to Leyton Station

temp_static <- leaflet(options = leafletOptions(attributionControl=F)) %>%
  addProviderTiles(providers$Esri.WorldTopoMap)%>%
  addFeatures(line_result[line_result$id == 1,])

mapshot(temp_static, file = "route_chingford_leyton.png",
        remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar"))
st_write(final_result, 'exposure_routes.geojson', delete_dsn=TRUE)
## Deleting source `/home/james/github/waltham_forest/exposure_routes.geojson' using driver `GeoJSON'
## Writing layer `exposure_routes' to data source `/home/james/github/waltham_forest/exposure_routes.geojson' using driver `GeoJSON'
## features:       1944
## fields:         7
## geometry type:  Point

Result available