Visualizing the Heatwave in Portland in 2021

portland_heat.R

library(tidyverse)
library(systemfonts)

## data
## via https://projects.oregonlive.com/weather/temps/
df <- 
  readr::read_csv("https://projects.oregonlive.com/weather/pdx_temps.csv") %>% 
  mutate(yday = lubridate::yday(date),
         year = lubridate::year(date),
         decade = year %/% 10 * 10,
         decade = if_else(decade == 2020, "2020–2021", paste0(decade, "–", decade + 9)),
         new_record = if_else(tmax > 100 & year == 2021, "Last 3 Days in Portland", ""), #New heat records in 2021
         desc = "The Pacific Northwest heat wave shatters temperature records in June 2021.")

months <- tibble(yday = c(1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335),
                 label = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))

## theming
theme_set(theme_minimal(base_size = 18, base_family = "Chivo"))
theme_update(
  panel.grid.major.x = element_blank(),
  panel.grid.minor = element_blank(),
  legend.position = "top",
  legend.title = element_text(size = 14, color = "grey20"),
  legend.text = element_text(size = 12, color = "grey50"),
  plot.title = element_text(size = 22, face = "bold", margin = margin(b = 15)),
  plot.caption = element_text(size = 10, color = "grey50", margin = margin(t = 25)),
  plot.title.position = "plot",
  plot.caption.position = "plot",
  axis.text = element_text(size = 14),
  axis.line.x = element_line(color = "grey20"),
  axis.ticks.x = element_line(color = "grey20"),
  plot.margin = margin(20, 20, 10, 20)
)

## all years
ga <- ggplot(df, aes(yday, tmax, color = year)) +
  ggforce::geom_mark_ellipse(
    aes(fill = new_record, label = new_record, filter = new_record != "", description = desc), 
    alpha = 0, color = "grey20", expand = unit(3, "mm"), con.cap = 0, show.legend = FALSE, label.buffer = unit(40, 'mm'),
    label.fontsize = c(12, 9), label.family = "Chivo"
  ) +
  geom_point(alpha = .5, size = .9) +
  geom_point(data = filter(df, new_record != ""), aes(color = year), size = 1.3) +
  coord_cartesian(clip = "off") +
  scale_x_continuous(breaks = months$yday, labels = months$label, expand = c(.001, .001)) +
  scale_y_continuous(labels = function(x) paste0(x, "°F"), breaks = seq(20, 120, by = 20)) +
  guides(color = guide_colorsteps(barwidth = unit(30, "lines"), barheight = unit(.4, "lines"))) +
  labs(x = NULL, y = NULL, title = "Daily maximum temperatures in Portland, 1938–2021",
       caption = "Data: National Oceanic and Atmospheric Administration via Oregon Live • Graphic: Cédric Scherer")

ga + scale_color_viridis_c(option = "turbo", direction = -1, name = NULL,
                           breaks = c(1938, seq(1950, 2010, by = 10), 2021))

ggsave("portland_heat_all_turbo.pdf", width = 12, height = 10, device = cairo_pdf)

ga + scico::scale_color_scico(palette = "hawaii", name = NULL, direction = -1,
                              breaks = c(1938, seq(1950, 2010, by = 10), 2021))

ggsave("portland_heat_all_hawaii.pdf", width = 12, height = 10, device = cairo_pdf)





## from 1979 onwards
gb <- df %>% 
  filter(year >= 1979) %>% 
  ggplot(aes(yday, tmax, color = year)) +
  ggforce::geom_mark_ellipse(
    aes(fill = new_record, label = new_record, filter = new_record != "", description = desc), 
    alpha = 0, color = "grey20", expand = unit(3, "mm"), con.cap = 0, show.legend = FALSE, label.buffer = unit(40, 'mm'),
    label.fontsize = c(12, 9), label.family = "Chivo"
  ) +
  geom_point(alpha = .5, size = 1.2) +
  geom_point(data = filter(df, new_record != ""), aes(color = year), size = 1.6) +
  coord_cartesian(clip = "off") +
  scale_x_continuous(breaks = months$yday, labels = months$label, expand = c(.001, .001)) +
  scale_y_continuous(labels = function(x) paste0(x, "°F"), breaks = seq(20, 120, by = 20)) +
  guides(color = guide_colorsteps(barwidth = unit(30, "lines"), barheight = unit(.4, "lines"))) +
  labs(x = NULL, y = NULL, title = "Daily maximum temperatures in Portland, 1979–2021",
       caption = "Data: National Oceanic and Atmospheric Administration via Oregon Live • Graphic: Cédric Scherer")

