knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE
)

1 Learning objectives

In this lab, we will treat stream gages as windows into integrated watershed behavior.

You will learn how to:

  1. clone the course GitHub repository;
  2. locate monitoring stations in USGS Water Data for the Nation;
  3. query the modern USGS Water Data APIs from R;
  4. download discharge, water temperature, and turbidity data;
  5. calculate annual water yield;
  6. compare turbidity among watersheds; and
  7. evaluate whether differences in data availability affect the conclusions we can draw.

2 Clone the GitHub repository

Our course repository is:

https://github.com/kellyloria/GPHS_782_watershed_sci_seminar.git

2.1 Using RStudio

Go to:

File → New Project → Version Control → Git

Then enter:

https://github.com/kellyloria/GPHS_782_watershed_sci_seminar.git

Choose a location on your computer and click Create Project.

2.2 Using the terminal

git clone https://github.com/kellyloria/GPHS_782_watershed_sci_seminar.git
cd GPHS_782_watershed_sci_seminar

Checkpoint: Look in the Files pane in RStudio. Are you working inside the cloned repository?

3 Explore USGS Water Data

Before writing code, open:

The dataRetrieval package provides R functions that construct requests to the modern USGS Water Data APIs.

For this class, we will use functions beginning with:

read_waterdata_*

4 Install and load packages

Run these installation lines once if needed.

install.packages(c(
  "dataRetrieval",
  "tidyverse",
  "lubridate",
  "scales"
))

Load packages:

library(dataRetrieval)
library(tidyverse)
library(lubridate)
library(scales)

5 Our watersheds

We will compare two Lake Tahoe tributaries with the much larger Truckee River.

sites <- tribble(
  ~site_name,        ~monitoring_location_id, ~drainage_area_mi2,
  "Blackwood Creek", "USGS-10336660",          11.2,
  "Incline Creek",   "USGS-10336700",           6.74,
  "Truckee River",   "USGS-10348000",        1067
)

sites
## # A tibble: 3 × 3
##   site_name       monitoring_location_id drainage_area_mi2
##   <chr>           <chr>                              <dbl>
## 1 Blackwood Creek USGS-10336660                      11.2 
## 2 Incline Creek   USGS-10336700                       6.74
## 3 Truckee River   USGS-10348000                    1067

The sites differ greatly in drainage area, so total discharge volume and watershed-normalized runoff answer different questions.

6 Look up station metadata with the API

site_metadata <- read_waterdata_monitoring_location(
  monitoring_location_id = sites$monitoring_location_id
)

site_metadata %>%
  as_tibble() %>%
  select(
    any_of(c(
      "id",
      "monitoring_location_id",
      "monitoring_location_name",
      "state_name",
      "county_name",
      "hydrologic_unit_code",
      "drainage_area"
    ))
  )
## # A tibble: 3 × 6
##   monitoring_location_id monitoring_location_name      state_name county_name  
##   <chr>                  <chr>                         <chr>      <chr>        
## 1 USGS-10336660          BLACKWOOD C NR TAHOE CITY CA  California Placer County
## 2 USGS-10336700          INCLINE CK NR CRYSTAL BAY, NV Nevada     Washoe County
## 3 USGS-10348000          TRUCKEE RV AT RENO, NV        Nevada     Washoe County
## # ℹ 2 more variables: hydrologic_unit_code <chr>, drainage_area <dbl>

7 USGS parameter codes

parameter_codes <- tribble(
  ~parameter_code, ~variable,
  "00060", "Discharge",
  "00010", "Water temperature",
  "63680", "Turbidity"
)

parameter_codes
## # A tibble: 3 × 2
##   parameter_code variable         
##   <chr>          <chr>            
## 1 00060          Discharge        
## 2 00010          Water temperature
## 3 63680          Turbidity

Ask what the parameter codes mean:

parameter_metadata <- read_waterdata_parameter_codes(
  parameter_code = parameter_codes$parameter_code
)

