What is dataRetrieval?

• R-package to get USGS/EPA water data into R

Where does the data come from?

• US Geological Survey water data
  • National Water Information System (NWIS)
• Water Quality Portal
  • USGS
  • EPA (EPA Storage and Retrieval Data Warehouse = STORET)
  • USDA
  • more being added….

What does dataRetrieval do to the data?

How to discover data?

• Examples will be provided

dataRetrieval is available on the CRAN repository. The CRAN version is the most stable and user-tested:

install.packages("dataRetrieval")

Bug fixes and feature upgrades are vetted through a version of dataRetrieval that is available on a USGS-maintained R repository. To install from that repository:

install.packages("dataRetrieval", 
                 repos=c("http://owi.usgs.gov/R",
                         getOption("repos")))

More information can be found at http://owi.usgs.gov/R/gran.html.

Finally, the absolute cutting-edge version of dataRetrieval can be installed using the devtools package which pulls from GitHub:

library(devtools)
install_github("USGS-R/dataRetrieval")

Once the dataRetrieval package has been installed, it needs to be loaded in order to use any of the functions:

library(dataRetrieval)

There is a vignette that covers the full scope of the dataRetrieval package. It can be accessed with the following command:

vignette("dataRetrieval",package = "dataRetrieval")

Additionally, each function has a help file. These can be accessed by typing a question mark, followed by the function name in the R console:

?readNWISuv

Each function's help file has working examples to demonstrate the usage. The examples may have comments "## Not run". These examples CAN be run, they just are not run by the CRAN maintainors due to the external service calls.

Finally, if there are still questions that the vignette and help files don't answer, please post an issue on the dataRetrieval GitHub page:

The USGS uses various codes for basic retrievals

• Site ID (often 8 or 15-digits)
• Parameter Code (5 digits)
  • Full list:
  • http://help.waterdata.usgs.gov/code/parameter_cd_query?fmt=rdb&inline=true&group_cd=%
• Statistic Code (for daily values)
  • Full list:
  • http://help.waterdata.usgs.gov/code/stat_cd_nm_query?stat_nm_cd=%25&fmt=html

dataRetrieval includes a data set parameterCdFile that allows you to explore the USGS parameter codes. For example, to find all the parameter codes with the word "phosphorus" in the name:

parameterCdFile <- parameterCdFile
names(parameterCdFile)

## [1] "parameter_cd"       "parameter_group_nm" "parameter_nm"      
## [4] "casrn"              "srsname"            "parameter_units"

phosCds <- parameterCdFile[grep("phosphorus",
                                parameterCdFile$parameter_nm,
                                ignore.case=TRUE),]

nrow(phosCds)

## [1] 93

unique(phosCds$parameter_units)

##  [1] "ml"         "nu"         "g"          "ug/l as P"  "mg/l as P" 
##  [6] "mg/kg as P" "ug/l"       "mg/kg"      "%"          "mg/l"      
## [11] "mg/l PO4"   "lb/day"     "mg/kg PO4"  "mg/m2 as P" "lb/d as P" 
## [16] "ug/L as P"  "tons/day"

Knowing a site number (or site numbers), paremeter code (or codes), and start and end date. Let's start by asking for discharge (parameter code = 00060) data for the Yahara River at Windsor, WI (an inlet to Lake Mendota).

siteNo <- "05427718"
pCode <- "00060"
start.date <- "2014-10-01"
end.date <- "2015-09-30"

yahara <- readNWISuv(siteNumbers = siteNo,
                     parameterCd = pCode,
                     startDate = start.date,
                     endDate = end.date)

From the Yahara example, let's look at the data. The column names are:

names(yahara)

## [1] "agency_cd"        "site_no"          "dateTime"        
## [4] "X_00060_00011"    "X_00060_00011_cd" "tz_cd"

The names of the columns are based on the parameter and statistic codes. In many cases, you can clean up the names with a convenience function renameNWISColumns:

yahara <- renameNWISColumns(yahara)
names(yahara)

## [1] "agency_cd"    "site_no"      "dateTime"     "Flow_Inst"   
## [5] "Flow_Inst_cd" "tz_cd"

head(yahara)

##   agency_cd  site_no            dateTime Flow_Inst Flow_Inst_cd tz_cd
## 1      USGS 05427718 2014-10-01 05:00:00        16            A   UTC
## 2      USGS 05427718 2014-10-01 05:15:00        16            A   UTC
## 3      USGS 05427718 2014-10-01 05:30:00        16            A   UTC
## 4      USGS 05427718 2014-10-01 05:45:00        16            A   UTC
## 5      USGS 05427718 2014-10-01 06:00:00        16            A   UTC
## 6      USGS 05427718 2014-10-01 06:15:00        16            A   UTC