gb + scale_color_viridis_c(option = "turbo", direction = -1, name = NULL,
                           breaks = seq(1980, 2020, by = 10))

ggsave("portland_heat_recent_turbo.pdf", width = 12, height = 10, device = cairo_pdf)

gb + scico::scale_color_scico(palette = "hawaii", name = NULL, direction = -1,
                              breaks = seq(1980, 2020, by = 10))

ggsave("portland_heat_recent_hawaii.pdf", width = 12, height = 10, device = cairo_pdf)



## as facet by decade
gf <- df %>% 
  filter(year >= 1940) %>% 
  ggplot(aes(yday, tmax, color = year)) +
  ggforce::geom_mark_ellipse(
    aes(fill = new_record, label = new_record, filter = new_record != ""), 
    alpha = 0, color = "00000000", label.colour = "grey20", con.colour = "grey20", 
    expand = unit(2, "mm"), con.cap = 0, show.legend = FALSE, 
    label.fontsize = 7, label.family = "Chivo"
  ) +
  geom_point(alpha = .3, size = .6) +
  geom_point(data = filter(df, new_record != ""), aes(color = year)) +
  coord_cartesian(clip = "off") +
  scale_x_continuous(breaks = months$yday, labels = months$label, expand = c(.001, .001)) +
  scale_y_continuous(labels = function(x) paste0(x, "°F"), breaks = seq(20, 120, by = 20)) +
  guides(color = guide_colorsteps(barwidth = unit(30, "lines"), barheight = unit(.4, "lines"))) +
  facet_wrap(~decade) +
  labs(x = NULL, y = NULL, title = "Daily maximum temperatures in Portland, 1940–2021",
       caption = "Data: National Oceanic and Atmospheric Administration via Oregon Live • Graphic: Cédric Scherer")

gf + scale_color_viridis_c(option = "turbo", direction = -1) + 
  theme(legend.position = "none", axis.text = element_text(size = 11))

ggsave("portland_heat_facet_turbo.pdf", width = 13, height = 10, device = cairo_pdf)

gf + scico::scale_color_scico(palette = "hawaii", direction = -1) +
  theme(legend.position = "none", axis.text = element_text(size = 11))

ggsave("portland_heat_facet_hawai.pdf", width = 13, height = 13, device = cairo_pdf)



## as line graph
gs <- df %>% 
  filter(year >= 1979) %>% 
  ggplot(aes(yday, tmax, color = year)) +
  ggforce::geom_mark_ellipse(
    aes(fill = new_record, label = new_record, filter = new_record != "", description = desc), 
    alpha = 0,color = "00000000", label.colour = "grey20", con.colour = "grey20",
    expand = unit(3, "mm"), con.cap = 0, show.legend = FALSE, label.buffer = unit(40, 'mm'),
    label.fontsize = c(12, 9), label.family = "Chivo"
  ) +
  geom_line(alpha = .4) +
  coord_cartesian(clip = "off") +
  scale_x_continuous(breaks = months$yday, labels = months$label, expand = c(.001, .001)) +
  scale_y_continuous(labels = function(x) paste0(x, "°F"), breaks = seq(20, 120, by = 20)) +
  guides(color = guide_colorsteps(barwidth = unit(30, "lines"), barheight = unit(.4, "lines"))) +
  labs(x = NULL, y = NULL, title = "Daily maximum temperatures in Portland, 1979–2021",
       caption = "Data: National Oceanic and Atmospheric Administration via Oregon Live • Graphic: Cédric Scherer")

gs + scale_color_viridis_c(option = "turbo", direction = -1, name = NULL,
                           breaks = seq(1980, 2020, by = 10))