parameter_metadata %>%
  as_tibble() %>%
  select(any_of(c(
    "parameter_code",
    "parameter_name",
    "parameter_description",
    "unit_of_measure"
  )))
## # A tibble: 3 × 4
##   parameter_code parameter_name       parameter_description      unit_of_measure
##   <chr>          <chr>                <chr>                      <chr>          
## 1 00010          Temperature, water   Temperature, water, degre… deg C          
## 2 00060          Discharge            Discharge, cubic feet per… ft3/s          
## 3 63680          Turbidity, Form Neph Turbidity, water, unfilte… FNU

8 Check data availability before downloading

A common mistake is assuming that every USGS gage measures every variable.

ts_metadata <- read_waterdata_ts_meta(
  monitoring_location_id = sites$monitoring_location_id
) %>%
  as_tibble()

names(ts_metadata)
##  [1] "unit_of_measure"               "parameter_name"               
##  [3] "parameter_code"                "statistic_id"                 
##  [5] "hydrologic_unit_code"          "state_name"                   
##  [7] "last_modified"                 "begin"                        
##  [9] "end"                           "begin_utc"                    
## [11] "end_utc"                       "computation_period_identifier"
## [13] "computation_identifier"        "thresholds"                   
## [15] "sublocation_identifier"        "primary"                      
## [17] "monitoring_location_id"        "web_description"              
## [19] "parameter_description"         "parent_time_series_id"        
## [21] "data_gap_interval"             "geometry"                     
## [23] "time_series_id"

Filter for the variables we care about:

availability <- ts_metadata %>%
  filter(parameter_code %in% parameter_codes$parameter_code) %>%
  left_join(parameter_codes, by = "parameter_code") %>%
  left_join(sites, by = "monitoring_location_id")

availability %>%
  select(
    any_of(c(
      "site_name",
      "monitoring_location_id",
      "variable",
      "parameter_code",
      "computation_period_identifier",
      "begin",
      "end",
      "begin_date",
      "end_date",
      "unit_of_measure"
    ))
  ) %>%
  arrange(site_name, variable)
## # A tibble: 50 × 8
##    site_name       monitoring_location_id variable  parameter_code
##    <chr>           <chr>                  <chr>     <chr>         
##  1 Blackwood Creek USGS-10336660          Discharge 00060         
##  2 Blackwood Creek USGS-10336660          Discharge 00060         
##  3 Blackwood Creek USGS-10336660          Discharge 00060         
##  4 Blackwood Creek USGS-10336660          Turbidity 63680         
##  5 Blackwood Creek USGS-10336660          Turbidity 63680         
##  6 Blackwood Creek USGS-10336660          Turbidity 63680         
##  7 Blackwood Creek USGS-10336660          Turbidity 63680         
##  8 Blackwood Creek USGS-10336660          Turbidity 63680         
##  9 Blackwood Creek USGS-10336660          Turbidity 63680         
## 10 Blackwood Creek USGS-10336660          Turbidity 63680         
## # ℹ 40 more rows
## # ℹ 4 more variables: computation_period_identifier <chr>, begin <dttm>,
## #   end <dttm>, unit_of_measure <chr>

Because metadata fields can evolve, any_of() keeps the display step from failing if a column name changes.

9 Define a water year

We will compare Water Year 2024.

start_date <- as.Date("2023-10-01")
end_date   <- as.Date("2024-09-30")

start_date
## [1] "2023-10-01"
end_date
## [1] "2024-09-30"

Water Year 2024 begins October 1, 2023 and ends September 30, 2024.

10 Download daily mean discharge

For discharge, we do not need every 15-minute observation.

  • Parameter 00060 = discharge
  • Statistic 00003 = mean
q_daily <- read_waterdata_daily(
  monitoring_location_id = sites$monitoring_location_id,
  parameter_code = "00060",
  statistic_id = "00003",
  time = c(start_date, end_date)
) %>%
  as_tibble()

