This post shows some examples of how to find and visualize USGS stream gauge data using R. The {dataRetrieval} package (DeCicco et al. 2024) is used to find and download USGS stream gauge data.
This tutorial on {dataRetrieval} was very helpful to get started using the package.
I will look at data from stream gauges in North Carolina during hurricane Helene in late September, 2024. Helene brought extreme amounts of rainfall to western North Carolina, causing extremely destructive flooding and unfortunately killing and injuring many people.
This is not meant to be an exhaustive analysis of the storm or flooding; rather it is meant to demonstrate some ways you can visualize and analyze stream gauge data in R. I hope it inspires you to create you own analyses.
library(dataRetrieval)
library(leaflet)
library(tidyverse)
library(DT)
library(plotly)
library(tidycensus)
options(tigris_use_cache = TRUE)
library(sf)I used the {tidycensus} (Walker and Herman 2024) package to get county boundaries for the state. I wanted to plot the county boundaries on the map, and have the option to filter by county if I want.
{tidycensus} uses {tigris} to get the geometries if requested (geometry = TRUE). If you just want the geometries, you can use the {tigris} package (Walker 2024) directly. By using {tidycensus} you can get other census variables such as population in the same dataframe.
nc_counties <- get_decennial(geography = "county",
state = "NC",
variables = "P001001", # total population,
year = 2010,
geometry = TRUE)
DT::datatable(nc_counties, options = list(pageLength = 5), rownames = FALSE)First we need to find what stream gauge sites are available. We can use the whatNWISsites() function from {dataRetrieval} to find all sites in North Carolina with a discharge parameter.
nc_sites <- whatNWISsites(stateCd = "NC",
parameterCd = "00060" # discharge parameter code
)
DT::datatable(nc_sites, options = list(pageLength = 5), rownames = FALSE)We can use {leaflet} (Cheng et al. 2024) to make an interactive map of stream gauge locations. Clicking on a marker will display the name and id of a site.
leaflet(data = nc_sites) |>
addProviderTiles(provider = providers$CartoDB.Voyager) |>
addPolygons(data = nc_counties |> sf::st_transform(4326), color = "black", weight = 1,
fillOpacity = 0.1, label = ~NAME) |>
addMarkers(lat = ~dec_lat_va, lng = ~dec_long_va,
popup = paste0(nc_sites$station_nm, "<br>", nc_sites$site_no ) ,
clusterOptions = markerClusterOptions())Once we identify a site to look at, we can get the timeseries data for that site using the readNWISuv() function from {dataRetrieval}.
df_inst <- readNWISuv(siteNumbers = "03451500",
parameterCd = "00060",
startDate = "2024-09-20",
endDate = "2024-10-05",tz = "US/Eastern" ) |>
renameNWISColumns()
df_inst |>
head(5) |>
DT::datatable(rownames = FALSE)Now we can plot the station data. Figure Figure 2 shows the timeseries of discharge at the selected stream gauge. Note that there is data missing during the peak flooding; I’m not sure if the measurement capabilities of the instrument were exceeded.
I wrote a helper function called add_labels_nwis() to rename the plot labels using the variableInfo attributes from the data frame. This replaces the variable code with the actual name of the variable.
##
# function to add better labels to plot using attributes of NWIS data frame
add_labels_nwis <- function(df, g){
parameterInfo <- attr(df, "variableInfo")
siteInfo <- attr(df, "siteInfo")
g <- g +
ylab(parameterInfo$variableDescription) +
ggtitle(siteInfo$station_nm)
return(g)
}
##
g <- df_inst |>
ggplot(aes(dateTime, Flow_Inst)) +
geom_point(color = "lightblue") +
geom_vline(xintercept = as.numeric(as.POSIXct("2024-09-26"))) +
labs(x = NULL) +
theme_minimal()
g <- add_labels_nwis(df_inst,g)
g
We can use the {plotly} package (Sievert 2020) to make an interactive plot that allows us to zoom in/out and hover over data points. Figure Figure 3 shows an interactive version of the previous figure. You can zoom and hover over data to display the values.
ggplotly(g)You can also make the plot directly with {plotly} instead of converting it from a ggplot figure:
parameterInfo <- attr(df, "variableInfo")
p <- df_inst |>
plot_ly(x = ~dateTime, y = ~Flow_Inst, type = "scatter", name = "Discharge") |>
add_lines(x = lubridate::ymd("2024-09-26"), y = range(df_inst$Flow_Inst, na.rm = TRUE),
line = list(color = "red", dash = "dash"), name = "Landfall") |>
layout(
xaxis = list(title = parameterInfo$variableDescription),
yaxis = list(title = "Streamflow [ft^3/s]")
)
pThis post showed how you can use the {dataRetrieval} and other packages to find, download, and analyze USGS stream gauge data in R. Data for stream gauges in North Carolina during hurricane Helene were used as an example. I hope you find this useful getting started with your own analysis. My thoughts go out to everyone affected by Helene and the resulting flooding.
R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
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/Denver
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] sf_1.0-18 tidycensus_1.6.6 plotly_4.10.4
[4] DT_0.33 lubridate_1.9.3 forcats_1.0.0
[7] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
[10] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[13] ggplot2_3.5.1 tidyverse_2.0.0 leaflet_2.2.2
[16] dataRetrieval_2.7.16
loaded via a namespace (and not attached):
[1] gtable_0.3.5 bslib_0.8.0 xfun_0.48
[4] tigris_2.1 htmlwidgets_1.6.4 tzdb_0.4.0
[7] leaflet.providers_2.0.0 vctrs_0.6.5 tools_4.4.1
[10] crosstalk_1.2.1 generics_0.1.3 curl_5.2.3
[13] proxy_0.4-27 fansi_1.0.6 pkgconfig_2.0.3
[16] KernSmooth_2.23-24 data.table_1.16.2 uuid_1.2-1
[19] lifecycle_1.0.4 farver_2.1.2 compiler_4.4.1
[22] munsell_0.5.1 sass_0.4.9 htmltools_0.5.8.1
[25] class_7.3-22 yaml_2.3.10 lazyeval_0.2.2
[28] jquerylib_0.1.4 crayon_1.5.3 pillar_1.9.0
[31] classInt_0.4-10 cachem_1.1.0 tidyselect_1.2.1
[34] rvest_1.0.4 digest_0.6.37 stringi_1.8.4
[37] labeling_0.4.3 fastmap_1.2.0 grid_4.4.1
[40] colorspace_2.1-1 cli_3.6.3 magrittr_2.0.3
[43] utf8_1.2.4 e1071_1.7-16 withr_3.0.1
[46] scales_1.3.0 rappdirs_0.3.3 timechange_0.3.0
[49] rmarkdown_2.28 httr_1.4.7 hms_1.1.3
[52] evaluate_1.0.1 knitr_1.48 viridisLite_0.4.2
[55] rlang_1.1.4 Rcpp_1.0.13 glue_1.8.0
[58] DBI_1.2.3 xml2_1.3.6 renv_1.0.9
[61] rstudioapi_0.16.0 jsonlite_1.8.9 R6_2.5.1
[64] units_0.8-5