README.qmd

---
title: netvis
# format: html
format: gfm
number-sections: true
toc: true
editor:
  render-on-save: true
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  warning = FALSE,
  message = FALSE,
  cache = TRUE
)
library(tidyverse)
library(tmap)
tmap_mode("view")
```

<!-- badges: start -->

<!-- badges: end -->

The goal of this repo is to demonstrate different visualisation techniques for communicating information about transport networks.

```{r, eval = FALSE, echo=FALSE}
# install.packages("devtools")
# devtools::install_github("ITSLeeds/netvis")
library(tidyverse)
```

# Example data

We'll use a route network dataset from the Propensity to Cycle Tool (PCT) to demonstrate the package.
The PCT is a web application that allows users to explore cycling potential across England and Wales.
The PCT is available at [www.pct.bike](https://www.pct.bike/).

We also provide an example dataset from Limerick, Ireland, from the CRUSE project.
This is imported and plotted with code shown below.

```{r}
#| eval: false
#| echo: false
u = "https://cruse.bike/limerick/rnet_combined_balanced.gpkg"
rnet_limerick_balanced_full = sf::read_sf(u)
# Make numeric cols numerics not chars:
# To do this, values of "1-10" need to be changed to 3 (or another integer)
rnet_limerick_balanced_full = rnet_limerick_balanced_full |> 
  mutate(across(starts_with("Bic"), as.numeric)) |> 
  mutate(across(starts_with("Bic"), replace_na, replace = 3)) 
zones_limerick = zonebuilder::zb_zone("limerick")
central_limerick = zones_limerick |>
  dplyr::filter(circle_id == 1) |>
  sf::st_union()
rnet_limerick = sf::st_intersection(rnet_limerick_balanced_full, central_limerick)

rnet_limerick = rnet_limerick |>
  mutate(across(c(`Bicycle (Baseline)`:`Bicycle (Ebike)`), ~as.numeric(.)))

sf::st_write(rnet_limerick, "test-data/rnet_limerick.geojson", delete_dsn = TRUE)
usethis::use_data(rnet_limerick)
```

```{r}
rnet_limerick = sf::read_sf("test-data/rnet_limerick.geojson")
plot(rnet_limerick)
```

# Static data visualisation

A simple visualisation of the data in a multi-panel static map is shown below.

```{r}
#| eval: false
#| echo: false
rnet_leeds = pct::get_pct_rnet("west-yorkshire")
rnet_leeds = rnet_leeds |>
  select(-1)
zones_leeds = zonebuilder::zb_zone("Leeds")
central_leeds = zones_leeds |>
  dplyr::filter(circle_id == 1) |>
  sf::st_union()
rnet_central = sf::st_intersection(rnet_leeds, central_leeds)
# sf::st_write(rnet_central, "test-data/rnet_central.geojson", delete_dsn = TRUE)
rnet_minimal = rnet_central |> 
  slice_max(ebike_slc, n = 2) |>
  mutate(across(where(is.numeric), function(x) x * 1:2))
rnet_minimal
usethis::use_data(rnet_minimal, overwrite = TRUE)
```

Read in data for Leeds with the following:

```{r}
rnet_central = sf::read_sf("test-data/rnet_central.geojson")
plot(rnet_central)
usethis::use_data(rnet_central, overwrite = TRUE)
```

In Python this would look something like this:

```{python}
import geopandas as gpd
rnet_central = gpd.read_file("test-data/rnet_central.geojson")
rnet_central.plot()
```

# Interactive maps

Interactive maps are more engaging and policy-relevant, allowing people to zoom in and explore results.
The images below illustrate interactive maps of the route network data focussed on the `bicycle` and `dutch_slc` variables.

```{r, echo=FALSE}
rnet_central = rnet_central |>
  select(bicycle, dutch_slc) 
m1 = tm_shape(rnet_central) +
  tm_lines(
    lwd = "bicycle",
    scale = 9
  )
tmap_save(m1, "maps/m1.html")
# browseURL("maps/m1.html")
# webshot2::webshot("maps/m1.html")
width_multiplier = round(max(rnet_central$dutch_slc) / max(rnet_central$bicycle))
```

The maximum level of flow in the Go Dutch scenario is `r round(width_multiplier)` times the maximum level of flow in the baseline scenario.
We can use this to scale the line widths as illustrated below.

```{r}
m2 = tm_shape(rnet_central) +
  tm_lines(
    lwd = "dutch_slc",
    scale = 9 * width_multiplier
  )