glimpse(q_daily)
## Rows: 1,098
## Columns: 11
## $ monitoring_location_id <chr> "USGS-10336660", "USGS-10336700", "USGS-1034800…
## $ parameter_code         <chr> "00060", "00060", "00060", "00060", "00060", "0…
## $ statistic_id           <chr> "00003", "00003", "00003", "00003", "00003", "0…
## $ time                   <date> 2023-10-01, 2023-10-01, 2023-10-01, 2023-10-02…
## $ value                  <dbl> 3.20, 6.81, 363.00, 2.91, 6.75, 350.00, 2.86, 6…
## $ unit_of_measure        <chr> "ft^3/s", "ft^3/s", "ft^3/s", "ft^3/s", "ft^3/s…
## $ approval_status        <chr> "Approved", "Approved", "Approved", "Approved",…
## $ last_modified          <dttm> 2025-03-10 21:59:45, 2025-03-10 23:06:23, 2025…
## $ qualifier              <chr> "", "", "", "", "ESTIMATED", "", "", "", "", ""…
## $ geometry               <POINT [°]> POINT (-120.1621 39.10739), POINT (-119.9…
## $ time_series_id         <chr> "5dacc3ee13924b11ad30295200f3983f", "e94627f4cc…

Clean the discharge data:

q_daily_clean <- q_daily %>%
  transmute(
    monitoring_location_id,
    date = as.Date(time),
    discharge_cfs = as.numeric(value)
  ) %>%
  left_join(sites, by = "monitoring_location_id")

head(q_daily_clean)
## # A tibble: 6 × 5
##   monitoring_location_id date       discharge_cfs site_name    drainage_area_mi2
##   <chr>                  <date>             <dbl> <chr>                    <dbl>
## 1 USGS-10336660          2023-10-01          3.2  Blackwood C…             11.2 
## 2 USGS-10336700          2023-10-01          6.81 Incline Cre…              6.74
## 3 USGS-10348000          2023-10-01        363    Truckee Riv…           1067   
## 4 USGS-10336660          2023-10-02          2.91 Blackwood C…             11.2 
## 5 USGS-10336700          2023-10-02          6.75 Incline Cre…              6.74
## 6 USGS-10348000          2023-10-02        350    Truckee Riv…           1067

11 Check discharge completeness

q_completeness <- q_daily_clean %>%
  group_by(site_name) %>%
  summarise(
    first_date = min(date, na.rm = TRUE),
    last_date = max(date, na.rm = TRUE),
    days_with_data = n_distinct(date[!is.na(discharge_cfs)]),
    missing_q_values = sum(is.na(discharge_cfs)),
    .groups = "drop"
  )

q_completeness
## # A tibble: 3 × 5
##   site_name       first_date last_date  days_with_data missing_q_values
##   <chr>           <date>     <date>              <int>            <int>
## 1 Blackwood Creek 2023-10-01 2024-09-30            366                0
## 2 Incline Creek   2023-10-01 2024-09-30            366                0
## 3 Truckee River   2023-10-01 2024-09-30            366                0

11.0.1 Question

Would you calculate annual water yield if a site were missing several months of discharge data? Why or why not?

12 Plot the hydrographs

ggplot(q_daily_clean, aes(x = date, y = discharge_cfs)) +
  geom_line() +
  facet_wrap(~ site_name, scales = "free_y", ncol = 1) +
  labs(
    x = NULL,
    y = expression("Daily mean discharge (ft"^3*"/s)"),
    title = "Water Year 2024 hydrographs"
  ) +
  theme_minimal()

Consider controls on hydrograph shape:

  • watershed size,
  • snow accumulation and melt,
  • groundwater contributions,
  • channel storage,
  • reservoirs and regulation,
  • diversions,
  • elevation,
  • precipitation timing.

13 Calculate total annual streamflow volume

Daily discharge is reported in cubic feet per second.

