dataRetrieval.Rmd 28.9 KB
Newer Older
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
1
2
3
4
---
title: "dataRetrieval Tutorial"
author: "<b>USGS</b>: Laura DeCicco"
date: "`r format(Sys.time(), '%d %B, %Y')`"
Laura A DeCicco's avatar
Laura A DeCicco committed
5
6
subtitle: "Advance slides using right and left arrow keys"

Laura A DeCicco's avatar
update  
Laura A DeCicco committed
7
8
9
10
knit: (function(inputFile, encoding) { 
      out_dir <- 'test';
      rmarkdown::render(inputFile,
                        encoding=encoding, 
Laura A DeCicco's avatar
Laura A DeCicco committed
11
                        output_file='dataRetrieval.html') })
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
12
13
output:
  ioslides_presentation:
Laura A DeCicco's avatar
Laura A DeCicco committed
14
15
    css: css/stylesSlides.css
    logo: img/simple-shadow.png
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
16
    smaller: yes
17
    widescreen: yes
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## Introduction

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



## Overview
Laura A DeCicco's avatar
Laura A DeCicco committed
51
<img src="img/2015-12-15_13-52-29.png" alt="Overview" style="width: 800px;"/>
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
52
53
54
55
56
57
58
59
60
61
62
63
64

## Installation

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

```{r echo=TRUE, eval=FALSE}
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:

```{r echo=TRUE, eval=FALSE}
install.packages("dataRetrieval", 
65
                 repos=c("https://owi.usgs.gov/R",
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
66
67
68
                         getOption("repos")))
```

69
More information can be found at [https://owi.usgs.gov/R/gran.html](https://owi.usgs.gov/R/gran.html).
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

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

```{r echo=TRUE, eval=FALSE}
library(devtools)
install_github("USGS-R/dataRetrieval")
```

## dataRetrieval Help

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

```{r echo=TRUE, eval=TRUE}
library(dataRetrieval)
```

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

```{r echo=TRUE, eval=FALSE}
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:

```{r echo=TRUE, eval=FALSE}
?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:

<center>[https://github.com/USGS-R/dataRetrieval/issues](https://github.com/USGS-R/dataRetrieval/issues)</center>

Laura A DeCicco's avatar
Laura A DeCicco committed
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
## US Geological Survey Water National Water Information System (NWIS) Overview 

|Unit Data|
|---------|
|Data reported at the frequency it is collected – e.g. 15 minute|
|Includes current real-time data|


|Daily Data|
|----------|
|Data aggregated to a daily statistic such as mean, min, or max|
|Includes, streamflow, groundwater levels, and water quality sensors|
|These data can can go back many decades|



|Discrete Data|Meta Data|
|-------------|---------|
|Water quality data|Site information|
|Groundwater level|Parameter information|
|Rating curves| |
|Surfacewater measurements| |
|Peak flow| |

128

Laura A DeCicco's avatar
update  
Laura A DeCicco committed
129
130
131
132
133
134
## USGS Basic Web Retrievals

The USGS uses various codes for basic retrievals

* Site ID (often 8 or 15-digits)
* Parameter Code (5 digits)
Laura A DeCicco's avatar
Laura A DeCicco committed
135
136
    + Full list:
    + [http://help.waterdata.usgs.gov/code/parameter_cd_query?fmt=rdb&inline=true&group_cd=%](http://help.waterdata.usgs.gov/code/parameter_cd_query?fmt=rdb&inline=true&group_cd=%)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
137
* Statistic Code (for daily values)
Laura A DeCicco's avatar
Laura A DeCicco committed
138
139
    + Full list:
    + [http://help.waterdata.usgs.gov/code/stat_cd_nm_query?stat_nm_cd=%25&fmt=html](http://help.waterdata.usgs.gov/code/stat_cd_nm_query?stat_nm_cd=%25&fmt=html)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170


Here are some examples of common codes:
<div class="columns-2">

```{r echo=FALSE, eval=TRUE}
library(knitr)
library(pander)

df <- data.frame(pCode = c("00060","00065","00010","00400"),
                 shName = c("Discharge","Gage Height","Temperature","pH"))

