README.qmd

---
format: gfm
title: "Voronoiesque polygons based on travel times isochrones"
warning: false
echo: false
---

# Pre-requisites

The code underlying this paper requires R to be installed.

```{r}
#| label: setup
# #| include: false
library(sf)
library(osmextract)
library(dplyr)
library(tmap)
```

# Input data

The input datasets for the example data are as follows:

- Street network in a 1 km buffer around central Oldenburg
- 4 points in Oldenburg

```{r}
#| label: extract-osm-data
centroid = osmextract:::oe_search("Oldenburg, Germany")
poly = zonebuilder::zb_zone(centroid, n_circles = 1)
# mapview::mapview(poly)
# walking_network = oe_get_network(poly, "walking", boundary = poly, boundary_type = "clipsrc")
# plot(walking_network$geometry)
# sf::write_sf(walking_network, "oldenburg_walking_network.geojson")
walking_network = sf::read_sf("oldenburg_walking_network.geojson")
# points = oe_get(poly, extra_tags = c("amenity"), query = "SELECT * FROM points WHERE amenity = 'pub'", boundary = poly, boundary_type = "clipsrc")
# points = points |>
#   select(name, osm_id)
# mapview::mapview(points)
# points = points |>
#   filter(stringr::str_detect(name, "Ben|Gast|AU|Kar"))
# sf::write_sf(points, "oldenburg_pubs.geojson", delete_dsn = TRUE)
points = sf::read_sf("oldenburg_pubs.geojson")
tm_shape(walking_network) + tm_lines() + tm_shape(points) + tm_dots(col = "red", size = 5)
```

# Voronoi polygons

```{r}
#| label: voronois
# voronoi = stplanr::geo_projected(points, sf::st_voronoi) # fails
# crsuggest::suggest_crs(points)
# 5652
local_crs = "EPSG:5652"
points_projected = st_transform(points, local_crs)
poly_projected = st_transform(poly, local_crs)
voronoi_projected = st_voronoi(st_union(points_projected), poly_projected$geometry)
voronoi_projected_polygons = st_collection_extract(voronoi_projected, type = "POLYGON")
voronoi = st_transform(voronoi_projected_polygons, "EPSG:4326")
voronoi = sf::st_join(st_as_sf(voronoi), points |> select(osm_id))
voronoi = left_join(st_drop_geometry(points), voronoi) |> 
  st_as_sf()
voronoi = st_intersection(voronoi, poly)
tm_shape(voronoi) + tm_polygons() + tm_shape(points) + tm_dots(col = "red", size = 0.8)
```

# Isochrones

```{r}
# osrm_iso = osrm::osrmIsochrone(loc = points[2, ], breaks = c(2, 4, 6, 8), osrm.profile = "foot")
# sf::write_sf(osrm_iso, "osrm_iso_1.geojson", delete_dsn = TRUE)
osrm_iso = read_sf("osrm_iso_1.geojson")
osrm_iso
tmap_options(check.and.fix = TRUE)
tm_shape(voronoi) +
  tm_borders(lwd = 5) +
  tm_shape(osrm_iso) +
  tm_polygons(col = "isomax", alpha = 0.3) +
  tm_shape(points) +
  tm_dots()

```

# Next steps with isochrone polygon intersection approach

The example above demonstrates the calculation of voronoi polygons and isochrone polygons associated with points.
To get from this example to catchment areas associated with travel times, building on the approach of calculating multiple isochrones, a number of problems need to be solved:

- Iterative union of isochrone polygons associated with each point for which there are no 'collissions'
- In cases where there are 'collisions' between isochrone polygons, erase polygons with larger travel times with polygons associated with a different point that have lower travel times
- Where isochrone polygons of equal travel time intersect, find the centreline of the intersection and partition polygons according, as outlined [here](https://gis.stackexchange.com/questions/217151/how-to-align-edges-of-overlapping-polygons-in-the-middle-line)

# Alternative approaches

Another approach would be to iteratively sample points located between points to find locations that have roughly equal travel times.
From these 'equal travel time points' polygons can be constructed.

# Nearest hex cells

```{r}
hex_grid = stplanr::geo_projected(
  voronoi,
  st_make_grid,
  cellsize = 200,
  square = FALSE
)
hex_grid = hex_grid[poly]
tm_shape(hex_grid) + tm_polygons() +
  tm_shape(points) + tm_dots(col = "red", size = 0.8)
```

We'll iterate over every hex cell to find the nearest pub, first using nearest distances:

```{r}
hex_df = data.frame(name = NA)
hex_centroids = st_as_sf(st_centroid(hex_grid))
nearest_points = st_join(hex_centroids, points, join = nngeo::st_nn, k = 1, progress = FALSE)
hex_joined = st_sf(
  st_drop_geometry(nearest_points),
  geometry = hex_grid
)
hex_joined_centroids = st_centroid(hex_joined)
voronoi_hex = hex_joined |>
  group_by(name) |>
  summarise(n = n())
tm_shape(voronoi_hex, bb = st_bbox(voronoi)) + tm_polygons(col = "name") +
  tm_shape(points) + tm_dots(col = "red", size = 0.8) +
  tm_shape(voronoi) + tm_borders(col = "blue", lwd = 5) +
  tm_layout(legend.outside = TRUE)
```

Next, we'll use travel times to find the nearest pub.
To minimise the number of requests, the strategy will be as follows: we will identify hex cells that touch the boundary between two or more territories.

```{r, echo=TRUE}
inner_lines = rmapshaper::ms_innerlines(voronoi_hex)
hex_boundary = hex_joined[inner_lines, ]
hex_boundary_centroids = st_centroid(hex_boundary)
tm_shape(voronoi_hex, bb = st_bbox(voronoi)) + tm_polygons(col = "name") +
  tm_shape(hex_boundary) + tm_fill(col = "grey", alpha = 0.8) +
  tm_shape(points) + tm_dots(col = "red", size = 0.8) +
  tm_shape(voronoi) + tm_borders(col = "blue", lwd = 5) +
  tm_layout(legend.outside = TRUE)
```

We'll prepare the OSM network for routing.

# Routing with sfnetworks

We'll start by demonstrating how the package works with the sample dataset.

```{r}
library(sfnetworks)
roxel
net = as_sfnetwork(roxel, directed = FALSE) |>
  activate("edges") |>
  mutate(weight = edge_length())
study_area = st_convex_hull(roxel)
set.seed(2023)
rpoints = st_sample(study_area, 10)
tm_shape(roxel) + tm_lines() +
  tm_shape(rpoints) + tm_dots(col = "red", size = 0.8)
```

We'll calculate the route from point 1 to point 2:

```{r}
path_1_2 = st_network_paths(net, rpoints[1], rpoints[2], weights = "weight")
path_1_2_sf = net |> 
  activate("edges") |> 
  slice(path_1_2$edge_paths[[1]]) |>
  sf::st_as_sf()
tm_shape(roxel) + tm_lines() +
  tm_shape(rpoints[1:2]) + tm_dots(col = "red", size = 0.8) +
  tm_shape(path_1_2_sf) + tm_lines(lwd = 5, col = "blue")
```

```{r}
#| echo: false
#| eval: false
#| label: with-nodes-on-network
net_nodes = net |> 
  activate("nodes") |>
  st_as_sf()
from_graph = nngeo::st_nn(
  rpoints[1],
  net_nodes, k = 1,
  progress = FALSE
)[[1]]
to_graph = nngeo::st_nn(
  rpoints[2],
  net_nodes, k = 1,
  progress = FALSE
)[[1]]
path_1_2 = st_network_paths(net, from_graph, to_graph)
path_1_2_sf = net |> 
  activate("edges") |> 
  slice(path_1_2$edge_paths[[1]]) |>
  sf::st_as_sf()
tm_shape(roxel) + tm_lines() +
  tm_shape(rpoints[1:2]) + tm_dots(col = "red", size = 0.8) +
  tm_shape(path_1_2_sf) + tm_lines(lwd = 5, col = "blue")
```

We can calculate many routes as follows:

```{r}
net_nodes = net |> 
  activate("nodes") |>
  st_as_sf()
point_ids = nngeo::st_nn(
  rpoints,
  net_nodes, k = 1,
  progress = FALSE
) |> unlist()
point_df = data.frame(
  from = rep(point_ids, each = length(point_ids)),
  to = rep(point_ids, length(point_ids))
) |>
  filter(from != to)

paths_all = st_network_paths(
  net,
  from = point_df$from,
  to = point_df$to,
  weights = "weight"
)
class(paths_all)
routes_list = lapply(seq(nrow(paths_all)), function(i) {
  net |> 
    activate("edges") |> 
    slice(paths_all$edge_paths[[i]]) |>
    mutate(route_number = i) |>
    sf::st_as_sf()
})
routes_list[[1]]
routes_sf = do.call(rbind, routes_list)
tm_shape(roxel) + tm_lines() +
  tm_shape(rpoints) + tm_dots(col = "red", size = 0.8) +
  tm_shape(routes_sf) + tm_lines(lwd = 5, col = "blue", alpha = 0.05)
```

We can calculate the amount of travel on each link as follows:

```{r}
routes_sf$n = 1
rnet = stplanr::overline(routes_sf, "n")
tm_shape(rnet) + tm_lines(lwd = "n", scale = 9)
```

```{r, echo=TRUE}
net_linestrings = sf::st_cast(walking_network, "LINESTRING")
net = sfnetworks::as_sfnetwork(net_linestrings, directed = FALSE)
library(tidygraph)
with_graph(net, graph_component_count())
net = net |>
  activate("edges") |>
  mutate(weight = edge_length()) |>
  activate("nodes") |>
  filter(group_components() == 1)
with_graph(net, graph_component_count())

net_sf = net |> 
  sfnetworks::activate("edges") |> 
  sf::st_as_sf() |> 
  dplyr::select(from, to, weight)
nrow(net_sf)
nrow(walking_network)
tm_shape(walking_network) + tm_lines("grey", lwd = 5) +
  tm_shape(net_sf) + tm_lines("blue", lwd = 2) 
```

We'll start by calculating routes from the first `hex_boundary` cell to the nearest point.

```{r}
net_nodes = net |> 
  activate("nodes") |>
  st_as_sf()
from_point = hex_joined_centroids[1, ]
to_point = points[1, ]
path = sfnetworks::st_network_paths(net, from_point, to_point)
path_sf = net |> 
  activate("edges") |> 
  slice(path$edge_paths[[1]]) |>
  sf::st_as_sf()

tm_shape(voronoi_hex, bb = st_bbox(voronoi)) + tm_polygons(col = "name") +
  tm_shape(hex_joined[1, ]) + tm_fill(col = "black") +
  tm_shape(hex_boundary) + tm_fill(col = "grey", alpha = 0.8) +
  tm_shape(points) + tm_dots(col = "red", size = 0.8) +
  tm_shape(voronoi) + tm_borders(col = "blue", lwd = 5) +
  tm_shape(path_sf) + tm_lines()
  tm_layout(legend.outside = TRUE)
```

## Calculation of shortest paths in boundary cell

A logical next step is to calculate the shortest path to n nearest (in Euclidean distance) destinations for 'boundary cells'.
We do this for the first boundary cell as follows:

```{r}
n = 3
first_boundary_point = hex_boundary_centroids[1, ]
nearest_point_ids = nngeo::st_nn(
  first_boundary_point,
  points, k = n,
  progress = FALSE
)[[1]]
nearest_points = points[nearest_point_ids, ]
# plot the result
tm_shape(voronoi_hex, bb = st_bbox(voronoi)) + tm_polygons(col = "name") +
  tm_shape(hex_boundary[1, ]) + tm_fill(col = "black") +
  tm_shape(hex_boundary) + tm_fill(col = "grey", alpha = 0.8) +
  tm_shape(points) + tm_dots(col = "red", size = 0.8) +
  tm_shape(voronoi) + tm_borders(col = "blue", lwd = 5) +
  tm_shape(nearest_points) + tm_dots(col = "green", size = 0.8) +
  tm_layout(legend.outside = TRUE)
```

Next we'll calculate the paths, keeping the total length of each path:

```{r}
paths = st_network_paths(
  net,
  from = first_boundary_point,
  to = nearest_points,
  weights = "weight"
)
path_1 = net |> 
  activate("edges") |> 
  slice(paths$edge_paths[[1]]) |>
  mutate(route_number = 1) |>
  sf::st_as_sf()
sum(path_1$weight)

path_weights = sapply(seq(nrow(paths)), function(i) {
  net |> 
    activate("edges") |> 
    slice(paths$edge_paths[[i]]) |>
    mutate(route_number = i) |>
    sf::st_as_sf() |>
    summarise(length = sum(weight)) |>
    pull(length)
})
point_shortest_id = which.min(path_weights)
point_shortest = nearest_points[point_shortest_id, ]
cell_value_original = first_boundary_point$name
cell_value_new = point_shortest$name
cell_value_original
cell_value_new
```

As shown, the pub associated with the shortest path is different from the pub associated with the original cell.
We will update the cell value to reflect this:

```{r}
which_hex = which(lengths(st_intersects(hex_joined, first_boundary_point)) > 0)
hex_iso = hex_joined
hex_iso$name[which_hex] 
hex_iso$name[which_hex] = cell_value_new
m1 = tm_shape(hex_joined) + tm_polygons(col = "name")
m2 = tm_shape(hex_iso) + tm_polygons(col = "name")
tmap_arrange(m1, m2)
```

We'll now repeat this process for all boundary cells:

```{r}
i = 2
for(i in seq(nrow(hex_boundary))) {
  first_boundary_point = hex_boundary_centroids[i, ]
  nearest_point_ids = nngeo::st_nn(
    first_boundary_point,
    points, k = n,
    progress = FALSE
  )[[1]]
  nearest_points = points[nearest_point_ids, ]
  paths = st_network_paths(
    net,
    from = first_boundary_point,
    to = nearest_points,
    weights = "weight"
  )
  path_weights = sapply(seq(nrow(paths)), function(i) {
    net |> 
      activate("edges") |> 
      slice(paths$edge_paths[[i]]) |>
      mutate(route_number = i) |>
      sf::st_as_sf() |>
      summarise(length = sum(weight)) |>
      pull(length)
  })
  point_shortest_id = which.min(path_weights)
  point_shortest = nearest_points[point_shortest_id, ]
  cell_value_original = first_boundary_point$name
  cell_value_new = point_shortest$name
  which_hex = which(lengths(st_intersects(hex_joined, first_boundary_point)) > 0)
  hex_iso$name[which_hex] = cell_value_new
}
# Plot the results next to the original:
tm_points = tm_shape(points) + tm_dots(col = "red", size = 0.8)
tm_net = tm_shape(net_sf) + tm_lines(col = "black")
m1_net = m1 + tm_net + tm_points
m2_net = tm_shape(hex_iso) + tm_polygons(col = "name") + tm_net + tm_points
tmap_arrange(m1_net, m2_net)
```

# Routing with cppRouting

```{r}
#| eval: false
net_df = net_sf |> 
  sf::st_drop_geometry()
head(net_df)
net_from = lwgeom::st_startpoint(net_sf)
net_to = lwgeom::st_endpoint(net_sf)
graph = cppRouting::makegraph(net_df, directed = FALSE)
names(graph)
names(graph$data)
summary(graph$data$dist)

#| label: single-path
head(graph$data)
head(graph$dict)
str(graph)
# calculate route from A to B with cppRouting
from_graph = nngeo::st_nn(
  hex_joined_centroids[1, ],
  net_from, k = 1,
  progress = FALSE
)[[1]]
to_graph = nngeo::st_nn(
  points[1, ],
  net_to, k = 1,
  progress = FALSE
)[[1]]
# cppRouting::get_path_pair(graph, from[rep(1, nrow(to))], to)
route_cpp = cppRouting::get_path_pair(graph, 1, 2, long = TRUE)
cppRouting::get_path_pair(graph, from_graph[1], to_graph[1])
route_cpp_df = cppRouting::get_path_pair(graph, from_graph[1], to_graph[1], long = TRUE)
str(route_cpp_df)
route_cpp_df_to_join = route_cpp_df |> 
  transmute(from = as.numeric(node))
# route_cpp1 = net_sf[as.numeric(route_cpp[[1]]), ]
route_cpp = inner_join(route_cpp_df_to_join, net_sf) |>
  st_as_sf()
tm_shape(voronoi_hex, bb = st_bbox(voronoi)) + tm_polygons(col = "name") +
  tm_shape(hex_joined[1, ]) + tm_fill(col = "black") +
  tm_shape(hex_boundary) + tm_fill(col = "grey", alpha = 0.8) +
  tm_shape(points) + tm_dots(col = "red", size = 0.8) +
  tm_shape(voronoi) + tm_borders(col = "blue", lwd = 5) +
  tm_shape(route_cpp) + tm_lines() +
  tm_layout(legend.outside = TRUE)
```

# Next steps

- Debug the results of the `sfnetworks` approach
- Get the `cppRouting` approach working
- Test out different routing backends, to overcome issues with local routing
- Improve networks used for routing with network cleaning approaches