evaluation.qmd

---
title: "Evaluating network-based interaction zones with real-world trajectory data"
format:
  html:
    self-contained: true
editor: visual
---

In the scope of the Bike2CAV project we have developed a workflow that automatically detects spatial areas on a traffic intersection in which it is expected that interactions between bicycles and cars will be taking place. This detection is purely based on the layout of the intersection as originally planned, with dedicated lanes of movement in different directions and for different vehicles. In this notebook we evaluate how accurate these detected zones are when considering how people actually move through the intersection, which might not always align to the original plan. We base this evaluation on trajectory data gathered at the intersection.

## Prepare

Load required libraries.

```{r}
#| message: false

library(sf)
library(tidyverse)
library(spatstat)
library(stars)
library(viridis)
library(units)
library(here)
```

Load utility functions.

```{r}
#| message: false

source(here("utils.R"))
```

Fetch all required data. They are stored in a single zip file. Some of the input datasets are already pre-processed. To get the zip file, please contact lucas.vandermeer\@plus.ac.at.

```{r}
unzip(here("data", "data.zip"), overwrite = TRUE, exdir = here("data"))
```

Read network-based interaction zones, as detected by earlier work in the project.

```{r}
zones = read_sf(here("data", "zones.geojson"))
```

Read trajectory data. These come in two separate datasets. The first one contains all observed sample points, with the time and location of the observation, and the index of the object being observed. The second one contains the lines that are drawn by chronologically connecting all the points belonging to the same object.

```{r}
all_points = read_sf(here("data", "trajectories.gpkg"), layer = "points")
all_lines = read_sf(here("data", "trajectories.gpkg"), layer = "lines")
```

Read other data. These include the polygonized lanelets of the intersection to be analyzed, and the focus area of that intersection, which we will limit our analysis to.

```{r}
lanes = read_sf(here("data", "lanes.geojson")) # Lanelets as polygons.
focus = read_sf(here("data", "focus.geojson")) # Focus area of the intersection.
```

## Pre-process

### Transform data

Transform all data into the same coordinate reference system.

```{r}
proj = 31287

all_points = st_transform(all_points, proj)
all_lines = st_transform(all_lines, proj)

lanes = st_transform(lanes, proj)
zones = st_transform(zones, proj)
focus = st_transform(focus, proj)
```

### Tidy trajectory data

Tidy up the trajectory data, by renaming columns and updating column types.

```{r}
names(all_points) = gsub(":", "_", names(all_points))
names(all_points) = gsub("[.]", "_", names(all_points))
st_agr(all_points) = setNames(st_agr(all_points), setdiff(names(all_points), "geom"))

all_points$track_id = as.character(all_points$track_id)
all_points$t = as.POSIXct(all_points$t / 1000, origin = "1970-01-01")
```

```{r}
all_lines = all_lines |>
  rename(t_start = startTimestamp, t_end = endTimestamp, t_update = updateTimestamp)

names(all_lines) = gsub(":", "_", names(all_lines))
names(all_lines) = gsub("[.]", "_", names(all_lines))
st_agr(all_lines) = setNames(st_agr(all_lines), setdiff(names(all_lines), "geom"))

all_lines$t_start = as.POSIXct(all_lines$t_start / 1000, origin = "1970-01-01")
all_lines$t_end = as.POSIXct(all_lines$t_end / 1000, origin = "1970-01-01")
all_lines$t_update = as.POSIXct(all_lines$t_update / 1000, origin = "1970-01-01")
```

### Filter trajectories

Calculate the length and displacement of each individual trajectory. The displacement is the straight-line distance from the start to the end of the trajectory.

```{r}
all_lines$length = st_length(all_lines)
all_lines$displacement = st_displacement(all_lines)
```

Remove all trajectories with a displacement lower than 10 or higher than 70 meters.

```{r}
lower = set_units(10, "m")
upper = set_units(70, "m")
sub_lines = filter(all_lines, displacement > lower & displacement < upper)
sub_points = filter(all_points, track_id %in% sub_lines$id)
```

Remove all trajectories that do not intersect with the focus area.

```{r}
sub_lines = st_filter(sub_lines, focus, .predicate = st_intersects)
sub_points = filter(sub_points, track_id %in% sub_lines$id)
```

```{r}
ggplot() + 
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(all_lines)), color = "steelblue", alpha = 0.5) +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank()) +
  ggtitle("All trajectories")
```

