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

## ---- 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 <- readr::read_csv("fun_table.csv")
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 <- buff_geo(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_cyclestreet(from = "Bradford", to = "Leeds")
#  route_c1_c2 <- route_cyclestreet(cents[1,], cents[2,])

## ---- eval=FALSE---------------------------------------------------------
#  route_bl_raw <- route_cyclestreet(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'---------------------------------------
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")

