Rough Seas Suppress Feeding — They Don’t Enhance It

Seabird feeding rates decline sharply as sea states worsen — calm-to-slight conditions yield the highest foraging activity.

Analysis of seabird feeding behavior across ocean conditions using 1969–1990 ship log data from New Zealand waters. A two-panel visualization reveals that feeding rates peak under calm-to-slight sea states and decline steadily as conditions worsen — evidence of opportunistic rather than forced foraging. Built with ggplot2 and patchwork, using Wilson confidence intervals for proportion estimation.

Figure 1: Two-panel data visualization titled “Rough Seas Suppress Feeding — They Don’t Enhance It.” Panel A is a heatmap showing survey effort by wind condition (Beaufort scale, binned Calm to Gale+) and sea state (SS1–SS6). Observation density peaks at slight-to-moderate seas with light-to-moderate winds, shown in deep navy. Panel B is a dot plot with Wilson 95% confidence intervals showing seabird feeding rates by sea state. Feeding peaks at 11.9% under calm, rippled conditions (SS1) and declines steadily through rough (SS5, ~4%) and very rough seas (SS6, ~2%). A dashed reference line marks the peak feeding rate. Data from Te Papa Tongarewa, Museum of New Zealand, 1969–1990.

Steps to Create this Graphic

1. Load Packages & Setup

```{r}
#| label: load
#| warning: false
#| message: false      
#| results: "hide"     

## 1. LOAD PACKAGES & SETUP ----
suppressPackageStartupMessages({
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
    tidyverse, ggtext, showtext, janitor, ggrepel,      
    scales, glue, skimr, patchwork, binom  
    )
})

### |- figure size ----
camcorder::gg_record(
  dir    = here::here("temp_plots"),
  device = "png",
  width  = 10,
  height = 8,
  units  = "in",
  dpi    = 320
)

# Source utility functions
suppressMessages(source(here::here("R/utils/fonts.R")))
source(here::here("R/utils/social_icons.R"))
source(here::here("R/utils/image_utils.R"))
source(here::here("R/themes/base_theme.R"))
```

2. Read in the Data

```{r}
#| label: read
#| include: true
#| eval: true
#| warning: false

tt <- tidytuesdayR::tt_load(2026, week = 15)
beaufort_scale <- tt$beaufort_scale |> clean_names()
birds <- tt$birds |> clean_names()
sea_states <- tt$sea_states |> clean_names()
ships <- tt$ships |> clean_names()
rm(tt)
```

3. Examine the Data

```{r}
#| label: examine
#| include: true
#| eval: true
#| results: 'hide'
#| warning: false

glimpse(beaufort_scale)
glimpse(birds)
glimpse(sea_states)
glimpse(ships)
```

4. Tidy Data