\[ V_{day} = Q \times 86,400 \; \frac{seconds}{day} \]

and

\[ 1 \; acre\!-\!foot = 43,560 \; ft^3 \]

Therefore:

cfs_day_to_acre_feet <- 86400 / 43560
cfs_day_to_acre_feet
## [1] 1.983471

Approximately 1 cfs sustained for one day equals 1.9835 acre-feet.

annual_yield <- q_daily_clean %>%
  mutate(
    daily_volume_af = discharge_cfs * cfs_day_to_acre_feet
  ) %>%
  group_by(site_name, drainage_area_mi2) %>%
  summarise(
    days_with_q = sum(!is.na(discharge_cfs)),
    annual_volume_af = sum(daily_volume_af, na.rm = TRUE),
    .groups = "drop"
  )

annual_yield
## # A tibble: 3 × 4
##   site_name       drainage_area_mi2 days_with_q annual_volume_af
##   <chr>                       <dbl>       <int>            <dbl>
## 1 Blackwood Creek             11.2          366           21529.
## 2 Incline Creek                6.74         366            5070.
## 3 Truckee River             1067            366          402575.

Important: na.rm = TRUE does not fill missing days. Interpret annual volume together with days_with_q.

14 Normalize annual flow by watershed area

One square mile contains 640 acres.

annual_yield <- annual_yield %>%
  mutate(
    watershed_area_acres = drainage_area_mi2 * 640
  )

annual_yield
## # A tibble: 3 × 5
##   site_name  drainage_area_mi2 days_with_q annual_volume_af watershed_area_acres
##   <chr>                  <dbl>       <int>            <dbl>                <dbl>
## 1 Blackwood…             11.2          366           21529.                7168 
## 2 Incline C…              6.74         366            5070.                4314.
## 3 Truckee R…           1067            366          402575.              682880

Plot total volume:

ggplot(annual_yield, aes(x = reorder(site_name, annual_volume_af), y = annual_volume_af)) +
  geom_col() +
  coord_flip() +
  scale_y_continuous(labels = comma) +
  labs(
    x = NULL,
    y = "Annual discharge volume (acre-feet)",
    title = "Total annual streamflow volume"
  ) +
  theme_minimal()

14.0.1 Questions

  1. Which site has the largest total annual discharge?
  2. Which site has the largest area-normalized runoff depth?
  3. Why are these rankings not necessarily the same?
  4. Is runoff depth a direct measurement of precipitation? Why not?

15 Download high-frequency water temperature and turbidity

Use the continuous-values API for sensor observations.

sensor_data <- read_waterdata_continuous(
  monitoring_location_id = sites$monitoring_location_id,
  parameter_code = c("00010", "63680"),
  time = c(start_date, end_date)
) %>%
  as_tibble()

glimpse(sensor_data)
## Rows: 503,773
## Columns: 10
## $ monitoring_location_id <chr> "USGS-10336660", "USGS-10336660", "USGS-1033670…
## $ parameter_code         <chr> "00010", "63680", "00010", "63680", "00010", "6…
## $ statistic_id           <chr> "00011", "00011", "00011", "00011", "00011", "0…
## $ time                   <dttm> 2023-10-01 00:00:00, 2023-10-01 00:00:00, 2023…
## $ value                  <dbl> 10.2, 0.4, 7.2, 0.9, 10.2, 0.5, 7.2, 0.8, 10.1,…
## $ unit_of_measure        <chr> "degC", "_FNU", "degC", "_FNU", "degC", "_FNU",…
## $ approval_status        <chr> "Approved", "Approved", "Approved", "Approved",…
## $ qualifier              <chr> "", "", "", "", "", "", "", "", "", "", "", "",…
## $ last_modified          <dttm> 2025-08-11 20:06:32, 2025-08-11 23:13:45, 2025…
## $ time_series_id         <chr> "3de52df18c0d46cc8acc46936cbed148", "4396be785c…

This request may take longer because continuous records contain many observations.