tmap_save(m2, "maps/m2.html")
# browseURL("maps/m2.html")
if(!file.exists("maps/m1.png")) {
  webshot2::webshot("maps/m1.html", "maps/m1.png")
  webshot2::webshot("maps/m2.html", "maps/m2.png")
}
```

::: {#m1-2 layout-ncol="2"}
![](maps/m1.png){#fig-surus}

![](maps/m2.png){#fig-hanno}
:::

There are two problems with line widths in the maps shown above:

-   The thinnes lines are too thin
-   The thickest lines are too thick

```{r, echo=FALSE}
max_width = 15
width_multiplier = 5
```

Given that the maximum width is determined by the `scale` argument, we can solve the first problem by increasing the value of `scale`.
The second problem can be solved with a multiplier that prevents lines being x times thicker than the thinnest lines, `r width_multiplier` times in this case.

```{r, echo=FALSE}
# get 95th percentile of line widths
combined_values = rnet_central |>
  sf::st_drop_geometry() |>
  sapply(function(x) {
    quantile(x, 0.95)
  }
  )
minimum_value_allowed = max(combined_values) / width_multiplier
```

```{r, echo=FALSE, out.width="50%", fig.show="hold"}
scale_bicycle = max_width / (quantile(rnet_central$dutch_slc, 0.95) / quantile(rnet_central$bicycle, 0.95))
# summary(rnet_central$bicycle)
rnet_bicycle = rnet_central |>
    mutate(lwd = case_when(
      bicycle < minimum_value_allowed ~ minimum_value_allowed,
      TRUE ~ bicycle
    ))
# summary(rnet_bicycle$lwd)
m3 = tm_shape(
  rnet_bicycle
 ) +
  tm_lines(
    lwd = "lwd",
    scale = scale_bicycle
  )
# tmap_save(m3, "maps/m3.html")
# browseURL("maps/m3.html")

# For the Go Dutch scenario
scale_dutch = max_width
m4 = tm_shape(
  rnet_central |>
    mutate(lwd = case_when(
      dutch_slc < minimum_value_allowed ~ minimum_value_allowed,
      TRUE ~ dutch_slc
    ))
 ) +
  tm_lines(
    lwd = "lwd",
    scale = scale_dutch
  )
# tmap_save(m4, "maps/m4.html")
# browseURL("maps/m4.html")
if(!file.exists("maps/m3.png")) {
  webshot2::webshot("maps/m3.html", "maps/m3.png")
  webshot2::webshot("maps/m4.html", "maps/m4.png")
}
```

::: {layout-ncol="2"}
![](maps/m3.png)

![](maps/m4.png)
:::

A function that does this would look something like this:

```{r}
scale_line_widths = function(x, max_width, min_width) {
  # ...
}
```

In the plot below, for example, we increase the thickness of the thinnest lines by setting `min_width` to 3, making the lines with less flow more visible.

```{r, echo=FALSE}
x = rnet_central |>
  select(bicycle, dutch_slc)
