Scrape and track Brooklyn Botanical Garden cherry blossom blooms

bbg_tree_bloom_2025.csv

brooklyn_botanical_garden_tree_tracker.R

# Code to scrape the Brooklyn Botanical Garden's website and track tree bloom
# https://www.bbg.org/collections/cherries

# Questions:
# - What proportion of trees are blooming?
# - Which flowers are in bloom?
# - Can you normalize the counts to compare trees?
# - Can you remove Kanzan trees? There are a lot of them
# - How do proportions of bloom change over time?
# - Can I make a moving GIF of proportion of trees blooming with the map?


# Setup ----

# Web
library(rvest) # Easily Harvest (Scrape) Web Pages CRAN v1.0.4
library(chromote) # Headless Chrome Web Browser Interface CRAN v0.5.0

# Data wrangling
library(dplyr) # A Grammar of Data Manipulation CRAN v1.1.4
library(tidyr) # Tidy Messy Data CRAN v1.3.1
library(stringr) # Simple, Consistent Wrappers for Common String Operations CRAN v1.5.1

# Utils
library(purrr) # Functional Programming Tools CRAN v1.0.2
library(readr) # Read Rectangular Text Data CRAN v2.1.5
library(forcats) # Tools for Working with Categorical Variables (Factors) CRAN v1.0.0
library(glue) # Interpreted String Literals CRAN v1.7.0
library(forecast) # Forecasting Functions for Time Series and Linear Models CRAN v8.22.0

# Dataviz
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics CRAN v3.5.1

# Image wrangling
library(magick)


# Read current records -----
# Read in records so far
message("Reading in records so far.")
f <- "~/Downloads/bbg_tree_bloom_2025.csv"
if (file.exists(f)) {
  records <- read_csv(f)
}

print(date <- Sys.Date())


# Scrape website ----
# Get live webpage and let JavaScript code load elements
flowers <- read_html_live("https://www.bbg.org/collections/cherries")

message(glue("Current date: {date}"))
message(glue("{html_text(html_elements(flowers, 'figcaption'))} <- Should be same"))

# Parse results ----
message("Get all annotated trees and form them into a nice data frame")
# Inspiration: https://stackoverflow.com/a/34513555/2468369
df_flowers <-
  flowers %>%
  html_nodes("span.location") %>%
  map(html_attrs) %>%
  map_df(~ as.list(.))

# Pre-process data
bloom_lvls <- c("Prebloom", "First Bloom", "Peak Bloom", "Post-Peak Bloom")
new_records <-
  df_flowers %>%
  separate_wider_delim(
    cols = class,
    delim = " ",
    names = c("tree", "location", "bloom", "tooltip")
  ) %>%
  mutate(bloom = case_when(
    bloom == "bloom0" ~ "Prebloom",
    bloom == "bloom1" ~ "First Bloom",
    bloom == "bloom2" ~ "Peak Bloom",
    bloom == "bloom3" ~ "Post-Peak Bloom",
    TRUE ~ "N/A"
  )) %>%
  separate_wider_regex(
    cols = tree,
    patterns = c(tree = "^[a-z_]*_", id = "[0-9]*$")
  ) %>%
  mutate(
    tree = str_remove(tree, pattern = regex("_$")),
    date = date,
    id = as.integer(id),
    bloom = bloom %>% fct_expand(bloom_lvls) %>% fct_relevel(bloom_lvls),
  ) %>%
  select(tree, id, bloom, date)


# Get positions of trees in the Brooklyn Botanical Garden ----
# Only needs to be run once

flower_pos <-
  flowers %>%
  html_nodes("div") %>%
  map(html_attrs) %>%
  map_df(~ as.list(.))

flower_pos %>%
  filter(!is.na(id), !is.na(style), is.na(class)) %>%
  rename(id_full = id) %>%
  mutate(
    top = str_extract(style, "top: ([0-9.]+)%;", group = 1),
    left = str_extract(style, "left: ([0-9.]+)%;", group = 1),
    tree = str_extract(id_full, "^[a-z_]+[a-z]"),
    id = str_extract(id_full, "[0-9]+$")
  ) %>%
  select(id_full, style, tree, id, top, left) %>%
  write_csv(file = "~/Downloads/bbg_flower_positions.csv")


# Colors ----

# Colors for plotting
prebloom <- "#CFC673"
first_bloom <- "#EC008C"
peak_bloom <- "#CD9CC7"
postpeak_bloom <- "#007647"