16 Clean the sensor data

sensor_clean <- sensor_data %>%
  transmute(
    monitoring_location_id,
    datetime = as.POSIXct(time),
    parameter_code,
    value = as.numeric(value),
    unit_of_measure = unit_of_measure
  ) %>%
  left_join(parameter_codes, by = "parameter_code") %>%
  left_join(sites, by = "monitoring_location_id")

sensor_clean %>%
  count(site_name, variable)
## # A tibble: 5 × 3
##   site_name       variable               n
##   <chr>           <chr>              <int>
## 1 Blackwood Creek Turbidity          99750
## 2 Blackwood Creek Water temperature 108078
## 3 Incline Creek   Turbidity          98926
## 4 Incline Creek   Water temperature 109388
## 5 Truckee River   Turbidity          87631

17 Summarize water temperature

Calculate daily mean temperature:

temp_daily <- sensor_clean %>%
  filter(parameter_code == "00010") %>%
  mutate(date = as.Date(datetime)) %>%
  group_by(site_name, date) %>%
  summarise(
    mean_temp_c = mean(value, na.rm = TRUE),
    n_obs = sum(!is.na(value)),
    .groups = "drop"
  )

temp_daily
## # A tibble: 730 × 4
##    site_name       date       mean_temp_c n_obs
##    <chr>           <date>           <dbl> <int>
##  1 Blackwood Creek 2023-10-01        8.97   288
##  2 Blackwood Creek 2023-10-02        9.12   288
##  3 Blackwood Creek 2023-10-03        9.71   288
##  4 Blackwood Creek 2023-10-04       10.8    288
##  5 Blackwood Creek 2023-10-05       10.6    288
##  6 Blackwood Creek 2023-10-06       10.5    288
##  7 Blackwood Creek 2023-10-07       10.6    288
##  8 Blackwood Creek 2023-10-08       10.6    288
##  9 Blackwood Creek 2023-10-09       10.5    288
## 10 Blackwood Creek 2023-10-10        9.48   288
## # ℹ 720 more rows

17.0.1 Question

Is mean daily the most helpful response? What about variance from solar noon? What about extremes related to short events?

Plot available records:

ggplot(temp_daily, aes(x = date, y = mean_temp_c, color = site_name)) +
  geom_line() +
  labs(
    x = NULL,
    y = "Daily mean water temperature (°C)",
    color = "Site",
    title = "Available water-temperature records"
  ) +
  theme_minimal()

The Truckee River at Reno is useful for discharge and turbidity, but it does not currently provide the same continuous temperature record as Blackwood and Incline Creeks.

This is a real-world feature of environmental databases: the best station depends on the question and variable.

Optional nearby Truckee River temperature site:

clark_temp <- read_waterdata_continuous(
  monitoring_location_id = "USGS-10350500",
  parameter_code = "00010",
  time = c(start_date, end_date)
)

18 Summarize turbidity carefully

Turbidity parameter 63680 is reported in FNU.

Rather than simply averaging every 15-minute value, first calculate daily summaries.

turbidity_daily <- sensor_clean %>%
  filter(parameter_code == "63680") %>%
  mutate(date = as.Date(datetime)) %>%
  group_by(site_name, date) %>%
  summarise(
    mean_turbidity_fnu = mean(value, na.rm = TRUE),
    median_turbidity_fnu = median(value, na.rm = TRUE),
    max_turbidity_fnu = max(value, na.rm = TRUE),
    n_obs = sum(!is.na(value)),
    .groups = "drop"
  )