```{r}
ggplot() + 
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(sub_lines)), color = "steelblue", alpha = 0.5) +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank()) +
  ggtitle("Filtered trajectories")
```

### Group trajectories by transport mode

Define the transport mode of each trajectory. This workflow utilizes two attributes in the trajectory data: the detected *object type*, which makes a distinction between vehicles and pedestrians, and the detected *vehicle type*, which makes a distinction between different types of vehicles (e.g. car, bike) whenever the object type is vehicle.

First, extract the last detected object and vehicle type of each trajectory.

```{r}
sub_lines = sub_lines |>
  rename(object_type_last = object_type, vehicle_type_last = vehicle_type)

last = sub_lines |>
  st_drop_geometry() |>
  select(id, object_type_last, vehicle_type_last)
```

Assign each individual trajectory point an object and vehicle type.

```{r}
sub_points = sub_points |>
  rename(object_type = tag_changed_object_type, vehicle_type = tag_changed_vehicle_type)

is_last = !duplicated(sub_points$track_id, fromLast = TRUE)

last_types = sub_points |>
  select(track_id) |>
  left_join(last, by = c(track_id = "id"))

sub_points$object_type[is_last] = last_types$object_type_last[is_last]
sub_points$vehicle_type[is_last] = last_types$vehicle_type_last[is_last]

sub_points = fill(sub_points, object_type, vehicle_type, .direction = "up")
```

Define the dominating object and vehicle type of each trajectory. This is the type that is most often assigned to a point within a trajectory. However, if less than 75% of the points in a trajectory have this type assigned, the dominating type is instead set to "unclear".

```{r}
dominator = function(x, threshold = 0.75) {
  if (all(is.na(x))) return(NA)
  tab = table(x)
  mod = names(tab)[which(tab == max(tab))[1]]
  if ((tab[mod] / sum(tab)) < threshold) {
    out = "unclear"
  } else {
    out = mod
  }
  out
}

dominant = sub_points |>
  st_drop_geometry() |>
  group_by(track_id) |>
  summarise(object_type_dominant = dominator(object_type), vehicle_type_dominant = dominator(vehicle_type))

sub_lines = left_join(sub_lines, dominant, by = c(id = "track_id"))
```

Define the transport mode of each trajectory by mapping dominant object and vehicle types to a specific transport mode.

```{r}
sub_lines$mode = case_when(
  sub_lines$object_type_dominant == "pedestrian" ~ "pedestrian",
  sub_lines$object_type_dominant == "unclear" ~ "unclear",
  sub_lines$vehicle_type_dominant == "passengerCar" ~ "car",
  sub_lines$vehicle_type_dominant == "bike" ~ "bike",
  sub_lines$vehicle_type_dominant == "motorcycle" ~ "motorcycle",
  sub_lines$vehicle_type_dominant == "bus" ~ "bus",
  sub_lines$vehicle_type_dominant == "heavyTruck" ~ "truck",
  sub_lines$vehicle_type_dominant == "unclear" ~ "unclear",
  TRUE ~ "unknown"
)

modes = sub_lines |>
  st_drop_geometry() |>
  select(id, mode)

sub_points = left_join(sub_points, modes, by = c(track_id = "id"))
```

```{r}
ggplot(modes |> group_by(mode) |> summarise(n = n())) + 
  geom_bar(aes(x = n, y = mode), stat = "identity", fill = "steelblue") +
  ggtitle("Number of trajectories per transport mode")
```

Group trajectories into car trajectories and bike trajectories.

```{r}
bike_points = filter(sub_points, mode == "bike")
bike_lines = filter(sub_lines, mode == "bike")

car_points = filter(sub_points, mode == "car")
car_lines = filter(sub_lines, mode == "car")
```

```{r}
ggplot() + 
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(car_lines)), color = "steelblue", alpha = 0.5) + 
  geom_sf(aes(geometry = st_geometry(bike_lines)), color = "violetred", alpha = 0.5) + 
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank()) +
  ggtitle("Trajectories for cars (blue) and bikes (red)")
```

## Evaluate

Map the network-based interaction zones to be evaluated. The darker grey area is the focus area of the intersection, in which we will conduct our analysis.

```{r}
ggplot() + 
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(focus)), fill = "black", color = NA, alpha = 0.3) +
  geom_sf(aes(geometry = st_geometry(zones)), fill = "steelblue", alpha = 0.8) +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), legend.position = "none") +
  ggtitle("Network-based interaction zones")
```

