load.Rmd

---
title: "National Propensity to Cycle Tool - local results"
author: "Created by the NPCT team"
output:
  word_document: 
    fig_caption: yes
  html_document:
    fig_caption: yes
---

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

This document was produced automatically at `r start_time`.

## Introduction

The results of NPCT scenarios are based on a model.
This document presents information about the input data, model diagnostics,
run time and information about its outputs.

The aim is to provide further information for transport planners and
researchers on the origin of results presented in the NPCT's interactive map.
This document is designed for use by advanced users: some technical knowledge
is needed to understand all its outputs.

The model driving the NPCT is licensed under the open source MIT License
and can be modified by others provided attribution to the original.

## Initial parameters

The preset values used to select the study area and filter the flow data
were as follows:

```{r, echo=FALSE}
if(!exists("la"))
  la <- "manchester" # create default LA name if none exists
```


```{r, warning=FALSE}
# Set local authority and ttwa zone names
la # print the name of the local authority
la_path <- file.path("pct-data", la)

if(!dir.exists(la_path)) dir.create(la_path) # on a unix machine

# Minimum flow between od pairs, subsetting lines. High means fewer lines.
mflow <- 30
mdist <- 15 # maximum euclidean distance (km) for subsetting lines
min_zones <- 60 # minumum number of zones in study area before buffer used
buff_dist <- 3 # buffer (km) used to select additional zones (often redundant)
```

## Input zone data

The input zones area are summarised in this section.

```{r plotzones, message=FALSE, warning=FALSE, results='hide', echo=FALSE}
ukmsoas <- shapefile("pct-bigdata/national/msoas.shp")

# Load population-weighted centroids
cents <- readOGR("pct-bigdata/national/cents.geojson", layer = "OGRGeoJSON")
cents$geo_code <- as.character(cents$geo_code)

# Load local authorities and districts
las <- readOGR(dsn = "pct-bigdata/national/las-pcycle.geojson", layer = "OGRGeoJSON")
sel_la <- grepl(pattern = la, x = las$NAME, ignore.case = T)
```


```{r, include=FALSE}
# lasdat <- SpatialPointsDataFrame(coords = coordinates(las), data = las@data)

# Load counties and unitary authorities
# cuas <- readOGR(dsn = "pct-bigdata/national/cuas.geojson", layer = "OGRGeoJSON")
# proj4string(lasdat) <- proj4string(cuas)
# cuas <- aggregate(lasdat, cuas, mean, na.action = na.omit()) # todo: fix data
# cua_shape <- cuas[grep(pattern = la, x = cuas@data$NAME)] # todo: fix
# # tmap::qtm(cuas2,fill = "clc")

la_shape <- las[grep(pattern = la, x = las@data$NAME, ignore.case = T),]

# select the cents of interest: by char string (above) or geography (below)
if(grepl(pattern = "cov", x = la) > 0){
  cents_la <- cents[grepl(la, cents$MSOA11NM, ignore.case = T), ]
} else{
  cents_la <- cents[la_shape,]
}

zones <- ukmsoas[ukmsoas@data$geo_code %in% cents_la@data$geo_code, ]
zbuf <- gBuffer(zones, width = buff_dist * 1000) # create la zone outline

# convert to osbg
cents <- spTransform(cents, CRSobj = proj4string(ukmsoas))
la_shape <- spTransform(la_shape, CRSobj = proj4string(ukmsoas))

# Extract zones to plot
# zones <- ukmsoas[ grep(la, ukmsoas$geo_label), ] # extract by name

if(nrow(zones) < min_zones){
  cents_la <- cents[zbuf, ]
  zones <- ukmsoas[ukmsoas@data$geo_code %in% cents_la@data$geo_code , ]
}

cents <- cents_la
```

The characteristics of zones are as follows:

```{r, echo=FALSE}
nzones <- nrow(zones) # how many zones?
mzarea <- round(median(gArea(zones, byid = T) / 10000), 2) # 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'}
flow <- readRDS("pct-bigdata/national/flow_eng_avlslope.Rds")

# Load flow-sex data. This is a local file
if(file.exists("private-data/wu01bew_msoa_v1.csv")){
  flowsex <- read_csv("private-data/wu01bew_msoa_v1.csv", col_names=F, skip=12)
} else{
  print("Warning: data for gender equality scenario missing")
  print("See loading-data/load-flowsex.R")
}

names(flowsex) <- c("Area.of.residence", "Area.of.workplace", "All.fs", "Male", "Female")
flowsex$id <- paste(flowsex$Area.of.residence, flowsex$Area.of.workplace)

# Subset by zones in the study area
o <- flow$Area.of.residence %in% cents$geo_code
d <- flow$Area.of.workplace %in% cents$geo_code
flow <- flow[o & d, ] # subset flows with o and d in study area
sel <- flow$All > mflow # subset flows by n. people using it
```

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

There are **`r nrow(flow)`** flows with origins and destinations in the study
area. Of these, **`r sum(sel)`** meet the criteria that `r mflow` people 
reported travelling between these two areas in the 2011 Census (a modifiable parameter used
to reduce the computational requirements
of the model). That is  **`r pmflow`%** of flows, accounting for
**`r pmflowa`%** of commutes in the study area.

```{r distance-dist, echo=FALSE, fig.cap="The study area (dark line), selected zones for training the model (grey) and the administrative zone of interest (red line). The black straight lines represent the most intensive commuting flows.", echo=FALSE, message=FALSE, warning=FALSE}
flow <- flow[sel, ]
# nrow(flow) # new flow rate

flow$id <- paste(flow$Area.of.residence, flow$Area.of.workplace)
flowsex <- dplyr::select(flowsex, Male, Female, id)
flow <- left_join(flow, flowsex, by = "id")

l <- gFlow2line(flow = flow, zones = cents)
plot(zones, col = "lightgrey")
plot(las, add = T)
plot(zbuf, lwd = 5, add = T)
plot(la_shape, border = "red", add = T, lwd = 3)
lines(l[l$All > 100,])
```

```{r, echo=FALSE}
proj4string(l) <- proj4string(cents)
l$dist <- gLength(l, byid = T) / 1000 # Euclidean distance
dsel <- l$dist < mdist

l <- l[dsel,]
l <- l[l$dist > 0, ] # to remove flows of 0 length
```

## Hilliness of flows

The average hilliness of zones in the study area is
`r round(mean(zones$avslope), 2)`
degrees,
compared with the national average of 
`r round(mean(ukmsoas$avslope, na.rm = T), 2)`. This data is displayed in the
figure below.

```{r, echo = FALSE}
tm_shape(zones) +
  tm_fill("avslope", n = 3, palette = "Oranges")
```

```{r, echo=FALSE}
# Hilliness of flows
# (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}
# # # # # # # # # # # # # # # # # #
# Calculate flow-level variables: #
# distances and olc for ag. model #
# # # # # # # # # # # # # # # # # #

# Calculate distances (eventually use route distance)


# Transform CRS to WGS84 for plotting
zones <- spTransform(zones, CRS("+init=epsg:4326"))
cents <- spTransform(cents, CRS("+init=epsg:4326"))
l <- spTransform(l, CRS("+init=epsg:4326"))

# # # # # # # # # # # # # # #
# Allocate flows to network #
# Warning: time-consuming!  #
# Needs CycleStreet.net API #
# # # # # # # # # # # # # # #

# Create local version of lines; if there are too many in the TTWA, sample!
l_local_sel <- l@data$Area.of.residence %in% zones$geo_code &
  l@data$Area.of.workplace %in% zones$geo_code
if(nrow(l) > 2 * sum(l_local_sel) & nrow(l) > 5000){ # sample if too many lines
  l_all <- l
  f <- list.files(paste0("pct-data/", la, "/"))
  if(sum(grepl("l_all", f)) == 0) saveRDS(l, paste0("pct-data/", la, "/l_all.Rds"))
#   l <- readRDS(paste0("pct-data/", la, "/l_all.Rds")) # restart point
  set.seed(2050)
  # sample from all routes in the TTWZ - change 1 for different % outside zone
  lsel <- sample(which(!l_local_sel), size = sum(l_local_sel) * 1)
  lsel <- c(lsel, which(l_local_sel))
  length(lsel)
  l <- l_all[lsel, ] # subset the lines
  # plot(l)
  # lines(l[2000:2600,], col = "blue") # ensure we have all the local ones
}

# Create route allocated lines
if(length(grep("rf.Rds|rq.Rds", list.files(paste0("pct-data/", la)))) >= 2){
  rf <- readRDS(paste0("pct-data/", la, "/rf_ttwa.Rds")) # if you've loaded them
  rq <- readRDS(paste0("pct-data/", la, "/rq_ttwa.Rds"))
  # l <- readRDS(paste0("pct-data/", la, "/l.Rds"))
  nrow(rf) == nrow(l)
} else{
  rf <- gLines2CyclePath(l[ l$dist > 0, ])
  rq <- gLines2CyclePath(l[ l$dist > 0, ], plan = "quietest")

  # Process route data
  rf$length <- rf$length / 1000
  rq$length <- rq$length / 1000
  rf$id <- l$id
  rq$id <- l$id
  saveRDS(rf, paste0("pct-data/", la, "/rf_ttwa.Rds")) # save the routes
  saveRDS(rq, paste0("pct-data/", la, "/rq_ttwa.Rds"))
  saveRDS(l, paste0("pct-data/", la, "/l.Rds"))
  }
```