# Distribution of blooms for today ----

new_records %>%
  count(bloom) %>%
  full_join(data.frame(bloom_lvls), by = c("bloom" = "bloom_lvls")) %>%
  mutate(
    bloom = bloom %>% fct_expand(bloom_lvls) %>% fct_relevel(bloom_lvls),
    n = if_else(is.na(n), 0, n),
    frac = n / sum(n)
  ) %>%
  ggplot(aes(bloom, frac, fill = bloom)) +
  geom_col() +
  theme_minimal() +
  scale_fill_manual(
    values = setNames(c(prebloom, first_bloom, peak_bloom, postpeak_bloom), bloom_lvls),
    guide = "none",
    drop = FALSE
  ) +
  ylim(c(0, 1)) +
  labs(
    x = glue("Bloom status ({date})"),
    y = "Fraction of trees bloomed",
    caption = "Data: https://www.bbg.org/collections/cherries"
  )
ggsave(filename = glue("~/Downloads/bbg_bloom_frac.jpg"), width = 7, height = 5)


# Distribution of flower types, raw ----
new_records %>%
  left_join(records %>% distinct(alt, tree, id), by = join_by(tree, id)) %>%
  count(alt, tree, bloom) %>%
  group_by(alt) %>%
  mutate(alt = factor(alt)) %>%
  ungroup() %>%
  ggplot(aes(n, fct_rev(alt), fill = bloom)) +
  geom_col(show.legend = TRUE) +
  labs(
    x = glue("Number of trees bloomed ({date})"),
    y = "Tree species",
    fill = "Bloom status",
    caption = "Data: https://www.bbg.org/collections/cherries"
  ) +
  theme_minimal() +
  scale_fill_manual(
    values = setNames(c(prebloom, first_bloom, peak_bloom, postpeak_bloom), bloom_lvls),
    drop = FALSE
  ) +
  theme(legend.position = "top", legend.justification = 'left')
ggsave(
  filename = glue("~/Downloads/bbg_bloom_tree_raw.jpg"),
  width = 9,
  height = 5
)


# Distribution of flower types, fraction ----
new_records %>%
  left_join(records %>% distinct(alt, tree, id), by = join_by(tree, id)) %>%
  count(alt, tree, bloom) %>%
  group_by(alt) %>%
  mutate(
    frac = n / sum(n),
    alt = factor(alt)
  ) %>%
  ungroup() %>%
  ggplot(aes(frac, fct_rev(alt), fill = bloom)) +
  geom_col(show.legend = TRUE) +
  labs(
    x = glue("Fraction of trees bloomed ({date})"),
    y = "Tree species",
    fill = "Bloom status",
    caption = "Data: https://www.bbg.org/collections/cherries"
  ) +
  theme_minimal() +
  scale_fill_manual(
    values = setNames(c(prebloom, first_bloom, peak_bloom, postpeak_bloom), bloom_lvls),
    drop = FALSE
  ) +
  theme(legend.position = "top", legend.justification = 'left')
ggsave(
  filename = glue("~/Downloads/bbg_bloom_tree_frac.jpg"),
  width = 9,
  height = 5
)


# Distribution of flower types, raw without Kanzan ----
new_records %>%
  left_join(records %>% distinct(alt, tree, id), by = join_by(tree, id)) %>%
  filter(tree != 'kanzan') %>%
  count(alt, tree, bloom) %>%
  group_by(alt) %>%
  mutate(alt = factor(alt)) %>%
  ungroup() %>%
  ggplot(aes(n, fct_rev(alt), fill = bloom)) +
  geom_col(show.legend = TRUE) +
  scale_x_continuous(breaks = scales::breaks_pretty()) +
  labs(
    x = glue("Number of trees bloomed ({date})"),
    y = "Tree species",
    fill = "Bloom status",
    caption = "Data: https://www.bbg.org/collections/cherries"
  ) +
  theme_minimal() +
  scale_fill_manual(
    values = setNames(c(prebloom, first_bloom, peak_bloom, postpeak_bloom), bloom_lvls),
    drop = FALSE
  ) +
  theme(legend.position = "top", legend.justification = 'left')
ggsave(
  filename = glue("~/Downloads/bbg_bloom_tree_raw_no_kanzan.jpg"),
  width = 9,
  height = 5
)


# Get all the data to see over time ----
new_records <-
  new_records %>%
  select(date, tree, id, bloom) %>%
  left_join(records %>% distinct(alt, tree, id), by = join_by(tree, id)) %>%
  select(date, alt, tree, id, bloom)

