stplanr-paper.R
In stplanr: Sustainable Transport Planning

## ---- echo=FALSE--------------------------------------------------------------
knitr::opts_chunk$set(fig.width = 7, fig.height = 5, eval = FALSE)

## ---- eval=FALSE--------------------------------------------------------------
#  install.packages("stplanr")

## -----------------------------------------------------------------------------
#  library(stplanr)

## ---- echo=FALSE, results='asis', message=FALSE-------------------------------
#  # stplanr_funs = ls("package:stplanr")
#  # sel_core = grep(pattern = "od_|^line_|route_", x = stplanr_funs)
#  # core_funs = stplanr_funs[sel_core]
#  # args(name = core_funs[1])
#  fun_table <- read.csv("fun_table.csv", stringsAsFactors = FALSE, check.names = FALSE)
#  knitr::kable(fun_table, caption = "Selection of functions for working with or generating OD, line and route data types.")

## ---- eval=FALSE--------------------------------------------------------------
#  dl_stats19() # download and extract stats19 road traffic casualty data

## ---- eval=FALSE--------------------------------------------------------------
#  ac <- read_stats19_ac()
#  ca <- read_stats19_ca()
#  ve <- read_stats19_ve()

## ---- eval=FALSE--------------------------------------------------------------
#  library(dplyr)
#  ca_ac <- inner_join(ca, ac)
#  ca_cycle <- ca_ac %>%
#    filter(Casualty_Severity == "Fatal" & !is.na(Latitude)) %>%
#    select(Age = Age_of_Casualty, Mode = Casualty_Type, Longitude, Latitude)
#  ca_sp <- SpatialPointsDataFrame(coords = ca_cycle[3:4], data = ca_cycle[1:2])

## ---- eval=FALSE--------------------------------------------------------------
#  data("route_network") # devtools::install_github("ropensci/splanr")version 0.1.7
#  proj4string(ca_sp) <- proj4string(route_network)
#  bb <- bb2poly(route_network)
#  proj4string(bb) <- proj4string(route_network)
#  ca_local <- ca_sp[bb, ]

## ---- echo=FALSE--------------------------------------------------------------
#  bb <- bb2poly(route_network)
#  load("reqfiles.RData")

## ---- message=FALSE-----------------------------------------------------------
#  rnet_buff_100 <- geo_buffer(route_network, width = 100)
#  ca_buff <- ca_local[rnet_buff_100, ]

## ----fats, fig.cap="Road traffic fatalities in the study area downloaded with with stplanr (crosses). Deaths that happened within 100 m of the route network are represented by circles.", out.width="50%", fig.align="center"----
#  plot(bb, lty = 4)
#  plot(rnet_buff_100, col = "grey", add = TRUE)
#  points(ca_local, pch = 4)
#  points(ca_buff, cex = 3)

## -----------------------------------------------------------------------------
#  data("flow", package = "stplanr")
#  head(flow[c(1:3, 12)])

## -----------------------------------------------------------------------------
#  data("cents", package = "stplanr")
#  as.data.frame(cents[1:3, -c(3, 4)])

## ---- warning=FALSE-----------------------------------------------------------
#  l <- od2line(flow = flow, zones = cents)

## ---- eval=FALSE--------------------------------------------------------------
#  route_bl <- route_cyclestreets(from = "Bradford", to = "Leeds")
#  route_c1_c2 <- route_cyclestreets(cents[1, ], cents[2, ])

## ---- eval=FALSE--------------------------------------------------------------
#  route_bl_raw <- route_cyclestreets(from = "Bradford", to = "Leeds", save_raw = TRUE)

## ----lines_routes, out.width='50%', fig.cap='Visualisation of travel desire lines, with width proportional to number of trips between origin and destination (black) and routes allocated to network  (red) in the left-hand panel. The right hand panel shows the route network dataset generated by overline().', fig.show='hold'----
#  plot(route_network, lwd = 0)
#  plot(l, lwd = l$All / 10, add = TRUE)
#  lines(routes_fast, col = "red")
#  routes_fast$All <- l$All
#  rnet <- overline(routes_fast, "All", fun = sum)
#  rnet$flow <- rnet$All / mean(rnet$All) * 3
#  plot(rnet, lwd = rnet$flow / mean(rnet$flow))

## ---- eval=FALSE, out.width='\\textwidth'-------------------------------------
#  ny2oaxaca1 <- route_graphhopper("New York", "Oaxaca", vehicle = "bike")
#  ny2oaxaca2 <- route_graphhopper("New York", "Oaxaca", vehicle = "car")
#  rbind(ny2oaxaca1@data, ny2oaxaca2@data)

## ---- eval=FALSE, echo=FALSE--------------------------------------------------
#  nytab <- rbind(ny2oaxaca1@data, ny2oaxaca2@data)
#  nytab <- cbind(Mode = c("Cycle", "Car"), nytab)
#  xtnyoa <- xtable(nytab, caption = "Attribute data from the route\\_graphhopper function, from New York to Oaxaca, by cycle and car.", label = "tab:xtnyoa")
#  print.xtable(xtnyoa, include.rownames = FALSE)
#  plot(ny2oaxaca1)
#  plot(ny2oaxaca2, add = TRUE, col = "red")
#  
#  ny2oaxaca1@data
#  ny2oaxaca2@data

