In USGS-R/dataRetrieval: Retrieval Functions for USGS and EPA Hydrology and Water Quality Data

library(knitr)
library(dplyr)

knitr::opts_chunk$set(echo = TRUE, 
                      eval = TRUE,
                      warning = FALSE, 
                      message = FALSE,
                      fig.height=5,
                      fig.width=7)

Increasingly there are high frequency sensor data available for water quality data. There's a need to join the sensor and discrete data by the closest time.

This article will discuss how to do that using the data.table package.

Let's look at site "01646500", and a nearby site with a real-time nitrate-plus-nitrite sensor. Our goal is to get the discharge and nitrate-plus-nitrite sensor data aligned with the discrete water quality samples.

library(dataRetrieval)

site_uv <- "01646500"
site_qw <- "USGS-01646580"
pcode_uv <- "99133"
pcode_qw <- "00631"
start_date <- as.Date("2018-01-01")
end_date <- as.Date("2020-01-01")
qw_data <- readWQPqw(site_qw, pcode_qw, 
                     startDate = start_date,
                     endDate = end_date)

uv_data <- readNWISuv(siteNumbers = site_uv, 
                      parameterCd = c(pcode_uv, "00060"),
                      startDate = start_date,
                     endDate = end_date)

The sensor data ("uv" data) at this particular site has 2 columns of data that are important. The first task is to combine those columns. This is rather unique to this particular site and probably won't need to be done generally.

library(dplyr)

uv_trim <- uv_data |> 
  select(uv_date = dateTime, 
         val1 = X_SUNA...Discontinued._99133_00000,
         val2 = X_SUNA_99133_00000,
         flow = X_00060_00000) |> 
  mutate(val_uv = if_else(is.na(val1), val2, val1)) |> 
  select(-val1, -val2)

knitr::kable(head(uv_trim))

Next we'll clean up the discrete water quality "qw" data to make it easy to follow in this tutorial.

qw_trim <- qw_data |> 
  filter(ActivityTypeCode == "Sample-Routine",
         !is.na(ActivityStartDateTime)) |> 
  select(qw_date = ActivityStartDateTime,
         val_qw = ResultMeasureValue,
         det_txt = ResultDetectionConditionText)

knitr::kable(head(qw_trim))

Finally, we'll use the data.table package to do a join to the nearest date. The code to do that is here:

library(data.table)

setDT(qw_trim)[, join_date := qw_date]

setDT(uv_trim)[, join_date := uv_date]

closest_dt <- uv_trim[qw_trim, on = .(join_date), roll = "nearest"]

closest_dt is a data.table object. It similar to a data.frame, but not identical. We can convert it to a data.frame and then use dplyr commands. Note: the whole analysis can be done via data.table, but most examples in dataRetrieval have used dplyr, which is why we bring it back to data.frame. dplyr also has a join_by(closest()) option, but it is more complicated because you can only specify the closeness in either the forward or backwards direction (and we want either direction).

We can calculate "delta_time" - the difference in time between the uv and qw data. We'll probably want to add a threshold that we don't join values if they are too far apart in time. In this example, if the difference is greater than 24 hours, we'll substitute NA.

qw_closest <- data.frame(closest_dt) |> 
  mutate(delta_time = difftime(qw_date, uv_date, 
                               units = "hours"),
         val_uv = if_else(abs(as.numeric(delta_time)) >= 24, NA, val_uv)) |> 
  select(-join_date)

knitr::kable(head(qw_closest))

Here are a few plots to show the applications of these joins:

library(ggplot2)

ggplot(data = qw_closest) +
  geom_point(aes(x = val_uv, y = val_qw)) +
  theme_bw() +
  xlab("Sensor") +
  ylab("Discrete")

ggplot(data = qw_closest) +
  geom_point(aes(x = flow, y = val_qw)) +
  theme_bw() +
  xlab("Discharge") +
  ylab("Concentration") +
  scale_x_log10() +
  scale_y_log10()

ggplot() +
  geom_line(data = uv_trim,
            aes(x = uv_date, val_uv),
            color = "lightgrey") +
  geom_point(data = qw_closest,
             aes(x = qw_date, y = val_qw),
             color= "red") +
  theme_bw() +
  ggtitle("Red dots = discrete samples, grey lines = continuous sensor") +
  xlab("") +
  ylab("Concentration")

USGS-R/dataRetrieval documentation built on Nov. 5, 2024, 9:18 p.m.