names(df) <- c("Parameter Codes", "Short Name")

pander(df)

```

```{r echo=FALSE, eval=TRUE}
df <- data.frame(pCode = c("00001","00002","00003","00008"),
                 shName = c("Maximum","Minimum","Mean","Median"))

names(df) <- c("Statistic Codes", "Short Name")

pander(df)
```

</div>

## USGS Basic Web Retrievals: Parameter Code

Laura A DeCicco's avatar
Laura A DeCicco committed
171
`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:
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
172
173
174

```{r echo=TRUE, eval=TRUE}
parameterCdFile <- parameterCdFile
Laura A DeCicco's avatar
Laura A DeCicco committed
175
names(parameterCdFile)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
176
177
178
179
180
181

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

nrow(phosCds)
Laura A DeCicco's avatar
Laura A DeCicco committed
182
unique(phosCds$parameter_units)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
183
184
185
186
187
188
189

```

## USGS Basic Web Retrievals: Parameter Code (cont.)

```{r echo=FALSE, eval=TRUE}
library(DT)
Laura A DeCicco's avatar
Laura A DeCicco committed
190
datatable(phosCds, rownames = FALSE,options = list(pageLength = 8))
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
```

## USGS Basic Web Retrievals: readNWISuv

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). 

```{r echo=TRUE, eval=TRUE}
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)

```

## USGS Basic Web Retrievals: renameNWISColumns

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

```{r echo=TRUE, eval=TRUE}
names(yahara)
```


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:

```{r echo=TRUE, eval=TRUE}
yahara <- renameNWISColumns(yahara)
names(yahara)
```

## Explore Data

```{r echo=TRUE, eval=TRUE}
head(yahara)
Laura A DeCicco's avatar
Laura A DeCicco committed
230
summary(yahara)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
```

## Explore Data (cont.)

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

```{r echo=TRUE, eval=TRUE}
names(attributes(yahara))

```

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

To access the attributes:

```{r echo=TRUE, eval=TRUE}
url <- attr(yahara, "url")
Laura A DeCicco's avatar
Laura A DeCicco committed
248
url
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
249
250
251
252
253
254
255
256
```

[Raw Data](`r url`)

##  Explore Data (cont.)

```{r echo=TRUE, eval=TRUE, fig.height=3.5}
library(ggplot2)
Laura A DeCicco's avatar
Laura A DeCicco committed
257
ts <- ggplot(data = yahara,
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
258
259
260
261
262
263
264
265
266
267
268
269
270
             aes(dateTime, Flow_Inst)) +
      geom_line()
ts
```
    
##  Use attributes for metadata:

```{r echo=TRUE, eval=TRUE, fig.height=3.5}
parameterInfo <- attr(yahara, "variableInfo")
siteInfo <- attr(yahara, "siteInfo")
  
ts <- ts +
      xlab("") +
271
      ylab(parameterInfo$variableDescription) +
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
272
273
274
275
276
      ggtitle(siteInfo$station_nm)
ts
```


Laura A DeCicco's avatar
Laura A DeCicco committed
277
## USGS Basic Web Retrievals (additional functions)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
278

Laura A DeCicco's avatar
Laura A DeCicco committed
279
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:
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
280

Laura A DeCicco's avatar
Laura A DeCicco committed
281
```{r echo=FALSE, eval=TRUE}
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
282

Laura A DeCicco's avatar
Laura A DeCicco committed
283
284
285
286
287
288
289
290
df <- data.frame(functionName = c("readNWISuv", "readNWISdv",                              "readNWISgwl", "readNWISmeas","readNWISpeak", 
                    "readNWISqw","readNWISrating","readNWISpCode",
                    "readNWISsite"),
                 service = c("Unit", "Daily", "Groundwater Level",
                             "Surface-water", "Peak Flow", 
                             "Water Quality", "Rating Curves",
                             "Parameter Code", "Site"),
                 stringsAsFactors = FALSE)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
291

Laura A DeCicco's avatar
Laura A DeCicco committed
292
names(df) <- c("Function Name", "Data")
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
293

Laura A DeCicco's avatar
Laura A DeCicco committed
294
pander(df)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
295

Laura A DeCicco's avatar
Laura A DeCicco committed
296
```
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
297