ggsave("portland_heat_lines_turbo.pdf", width = 12, height = 9, device = cairo_pdf)

ggsave("portland_heat_lines_inferno.pdf", width = 12, height = 9, device = cairo_pdf)

gs + rcartocolor::scale_color_carto_c(palette = "Sunset", name = NULL,
                                      breaks = seq(1980, 2020, by = 10))

ggsave("portland_heat_lines_sunset.pdf", width = 12, height = 9, device = cairo_pdf)


## max temperature per year
df %>% 
  group_by(year) %>% 
  filter(tmax == max(tmax)) %>% 
  arrange(yday) %>% 
  slice(1) %>% 
  ggplot(aes(year, tmax)) +
    geom_col(aes(fill = yday, fill = after_scale(colorspace::desaturate(fill, .3)))) +
    scale_x_continuous(breaks = c(1938, seq(1950, 2010, by = 10), 2021), expand = c(.001, .001)) +
    scale_y_continuous(labels = function(x) paste0(x, "°F"), breaks = seq(0, 120, by = 20), 
                       limits = c(0, 120), expand = c(.01, .01)) +
    # scale_fill_viridis_c(option = "rocket", end = .9, direction = -1, 
    #                      name = "Day of the year maximum temperature was reached") +
    scico::scale_fill_scico(palette = "batlow", direction = -1, 
                            name = "Day of the year maximum temperature was reached") +
    coord_cartesian(clip = "off") +
    guides(fill = guide_colorbar(barwidth = unit(30, "lines"), barheight = unit(.4, "lines"),
                                 title.position = "top", title.hjust = .5)) +
    labs(x = NULL, y = NULL, title = "Daily maximum temperatures per year in Portland, 1938–2021",
         caption = "Data: National Oceanic and Atmospheric Administration via Oregon Live • Graphic: Cédric Scherer")

ggsave("portland_heat_bars_batlow.pdf", width = 12, height = 8, device = cairo_pdf)

df %>% 
  mutate(last = if_else(year == 2021, "2021 so far", ""),
         desc = if_else(year == 2021, "Even though 2021 is only 116 days old, we've already seen the hottest days ever since the beginning of weather recording.", "")) %>% 
  group_by(year) %>% 
  ggplot(aes(year, tmax)) +
    ggforce::geom_mark_ellipse(
      aes(fill = year, label = last, filter = last != "", description = desc), 
      alpha = 0, color = "00000000", label.colour = "grey20", con.colour = "grey20", 
      expand = unit(0, "mm"), con.cap = 0, show.legend = FALSE, 
      label.fontsize = c(14, 8), label.family = "Chivo"
    ) +
    stat_summary(geom = "col", fun = max, fill = "white", color = "black", width = 1.01) +
    geom_col(alpha = 1 / 365, fill = "grey20", position = position_dodge(width = 0), width = 1.01) +
    scale_x_continuous(breaks = c(1938, seq(1950, 2020, by = 10)), expand = c(.001, .001)) +
    scale_y_continuous(labels = function(x) paste0(x, "°F"), breaks = seq(0, 120, by = 20), 
                       limits = c(0, 127), expand = c(.01, .01)) +
    coord_cartesian(clip = "off") +
    guides(fill = guide_colorbar(barwidth = unit(30, "lines"), barheight = unit(.4, "lines"),
                                 title.position = "top", title.hjust = .5)) +
    labs(x = NULL, y = NULL, title = "Distribution of daily maximum temperatures per year in Portland, 1938–2021",
         caption = "Data: National Oceanic and Atmospheric Administration via Oregon Live • Graphic: Cédric Scherer")

ggsave("portland_heat_bars_gradient.pdf", width = 12, height = 8, device = cairo_pdf)


## convert PDFs
pdfs <- list.files(here::here(), pattern = "*.pdf", recursive = TRUE)
for(pdf in pdfs) {
  pdftools::pdf_convert(
    pdf = glue::glue("{here::here()}/{pdf}"), 
    filenames = glue::glue("{here::here()}/{str_remove(pdf, '.pdf')}.png"),
    format = "png", dpi = 500
  )
}