```{r plot-rlines, echo = FALSE, fig.cap="Sample of the straight and route-lines allocated to the travel network"}
library(sp)
plot(l[1:300,])
lines(rf[1:300,], col = "red")
lines(rq[1:300,], col = "green")
```

```{r, echo=FALSE}

# Allocate route factors to flows
# nz <- which(l$dist > 0) # non-zero lengths = nz
l$dist_quiet <- l$dist_fast <- l$cirquity <- l$distq_f <- NA
if(nrow(rf) == nrow(l)) print("Warning, lines and routes are different lengths")
  l$dist_fast <- rf$length
  l$dist_quiet <- rq$length
  l$cirquity <- rf$length / l$dist
  l$distq_f <- rq$length / rf$length
```

## 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("pct-bigdata/national/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, 2)
pl5km <- round(sum(l$All[l$dist_fast < 5]) /
    sum(l$All) * 100, 2)
```

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 flows
further reduces the number of flows from `r sum(sel)`
to `r nrow(l)`.

```{r, echo=FALSE}
# a = 11
# plot(l[a,])
# lines(rf[a,], col = "red")
# lines(rq[a,], col = "green")

# # # # # # # # # # # # # #
# Estimates slc from olc  #
# # # # # # # # # # # # # #

l$clc <- l$Bicycle / l$All
flow_ttwa <- flow # save flows for the ttwa
flow <- l@data
```

## 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}
source("models/aggregate-model.R") # this model creates the variable 'slc'
cormod <- cor(flow$clc, mod_logsqr$fitted.values) # crude indication of goodness-of-fit
# summary(mod_logsqr)

mod_nat <- readRDS("pct-bigdata/national/mod_logsqr_national_8.Rds")
```

## Cycling in the study area

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

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

## Scenarios

```{r, include=FALSE}
l$slc <- flow$plc
l$base_olc <- l$Bicycle
l$base_slc <- l$slc * l$All
l$base_sic <- l$base_slc - l$base_olc
# l$sic2 <- l$slc * l$All - l$Bicycle # identical sic result

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

# Additional scenarios
# Replace with source("models/aggregate-model-dutch|gendereq|ebike.R"))
set.seed(2015)
l$npred <- exp(predict(mod_nat, flow))

l$cdp_slc <-l$All * (l$clc + l$npred) # may be more that 1
l$cdp_sic <- l$cdp_slc - l$Bicycle

# gendereq scenario
area_pcycle <- las@data[sel_la,]
p_trips_male <- area_pcycle$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$base_olc

# Dutch scenario - coefficients calculated from Dutch NTS by A. Goodman
mod_dutch <- mod_nat
mod_dutch$coefficients[1] <- -0.3253
mod_dutch$coefficients[2] <- -0.3543
mod_dutch$coefficients[3] <- 0.7750
l$dutch_slc <- l$All * exp(predict(mod_dutch, flow))
l$dutch_sic <- l$dutch_slc - l$base_olc

mod_ebike <- mod_dutch
mod_ebike$coefficients[2] <- -0.3

l$ebike_slc <- l$All * exp(predict(mod_ebike, flow))
l$ebike_sic <- l$ebike_slc - l$base_olc

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:7])* 100, 2))
dfspread[3:7] <- do.call(cbind, apply(dfspread[3:7], 2, function(x) round(x / dfspread[2] * 100, 2)))
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
```