Laura A DeCicco's avatar
Laura A DeCicco committed
298
##  USGS Advanced Retrievals: readNWISdata
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
299

Laura A DeCicco's avatar
Laura A DeCicco committed
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
* 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")

```{r eval=FALSE, echo=TRUE}
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")
```                
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
320
321
322
323
324
325
326

##  USGS Advanced Retrievals: Discover Data

This is all great when you know your site numbers.

What do you do when you don't?

Laura A DeCicco's avatar
Laura A DeCicco committed
327
328
329
330
NWIS Mapper:
[http://maps.waterdata.usgs.gov/mapper/index.html](http://maps.waterdata.usgs.gov/mapper/index.html)

Environmental Data Discovery and Transformation tool: [http://cida.usgs.gov/enddat/dataDiscovery.jsp](http://cida.usgs.gov/enddat/dataDiscovery.jsp)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
331

Laura A DeCicco's avatar
Laura A DeCicco committed
332
For dataRetrieval:
Laura A DeCicco's avatar
Laura A DeCicco committed
333
Become familiar with the possibilities of the web services [http://waterservices.usgs.gov/](http://waterservices.usgs.gov/)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
334

Laura A DeCicco's avatar
Laura A DeCicco committed
335
<img src="img/waterservices.png" alt="Overview" title="dataRetrieval Overview" style="width: 450px;"/>
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
336
337
338
339
340
341
342
343
344

##  USGS Advanced Retrievals: readNWISdata

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

```{r echo=TRUE, eval=FALSE}
?readNWISdata

pCode <- c("00662","00665")
Laura A DeCicco's avatar
Laura A DeCicco committed
345
346
phWI <- readNWISdata(stateCd="WI", parameterCd=pCode,
                     service="site", seriesCatalogOutput=TRUE)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
347

Laura A DeCicco's avatar
Laura A DeCicco committed
348
library(dplyr)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
349
350
351
352
353
354
355
356
357
358
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)

```

##  readNWISdata

```{r echo=FALSE, eval=TRUE, message=FALSE}
Laura A DeCicco's avatar
Laura A DeCicco committed
359
library(dplyr)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
360
pCode <- c("00662","00665")
Laura A DeCicco's avatar
Laura A DeCicco committed
361
362
phWI <- readNWISdata(stateCd="WI", parameterCd=pCode,
                     service="site",
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
                     seriesCatalogOutput=TRUE)

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) %>%
            arrange(-count_nu)

datatable(phWI.1[,c("site_no","station_nm","data_type_cd","begin_date","end_date","count_nu")],rownames = FALSE, 
            options = list(pageLength = 10))
```
 
##  Let's look on a map 

```{r echo=TRUE, eval=TRUE, message=FALSE, fig.height=4}
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)
        
```
    
    
##  More readNWISdata examples     
    
```{r eval=FALSE, echo=TRUE}

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 (WQP)

[Water Quality Portal](http://www.waterqualitydata.us/)

* 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.


## Water Quality Portal Basic Retrieval: readWQPqw

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.

```{r echo=TRUE, eval=TRUE}
sitePH <- phWI.1$site_no[1]

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

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

```

Explore these results in RStudio.

## Censored Data: NWIS

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. 

```{r eval=TRUE, echo=FALSE}
Laura A DeCicco's avatar
Laura A DeCicco committed
456
subNWIS <- select(nwisQW, startDateTime, result_va,remark_cd)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
457
458
459
460
461
462
463
464
465
466
datatable(subNWIS,rownames = FALSE, 
            options = list(pageLength = 5))

```

## Censored Data: WQP

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

```{r eval=TRUE, echo=FALSE}
Laura A DeCicco's avatar
Laura A DeCicco committed
467
subWQP <- select(wqpQW, ActivityStartDateTime, DetectionQuantitationLimitMeasure.MeasureValue,ResultMeasureValue,ResultDetectionConditionText) %>%
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
  rename(value=ResultMeasureValue) %>%
  rename(dateTime = ActivityStartDateTime) %>%
  rename(detectionText = ResultDetectionConditionText) %>%
  rename(detectionLimit = DetectionQuantitationLimitMeasure.MeasureValue)
datatable(subWQP,rownames = FALSE, 
            options = list(pageLength = 5))

```

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

## Water Quality Portal Retrieval: readWQPdata

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

Laura A DeCicco's avatar
Laura A DeCicco committed
483
Become familiar with the possibilities of the web services [http://www.waterqualitydata.us/](http://www.waterqualitydata.us/)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
484

485
486
487
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.
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505

```{r echo=TRUE, eval=FALSE}

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

## Water Quality Portal Retrieval: readWQPdata

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

```{r echo=TRUE, eval=FALSE}
phosData <- readWQPdata(statecode="WI",characteristicName="Phosphorus")
unique(phosData$ResultMeasure.MeasureUnitCode)
  
```

```{r echo=FALSE, eval=TRUE}
Laura A DeCicco's avatar
Laura A DeCicco committed
506
phosData <- readRDS("img/phosData.rds")
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
unique(phosData$ResultMeasure.MeasureUnitCode)
```

Use a little `dplyr` to get some information:

```{r eval=TRUE, echo=TRUE}
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) %>%
Laura A DeCicco's avatar
Laura A DeCicco committed
524
  left_join(siteInfo, by = "MonitoringLocationIdentifier")
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
525
526
527

```

Laura A DeCicco's avatar
Laura A DeCicco committed
528

Laura A DeCicco's avatar
update  
Laura A DeCicco committed
529
530
## Water Quality Portal Retrieval: readWQPdata (cont.)

Laura A DeCicco's avatar
Laura A DeCicco committed
531
532
533
Code for map on next slide:

```{r echo=FALSE, eval=TRUE}
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
534
col_types <- c("darkblue","dodgerblue","green4","gold1","orange","brown","red")
Laura A DeCicco's avatar
Laura A DeCicco committed
535
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)))
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
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,
Laura A DeCicco's avatar
Laura A DeCicco committed
551
            labFormat = labelFormat(digits = 1), #transform = function(x) as.integer(x)),
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
552
553
554
555
556
557
            title = 'Max Value')

```

## Water Quality Portal Retrieval: readWQPdata (cont.)

Laura A DeCicco's avatar
Laura A DeCicco committed
558
```{r echo=TRUE, eval=FALSE}
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
559
col_types <- c("darkblue","dodgerblue","green4","gold1","orange","brown","red")
Laura A DeCicco's avatar
Laura A DeCicco committed
560
561
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)))
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
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,
Laura A DeCicco's avatar
Laura A DeCicco committed
577
            labFormat = labelFormat(digits = 1), 
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
            title = 'Max Value')

```


##  Time/Time zone discussion

* 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

## Allowable timezones

```{r echo=FALSE, eval=TRUE}
df <- data.frame(tz = c("America/New_York","America/Chicago","America/Denver",
                        "America/Los_Angeles","America/Anchorage","America/Honolulu",
                        "America/Jamaica","America/Managua","America/Phoenix", "America/Metlakatla"),
                 Name = c("Eastern Time","Central Time","Mountain Time","Pacific Time",
                          "Alaska Time","Hawaii Time","Eastern Standard Time(year-round)","Central Standard Time (year-round)","Mountain Standard Time(year-round)","Pacific Standard Time(year-round"),
                 UTC_offset = c("-05:00","-06:00","-07:00","-08:00","-09:00","-10:00","-05:00","-06:00","-07:00","-08:00"),
                 UTC_DST = c("-04:00","-05:00","-06:00","-07:00","-08:00","-09:00","-05:00","-06:00","-07:00","-08:00"))



pander(df)
```


## Verbose and Progress

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

```{r echo=TRUE, eval=FALSE}
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
```



## More information:

Laura A DeCicco's avatar
Laura A DeCicco committed
647
648
649
650
651
652
653
654
dataRetrieval

* These slides:
  + [http://usgs-r.github.io/dataRetrieval](http://usgs-r.github.io/dataRetrieval)

* Report issues and ask questions:
  + [https://github.com/USGS-R/dataRetrieval/issues](https://github.com/USGS-R/dataRetrieval/issues)

Laura A DeCicco's avatar
update  
Laura A DeCicco committed
655
656
NWIS

Laura A DeCicco's avatar
Laura A DeCicco committed
657
658
* Water Services:
  + [http://waterservices.usgs.gov/](http://waterservices.usgs.gov/)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
659

Laura A DeCicco's avatar
Laura A DeCicco committed
660
661
* Help:
  + [http://help.waterdata.usgs.gov/](http://help.waterdata.usgs.gov/)
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
662
663
664

Water Quality Portal

Laura A DeCicco's avatar
Laura A DeCicco committed
665
666
  + [http://www.waterqualitydata.us/](http://www.waterqualitydata.us/)

Laura A DeCicco's avatar
Laura A DeCicco committed
667
668
669
670
671

## Supplemental Slides: Example Data Exploration

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

Laura A DeCicco's avatar
Laura A DeCicco committed
672
```{r echo=TRUE, eval=TRUE}
Laura A DeCicco's avatar
Laura A DeCicco committed
673
pCode <- c("00665","80154") #Phos and Sed
Laura A DeCicco's avatar
Laura A DeCicco committed
674
675
sitesVA <- readNWISdata(stateCd="VA", parameterCd=pCode,
                        service="site", seriesCatalogOutput=TRUE)
Laura A DeCicco's avatar
Laura A DeCicco committed
676
677
678
sitesVA.1 <- filter(sitesVA, parm_cd %in% pCode) %>%
             filter(data_type_cd == "qw") %>%
             filter(count_nu > 500) %>%
Laura A DeCicco's avatar
Laura A DeCicco committed
679
680
             filter(duplicated(site_no)) %>%
             arrange(desc(count_nu))
Laura A DeCicco's avatar
Laura A DeCicco committed
681
682
683
684
685
686
687
688
689
690
691
692
693

```
 
```{r echo=FALSE}

pander(sitesVA.1[,c("site_no","station_nm","begin_date","end_date")])
```
 
 
## Supplemental Slides: Example Data Exploration (cont.)

Get samples that have both parameters:

Laura A DeCicco's avatar
Laura A DeCicco committed
694
```{r echo=TRUE, eval = FALSE}
Laura A DeCicco's avatar
Laura A DeCicco committed
695
696
comboData <- readNWISqw(sitesVA.1$site_no[1], 
                        pCode, reshape = TRUE) %>% 
Laura A DeCicco's avatar
Laura A DeCicco committed
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
                  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

```

## Supplemental Slides: Example Data Exploration (cont.)

Laura A DeCicco's avatar
Laura A DeCicco committed
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
Get samples that have both parameters:

```{r echo=FALSE, eval = TRUE}
comboData <- readNWISqw(sitesVA.1$site_no[1], pCode, reshape = TRUE) %>% 
                  renameNWISColumns() %>%
                  filter(!is.na(result_va_00665)) %>%
                  filter(!is.na(result_va_80154)) %>%
                  filter(result_va_00665 > 0) %>%
                  filter(result_va_80154 > 0)

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

```


## Supplemental Slides: Example Data Exploration (cont.)

Run a simple linear regression

```{r echo = TRUE, eval = TRUE}
mod <- lm(log(result_va_00665)~log(result_va_80154),data=comboData)
print(summary(mod))
```

## Supplemental Slides: Example Data Exploration (cont.)

Plot residuals:

```{r echo=TRUE, eval=FALSE}

plotResid <- ggplot() +
        geom_point(aes(x=comboData$sample_dt,
                       y=mod$residuals)) +
        geom_hline(yintercept = 0) +
        xlab("") + ylab("Residuals")
plotResid
Laura A DeCicco's avatar
Laura A DeCicco committed
757
758
759

```

Laura A DeCicco's avatar
Laura A DeCicco committed
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
## Supplemental Slides: Example Data Exploration (cont.)

Plot residuals:

```{r echo=FALSE, eval=TRUE}

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

```


Laura A DeCicco's avatar
Laura A DeCicco committed
776
777
778
##  Supplemental Slides: USGS Basic Web Retrievals: readNWISdv


Laura A DeCicco's avatar
Laura A DeCicco committed
779
780
781
782
783
784
785
786
```{r echo=TRUE, eval=FALSE}
siteNo <- "05427718"
pCode <- "00060"

daily <- readNWISdv(siteNo, pCode, "1990-01-01","2014-12-31")
daily <- renameNWISColumns(daily)
dd <- ggplot(daily, aes(Date,Flow)) +
  xlab("") +
787
  ylab(attr(daily,"variableInfo")$variableName) +
Laura A DeCicco's avatar
Laura A DeCicco committed
788
789
790
791
792
793
794
795
796
  geom_line()
dd
```
   
##  Supplemental Slides: USGS Basic Web Retrievals: readNWISdv



```{r echo=FALSE, eval=TRUE}
797
798


Laura A DeCicco's avatar
Laura A DeCicco committed
799
800
801
802
803
siteNo <- "05427718"
pCode <- "00060"

daily <- readNWISdv(siteNo, pCode, "1990-01-01","2014-12-31")
daily <- renameNWISColumns(daily)
804

Laura A DeCicco's avatar
Laura A DeCicco committed
805
806
dd <- ggplot(daily, aes(Date,Flow)) +
  xlab("") +
807
  ylab(attr(daily,"variableInfo")$variableName) +
Laura A DeCicco's avatar
Laura A DeCicco committed
808
809
810
811
  geom_line()
dd
```
    
Laura A DeCicco's avatar
Laura A DeCicco committed
812
 
Laura A DeCicco's avatar
Laura A DeCicco committed
813
814
    
    
Laura A DeCicco's avatar
Laura A DeCicco committed
815
816
## Supplemental Slides: Multiple site and parameters

Laura A DeCicco's avatar
Laura A DeCicco committed
817
818
The NWIS functions can be used to get multiple sites and multiple parameters in one call.

Laura A DeCicco's avatar
Laura A DeCicco committed
819
820
821
822
```{r echo=TRUE, eval=TRUE}
sites <- c("05427850","05427718") # 2 stations on Yahara
pCodes <- c("00060","00010") #Temperature and discharge

Laura A DeCicco's avatar
Laura A DeCicco committed
823
824
825
826
today <- Sys.Date()
last.week <- Sys.Date()-7

multi <- readNWISuv(sites,pCodes, last.week, today)
Laura A DeCicco's avatar
Laura A DeCicco committed
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
multi <- renameNWISColumns(multi)

names(multi)

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

```

## Supplemental Slides: Multiple site and parameters (cont.)

Base R plots:

```{r echo=TRUE, eval=FALSE}
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()


```

## Supplemental Slides: Multiple site and parameters (cont.)

Base R plots:

```{r echo=FALSE, eval=TRUE}
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()

```


## Supplemental Slides: Multiple site and parameters (cont.)

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

```{r echo=TRUE, eval=TRUE, message=FALSE}
library(tidyr)
library(dplyr)
multi.data <- select(multi, site_no, dateTime, Wtemp_Inst, Flow_Inst)
head(multi.data)

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


## Supplemental Slides: Multiple site and parameters (cont.)

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

```{r echo=TRUE, eval=TRUE, fig.height=4}

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

```

## Supplemental Slides: Multiple site and parameters (cont.)

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

```{r echo=TRUE, eval=FALSE}
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

```

## Supplemental Slides: Multiple site and parameters (cont.)

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

```{r echo=FALSE, eval=TRUE, warning=FALSE, fig.height=6}
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"))
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
956

Laura A DeCicco's avatar
Laura A DeCicco committed
957
958
959
960
961
962
963
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))

964
965
966
levels(longMulti$variable) <- paramInfo$unit[which(paste0("X_",paramInfo$variableCode,"_00000") %in% levels(longMulti$variable))]

levels(longMulti$site_no) <- siteInfo$station_nm[which(siteInfo$site_no %in% levels(longMulti$site_no))]
Laura A DeCicco's avatar
Laura A DeCicco committed
967
968
969
970
971
972
973
974
975
976

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

```
Laura A DeCicco's avatar
update  
Laura A DeCicco committed
977