Exercise 4: Compute Degree Days for the Entire Season

Summary

In this Notebook, we import the weather data tables we imported in Exercise 2 & 3 for three time periods:

  1. Recent past (Exercise 2)

  2. Short-term forecast (Exercise 3)

  3. Seasonal forecast (Exercise 3):
    Open-Meteo’s Seasonal Forecast service
    Long-term daily averages from gridMet

We will then:

  1. append the tables together to get a time series for the entire season
  2. compute degrees for the entire season
  3. make predictions about the best time to till to disrupt the life cycle of the three-cornered alfalfa hopper (exercise 1)
  4. compare the predictions for the two proxies for the seasonal forecast

1 Import the individual weather data files

Begin by importing the weather data tables we created in exercises 2 and 3:

1.1 Recent past:

data_dir <- here::here("exercises/data")
dir.exists(data_dir)

## For the recent past data for CIMIS Station 077, we will use a saved copy of the station
## data download using the CIMIS API (see import-data-cimis-api.qmd), 
## rather than the copy we imported thru Syntopic,
## We do this because there is ~6 week gap in the Synoptic data during Jan & Feb 2024 
## when CIMIS data were not ingested into Synoptic.

stn_rctpast_dlytas_tbl <- readRDS(file.path(data_dir, "stn_rctpast_cimis_dlytas_tbl.Rds"))
stn_rctpast_dlytas_tbl |> head()
stn_rctpast_dlytas_tbl$date |> range()
[1] TRUE
[1] "2024-01-01" "2024-05-21"


1.2 Short-term forecast

stn_shrtfrcst_dlytas_tbl <- readRDS(file.path(data_dir, "stn_shrtfrcst_dlytas_tbl.Rds"))
stn_shrtfrcst_dlytas_tbl |> head()
stn_shrtfrcst_dlytas_tbl$date |> range()
[1] "2024-05-22" "2024-05-29"


1.3 Seasonal Forecast

stn_ssnfcast_dlytas_tbl <- readRDS(file.path(data_dir, "stn_ssnfcast_dlytas_tbl.Rds"))
stn_ssnfcast_dlytas_tbl |> head()
stn_ssnfcast_dlytas_tbl$date |> range()
[1] "2024-05-30" "2024-08-28"


1.4 Long-term daily averages

stn_histavg_dlytas_tbl <- readRDS(file.path(data_dir, "stn_histavg_dlytas_tbl.Rds"))
stn_histavg_dlytas_tbl |> head()
stn_histavg_dlytas_tbl$date |> range()
[1] "2024-05-30" "2024-08-28"


1.5 Append the tables

First, load packages:

library(dplyr) |> suppressPackageStartupMessages()

# Set preferences for functions with common names 
library(conflicted)
conflict_prefer("filter", "dplyr", quiet = TRUE)
conflict_prefer("count", "dplyr", quiet = TRUE)
conflict_prefer("select", "dplyr", quiet = TRUE)
conflict_prefer("arrange", "dplyr", quiet = TRUE)


Append the tables:

Note: after combining the tables we convert the period column to a factor with manually specified levels. This is how we control the order of the three periods in the legend in ggplot.

stn_all_dlytas_ssnfcast_tbl <- stn_rctpast_dlytas_tbl |> 
  bind_rows(stn_shrtfrcst_dlytas_tbl) |> 
  bind_rows(stn_ssnfcast_dlytas_tbl) |> 
  arrange(date) |> 
  mutate(period = factor(period, levels = c("rp", "stf", "seasn")))


1.6 Plot time series

First, we plot the Open-Meteo Seasonal Forecast as our best guess for the rest of the season:

library(ggplot2)