turbidity_daily
## # A tibble: 1,035 × 6
##    site_name       date       mean_turbidity_fnu median_turbidity_fnu
##    <chr>           <date>                  <dbl>                <dbl>
##  1 Blackwood Creek 2023-10-01              0.564                  0.5
##  2 Blackwood Creek 2023-10-02              0.735                  0.7
##  3 Blackwood Creek 2023-10-03              0.675                  0.7
##  4 Blackwood Creek 2023-10-04              0.715                  0.7
##  5 Blackwood Creek 2023-10-05              0.722                  0.7
##  6 Blackwood Creek 2023-10-06              0.714                  0.7
##  7 Blackwood Creek 2023-10-07              0.649                  0.6
##  8 Blackwood Creek 2023-10-08              0.687                  0.6
##  9 Blackwood Creek 2023-10-09              0.657                  0.6
## 10 Blackwood Creek 2023-10-10              0.669                  0.6
## # ℹ 1,025 more rows
## # ℹ 2 more variables: max_turbidity_fnu <dbl>, n_obs <int>

Check completeness:

turbidity_completeness <- turbidity_daily %>%
  group_by(site_name) %>%
  summarise(
    days_with_turbidity = n_distinct(date),
    first_day = min(date),
    last_day = max(date),
    .groups = "drop"
  )

turbidity_completeness
## # A tibble: 3 × 4
##   site_name       days_with_turbidity first_day  last_day  
##   <chr>                         <int> <date>     <date>    
## 1 Blackwood Creek                 364 2023-10-01 2024-09-30
## 2 Incline Creek                   348 2023-10-01 2024-09-30
## 3 Truckee River                   323 2023-10-03 2024-09-30

19 Compare mean turbidity among watersheds

turbidity_summary <- turbidity_daily %>%
  group_by(site_name) %>%
  summarise(
    days_with_data = n(),
    mean_daily_turbidity_fnu = mean(mean_turbidity_fnu, na.rm = TRUE),
    median_daily_turbidity_fnu = median(mean_turbidity_fnu, na.rm = TRUE),
    maximum_daily_turbidity_fnu = max(max_turbidity_fnu, na.rm = TRUE),
    .groups = "drop"
  )

turbidity_summary
## # A tibble: 3 × 5
##   site_name       days_with_data mean_daily_turbidity_fnu median_daily_turbidi…¹
##   <chr>                    <int>                    <dbl>                  <dbl>
## 1 Blackwood Creek            364                     1.37                  0.822
## 2 Incline Creek              348                     3.85                  3.42 
## 3 Truckee River              323                     5.99                  5.59 
## # ℹ abbreviated name: ¹​median_daily_turbidity_fnu
## # ℹ 1 more variable: maximum_daily_turbidity_fnu <dbl>
ggplot(
  turbidity_summary,
  aes(
    x = reorder(site_name, mean_daily_turbidity_fnu),
    y = mean_daily_turbidity_fnu
  )
) +
  geom_col() +
  coord_flip() +
  labs(
    x = NULL,
    y = "Mean daily turbidity (FNU)",
    title = "Mean turbidity during Water Year 2024"
  ) +
  theme_minimal()

20 Turbidity through time

ggplot(turbidity_daily, aes(x = date, y = mean_turbidity_fnu)) +
  geom_line() +
  facet_wrap(~ site_name, scales = "free_y", ncol = 1) +
  labs(
    x = NULL,
    y = "Daily mean turbidity (FNU)",
    title = "Turbidity is episodic"
  ) +
  theme_minimal()

20.0.1 Questions

  1. Are the highest-turbidity days associated with high-flow periods?
  2. Would the annual mean or event maximum better describe sediment transport?
  3. Why might the answer depend on the ecological or management question?

21 Put discharge and turbidity together

q_turbidity <- q_daily_clean %>%
  select(site_name, date, discharge_cfs) %>%
  inner_join(
    turbidity_daily %>%
      select(site_name, date, mean_turbidity_fnu),
    by = c("site_name", "date")
  )

head(q_turbidity)
## # A tibble: 6 × 4
##   site_name       date       discharge_cfs mean_turbidity_fnu
##   <chr>           <date>             <dbl>              <dbl>
## 1 Blackwood Creek 2023-10-01          3.2               0.564
## 2 Incline Creek   2023-10-01          6.81              1.30 
## 3 Blackwood Creek 2023-10-02          2.91              0.735
## 4 Incline Creek   2023-10-02          6.75              0.929
## 5 Blackwood Creek 2023-10-03          2.86              0.675
## 6 Incline Creek   2023-10-03          6.66              0.579

