Riparian Hydro‑Ecology

Author

Inyo County Water Department

Original code flattened the hydrograph It’s a limitation of the ggplot2 package by design. ggplot2 deliberately discourages dual‑axis plots. Hadley Wickham’s design rationale is that a figure should have one visual scale per aesthetic; two different units mapped to the same aesthetic (y‑position) can mislead. If dual axes are unavoidable, the relationship must be linear so viewers can mentally transform one to the other.

Because of that philosophy, ggplot2 offers only a re‑labelling tool (sec_axis), not a second independent scale.
Hence the manual work you have to do with when forcing ggplot to do this:

  1. Calculate a scale factor k so the secondary data fit on the primary axis.
  2. Multiply the data in the layer (flow_acft * k).
  3. Provide the inverse transformation (~ . / k) in sec.axis.
  4. Format tick labels yourself (scales::comma_format()).

It’s easier to do this in python - letting you call ax.twinx() and be done, but in ggplot2 the extra arithmetic keeps the core grammar simple and guards against misleading charts—at the cost of annoyances and manual work when forcing it to do dual axes.

  • Before, the histogram’s axis stopped around 15 trees, so flows > 100 000 acre‑ft were squeezed into 0 – 15 so the line looked flat.

In a nutshell, you have to first rescale the flow values

k <- max(tree_cnt$n_trees) / max(flow_df$flow_acft)
geom_line(aes(y = flow_acft * k, colour = station))

Then undo that scaling on the secondary axis

scale_y_continuous(
  name     = "Number of Trees",
  sec.axis = sec_axis(
    ~ . / k,
    name   = "Flow (acre‑ft)",
    labels = scales::comma_format()
  )
)

here’s a reproducible code below. Note - I wrote two functions which can be updated and called more easily than copying code over and over to iterate. Hope this helps!

libraries 
library(tidyverse)
library(here)
library(scales)
functions 
#────────────────────────────────────────────────────────────
# 1) multi‑station flow
#────────────────────────────────────────────────────────────
plot_tree_vs_flow <- function(tree_df,
                              flow_df,
                              binwidth   = 1,
                              fill_cols  = c(off      = "brown4",
                                             on       = "turquoise4",
                                             onfire   = "darkorange1",
                                             `onfire C` = "gold",
                                             onsuppl  = "blue2")) {

  ## histogram counts ---------------------------------------------------------
  tree_cnt <- tree_df %>%
    count(strata_clone, comb_date, name = "n_trees")

  facets  <- unique(tree_cnt$strata_clone)
  # min_yr  <- min(tree_cnt$comb_date, na.rm = TRUE)
  min_yr  <- 1885
  # max_yr  <- max(tree_cnt$comb_date, na.rm = TRUE)
  max_yr <- 2025

  ## flow data (drop stray strata column if present) --------------------------
  flow_base <- flow_df %>%
    select(-any_of("strata_clone")) %>%
    mutate(comb_date = as.integer(comb_date))

  k <- max(tree_cnt$n_trees, na.rm = TRUE) /
       max(flow_base$flow_acft, na.rm = TRUE)

  # flow_plot <- flow_base %>%
  #   mutate(flow_scaled = flow_acft * k) %>%
  #   tidyr::crossing(strata_clone = facets)


  ## figure -------------------------------------------------------------------
  ggplot() +
    geom_col(data = tree_cnt,
             aes(comb_date, n_trees, fill = strata_clone),
             width = binwidth, colour = "black") +

    # geom_line(data = flow_plot,
    #           aes(comb_date, flow_scaled, colour = station),
    #           size = .7, alpha = .8) +
      geom_line(data = flow_base,
                aes(comb_date, flow_acft * k, colour = station),
                size = .7, alpha = .8)+
    scale_y_continuous(
      name     = "Number of Trees",
      sec.axis = sec_axis(~ . / k,
                          name   = "Flow (acre‑ft)",
                          labels = comma_format())
    ) +
    scale_x_continuous(
      breaks       = seq(min_yr, max_yr, by = 5),  # label every 5 yrs
      minor_breaks = seq(min_yr, max_yr, by = 1),   # grid every year
      name         = "Recruitment year",
      labels       = number_format(accuracy = 1)
    ) +
    scale_fill_manual(values = fill_cols, name = NULL) +
    scale_colour_brewer(palette = "Dark2", name = "Station") +
    facet_grid(strata_clone ~ .) +
    theme_light() +
    theme(axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6),
          legend.position = "bottom")
}

#────────────────────────────────────────────────────────────
# 2) OV runoff
#────────────────────────────────────────────────────────────
plot_tree_vs_runoff <- function(tree_df,
                                runoff_df,
                                binwidth   = 1,
                                fill_cols  = c(off      = "brown4",
                                               on       = "turquoise4",
                                               onfire   = "darkorange1",
                                               `onfire C` = "gold",
                                               onsuppl  = "blue2")) {

  tree_cnt <- tree_df %>%
    count(strata_clone, comb_date, name = "n_trees")

  facets  <- unique(tree_cnt$strata_clone)
  # min_yr  <- min(tree_cnt$comb_date, na.rm = TRUE)
  min_yr  <- 1885
  # max_yr  <- max(tree_cnt$comb_date, na.rm = TRUE)
  max_yr <- 2025

  runoff_base <- runoff_df %>%
    select(-any_of("strata_clone")) %>%
    mutate(comb_date = as.integer(comb_date))

  k <- max(tree_cnt$n_trees, na.rm = TRUE) /
       max(runoff_base$ovr_acft, na.rm = TRUE)

  # runoff_plot <- runoff_base %>%
  #   mutate(run_scaled = ovr_acft * k) %>%
  #   tidyr::crossing(strata_clone = facets)

  ggplot() +
    geom_col(data = tree_cnt,
             aes(comb_date, n_trees, fill = strata_clone),
             width = binwidth, colour = "black") +

    # geom_line(data = runoff_plot,
    #           aes(comb_date, run_scaled),
    #           colour = "black", size = .9) +
    geom_line(data = runoff_base,
                aes(comb_date, ovr_acft * k), colour = "black", size = .9)+
    scale_y_continuous(
      name     = "Number of Trees",
      sec.axis = sec_axis(~ . / k,
                          name   = "Total Runoff (acre‑ft)",
                          labels = comma_format())
    ) +
    scale_x_continuous(
      breaks       = seq(min_yr, max_yr, by = 5),
      minor_breaks = seq(min_yr, max_yr, by = 1),
      name         = "Recruitment year",
      labels       = number_format(accuracy = 1)
    ) +
    scale_fill_manual(values = fill_cols, name = NULL) +
    facet_grid(strata_clone ~ .) +
    theme_light() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6),
          legend.position = "bottom")
}
load-data 
## load data
trees   <- read_csv(here("data","cored_trees_plt.csv"),  show_col_types = FALSE)
flows   <- read_csv(here("data","cored_trees_hy_test.csv"), show_col_types = FALSE)
runoff  <- read_csv(here("data","cored_trees_run.csv"), show_col_types = FALSE)
call-multi-station 
## multi‑station figure
plot_tree_vs_flow(trees, flows)

call-total-runoff 
## OV‑runoff figure
plot_tree_vs_runoff(trees, runoff)