## ----loadshapefiles, results='hide',message='hide'----------------------------
#  data_dir <- system.file("extdata", package = "stplanr")
#  unzip(file.path(data_dir, "smallsa1.zip"))
#  unzip(file.path(data_dir, "testcycleway.zip"))
#  sa1income <- as(sf::read_sf("smallsa1.shp"), "Spatial")
#  testcycleway <- as(sf::read_sf("testcycleway.shp"), "Spatial")
#  # Remove unzipped files
#  file.remove(list.files(pattern = "^(smallsa1|testcycleway).*"))

## ----calccatchment, results='hide', eval=FALSE--------------------------------
#  remotes::install_github("ropensci/stplanr")
#  catch800m <- calc_catchment(
#    polygonlayer = sa1income,
#    targetlayer = testcycleway,
#    calccols = c("Total"),
#    distance = 800,
#    projection = "austalbers",
#    dissolve = TRUE
#  )

## ----catchmentplot, fig.cap='An 800 metre catchment area (red) associated with a cycle path (green) using straight-line distance in Sydney.'----
#  plot(sa1income, col = "light grey")
#  plot(catch800m, col = rgb(1, 0, 0, 0.5), add = TRUE)
#  plot(testcycleway, col = "green", add = TRUE)

## ---- echo=TRUE, message=FALSE, warning=FALSE, results='hide'-----------------
#  unzip(file.path(data_dir, "sydroads.zip"))
#  sydroads <- as(sf::read_sf(".", "roads"), "Spatial")
#  file.remove(list.files(pattern = "^(roads).*"))
#  sydnetwork <- SpatialLinesNetwork(sydroads)

## ---- warning=FALSE-----------------------------------------------------------
#  netcatch800m <- calc_network_catchment(
#    sln = sydnetwork,
#    polygonlayer = sa1income,
#    targetlayer = testcycleway,
#    calccols = c("Total"),
#    maximpedance = 800,
#    distance = 100,
#    projection = "austalbers"
#  )

## ----netcatchplot, fig.cap='A 800 metre network catchment are (blue) compared with a catchment area based on Euclidean distance (red) associated with a cycle path (green).'----
#  plot(sa1income, col = "light grey")
#  plot(catch800m, col = rgb(1, 0, 0, 0.5), add = TRUE)
#  plot(netcatch800m, col = rgb(0, 0, 1, 0.5), add = TRUE)
#  plot(testcycleway, col = "green", add = TRUE)

## ---- echo=FALSE, message=FALSE-----------------------------------------------
#  l$d_euclidean <- line_length(l)
#  l$d_rf <- routes_fast$length

## ---- eval=FALSE--------------------------------------------------------------
#  routes_slow <- line2route(l, route_cyclestreet, plan = "quietest")

## -----------------------------------------------------------------------------
#  l$d_rq <- routes_slow$length # quietest route distance
#  Q <- mean(l$d_rf / l$d_euclidean, na.rm = TRUE)
#  QDF <- mean(l$d_rq / l$d_rf, na.rm = TRUE)
#  Q
#  QDF

## -----------------------------------------------------------------------------
#  (QDFt <- mean(routes_slow$time / routes_fast$time, na.rm = TRUE))

## ----euclidwalking1, fig.cap='Euclidean distance and walking trips', eval=FALSE----
#  l$pwalk <- l$On.foot / l$All
#  plot(l$d_euclidean, l$pwalk,
#    cex = l$All / 50,
#    xlab = "Euclidean distance (m)", ylab = "Proportion of trips by foot"
#  )

## ----euclidfastest, out.width='100%', fig.cap='Euclidean and fastest route distance of trips in the study area (left) and Euclidean distance vs the proportion of trips made by walking (right).', echo=FALSE----
#  par(mfrow = c(1, 2))
#  lgb <- sp::spTransform(l, CRSobj = sp::CRS("+init=epsg:27700"))
#  l$d_euclidean <- rgeos::gLength(lgb, byid = T)
#  l$d_rf <- routes_fast@data$length
#  plot(l$d_euclidean, l$d_rf,
#    xlab = "Euclidean distance", ylab = "Route distance"
#  )
#  abline(a = 0, b = 1)
#  abline(a = 0, b = 1.2, col = "green")
#  abline(a = 0, b = 1.5, col = "red")
#  l$pwalk <- l$On.foot / l$All
#  plot(l$d_euclidean, l$pwalk,
#    cex = l$All / 50,
#    xlab = "Euclidean distance (m)", ylab = "Proportion of trips by foot"
#  )

## -----------------------------------------------------------------------------
#  lm1 <- lm(pwalk ~ d_euclidean, data = l@data, weights = All)
#  lm2 <- lm(pwalk ~ d_rf, data = l@data, weights = All)
#  lm3 <- glm(pwalk ~ d_rf + I(d_rf^0.5),
#    data = l@data, weights = All, family = quasipoisson(link = "log")
#  )

## ---- echo=FALSE, eval=FALSE--------------------------------------------------
#  summary(lm1)
#  summary(lm2)
#  summary(lm3)

## ----euclidwalking2, fig.cap='Relationship between euclidean distance and walking', out.width="75%", fig.align="center"----
#  plot(l$d_euclidean, l$pwalk,
#    cex = l$All / 50,
#    xlab = "Euclidean distance (m)", ylab = "Proportion of trips by foot"
#  )
#  l2 <- data.frame(d_euclidean = 1:5000, d_rf = 1:5000)
#  lm1p <- predict(lm1, l2)
#  lm2p <- predict(lm2, l2)
#  lm3p <- predict(lm3, l2)
#  lines(l2$d_euclidean, lm1p)
#  lines(l2$d_euclidean, exp(lm2p), col = "green")
#  lines(l2$d_euclidean, exp(lm3p), col = "red")