Some R code to pull down a time series of precipitation data for a HUC12 and plot annual maximums.

geoknife_zoo_analysis.R

library(geoknife)
# Use the WBD feature service available from the National Water Census.
stencil <- webgeom(url = "http://cida.usgs.gov/nwc/geoserver/WBD/ows")
query(stencil,'geoms')
# [1] "WBD:huc08"    "WBD:huc12"    "WBD:huc12agg" "WBD:huc12all"
geom(stencil) <- 'WBD:huc12'
query(stencil, 'attributes')
# [1] "ogc_fid"    "tnmid"      "metasource" "sourcedata" "sourceorig" "sourcefeat" "loaddate"   "gnis_id"   
# [9] "areaacres"  "areasqkm"   "states"     "huc12"      "hutype"     "humod"      "tohuc"      "noncontrib"
# [17] "noncontr_1" "shape_leng" "shape_area"
attribute(stencil) <- 'huc12'
# See: http://cida.usgs.gov/nwc/#!waterbudget/achuc/020600031101 or other lookup for where to find a HUC.
values(stencil) <- '020600031101'
# Run a short time when flash flooding happened.
startDate <- "2016-07-30 22:00:00"
endDate <- "2016-07-31 01:00:00"
fabric <- webdata(url = 'http://cida.usgs.gov/thredds/dodsC/stageiv_combined',
                  variables = "Total_precipitation_surface_1_Hour_Accumulation",
                  times = c(as.POSIXct(startDate),
                            as.POSIXct(endDate)))
job <- geoknife(stencil, fabric, wait=TRUE)
data.out <- result(job, with.units=TRUE)
data.out <- result(job, with.units=TRUE)
data.out
# DateTime 020600031101                                        variable statistic   units
# 1 2016-07-30 22:00:00    0.3466465 Total_precipitation_surface_1_Hour_Accumulation      MEAN kg m^-2
# 2 2016-07-30 23:00:00   30.0033680 Total_precipitation_surface_1_Hour_Accumulation      MEAN kg m^-2
# 3 2016-07-31 00:00:00   65.0995800 Total_precipitation_surface_1_Hour_Accumulation      MEAN kg m^-2
# 4 2016-07-31 01:00:00   24.0071700 Total_precipitation_surface_1_Hour_Accumulation      MEAN kg m^-2

# Now run a long time:
query(fabric,'times')
# "2002-01-01 UTC" "2016-08-02 UTC"
times(fabric) <- c('2002-01-01', '2016-08-01')
job_long <- geoknife(stencil, fabric, wait=TRUE)
data.out.long <- result(job_long, with.units=TRUE)

# Use zoo to calculate annual maximum one and three-hour.
library(zoo)
data.out.zoo<-zoo(data.out.long$`020600031101`,data.out.long$DateTime)
maxann.one.hour<-aggregate(data.out.zoo,cut(time(data.out.zoo),"y"),max,na.rm=TRUE)

maxann.one.hour
# 2002-01-01 2003-01-01 2004-01-01 2005-01-01 2006-01-01 2007-01-01 2008-01-01 2009-01-01 2010-01-01 2011-01-01 2012-01-01 2013-01-01 2014-01-01 2015-01-01 2016-01-01 
# 13.69484   20.33760   22.81684   19.75270   27.62504   11.92129   18.48002   25.62986   17.39823   34.54271   18.38567   15.26099   24.77284   39.66375   65.0995
data.out.zoo.threehour<-zoo(rollapply(data.out.long$`020600031101`, 3, sum, fill=NA,na.rm = TRUE),data.out.long$DateTime)
maxann.three.hour<-aggregate(data.out.zoo.threehour,cut(time(data.out.zoo.threehour),"y"),max,na.rm=TRUE)

maxann.three.hour
# 2002-01-01 2003-01-01 2004-01-01 2005-01-01 2006-01-01 2007-01-01 2008-01-01 2009-01-01 2010-01-01 2011-01-01 2012-01-01 2013-01-01 2014-01-01 2015-01-01 2016-01-01 
#   28.37660   39.83037   43.88621   32.82724   48.55959   22.09352   26.18577   36.85747   32.00337   69.22680   30.07654   27.73321   39.49855   51.50629  119.11012 

df_maxxann.one.hour <- data.frame(value = as.vector(maxann.one.hour),
                     time = as.POSIXct(time(maxann.one.hour)))
df_maxxann.three.hour <- data.frame(value = as.vector(maxann.three.hour),
                                  time = as.POSIXct(time(maxann.three.hour)))
plot(df_maxxann.one.hour$time,df_maxxann.one.hour$value,
     ylab = "Precip (mm)", 
     xlab="Year",
     type="p", 
     las=1,tcl=0.3,col="red",
     main="Maximum Annual One Hour Precipitation \n Over Brice Run-Patapsco River HUC12")

plot(df_maxxann.three.hour$time,df_maxxann.three.hour$value,
     ylab = "Precip (mm)",
     xlab="Year",
     type="p",
     las=1,tcl=0.3,col="red",
     main="Maximum Annual Three Hour Precipitation \n Over Brice Run-Patapsco River HUC12")

The rainfall data source for flash flooding in Maryland:

County based analysis for the three hours of intense rainfall using code from here.

Resulting plots showing how much of an outlier this event really was!