4. Routing

# Select a random set of sites and facilities.
# We select random locations within the bounding box of the nodes.
set.seed(128)
hull = net %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_combine() %>%
  st_convex_hull()

sites = st_sample(hull, 50, type = "random")
facilities = st_sample(hull, 5, type = "random")

# Blend the sites and facilities into the network to get better results.
# Also select only the largest connected component.
new_net = net %>%
  activate("nodes") %>%
  filter(group_components() == 1) %>%
  st_network_blend(c(sites, facilities))
#> Warning: st_network_blend assumes attributes are constant over geometries

# Calculate the cost matrix.
cost_matrix = st_network_cost(new_net, from = sites, to = facilities, weights = "weight")

# Find for each site which facility is closest.
closest = facilities[apply(cost_matrix, 1, function(x) which(x == min(x))[1])]

# Create a line between each site and its closest facility, for visualization.
draw_lines = function(sources, targets) {
  lines = mapply(
    function(a, b) st_sfc(st_cast(c(a, b), "LINESTRING"), crs = st_crs(net)),
    sources,
    targets,
    SIMPLIFY = FALSE
  )
  do.call("c", lines)
}

connections = draw_lines(sites, closest)

# Plot the results.
plot(new_net, col = "grey")
plot(connections, col = colors[2], lwd = 2, add = TRUE)
plot(facilities, pch = 8, cex = 2, lwd = 2, col = colors[3], add = TRUE)
plot(sites, pch = 20, cex = 2, col = colors[2], add = TRUE)

set.seed(403)
rdm = net %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_sample(10, type = "random")

net = activate(net, "nodes")
cost_matrix = st_network_cost(net, from = rdm, to = rdm, weights = "weight")

# Use nearest node indices as row and column names.
rdm_idxs = st_nearest_feature(rdm, net)
row.names(cost_matrix) = rdm_idxs
colnames(cost_matrix) = rdm_idxs

round(cost_matrix, 0)
#> Units: [m]
#>      252  545  698   15  224  323   31  510   71   61
#> 252    0  914 1335 1947 1010 1122  776 1052 1365 1204
#> 545  914    0  421 1575  503  358  962  680  451  291
#> 698 1335  421    0 1552  590  241 1383 1031  227  130
#> 15  1947 1575 1552    0 2046 1408 1314 1015 1616 1620
#> 224 1010  503  590 2046    0  767 1155 1151  551  459
#> 323 1122  358  241 1408  767    0 1170  800  304  308
#> 31   776  962 1383 1314 1155 1170    0  796 1413 1252
#> 510 1052  680 1031 1015 1151  800  796    0 1094  970
#> 71  1365  451  227 1616  551  304 1413 1094    0  161
#> 61  1204  291  130 1620  459  308 1252  970  161    0

tour = solve_TSP(TSP(units::drop_units(cost_matrix)))
tour_idxs = as.numeric(names(tour))
tour_idxs
#>  [1]  71  61 545 224 252  31  15 510 323 698

# Approximate length of the route.
# In meters, since that was the unit of our cost values.
round(tour_length(tour), 0)
#> [1] 6337

# Define the nodes to calculate the shortest paths from.
# Define the nodes to calculate the shortest paths to.
# All based on the calculated order of visit.
from_idxs = tour_idxs
to_idxs = c(tour_idxs[2:length(tour_idxs)], tour_idxs[1])

# Calculate the specified paths.
tsp_paths = mapply(st_network_paths,
    from = from_idxs,
    to = to_idxs,
    MoreArgs = list(x = net, weights = "weight")
  )["node_paths", ] %>%
  unlist(recursive = FALSE)

# Plot the results.
plot(net, col = "grey")
plot(rdm, pch = 20, col = colors[2], add = TRUE)

walk(tsp_paths, plot_path) # Reuse the plot_path function defined earlier.

plot(
  st_as_sf(slice(net, rdm_idxs)),
  pch = 20, col = colors[3], add = TRUE
)
plot(
  st_as_sf(slice(net, tour_idxs[1])),
  pch = 8, cex = 2, lwd = 2, col = colors[3], add = TRUE
)
text(
  st_coordinates(st_as_sf(slice(net, tour_idxs[1]))) - c(200, 90),
  labels = "start/end\npoint"
)