scale_line_widths = function(
  x,
  max_width = 15,
  min_width = 1,
  ptile = 0.99,
  width_regex = "bicycle|dutch_slc",
  width_multiplier = NULL
) {
  names_width = names(x)[grepl(width_regex, names(x))]
  names_width_lwd = paste0(names_width, "_lwd")
  names(names_width_lwd) = names_width
  if(is.null(width_multiplier)) {
    width_multiplier = max_width / min_width
  }
  x_no_outliers = x |>
    sf::st_drop_geometry() |>
    select(matches(width_regex)) |>
    mutate_all(function(x) {
      x[x > quantile(x, ptile)] = quantile(x, ptile)
      x
    })
  # waldo::compare(x, x_no_outliers) # large values gone
  max_value = max(x_no_outliers)
  x_normalized = x_no_outliers |>
    mutate_all(function(x) x / max_value)
  # summary(x_normalized)
  minimum_value_allowed = 1 / width_multiplier
  x_scaled = x_normalized |>
    mutate_all(function(x) case_when(
      x < minimum_value_allowed ~ minimum_value_allowed,
      TRUE ~ x
    ))
  summary(x_scaled)
  maximums = x_scaled |> sapply(max)
  max_widths = maximums * max_width
  names(x_scaled) = names_width_lwd
  x_to_plot = cbind(x, x_scaled)
  nm = names_width[1]
  map_list = lapply(names_width, function(nm) {
    tm_shape(x_to_plot) +
      tm_lines(
        lwd = names_width_lwd[nm],
        scale = max_widths[[nm]]
      )
      })
  # map_list[[1]]
  map_list
}
```

```{r}
m5 = scale_line_widths(x, max_width = 15, min_width = 3)
# m5[[1]]
# m5[[2]]
```

```{r, echo=FALSE, out.width="50%", fig.show="hold"}
tmap_save(m5[[1]], "maps/m5.html")
# browseURL("maps/m5.html")
# webshot2::webshot("maps/m5.html", "maps/m5.png")
tmap_save(m5[[2]], "maps/m6.html")
# browseURL("maps/m6.html")
# webshot2::webshot("maps/m6.html", "maps/m6.png")
```

::: {layout-ncol="2"}
![](maps/m5.png)

![](maps/m6.png)
:::

Instead of showing different layers side-by-side, it's useful to be able to select them as different layers interactively.
The `netvis()` function defined in this packages does that, as illustrated below.

![](https://user-images.githubusercontent.com/1825120/253802338-83fdcbea-3a43-4ae7-b2d6-db6810e9c89c.png)

<!-- Python attempt: -->

```{python}
#| eval: false
#| echo: false
def scale_line_widths(
    data,
    baseline_column,
    dutch_column,
    max_width = 15,
    width_multiplier = 5
):
    """
    Scale line widths for interactive maps.

    Parameters
    ----------
    data : pandas.DataFrame
        Dataframe containing the data to be plotted.
    baseline_column : str
        Name of the column containing the baseline values.
    dutch_column : str
        Name of the column containing the Go Dutch values.
    max_width : int, optional
        Maximum width of the lines in the interactive map. The default is 15.
    width_multiplier : int, optional
        The minimum width of the lines is the maximum value in the baseline
        data divided by this value. The default is 5.

    Returns
    -------
    data : pandas.DataFrame
        Dataframe with a new column containing the scaled line widths.

    """
    # get 95th percentile of line widths
    combined_values = data[[baseline_column, dutch_column]].quantile(0.95)
    # Set mininum value to the maximum value in the combined data divided by the minimum width multiplier
    minumum_value_allowed = combined_values.max() / width_multiplier
    # scale line widths
    scale = max_width / (combined_values[dutch_column] / combined_values[baseline_column])
    data["lwd_dutch"] = data[dutch_column].apply(lambda x: minimum_value_allowed if x < minimum_value_allowed else x)
    data["lwd_bicycle"] = data[baseline_column].apply(lambda x: minimum_value_allowed if x < minimum_value_allowed else x)
    return data, scale
```

<!-- Starting with the data in `test-data/rnet_central.geojson` let's  test this function on a small subset of the data. -->

```{python}
#| eval: false
#| echo: false
import pandas as pd
import geopandas as gpd
import numpy as np
data = gpd.read_file("test-data/rnet_central.geojson")
data = data[["bicycle", "dutch_slc"]]
data = data.iloc[0:10, :]
dw, scale = scale_line_widths(data, "bicycle", "dutch_slc")
print(dw)
```

```{r}
#| eval: false
#| echo: false

# https://github.com/nptscot/outputdata/releases/tag/v2023-07-11-00-00-48.247456_commit_282d4661618e5e60ad9c34a6368db84745473d18
rnet_edinburgh = readRDS("/tmp/combined_network.Rds")
sf::write_sf(rnet_edinburgh, "~/github/nptscot/networkmerge/large_route_network_example_edingurgh.geojson")
zones_edinburgh = zonebuilder::zb_zone("princes street edinburgh")
mapview::mapview(zones_edinburgh)
# library(osmdata)
# princes_street_all = opq("edinburgh") %>%
#   add_osm_feature(key = "name", value = "Princes Street") %>% 
#   osmdata_sf()
osm_edinburgh = osmextract::oe_get("Edinburgh")
princes_street = osm_edinburgh |> 
  filter(name == "Princes Street")
plot(princes_street$geometry)
princes_street_buffer = princes_street |> 
  sf::st_union() |> 
  sf::st_buffer(dist = 100)
rnet_princes_street = rnet_edinburgh[princes_street_buffer, , op = sf::st_within]
plot(rnet_princes_street$geometry)
rnet_princes_street = rnet_princes_street |> 
  select(commute_fastest_bicycle_go_dutch, Quietness)
rnet_princes_street |> plot()
sf::write_sf(rnet_princes_street, "data/rnet_princes_street_minimal.geojson")
sf::write_sf(rnet_princes_street, "~/github/nptscot/networkmerge/data/rnet_princes_street_minimal.geojson")

```