summary(yahara)

##   agency_cd           site_no             dateTime                  
##  Length:33449       Length:33449       Min.   :2014-10-01 05:00:00  
##  Class :character   Class :character   1st Qu.:2014-12-29 20:15:00  
##  Mode  :character   Mode  :character   Median :2015-05-09 19:00:00  
##                                        Mean   :2015-04-21 17:06:32  
##                                        3rd Qu.:2015-07-19 14:05:00  
##                                        Max.   :2015-10-01 04:45:00  
##    Flow_Inst      Flow_Inst_cd          tz_cd          
##  Min.   :  5.10   Length:33449       Length:33449      
##  1st Qu.: 15.00   Class :character   Class :character  
##  Median : 20.00   Mode  :character   Mode  :character  
##  Mean   : 26.77                                        
##  3rd Qu.: 25.00                                        
##  Max.   :253.00

The returned data also has several attributes attached to the data frame. To see what the attributes are:

names(attributes(yahara))

## [1] "names"         "row.names"     "class"         "url"          
## [5] "siteInfo"      "variableInfo"  "disclaimer"    "statisticInfo"
## [9] "queryTime"

Each dataRetrieval return should have the attributes url, siteInfo, and variableInfo. Additional attributes are available depending on the data.

To access the attributes:

url <- attr(yahara, "url")
url

## [1] "http://waterservices.usgs.gov/nwis/iv/?site=05427718&format=waterml,1.1&ParameterCd=00060&startDT=2014-10-01&endDT=2015-09-30"

Raw Data

library(ggplot2)
ts <- ggplot(data = yahara,
             aes(dateTime, Flow_Inst)) +
      geom_line()
ts

parameterInfo <- attr(yahara, "variableInfo")
siteInfo <- attr(yahara, "siteInfo")
  
ts <- ts +
      xlab("") +
      ylab(parameterInfo$parameter_desc) +
      ggtitle(siteInfo$station_nm)
ts

Aside from readNWISuv, there are other functions in dataRetrieval that take essentially the same 4 inputs (sites, parameter codes, start date, end date), but deliver data from different NWIS services:

• A single function exists: readNWISdata
• Allow users to specify the service type in the function call
• Available services:
  • Unit (service = "iv")
  • Daily (service = "dv")
  • Groundwater levels (service = "gwlevels")
  • Surface-water measurements (service = "measurement")
  • Water Quality (service = "qw")
  • Site (service = "site")

siteNo <- "05427718"
pCode <- "00060"
start.date <- "2014-10-01"
end.date <- "2015-09-30"

yahara2 <- readNWISdata(siteNumbers = siteNo, parameterCd = pCode,
                     startDate = start.date, endDate = end.date,
                     service = "uv")

This is all great when you know your site numbers.

What do you do when you don't?

NWIS Mapper: http://maps.waterdata.usgs.gov/mapper/index.html

Environmental Data Discovery and Transformation tool: http://cida.usgs.gov/enddat/dataDiscovery.jsp

For dataRetrieval: Become familiar with the possibilities of the web services http://waterservices.usgs.gov/

Let's look for long-term USGS phosphorous data in Wisconsin:

?readNWISdata

pCode <- c("00662","00665")
phWI <- readNWISdata(stateCd="WI", parameterCd=pCode,
                     service="site", seriesCatalogOutput=TRUE)

library(dplyr)
phWI.1 <- filter(phWI, parm_cd %in% pCode) %>%
            filter(count_nu > 300) %>%
            mutate(period = as.Date(end_date) - as.Date(begin_date)) %>%
            filter(period > 15*365)

library(leaflet)
leaflet(data=phWI.1) %>% 
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(~dec_long_va,~dec_lat_va,
                   color = "red", radius=3, stroke=FALSE,
                   fillOpacity = 0.8, opacity = 0.8,
                   popup=~station_nm)

dataTemp <- readNWISdata(stateCd="OH", parameterCd="00010", service="dv")

instFlow <- readNWISdata(sites="05114000", service="iv", 
                   parameterCd="00060", 
                   startDate="2014-05-01T00:00Z",endDate="2014-05-01T12:00Z")
                                                   
