First look: USGS Stage/Discharge Gages¶
Overview: Download Data using the USGS Data Retrieval (R Package )¶
USGS Package Overview¶
GitHub Repo¶
CRAN¶
In [1]:
require(dataRetrieval)
require(curl)
library(repr)
library(maps)
library(ggplot2)
library(leaflet)
#library(skimr)
library(htmlwidgets)
options(repr.matrix.max.cols = 40)
Select states of interest & search for available USGS gages¶
- Pull metadata from stations with flow & stage records
- Plot
In [2]:
# Select states to refine search
states <- c("TX")
#--- Assign parameter of interest
stageCode <- "00065"
flowCode <- "00060"
#bBox not working directly
#stage_sites <- whatNWISsites(bBox= bbox, parameterCd=stageCode)
stage_sites <- data.frame()
flow_sites <- data.frame()
for(state in states){
print(state)
tmp <- whatNWISsites(stateCd=state ,parameterCd=stageCode)
stage_sites <- rbind(tmp, stage_sites)
tmp <- whatNWISsites(stateCd=state ,parameterCd=flowCode)
flow_sites <- rbind(tmp, flow_sites)
}
stage_sites$gtype = stageCode
flow_sites$gtype = flowCode
map('state', region = c('texas'))
#map('state', region = c('mississippi'))
title(main="USGS Gauges found in the NWIS")
points(x=flow_sites$dec_long_va,
y=flow_sites$dec_lat_va,
col='black')
points(x=stage_sites$dec_long_va,
y=stage_sites$dec_lat_va,
col='blue')
head(flow_sites)
Refine gage list using bounding box¶
- Select gages within a refined areas of interest
- Filter out gages with limited data
In [3]:
# Bounding box
bbox <- c(29.15, -96.00, 30.5, -94.5)
df <- rbind(stage_sites, flow_sites)
df <- df[df$dec_lat_va <= bbox[3] & df$dec_lat_va >= bbox[1],]
df <- df[df$dec_long_va >= bbox[2] & df$dec_long_va <= bbox[4],]
s <- df[df['gtype'] == stageCode,]
q <- df[df['gtype'] == flowCode,]
head(q, 5)
Filter by # of records available¶
In [4]:
nrecords = 10000
s_meta <- whatNWISdata(siteNumber=q$site_no, parameterCd=stageCode)
q_meta <- whatNWISdata(siteNumber=q$site_no, parameterCd=flowCode)
mysdata <- s_meta[s_meta['count_nu'] > nrecords,]
myqdata <- q_meta[q_meta['count_nu'] > nrecords,]
head(myqdata)
Let's look at data by location¶
In [5]:
pal <- colorFactor(c("black","blue"), domain = c("00060", "00065"))
m <- leaflet(
data=na.omit(df), height=500, width=1000) %>%
addTiles('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
addCircleMarkers(~dec_long_va,~dec_lat_va, color = ~pal(gtype),
radius=5, stroke=FALSE,fillOpacity = 0.9, opacity = 0.8,
popup=~c(station_nm)) %>%
addLegend("topright", pal = pal, values = ~gtype,
title = "USGS Gages in Study Area",
labFormat = labelFormat(prefix = "")
)
head(df,2)
saveWidget(m, 'HoustonGages.html')
m
In [10]:
regex <- 'Bayou'
idx <- grep(regex,q_meta$station_nm)
dc_q_gages <- q_meta[idx,]#$station_nm
idx <- grep(regex,s_meta$station_nm)
dc_s_gages <- s_meta[idx,]#$station_nm
m2 <- (leaflet(data=dc_s_gages, height=500, width=1000) %>%
addTiles('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
addCircleMarkers(~dec_long_va,~dec_lat_va, color = "red",
radius=8, stroke=FALSE, fillOpacity = 0.9,
opacity = 0.8, popup=~c(site_no)) )
head(dc_q_gages,3)
saveWidget(m2, 'BayouGages.html')
m2
Click HERE if points don't plot
Pull Metadata¶
In [7]:
s_meta <- whatNWISdata(siteNumber=s$site_no, parameterCd=stageCode)
q_meta <- whatNWISdata(siteNumber=q$site_no, parameterCd=flowCode)
options(repr.plot.width=7, repr.plot.height=2)# Set plotsize
qplot(s_meta$begin_date, geom="histogram", binwidth = 50) # Plot stage counts by starting year
qplot(q_meta$begin_date, geom="histogram", binwidth = 75) # Plot flow counts by starting year
Download Flow Data & Plot¶
In [8]:
sites <- q$site_no
print(sites)
siteNo <- sites[3]
gage_data <- readNWISdv(siteNumbers=siteNo, parameterCd=stageCode, startDate = "1920-01-01", endDate = "2018-01-01")
gage_data <- renameNWISColumns(gage_data)
head(gage_data,2)
In [9]:
siteNo <- sites[15]
gage_data <- readNWISdv(siteNumbers=siteNo, parameterCd=flowCode, startDate = "1920-01-01", endDate = "2018-01-01")
gage_data <- renameNWISColumns(gage_data)
head(gage_data,2)
parameterInfo <- attr(gage_data, "variableInfo")
siteInfo <- attr(gage_data, "siteInfo")
options(repr.plot.width=8, repr.plot.height=3)
ts <- ggplot(data = gage_data,
aes(Date, Flow)) +
geom_line()+
xlab(siteInfo$station_nm) +
ylab(parameterInfo$parameter_desc) +
ggtitle(siteNo)
ts
#skim(gage_data) %>% skimr::kable()
