Calculating normal grazing periods using historical nClimGrid-daily data

In several US Department of Agriculture (USDA) Farm Service Agency (FSA) drought and disaster relief programs — including the Noninsured Crop Disaster Assistance Program (NAP), Livestock Forage Disaster Program (LFP), and Emergency Assistance for Livestock, Honey Bees, and Farm-raised Fish (ELAP) programs — program eligibility is determined in part by whether the qualifying disaster (and resulting loss) occurs during the “normal grazing period” for a particular type of pasture or grazing crop. State FSA offices have broad discretion to define “normal grazing period” for their regions based on local conditions.

In April 2018, USDA FSA issued Notice NAP-190, Review of Normal Grazing Periods and Carrying Capacities for NAP, LFP, and ELAP, which sought to provide a more consistent definition of growing, coverage, and grazing periods across these programs. According to NAP-190, “the normal grazing periods established for each pasture or grazing crop type for LFP are to be consistent with grazing periods established for NAP purposes according to 1-LDAP (Rev. 1), subparagraph 411 M. DAFP has determined the grazing periods are to be consistent with ELAP as well.” Also, “for program integrity purposes, growing periods, coverage periods, and grazing periods are analogous to each other and need to be consistent for each crop and type intended for forage or grazing” (Farm Service Agency 2018, 1).

NAP-190 defines a process for proposing normal grazing periods by county across the US based on daily temperature recordings (Farm Service Agency 2018, 1). The USDA Deputy Administrator for Farm Programs office initially calculated normal grazing periods across the US. These county-level grazing period definitions were then shared with the state FSA offices, which were invited to review and recommend amendments. Amended grazing periods were returned to the National Office for review and approval.

Here, we use the procedures outlined in NAP-190 (Farm Service Agency 2018, 4–7) to re-calculate normal grazing periods for all counties in the US, using the latest thirty-year climate normals (1991–2020).

Guidelines for Establishing of Normal Grazing Periods and Carrying Capacities

NAP-190 descripbes a probabilistic, temperature-based definition of the start and end dates of the grazing period (Farm Service Agency 2018, 4):

Livestock receive daily nutrients and obtain net energy requirements during the spring and early summer months of the forage growth period without supplemental feed. DAFP has reviewed available information and determined:

• full season and/or warm season grazed forage crops’ types plant growth historically:
  • begins on or after the last date of an 80 percent probability that a 28 degree Fahrenheit temperature (freeze date) occurs in a county, known as the beginning date of the grazing period
  • ends on or before the first date of an 80 percent probability of a 28 degree Fahrenheit temperature (freeze date) occurring in a county, known as the normal harvest date (grazing period end date).
• cool season grazed forage crops’ plant growth development for nutritional value historically begins on or after the first date the minimum average temperature is approximately 50 degrees Fahrenheit and transitions to a semi-dormant stage when temperatures reach an average maximum of 90 degrees Fahrenheit. This date is the normal harvest date (grazing period end date) for cool season grasses.

The definition used for cool season crops is unfortunately ambiguous. The start date used the phrasing “minimum average temperature” (suggesting the minimum of averages, presumably across the county), while the end date uses “average maximum” (suggesting the average maximum temperature that occurs in a county, analogous to the definition for full/warm season crops). Here, we apply the phrasing for the end date to both the start and end dates, and use the mode of the distribution instead of the average (due to skewness) — we take the climatological modes of minimum and maximum temperatures that occur in each county, and set start and end dates using the 50º and 90º thresholds, as described.

NAP-190 offers further guidance on the determination of warm and cool season grazing periods (Farm Service Agency 2018, 3):

STC’s may seek approval from DAFP to establish both warm and cool season grazing periods. As of the date of this notice, only Arkansas, Mississippi, and Oklahoma have been approved by DAFP for establishing warm season and cool season grazing periods. States where both warm season and cool season grazing periods may also apply include: Alabama, Florida, Georgia, Louisiana, North Carolina, South Carolina, and Texas. Before a State can establish both warm and cool season grazing periods, DAFP approval must be obtained. Note: The establishment of both warm and cool season grazing periods:

• does not require the grazing periods to be equal in length of time
• combined will not exceed 12 months
• cannot result in an overlap of the grazing periods.

Dates are customarily rounded to the 1st or 15th of the month (Farm Service Agency 2018, 7):

The proposed grazing period, established as discussed in subparagraph 2 B, has the:

• beginning date of the grazing period (start date) rounded to the 1st or 15th of the month
• normal harvest date (end date) rounded to the 15th or last day of the month.

Establishing normal grazing periods using a 1991–2020 climatology

“Climate normals” are statistical descriptions of typical climate conditions. Normals are often calculated using past weather data, and can be predictive of conditions in the near future. Conventional practice since the 1930s is to calculate climate normals using 30 years of data, updated decadally; the last decadal update of the NOAA climate normals is for the 1991–2020 period.

Traditionally, temperature normals (for daily minimum and maximum temperature) are represented by the average value of thirty years of data. However, it is useful to represent the expected range of variability as well, and the methods described in NAP-190 require probabilistic estimates of temperatures to define grazing period start and end dates. Here, we develop a method for estimating the parameters of statistical distributions describing daily temperature extremes. Our methods are similar in principle to those described by Rennie, Durre, and Palecki (2022), except that we derive daily distribution parameters which we smooth annually, and from which we calculate values such as the mode and exceedance probabilities.

Method summary and example

In order to establish normal grazing periods for each county, we must first define normal climatologies for the US. The basic steps of the analysis are:

1. Find the minimum and maximum daily temperature in each county
2. Estimate Gumbel distribution parameters using L-moments
3. Smooth annual series of daily Gumbel distribution parameters annually using a cyclic smoother
4. Identify and round first and last days of grazing periods

