knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE
)
In this lab, we will treat stream gages as windows into integrated watershed behavior.
You will learn how to:
Our course repository is:
https://github.com/kellyloria/GPHS_782_watershed_sci_seminar.git
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.
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?
Before writing code, open:
dataRetrieval documentation: https://water.code-pages.usgs.gov/dataRetrieval/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_*
Run these installation lines once if needed.
install.packages(c(
"dataRetrieval",
"tidyverse",
"lubridate",
"scales"
))
Load packages:
library(dataRetrieval)
library(tidyverse)
library(lubridate)
library(scales)
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.
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>
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
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.
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.
For discharge, we do not need every 15-minute observation.
00060 = discharge00003 = meanq_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
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
Would you calculate annual water yield if a site were missing several months of discharge data? Why or why not?
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:
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.
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()
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.
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
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
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)
)
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
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()
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()
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.
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>
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
¶meter_code=00060
&statistic_id=00003
&time=2023-10-01/2024-09-30
The important parts are:
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