# NOTE: On 2025-04-20, Prunus pendula ‘Pendula Plena Rosea’ ID 81 went missing

# Basic checks before joining
dim(new_records) # Should be 152 (sometimes 151)
dim(records)
unique(records$date) # Check that it should be days up until yesterday
records %>% count(alt, tree, id) %>% count(n) # Check counts for days

if (exists("records")) {
  records <- bind_rows(records, new_records)
} else {
  records <- new_records
}
message("Unique dates now:")
print(glue("{unique(records$date)}"))
message(glue("Current date: {date} <-- Should be last date"))


# Bloom fraction over time ----

records %>%
  mutate(
    bloom = bloom %>% fct_expand(bloom_lvls) %>% fct_relevel(bloom_lvls)
  ) %>%
  group_by(date, bloom) %>%
  count(bloom) %>%
  full_join(records %>% expand(date, bloom), by = c("date", "bloom")) %>%
  mutate(n = if_else(is.na(n), 0, n)) %>%
  group_by(date) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x = date, y = frac, color = bloom)) +
  geom_line(linewidth = 1.5, show.legend = TRUE) +
  ylim(c(0, 1)) +
  labs(
    x = "Date",
    y = "Fraction of trees",
    color = "Bloom status",
    caption = "Data: https://www.bbg.org/collections/cherries"
  ) +
  theme_minimal() +
  scale_color_manual(
    values = setNames(c(prebloom, first_bloom, peak_bloom, postpeak_bloom), bloom_lvls),
    drop = FALSE
  ) +
  theme(legend.position = "top", legend.justification = 'top')
ggsave(
  filename = glue("~/Downloads/bbg_bloom_time.jpg"),
  width = 7,
  height = 5
)


# Create GIF of bloom fraction over time ----

dates <- unique(records$date)
for (i in 2:length(dates)) {
  tmp_counts <-
    records %>%
    filter(date <= dates[i]) %>%
    mutate(
      bloom = bloom %>% fct_expand(bloom_lvls) %>% fct_relevel(bloom_lvls)
    ) %>%
    group_by(date, bloom) %>%
    count(bloom) %>%
    ungroup()

  tmp_counts %>%
    full_join(tmp_counts %>% expand(date, bloom), by = c("date", "bloom")) %>%
    mutate(n = if_else(is.na(n), 0, n)) %>%
    group_by(date) %>%
    mutate(frac = n / sum(n)) %>%
    ungroup() %>%
    ggplot(aes(x = date, y = frac, color = bloom)) +
    geom_line(linewidth = 1.5) +
    ylim(c(0, 1)) +
    labs(
      x = "Date",
      y = "Fraction of trees",
      color = "Bloom status",
      caption = "Data: https://www.bbg.org/collections/cherries"
    ) +
    theme_minimal() +
    scale_color_manual(values = c(prebloom, first_bloom, peak_bloom, postpeak_bloom)) +
    theme(legend.position = "top", legend.justification = 'top')

  ggsave(
    filename = glue("~/Downloads/bbg_over_time/bbg_frac_{dates[i]}.jpg"),
    width = 7,
    height = 5
  )
}
list.files(
  path = "~/Downloads/bbg_over_time/",
  pattern = "*.jpg",
  full.names = TRUE
) %>%
  image_read() %>% # reads each path file
  image_join() %>% # joins image
  image_animate(fps = 2) %>% # animates, can opt for number of loops
  image_write("~/Downloads/bbg_over_time/bbg_frac_timelapse.gif")


# Percent of trees changed in the past day ----

# Get basic counts of tree changes per day
# Also need this to join later in the even there is a day without change
tree_change_counts <-
  records %>%
  group_by(tree, id) %>%
  mutate(
    bloom_prev = lag(bloom),
    bloom_diff = if_else(bloom != bloom_prev, 1, 0)
  ) %>%
  ungroup() %>%
  group_by(date) %>%
  count(bloom_diff)