To illustrate the method, let’s derive annual daily climatologies for minimum and maximum temperature for Missoula County, Montana. All analyses performed below are scripted in the R scientific programming language for reproducibility (R Core Team 2023). The data processing described here relies heavily on functions in the dplyr (Wickham et al. 2023), tidyr (Wickham, Vaughan, and Girlich 2023), and magrittr (Bache and Wickham 2022) packages for R.

Find the minimum daily temperature in each county

We start with the daily timeseries of minimum and maximum temperatures from 1991–2020, extracted from the nClimGrid-daily dataset.

## NCLIMGRID Grazing Periods, 1991-2020
dir.create("data-raw/nclimgrid",
           recursive = TRUE,
           showWarnings = FALSE
)

ncores <- max(1, future::availableCores() - 2)
future::plan(multisession,
             workers = ncores)


get_nclimgrid <- function(x, y){
  if(file.exists(y))
    return(NULL)
  
  httr::GET(url = x, httr::write_disk(path = y))
}

nclimgrid_tmin <-
  tidyr::crossing(year = 1991:2020,
                  month = stringr::str_pad(1:12, width = 2, pad = 0)) %>%
  dplyr::mutate(url = paste0("https://noaa-nclimgrid-daily-pds.s3.amazonaws.com/beta/by-month/",year,"/",month,"/tmin-",year,month,"-grd-scaled.nc"),
                file = paste0("data-raw/nclimgrid/tmin-",year,month,"-grd-scaled.nc"),
                month = as.integer(month))

# Download data
nclimgrid_tmin %$%
  furrr::future_walk2(.x = url, .y = file, .f = ~get_nclimgrid(x = .x, y = .y))

nclimgrid_tmax <-
  tidyr::crossing(year = 1991:2020,
                  month = stringr::str_pad(1:12, width = 2, pad = 0)) %>%
  dplyr::mutate(url = paste0("https://noaa-nclimgrid-daily-pds.s3.amazonaws.com/beta/by-month/",year,"/",month,"/tmax-",year,month,"-grd-scaled.nc"),
                file = paste0("data-raw/nclimgrid/tmax-",year,month,"-grd-scaled.nc"),
                month = as.integer(month))

# Download data
nclimgrid_tmax %$%
  furrr::future_walk2(.x = url, .y = file, .f = ~get_nclimgrid(x = .x, y = .y))

missoula_county <- 
  sf::read_sf("data-raw/counties.fgb") %>%
  dplyr::filter(county == "Missoula")

set_year <- function(x, y){
  lubridate::year(x) <- y
  x
}

ncores <- max(1, future::availableCores() - 2)
future::plan(multisession,
             workers = ncores)

missoula_tmin <-
  list.files("data-raw/nclimgrid/",
             pattern = "tmin",
             full.names = TRUE) %>%
  {
    split(., f = ggplot2::cut_number(1:length(.), n = ncores))
  } %>%
  furrr::future_map_dfr(
    .f = 
      function(x){
        terra::rast(x) %>%
          magrittr::set_names(., terra::time(.)) %>%
          {
            tibble::tibble(
              date = lubridate::as_date(names(.)),
              yday = lubridate::as_date(names(.)) %>%
                # We use 1993-1994 as generic years throughout this 
                # analysis because they are not leap years.
                set_year(1993) %>%
                lubridate::yday(),
              tmin = exactextractr::exact_extract(.,
                                                  missoula_county,
                                                  fun = "min") %>%
                as.numeric()
            )
          } %>%
          dplyr::filter(!is.na(yday))
      },
    .options = furrr::furrr_options(seed = NULL))

missoula_tmax <-
  list.files("data-raw/nclimgrid/",
             pattern = "tmax",
             full.names = TRUE) %>%
  {
    split(., f = ggplot2::cut_number(1:length(.), n = ncores))
  } %>%
  furrr::future_map_dfr(
    .f = 
      function(x){
        terra::rast(x) %>%
          magrittr::set_names(., terra::time(.)) %>%
          {
            tibble::tibble(
              date = lubridate::as_date(names(.)),
              yday = lubridate::as_date(names(.)) %>%
                set_year(1993) %>%
                lubridate::yday(),
              tmax = exactextractr::exact_extract(.,
                                                  missoula_county,
                                                  fun = "max") %>%
                as.numeric()
            )
          } %>%
          dplyr::filter(!is.na(yday))
      },
    .options = furrr::furrr_options(seed = NULL))

missoula <- 
  dplyr::full_join(missoula_tmin, missoula_tmax)%>%
  na.omit() %>%
  dplyr::mutate(
    dplyr::across(
      tmin:tmax, 
      ~(
        units::set_units(.x, "degC") %>%
          units::set_units("degF") %>%
          units::drop_units()
      )
    )
  )