Plot turbidity against discharge:

ggplot(
  q_turbidity,
  aes(x = discharge_cfs, y = mean_turbidity_fnu)
) +
  geom_point(alpha = 0.5) +
  facet_wrap(~ site_name, scales = "free") +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    x = expression("Daily mean discharge (ft"^3*"/s; log scale)"),
    y = "Daily mean turbidity (FNU; log scale)",
    title = "Does turbidity increase with streamflow?"
  ) +
  theme_minimal()

Note: Log scales cannot display zero values. Inspect zeros before deciding how to handle them.

22 Summary table

final_summary <- annual_yield %>%
  left_join(
    turbidity_summary,
    by = "site_name"
  ) %>%
  arrange(desc(drainage_area_mi2))

final_summary
## # A tibble: 3 × 9
##   site_name  drainage_area_mi2 days_with_q annual_volume_af watershed_area_acres
##   <chr>                  <dbl>       <int>            <dbl>                <dbl>
## 1 Truckee R…           1067            366          402575.              682880 
## 2 Blackwood…             11.2          366           21529.                7168 
## 3 Incline C…              6.74         366            5070.                4314.
## # ℹ 4 more variables: days_with_data <int>, mean_daily_turbidity_fnu <dbl>,
## #   median_daily_turbidity_fnu <dbl>, maximum_daily_turbidity_fnu <dbl>

23 Inspect the API request itself

The dataRetrieval functions construct API requests for us.

A daily mean discharge request follows this logic:

https://api.waterdata.usgs.gov/ogcapi/v0/collections/daily/items
    ?monitoring_location_id=USGS-10336660
    &parameter_code=00060
    &statistic_id=00003
    &time=2023-10-01/2024-09-30

The important parts are:

  • where? monitoring location
  • what? parameter code
  • which statistic? mean
  • when? date range

24 Reproducibility

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.4.0         lubridate_1.9.5      forcats_1.0.1       
##  [4] stringr_1.6.0        dplyr_1.2.1          purrr_1.2.2         
##  [7] readr_2.2.0          tidyr_1.3.2          tibble_3.3.1        
## [10] ggplot2_4.0.3        tidyverse_2.0.0      dataRetrieval_2.7.25
## 
## loaded via a namespace (and not attached):
##  [1] utf8_1.2.6         sass_0.4.10        generics_0.1.4     class_7.3-24      
##  [5] KernSmooth_2.23-26 stringi_1.8.9      hms_1.1.4          digest_0.6.39     
##  [9] magrittr_2.0.5     timechange_0.4.0   evaluate_1.0.5     grid_4.6.1        
## [13] RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0     e1071_1.7-17      
## [17] DBI_1.3.0          httr2_1.3.0        jquerylib_0.1.4    cli_3.6.6         
## [21] rlang_1.3.0        units_1.0-1        withr_3.0.3        cachem_1.1.0      
## [25] yaml_2.3.12        otel_0.2.0         tools_4.6.1        tzdb_0.5.0        
## [29] curl_7.1.0         vctrs_0.7.3        R6_2.6.1           proxy_0.4-29      
## [33] lifecycle_1.0.5    classInt_0.4-11    pkgconfig_2.0.3    pillar_1.11.1     
## [37] bslib_0.12.0       gtable_0.3.6       data.table_1.18.4  glue_1.8.1        
## [41] Rcpp_1.1.2         sf_1.1-2           xfun_0.60          tidyselect_1.2.1  
## [45] rstudioapi_0.19.0  knitr_1.51         farver_2.1.2       htmltools_0.5.9   
## [49] labeling_0.4.3     rmarkdown_2.31     compiler_4.6.1     S7_0.2.2