```{r}
#| label: tidy
#| warning: false

# Join birds to ship conditions via record_id
# Filter sentinel count values (99999 = "over 100,000")
# Only use full 10-minute census records for comparability

bird_ship <- birds |>
    filter(!is.na(record_id)) |>
    inner_join(ships, by = "record_id") |>
    filter(
        census_method == "full",
        count < 99999,
        !is.na(sea_state_class),
        !is.na(wind_speed_class)
    )

### |- Panel A: Environmental Envelope ----
# Collapse 12 Beaufort classes into 5 meaningful bins for readability.

beaufort_bins <- tibble::tribble(
    ~wind_speed_class, ~wind_bin, ~wind_bin_order,
    0,  "Calm\n(Bf 0–1)",     1,
    1,  "Calm\n(Bf 0–1)",     1,
    2,  "Light\n(Bf 2–3)",    2,
    3,  "Light\n(Bf 2–3)",    2,
    4,  "Moderate\n(Bf 4–5)", 3,
    5,  "Moderate\n(Bf 4–5)", 3,
    6,  "Strong\n(Bf 6–7)",   4,
    7,  "Strong\n(Bf 6–7)",   4,
    8,  "Gale+\n(Bf 8–11)",   5,
    9,  "Gale+\n(Bf 8–11)",   5,
    10, "Gale+\n(Bf 8–11)",   5,
    11, "Gale+\n(Bf 8–11)",   5
)

envelope_data <- bird_ship |>
    distinct(record_id, wind_speed_class, sea_state_class) |>
    left_join(beaufort_bins, by = "wind_speed_class") |>
    left_join(
        sea_states |> select(sea_state_class, sea_state_description),
        by = "sea_state_class"
    ) |>
    count(wind_bin, wind_bin_order, sea_state_class, sea_state_description,
          name = "n_obs") |>
    mutate(
        wind_label = fct_reorder(wind_bin, wind_bin_order),
        sea_label  = glue("SS {sea_state_class}: {str_to_title(sea_state_description)}"),
        sea_label  = fct_reorder(sea_label, sea_state_class)
    )

### |- Panel B: Feeding Rate × Sea State ----
# For each sea state class: compute feeding rate with Wilson CI
# Feeding rate = proportion of bird observations where feeding = TRUE
# Restrict to sea states with >= 30 observations for reliability

feeding_rate_data <- bird_ship |>
    filter(!is.na(feeding)) |>
    group_by(sea_state_class) |>
    summarise(
        n_obs      = n(),
        n_feeding  = sum(feeding, na.rm = TRUE),
        .groups    = "drop"
    ) |>
    filter(n_obs >= 30) |>
    # Wilson confidence intervals (statistically appropriate for proportions)
    mutate(
        ci    = map2(n_feeding, n_obs, ~ binom.wilson(.x, .y)),
        rate  = map_dbl(ci, ~ .x$mean),
        lower = map_dbl(ci, ~ .x$lower),
        upper = map_dbl(ci, ~ .x$upper)
    ) |>
    left_join(
        sea_states |> select(sea_state_class, sea_state_description),
        by = "sea_state_class"
    ) |>
    # Complete sea states 0-6 after CI computation so ghost rows never hit binom.wilson.
    # Missing states join sea_states for labels, then bind cleanly as NA rows.
    (\(computed) {
        present <- unique(computed$sea_state_class)
        missing <- setdiff(0:6, present)
        if (length(missing) > 0) {
            ghost <- tibble::tibble(sea_state_class = missing) |>
                left_join(
                    sea_states |> select(sea_state_class, sea_state_description),
                    by = "sea_state_class"
                ) |>
                mutate(
                    n_obs = 0L, n_feeding = 0L,
                    ci = list(NULL), rate = NA_real_, lower = NA_real_, upper = NA_real_
                )
            dplyr::bind_rows(computed, ghost)
        } else {
            computed
        }
    })() |>
    mutate(
        sea_label = glue("SS {sea_state_class}\n{str_to_title(sea_state_description)}"),
        sea_label = fct_reorder(sea_label, sea_state_class),
        # na.rm = TRUE so ghost NA rows don't poison max()
        is_peak   = !is.na(rate) & rate == max(rate, na.rm = TRUE)
    )
```

5. Visualization Parameters

```{r}
#| label: params
#| include: true
#| warning: false

### |-  plot aesthetics ----
colors <- get_theme_colors(
    palette = list(
        "accent"   = "#1B6CA8",   
        "low"      = "#D9EAF7",  
        "high"     = "#0D3B6E",   
        "peak"     = "#722F37",   
        "gray"     = "gray75",
        "bg"       = "#F8F9FA"
    )
)

### |- titles and caption ----
title_text    <- str_glue("Rough Seas Suppress Feeding — They Don't Enhance It")

subtitle_text <- str_glue(
    "Seabird feeding rates **decline sharply as sea states worsen** — calm-to-slight conditions<br>",
    "yield the highest foraging activity. Surveys clustered in moderate wind and wave conditions<br>",
    "*(Panel A)*, yet feeding rates *(Panel B)* decline steadily beyond slight conditions — suggesting<br>",
    "birds feed opportunistically, not because rough water forces prey to the surface."
)

caption_text <- create_social_caption(
    tt_year = 2026,
    tt_week = 15,
    source_text = "Te Papa Tongarewa — Museum of New Zealand"
)

### |-  fonts ----
setup_fonts()
fonts <- get_font_families()

### |-  plot theme ----
base_theme <- create_base_theme(colors)

weekly_theme <- extend_weekly_theme(
    base_theme,
    theme(
        # Panel
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color = "gray92", linewidth = 0.3),
        panel.grid.minor   = element_blank(),
        
        # Axes
        axis.ticks = element_blank(),
        axis.text = element_text(family = fonts$text, size = 7.5, color = "gray40"),
        axis.title = element_text(family = fonts$text, size = 8.5, color = "gray30"),
        
        # Strip (for any facets)
        strip.text  = element_text(family = fonts$text, size = 8, face = "bold"),
        
        # Legend
        legend.position = "right",
        legend.title = element_text(family = fonts$text, size = 7.5),
        legend.text = element_text(family = fonts$text, size = 7),
        legend.key.size = unit(0.4, "cm")
    )
)

theme_set(weekly_theme)
```

