Agroclimate Metrics Notebook #2: Cummulative metrics and multi-year summaries

Summary

In this Notebook, we will see how to:

• download 10 years of daily historic observed temperature data for a single location from an interpolated dataset called gridMET hosted on Cal-Adapt

• convert units using the units package

• time slice multi-year data including custom time periods based on months or calendar dates

• create tabular and visual summaries of multi-year metrics

• compute accumulated degree days for specific crops & pests using the degday package

• find the date when a specific accumulated degree day threshold is reached

• interpolate hourly temperatures from the daily min and max

• compute accumulated chill portions using the chillR package

Load Packages

First load the packages we’ll be using:

library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)

Import gridMET data

For this exercise, we’ll work with 10 years (2011-2020) of daily observed historical data from gridMET for a location in the northern San Joaquin Valley. gridMET weather data come as 4km rasters interpolated from weather stations. It is based on PRISM with additional regional reanalysis using climatically aided interpolation to better capture microclimates.

We can get gridMET data from Cal-Adapt, which hosts gridMET data up through 2020. The caladaptR package allows us to import data through the Cal-Adapt API. The first step is to create the API request object (more info). For gridMET data, we need to identify the dataset by its ‘slug’, which we can find by searching the Cal-Adapt API data catalog:

library(caladaptr)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
caladaptr (version 0.6.8)
URL: https://ucanr-igis.github.io/caladaptr
Bug reports: https://github.com/ucanr-igis/caladaptr/issues

## Search the data catalog for gridMET:  
ca_catalog_search("gridmet")

pr_day_gridmet
  name: gridMET daily precipitation historical
  url: https://api.cal-adapt.org/api/series/pr_day_gridmet/
  tres: daily
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: mm
  num_rast: 1
  id: 338
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

pr_year_gridmet
  name: gridMET yearly precipitation historical
  url: https://api.cal-adapt.org/api/series/pr_year_gridmet/
  tres: annual
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: mm
  num_rast: 42
  id: 381
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

tmmn_day_gridmet
  name: gridMET daily minimum temperature historical
  url: https://api.cal-adapt.org/api/series/tmmn_day_gridmet/
  tres: daily
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: K
  num_rast: 1
  id: 919
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

tmmn_year_gridmet
  name: gridMET yearly minimum temperature historical
  url: https://api.cal-adapt.org/api/series/tmmn_year_gridmet/
  tres: annual
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: K
  num_rast: 42
  id: 920
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

tmmx_day_gridmet
  name: gridMET daily maximum temperature historical
  url: https://api.cal-adapt.org/api/series/tmmx_day_gridmet/
  tres: daily
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: K
  num_rast: 1
  id: 921
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

tmmx_year_gridmet
  name: gridMET yearly maximum temperature historical
  url: https://api.cal-adapt.org/api/series/tmmx_year_gridmet/
  tres: annual
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: K
  num_rast: 42
  id: 922
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

Once you find the ‘slugs’ of the dataset(s) of interest, you can construct a API request object:

## Define an object to hold longitude & latitude coordinates
pt1_coords <- c(-121.171, 37.730)

pt1_cap <- ca_loc_pt(coords = pt1_coords) |>
  ca_slug(c("tmmn_day_gridmet", "tmmx_day_gridmet")) |> 
  ca_dates(start = as.Date("2010-10-01"), end = as.Date("2020-09-30"))

pt1_cap

Cal-Adapt API Request
Location(s): 
  x: -121.171
  y: 37.73
Slug(s): tmmn_day_gridmet, tmmx_day_gridmet
Dates: 2010-10-01 to 2020-09-30
 

You can plot a caladaptR API request object to see exactly where it is:

plot(pt1_cap)

To actually retrieve data, you feed the API request object into a function that communicates with the server and retrieves data:

## Uncomment the following to retreive data from the server
# pt1_tbl <- pt1_cap |> ca_getvals_tbl()

pt1_tbl <- readRDS("./data/pt1_tbl.Rds")

glimpse(pt1_tbl)