```{r, echo=FALSE}
kable(dfspread)
```


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

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

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

  cents$base_olc[i] <- sum(l$Bicycle[j])
  cents$base_slc[i] <- sum(l$base_slc[j])
  cents$base_sic[i] <- sum(l$base_sic[j])

  # values for scenarios
  cents$cdp_slc[i] <- sum(l$cdp_slc[j])
  cents$cdp_sic[i] <- sum(l$cdp_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])
}

# names(l) # which line names can be added for non-directional flows?
# dput(c(names(l)[addids], c("cdp_slc", "cdp_sic")))
addidsn <- c("All", "Work.mainly.at.or.from.home", "Underground..metro..light.rail..tram", 
"Train", "Bus..minibus.or.coach", "Taxi", "Motorcycle..scooter.or.moped", 
"Driving.a.car.or.van", "Passenger.in.a.car.or.van", "Bicycle", 
"On.foot", "Other.method.of.travel.to.work", "base_olc", "base_slc", 
"base_sic", "gendereq_slc", "gendereq_sic", "dutch_slc", "dutch_sic", 
"ebike_slc", "ebike_sic", "cdp_slc", "cdp_sic")
addids <- which(names(l) %in% addidsn)
# addids <- c(3:14, 23:31)
# summary(l[addids])


# Aggregate bi-directional flows

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

# 4: by aggregating 2 way flows
l <- gOnewayid(l, attrib = c(addids))

l$clc <- l$Bicycle / l$All
l$slc <- l$base_slc / l$All

nrow(l)
idsel <- l$id
rf <- rf[rf@data$id %in% idsel,]
rq <- rq[rq@data$id %in% idsel,]

# Sanity test
head(l@data[1:5])
cents_ttwa <- cents # copy cents data (we'll overwrite cents)

# # Subset to zone
# cents <- cents_ttwa[zone,] # subset centroids geographically
# zones <- zones[cents,]
```

```{r, echo=FALSE, results='hide', fig.cap="Illustration of flows on travel network"}
zbuf <- spTransform(zbuf, CRS("+init=epsg:4326"))
plot(zbuf)
plot(zones, add = T)
points(cents_ttwa, 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"}
source("models/aggregate-model.R") # this model creates the variable 'slc'
cormod <- cor(flow$clc, mod_logsqr$fitted.values) # crude indication of goodness-of-fit
# summary(mod_logsqr)

mod_nat <- readRDS("pct-bigdata/national/mod_logsqr_national_8.Rds")

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 <- exp(predict(mod_logsqr, justdist))
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), coef(mod_logsqr)), 3)
dfcos <- cbind(Model = c("National", "Local"), 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}
kable(dfcos, digits = 3)
```

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

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

# 
# # Save objects
# Save objects # uncomment these lines to save model output
saveRDS(zones, paste0("pct-data/", la, "/z.Rds"))
saveRDS(cents, paste0("pct-data/", la, "/c.Rds"))
saveRDS(l, paste0("pct-data/", la, "/l.Rds"))
saveRDS(rf, paste0("pct-data/", la, "/rf.Rds"))
saveRDS(rq, paste0("pct-data/", la, "/rq.Rds"))
saveRDS(mod_logsqr, paste0("pct-data/", la, "/model.Rds"))
# 
# # # Save data for wider ttwz area
# # saveRDS(ttwa_zone, paste0("pct-data/", la, "/ttw_zone.Rds"))
# # saveRDS(cents_ttwa, paste0("pct-data/", la, "/c_ttwa.Rds"))
# # saveRDS(l_ttwa, paste0("pct-data/", la, "/l_ttwa.Rds"))
# 
# # Create new folder in pct-shiny repo
# rname <- tolower(la)
# dname <- paste0("~/repos/pct-shiny/", rname, "/")
# dir.create(dname)
# files <- list.files("~/repos/pct-shiny/manchester/", full.names = T)
# file.copy(files, dname)
# server <- readLines(paste0(dname, "server.R"))
# server <- gsub("manchester", la, server)
# writeLines(server, paste0(dname, "server.R"))
# 
# # Save the script that loaded the lines into the data directory
file.copy("load.Rmd", paste0("pct-data/", la, "/load.Rmd"))
```

## Time taken

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

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

end_time - start_time
```