6. Plot

```{r}
#| label: plot
#| warning: false

### |-  Panel A: Environmental Envelope heatmap ----

p_envelope <- envelope_data |>
  ggplot(aes(x = wind_label, y = sea_label, fill = n_obs)) +
  # Geoms
  geom_tile(color = "white", linewidth = 0.6) +
  scale_fill_gradient(
    low = colors$palette$low,
    high = colors$palette$high,
    name = "Observations",
    labels = comma,
    guide = guide_colorbar(barwidth = 0.4, barheight = 4)
  ) +
  # Scales
  scale_x_discrete(position = "bottom") +
  # Labs
  labs(
    title = "Panel A — Survey Conditions",
    subtitle = "Surveys concentrate in slight–moderate seas with light–moderate winds — rougher conditions were rarely sampled",
    x = "Wind Condition",
    y  = "Sea State"
  ) +
  # Them
  theme(
    plot.title = element_text(
      family = fonts$title, size = 10, face = "bold", color = "gray20",
      margin = margin(b = 3)
    ),
    plot.subtitle = element_text(
      family = fonts$text, size = 7.5, color = "gray45",
      margin = margin(b = 8)
    ),
    axis.text.x = element_text(size = 7.5, lineheight = 1.3),
    axis.text.y = element_text(size = 7.5),
    panel.grid = element_blank()
  )

### |-  Panel B: Feeding Rate × Sea State ----

# Identify peak for annotation
peak_row <- feeding_rate_data |> filter(is_peak)

p_feeding <- feeding_rate_data |>
  ggplot(aes(x = sea_label, y = rate)) +

  # Geoms
  geom_pointrange(
    aes(ymin = lower, ymax = upper, color = is_peak),
    linewidth = 0.7,
    size = 0.2,
    show.legend = FALSE
  ) +
  geom_hline(
    yintercept = peak_row$rate,
    color      = "gray80",
    linewidth  = 0.4,
    linetype   = "dashed"
  ) +
  geom_text(
    aes(y = lower - 0.007, label = glue("n={comma(n_obs)}")),
    size = 2.3,
    color = "gray55",
    family = fonts$text
  ) +
  # Annotate
  annotate(
    "text",
    x = as.character(peak_row$sea_label),
    y = peak_row$upper + 0.009,
    label = glue("Peak feeding\n({percent(peak_row$rate, accuracy = 0.1)})"),
    size = 2.8,
    color = colors$palette$peak,
    family = fonts$text,
    lineheight = 1.1,
    vjust = 0
  ) +
  # Scales
  scale_y_continuous(
    labels = percent_format(accuracy = 0.1),
    expand = expansion(mult = c(0.14, 0.18))
  ) +
  scale_x_discrete(drop = FALSE) +
  scale_color_manual(
    values = c("FALSE" = colors$palette$accent, "TRUE" = colors$palette$peak)
  ) +
  # Labs
  labs(
    title = "Panel B — Feeding Rate by Sea State",
    subtitle = "Proportion of observations with active feeding · Wilson 95% CI · full censuses only · n ≥ 30",
    x = "Sea State",
    y = "Feeding Rate"
  ) +
  # Theme
  theme(
    plot.title = element_text(
      family = fonts$title, size = 10, face = "bold", color = "gray20",
      margin = margin(b = 3)
    ),
    plot.subtitle = element_text(
      family = fonts$text, size = 7.5, color = "gray45",
      margin = margin(b = 8)
    ),
    axis.text.x = element_text(size = 7.5, lineheight = 1.3),
    panel.grid.major.x = element_blank()
  )

### |-  Combined layout ----
combined_plot <- p_envelope / p_feeding +
  plot_layout(heights = c(1, 1.4)) +
  plot_annotation(
    title = title_text,
    subtitle = subtitle_text,
    caption = caption_text,
    theme = theme(
      plot.title = element_text(
        family = fonts$title,
        size = 18,
        face = "bold",
        color  = "gray10",
        lineheight = 1.1,
        margin = margin(b = 6)
      ),
      plot.subtitle = element_markdown(
        family = fonts$text,
        size = 10,
        color = "gray35",
        lineheight = 1.5,
        margin = margin(b = 16)
      ),
      plot.caption = element_markdown(
        family = fonts$text,
        size = 7,
        color = "gray55",
        hjust = 0,
        margin = margin(t = 12)
      ),
      plot.margin = margin(20, 20, 12, 20),
      plot.background = element_rect(fill = colors$palette$bg, color = NA),
      panel.background = element_rect(fill = colors$palette$bg, color = NA)
    )
  )
```