ggplot2::ggplot(
  data = missoula %>%
    tidyr::pivot_longer(tmin:tmax,
                        names_to = "type"),
  mapping = aes(x = date, y = value, color = type)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::scale_y_continuous(limits = c(-45, 105),
                              breaks = seq(-20,100,20),
                              name = "Temperature (ºF)",
                              expand = expansion(0,0)
  ) +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  scale_x_date(date_labels = "%Y"
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank()
  )

Daily minimum and maximum temperatures for Missoula County, Montana, 1991–2020. Data are from the NOAA nClimGrid-daily dataset.

Each day of the year is represented by a series of minimum and maximum temperatures that occur in the county, such as these for October 23 of each year:

missoula_oct23 <- 
  missoula %>%
  dplyr::filter(lubridate::month(date) == 10,
                lubridate::day(date) == 23)

ggplot2::ggplot(
  data = missoula_oct23 %>%
    tidyr::pivot_longer(tmin:tmax,
                        names_to = "type"),
  mapping = aes(x = lubridate::year(date), y = value, color = type)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::scale_y_continuous(limits = c(0, 100),
                              breaks = seq(0,100,20),
                              name = "Temperature (ºF)",
                              expand = expansion(0,0)
  ) +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  scale_x_continuous(name = "Year"
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        # panel.grid = element_blank()
  )

October 23 daily minimum and maximum temperatures for Missoula County, Montana, 1991–2020. Daily low temperatures range between 20ºF and 40ºF; highs between 40ºF and 75ºF.

Estimate Gumbel distribution parameters using L-moments

To calculate the normals, we fit a probability distribution to each of these series. Gaussian (normal) distributions are traditionally used for average temperature, but minimum and maximum temperature distributions are often skewed. We therefore fit Gumbel distributions to the data, which are appropriate for modeling extreme values such as temperature minima and maxima. Here, we fit Gumbel distributions to daily minimum and maximum temperatures in each county using L-moments. L-moments are measures of the location, scale, and shape of probability distributions or data samples (J. Hosking 1990; J. Hosking and Wallis 1997), and sample L-moments provide a robust estimate of distribution parameters. The first L-moment, \(\lambda_1\), is the sample mean.

Fitting for our data from October 23 yields these distributions:

missoula_oct23_dists <-
  missoula_oct23 %>%
  dplyr::filter(lubridate::month(date) == 10,
                lubridate::day(date) == 23) %>%
  dplyr::summarise(
    tmin = list(
      
      lmomco::pargum(
        lmomco::lmoms(-1 * tmin)
      )
      
    ),
    tmax = list(
      
      lmomco::pargum(
        lmomco::lmoms(tmax)
      )
      
    )
  )

ggplot2::ggplot() +
  stat_function(
    fun = ~lmomco::pdfgum(.x, 
                          para = missoula_oct23_dists$tmax[[1]]), 
    xlim = lmomco::quagum(f = c(0.2,0.8),
                          para = missoula_oct23_dists$tmax[[1]]),
    fill = "red",
    alpha = 0.5,
    geom = "area") +
  stat_function(
    fun = ~lmomco::pdfgum(.x, 
                          para = missoula_oct23_dists$tmax[[1]]), 
    xlim = c(-10,100),
    color = "red",
    geom = "line",
    linewidth = 1,
    n = 1000) +
  stat_function(
    fun = ~lmomco::pdfgum(-1 * .x, 
                          para = missoula_oct23_dists$tmin[[1]]), 
    xlim = -1 * lmomco::quagum(f = c(0.2,0.8),
                               para = missoula_oct23_dists$tmin[[1]]),
    fill = "dodgerblue",
    alpha = 0.5,
    geom = "area") +
  stat_function(
    fun = ~lmomco::pdfgum(-1  * .x, 
                          para = missoula_oct23_dists$tmin[[1]]), 
    xlim = c(-10,100),
    color = "dodgerblue",
    geom = "line",
    linewidth = 1,
    n = 1000) +
  ggplot2::geom_point(data = missoula_oct23 %>%
                        tidyr::pivot_longer(tmin:tmax,
                                            names_to = "type"),
                      mapping = aes(x = value,
                                    y = 0,
                                    color = type)) +
  ggplot2::scale_x_continuous(limits = c(-10, 100),
                              name = "Temperature (ºF)",
                              expand = expansion(0,0)) +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  ggplot2::scale_y_continuous(limits = c(0,NA),
                              expand = expansion(c(0,0.05),c(0,0))) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  ) +
  ggplot2::coord_flip()

Gumbel distributions fit to 30-year minimum and maximum temperature for October 23 of each year. The shaded regions show the 20th to 80th percentiles of each distribution.

We’ve also shaded the 20th to 80th percentiles, 60% of the data fall within these boundaries. Notice how the distributions for minimum temperature on October 23 is narrower than that for maximum temperature.

We can perform a similar process for every day of the year—fitting Gumbel distributions to the observations for each day. If we plot just the raw data by day, we start to get a sense of the variability:

background <-
  tibble::tibble(
    start = seq(lubridate::as_date("1993-01-01"),
                lubridate::as_date("1993-11-01"),
                "2 months"),
    end = seq(lubridate::as_date("1993-02-01"),
              lubridate::as_date("1993-12-01"),
              "2 months"),
    yday = lubridate::yday(start))


ggplot2::ggplot(
  data = missoula %>%
    tidyr::pivot_longer(tmin:tmax,
                        names_to = "type"),
  mapping = aes(x = yday + lubridate::as_date("1992-12-31"))) +
  ggplot2::geom_rect(data = background,
                     mapping = aes(xmin = start,
                                   xmax = end),
                     ymin = -Inf,
                     ymax = Inf,
                     fill = "gray95") +
  ggplot2::geom_point(
    data = missoula %>%
      tidyr::pivot_longer(tmin:tmax,
                          names_to = "type"),
    mapping = ggplot2::aes(y = value,
                           color = type,
                           group = type), 
    size = 0.1) +
  ggplot2::scale_y_continuous(limits = c(-45, 110),
                              breaks = seq(-20,100,20),
                              name = "Temperature (ºF)",
                              expand = expansion(0,0)
  ) +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  scale_x_date(date_breaks = "month",
               date_labels = "%b",
               expand = expansion(0,0)
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        panel.grid = element_blank()
  )

Daily minimum and maximum temperatures for Missoula County, Montana, 1991–2020, organized by day of the year. Every day has thirty points.

Notice how the spread of distribution of maximum temperature seems to be relatively stable through time, while minimum temperature is narrower in the summer and wider in the winter.

Now, let’s do the same plot, but add the daily 20% to 80% distribution range.

missoula_dists <-
  missoula %>%
  dplyr::group_by(yday) %>%
  dplyr::summarise(
    tmin = list(
      
      lmom::pelgum(
        lmom::.samlmu(-1 * tmin)
      )
      
    ),
    tmax = list(
      
      lmom::pelgum(
        lmom::.samlmu(tmax)
      )
      
    ),
    .groups = "drop"
  )

missoula_quants <-
  missoula_dists %>%
  dplyr::rowwise() %>%
  dplyr::mutate(
    tmin = 
      list(
        -1 * lmom::quagum(tmin, f = c(0.8,0.5,0.2)) %>%
          magrittr::set_names(c("20%","50%","80%")) %>%
          as.list() %>%
          tibble::as_tibble()
      ),
    tmax = 
      list(
        lmom::quagum(tmax, f = c(0.2,0.5,0.8)) %>%
          magrittr::set_names(c("20%","50%","80%")) %>%
          as.list() %>%
          tibble::as_tibble()
      )
  ) %>%
  tidyr::pivot_longer(tmin:tmax,
                      names_to = "type") %>%
  tidyr::unnest(value)

ggplot2::ggplot(
  data = missoula %>%
    tidyr::pivot_longer(tmin:tmax,
                        names_to = "type"),
  mapping = aes(x = yday + lubridate::as_date("1992-12-31"))) +
  ggplot2::geom_rect(data = background,
                     mapping = aes(xmin = start,
                                   xmax = end),
                     ymin = -Inf,
                     ymax = Inf,
                     fill = "gray95") +
  ggplot2::geom_point(
    data = missoula %>%
      tidyr::pivot_longer(tmin:tmax,
                          names_to = "type"),
    mapping = ggplot2::aes(y = value,
                           color = type,
                           group = type), 
    size = 0.1) +
  geom_ribbon(
    data = 
      missoula_quants, 
    mapping = ggplot2::aes(ymin = `20%`, 
                           ymax = `80%`,
                           group = type), 
    fill = "black",
    color = "black",
    alpha = 0.5) +
  ggplot2::geom_line(
    data = missoula_quants,
    mapping = ggplot2::aes(y = `50%`,
                           group = type),
    color = "black",
    linewidth = 1) +
  ggplot2::scale_y_continuous(limits = c(-45, 110),
                              breaks = seq(-20,100,20),
                              name = "Temperature (ºF)",
                              expand = expansion(0,0)
  ) +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  scale_x_date(date_breaks = "month",
               date_labels = "%b",
               expand = expansion(0,0)
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        panel.grid = element_blank()
  )

Daily minimum and maximum temperatures for Missoula County, Montana, with 20th, 50th, and 80th percentiles plotted from the fitted Gumbel distributions.

Smooth annual series of daily Gumbel distribution parameters annually using a cyclic smoother

These jagged representations of the distributions are typical of daily normals; they ought to be smoothed to better represent the continuous nature of expected temperatures from day to day. We apply a cyclical generalized additive model to the parameters of the daily Gumbel distributions in order to estimate smoothed parameters. We used 12 “knots”, or hinge points, in the smoothing function in order to allow for monthly variability. Doing so yields smooth daily climate normals that capture variability in both the magnitude and spread of temperature through the year.

cyclic_smooth <-
  function(y){
    if(is.na(y[1]))
      return(rep(NA,length(y)))
    
    as.numeric(
      predict(
        mgcv::gam(
          formula = y ~ s(x, k = 12, bs = "cc"),
          data = data.frame(x = 1:365,
                            y = y)
        ),
        type = "response")
    )
  }


smooth_dist <-
  function(x){
    
    x %>%
      dplyr::bind_rows() %>%
      dplyr::mutate(dplyr::across(dplyr::everything(), cyclic_smooth)) %>%
      dplyr::rowwise() %>%
      dplyr::group_split() %>%
      purrr::map(unlist)
  }


missoula_quants_smooth <-
  missoula_dists %>%
  dplyr::mutate(dplyr::across(tmin:tmax, smooth_dist)) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(
    tmin = 
      list(
        -1 * lmom::quagum(tmin, f = c(0.8,0.5,0.2)) %>%
          magrittr::set_names(c("20%","50%","80%")) %>%
          as.list() %>%
          tibble::as_tibble()
      ),
    tmax = 
      list(
        lmom::quagum(tmax, f = c(0.2,0.5,0.8)) %>%
          magrittr::set_names(c("20%","50%","80%")) %>%
          as.list() %>%
          tibble::as_tibble()
      )
  ) %>%
  tidyr::pivot_longer(tmin:tmax,
                      names_to = "type") %>%
  tidyr::unnest(value)




ggplot2::ggplot(
  data = missoula%>%
    tidyr::pivot_longer(tmin:tmax,
                        names_to = "type"),
  mapping = aes(x = yday + lubridate::as_date("1992-12-31"))) +
  ggplot2::geom_rect(data = background,
                     mapping = aes(xmin = start,
                                   xmax = end),
                     ymin = -Inf,
                     ymax = Inf,
                     fill = "gray95") +
  ggplot2::geom_point(
    data = missoula%>%
      tidyr::pivot_longer(tmin:tmax,
                          names_to = "type"),
    mapping = ggplot2::aes(y = value,
                           color = type,
                           group = type), 
    size = 0.1) +
  geom_ribbon(
    data = missoula_quants_smooth, 
    mapping = ggplot2::aes(ymin = `20%`, ymax = `80%`,
                           # color = type,
                           group = type), 
    fill = "black",
    color = "black",
    alpha = 0.5) +
  ggplot2::geom_line(
    data = missoula_quants_smooth,
    mapping = ggplot2::aes(y = `50%`,
                           # color = type,
                           group = type),
    color = "black",
    linewidth = 1) +
  ggplot2::scale_y_continuous(limits = c(-45, 110),
                              breaks = seq(-20,100,20),
                              name = "Temperature (ºF)",
                              expand = expansion(0,0)
  ) +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  scale_x_date(date_breaks = "month",
               date_labels = "%b",
               expand = expansion(0,0)
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        panel.grid = element_blank()
  )

Daily minimum and maximum temperatures for Missoula County, Montana, with 20th, 50th, and 80th percentiles plotted from the smoothed Gumbel distributions.

Identify and round first and last days of grazing periods

The guidance for establishing normal grazing periods (in areas without cool-season forage) states that the grazing period “begins on or after the last date of an 80 percent probability that a 28 degree Fahrenheit temperature (freeze date) occurs in a county” and “ends on or before the first date of an 80 percent probability of a 28 degree Fahrenheit temperature (freeze date) occurring in a county” (Farm Service Agency 2018, 4). So, for most of the US (including Missoula), grazing period only depends on the 80th percentile of the minimum temperature; it is defined as the period of time when the 80th percentile of the minimum temperature is 28ºF or higher.

round_season <-
  function(x){
    if(is.na(x[1]))
      return(x)
    season <- as.Date(x, origin = '2017-12-31')
    if(x[1] >= x[2])
      lubridate::year(season[2]) <- lubridate::year(season[2]) + 1
    
    if(lubridate::mday(season[1]) < 15){
      lubridate::day(season[1]) <- 1
    }else{
      lubridate::day(season[1]) <- 15
    }
    
    if(lubridate::mday(season[2]) < 15){
      lubridate::day(season[2]) <- 14
    }else{
      lubridate::day(season[2]) <- lubridate::days_in_month(season[2])
    }
    
    if(diff(season) > 365)
      return(c(1,365))
    
    return(lubridate::yday(season))
    
  }

missoula_start_end <-
  missoula_quants_smooth %>%
  dplyr::filter(type == "tmin") %$%
  which(`80%` >= 28) %>%
  range() %>%
  magrittr::add(lubridate::as_date("1992-12-31"))

missoula_start_end_rounded <-
  round_season(missoula_start_end) %>%
  magrittr::add(lubridate::as_date("1992-12-31"))


ggplot2::ggplot(
  data = missoula %>%
    tidyr::pivot_longer(tmin:tmax,
                        names_to = "type"),
  mapping = aes(x = yday + lubridate::as_date("1992-12-31"))) +
  ggplot2::geom_rect(data = background,
                     mapping = aes(xmin = start,
                                   xmax = end),
                     ymin = -Inf,
                     ymax = Inf,
                     fill = "gray95") +
  geom_ribbon(data = missoula_quants_smooth %>%
                dplyr::filter(type == "tmin") %>%
                dplyr::mutate(`80up` = ifelse(`80%` >= 28, `80%`, NA)), 
              mapping = aes(ymin = 28,
                            ymax = `80up`),
              fill = "dodgerblue",
              alpha = 0.5
  ) + 
  geom_line(
    data = missoula_quants_smooth %>%
      dplyr::filter(type == "tmin"), 
    mapping = ggplot2::aes(y = `80%`,
                           # color = type,
                           group = type), 
    color = "dodgerblue",
    linewidth = 1) +
  ggplot2::geom_hline(yintercept = 28, color = "black") +
  ggplot2::scale_y_continuous(limits = c(15, 50),
                              breaks = c(20, 28, 30, 40, 50),
                              name = "Temperature (ºF)") +
  ggplot2::scale_color_manual(values = c(tmin = "dodgerblue",
                                         tmax = "red")) +
  scale_x_date(date_breaks = "month",
               date_labels = "%b",
               expand = expansion(0,0)
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        panel.grid = element_blank()
  )

The 80th percentile of the normal daily minimum temperatures in Missoula County, Montana (blue line), compared to the 28ºF threshold for defining the normal grazing period (black line). The shaded area represents the un-rounded normal grazing period for Missoula County.

Defined as above, the normal grazing period calculated for Missoula County begins on April 20 and ends on October 24. NAP-190 also allows for rounding the season to the 1st, 15th, or last day of the month. We feel that the grazing period should only be rounded so as to expand the grazing period; thus, the start day is rounded earlier to the 1st or 15th of the month, and the end day is rounded later to the 14th or last day of the month. Thus, the rounded normal grazing period for Missoula County is April 15 – October 31.

1991–2020 Normal Grazing Periods across the contiguous United States

We can now calculate grazing periods using the methods described above for every county in thie contiguous United States (CONUS). Grazing periods could, in principal, be calculated for areas outside the CONUS as well (e.g., Alaska, Hawaii, and other US territories). However, here we use a dataset — NOAA’s nClimGrid-daily gridded fields of daily weather observations from January 1, 1951–present (Durre et al. 2022) — that is only available within the CONUS. future work will identify appropriate gridded datasets for calculating normal grazing periods elsewhere.

nClimGrid-daily data

The NOAA nClimGrid-daily dataset is NOAA’s operational gridded weather dataset. The data are derived from in situ weather station observations from the Global Historical Climatology Network daily data (GHCNd), and processed to capture weather variability due to topographic and other landscape effects, while addressing the temporal and spatial idiosyncrasies in the GHCNd data. The nClimGrid-daily dataset is produced at a resolution of 1/24 of a degree (0.0417°), or roughly 3.75 km per grid cell (a ~14 km\(^2\) area). Data become available two to three days after the observation date.

For this analysis, we downloaded the nClimGrid-daily minimum (tmin) and maximum (tmax) for all days from January 1, 1991 through December 31, 2020. These data are hosted in an Amazon S3 data bucket, and are free and publicly available in the NetCDF version 4 data format. In the code below, we use the furrr (Vaughan and Dancho 2022) and future (Bengtsson 2021) packages for R to perform the downloads in parallel, which greatly speeds up the processing time.

Calculating normal grazing periods for each county

We then calculated daily summaries for each county in the CONUS, recording the daily minimum and maximum temperature occurrences in each county. We use generalized (1:500k) county boundaries from the US Census TIGER/Line Shapefiles database, accessed through the tigris package for R (Walker 2023). We use the exactextractr package for R (Daniel Baston 2022), which enables fast extraction and zonal summaries for gridded datasets. Daily county summaries are stored on disk in the Apache Arrow IPC file format using the arrow package for R (Richardson et al. 2023).

counties <-
  sf::read_sf("data-raw/counties.fgb") %>%
  dplyr::mutate(state = factor(state,
                               ordered = TRUE),
                county = factor(county,
                                ordered = TRUE))

if(!file.exists("data-raw/counties-tmin.arrow")){
  
  ncores <- max(1, future::availableCores() - 2)
  future::plan(multisession,
               workers = ncores)
  
  list.files("data-raw/nclimgrid/",
             pattern = "tmin",
             full.names = TRUE) %>%
    {
      split(., f = ggplot2::cut_number(1:length(.), n = ncores))
    } %>%
    furrr::future_map_dfr(
      .f =
        function(x){
          terra::rast(x) %>%
            magrittr::set_names(., terra::time(.)) %>%
            exactextractr::exact_extract(
              counties,
              fun = "min",
              append_cols = TRUE,
              progress = FALSE,
              colname_fun = function(values, ...){values}) %>%
            tidyr::pivot_longer(-state:-fips,
                                names_to = "date",
                                values_to = "tmin")
        },
      .options = furrr::furrr_options(seed = NULL)) %>%
    dplyr::transmute(state,
                     county,
                     fips,
                     date = lubridate::as_date(date),
                     yday = lubridate::as_date(date) %>%
                       set_year(1993) %>%
                       lubridate::yday(),
                     tmin
    ) %>%
    na.omit() %>%
    dplyr::mutate(tmin =
                    tmin %>%
                    units::set_units("degC") %>%
                    units::set_units("degF") %>%
                    units::drop_units()) %>%
    dplyr::arrange(state, county, date) %>%
    arrow::write_feather("data-raw/counties-tmin.arrow")
  
}

counties_tmin <-
  arrow::read_feather("data-raw/counties-tmin.arrow")

if(!file.exists("data-raw/counties-tmax.arrow")){
  
  ncores <- max(1, future::availableCores() - 2)
  future::plan(multisession,
               workers = ncores)
  
  list.files("data-raw/nclimgrid/",
             pattern = "tmax",
             full.names = TRUE) %>%
    {
      split(., f = ggplot2::cut_number(1:length(.), n = ncores))
    } %>%
    furrr::future_map_dfr(
      .f =
        function(x){
          terra::rast(x) %>%
            magrittr::set_names(., terra::time(.)) %>%
            exactextractr::exact_extract(
              counties,
              fun = "max",
              append_cols = TRUE,
              progress = FALSE,
              colname_fun = function(values, ...){values}) %>%
            tidyr::pivot_longer(-state:-fips,
                                names_to = "date",
                                values_to = "tmax")
        },
      .options = furrr::furrr_options(seed = NULL)) %>%
    dplyr::transmute(state,
                     county,
                     fips,
                     date = lubridate::as_date(date),
                     yday = lubridate::as_date(date) %>%
                       set_year(1993) %>%
                       lubridate::yday(),
                     tmax
    ) %>%
    na.omit() %>%
    dplyr::mutate(tmax =
                    tmax %>%
                    units::set_units("degC") %>%
                    units::set_units("degF") %>%
                    units::drop_units()) %>%
    dplyr::arrange(state, county, date) %>%
    arrow::write_feather("data-raw/counties-tmax.arrow")
  
}

counties_tmax <-
  arrow::read_feather("data-raw/counties-tmax.arrow")

We calculate sample L-moments as above and estimate parameters for the Gumbel distribution using the lmom package for R (J. R. M. Hosking 2022), which includes routines for very fast parameter estimation for several families of distributions. We smooth the parameter distributions annually using a generalized additive model with integrated smoothness estimation using a 12-dimension (e.g., monthly) cyclic cubic regression spline, using functions provided by the mgcv package for R (Wood 2011).

cyclic_smooth <-
  function(y){
    if(is.na(y[1]))
      return(rep(NA,length(y)))
    
    as.numeric(
      predict(
        mgcv::gam(
          formula = y ~ s(x, k = 12, bs = "cc"),
          data = data.frame(x = 1:365,
                            y = y)
        ),
        type = "response")
    )
  }


smooth_dist <-
  function(x){
    
    x %>%
      dplyr::bind_rows() %>%
      dplyr::mutate(dplyr::across(dplyr::everything(), cyclic_smooth)) %>%
      dplyr::rowwise() %>%
      dplyr::group_split() %>%
      purrr::map(unlist)
  }

counties_data <-
  dplyr::full_join(counties_tmin, counties_tmax,
                   by = dplyr::join_by(state, county, fips, date, yday))

counties_dists <-
  counties_data %>%
  dplyr::group_by(state, county, fips, yday) %>%
  dplyr::summarise(
    tmin = list(
      lmom::pelgum(
        lmom::.samlmu(-1 * tmin)
      )
    ),
    tmax = list(
      lmom::pelgum(
        lmom::.samlmu(tmax)
      )
    ),
    .groups = "drop"
  )

counties_dists_smooth <-
  counties_dists %>%
  dplyr::group_by(state, county, fips) %>%
  dplyr::mutate(dplyr::across(tmin:tmax, smooth_dist)) %>%
  dplyr::ungroup()

We then calculate full-, cool-, and warm-season normal grazing periods. Full-season grazing periods are defined as described above for Missoula County. In areas of the country that never have above an 80% probability of a 28ºF minimum temperature, the normal grazing period is the full year. Cool-season grazing periods are defined using the modal minimum and maximum daily temperatures, and are rounded earlier and later as described above; warm-season grazing periods are then defined as the period after cool-season grazing within the full-season normal grazing period for the county (e.g., the warm-season starts the day after the cool-season ends, and continues through the end of the full season or until the next cool season starts, whichever comes first). Note that seasons are specified climatologically from the nadir of the annual daily minimum temperature cycle; in rare circumstances, grazing periods begin prior to January 1. Warm- and cool-season grazing periods are only calculated for states specified in NAP-190 (Farm Service Agency 2018), and only for counties within those states where valid grazing periods can be calculated.

calc_full_season <-
  function(x, y){
    if(is.na(x[1]))
      return(c(NA,NA))
    
    x_min <- which.min(x)
    
    if(x[x_min] > y)
      return(c(1, length(x)))
    
    if(x_min != 1)
      x <- x[c(x_min:length(x),1:(x_min - 1))]
    
    thresh <- which(x >= y)
    
    out <- range(thresh) - (length(x) - (x_min - 1))
    
    out[out < 1] <- out[out < 1] + length(x)
    
    out
    
  }

calc_cool_season_start <-
  function(x, y){
    if(is.na(x[1]))
      return(NA)
    
    x_min <- which.min(x)
    
    if(x[x_min] > y)
      return(NA)
    
    if(x_min != 1)
      x <- x[c(x_min:length(x),1:(x_min - 1))]
    
    thresh <- which(x >= y)[1]
    
    out <- thresh - (length(x) - (x_min - 1))
    
    out[out < 1] <- out[out < 1] + length(x)
    
    out
    
  }

calc_cool_season_end <-
  function(x, y){
    if(is.na(x[1]))
      return(NA)
    
    x_min <- which.min(x)
    
    if(x[which.max(x)] < y)
      return(NA)
    
    if(x_min != 1)
      x <- x[c(x_min:length(x),1:(x_min - 1))]
    
    thresh <- which(x >= y)[1]
    
    out <- thresh - (length(x) - (x_min - 1))
    
    out[out < 1] <- out[out < 1] + length(x)
    
    out
    
  }

calc_cool_season <- 
  function(min, max){
    return(
      c(
        calc_cool_season_start(min, 50),
        calc_cool_season_end(max, 90)
      )
    )
  }

round_season <-
  function(x){
    if(any(is.na(x)))
      return(c(NA, NA))
    season <-
      x + lubridate::as_date("1992-12-31")
    if(x[1] >= x[2])
      lubridate::year(season[2]) <- lubridate::year(season[2]) + 1
    
    if(lubridate::mday(season[1]) < 15){
      lubridate::day(season[1]) <- 1
    }else{
      lubridate::day(season[1]) <- 15
    }
    
    if(lubridate::mday(season[2]) < 15){
      lubridate::day(season[2]) <- 14
    }else{
      lubridate::day(season[2]) <- lubridate::days_in_month(season[2])
    }
    
    if(diff(season) > 365)
      return(c(1,365) %>%
               magrittr::add(lubridate::as_date("1992-12-31")))
    
    return(season)
    
  }

counties_full_season <-
  counties_dists_smooth %>%
  dplyr::rowwise() %>%
  dplyr::mutate(threshold = 
                  (-1 * lmom::quagum(tmin, f = c(0.2)))) %>%
  dplyr::group_by(state, county, fips) %>%
  dplyr::summarise(full_season = 
                     list(
                       threshold %>%
                         calc_full_season(28) %>%
                         round_season() %>%
                         # magrittr::add(lubridate::as_date("1992-12-31")) %>%
                         magrittr::set_names(c("start","end")) %>%
                         as.list() %>%
                         tibble::as_tibble()
                     ),
                   .groups = "drop"
  ) %>%
  tidyr::unnest(full_season)

counties_cool_season <-
  counties_dists_smooth %>%
  dplyr::rowwise() %>%
  dplyr::mutate(tmin = 
                  -1 * tmin["xi"],
                tmax = tmax["xi"]) %>%
  dplyr::group_by(state, county, fips) %>%
  dplyr::summarise(cool_season = 
                     list(
                       calc_cool_season(tmin, tmax) %>%
                         round_season() %>%
                         # magrittr::add(lubridate::as_date("1992-12-31")) %>%
                         magrittr::set_names(c("start","end")) %>%
                         as.list() %>%
                         tibble::as_tibble()
                     ),
                   .groups = "drop"
  ) %>%
  tidyr::unnest(cool_season)


counties_warm_season <-
  dplyr::full_join(counties_full_season, 
                   counties_cool_season, 
                   by = c("state", "county", "fips")) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(start = max(start.x, end.y + 1),
                end = min(end.x, start.y - 1),
                end = lubridate::as_date(ifelse(start >= end, end + years(1), end))) %>%
  dplyr::select(state, county, fips, start, end)

grazing_periods <-
  list(`Full Season` = counties_full_season,
       `Cool Season` = counties_cool_season,
       `Warm Season` = counties_warm_season) %>%
  dplyr::bind_rows(.id = "Pasture Grazing Type") %>%
  dplyr::transmute(`Pasture Grazing Type` = factor(`Pasture Grazing Type`, 
                                                   levels = c("Full Season", 
                                                              "Cool Season",
                                                              "Warm Season"),
                                                   ordered = TRUE),
                   State = state,
                   County = county,
                   FIPS = fips,
                   `Normal Grazing Period Start Date` = start,
                   `Normal Grazing Period End Date` = end,
                   `Number of Days in Grazing Period` = as.integer(end - start + 1),
                   dplyr::across(`Normal Grazing Period Start Date`:`Normal Grazing Period End Date`,
                                 ~format(.x, "%B %e") %>% stringr::str_replace("  "," "))
  ) %>%
  dplyr::arrange(State, County, `Pasture Grazing Type`) %>%
  arrow::write_feather("nclimgrid-normal-grazing-period.arrow")

Results

The output of this analysis is a table of CONUS counties and their respective full-, cool-, and warm-season normal grazing periods, as derived for the 1991–2020 climatology period.

arrow::read_feather("nclimgrid-normal-grazing-period.arrow") %>%
  na.omit() %>%
  datatable(filter = 'top', 
            extensions = 'Buttons', 
            options = list(
              dom = 'Bfrtip',
              buttons = list(
                list(extend = 'csv', 
                     title = "nclimgrid-normal-grazing-period"),
                list(extend = 'excel', 
                     title = "nclimgrid-normal-grazing-period")
              ),
              initComplete = JS(
                "function(settings, json) {",
                "$(this.api().table().header()).css({'font-size': '75%'});",
                
                "}")
            ),
            rownames = FALSE)  %>%
  DT::formatStyle(columns = 1:7, fontSize = '75%')

These data may be mapped to visualize the start, end, and duration of each grazing period.

# Install mapview if necessary
if (!require("mapview")) {
  install.packages("mapview")
}
library(mapview)
mapviewOptions(
  basemaps = c(),
  homebutton = FALSE,
  query.position = "topright",
  query.digits = 2,
  query.prefix = "",
  legend.pos = "bottomright"
)

(arrow::read_feather("nclimgrid-normal-grazing-period.arrow") %>%
  dplyr::left_join(counties,
                    by = c("State" = "state", 
                           "County" = "county",
                           "FIPS" = "fips")) %>%
  sf::st_as_sf() %>%
  dplyr::mutate(
    # dplyr::across(`Normal Grazing Period Start Date`:`Normal Grazing Period End Date`,
    #               ~paste0(.x, ", 1993") %>% 
    #                 lubridate::as_date(format = "%B %e, %Y")),
    # `Normal Grazing Period End Date` = 
    #   lubridate::as_date(
    #     ifelse(`Normal Grazing Period End Date` <= `Normal Grazing Period Start Date`,
    #            `Normal Grazing Period End Date` + years(1),
    #            `Normal Grazing Period End Date`)
    #   ),
    label = paste0(County, " County, ", State, ": ", `Number of Days in Grazing Period`, " days")
  ) %>%
  dplyr::filter(`Pasture Grazing Type` == "Full Season") %>%
  mapview::mapview(
    .,
    zcol = "Number of Days in Grazing Period",
    layer.name = "Grazing Period<br>(days)",
    label = "label",
    popup = 
      leafpop::popupTable(
        ., 
        zcol = c("County", "State", 
                 "Normal Grazing Period Start Date", 
                 "Normal Grazing Period End Date",
                 "Number of Days in Grazing Period"),
        row.numbers = FALSE, 
        feature.id = FALSE
      )
  ))@map %>%
      leaflet::removeLayersControl() %>%
      leaflet::addTiles(
        urlTemplate = "https://basemap.nationalmap.gov/ArcGIS/rest/services/USGSShadedReliefOnly/MapServer/tile/{z}/{y}/{x}"
      ) %>%
      leaflet::addProviderTiles("Stamen.TonerLines") %>%
      leaflet::addProviderTiles("Stamen.TonerLabels")

References

Bache, Stefan Milton, and Hadley Wickham. 2022. magrittr: A Forward-Pipe Operator for R. https://CRAN.R-project.org/package=magrittr.

Bengtsson, Henrik. 2021. “A Unifying Framework for Parallel and Distributed Processing in R Using Futures.” The R Journal 13 (2): 208–27. https://doi.org/10.32614/RJ-2021-048.

Daniel Baston. 2022. exactextractr: Fast Extraction from Raster Datasets Using Polygons. https://CRAN.R-project.org/package=exactextractr.

Durre, Imke, Anthony Arguez, Carl J Schreck III, Michael F Squires, and Russell S Vose. 2022. “Daily High-Resolution Temperature and Precipitation Fields for the Contiguous United States from 1951 to Present.” Journal of Atmospheric and Oceanic Technology 39 (12): 1837–55.

Farm Service Agency. 2018. Review of Normal Grazing Periods and Carrying Capacities for NAP, LFP, and ELAP. Washington, DC: US Department of Agriculture. https://www.fsa.usda.gov/Internet/FSA_Notice/nap_190.pdf.

Hosking, J. R. M. 2022. L-Moments. https://CRAN.R-project.org/package=lmom.

Hosking, JRM. 1990. “L-Moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics.” Journal of the Royal Statistical Society: Series B (Methodological) 52 (1): 105–24.

Hosking, JRM, and JR Wallis. 1997. Regional Frequency Analysis: An Approach Based on l-Moments. Cambridge University Press, New York.

R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rennie, J, I Durre, and M Palecki. 2022. 2020 U.S. Climate Normals: Daily Gridded Normals, Version 1.0. Washington, DC: National Oceanic; Atmospheric Administration. https://www.ncei.noaa.gov/sites/default/files/2022-09/Documentation_Daily_Gridded_Normals%20V1.0.pdf.

Richardson, Neal, Ian Cook, Nic Crane, Dewey Dunnington, Romain François, Jonathan Keane, Dragoș Moldovan-Grünfeld, Jeroen Ooms, and Apache Arrow. 2023. arrow: Integration to Apache Arrow. https://CRAN.R-project.org/package=arrow.

Vaughan, Davis, and Matt Dancho. 2022. furrr: Apply Mapping Functions in Parallel Using Futures. https://CRAN.R-project.org/package=furrr.

Walker, Kyle. 2023. tigris: Load Census TIGER/Line Shapefiles. https://CRAN.R-project.org/package=tigris.

Wickham, Hadley, Romain François, Lionel Henry, Kirill Müller, and Davis Vaughan. 2023. dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, Davis Vaughan, and Maximilian Girlich. 2023. tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

Wood, S. N. 2011. “Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models.” Journal of the Royal Statistical Society (B) 73 (1): 3–36.

Reuse

Copyright

R. Kyle Bocinsky, Montana Climate Office