---
title: "Desire Lines Extended"
author: "Robin Lovelace, Jakub Nowosad, Jannes Muenchow"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Desire Lines Extended}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

This vignette builds on the [transport chapter](https://geocompr.robinlovelace.net/transport.html) of [the Geocomputation with R book](https://geocompr.github.io/) by showing how to create multi-stage desire lines, from the ground-up.

It depends on these packages and datasets:

```{r, message=FALSE, warning=FALSE}
library(sf)
library(stplanr)
library(dplyr)
library(spDataLarge)
desire_lines = od2line(bristol_od, bristol_zones)
desire_rail = top_n(desire_lines, n = 3, wt = train)
```

The first stage is to create matrices of coordinates that will subsequently be used to create matrices representing each leg:

```{r}
mat_orig = as.matrix(line2df(desire_rail)[c("fx", "fy")])
mat_dest = as.matrix(line2df(desire_rail)[c("tx", "ty")])
mat_rail = st_coordinates(bristol_stations)
```

The outputs are three matrices representing the starting points of the trips, their destinations and possible intermediary points at public transport nodes (named `orig`, `dest` and `rail` respectively).
But how to identify *which* intermediary points to use for each desire line?
The `knn()` function from the **nabor** package (which is used internally by **stplanr** so it should already be installed) solves this problem by finding *k nearest neighbors* between two sets of coordinates.
By setting the `k` parameter, one can define how many nearest neighbors should be returned. 
Of course, `k` cannot exceed the number of observations in the input (here: `mat_rail`).
We are interested in just one nearest neighbor, namely, the closest railway station:

```{r}
knn_orig = nabor::knn(mat_rail, query = mat_orig, k = 1)$nn.idx
knn_dest = nabor::knn(mat_rail, query = mat_dest, k = 1)$nn.idx
```

This results not in matrices of coordinates, but row indices that can subsequently be used to subset the `mat_rail`.
It is worth taking a look at the results to ensure that the process has worked properly, and to explain what has happened:

```{r}
as.numeric(knn_orig)
as.numeric(knn_dest)
```

The output demonstrates that each object contains three whole numbers (the number of rows in `desire_rail`) representing the rail station closest to the origin and destination of each desire line.
Note that while each 'origin station' is different, the destination (station `30`) is the same for all desire lines.
This is to be expected because rail travel in cities tends to converge on a single large station (in this case Bristol Temple Meads).
The indices can now be used to create matrices representing the rail station of origin and destination:

```{r}
mat_rail_o = mat_rail[knn_orig, ]
mat_rail_d = mat_rail[knn_dest, ]
```

The final stage is to convert these matrices into meaningful geographic objects, in this case simple feature 'multilinestrings' that capture the fact that each stage is a separate line, but part of the same overall trip:

```{r}
mats2line = function(mat1, mat2) {
  lapply(1:nrow(mat1), function(i) {
    rbind(mat1[i, ], mat2[i, ]) %>%
    st_linestring()
  }) %>% st_sfc()
}
desire_rail$leg_orig = mats2line(mat_orig, mat_rail_o)
desire_rail$leg_rail = mats2line(mat_rail_o, mat_rail_d)
desire_rail$leg_dest = mats2line(mat_rail_d, mat_dest)
```


```{r, eval=FALSE, echo=FALSE}
# Create equivalent geographic objects using multilinestrings:
tmp = lapply(1:nrow(mat_orig), function(i) {
      list(
        st_linestring(rbind(mat_orig[i, ], mat_rail_o[i, ])), 
        st_linestring(rbind(mat_rail_o[i, ], mat_rail_d[i, ])),
        st_linestring(rbind(mat_rail_d[i, ], mat_dest[i, ]))) %>%
  st_multilinestring
  })
tmp = st_sfc(tmp)
plot(tmp, type = "n")
lapply(tmp, function(x) plot(x[[1]], col = "red", add = TRUE))
lapply(tmp, function(x) plot(x[[2]], col = "blue", add = TRUE))
lapply(tmp, function(x) plot(x[[3]], col = "green", add = TRUE))
```

The results are visualised below:

```{r stations, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Station nodes (red dots) used as intermediary points that convert straight desire lines with high rail usage (black) into three legs: to the origin station (red) via public transport (grey) and to the destination (a very short blue line)."}
zone_cents = st_centroid(bristol_zones)
zone_cents_rail = zone_cents[desire_rail,]
plot(desire_rail$geometry, expandBB = c(.1, .1, .1, .1))
plot(desire_rail$leg_orig, add = TRUE, col = "red", lwd = 3)
plot(desire_rail$leg_rail, add = TRUE, col = "grey", lwd = 2)
plot(bristol_stations, add = TRUE, col = "red")
plot(desire_rail$leg_dest, add = TRUE, col = "blue", lwd = 5)
plot(zone_cents_rail, add = TRUE, col = "black")
# library(tmap)
# tmap_mode("plot")
# qtm(bristol_stations, basemaps = "https://{s}.tile.thunderforest.com/transport/{z}/{x}/{y}.png?apikey=feae177da543411c9efa64160305212d", dots.col = "red", symbols.size = 2) +
#   tm_shape(desire_rail) +
#   tm_lines(col = "black", lwd = 4) +
#   tm_shape(legs) +
#   tm_lines()
```