7. Save

```{r}
#| label: save
#| warning: false

### |-  plot image ----  
save_plot_patchwork(
  plot = combined_plot, 
  type = "tidytuesday", 
  year = 2026, 
  week = 15, 
  width  = 10,
  height = 8
  )
```

8. Session Info

TipExpand for Session Info

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] here_1.0.2      binom_1.1-1.1   patchwork_1.3.2 skimr_2.2.2    
 [5] glue_1.8.0      scales_1.4.0    ggrepel_0.9.8   janitor_2.2.1  
 [9] showtext_0.9-7  showtextdb_3.0  sysfonts_0.8.9  ggtext_0.1.2   
[13] lubridate_1.9.5 forcats_1.0.1   stringr_1.6.0   dplyr_1.2.0    
[17] purrr_1.2.1     readr_2.2.0     tidyr_1.3.2     tibble_3.2.1   
[21] ggplot2_4.0.2   tidyverse_2.0.0 pacman_0.5.1   

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1   farver_2.1.2       S7_0.2.0           fastmap_1.2.0     
 [5] gh_1.4.1           digest_0.6.39      timechange_0.4.0   lifecycle_1.0.5   
 [9] rsvg_2.6.2         magrittr_2.0.3     compiler_4.3.1     rlang_1.1.7       
[13] tools_4.3.1        yaml_2.3.12        knitr_1.51         labeling_0.4.3    
[17] htmlwidgets_1.6.4  bit_4.6.0          curl_7.0.0         xml2_1.5.2        
[21] camcorder_0.1.0    repr_1.1.7         RColorBrewer_1.1-3 tidytuesdayR_1.2.1
[25] withr_3.0.2        grid_4.3.1         gitcreds_0.1.2     cli_3.6.5         
[29] rmarkdown_2.30     crayon_1.5.3       generics_0.1.4     otel_0.2.0        
[33] rstudioapi_0.18.0  tzdb_0.5.0         commonmark_2.0.0   parallel_4.3.1    
[37] ggplotify_0.1.3    base64enc_0.1-6    vctrs_0.7.1        yulab.utils_0.2.4 
[41] jsonlite_2.0.0     litedown_0.9       gridGraphics_0.5-1 hms_1.1.4         
[45] bit64_4.6.0-1      systemfonts_1.3.2  magick_2.8.6       gifski_1.32.0-2   
[49] codetools_0.2-19   stringi_1.8.7      gtable_0.3.6       pillar_1.11.1     
[53] rappdirs_0.3.4     htmltools_0.5.9    R6_2.6.1           httr2_1.2.2       
[57] rprojroot_2.1.1    vroom_1.7.0        evaluate_1.0.5     markdown_2.0      
[61] gridtext_0.1.6     snakecase_0.11.1   Rcpp_1.1.1         svglite_2.1.3     
[65] xfun_0.56          fs_1.6.7           pkgconfig_2.0.3   

9. GitHub Repository

TipExpand for GitHub Repo

The complete code for this analysis is available in tt_2026_15.qmd.

For the full repository, click here.

10. References

TipExpand for References

1. Data Source:
   • TidyTuesday 2026 Week 15: Repair Cafes Worldwide

11. Custom Functions Documentation

Note📦 Custom Helper Functions

This analysis uses custom functions from my personal module library for efficiency and consistency across projects.

Functions Used:

• fonts.R: setup_fonts(), get_font_families() - Font management with showtext
• social_icons.R: create_social_caption() - Generates formatted social media captions
• image_utils.R: save_plot() - Consistent plot saving with naming conventions
• base_theme.R: create_base_theme(), extend_weekly_theme(), get_theme_colors() - Custom ggplot2 themes

Why custom functions?
These utilities standardize theming, fonts, and output across all my data visualizations. The core analysis (data tidying and visualization logic) uses only standard tidyverse packages.

Source Code:
View all custom functions → GitHub: R/utils