# How many edge types are there?
types = net %>%
  activate("edges") %>%
  pull(type) %>%
  unique()

types
#> [1] residential  secondary    path         footway      track       
#> [6] unclassified pedestrian   service      cycleway    
#> 9 Levels: cycleway footway path pedestrian residential secondary ... unclassified

# Randomly define a walking speed in m/s for each type.
# With values between 3 and 7 km/hr.
set.seed(1)
speeds = runif(length(types), 3 * 1000 / 60 / 60, 7 * 1000 / 60 / 60)

# Assign a speed to each edge based on its type.
# Calculate travel time for each edge based on that.
net = net %>%
  activate("edges") %>%
  group_by(type) %>%
  mutate(speed = units::set_units(speeds[cur_group_id()], "m/s")) %>%
  mutate(time = weight / speed) %>%
  ungroup()

net

net = activate(net, "nodes")

p = net %>%
  st_geometry() %>%
  st_combine() %>%
  st_centroid()

iso = net %>%
  filter(node_distance_from(st_nearest_feature(p, net), weights = time) <= 600)

iso_poly = iso %>%
  st_geometry() %>%
  st_combine() %>%
  st_convex_hull()

plot(net, col = "grey")
plot(iso_poly, col = NA, border = "black", lwd = 3, add = TRUE)
plot(iso, col = colors[2], add = TRUE)
plot(p, col = colors[1], pch = 8, cex = 2, lwd = 2, add = TRUE)

# Define the threshold values (in seconds).
# Define also the colors to plot the neighborhoods in.
thresholds = rev(seq(60, 600, 60))
palette = sf.colors(n = 10)

# Plot the results.
plot(net, col = "grey")

for (i in c(1:10)) {
  nbh = convert(net, to_spatial_neighborhood, p, thresholds[i], weights = "time")
  plot(nbh, col = palette[i], add = TRUE)
}

plot(p, pch = 8, cex = 2, lwd = 2, add = TRUE)

table(roxel$type)
#> 
#>     cycleway      footway         path   pedestrian  residential    secondary 
#>           42           67           95            1          394           62 
#>      service        track unclassified 
#>          160            5           25

paths_sf = net %>%
  activate("edges") %>%
  slice(unlist(paths$edge_paths)) %>%
  st_as_sf()

table(paths_sf$type)
#> 
#>     cycleway      footway         path   pedestrian  residential    secondary 
#>            2            0            6            0           38            1 
#>      service        track unclassified 
#>            2            0            6

plot(paths_sf["type"], lwd = 4, key.pos = 4, key.width = lcm(4))

weighting_profile = c(
  cycleway = Inf,
  footway = Inf,
  path = Inf,
  pedestrian = Inf,
  residential = 3,
  secondary = 1,
  service = 1,
  track = 10,
  unclassified = 1
)

weighted_net = net %>%
  activate("edges") %>%
  mutate(multiplier = weighting_profile[type]) %>%
  mutate(weight = edge_length() * multiplier)

weighted_paths = st_network_paths(weighted_net, from = 495, to = c(458, 121), weights = "weight")

weighted_paths_sf = weighted_net %>%
  activate("edges") %>%
  slice(unlist(weighted_paths$edge_paths)) %>%
  st_as_sf()

table(weighted_paths_sf$type)
#> 
#>     cycleway      footway         path   pedestrian  residential    secondary 
#>            0            0            0            0            9           43 
#>      service        track unclassified 
#>            3            0            0

plot(st_as_sf(net, "edges")["type"], lwd = 4,
     key.pos = NULL, reset = FALSE, main = "Distance weights")
plot(st_geometry(paths_sf), add = TRUE)
plot(st_as_sf(net, "edges")["type"], lwd = 4,
     key.pos = NULL, reset = FALSE, main = "Custom weights")
plot(st_geometry(weighted_paths_sf), add = TRUE)