load.Rmd

---
title: "National Propensity to Cycle Tool - local results"
author: "Created by the NPCT team"
output:
  html_document:
    fig_caption: yes
    highlight: pygments
    theme: null
    toc: yes
---

```{r, include=FALSE}
start_time <- Sys.time() # for timing the script
source("set-up.R") # load packages needed
```

This document was produced automatically at `r start_time`.

## Introduction

The results of the National Propensity to Cycle Tool (NPCT) scenarios are based on a model.
This document presents information about the input data, model diagnostics,
run time and the key outputs for each region.
The aim is to inform users of the NPCT's interactive map how the results were generated.
This document assumes some technical knowledge in the area of transport planning.

The code underlying the NPCT is open source, under the conditions of the MIT License.
Both the [pct](https://github.com/npct/pct) and [pct-shiny](https://github.com/npct/pct-shiny) and be modified by others provided attribution to the original.

## Initial parameters

The preset values used to select the study area and filter the origin destination data are described in this section.

```{r, echo=FALSE}
# Create default LA name if none exists
if(!exists("region")) region <- "liverpool-city-region"
pct_data <- file.path("..", "pct-data")
pct_bigdata <- file.path("..", "pct-bigdata")
pct_privatedata <- file.path("..", "pct-privatedata")
pct_shiny_regions <- file.path("..", "pct-shiny", "regions_www")
if(!file.exists(pct_data)) stop(paste("The pct-data repository has cannot be found.  Please clone https://github.com/npct/pct-data in", dirname(getwd())))
if(!file.exists(pct_bigdata)) stop(paste("The pct-bigdata repository has cannot be found.  Please clone https://github.com/npct/pct-bigdata in", dirname(getwd())))
if(!file.exists(pct_privatedata)){
  message("The pct_privatedata repository cannot be found")
  scens <- c("govtarget_slc", "dutch_slc", "ebike_slc")
}else{
  scens <- c("govtarget_slc", "gendereq_slc", "dutch_slc", "ebike_slc")
}
if(Sys.getenv("CYCLESTREET") == "")
  warning("No CycleStreet API key. See ?cyclestreet_pat")
```

```{r, warning=FALSE}
# Set local authority and ttwa zone names
region # name of the region
region_path <- file.path(pct_data, region)
if(!dir.exists(region_path)) dir.create(region_path) # create data directory

# Minimum flow between od pairs, subsetting lines. High means fewer lines
mflow <- 20
mflow_short <- 1
mdist <- 20 # maximum euclidean distance (km) for subsetting lines
max_all_dist <- 7 # maximum distance (km) below which more lines are selected
buff_dist <- 0 # buffer (km) used to select additional zones (often zero = ok)
buff_geo_dist <- 250 # buffer (m) for removing line start and end points for network
# save the initial parameters to reproduce results
save(region, mflow, mflow_short, mdist, max_all_dist, buff_dist, buff_geo_dist, file = file.path(region_path, "params.RData"))
```

## Input zone data

The input zones area are summarised in this section.

```{r plotzones, message=FALSE, warning=FALSE, results='hide', echo=FALSE}
if(!exists("ukmsoas")) ukmsoas <- shapefile(file.path(pct_bigdata, "msoas.shp"))
# Load population-weighted centroids
if(!exists("centsa"))
  centsa <- readOGR(file.path(pct_bigdata, "cents.geojson"), "OGRGeoJSON")
centsa$geo_code <- as.character(centsa$geo_code)
las <- readOGR(dsn = file.path(pct_bigdata, "las-pcycle.geojson"), layer = "OGRGeoJSON")
las_cents <- SpatialPoints(coordinates(las))

# Load local authorities and districts
if(!exists("geo_level")) geo_level <- "regions"
if(geo_level == "cua"){
  regions <- readOGR(file.path(pct_bigdata, "regions.geojson"), layer = "OGRGeoJSON")
} else if(geo_level == "custom"){
  # if you use a custom geometry, regions should already be saved from buildmaster.R
} else{
    regions <- readOGR(dsn = file.path(pct_bigdata, "cuas-mf.geojson"), layer = "OGRGeoJSON")
  regions$Region <- regions$CTYUA12NM
}

# create region shape (and add buffer in m)
region_shape <- region_orig <-
  regions[grep(pattern = region, x = regions$Region, ignore.case = T),]

# Only transform if needed
if(buff_dist > 0){
  region_shape <- spTransform(region_shape, CRS("+init=epsg:27700"))
  region_shape <- gBuffer(region_shape, width = buff_dist * 1000)
  region_shape <- spTransform(region_shape, proj4string(centsa))
}

las_in_region <- rgeos::gIntersects(las_cents, region_shape, byid = T)
las_in_region <- las_in_region[1,]
las_in_region <- las[las_in_region,]
```

The selected region is `r region`.

```{r, include=FALSE}
# select msoas of interest
if(proj4string(region_shape) != proj4string(centsa))
  region_shape <- spTransform(region_shape, proj4string(centsa))
cents <- centsa[region_shape,]
zones <- ukmsoas[ukmsoas@data$geo_code %in% cents$geo_code, ]
```

The characteristics of zones are as follows:

```{r, echo=FALSE}
nzones <- nrow(zones) # how many zones?
zones_wgs <- spTransform(zones, CRS("+init=epsg:27700"))
mzarea <- round(median(gArea(zones_wgs, byid = T) / 10000), 1) # average area of zones, sq km
```

- Number of zones: `r nzones`, compared with 6791 in England
- Median area of zones: `r mzarea` ha, compared with 300 ha across England

## Input flow data

```{r, echo=FALSE, results='hide'}
# load flow dataset, depending on availability
if(!exists("flow_nat")){
if(dir.exists(pct_privatedata)){
  flow_nat <- readRDS(file.path(pct_privatedata, "flowsex-merged.Rds"))
  } else{
  flow_nat <- readRDS(file.path(pct_bigdata, "flow.Rds"))
  } # load open data
}

# Subset by zones in the study area
o <- flow_nat$Area.of.residence %in% cents$geo_code
d <- flow_nat$Area.of.workplace %in% cents$geo_code
flow <- flow_nat[o & d, ] # subset OD pairs with o and d in study area
n_flow_region <- nrow(flow)
n_commutes_region <- sum(flow$All)
```

```{r distance-dist, echo=FALSE, fig.cap="The study region (thick black border), selected zones (grey), the administrative zone region (red line) and local authorities (blue line). The black straight green represent the most intensive commuting OD pairs.", echo=FALSE, message=FALSE, warning=FALSE}
# Calculate line lengths (in km)
coord_from <- coordinates(cents[match(flow$Area.of.residence, cents$geo_code),])
coord_to <- coordinates(cents[match(flow$Area.of.workplace, cents$geo_code),])
# Euclidean distance (km)
flow$dist <- geosphere::distHaversine(coord_from, coord_to) / 1000
flow <- flow[flow$dist > 0,] # remove 0 dist flows

# Subset lines
# subset OD pairs by n. people using it
sel_long <- flow$All > mflow & flow$dist < mdist 
sel_short <- flow$dist < max_all_dist & flow$All > mflow_short
sel <- sel_long | sel_short
flow <- flow[sel, ]
# summary(flow$dist)
l <- od2line(flow = flow, zones = cents)

plot(zones, col = "lightgrey")
plot(regions, add = T)
plot(las_in_region, border = "blue", add = T, lwd = 2)
plot(region_orig, lwd = 5, add = T)
plot(region_shape, border = "red", add = T, lwd = 2)
lines(l[l$All > 100,], col = "green")
```

```{r, echo=FALSE}
# nrow(flow) # how many OD pairs in the study area?
# proportion of OD pairs in min-flow based subset
pmflow <- round(nrow(l) / n_flow_region * 100, 1)
# % all trips covered
pmflowa <- round(sum(l$All) / n_commutes_region * 100, 1)
```

There are **`r n_flow_region`** OD pairs with origins and destinations in the study
area. Of these, `r sum(sel_long)` meet the criteria that at least `r mflow` people travelled to work along OD pairs up to `r mdist` km in the 2011 Census. The additional selection criteria that at least `r mflow_short` people travelled to work along OD pairs up to `r max_all_dist` km was met by `r sum(sel_short)` OD pairs.
Adding those (overlapping) selection criteria resulted in
**`r nrow(flow)`** or **`r pmflow`%** of all inter-zone OD pairs were selected in the region, accounting for
**`r pmflowa`%** of inter-zone commutes in the study area.

## Hilliness of OD pairs

The average hilliness of zones in the study area is
`r round(100 * mean(rf$av_incline), 1)`
percent.

```{r, echo = FALSE}
# # It used to say - see below. What to replace this with?
# compared with the national average of
# `r round(mean(ukmsoas$avslope, na.rm = T), 1)`. This data is displayed in the
# figure below.
tm_shape(zones) +
  tm_fill("avslope", n = 3, palette = "Oranges")
```

```{r, echo=FALSE}
# # Hilliness of OD pairs - calculated the old way (depreciated)
# (calculated as the average gradient of the zone
# of the flow's origin and destination, in degrees)
# is
# `r round(mean(flow$avslope * flow$All / mean(flow$All), na.rm = T), 2)`.
# The UK
# average is xx degrees
```

## Lines allocated to the road network

We use CycleStreets.net to estimate optimal routes.
An illustration of these routes is presented below.

```{r flow-vars, echo=FALSE, message=FALSE}
# # # # # # # # # # # # # # # # # #
# Calculate flow-level variables: #
# distances and olc for ag. model #
# # # # # # # # # # # # # # # # # #

# Calculate distances (eventually use route distance)

# # # # # # # # # # # # # # #
# Allocate OD pairs2network #
# Warning: time-consuming!  #
# Needs CycleStreet.net API #
# # # # # # # # # # # # # # #
rf <- line2route(l, silent = TRUE)
rq <- line2route(l, plan = "quietest", silent = T)
rf$length <- rf$length / 1000 # set length correctly
rq$length <- rq$length / 1000

# debug lines which failed
if(!(nrow(l) == nrow(rf) & nrow(l) == nrow(rq))){
  # which paths succeeded
  path_ok <- row.names(l) %in% row.names(rf) &
                   row.names(l) %in% row.names(rq)
  # summary(path_ok)
  l <- l[path_ok,]
  path_ok <- row.names(rf) %in% row.names(l)
  rf <- rf[path_ok,]
  path_ok <- row.names(rq) %in% row.names(l)
  rq <- rq[path_ok,]
}
# saveRDS(rf, file.path(pct_data, region, "rf.Rds")) # save the routes
# saveRDS(rq, file.path(pct_data, region, "rq.Rds"))

# add line id
l$id <- paste(l@data$Area.of.residence, l@data$Area.of.workplace)
```

```{r plot-rlines, warning=FALSE, echo = FALSE, fig.cap="Straight and route-lines allocated to the travel network"}
plot(l)
lines(rf, col = "red")
lines(rq, col = "green")
```

```{r check-line, warning=FALSE, echo = FALSE, fig.cap="Check the final line"}
plot(l[nrow(l),])
lines(rf[nrow(l),], col = "red")
lines(rq[nrow(l),], col = "green")
```

Of the lines allocated to the route network, `r sum(!path_ok)` failed.

```{r, echo=FALSE}
# Allocate route characteristics to OD pairs
l$dist_fast <- rf$length
l$dist_quiet <- rq$length
l$time_fast <- rf$time
l$time_quiet <- rq$time
l$cirquity <- rf$length / l$dist
l$distq_f <- rq$length / rf$length
l$avslope <- rf$av_incline
l$co2_saving <- rf$co2_saving
l$calories <- rf$calories
l$busyness <- rf$busyness
l$avslope_q <- rq$av_incline
l$co2_saving_q <- rq$co2_saving
l$calories_q <- rq$calories
l$busyness_q <- rq$busyness
```

## Distance distributions

The distance distribution of trips in the study area is displayed in the figure below, which compares the result with the distribution of trips nationwide.

```{r, echo=FALSE, fig.cap="Distance distribution of all trips in study lines (blue) compared with national average (dotted bars)"}
luk <- readRDS(file.path(pct_bigdata, "l_sam8.Rds"))

hdfl <- dplyr::select(l@data, All, dist_fast)
hdfl$Scope <- "Local"
hdfl$All <- hdfl$All / sum(hdfl$All)

hdfu <- dplyr::select(luk@data, All, dist_fast)
hdfu$Scope <- "National"
hdfu$All <- hdfu$All / sum(hdfu$All)

histdf <- rbind(hdfl, hdfu)

ggplot(histdf) +
  geom_histogram(aes(dist_fast, weight = All, fill = Scope, linetype = Scope),
    position = "identity", colour = "black", binwidth = 0.5) +
  scale_fill_manual(values = c("lightblue", NA)) +
  scale_linetype(c(1, 2), guide = "none") +
  scale_y_continuous() +
  # scale_y_continuous(labels = percent) +
  xlab("Route distance (km)") +
  ylab("Proportion of trips in each band") +
  xlim(c(0,13)) +
  theme_bw()

pl5kmuk <- round(sum(luk$All[luk$dist_fast < 5]) /
    sum(luk$All) * 100, 1)
pl5km <- round(sum(l$All[l$dist_fast < 5]) /
    sum(l$All) * 100, 1)
```

From the nationwide sample of trips, `r pl5kmuk`% of trips are less than 5km.

In the case study area
`r pl5km`% of sampled trips are less than 5km.

Subsetting by distance (set
to `r mdist` km) and removing inter-zone OD pairs
further reduces the number of OD pairs from `r sum(sel)`
to `r nrow(l)`.

```{r, echo=FALSE}
# # # # # # # # # # # # # #
# Estimates slc from olc  #
# # # # # # # # # # # # # #

l$clc <- l$Bicycle / l$All
l$clcar <- (l$Car_driver + l$Car_passenger) / l$All
l$clfoot <- l$Foot / l$All   #% walk population
```

## The flow model

To estimate the potential rate of cycling under different scenarios
regression models operating at the flow level are used.
These can be seen in the model script which is available
[online](https://github.com/npct/pct/blob/master/models/aggregate-model.R).

```{r, echo=FALSE}
# Parameter creation for model
l$clc[l$clc < 0.001] <- 0.001 # for logit link to work
l$clc_logit <- log(l$clc / (1 - l$clc))

source("models.R")
mod_local <- lm(formula = logit_interact_formula, data = l@data)
logit_loc <- predict(object = mod_local, data = l@data)
l$local_slc <- exp(logit_loc) / (1 + exp(logit_loc))
mod_nat <- mod_local
mod_nat$coefficients <- c(-3.2941, -0.4659, 1.3430, 0.0062, -31.9249, 1.0989, -6.6162)
logit_nat <- predict(object = mod_nat, data = l@data)
l$slc <- exp(logit_nat) / (1 + exp(logit_nat))
cormod <- cor(l$clc, l$slc) # crude indication of goodness-of-fit
```

## Cycling in the study area

```{r, echo=FALSE}
rcycle <- round(100 * sum(l$Bicycle) / sum(l$All), 1)
rcarusers <- round (100 * sum(l$Car_driver+l$Car_passenger) / sum(l$All), 1)
natcyc <- sum(luk$Bicycle) / sum(luk$All)
```

The overall rate of cycling in the OD pairs in the study area
(after subsetting for distance) is `r rcycle`%, compared a
rate from the national data (of equally short OD pairs)
of 5.0%.

## Scenarios

```{r, include=FALSE}
l$base_slc <- l$slc * l$All

# # # # # # # # # # # # #
# Additional scenarios  #
# # # # # # # # # # # # #

# Replace with source("models/aggregate-model-dutch|gendereq|ebike.R"))
l$govtarget_slc <- l$Bicycle + l$base_slc # may be more that 1
l$govtarget_sic <- l$govtarget_slc - l$Bicycle

# gendereq scenario
has_private_data <- file.exists(pct_privatedata)
if(has_private_data){
  p_trips_male <- mean(las_in_region$clc_m) # proportion of bicycle trips by males
  
  clc_m <- l$Bicycle * p_trips_male
  pmale_c <- clc_m / l$Male
  slc_gendereq_f <- l$Female * pmale_c
  slc_gendereq <- clc_m + slc_gendereq_f

  l$gendereq_slc <- slc_gendereq
  l$gendereq_sic <- l$gendereq_slc - l$Bicycle
}else{
  l$gendereq_slc <- NA
  l$gendereq_sic <- NA
}
# Dutch scenario - coefficients calculated from Dutch NTS by A. Goodman
logit_dutch <- logit_nat + 2.4783 -0.0706 * l$dist_fast
l$dutch_slc <- exp(logit_dutch) / (1 + exp(logit_dutch)) * l$All
l$dutch_sic <- l$dutch_slc - l$Bicycle

# Ebike scenario
logit_ebike <- logit_dutch + 0.0286 + 0.0902 * l$dist_fast - 0.0014 * l$dist_fast^2

l$ebike_slc <- exp(logit_ebike) / (1 + exp(logit_ebike)) * l$All
l$ebike_sic <- l$ebike_slc - l$Bicycle

dfscen <- dplyr::select(l@data, contains("slc"), All, olc = Bicycle, dist_fast)
dfscen <- dfscen[-which(names(dfscen) == "slc")]
dfscen <- dfscen[-which(names(dfscen) == "base_slc")]
# head(dfscen)

dfsp <- gather(dfscen, key = scenario, value = slc, -dist_fast)
# head(dfsp)
dfsp$scenario <- factor(dfsp$scenario)
 summary(dfsp$scenario)
dfsp$scenario <-
  factor(dfsp$scenario, levels = levels(dfsp$scenario)[c(5, 4, 1, 2, 3, 6)])
scalenum <- sum(l$All)
```

```{r, echo=FALSE, warning=FALSE, fig.cap="Rate of cycling in model scenarios. Note the total percentage cycling is equal to the area under each line."}
ggplot(dfsp) +
  geom_freqpoly(aes(dist_fast, weight = slc,
    color = scenario), binwidth = 1) +
  ylab("Total number of trips") +
  xlab("Route distance (km)") +
  scale_color_discrete(name = "Mode and\nscenario\n(cycling)") +
  xlim(c(0,12)) +
  theme_bw()

dfsp$dist_band <- cut(dfsp$dist_fast, c(0, 2, 5, 10, 20))
dfsum <- summarise(group_by(dfsp, scenario, dist_band), Percent = sum(slc) / sum(l$All))
dfsum$Percent <- dfsum$Percent
dfspread <- spread(dfsum, scenario, Percent)
dfspread$dist_band <- as.character(dfspread$dist_band)
dfspreadf <- c("Total", round(colSums(dfspread[2:ncol(dfspread)])* 100, 1))
dfspread[3:ncol(dfspread)] <- do.call(cbind, apply(dfspread[3:ncol(dfspread)], 2, function(x) round(x / dfspread[2] * 100, 1)))
dfspread <- rbind(dfspread, dfspreadf)
# dfspread <- dfspread[c(1, 2, 7, 3, 4, 5, 6)]
dfspread$All <- round(as.numeric(dfspread$All) * 100, 1)
dfspread$All[nrow(dfspread)] <- dfspread$All[nrow(dfspread)] / 100
```

The table below illustrates the same information by distance band.

```{r, echo=FALSE}
# names(dfspread)[1:3] <- c("Distance band", 	"All modes", "Observed level (OLC)")
kable(dfspread, format = "html", digits = 1)
```


```{r, include=FALSE}
# # # # # # # # # # # # # # # # # #
# Extract area-level commute data #
# # # # # # # # # # # # # # # # # #

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

  # all OD pairs originating from centroid i
  j <- which(l$Area.of.residence == cents$geo_code[i])

  cents$Bicycle[i] <- sum(l$Bicycle[j])
  cents$base_slc[i] <- sum(l$base_slc[j])
  cents$base_sic[i] <- sum(l$base_sic[j])
  
  # car users % added to centroids
  cents$All[i] <- sum(l$All[j])
  cents$Car[i] <- sum(l$Car_driver[j])
  cents$base_olcarusers[i] <- cents$Car[i] / cents$All[i]

  # values for scenarios
  cents$govtarget_slc[i] <- sum(l$govtarget_slc[j])
  cents$govtarget_sic[i] <- sum(l$govtarget_sic[j])

  cents$gendereq_slc[i] <- sum(l$gendereq_slc[j])
  cents$gendereq_sic[i] <- sum(l$gendereq_sic[j])

  cents$dutch_slc[i] <- sum(l$dutch_slc[j])
  cents$dutch_sic[i] <- sum(l$dutch_sic[j])

  cents$ebike_slc[i] <- sum(l$ebike_slc[j])
  cents$ebike_sic[i] <- sum(l$ebike_sic[j])

  cents$av_distance[i] <- sum(l$dist[j] * l$All[j])  / sum(l$All[j])
  cents$cirquity[i] <- sum(l$cirquity[j] * l$All[j], na.rm = T )  / sum(l$All[j])
  cents$distq_f[i] <- sum(l$distq_f[j] * l$All[j], na.rm = T )  / sum(l$All[j])
}

l$Rail <- l$Light_rail + l$Train
l$Other <- l$Other + l$Taxi + l$Motorbike
# Remove aggregated categories
l$Taxi <- l$Light_rail <- l$Train <- l$Motorbike <- NULL 

addidsn <- c("All", "From_home", "Bus",
"Rail", "Taxi", "Motorbike", "Car_driver", "Car_passenger", "Bicycle",
"Foot", "Other", "base_slc",
"base_sic", "gendereq_slc", "gendereq_sic", "dutch_slc", "dutch_sic",
"ebike_slc", "ebike_sic", "govtarget_slc", "govtarget_sic")
addids <- which(names(l) %in% addidsn)
# addids <- c(3:14, 23:31)
# summary(l[addids])


# Aggregate bi-directional OD pairs

# Subset by zone bounding box
# l <- l[as.logical(gContains(zone, l, byid = T)),]
# nrow(l)

# 4: by aggregating 2 way OD pairs
l <- onewayid(l, attrib = addids)

# remove origin and destinations now they are redundant
l$Area.of.residence <- l$Area.of.workplace <- NULL
l$clc <- l$Bicycle / l$All
l$slc <- l$base_slc / l$All

nrow(l)
idsel <- row.names(l)
rf <- rf[row.names(rf) %in% idsel,]
rq <- rq[row.names(rq) %in% idsel,]

# Sanity test
head(l@data[1:5])
```

```{r, echo=FALSE, results='hide', fig.cap="Illustration of OD pairs on travel network"}
plot(region_shape)
plot(zones, add = T)
points(cents, col = "red")
lines(l, col = "black")
lines(rq, col = "green")
lines(rf, col = "blue")
```

## Flow model results

To estimate the potential rate of cycling under different scenarios
regression models operating at the flow level are used.
These can be seen in the model script which is available
[online](https://github.com/npct/pct/blob/master/models/aggregate-model.R).

```{r, echo=FALSE, fig.cap="National vs local cycling characteristics with hilliness, captured in the model results"}
justdist1 <- data.frame(
  dist_fast = 1:20,
  avslope = 1,
  type = "Flat")

justdist2 <- justdist1
justdist2$avslope <- 1.5
justdist2$type <- "Hilly"

justdist <- rbind(justdist1, justdist2) # for prediction
justdist$model <- "National"

justdist5 <- justdist6 <- justdist7 <- justdist # replicate
justdist5$model <- "Local"
justdist6$model <- "Dutch"
justdist7$model <- "Ebike"

justdist$npred <- exp(predict(mod_nat, justdist))
justdist5$npred <- NA # todo: fix this
justdist6$npred <- exp(predict(mod_dutch, justdist))
justdist7$npred <- exp(predict(mod_ebike, justdist))

justdist <- rbind(justdist, justdist5, justdist6, justdist7)

ggplot(justdist) +
  geom_line(aes(dist_fast, npred, color = model, linetype = type),
    size = 1.5) +
  xlab("Route distance (km)") + ylab("Expected proportion cycling") +
  theme_bw()

dfcos <- round(rbind(coef(mod_nat)), 3)
dfcos <- cbind(Model = c("National"), dfcos)
```

The correlation between fitted and observed cycling in the model is
`r round(cormod, 2)`, compared with 0.39 nationally.

The values for the coefficients are presented in the table below.

```{r, echo=FALSE}
dfcos <- data.frame(dfcos)
names(dfcos) <- c("Model", "Alpha", "Distance", "Dist^0.5", "Hills", "Dist/Hills")
kable(dfcos, digits = 3)
```

## Network analysis

Now we aggregate the overlapping routes to create a route network.
The value of each segment in the network corresponds to the total number of cyclists who we estimate to use the segment.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="The route network, with widths proportional to the current estimated number of commuter cyclists"}
# Scenario names - may need updating
# nrow(l) == nrow(rf)

rft <- rf
# Stop rnet lines going to centroid (optional)
rft <- toptail(rf, toptail_dist = buff_geo_dist)
if(length(rft) == length(rf)){
  row.names(rft) <- row.names(rf)
  rft <- SpatialLinesDataFrame(rft, rf@data)
} else print("Error: toptailed lines do not match lines")
rft$Bicycle <- l$Bicycle

# Simplify line geometries (if mapshaper is available)
# this greatly speeds up the build (due to calls to overline)
# needs mapshaper installed and available to system():
# see https://github.com/mbloch/mapshaper/wiki/
rft <- ms_simplify(rft, keep = 0.1)
rnet <- overline(rft, "Bicycle")
# object.size(rnet)
# test the resulting plot
plot(rnet, lwd = rnet$Bicycle / mean(rnet$Bicycle))
for(i in scens){
  rft@data[i] <- l@data[i]
  rnet_tmp <- overline(rft, i)
  rnet@data[i] <- rnet_tmp@data[i]
  rft@data[i] <- NULL
}
if(!"gendereq_slc" %in% scens)
  rnet$gendereq_slc <- NA
# plot(rnet, lwd = rnet$govtarget_slc / mean(rnet$Bicycle))
```


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}
# # # # # # # # #
# Save the data #
# # # # # # # # #

# Remove/change private/superfluous variables
l$Male <- l$Female <- l$From_home <- l$calories <-
  l$co2_saving_q <-l$calories_q <- l$busyness_q <- 
  # data used in the model - superflous for pct-shiny
  l$dist_fastsq <- l$dist_fastsqrt <- l$ned_avslope <- 
  l$interact <- l$interactsq <- l$interactsqrt <- NULL

# Transfer cents data to zones
c_in_z <- names(cents) == "avslope"
zones@data <- left_join(zones@data, cents@data[,!c_in_z])

# Remove NAN numbers (cause issues with geojson_write)
na_cols  <- which(names(zones) %in%
  c("av_distance", "cirquity", "distq_f", "base_olcarusers", "gendereq_slc", "gendereq_sic"))
for(ii in na_cols){
  zones@data[[ii]][is.nan(zones@data[[ii]])] <- NA
}

# # Save objects
# Save objects # uncomment these lines to save model output
if(exists("isolated")) file.create(file.path(pct_data, region, "isolated"))
saveRDS(zones, file.path(pct_data, region, "z.Rds"))
geojson_write( ms_simplify(zones, keep = 0.1), file = file.path(pct_data, region, "z"))
saveRDS(cents, file.path(pct_data, region, "c.Rds"))
saveRDS(l, file.path(pct_data, region, "l.Rds"))
saveRDS(rf, file.path(pct_data, region, "rf.Rds"))
saveRDS(rq, file.path(pct_data, region, "rq.Rds"))
saveRDS(rnet, file.path(pct_data, region, "rnet.Rds"))
saveRDS(mod_logsqr, file.path(pct_data, region, "model.Rds"))

# # Save the script that loaded the lines into the data directory
file.copy("load.Rmd", file.path(pct_data, region, "load.Rmd"))

# Create folder in shiny app folder
region_dir <- file.path(file.path(pct_shiny_regions, region))
dir.create(region_dir)
ui_text <- 'source("../../ui-base.R", local = T, chdir = T)'
server_text <- paste0('startingCity <- "', region, '"', "\n",
                      'source("../../server-base.R", local = T, chdir = T)')
write(ui_text, file = file.path(region_dir, "ui.R"))
write(server_text, file = file.path(region_dir, "server.R"))
file.symlink(file.path("..", "..","www"), region_dir)
```

## Time taken

The time taken to run the analysis for this area is presented below:

```{r}
end_time <- Sys.time()

end_time - start_time
```