# Join expanded data so that even days with no change are plotted
tree_change_counts %>%
  full_join(
    tree_change_counts %>% expand(date, bloom_diff) %>% na.omit(),
    by = join_by(date, bloom_diff)
  ) %>%
  distinct() %>%
  mutate(n = if_else(is.na(n), 0, n)) %>%
  mutate(total = sum(n)) %>%
  filter(bloom_diff == 1) %>%
  mutate(diff_frac = n / total) %>%
  ggplot(aes(date, n)) +
  geom_col(fill = "#009AB5") +
  # ylim(c(0, 0.25)) +
  # ylim(c(0,25)) +
  labs(
    title = "Number of trees changing their bloom status",
    x = "Date",
    # y = "Fraction of trees changed"
    y = "Number of trees changed",
    caption = "Data: https://www.bbg.org/collections/cherries"
  ) +
  theme_minimal()
  # scale_x_date(date_breaks = "3 days", date_labels = "%b %d")
ggsave(filename = glue("~/Downloads/bbg_bloom_change.jpg"), width = 7, height = 5)


# Forecasting ----

records %>%
  mutate(bloom = fct_relevel(bloom, "Prebloom", "First Bloom", "Peak Bloom", "Post-Peak Bloom")) %>%
  group_by(date, bloom) %>%
  count(bloom) %>%
  group_by(date) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup() %>%
  filter(bloom == "Prebloom") %>%
  pull(frac) %>%
  ts(start = c(2025, 103), frequency = 365) -> ts_data
forecast(ts_data, h = 12)


# Make GIF of BBG map ----

# Source: https://stackoverflow.com/a/56389538/2468369
if (dir.exists("~/Downloads/bbg_map")) {
  # Fast
  list.files(
    path = "~/Downloads/bbg_map/",
    pattern = "*.png",
    full.names = TRUE
  ) %>%
    image_read() %>% # reads each path file
    image_join() %>% # joins image
    image_animate(fps = 2) %>% # animates, can opt for number of loops
    image_write("~/Downloads/bbg_map/bbg_map_timelapse.gif")

  # Slower
  list.files(
    path = "~/Downloads/bbg_map/",
    pattern = "*.png",
    full.names = TRUE
  ) %>%
    image_read() %>% # reads each path file
    image_join() %>% # joins image
    image_animate(fps = 0.5) %>% # animates, can opt for number of loops
    image_write("~/Downloads/bbg_map/bbg_map_timelapse_slow.gif")
}


# Create map + graph combo ----

fraction_imgs <- list.files(path = "~/Downloads/bbg_over_time/", pattern = "*.jpg")
map_imgs <- list.files(path = "~/Downloads/bbg_map/", pattern = "*.png")
for (i in 1:length(fraction_imgs)) {
  c(
    image_read(glue("~/Downloads/bbg_over_time/{fraction_imgs[i]}")),
    image_read(glue("~/Downloads/bbg_map/{map_imgs[i]}"))
  ) %>%
    image_scale("x750") %>%
    image_append() %>%
    image_write(
      path = glue("~/Downloads/bbg_map_time/{str_pad(i, 2, pad = '0')}.jpg"),
      format = "jpg"
    )
}
list.files(
  path = "~/Downloads/bbg_map_time/",
  pattern = "*.jpg",
  full.names = TRUE
) %>%
  image_read() %>% # reads each path file
  image_join() %>% # joins image
  image_animate(fps = 2) %>% # animates, can opt for number of loops
  image_write("~/Downloads/bbg_map_time/bbg_map_time_timelapse.gif")


# Make small inset animation ----

if (!dir.exists("~/Downloads/bbg_map_time_inset")) {
  dir.create("~/Downloads/bbg_map_time_inset")
}
for (i in 1:length(fraction_imgs)) {
  # Make and scale images
  small_map <-
    glue("~/Downloads/bbg_map/{map_imgs[i]}") %>%
    image_read() %>%
    image_background("none") %>%
    image_scale("x175")
  big_plot <-
    glue("~/Downloads/bbg_over_time/{fraction_imgs[i]}") %>%
    image_read() %>%
    image_scale("x800")

  # Put everything together and arrange images together
  image_composite(big_plot, small_map, offset = "+875+80") %>% # +XXX+YYY
    image_write(
      glue("~/Downloads/bbg_map_time_inset/map_inset_{str_pad(i, 2, pad = '0')}.jpg")
    )
}
list.files(
  path = "~/Downloads/bbg_map_time_inset/",
  pattern = "*.jpg",
  full.names = TRUE
) %>%
  image_read() %>% # reads each path file
  image_join() %>% # joins image
  image_animate(fps = 2) %>% # animates, can opt for number of loops
  image_write("~/Downloads/bbg_map_time_inset/bbg_map_time_inset_timelapse.gif")


# Write out results ----

write_csv(x = records, file = f)