### Detect interactions in trajectory data

Find locations where bike and car trajectories interacted inside the focus area of the intersection.

#### Approach I: Post-encroachment times

This approach finds all locations where a segment of a bike trajectory crosses a segment of a car trajectory, and then calculated the post-encroachment time of that crossing. The post-encroachment time in this case is defined as the time between the moment that the cyclist leaves the trajectory of the car driver and the moment that the car driver reaches the trajectory of the cyclist. In that way, the post-encroachment indicates the extent to which they missed each other. See [here](https://www.youtube.com/watch?v=rSABRwcVvUU) for a visual explanation. To define an interaction, the approach uses a threshold of 2 seconds (i.e. a crossing of trajectories is considered an interaction if the post-encroachment time is 2 seconds or less).

```{r}
bike_segments = st_points_to_segments(bike_points, linestring_id = "track_id")

bike_times = bike_points |> 
  st_drop_geometry() |> 
  select(track_id, t) |> 
  group_by(track_id) |> 
  group_split() |> 
  lapply(\(x) tibble(t0 = x$t[1:(nrow(x) - 1)], t1 = x$t[2:nrow(x)])) |> 
  bind_rows()

bike_segments = bind_cols(bike_segments, bike_times)
```

```{r}
car_segments = st_points_to_segments(car_points, linestring_id = "track_id")

car_times = car_points |> 
  st_drop_geometry() |> 
  select(track_id, t) |> 
  group_by(track_id) |> 
  group_split() |> 
  lapply(\(x) tibble(t0 = x$t[1:(nrow(x) - 1)], t1 = x$t[2:nrow(x)])) |> 
  bind_rows()

car_segments = bind_cols(car_segments, car_times)
```

```{r}
#| warning: false

segment_crossings = st_crossings(bike_segments, car_segments) |>
  remove_rownames() |>
  st_as_tibble() |>
  rename(segment_id_bike = segment_id, track_id_bike = track_id) |>
  rename(segment_id_car = segment_id.1, track_id_car = track_id.1) |>
  mutate(t_bike = t0 + (t1 - t0) / 2, t_car = t0.1 + (t1.1 - t0.1) / 2) |>
  mutate(pet = abs(t_car - t_bike)) |>
  select(ends_with("_bike"), ends_with("car"), pet)
```

```{r}
pet_interactions = segment_crossings |>
  filter(pet <= 2) |>
  st_filter(focus)
```

```{r}
ggplot() +
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(focus)), fill = "black", color = NA, alpha = 0.3) +
  geom_sf(aes(geometry = st_geometry(pet_interactions)), cex = 0.8, alpha = 0.8) +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank()) +
  ggtitle("Trajectory-based interaction points", subtitle = "Approach I: Post-encroachment times")
```

#### Approach II: Space-time prisms

This approach utilizes the concept of [space-time prisms](https://link.springer.com/referenceworkentry/10.1007/978-3-319-17885-1_1599). For each individual bike trajectory point, it finds those car trajectory points that are close in space (within 2 meters), and those car trajectory points that are close in time (within 2 seconds). A bike trajectory point becomes an interaction point if the same car trajectory point occurs in both the set of spatial neighbors and the set of temporal neighbors.

```{r}
spatial_neighbors = function(x, y = x, threshold = 2) {
  st_contains(st_buffer(x, threshold), y)
}

temporal_neighbors = function(x, y = x, threshold = 2, time_column = "time") {
  xtime = x[[time_column]]
  ytime = y[[time_column]]
  lapply(xtime, \(x) which(ytime > (x - threshold) & ytime < (x + threshold)))
}

sn = spatial_neighbors(bike_points, car_points)
tn = temporal_neighbors(bike_points, car_points, time_column = "t")

bike_points$interaction = mapply(\(x, y) any(x %in% y), sn, tn)
```

```{r}
prism_interactions = bike_points |>
  filter(interaction) |>
  st_filter(focus)
```

```{r}
ggplot() +
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(focus)), fill = "black", color = NA, alpha = 0.3) +
  geom_sf(aes(geometry = st_geometry(prism_interactions)), cex = 0.8, alpha = 0.8) +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank()) +
  ggtitle("Trajectory-based interaction points", subtitle = "Approach II: Space-time prisms")
```

### Compare interaction points with interaction zones

First, choose an approach for the detection of interaction points within trajectory data.

```{r}
approach = "pet" # Can be set to "pet" or "prism"
interactions = get(paste0(approach, "_interactions"))
```

Compare the number of trajectory-based interaction points that are located inside a network-based interaction zone to those that are located outside of a network-based interaction zone.

We see here that the majority of trajectory-based interaction points are located outside of any network-based interaction zone.

```{r}
interactions = st_join(interactions, select(zones, id)) |>
  rename(zone_id = id)

interactions$location = case_when(is.na(interactions$zone_id) ~ "out", TRUE ~ "in")
```

```{r}
ggplot() + 
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") + 
  geom_sf(aes(geometry = st_geometry(focus)), fill = "black", color = NA, alpha = 0.3) +
  geom_sf(aes(geometry = st_geometry(zones)), fill = "steelblue", alpha = 0.8) +
  geom_sf(aes(geometry = st_geometry(interactions)), cex = 0.8, alpha = 0.8) +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), legend.position = "none") +
  ggtitle("Network-based interaction zones and trajectory-based interaction points")
```

```{r}
counts = interactions |>
  st_drop_geometry() |>
  group_by(location) |>
  summarise(count = n()) |>
  mutate(percentage = count / nrow(interactions) * 100)

counts
```

```{r}
ggplot(counts) + 
  geom_bar(aes(x = count, y = location), fill = "steelblue", stat = "identity") +
  ggtitle("Number of interaction points inside and outside of interaction zones")
```

Using quadrat counts of interaction occurrences on a 10x10 regular grid covering the focus area, we can assess if the point pattern of interaction locations is completely spatially random. We see here that this hypothesis can be discarded (given the small p-value), i.e. the observed pattern is likely not created by a random process.

```{r}
#| warning: false

pattern = as.ppp(c(st_geometry(focus), st_geometry(interactions))) |>
  rjitter()

quadrat.test(pattern, nx = 10, ny = 10)
```

Using basic kernel density estimation, we can get a better idea of where high density clusters of trajectory-based interaction points are located, and compare that to the network-based interaction zones. We see here that many high-density clusters are located outside of the network-based zones, and only few inside of them.

```{r}
density = density(pattern, sigma = 1.5) |>
  st_as_stars() |>
  st_as_sf(as_points = FALSE) |>
  rename(density = v)
```

```{r}
ggplot() + 
  geom_sf(aes(geometry = st_geometry(lanes)), fill = "grey90", color = "grey85") +
  geom_sf(aes(geometry = st_geometry(density), fill = density$density), lwd = 0.1) +
  geom_sf(aes(geometry = st_geometry(zones)), fill = NA, color = "white") +
  labs(fill = "density") +
  scale_fill_viridis(option = "plasma") +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), legend.position = c(0.9, 0.7), legend.title = element_text(size = 8), legend.text = element_text(size = 6), legend.background = element_rect(fill = NA)) +
  ggtitle("Gridded interaction point density overlayed with interaction zone borders")
```

## Discussion and conclusion

By far most trajectory-based interaction points are found outside of the network-based interaction zones. Many areas with a high density of trajectory-based interaction points are not located inside any network-based interaction zone. This implies that the detected interaction zones do not accurately reflect the real-world situation. The network-based interaction zones expect users of specific modes to follow specific lanes, and interact only in small areas where these lanes cross. The trajectory-based interaction points however show a pattern in which many high-density clusters follow the shape of certain lanes, meaning that there is more mixed-used of lanes than we would expect from the network design. An explanation for this could be that especially cyclists do in many cases not behave according to what the traffic planners had in mind when designing the intersection.

A few notes need to me made:

-   The analyzed intersection has traffic lights, which minimize the possibility for bike-car interactions in many directions. These traffic lights where *not* taken into account in the detection of the interaction zones (this was based solely on the HD map of the intersection). It would be interesting to analyze an uncontrolled intersection, i.e. without traffic lights, instead.
-   Many of the inferred interaction points are part of bike trajectories that move into a tunnel where bikes normally would not go. It is not clear if these trajectories are real bike trajectories, or that they are misclassified.

To obtain more accurate interaction zones, a solution could be to run micro-simulation models that simulate movements of bikes and cars, trained with collected trajectory data. However, transferring calibrated model parameters to other intersections should be handled with care, as shown by [previous research](https://www.sciencedirect.com/science/article/abs/pii/S0001457515300385). Furthermore, it should be critically evaluated if all interactions are by definition bad and should always be avoided. A forced separation between traffic modes to minimize interactions may lead to a public space design that negatively affects well-being of the people making use of it.