instFlowCDT <- readNWISdata(sites="05114000", service="iv", 
                   parameterCd="00060", 
                   startDate="2014-05-01T00:00",endDate="2014-05-01T12:00",
                   tz="America/Chicago")

bBoxEx <- readNWISdata(bBox=c(-83,36.5,-81,38.5), parameterCd="00010")

waterYear <- readNWISdata(bBox=c(-83,36.5,-81,38.5), parameterCd="00010", 
                  service="dv", startDate="2013-10-01", endDate="2014-09-30")

wiGWL <- readNWISdata(stateCd="WI",service="gwlevels")
meas <- readNWISdata(state_cd="WI",service="measurements",format="rdb_expanded")

Water Quality Portal

• Multiple agencies
  • USGS data comes from the NWIS database
  • EPA data comes from the STORET database (this includes many state, tribal, NGO, and academic groups)

• WQP brings data from all these oranizations together and provides it in a single format

• Much more verbose output than NWIS

• To get non-NWIS data, need to use CharacteristicName instead of parameter code.

Much like the convenience functions for NWIS, there's a simple function for retrievals if the site number and parameter code or characteristic name is known.

sitePH <- phWI.1$site_no[1]

nwisQW <- readNWISqw(sitePH, "00665")
ncol(nwisQW)

## [1] 36

wqpQW <- readWQPqw(paste0("USGS-",sitePH),"00665")
ncol(wqpQW)

## [1] 65

Explore these results in RStudio.

Censored data is particularly important for water quality data. Two examples of censored data are: * Left-censored - the data is less than the detection level of the measurement technique * Right-censored - the data is greater than the upper-limit of the measurement technique

NWIS data makes identifing this data easy using the column result_cd.

There's a little more work for WQP data. The following table has renamed the columns for space considerations.

Non-NWIS data might have different ways to indicate censoring.

The true value of the Water Quality Portal is to explore water quality data from different sources.

Become familiar with the possibilities of the web services http://www.waterqualitydata.us/

There's not a function in WQP that returns period of record information like we did above via NWIS data…(that feature may be implemented in the future)

The following function returns sites that have collected phosphorus data in Wisconsin. There's no way to know if that site has collected one sample, or thousands.

phosSites <- whatWQPsites(statecode="WI",characteristicName="Phosphorus")

So, to get that information, you need to actually get that data.

phosData <- readWQPdata(statecode="WI",characteristicName="Phosphorus")
unique(phosData$ResultMeasure.MeasureUnitCode)

##  [1] "mg/l"       "ug/l"       NA           "mg/kg"      "ppb"       
##  [6] "mg/l as P"  "mg/kg as P" "mg/l PO4"   "lb/day"     "mg/kg PO4" 
## [11] "%"

Use a little dplyr to get some information:

siteInfo <- attr(phosData, "siteInfo")

wiSummary <- phosData %>%
  filter(ResultMeasure.MeasureUnitCode %in% c("mg/l","mg/l as P")) %>%
  group_by(MonitoringLocationIdentifier) %>%
  summarise(count=n(),
            start=min(ActivityStartDateTime),
            end=max(ActivityStartDateTime),
            max = max(ResultMeasureValue, na.rm = TRUE)) %>%
  filter(count > 300) %>%
  arrange(-count) %>%
  left_join(siteInfo, by = "MonitoringLocationIdentifier")

Code for map on next slide:

col_types <- c("darkblue","dodgerblue","green4","gold1","orange","brown","red")
leg_vals <- unique(as.numeric(quantile(wiSummary$max, 
                probs=c(0,0.01,0.1,0.25,0.5,0.75,0.9,.99,1), na.rm=TRUE)))
pal = colorBin(col_types, wiSummary$max, bins = leg_vals)
rad <-3*seq(1,4,length.out = 16)
wiSummary$sizes <- rad[as.numeric(cut(wiSummary$count, breaks=16))]
          
leaflet(data=wiSummary) %>% 
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(~dec_lon_va,~dec_lat_va,
                   fillColor = ~pal(max),
                   radius = ~sizes,
                   fillOpacity = 0.8, opacity = 0.8,stroke=FALSE,
                   popup=~station_nm) %>%
  addLegend(position = 'bottomleft',
            pal=pal,
            values=~max,
            opacity = 0.8,
            labFormat = labelFormat(digits = 1), 
            title = 'Max Value')

• The arguments for all dataRetrieval functions concerning dates (startDate, endDate) can be R Date objects, or character strings, as long as the string is in the form "YYYY-MM-DD"

