Monthly USGS Obs and Precip

require(plyr)
require(stringr)
require(ggplot2)
require(dataRetrieval) #this library is published by USGS
require(geoknife)
require(dplyr)
require(zoo)
require(lubridate)

#####Enter vector USGS site numbers to be evaluated
site_no<-c("05435943")

#####Get data for site numbers
data_avail<-whatNWISdata(siteNumber = c(site_no), parameterCd = c("72019", "00060", "00065")) 
site_desc<-readNWISsite(c(site_no))
#if(any(data_avail$data_type_cd =="gw")){gw_data<-readNWISgwl(site_no)[,c(4, 7, 2)]}
if(any(data_avail$data_type_cd =="gw")){gw_data<-readNWISgwl(site_no, convertType = FALSE)[,c(4, 7, 2)]} 
if(exists("gw_data")){gw_data$lev_dt<-as.Date(gw_data$lev_dt, "%Y-%m-%d")}
if(any(data_avail$data_type_cd =="dv" & data_avail$parm_cd=="72019")){dv_gw_data<-readNWISdv(site_no, parameterCd =  "72019" , statCd = "00001")[,c(3, 4, 2)]}
if(any(data_avail$data_type_cd =="dv" & data_avail$parm_cd=="00060")){dv_cfs_data<-readNWISdv(site_no, parameterCd =  "00060" , statCd = "00003")[,c(3, 4, 2)]}
if(any(data_avail$data_type_cd =="dv" & data_avail$parm_cd=="00065")){dv_lvl_data<-readNWISdv(site_no, parameterCd =  "00065" , statCd = "00003")[,c(3, 4, 2)]} 
if(exists("gw_data")) {colnames(gw_data)<-c("Date", "OBS", "site_no")}
if(exists("dv_gw_data")) {colnames(dv_gw_data)<-c("Date", "OBS", "site_no")}
if(exists("dv_cfs_data")) {colnames(dv_cfs_data)<-c("Date", "OBS", "site_no")}
if(exists("dv_lvl_data")) {colnames(dv_lvl_data)<-c("Date", "OBS", "site_no")}

if(exists("gw_data") & exists("dv_gw_data")) {gwl_data<-rbind(gw_data, dv_gw_data)}
if(exists("gw_data") & exists("dv_gw_data")==FALSE) {gwl_data<-gw_data}
if(exists("gw_data")==FALSE & exists("dv_gw_data")) {gwl_data<-dv_gw_data}
if(exists("gw_data")) {colnames(gwl_data)<-c("Date", "OBS", "site_no")}

if(exists("gwl_data")){gwl_data$OBS<-as.numeric(gwl_data$OBS)}
if(exists("gwl_data")){gwl_data$OBS<-gwl_data$OBS*(-1)}
if(exists("gwl_data")){gwl_data$Unit<-"Ft_BLS"}
if(exists("gwl_data")){gwl_data<-join(gwl_data, site_desc, by= "site_no")[, c(1,2,3,4,6,7,10,11,23)]}
if(exists("gwl_data")){gwl_data$Wat_elev<-gwl_data$alt_va+gwl_data$OBS}
if(exists("dv_cfs_data")){dv_cfs_data$Unit<-"CFS"}
if(exists("dv_cfs_data")){dv_cfs_data<-join(dv_cfs_data, site_desc, by= "site_no")[, c(1,2,3,4,6,7,10,11,23)]}
if(exists("dv_cfs_data")){dv_cfs_data$Wat_elev<-NA}
if(exists("dv_lvl_data")){dv_lvl_data$Unit<-"Elev_Ft"}
if(exists("dv_lvl_data")){dv_lvl_data<-join(dv_lvl_data, site_desc, by= "site_no")[, c(1,2,3,4,6,7,10,11,23)]}
if(exists("dv_lvl_data")){dv_lvl_data$Wat_elev<-dv_lvl_data$alt_va+dv_lvl_data$OBS}

if(exists("dv_cfs_data")) {dv_cfs_data$OBS<-log(dv_cfs_data$OBS)}

lvl_obs<-rbind(if(exists("gwl_data"))gwl_data, if(exists("dv_cfs_data"))dv_cfs_data, if(exists("dv_lvl_data"))dv_lvl_data)
lvl_obs$DateYM<-str_sub(lvl_obs$Date, 1, 7)
if(any (lvl_obs$OBS=="-999999")) lvl_obs$OBS[lvl_obs$OBS=="-999999"]<-NA  

#####Combine level data into one monthly table
obs_data_monthly<-aggregate(lvl_obs$OBS ~ lvl_obs$DateYM +lvl_obs$site_no+lvl_obs$Unit , lvl_obs,  mean)
colnames(obs_data_monthly)<-c("DateYM", "site_no","Unit", "OBS")
obs_data_monthly$DateYMD<-as.Date(paste(obs_data_monthly$DateYM, "-01", sep=""))



obs_data_monthly<-as.data.frame(obs_data_monthly%>%
                                  group_by(site_no)%>%
                                  mutate(z_score=(OBS-mean(OBS))/sd(OBS)))

# obs_data_monthly<-as.data.frame(obs_data_monthly%>%
#                                   group_by(site_no)%>%
#                                   filter(DateYMD > "2002-01-01") %>%
#                                   mutate(z_score=(OBS-mean(OBS))/sd(OBS)))

#####Create data for potential stencils derived from site data
mean_site_long<-mean(site_desc$dec_long_va)
mean_site_lat<- mean(site_desc$dec_lat_va)
huc8_stencil<-paste("HUC8::", paste(site_desc$huc_cd, collapse = ', '))


#####Get monthly PPT data for selected stencil - comment out unused stencils
#stencil<-webgeom("HUC8:: 07070002") #selected huc8
#stencil<-webgeom(huc8_stencil) #returns the average for all huc8s in the site descriptions
stencil<-simplegeom(c(mean_site_long, mean_site_lat)) #takes the mean of the lat/longs selected
#stencil<-simplegeom(c(-89.4, 44.4))
#stencil<-webgeom("state::Wisconsin")
ppt_job<-geoknife(stencil, (list(times = as.POSIXct(c('1895-01-01','2018-01-01')),url = 'http://cida.usgs.gov/thredds/dodsC/prism_v2', variables = 'ppt')), wait = TRUE)
ppt_data = result(ppt_job)[,c(1,2)] 
colnames(ppt_data)<-c("Date", "PPT_MM")
ppt_data$Date<-as.Date(as.POSIXct(ppt_data$Date))

#####calculate and add cumulative deviation from 10yr mean variables
ppt_data$mean10yr<- rollmean(ppt_data$PPT_MM, k=120, align = "right", fill = NA, na.rm=TRUE)
ppt_data$mo_dev10<-(ppt_data$PPT_MM-ppt_data$mean10yr)
ppt_data$mo_dev10[is.na(ppt_data$mo_dev10)]<-0
ppt_data$sum_dev10<-cumsum(ppt_data$mo_dev10)
ppt_data$site<-"Precip"
ppt_data$data$Unit<-"MM"

#####Set start time and graph variable
startTime<- as.Date(as.POSIXct("1895-01-01 00:00"))
ggplot(ppt_data, aes(ppt_data$Date, ppt_data$PPT_MM)) + geom_line() +(scale_x_date(limits = c(startTime, NA)))
ggplot(ppt_data, aes(ppt_data$Date, ppt_data$sum_dev10)) + geom_line() +(scale_x_date(limits = c(startTime, NA)))


#####Integrate PPT and levels data, create monthly mean data table
ppt_data$sum_dev10_Z<-((ppt_data$sum_dev10-mean(ppt_data$sum_dev10))/sd(ppt_data$sum_dev10))
ppt_z<-ppt_data[, c("PPT_MM","Date","sum_dev10_Z","site")]
lvl_z<-obs_data_monthly[, c("OBS","DateYMD","z_score","site_no")]
colnames(ppt_z)<-c("obs", "month", "z_score", "site")
colnames(lvl_z)<-c("obs", "month", "z_score", "site")
ppt_lvls_z<-rbind(ppt_z, lvl_z)
ppt_lvls_z$mo_no<-month(ppt_lvls_z$month)
ppt_lvls_mo_mean<-aggregate(obs~mo_no+site, ppt_lvls_z, mean)

ppt_lvls_z %>% 
  ggplot(aes(month, z_score, color = site))+
  geom_line(size=1)+
  scale_x_date(limits = as.Date(c('1995-01-01', '2017-01-01')))+
  labs(color = "Observation Type"
       , title = "Monthly Z Score\nCumulative Precip Deviation\nand Water Levels Variation")


ppt_lvls_z <- ppt_lvls_z %>%
  group_by(site) %>%
  mutate(z_1yr=rollmean(z_score, k=12, fill = NA, align = "right"))

ppt_lvls_z %>% 
  ggplot(aes(month, z_1yr, color = site))+
  geom_line(size=1)+
  scale_x_date(limits = as.Date(c('1995-01-01', '2018-01-01')))+
  labs(color = "Observation Type"
       , title = "1 Year Average Monthly Z Score\nCumulative Precip Deviation\nand Water Levels Variation")

ppt_lvls_z<- ppt_lvls_z%>%
  group_by(site)%>%
  filter(n()>60)

ppt_lvls_z <- ppt_lvls_z %>%
  group_by(site) %>%
  mutate(z_5yr=rollmean(z_score, k=60, fill = NA, align = "right"))

ppt_lvls_z %>% 
  ggplot(aes(month, z_5yr, color = site))+
  geom_line(size=1)+
  scale_x_date(limits = as.Date(c('1905-01-01', '2018-01-01')))+
  labs(color = "Observation Type"
       , title = "5 Year Average Monthly Z Score\nCumulative Precip Deviation\nand Water Levels Variation")

ppt_lvls_z<- ppt_lvls_z%>%
  group_by(site)%>%
  filter(n()>120)

ppt_lvls_z <- ppt_lvls_z %>%
  group_by(site) %>%
  mutate(z_10yr=rollmean(z_score, 120, na.pad=TRUE, align = "right"))

ppt_lvls_z %>% 
  ggplot(aes(month, z_10yr, color = site))+
  geom_line(size=1)+
  scale_x_date(limits = as.Date(c('1905-01-01', '2018-01-01')))+
  labs(color = "Observation Type"
       , title = "10 Year Average Monthly Z Score\nCumulative Precip Deviation\nand Water Levels Variation")

ppt_lvls_z <-ppt_lvls_z%>% filter(1==1)

obs_data_monthly %>% 
  ggplot(aes(DateYMD, OBS, color = site_no))+
  geom_line(size=1)+
  scale_x_date(limits = as.Date(c('1950-01-01', '2018-01-01')))
####################################################