Rows: 7,306
Columns: 5
$ id   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ slug <chr> "tmmn_day_gridmet", "tmmn_day_gridmet", "tmmn_day_gridmet", "tmmn_day_gridmet", "tmmn_day_gridmet", "tmmn_day_gridme…
$ spag <fct> none, none, none, none, none, none, none, none, none, none, none, none, none, none, none, none, none, none, none, no…
$ dt   <chr> "2010-10-01", "2010-10-02", "2010-10-03", "2010-10-04", "2010-10-05", "2010-10-06", "2010-10-07", "2010-10-08", "201…
$ val  [K] 288.6 [K], 289.2 [K], 286.4 [K], 285.6 [K], 284.1 [K], 283.1 [K], 284.0 [K], 283.3 [K], 284.3 [K], 285.5 [K], 287.7 [K…

head(pt1_tbl)

Convert Units

The climate variable column returned by ca_getvals_tbl() (in this case temperature) is a units objects (i.e., numeric with the units encoded) from the units package. This makes it easy to convert units, for example Kelvin to Fahrenheit, with set_units().

We also need to convert the dt column from character values to Date values, so we can use it for filtering and sorting.

library(units)

udunits database from C:/Users/Andy/AppData/Local/R/win-library/4.2/units/share/udunits/udunits2.xml

pt1_degf_tbl <- pt1_tbl |> 
  mutate(dt_date = as.Date(dt),
         temp_f = set_units(val, degF)) |> 
  select(dt_date, slug, val, temp_f)

head(pt1_degf_tbl)

Plot the raw data:

ggplot(pt1_degf_tbl, mapping = aes(x = dt_date, y = temp_f, group = slug)) + 
  geom_line() + 
  scale_x_date(date_breaks = "1 years", date_labels = "%Y") +
  facet_wrap(~ slug, ncol = 1) 

Note: Because we are using modeled data, gridMET doesn’t contain any missing values. Normally, if we were using weather station data, we would first need to check for missing rows and/or NA values, and then deal with them using substitution and/or interpolation (details).

Time Filtering Multi-Year Datasets

To filter multi-year datasets by month or season, we can use functions from lubridate to pull out date parts. For example to add columns for the month and year:

pt1_month_year_tbl <- pt1_degf_tbl |> 
  mutate(month_num = lubridate::month(dt_date), 
         year = lubridate::year(dt_date)) |> 
  select(dt_date, month_num, year, slug, temp_f)

head(pt1_month_year_tbl)

Group and Summarize

To compute descriptive stats of groups of rows, you can use dplyr’s group_by() followed by summarise(). For example to compute the average daily high temp by year:

pt1_month_year_tbl |> 
  filter(slug == "tmmx_day_gridmet") |> 
  group_by(year) |> 
  summarise(mean_daily_high = mean(temp_f))

McBride and Lacan (2018) computed the average daily high temperature in July as a measure of heat stress for trees. Here’s how that would look for our historical data:

pt1_month_year_tbl |> 
  filter(slug == "tmmx_day_gridmet", month_num == 7) |> 
  group_by(year) |> 
  summarise(mean_daily_high_july = mean(temp_f))

If we wanted the mean daily high in July for this entire time period (not by year), simply remove the group_by() statement:

pt1_month_year_tbl |> 
  filter(slug == "tmmx_day_gridmet", month_num == 7) |> 
  summarise(mean_daily_high_july_2010s = mean(temp_f)) 

Custom Time Slices

Some metrics are computed for custom time periods defined by months or calendar days. You can create columns for custom time periods using vectorized conditional functions such as dplyr::if_else() and dplyr::case_when().

Water Year

The water year starts on Oct. 1st and goes through Sep. 30. It is designated by the calendar year on which it ends.

We can add a column for water year using a mutate expression that includes if_else():

pt1_water_year_tbl <- pt1_month_year_tbl |> 
  mutate(water_year = year + if_else(month_num >= 10, 1, 0))

head(pt1_water_year_tbl)

For more flexibility, lubridate::yday() returns the ‘Julian day’ (1..365) of a date. For example, to pull out April 15 - June 15 each year, we can use:

(jday_start <- lubridate::yday(as.Date("1970-04-15")))

[1] 105

(jday_end <- lubridate::yday(as.Date("1970-06-15")))

[1] 166

pt1_afterbloom_tbl <- pt1_month_year_tbl |> 
  mutate(jday = yday(dt_date)) |> 
  filter(jday >= jday_start, jday <= jday_end) |> 
  select(dt_date, year, jday, slug, temp_f)
  
head(pt1_afterbloom_tbl)

To display the minimum daily temperature during this period as box and whiskers plots:

ggplot(pt1_afterbloom_tbl |> filter(slug == "tmmn_day_gridmet"), 
       mapping = aes(x = year, y = temp_f, group = year)) +
  geom_boxplot() +
  labs(title = "Minimum Daily Temp April 15 - June 15",
       subtitle = "Point 1, 2011 - 2022",
       x = "year",
       y = "temp")

Degree Days

Many phenology events for trees (e.g., blooming) and insects (e.g., egg laying) can be predicted by degree days. Degree days can be thought of as the total amount of warmth, within a usable temperature range, that a plant or insect is exposed to over time. Accumulated degree days make good predictors because plants and insects are cold blooded, hence their development rates are influenced by the ambient temperature.

Some things to know about degree days:

• Degree day are not a real thing that you can measure with a sensor, like temperature or humidity. Rather they are an analytical construct that aims to mirror plant and insect physiology.

• There is no such thing as a ‘universal’ or ‘standard’ degree day. Degree days take into consider the usable range of temperature for a specific species, so they are always in reference to a specific insect, crop, or insect-crop combo.

• Degree days are computed from the daily minimum and maximum temperatures, with additional parameters specific to a crop and/or insect (degree hours are computed from hourly temperature).

• Degree days be can computed in Fahrenheit or Celsius.

• Degree days are not very useful by themselves. You need to use them with a phenology table that predicts when events will take place based on accumulated degree days.

• Phenology tables also tell you when to start counting degree days. This could be a calendar event, or the date of an observation such as when you see eggs in your bug traps.

• There are different formula for computing degree days, including the simple average method, single sine, single triangle, double-sine, and double-triangle. The phenology table will tell you which one to use.

• You can compute degree days using the degday package.

Example: Navel Orangeworm Degree Days

In this example, we’ll use degree days to explore the timing of generations of Navel Orangeworm living in an almond orchard. We begin by looking up which degree day formula to use, and the range of usable temperatures, from UC IPM website, which tells us:

Navel Orangeworm in Almonds

Lower/upper threshold: 55/94°F

Calculation/upper cutoff method: single sine/horizontal

Biofix: The first biofix is the beginning of a consistent increase in egg laying on egg traps. When at least 75% of the egg traps in a given location show increases in the number of eggs on two consecutive monitoring dates, the biofix is the first of those two dates.

Degree Day Events:

• the best time to spray is when 100 NOW-DD have accumulated after biofix

• the next generation of adults can be expected in 1056 NOW-DD after biofix

To compute degree days, we start by putting the daily min and max temperatures in separate columns:

pt1_minmax_tbl <- pt1_tbl |> 
  mutate(dt = as.Date(dt),
         temp_f = set_units(val, degF)) |> 
  pivot_wider(id_cols = dt, names_from = slug, values_from = temp_f) |> 
  rename(tmin = tmmn_day_gridmet, tmax = tmmx_day_gridmet) |> 
  mutate(tmax = if_else(tmax < tmin, tmin, tmax))

head(pt1_minmax_tbl)

Now we can compute NOW Almond degree days:

library(degday)
thresh_low <- 55
thresh_up <- 94

pt1_nowdd_tbl <- pt1_minmax_tbl |> 
  mutate(now_dd = dd_sng_sine(daily_min = tmin, 
                              daily_max = tmax, 
                              thresh_low = thresh_low, 
                              thresh_up = thresh_up))

 - using single sine method

  
head(pt1_nowdd_tbl)

Applying Degree Days: Pest Management

Suppose an almond grower sees a consistent increase in egg laying on egg traps on April 18, 2011 (i.e., the biofix event). The UC IPM website says the best day to spray is when 100 DD have accumulated after the biofix, and the next generation of adults can be expected 1056 DD after biofix. Find the dates for these events.

Step 1 is to filter the dates to begin with the day after biofix, and add a column for accumulated degree days:

pt1_nowdd_2011_tbl <- pt1_nowdd_tbl |> 
  filter(dt > as.Date("2011-04-18"), dt <= as.Date("2011-10-31")) |> 
  mutate(now_dd_acc = cumsum(now_dd)) 

head(pt1_nowdd_2011_tbl)  

Find the date when 100 DD have accumulated:

pt1_nowdd_2011_tbl |> 
  filter(now_dd_acc >= 100) |> 
  slice(1) 

And 1056 DD:

pt1_nowdd_2011_tbl |> 
  filter(now_dd_acc >= 1056) |> 
  slice(1) |> 
  pull(dt)

[1] "2011-07-17"

Applying Degree Days: Estimating the biofix when you don’t have observations

Pathak et al (2021) estimate the emergence of the first flight of Navel Orangeworm as occurring when 300 °F NOW DD have accumulated after January 1st. Compute when this threshold was reached for the historic period.

Step 1 is to remove incomplete years and compute accumulated degree days for each year:

pt1_nowdd_yr_acc_tbl <- pt1_nowdd_tbl |> 
  mutate(year = lubridate::year(dt)) |>    ## add a year column
  filter(year >= 2011) |>                  ## remove 2010 (incomplete year)
  group_by(year) |>                        ## group by years
  mutate(dd_acc_yr = cumsum(now_dd))       ## for each year, compute accumulated DD

head(pt1_nowdd_yr_acc_tbl)

Step 2: on what day each year did we reach 300 DD °F?

pt1_nowdd_yr_acc_tbl |> 
  filter(dd_acc_yr >= 300) |> 
  summarise(first_300dd_date = min(dt), first_300dd_jday = yday(min(dt)))

Adding a column with tomorrow’s temperature

Some degree day formulas (i.e., double-sine and double-triangle) require the next day’s minimum temp to be included. We can add this to our data frame using dplyr::lead().

For example, starting with:

head(pt1_minmax_tbl)

Add the next day’s minimum temperature:

pt1_minmax_nextmin_tbl <- pt1_minmax_tbl |> 
  mutate(tmin_next = lead(tmin, n = 1))

head(pt1_minmax_nextmin_tbl)

Note how lead() treats the last row:

tail(pt1_minmax_nextmin_tbl)

Now we can compute degree days using the double-sine method:

thresh_low <- 55
thresh_up <- 94

pt1_minmax_nextmin_tbl |>
  mutate(now_dd_dblsine = dd_dbl_sine(daily_min = tmin,
                                      daily_max = tmax, 
                                      nextday_min = tmin_next,
                                      thresh_low = thresh_low, 
                                      thresh_up = thresh_up)) |> 
  head()

 - using double sine method

Interpolating Hourly Temps

Some agroclimate metrics require hourly temperatures, such as chill hours and frost exposure (Parker et al, 2021). Ideally you would have hourly temperature data for these metrics, but if not you can interpolate hourly temps from the daily min and max.

One of the best algorithms for interpolating hourly temps comes from Linvill (1990). This method uses an idealized sine curve to describe daytime warming, and a logarithmic decay function for nighttime cooling. The transition between warming and cool is a function of the day length, which in turn is modeled by sunrise and sunset, which in turn is modeled by latitude.

The chillR package has a function make_hourly_temps() hat applies the Linvill method. You need to feed it a data frame that is formatted with specific columns, plus a latitude value (details).

The first step is to make the min and max temps separate columns:

pt1_minmax_tbl <- pt1_tbl |> 
  mutate(dt = as.Date(dt),
         temp_f = set_units(val, degF)) |> 
  pivot_wider(id_cols = dt, names_from = slug, values_from = temp_f) |> 
  rename(tmin = tmmn_day_gridmet, tmax = tmmx_day_gridmet) |> 
  mutate(tmax = if_else(tmax < tmin, tmin, tmax))

head(pt1_minmax_tbl)

Next, we have to add a couple of columns, and change the column names, to match the format expected by chillR::make_hourly_temps() (as described on the help page). This is an example of data wrangling to ‘work backwards’ from what you need to what you’ve got:

pt1_minmax_chillr_tbl <- pt1_minmax_tbl |> 
  mutate(Year = lubridate::year(dt),
         Month = lubridate::month(dt),
         Day = lubridate::day(dt),
         Tmax = as.numeric(tmax),
         Tmin = as.numeric(tmin)) |> 
  select(DATE = dt, Year, Month, Day, Tmax, Tmin)

head(pt1_minmax_chillr_tbl)

Now we can call chillR::make_hourly_temps(), also passing the latitude of our location:

library(chillR)

Attaching package: ‘chillR’

The following object is masked from ‘package:lubridate’:

    leap_year

pt1_coords[2]

[1] 37.73

pt1_hourtemps_wide_tbl <- make_hourly_temps(latitude = pt1_coords[2],
                                            year_file = pt1_minmax_chillr_tbl)
head(pt1_hourtemps_wide_tbl)

make_hourly_temps() gives us wide data. It’s generally easier to work with hourly temperature data in a long format. We can reshape the hourly temps with tidyr::pivot_longer().

pt1_hourtemps_long_tbl <- pt1_hourtemps_wide_tbl |>
  pivot_longer(cols = starts_with("Hour_"),
               names_to = "Hour",
               names_prefix = "Hour_",
               names_transform = list(Hour = as.integer),
               values_to = "temp_f") |>
  mutate(date_hour = lubridate::make_datetime(year = Year, month = Month,
                                              day = Day, hour = Hour,
                                              tz = "America/Los_Angeles")) |>
  select(date_hour, temp_f)

head(pt1_hourtemps_long_tbl)

To see what these hourly temperatures look like, let’s plot one week of them:

pt1_hourtemps_long_tbl |> 
  filter(date_hour >= as.Date("2011-01-01"), date_hour <= as.Date("2011-01-07")) |> 
  ggplot(aes(x = date_hour, y = temp_f)) +
  geom_line(aes(color="red"), show.legend = FALSE)

Chill Portions

Chill portions and chill hours are like degree days, but for cold temperatures. Certain phenology events, like blooming in fruit and nut trees, are correlated with the net amount of cold temperatures the trees been exposed to. This reflect an evolutionary adaptation trees have developed to prevent blooming too early and hence risking damage from a late frost. Different fruit and nut trees have developed different thresholds of chill portions that tell them when its time to ‘wake up’ from their winter dormancy.

We can compute accumulated chill portions by passing a vector of hourly temperature values to chillR::Dynamic_Model(). Note however that Dynamic_Model() expects the temperatures to be in Celsius, so the first step is to add a column with the temperature in °C using the units package:

pt1_hourtemps_degc_tbl <- pt1_hourtemps_long_tbl |> 
  mutate(temp_c = as.numeric(set_units(set_units(temp_f, degF), degC)))
         
head(pt1_hourtemps_degc_tbl)        

Now we can compute accumulated chill portions using Dynamic_Model():

pt1_chillport_tbl <- pt1_hourtemps_degc_tbl |> 
  mutate(chillport_acc = chillR::Dynamic_Model(temp_c, summ = TRUE))

head(pt1_chillport_tbl)

Plot accumulated chill portions for one season:

pt1_chillport_tbl |> 
  filter(date_hour >= as.Date("2010-10-01"), date_hour <= as.Date("2011-07-01")) |> 
  ggplot(aes(x = date_hour, y = chillport_acc)) +
  geom_line() +
  ggtitle("Chill Portions 2010-11")

Challenge Question

New Star Cherries require 54 chill portions to come out of their winter dormancy. Identify the date when this level of chill was reached each year of the observed dataset. For the purposes of the exercise, consider the chill season to go from Nov 1 thru June 30.

Hint: This will look a lot like the question about the date when degree days were reached. Answer

pt1_hourtemps_degc_tbl |> 
  mutate(year = year(date_hour), month = month(date_hour)) |> 
  filter(date_hour >= as.Date("2010-11-01") & (month >= 11 | month <= 6)) |> 
  mutate(chill_season = year + if_else(month >= 11, 1, 0)) |> 
  group_by(chill_season) |> 
  mutate(chillport_acc = chillR::Dynamic_Model(temp_c)) |> 
  filter(chillport_acc >= 54) |> 
  summarise(req_chill_reached = min(date_hour)) |> 
  mutate(month = month(req_chill_reached), day = day(req_chill_reached))

End!

Remember to save the Notebook to generate a HTML version that includes all executed code that you can save for keeps!

---