• In R, one vector (or column in a data frame) can only have ONE timezone attribute
  • Sometimes in a single state, some sites will acknowledge daylight savings and some don't
  • dataRetrieval queries could easily span timezones

• Therefore, dataRetrieval converts all date/times to UTC.

• The user can specify a single timezone to override UTC. The allowable tz arguments are listed on the next slide

Use tools from the httr package to get data transfer information, and/or a progress bar on long-running retrievals.

library(httr)
set_config(verbose())
set_config(progress())
daily <- readNWISdv(siteNo, pCode)
-> GET /nwis/dv/?site=05427718&format=waterml%2C1.1&
  ParameterCd=00060&StatCd=00003&startDT=1851-01-01 
  HTTP/1.1
-> Host: waterservices.usgs.gov
-> User-Agent: libcurl/7.43.0 httr/1.1.0 dataRetrieval/2.5.2
-> Accept-Encoding: gzip, deflate
-> Accept: application/json, text/xml, application/xml, */*
-> 
<- HTTP/1.1 200 OK
<- Date: Wed, 09 Mar 2016 18:10:29 GMT
<- Server: GlassFish Server Open Source Edition  4.1
<- Vary: Accept-Encoding
<- Content-Encoding: gzip
<- Content-Type: text/xml;charset=UTF-8
<- X-UA-Compatible: IE=edge,chrome=1
<- ResponseTime: D=258500 t=1457547029642140
<- Connection: close
<- Transfer-Encoding: chunked
<- 
Downloading: 53 kB

dataRetrieval

• These slides:

• http://usgs-r.github.io/dataRetrieval

• Report issues and ask questions:

• https://github.com/USGS-R/dataRetrieval/issues

NWIS

• Water Services:

• http://waterservices.usgs.gov/

• Help:

• http://help.waterdata.usgs.gov/

Water Quality Portal

• http://www.waterqualitydata.us/

Let's work through a problem. Phosphorus ("00665") and Sediment ("80154")

pCode <- c("00665","80154") #Phos and Sed
sitesVA <- readNWISdata(stateCd="VA", parameterCd=pCode,
                        service="site", seriesCatalogOutput=TRUE)
sitesVA.1 <- filter(sitesVA, parm_cd %in% pCode) %>%
             filter(data_type_cd == "qw") %>%
             filter(count_nu > 500) %>%
             filter(duplicated(site_no)) %>%
             arrange(desc(count_nu))

Get samples that have both parameters:

comboData <- readNWISqw(sitesVA.1$site_no[1], 
                        pCode, reshape = TRUE) %>% 
                  filter(!is.na(result_va_00665) & !is.na(result_va_80154))

variableInfo <- attr(comboData, "variableInfo")

plotCombo <- ggplot(data=comboData) +
        geom_point(aes(x=result_va_80154, y=result_va_00665)) +
        scale_y_log10() +
        scale_x_log10() +
        xlab(variableInfo$srsname[variableInfo$parameter_cd == "80154"]) +
        ylab(variableInfo$srsname[variableInfo$parameter_cd == "00665"])
plotCombo

Get samples that have both parameters:

Run a simple linear regression

mod <- lm(log(result_va_00665)~log(result_va_80154),data=comboData)
print(summary(mod))

## 
## Call:
## lm(formula = log(result_va_00665) ~ log(result_va_80154), data = comboData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.44867 -0.30395  0.00387  0.29651  3.06476 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -4.08483    0.04353  -93.84   <2e-16 ***
## log(result_va_80154)  0.47166    0.01285   36.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5241 on 666 degrees of freedom
## Multiple R-squared:  0.6691, Adjusted R-squared:  0.6686 
## F-statistic:  1346 on 1 and 666 DF,  p-value: < 2.2e-16

Plot residuals:

plotResid <- ggplot() +
        geom_point(aes(x=comboData$sample_dt,
                       y=mod$residuals)) +
        geom_hline(yintercept = 0) +
        xlab("") + ylab("Residuals")
plotResid

Plot residuals:

Unit value discharge data is 15 minute data availabe back to 2007.

Daily mean data is available for far longer.

siteNo <- "05427718"
pCode <- "00060"

daily <- readNWISdv(siteNo, pCode, "1990-01-01","2014-12-31")
daily <- renameNWISColumns(daily)
dd <- ggplot(daily, aes(Date,Flow)) +
  xlab("") +
  ylab(attr(daily,"variableInfo")$parameter_nm) +
  geom_line()
dd

Unit value discharge data is 15 minute data availabe back to 2007.

Daily mean data is available for far longer.

The NWIS functions can be used to get multiple sites and multiple parameters in one call.

sites <- c("05427850","05427718") # 2 stations on Yahara
pCodes <- c("00060","00010") #Temperature and discharge

today <- Sys.Date()
last.week <- Sys.Date()-7

multi <- readNWISuv(sites,pCodes, last.week, today)
multi <- renameNWISColumns(multi)

names(multi)

## [1] "agency_cd"     "site_no"       "dateTime"      "Wtemp_Inst"   
## [5] "Flow_Inst"     "Wtemp_Inst_cd" "Flow_Inst_cd"  "tz_cd"

variableInfo <- attr(multi, "variableInfo")
siteInfo <- attr(multi, "siteInfo")

Base R plots:

par(mfcol=c(2,1),oma=c(2,1,1,1), mar=c(1,4,0.1,0.1))

site1 <- multi[multi$site_no == sites[1],]
site2 <- multi[multi$site_no == sites[2],]
plot(site1$dateTime, site1$Flow_Inst, type="l",col="red",
     xlab = "",xaxt="n",
     ylab=variableInfo$param_units[variableInfo$parameterCd == "00060"])

lines(site2$dateTime, site2$Flow_Inst)
plot(site1$dateTime, site1$Wtemp_Inst, type="l",col="red",
     xlab = "", 
     ylab=variableInfo$param_units[variableInfo$parameterCd == "00010"])
lines(site2$dateTime, site2$Wtemp_Inst)
legend("topleft", legend = siteInfo$site_no, 
       cex=0.75, col = c("red","black"),lwd = 1)
box()

Base R plots:

ggplot2: Reshape data from "wide" to "long" first using tidyr package.

library(tidyr)
library(dplyr)
multi.data <- select(multi, site_no, dateTime, Wtemp_Inst, Flow_Inst)
head(multi.data)

##    site_no            dateTime Wtemp_Inst Flow_Inst
## 1 05427718 2016-08-15 05:00:00       19.9       8.1
## 2 05427718 2016-08-15 05:15:00       19.9       7.7
## 3 05427718 2016-08-15 05:30:00       19.8       8.1
## 4 05427718 2016-08-15 05:45:00       19.7       8.1
## 5 05427718 2016-08-15 06:00:00       19.6       7.7
## 6 05427718 2016-08-15 06:15:00       19.6       7.7

longMulti <- gather(multi.data, variable, value, -site_no, -dateTime)
head(longMulti)

##    site_no            dateTime   variable value
## 1 05427718 2016-08-15 05:00:00 Wtemp_Inst  19.9
## 2 05427718 2016-08-15 05:15:00 Wtemp_Inst  19.9
## 3 05427718 2016-08-15 05:30:00 Wtemp_Inst  19.8
## 4 05427718 2016-08-15 05:45:00 Wtemp_Inst  19.7
## 5 05427718 2016-08-15 06:00:00 Wtemp_Inst  19.6
## 6 05427718 2016-08-15 06:15:00 Wtemp_Inst  19.6

ggplot2: Reshape data from "wide" to "long" first using tidyr package.

gp <- ggplot(longMulti, 
             aes(dateTime, value, color=site_no)) +
  geom_line() +
  facet_grid(variable ~ .,scales= "free")
gp

To demonstrate the power of tidyr + dplyr + ggplot2 + dataRetrieval:

sites <- c("05427850","05427718","05428000","05427880") 
pCodes <- c("00060","00010","00055","00065") 

wideMulti <- readNWISuv(sites,pCodes, "2016-03-01","2016-03-08") %>%
         select(-ends_with("_cd"))

siteInfo <- attr(wideMulti, "siteInfo")
paramInfo <- attr(wideMulti, "variableInfo")

longMulti <- gather(wideMulti, variable, value, -site_no, -dateTime) %>%
         mutate(variable = as.factor(variable)) %>%
         mutate(site_no = as.factor(site_no))

levels(longMulti$variable) <- paramInfo$param_units
levels(longMulti$site_no) <- siteInfo$station_nm

gp <- ggplot(longMulti, 
             aes(dateTime, value, color=site_no)) +
  geom_line(size=1.5) + xlab("") + ylab("") +
  facet_grid(variable ~ .,scales= "free") + 
  theme(legend.title=element_blank(),
        legend.position='bottom',legend.direction = "vertical")
gp

To demonstrate the power of tidyr + dplyr + ggplot2 + dataRetrieval: