Heatmaps and other surfaces

Overview

This notebook provides examples of creating heatmaps and other colour surfaces to plot outcomes as colours with two driver axes (e.g. temperature and rainfall changes).

Demonstration setup

As usual, we need paths to the data. We use the ‘more scenarios’ examples for all the plots, with processing as in the website workflow.

project_dir <- file.path("more_scenarios")
hydro_dir <- file.path(project_dir, "hydrographs")
agg_dir <- file.path(project_dir, "aggregator_output")

Read in the data

We read in the example data we will use for all plots.

agged_data <- readRDS(file.path(agg_dir, "achievement_aggregated.rds"))

That has all the steps in the aggregation, but most of the plots here will only use a subset to demonstrate.

To make visualisation easier, the SDL units data is given a grouping column that puts the many env_obj variables in groups defined by their first two letters, e.g. EF for Ecosystem Function. These correspond to the ‘Target’ level, but it can be useful to have the two groupings together for some examples.

If we had used multiple aggregation functions at any step, we should filter down to the one we want here, but we only used one for this example.

scenarios <- yaml::read_yaml(file.path(hydro_dir, "scenario_metadata.yml")) |>
  tibble::as_tibble()

# join the scenario information to each of those sheets
agged_data <- agged_data |>
  purrr::map(\(x) left_join(x, scenarios, by = "scenario"))

Make surface plots

The main thing to remember when making heatmaps is to use plot_type = 'heatmap', and include x_col, y_col and outcome_col (the outcome plotted as colour) arguments.

Qualitative axes

If we have qualitative descriptions on x- and y-, the heatmap will have even steps. Note here the simple cleanup of the names and the need to make the adapt_code column non-numeric.

qual_heatmap <- agged_data$Target |>
  dplyr::filter(adapt_code %in% c(1, 2, 3, 4)) |>
  dplyr::filter(!is.na(Target)) |>
  # clean names
  mutate(
    SWSDLName = stringr::str_replace(SWSDLName, "–", "-\n"),
    Target = stringr::str_replace(Target, "Priority e", "E"),
    Target = stringr::str_wrap(Target, width = 12),
    adapt_code = as.character(adapt_code)
  ) |>
  plot_outcomes(
    outcome_col = "ewr_achieved",
    outcome_lab = "Condition",
    y_col = "adapt_code",
    y_lab = "Adaptation option",
    x_col = "climate_code",
    x_lab = "Climate scenario",
    plot_type = "heatmap",
    colorset = "ewr_achieved",
    pal_list = "grDevices::Viridis",
    facet_row = "Target",
    facet_col = "SWSDLName"
  )

qual_heatmap

Quantitative axes

Here, our scenarios are defined quantitatively, and so we can allow x and y to be quantitative. A key issue here is the spacing of the scenarios. We log-transform the axes to help a bit, but if this were a targetted output, even scenario spacing would be desirable.

# Quantitative axes- it's more informative, but uglier
quant_heatmap <- agged_data$Target |>
  dplyr::filter(adapt_code %in% c(1, 2, 3)) |>
  # clean names
  mutate(
    SWSDLName = stringr::str_replace(SWSDLName, "–", "-\n"),
    Target = stringr::str_replace(Target, "Priority e", "E"),
    Target = stringr::str_wrap(Target, width = 12)
  ) |>
  mutate(flow_addition = flow_addition + 1) |>
  plot_outcomes(
    outcome_col = "ewr_achieved",
    outcome_lab = "Condition",
    y_col = "flow_addition",
    y_lab = "Flow addition ('adaptation')",
    x_col = "flow_multiplier",
    x_lab = "Flow multiplier ('climate')",
    plot_type = "heatmap",
    transy = "log10",
    transx = "log10",
    colorset = "ewr_achieved",
    pal_list = "grDevices::Viridis",
    facet_row = "Target",
    facet_col = "SWSDLName"
  )

quant_heatmap

Interpolated raster

We can try to fill in the missing locations by interpolating a raster using the contour_arglist, which passes arguments to geom_raster (if it contains ‘interpolate’) or geom_contour (if it does not). That helps, but the scenarios are still really too far apart, and so closer scenario spacing would be better.

quant_interp_heatmap <- agged_data$Target |>
  dplyr::filter(adapt_code %in% c(1, 2, 3)) |>
  # clean names
  mutate(
    SWSDLName = stringr::str_replace(SWSDLName, "–", "-\n"),
    Target = stringr::str_replace(Target, "Priority e", "E"),
    Target = stringr::str_wrap(Target, width = 12)
  ) |>
  mutate(flow_addition = flow_addition + 1) |>
  plot_outcomes(
    outcome_col = "ewr_achieved",
    outcome_lab = "Condition",
    y_col = "flow_addition",
    y_lab = "Flow addition ('adaptation')",
    x_col = "flow_multiplier",
    x_lab = "Flow multiplier ('climate')",
    plot_type = "heatmap",
    transy = "log10",
    transx = "log10",
    colorset = "ewr_achieved",
    pal_list = "grDevices::Viridis",
    facet_row = "Target",
    facet_col = "SWSDLName",
    contour_arglist = list(interpolate = TRUE)
  )

quant_interp_heatmap
Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.
Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.

Contours

We can also just fit contour surfaces, done by passing a list of arguments to contour_arglist. If it is an empty list, as here, it simply triggers the contour. Again, these contours would be better if they were fit to tighter scenario spacing.

contour_heatmap <- agged_data$Target |>
  dplyr::filter(adapt_code %in% c(1, 2, 3)) |>
  # clean names
  mutate(
    scenario = paste0(climate_code, adapt_code),
    SWSDLName = stringr::str_replace(SWSDLName, "–", "-\n"),
    Target = stringr::str_replace(Target, "Priority e", "E"),
    Target = stringr::str_wrap(Target, width = 12),
    adapt_code = as.character(adapt_code)
  ) |>
  mutate(flow_addition = flow_addition + 1) |>
  plot_outcomes(
    outcome_col = "ewr_achieved",
    outcome_lab = "Proportion\nEWR achieved",
    y_col = "flow_addition",
    y_lab = "Flow addition ('adaptation')",
    x_col = "flow_multiplier",
    x_lab = "Flow multiplier ('climate')",
    plot_type = "heatmap",
    transy = "log10",
    transx = "log10",
    colorset = "ewr_achieved",
    pal_list = "grDevices::Viridis",
    facet_row = "Target",
    facet_col = "SWSDLName",
    contour_arglist = list()
  )

contour_heatmap

We can use contour_arglist to pass arguments to geom_contour, e.g. changing the bins, binwidth, or breaks.

contour_heatmap_adj <- agged_data$Target |>
  dplyr::filter(adapt_code %in% c(1, 2, 3)) |>
  # clean names
  mutate(
    scenario = paste0(climate_code, adapt_code),
    SWSDLName = stringr::str_replace(SWSDLName, "–", "-\n"),
    Target = stringr::str_replace(Target, "Priority e", "E"),
    Target = stringr::str_wrap(Target, width = 12),
    adapt_code = as.character(adapt_code)
  ) |>
  mutate(flow_addition = flow_addition + 1) |>
  plot_outcomes(
    outcome_col = "ewr_achieved",
    outcome_lab = "Proportion\nEWR achieved",
    y_col = "flow_addition",
    y_lab = "Flow addition ('adaptation')",
    x_col = "flow_multiplier",
    x_lab = "Flow multiplier ('climate')",
    plot_type = "heatmap",
    transy = "log10",
    transx = "log10",
    colorset = "ewr_achieved",
    pal_list = "grDevices::Viridis",
    facet_row = "Target",
    facet_col = "SWSDLName",
    contour_arglist = list(bins = 30)
  )

contour_heatmap_adj

Overplotting points

At present, we cannot plot points on top of the heatmap internally to plot_outcomes(), but it should be possible in future with overlay_list(). Currently, we can take advantage of the fact that these are ggplot objects and add another geom:

contour_heatmap +
  geom_point(
    data = scenarios |>
      filter(adapt_code %in% c(1, 2, 3)),
    color = "black"
  )
Warning in ggplot2::scale_y_continuous(trans = transy): log-10 transformation
introduced infinite values.

Heatmaps for EWRs

Heatmaps can also be useful for examining the EWRs that come in as input. Since they are unique by gauge, planning-unit, and sdl unit, the names can get ugly, but it can give a good sense of how the EWRs themselves are performing. Here, we’ll look over time, and only for one catchment. If the names were simpler, we could facet by catchment or adaptation option.

agged_data$all_time |>
  dplyr::filter(adapt_code %in% c(1) &
    climate_code %in% c("A", "E", "I") &
    SWSDLName == "Namoi") |>
  dplyr::mutate(unique_ewr = paste0(ewr_code_timing, gauge, planning_unit_name)) |>
  # clean names
  plot_outcomes(
    outcome_col = "ewr_achieved",
    y_col = "unique_ewr",
    x_col = "climate_code",
    x_lab = "Climate scenario",
    plot_type = "heatmap",
    colorset = "ewr_achieved",
    pal_list = "grDevices::Viridis",
    facet_wrapper = "SWSDLName",
    scales = "free_y"
  )