ggplot(stn_all_dlytas_ssnfcast_tbl, mapping = aes(x = date, y = tasmax, col = period)) +
  geom_line() +
  geom_line(mapping = aes(y = tasmin)) +
  ylab("air temperature") +
  xlab(NULL) +
  labs(title = "Minimum and Maximum Daily Temperature",
       subtitle = "CIMIS Station #077, Oakville. 2024",
       caption = "Sources:\nRecent past: CIMIS API\nShort-term forecast: Open-Meteo Weather Forecast\nSeasonal forecast: Open-Meteo Seasonal Forecast") +
  theme(plot.caption = element_text(hjust = 0))


Next, do the same using the historic daily average as the proxy for the seasonal forecast:

stn_all_dlytas_histavg_tbl <- stn_rctpast_dlytas_tbl |> 
  bind_rows(stn_shrtfrcst_dlytas_tbl) |> 
  bind_rows(stn_histavg_dlytas_tbl) |> 
  arrange(date) |> 
  mutate(period = factor(period, levels = c("rp", "stf", "histavg")))

ggplot(stn_all_dlytas_histavg_tbl, mapping = aes(x = date, y = tasmax, col = period)) +
  geom_line() +
  geom_line(mapping = aes(y = tasmin)) +
  ylab("air temperature") +
  xlab(NULL) +
  labs(title = "Minimum and Maximum Daily Temperature",
       subtitle = "CIMIS Station #077, Oakville. 2024",
       caption = "Sources:\nRecent past: CIMIS API\nShort-term forecast: Open-Meteo Weather Forecast\nSeasonal forecast: Long-term daily average, gridMet 1990-2020") +
  theme(plot.caption = element_text(hjust = 0))


2 Compute Degree Days

library(degday)

oakville_dd_ssnfast_tbl <- stn_all_dlytas_ssnfcast_tbl |> 
  mutate(dd_f = dd_sng_sine(daily_min = tasmin, 
                          daily_max = tasmax, 
                          thresh_low = 32,
                          thresh_up = 999),
         acc_dd_f = cumsum(dd_f))
 - using single sine method
oakville_dd_histavg_tbl <- stn_all_dlytas_histavg_tbl |> 
  mutate(dd_f = dd_sng_sine(daily_min = tasmin, 
                          daily_max = tasmax, 
                          thresh_low = 32,
                          thresh_up = 999),
         acc_dd_f = cumsum(dd_f))
 - using single sine method
oakville_dd_histavg_tbl


Let’s visualize the difference between accumulated degree days the seasonal forecast vs historic daily averages make in the accumulation of degree days?

oakville_dd_comb_tbl <- oakville_dd_ssnfast_tbl |> 
  bind_rows(oakville_dd_histavg_tbl |> filter(period == "histavg"))

ggplot(oakville_dd_comb_tbl, mapping = aes(x = date, y = acc_dd_f, colour = period)) +
  geom_line() +
  ylab("degree days (F)") +
  xlab(NULL) +
  labs(title = "Accumulated Degree Days",
       subtitle = "CIMIS Station #077, Oakville. 2024",
       caption = "Sources:\nRecent past: CIMIS API\nShort-term forecast: Open-Meteo Weather Forecast\nSeasonal forecast: Long-term daily average, gridMet 1990-2020; Open-Meteo Seasonal Forecast") +
  theme(plot.caption = element_text(hjust = 0))

2.1 Compare the predictions

Let’s compare the predictions regarding the best window to till the soil:

till_yn <- function(x) {x >= 2358 & x <= 2817}

## OM Seasonal Forecast:
oakville_dd_ssnfast_tbl |> filter(till_yn(acc_dd_f)) |> pull(date) |> range()
[1] "2024-04-26" "2024-05-11"
## Long-term historic daily averages
oakville_dd_histavg_tbl |> filter(till_yn(acc_dd_f)) |> pull(date) |> range()
[1] "2024-04-26" "2024-05-11"


Not surprisingly - they’re the same!

3 CHALLENGE

Compare the predictions by the seasonal forecast and long-term daily averages regarding when adults are expected (3137-3447 accumulated degree days Fahrenheit). [Solution]

## Your answer here

4 END!

We now have all the code we need to make a Shiny app!