In HARPgroup/cbp6:
#token = params$token
#riv.seg = params$riv.seg
# Loading Necessary Packages
library(dataRetrieval)
library(lubridate)
library(plyr)
library(zoo)
library(knitr)
library(ggplot2)
library(stats)
library(lfstat)
library(rstudioapi)
library(IHA)
library(PearsonDS)
library(lfstat)
library(scales)
library(readr)
library(httr)
library(dplyr)
library(stringr)
library(RCurl)
library(rgeos)
library(ggmap)
library(ggsn)
library(sp)
library(rlist)
library(tinytex)
options(tinytex.verbose = TRUE)
# # INPUTS
riv.seg <- 'PU3_4450_4440'
run.id1 <- '120'
run.id2 <- '122'
start.date <- '1984-01-01'
end.date <- '2000-12-31'
github_link <- "C:\\Users\\danie\\Documents\\HARP\\GitHub"
cbp6_link <- paste0(github_link, "/cbp6/code");
site_url <- "http://deq2.bse.vt.edu/d.dh"
site.or.server <- 'site'
site <- site_url
# obsolete inputs
mod.phase1 <- 'p6/p6_gb604' #or "p532c-sova" (phase 5)
mod.scenario1 <- 'CFBASE30Y20180615' #or "p532cal_062211" (phase 5)
mod.phase2 <- 'p6/p6_gb604' #or "p532c-sova" (phase 5)
mod.scenario2 <- 'CBASE1808L55CY55R45P50R45P50Y' #or "p532cal_062211" (phase 5)
site_number <- '00000000'
# GETTING MODEL DATA FROM VA HYDRO
# SETUP
# SETUP
# Setting active directory
# Setting working directory to the source file location
# Sourcing functions
#source(paste0(cbp6_link,"/cbp6_functions.R"))
#source(paste(github_link,"auth.private", sep = "/"));#load rest username and password, contained in
#token <- rest_token(site, token, rest_uname, rest_pw);
#options(timeout=1200); # set timeout to twice default level to avoid abort due to high traffic
source(paste0(cbp6_link,"/cbp6_functions.R"))
source(paste(github_link,"auth.private", sep = "/"));#load rest username and password, contained in
source(paste(cbp6_link, "/fn_vahydro-1.0.R", sep = ''))
token <- rest_token(site, token, rest_uname, rest_pw);
options(timeout=1200); # set timeout to twice default level to avoid abort due to high traffic
# Should new or original data be used?
new.or.original <- "new"
# CARRYOVER IF MASTER IS BEING RUN ----------------------------------------
if (exists("container.master") == TRUE) {
github_link <- container.master
riv.seg <- riv.seg.master
new.or.original <- new.or.original.master
}
# LOADING DATA ------------------------------------------------------------
if (site.or.server == 'site') {
#data1 <- model_import_data_cfs(riv.seg, mod.phase1, mod.scenario1, start.date, end.date)
data1 <- vahydro_import_data_cfs(riv.seg, run.id1, token, site = site_url, start.date, end.date)
#data2 <- model_import_data_cfs(riv.seg, mod.phase2, mod.scenario2, start.date, end.date)
data2 <- vahydro_import_data_cfs(riv.seg, run.id2, token, site = site_url, start.date, end.date)
} else if (site.or.server == 'server') {
data1 <- model_server_import_data_cfs(riv.seg, mod.phase1, mod.scenario1, start.date, end.date)
data2 <- model_server_import_data_cfs(riv.seg, mod.phase2, mod.scenario2, start.date, end.date)
}
data1 <- water_year_trim(data1)
data2 <- water_year_trim(data2)
title: "River Segment: r paste(riv.seg,"- VaHydro Run 120 (Base) vs. VAHydro Run 121 (CC)")"
header-includes:
- \usepackage{titling}
- \pretitle{\begin{flushleft}}
- \posttitle{\end{flushleft}}
output:
pdf_document
# #retrieve rest token
# source(paste0(github_link, "/auth.private"));
#
# #load rest username and password, contained in auth.private file
#
# token <- rest_token(site_url, token, rest_uname, rest_pw);
# options(timeout=120); # set timeout to twice default level to avoid abort due to high traffic
hydrocode = paste("vahydrosw_wshed_",riv.seg,sep="");
ftype = 'vahydro'; # nhd_huc8, nhd_huc10, vahydro
inputs <- list (
hydrocode = hydrocode,
bundle = 'watershed',
ftype = 'vahydro'
)
#property dataframe returned
feature = FALSE;
odata <- getFeature(inputs, token, site_url, feature);
hydroid <- odata[1,"hydroid"];
fname <- as.character(odata[1,]$name );
print(paste("Retrieved hydroid",hydroid,"for", fname,riv.seg, sep=' '));
# Getting the local drainage area feature
areainfo <- list(
varkey = "wshed_drainage_area_sqmi",
featureid = as.integer(as.character(hydroid)),
entity_type = "dh_feature"
)
model.area <- getPropertyALT(areainfo, site_url, model.area, token)
area <- model.area$propvalue
area <- area*27878400 #sq ft
# Loading written gage description
description <- try(read_file(paste0(cbp6_link, "/gage_descriptions/", site_number, ".txt")))
if (class(description) == "try-error") {
description <- ""
}
# Generating gage location maps
gis_img <- fn_gage_and_seg_mapperALT(riv.seg, site_number, site_url, cbp6_link, token)
ggsave("gis.png", plot = gis_img, device = 'png', width = 8, height = 5.5, units = 'in')
# Creating Data Frame with calculated metrics
metrics1 <- metrics_calc_all(data1)
metrics2 <- metrics_calc_all(data2)
percent_difference <- metrics_compare(metrics1, metrics2, riv.seg)
# initialize variables ----------------------------------------------------
# This section will create a hydrograph that will zoom in on 3 month segments where difference is high
# It does so for the top three highest difference periods
#all_data puts scenario flows and corresponding dates in one data frame
all_data <- data.frame(data1$date, data1$flow, data2$flow)
all_data$counter <- 1:length(all_data$data1.date) # counter fixes issues with row numbers later on in script
colnames(all_data) <- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
# find the first date for which data is collected, (in date format)
# and a date that is roughly one year and two months past the first date
YearStart <- as.character(as.Date(all_data$Date[1]))
fixer <- as.numeric(which(all_data$Date == paste0((year(YearStart)+1),"-11-30")))
YearEnd <- as.character(as.Date(all_data$Date[fixer]))
# YearStart_Row and YearEnd_Row are the rows corresponding the the YearStart and YearEnd dates
YearStart_Row <- as.numeric(which(all_data$Date==YearStart))
YearEnd_Row <- as.numeric(which(all_data$Date==YearEnd))
# initalize dataframes and counters, assign names for dataframe columns
#AvgMonthlyDifference: used within nested for loop to create a 1x12 matrix that holds 1 year of 3 month difference segments
#Timespan_Difference: used in large loop to store values from AvgMonthlyDifference; holds entire timespan of 3 month difference segments
AvgMonthlyDifference <- data.frame(matrix(nrow=1,ncol=1));
names(AvgMonthlyDifference)<-'Difference'
Timespan_Difference <- data.frame(matrix(nrow=1, ncol=1));
names(Timespan_Difference)<-'Difference'
i <- 1; # used for first for loop to advance a year
x <- 1 # x and y used to advance dataframes
y <- 12
# start loops used for yearly and monthly data -------------------------------------------------------------
loop <- as.numeric(round(length(data1$date)/365, digits = 0))-1
for (i in 1:loop){ # run loop for an entire data series
year <- all_data[YearStart_Row:YearEnd_Row,] # specify year: 10-01-year1 to 11-30-year2
m <- 1 # counter for nested loop
MonthStart <- YearStart # first date for 3 month timespan
doi <- as.Date(MonthStart)
doi <- doi + seq(0,365,31)
# doi= date of interest. dummy variable just to create the function next.month
next.month <- function(doi) as.Date(as.yearmon(doi) + 1/12) + as.numeric(as.Date(doi)-(as.Date(as.yearmon(doi))))
#(re)initalize variables for the nested loop
MonthEnd <-data.frame(next.month(doi)); # last date in 3 month timespan - used function to determine 3rd month
MonthEnd <- MonthEnd[3,1]-2 # specifies end of month 3 as last date
# (technically specifies 01 of month 4)
# row numbers corresponding with start and end dates, as a number. See note below
MonthStart_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthStart)))
MonthEnd_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthEnd)))
# Note: Counter column is used here to specify which row starts MonthStart and MonthEnd_Row.
# When rows are pulled from year row numbers are also pulled,
# so a counter must be used for proper row numbers.
Start_new <- as.numeric(which(year$Counter==MonthStart_Row))
End_new <- as.numeric(which(year$Counter==MonthEnd_Row))
# begin nested loop
for (m in 1:12){
month_time <- year[Start_new:End_new ,] #extract data for 3 month timespan within year of interest
avgmonth_scenario1 <- mean(month_time$`Scenario 1 Flow`) # find average of scenario 1 flow for 3 months
avgmonth_scenario2 <- mean(month_time$`Scenario 2 Flow`) # find average of scenario 2 flow for 3 months
AvgMonthlyDifference[m,1] <- (avgmonth_scenario1 - avgmonth_scenario2)/ avgmonth_scenario1 * 100 # percent difference between scenarios
MonthEnd<-as.Date(MonthEnd)
MonthEnd<-MonthEnd+1
MonthEndyear <- year(MonthEnd) # Year associated with last month of extracted data
MonthEndmonth <- month(MonthEnd) # Month associated with last month of extracted data
# the next three lines are for the difference calculations -- stop on 1st of month 4 (31 of month 3)
# Note: this DOES include the 1st of the next month in difference calculation
# Put a control on what date the script advances by - if date is not 1st of month, reset it
DateCheck <- as.Date(paste0(MonthEndyear,'-',MonthEndmonth,'-01'))
if (MonthEnd != DateCheck)
MonthEnd <- as.Date(paste0(MonthEndyear,'-', MonthEndmonth, '-01'))
MonthEnd <- MonthEnd-1
stopdate <- as.Date(MonthEnd)
AvgMonthlyDifference[m,2] <- stopdate
# Advance to next month or count
MonthStart <- next.month(MonthStart)
MonthEnd <- MonthEnd+1
MonthEnd <- next.month(MonthEnd)
MonthEnd <- MonthEnd-1
StartMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthStart));
StartMonth_Row <- as.numeric(StartMonth_Row)
EndMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthEnd));
EndMonth_Row <- as.numeric(EndMonth_Row)
Start_new <- which(year$Counter==StartMonth_Row)
End_new <- which(year$Counter==EndMonth_Row)
m <- m + 1
}
Timespan_Difference[x:y, 1] <- AvgMonthlyDifference[,1] # save the difference entries from AvgMonthlyDifference
Timespan_Difference[x:y, 2] <- AvgMonthlyDifference[,2] # save the dates
# advance Timespan_Difference for next run
x <- x + 12
y <- y + 12
YearStart <- as.Date(YearStart) + 365 # Advance 1 year
YearEnd <- as.Date(YearEnd) + 365 # Advance 1 year & 2 months (from 10-01 to 11-30)
# Put a control on what date the script advances by - if end date is not 11-30, reset it
# - if begin date is not -10-01, reset it
YearBeginyear <- year(YearStart) # pull year of beginning year
YearBeginCheck <- as.Date(paste0(YearBeginyear,'-10-01'))
if (YearBeginyear != YearBeginCheck)
YearStart <- as.Date(paste0(YearBeginyear,'-10-01'))
YearEndyear <- year(YearEnd) # pull year of ending year
YearEndCheck <- as.Date(paste0(YearEndyear,'-11-30'))
if (YearEnd != YearEndCheck)
YearEnd <- as.Date(paste0(YearEndyear,'-11-30'))
YearStart_Row <- which(as.Date(all_data$Date)== as.Date(YearStart))
YearEnd_Row <- which(as.Date(all_data$Date) == as.Date(YearEnd))
i <- i + 1
}
# This section of code will plot timeframes with high difference.
# count the number of 3 month periods over 20% difference, plot the highest 3 periods.
Timespan_Difference$Logic <- Timespan_Difference$Difference>=20 | Timespan_Difference$Difference<= -20
over20 <- Timespan_Difference[Timespan_Difference$Logic=='TRUE',]
HighDifference <- Timespan_Difference[order(abs(Timespan_Difference$Difference), decreasing = TRUE),]
names(HighDifference)<-c('Difference', 'Date', 'Logic')
# pull data for each of these 3 month segments.
HighestDifferences <- HighDifference[1:3,]
HighestDifferences$Date <- as.Date(HighestDifferences$Date)
# initalize variables for loop
differenceyear <- data.frame(matrix(nrow=1,ncol=6))
differencedates <- data.frame(matrix(nrow=1, ncol=2))
names(differenceyear)<- c('endyear', 'endmonth', 'enddate', 'startyear', 'startmonth', 'startdate')
names(differencedates)<- c('start date row', 'end date row')
storeplotdata1<- data.frame(matrix(nrow=1, ncol=4))
names(storeplotdata1)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
storeplotdata2<- data.frame(matrix(nrow=1, ncol=4))
names(storeplotdata2)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
storeplotdata3<- data.frame(matrix(nrow=1, ncol=4))
names(storeplotdata3)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
q <- 1
for (q in 1:length(HighestDifferences)){
differenceyear[q,1] <- year(HighestDifferences$Date[q]) # ending year
differenceyear[q,2]<- month(HighestDifferences$Date[q]) + 1 # ending month
differenceyear[q,4]<- year(HighestDifferences$Date[q]) #startyear
differenceyear[q,5]<- month(HighestDifferences$Date[q])-2 #startmonth
if (differenceyear[q,2] > 12) { # if end month is jan, must move year up
differenceyear[q,4] <- differenceyear[q,1]
differenceyear[q,1]<- differenceyear[q,1] + 1 # year for jan moves
differenceyear[q,2] <- 1
}else if (differenceyear[q,5] == -1) {
differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november
differenceyear[q,5] <- 11
}else if (differenceyear[q,5] == 0) {
differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november
differenceyear[q,5] <- 12
} else{
differenceyear[q,1]<- differenceyear[q,1] #endyear
differenceyear[q,2]<- differenceyear[q,2] #endmonth
differenceyear[q,4]<- differenceyear[q,4] #startyear
differenceyear[q,5]<- differenceyear[q,5] #startmonth
}
differenceyear[q,3]<- paste0(differenceyear[q,1], '-',differenceyear[q,2], '-01') #enddate
differenceyear$enddate <- as.Date(differenceyear$enddate)
differenceyear[q,6]<- as.Date(paste0(differenceyear[q,4], '-', differenceyear[q,5], '-01')) #startdate
differenceyear$startdate <- as.Date(differenceyear$startdate)
differencedates[q,1]<- as.character(differenceyear$startdate[q])
differencedates[q,2]<- as.character(differenceyear$enddate[q]-1)
differencedates[q,3]<- which(as.Date(all_data$Date)==as.Date(differencedates$`start date row`[q]))
differencedates[q,4]<- which(as.Date(all_data$Date)==as.Date(differencedates$`end date row`[q]))
plot1<-all_data[differencedates$V3[q]:differencedates$V4[q],]
if (q==1){
storeplotdata1<- plot1
}else if(q==2){
storeplotdata2<- plot1
}else if(q==3){
storeplotdata3<- plot1
}
q <- q+1
}
# # create and export 3 plots: \plot for info of row q
difference1 <- signif(HighestDifferences$Difference[1], digits=3) #Create difference variable to display on graph
difference2 <- signif(HighestDifferences$Difference[2], digits=3)
difference3 <- signif(HighestDifferences$Difference[3], digits=3)
# CREATES OUTPUT MATRIX -------------------------------------------------------
avg_scenario1 <- mean(data1$flow)
avg_scenario2 <- mean(data2$flow)
# also want to list the number of timespans that were over 20% difference.
over20 <- signif(nrow(over20)/nrow(Timespan_Difference)*100, digits=3)
OUTPUT_MATRIX <- matrix(c(avg_scenario1, avg_scenario2, over20), nrow=1, ncol=3)
rownames(OUTPUT_MATRIX) = c("Flow")
colnames(OUTPUT_MATRIX) = c('Scenario 1', 'Scenario 2', 'Difference>20 (%)')
overall_difference <- signif((OUTPUT_MATRIX[1,1]-OUTPUT_MATRIX[1,2])/OUTPUT_MATRIX[1,1]*100, digits=3)
OUTPUT_MATRIX <- matrix(c(over20, percent_difference[3,]), nrow=1, ncol=2)
rownames(OUTPUT_MATRIX) = c("Percent")
colnames(OUTPUT_MATRIX) = c('Difference > 20%', 'Overall Difference')
OUTPUT_MATRIX <- signif(as.numeric(OUTPUT_MATRIX, digits = 2))
OUTPUT_MATRIXsaver <- OUTPUT_MATRIX
OUTPUT_MATRIX <- kable(format(OUTPUT_MATRIX, digits = 3))
# plot for highest difference
# Max/min for y axis scaling
max <- max(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`)));
min <- min(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`)));
xpos1 <- min(storeplotdata1$Date)+20
xpos2 <- max(storeplotdata1$Date)-20
# Creating and exporting plot
df <- data.frame(as.Date(storeplotdata1$Date), storeplotdata1$`Scenario 1 Flow`, storeplotdata1$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
difference1plot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=14),
axis.text=element_text(size=14, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference1, '', '%'), size=3)+
annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '',
min(storeplotdata1$Date),': ', max(storeplotdata1$Date)), size=3)+
labs(y = "Flow (cfs)")
ggsave("fig6.png", plot = difference1plot, device = 'png', width = 8, height = 5.5, units = 'in')
# plot for second highest difference
# Max/min for y axis scaling
max <- max(c(max(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`)));
min <- min(c(max(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`)));
xpos1 <- min(storeplotdata2$Date)+20
xpos2 <- max(storeplotdata2$Date)-20
# Creating and exporting plot
df <- data.frame(as.Date(storeplotdata2$Date), storeplotdata2$`Scenario 1 Flow`, storeplotdata2$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
difference2plot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=14),
axis.text=element_text(size=14, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference2, '', '%'), size=3)+
annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '',
min(storeplotdata2$Date),': ', max(storeplotdata2$Date)), size=3)+
labs(y = "Flow (cfs)")
ggsave("fig7.png", plot = difference2plot, device = 'png', width = 8, height = 5.5, units = 'in')
# plot for third highest difference
# Max/min for y axis scaling
max <- max(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`)));
min <- min(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`)));
xpos1 <- min(storeplotdata3$Date)+20
xpos2 <- max(storeplotdata3$Date)-20
# Creating and exporting plot
df <- data.frame(as.Date(storeplotdata3$Date), storeplotdata3$`Scenario 1 Flow`, storeplotdata3$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
difference3plot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=14),
axis.text=element_text(size=14, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference3, '', '%'), size=3)+
annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '',
min(storeplotdata3$Date),': ', max(storeplotdata3$Date)), size=3)+
labs(y = "Flow (cfs)")
ggsave("fig8.png", plot = difference3plot, device = 'png', width = 8, height = 5.5, units = 'in')
# initialize variables ----------------------------------------------------
# This section will create a hydrograph that will zoom in on 3 month segments where difference is low
# It does so for the top three lowest difference periods
#all_data puts scenario flows and corresponding dates in one data frame
all_data <- data.frame(data1$date, data1$flow, data2$flow)
all_data$counter <- 1:length(all_data$data1.date) # counter fixes issues with row numbers later on in script
colnames(all_data) <- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
# find the first date for which data is collected, (in date format)
# and a date that is roughly one year and two months past the first date
YearStart <- as.character(as.Date(all_data$Date[1]))
fixer <- as.numeric(which(all_data$Date == paste0((year(YearStart)+1),"-11-30")))
YearEnd <- as.character(as.Date(all_data$Date[fixer]))
# YearStart_Row and YearEnd_Row are the rows corresponding the the YearStart and YearEnd dates
YearStart_Row <- as.numeric(which(all_data$Date==YearStart))
YearEnd_Row <- as.numeric(which(all_data$Date==YearEnd))
# initalize dataframes and counters, assign names for dataframe columns
#AvgMonthlyDifference: used within nested for loop to create a 1x12 matrix that holds 1 year of 3 month difference segments
#Timespan_Difference: used in large loop to store values from AvgMonthlyDifference; holds entire timespan of 3 month difference segments
AvgMonthlyDifference <- data.frame(matrix(nrow=1,ncol=1));
names(AvgMonthlyDifference)<-'Difference'
Timespan_Difference <- data.frame(matrix(nrow=1, ncol=1));
names(Timespan_Difference)<-'Difference'
i <- 1; # used for first for loop to advance a year
x <- 1 # x and y used to advance dataframes
y <- 12
# start loops used for yearly and monthly data -------------------------------------------------------------
loop <- as.numeric(round(length(data1$date)/365, digits = 0))-1
for (i in 1:loop){ # run loop for an entire data series
year <- all_data[YearStart_Row:YearEnd_Row,] # specify year: 10-01-year1 to 11-30-year2
m <- 1 # counter for nested loop
MonthStart <- YearStart # first date for 3 month timespan
doi <- as.Date(MonthStart)
doi <- doi + seq(0,365,31)
# doi= date of interest. dummy variable just to create the function next.month
next.month <- function(doi) as.Date(as.yearmon(doi) + 1/12) + as.numeric(as.Date(doi)-(as.Date(as.yearmon(doi))))
#(re)initalize variables for the nested loop
MonthEnd <-data.frame(next.month(doi)); # last date in 3 month timespan - used function to determine 3rd month
MonthEnd <- MonthEnd[3,1]-2 # specifies end of month 3 as last date
# (technically specifies 01 of month 4)
# row numbers corresponding with start and end dates, as a number. See note below
MonthStart_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthStart)))
MonthEnd_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthEnd)))
# Note: Counter column is used here to specify which row starts MonthStart and MonthEnd_Row.
# When rows are pulled from year row numbers are also pulled,
# so a counter must be used for proper row numbers.
Start_new <- as.numeric(which(year$Counter==MonthStart_Row))
End_new <- as.numeric(which(year$Counter==MonthEnd_Row))
# begin nested loop
for (m in 1:12){
month_time <- year[Start_new:End_new ,] #extract data for 3 month timespan within year of interest
avgmonth_scenario1 <- mean(month_time$`Scenario 1 Flow`) # find average of scenario 1 flow for 3 months
avgmonth_scenario2 <- mean(month_time$`Scenario 2 Flow`) # find average of scenario 2 flow for 3 months
AvgMonthlyDifference[m,1] <- (avgmonth_scenario1 - avgmonth_scenario2)/ avgmonth_scenario1 * 100 # percent difference between scenarios
MonthEnd<-as.Date(MonthEnd)
MonthEnd<-MonthEnd+1
MonthEndyear <- year(MonthEnd) # Year associated with last month of extracted data
MonthEndmonth <- month(MonthEnd) # Month associated with last month of extracted data
# the next three lines are for the difference calculations -- stop on 1st of month 4 (31 of month 3)
# Note: this DOES include the 1st of the next month in difference calculation
# Put a control on what date the script advances by - if date is not 1st of month, reset it
DateCheck <- as.Date(paste0(MonthEndyear,'-',MonthEndmonth,'-01'))
if (MonthEnd != DateCheck)
MonthEnd <- as.Date(paste0(MonthEndyear,'-', MonthEndmonth, '-01'))
MonthEnd <- MonthEnd-1
stopdate <- as.Date(MonthEnd)
AvgMonthlyDifference[m,2] <- stopdate
# Advance to next month or count
MonthStart <- next.month(MonthStart)
MonthEnd <- MonthEnd+1
MonthEnd <- next.month(MonthEnd)
MonthEnd <- MonthEnd-1
StartMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthStart));
StartMonth_Row <- as.numeric(StartMonth_Row)
EndMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthEnd));
EndMonth_Row <- as.numeric(EndMonth_Row)
Start_new <- which(year$Counter==StartMonth_Row)
End_new <- which(year$Counter==EndMonth_Row)
m <- m + 1
}
Timespan_Difference[x:y, 1] <- AvgMonthlyDifference[,1] # save the difference entries from AvgMonthlyDifference
Timespan_Difference[x:y, 2] <- AvgMonthlyDifference[,2] # save the dates
# advance Timespan_Difference for next run
x <- x + 12
y <- y + 12
YearStart <- as.Date(YearStart) + 365 # Advance 1 year
YearEnd <- as.Date(YearEnd) + 365 # Advance 1 year & 2 months (from 10-01 to 11-30)
# Put a control on what date the script advances by - if end date is not 11-30, reset it
# - if begin date is not -10-01, reset it
YearBeginyear <- year(YearStart) # pull year of beginning year
YearBeginCheck <- as.Date(paste0(YearBeginyear,'-10-01'))
if (YearBeginyear != YearBeginCheck)
YearStart <- as.Date(paste0(YearBeginyear,'-10-01'))
YearEndyear <- year(YearEnd) # pull year of ending year
YearEndCheck <- as.Date(paste0(YearEndyear,'-11-30'))
if (YearEnd != YearEndCheck)
YearEnd <- as.Date(paste0(YearEndyear,'-11-30'))
YearStart_Row <- which(as.Date(all_data$Date)== as.Date(YearStart))
YearEnd_Row <- which(as.Date(all_data$Date) == as.Date(YearEnd))
i <- i + 1
}
# This section of code will plot timeframes with high difference.
# count the number of 3 month periods over 20% difference, plot the highest 3 periods.
Timespan_Difference$Logic <- Timespan_Difference$Difference<=20 | Timespan_Difference$Difference>= -20
less20 <- Timespan_Difference[Timespan_Difference$Logic=='TRUE',]
HighDifference <- Timespan_Difference[order(abs(Timespan_Difference$Difference), decreasing = FALSE),]
names(HighDifference)<-c('Difference', 'Date', 'Logic')
# pull data for each of these 3 month segments.
HighestDifferences <- HighDifference[1:3,]
HighestDifferences$Date <- as.Date(HighestDifferences$Date)
# initalize variables for loop
differenceyear <- data.frame(matrix(nrow=1,ncol=6))
differencedates <- data.frame(matrix(nrow=1, ncol=2))
names(differenceyear)<- c('endyear', 'endmonth', 'enddate', 'startyear', 'startmonth', 'startdate')
names(differencedates)<- c('start date row', 'end date row')
storeplotdata1<- data.frame(matrix(nrow=1, ncol=4))
names(storeplotdata1)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
storeplotdata2<- data.frame(matrix(nrow=1, ncol=4))
names(storeplotdata2)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
storeplotdata3<- data.frame(matrix(nrow=1, ncol=4))
names(storeplotdata3)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter')
q <- 1
for (q in 1:length(HighestDifferences)){
differenceyear[q,1] <- year(HighestDifferences$Date[q]) # ending year
differenceyear[q,2]<- month(HighestDifferences$Date[q]) + 1 # ending month
differenceyear[q,4]<- year(HighestDifferences$Date[q]) #startyear
differenceyear[q,5]<- month(HighestDifferences$Date[q])-2 #startmonth
if (differenceyear[q,2] > 12) { # if end month is jan, must move year up
differenceyear[q,4] <- differenceyear[q,1]
differenceyear[q,1]<- differenceyear[q,1] + 1 # year for jan moves
differenceyear[q,2] <- 1
}else if (differenceyear[q,5] == -1) {
differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november
differenceyear[q,5] <- 11
}else if (differenceyear[q,5] == 0) {
differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november
differenceyear[q,5] <- 12
} else{
differenceyear[q,1]<- differenceyear[q,1] #endyear
differenceyear[q,2]<- differenceyear[q,2] #endmonth
differenceyear[q,4]<- differenceyear[q,4] #startyear
differenceyear[q,5]<- differenceyear[q,5] #startmonth
}
differenceyear[q,3]<- paste0(differenceyear[q,1], '-',differenceyear[q,2], '-01') #enddate
differenceyear$enddate <- as.Date(differenceyear$enddate)
differenceyear[q,6]<- as.Date(paste0(differenceyear[q,4], '-', differenceyear[q,5], '-01')) #startdate
differenceyear$startdate <- as.Date(differenceyear$startdate)
differencedates[q,1]<- as.character(differenceyear$startdate[q])
differencedates[q,2]<- as.character(differenceyear$enddate[q]-1)
differencedates[q,3]<- which(as.Date(all_data$Date)==as.Date(differencedates$`start date row`[q]))
differencedates[q,4]<- which(as.Date(all_data$Date)==as.Date(differencedates$`end date row`[q]))
plot1<-all_data[differencedates$V3[q]:differencedates$V4[q],]
if (q==1){
storeplotdata1<- plot1
}else if(q==2){
storeplotdata2<- plot1
}else if(q==3){
storeplotdata3<- plot1
}
q <- q+1
}
# # create and export 3 plots: \plot for info of row q
difference1 <- signif(HighestDifferences$Difference[1], digits=3) #Create difference variable to display on graph
difference2 <- signif(HighestDifferences$Difference[2], digits=3)
difference3 <- signif(HighestDifferences$Difference[3], digits=3)
# CREATES OUTPUT MATRIX -------------------------------------------------------
avg_scenario1 <- mean(data1$flow)
avg_scenario2 <- mean(data2$flow)
# also want to list the number of timespans that were over 20% difference.
less20 <- signif(nrow(less20)/nrow(Timespan_Difference)*100, digits=3)
OUTPUT_MATRIX_LESS <- matrix(c(avg_scenario1, avg_scenario2, less20), nrow=1, ncol=3)
rownames(OUTPUT_MATRIX_LESS) = c("Flow")
colnames(OUTPUT_MATRIX_LESS) = c('Scenario 1', 'Scenario 2', 'Difference<20 (%)')
overall_difference <- signif((OUTPUT_MATRIX_LESS[1,1]-OUTPUT_MATRIX_LESS[1,2])/OUTPUT_MATRIX_LESS[1,1]*100, digits=3)
OUTPUT_MATRIX_LESS <- matrix(c(less20, percent_difference[3,]), nrow=1, ncol=2)
rownames(OUTPUT_MATRIX_LESS) = c("Percent")
colnames(OUTPUT_MATRIX_LESS) = c('Difference < 20%', 'Overall Difference')
OUTPUT_MATRIX_LESS <- signif(as.numeric(OUTPUT_MATRIX_LESS, digits = 2))
OUTPUT_MATRIX_LESSsaver <- OUTPUT_MATRIX_LESS
OUTPUT_MATRIX_LESS <- kable(format(OUTPUT_MATRIX_LESS, digits = 3))
# plot for highest difference
# Max/min for y axis scaling
max <- max(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`)));
min <- min(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`)));
xpos1 <- min(storeplotdata1$Date)+20
xpos2 <- max(storeplotdata1$Date)-20
# Creating and exporting plot
df <- data.frame(as.Date(storeplotdata1$Date), storeplotdata1$`Scenario 1 Flow`, storeplotdata1$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
difference1plot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=14),
axis.text=element_text(size=14, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference1, '', '%'), size=3)+
annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '',
min(storeplotdata1$Date),': ', max(storeplotdata1$Date)), size=3)+
labs(y = "Flow (cfs)")
ggsave("fig11.png", plot = difference1plot, device = 'png', width = 8, height = 5.5, units = 'in')
# plot for second highest difference
# Max/min for y axis scaling
max <- max(c(max(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`)));
min <- max(c(min(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`)));
xpos1 <- min(storeplotdata2$Date)+20
xpos2 <- max(storeplotdata2$Date)-20
# Creating and exporting plot
df <- data.frame(as.Date(storeplotdata2$Date), storeplotdata2$`Scenario 1 Flow`, storeplotdata2$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
difference2plot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=14),
axis.text=element_text(size=14, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference2, '', '%'), size=3)+
annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '',
min(storeplotdata2$Date),': ', max(storeplotdata2$Date)), size=3)+
labs(y = "Flow (cfs)")
ggsave("fig12.png", plot = difference2plot, device = 'png', width = 8, height = 5.5, units = 'in')
# plot for third highest difference
# Max/min for y axis scaling
max <- max(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`)));
min <- min(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`)));
xpos1 <- min(storeplotdata3$Date)+20
xpos2 <- max(storeplotdata3$Date)-20
# Creating and exporting plot
df <- data.frame(as.Date(storeplotdata3$Date), storeplotdata3$`Scenario 1 Flow`, storeplotdata3$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
difference3plot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=14),
axis.text=element_text(size=14, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference3, '', '%'), size=3)+
annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '',
min(storeplotdata3$Date),': ', max(storeplotdata3$Date)), size=3)+
labs(y = "Flow (cfs)")
ggsave("fig13.png", plot = difference3plot, device = 'png', width = 8, height = 5.5, units = 'in')
# OUTPUT TABLES -----
# Table 1: Monthly Average Flow
Table1 <- matrix(c(percent_difference[1,1], percent_difference[1,14], percent_difference[1,15], percent_difference[1,16], percent_difference[1,17], percent_difference[1,18], percent_difference[1,19], percent_difference[1,20], percent_difference[1,21], percent_difference[1,22], percent_difference[1,23], percent_difference[1,24], percent_difference[1,25], percent_difference[2,1], percent_difference[2,14], percent_difference[2,15], percent_difference[2,16], percent_difference[2,17], percent_difference[2,18], percent_difference[2,19], percent_difference[2,20], percent_difference[2,21], percent_difference[2,22], percent_difference[2,23], percent_difference[2,24], percent_difference[2,25], percent_difference[3,1], percent_difference[3,14], percent_difference[3,15], percent_difference[3,16], percent_difference[3,17], percent_difference[3,18], percent_difference[3,19], percent_difference[3,20], percent_difference[3,21], percent_difference[3,22], percent_difference[3,23], percent_difference[3,24], percent_difference[3,25]),
nrow = 13, ncol = 3);
colnames(Table1) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference");
rownames(Table1) = c("Overall Mean Flow",
"Jan. Mean Flow", "Feb. Mean Flow",
"Mar. Mean Flow", "Apr. Mean Flow",
"May Mean Flow", "Jun. Mean Flow",
"Jul. Mean Flow", "Aug. Mean Flow",
"Sep. Mean Flow", "Oct. Mean Flow",
"Nov. Mean Flow", "Dec. Mean Flow");
Table1 <- signif(Table1, digits = 3)
Table1 <- round(Table1, digits = 2)
Table1 <- kable(format(Table1, digits = 3, drop0trailing = TRUE))
Table2 <- matrix(c(percent_difference[1,2], percent_difference[1,3],
percent_difference[1,4], percent_difference[1,5], percent_difference[1,6],
percent_difference[1,7], percent_difference[1,8], percent_difference[1,9],
percent_difference[1,10], percent_difference[1,11], percent_difference[1,12],
percent_difference[1,13], percent_difference[2,2],
percent_difference[2,3], percent_difference[2,4], percent_difference[2,5],
percent_difference[2,6], percent_difference[2,7], percent_difference[2,8],
percent_difference[2,9], percent_difference[2,10], percent_difference[2,11],
percent_difference[2,12], percent_difference[2,13],
percent_difference[3,2], percent_difference[3,3], percent_difference[3,4],
percent_difference[3,5], percent_difference[3,6], percent_difference[3,7],
percent_difference[3,8], percent_difference[3,9], percent_difference[3,10],
percent_difference[3,11], percent_difference[3,12], percent_difference[3,13]),
nrow = 12, ncol = 3);
colnames(Table2) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference");
rownames(Table2) = c("Jan. Low Flow", "Feb. Low Flow",
"Mar. Low Flow", "Apr. Low Flow",
"May Low Flow", "Jun. Low Flow",
"Jul. Low Flow", "Aug. Low Flow",
"Sep. Low Flow", "Oct. Low Flow",
"Nov. Low Flow", "Dec. Low Flow");
Table2 <- signif(Table2, digits = 3)
Table2 <- round(Table2, digits = 2)
Table2 <- kable(format(Table2, digits = 3, drop0trailing = TRUE))
# Table 3: Monthly High Flow
Table3 <- matrix(c(percent_difference[1,26], percent_difference[1,27], percent_difference[1,28],percent_difference[1,29], percent_difference[1,30], percent_difference[1,31],percent_difference[1,32], percent_difference[1,33], percent_difference[1,34],percent_difference[1,35], percent_difference[1,36], percent_difference[1,37],percent_difference[2,26],percent_difference[2,27], percent_difference[2,28], percent_difference[2,29],percent_difference[2,30], percent_difference[2,31], percent_difference[2,32],percent_difference[2,33], percent_difference[2,34], percent_difference[2,35],percent_difference[2,36], percent_difference[2,37],percent_difference[3,26], percent_difference[3,27], percent_difference[3,28],percent_difference[3,29], percent_difference[3,30], percent_difference[3,31],percent_difference[3,32], percent_difference[3,33], percent_difference[3,34],percent_difference[3,35], percent_difference[3,36], percent_difference[3,37]),
nrow = 12, ncol = 3);
colnames(Table3) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference");
rownames(Table3) = c("Jan. High Flow", "Feb. High Flow",
"Mar. High Flow", "Apr. High Flow",
"May High Flow", "Jun. High Flow",
"Jul. High Flow", "Aug. High Flow",
"Sep. High Flow", "Oct. High Flow",
"Nov. High Flow", "Dec. High Flow");
Table3 <- signif(Table3, digits = 3)
Table3 <- round(Table3, digits = 2)
Table3 <- kable(format(Table3, digits = 3, drop0trailing = TRUE))
Table4 <- matrix(c(percent_difference[1,38], percent_difference[1,43],
percent_difference[1,39], percent_difference[1,44],
percent_difference[1,40], percent_difference[1,45],
percent_difference[1,41], percent_difference[1,46],
percent_difference[1,42], percent_difference[1,47],
percent_difference[1,59], percent_difference[1,60],
percent_difference[1,58], percent_difference[1,67],
percent_difference[2,38], percent_difference[2,43],
percent_difference[2,39], percent_difference[2,44],
percent_difference[2,40], percent_difference[2,45],
percent_difference[2,41], percent_difference[2,46],
percent_difference[2,42], percent_difference[2,47],
percent_difference[2,59], percent_difference[2,60],
percent_difference[2,58], percent_difference[2,67],
percent_difference[3,38], percent_difference[3,43],
percent_difference[3,39], percent_difference[3,44],
percent_difference[3,40], percent_difference[3,45],
percent_difference[3,41], percent_difference[3,46],
percent_difference[3,42], percent_difference[3,47],
percent_difference[3,59], percent_difference[3,60],
percent_difference[3,58], percent_difference[3,67]),
nrow = 14, ncol = 3);
colnames(Table4) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference");
rownames(Table4) = c("Min. 1 Day Min", "Med. 1 Day Min",
"Min. 3 Day Min", "Med. 3 Day Min",
"Min. 7 Day Min", "Med. 7 Day Min",
"Min. 30 Day Min", "Med. 30 Day Min",
"Min. 90 Day Min", "Med. 90 Day Min",
"7Q10", "Year of 90-Day Min. Flow",
"Drought Year Mean", "Mean Baseflow");
Table4 <- signif(Table4, digits = 3)
Table4 <- round(Table4, digits = 2)
Table4 <- kable(format(Table4, digits = 3, drop0trailing = TRUE))
# Table 5: Period High Flows
Table5 <- matrix(c(percent_difference[1,48], percent_difference[1,53],
percent_difference[1,49], percent_difference[1,54],
percent_difference[1,50], percent_difference[1,55],
percent_difference[1,51], percent_difference[1,56],
percent_difference[1,52], percent_difference[1,57],
percent_difference[2,48], percent_difference[2,53],
percent_difference[2,49], percent_difference[2,54],
percent_difference[2,50], percent_difference[2,55],
percent_difference[2,51], percent_difference[2,56],
percent_difference[2,52], percent_difference[2,57],
percent_difference[3,48], percent_difference[3,53],
percent_difference[3,49], percent_difference[3,54],
percent_difference[3,50], percent_difference[3,55],
percent_difference[3,51], percent_difference[3,56],
percent_difference[3,52], percent_difference[3,57]),
nrow = 10, ncol = 3);
colnames(Table5) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference");
rownames(Table5) = c("Max. 1 Day Max", "Med. 1 Day Max",
"Max. 3 Day Max", "Med. 3 Day Max",
"Max. 7 Day Max", "Med. 7 Day Max",
"Max. 30 Day Max", "Med. 30 Day Max",
"Max. 90 Day Max", "Med. 90 Day Max");
Table5 <- signif(Table5, digits = 3)
Table5 <- round(Table5, digits = 2)
Table5 <- kable(format(Table5, digits = 3, drop0trailing = TRUE))
# Table 6: Non-Exceedance Flows
Table6 <- matrix(c(percent_difference[1,62], percent_difference[1,63], percent_difference[1,64],percent_difference[1,65], percent_difference[1,66], percent_difference[1,61],percent_difference[2,62], percent_difference[2,63], percent_difference[2,64],percent_difference[2,65], percent_difference[2,66], percent_difference[2,61],percent_difference[3,62], percent_difference[3,63], percent_difference[3,64],percent_difference[3,65], percent_difference[3,66], percent_difference[3,61]),
nrow = 6, ncol = 3);
colnames(Table6) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference");
rownames(Table6) = c("1% Non-Exceedance", "5% Non-Exceedance",
"50% Non-Exceedance", "95% Non-Exceedance",
"99% Non-Exceedance", "Sept. 10% Non-Exceedance");
Table6 <- signif(Table6, digits = 3)
Table6 <- round(Table6, digits = 2)
Table6 <- kable(format(Table6, digits = 3, drop0trailing = TRUE))
# Creating names for plot legends
name_scenario1 <- paste('Scenario 1', mod.scenario1);
name_scenario2 <- paste('Scenario 2', mod.scenario2);
# CREATE FUNCTIONS FOR PLOTTING ------------------------------------------
base_breaks <- function(n = 10){
function(x) {
axisTicks(log10(range(x, na.rm = TRUE)), log = TRUE, n = n)
}
}
scaleFUN <- function(x) sprintf("%.0f", x)
include_graphics("gis.png")
description.cont <- paste0(" The average daily discharge change between scenario 1 and scenario 2 for the 20 year timespan was ", OUTPUT_MATRIXsaver[2] ,"%, with ", OUTPUT_MATRIXsaver[1], "% of its rolling three month time spans above 20% difference.")
cat(description)
cat(description.cont)
\pagebreak
Table 1: Monthly Low Flows
Table2
Table 2: Monthly Average Flows
Table1
\pagebreak
Table 3: Monthly High Flows
Table3
Table 4: Period Low Flows
Table4
Table 5: Period High Flows
Table5
Table 6: Non-Exceedance Flows
Table6
# SETTING UP PLOTS
# Basic hydrograph -----
# Max/min for y axis scaling
max <- max(c(max(all_data$`Scenario 1 Flow`), max(all_data$`Scenario 2 Flow`)));
min <- min(c(min(all_data$`Scenario 1 Flow`), min(all_data$`Scenario 2 Flow`)));
if (max > 10000){
max <- 100000
}else if (max > 1000){
max <- 10000
}else if (max > 100){
max <- 1000
}else if (max > 10){
max <- 100
}
if (min>100){
min<-100
}else if (min>10){
min<-10
}else
min<-1
if (min==100){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(100, 1000, 10000, 100000),
limits=c(min,max))
}else if (min==10){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(10, 100, 1000, 10000),
limits=c(min,max))
}else if (min==1){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(1, 10, 100, 1000, 10000),
limits=c(min,max))
}else
fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(),
labels=scaleFUN, limits=c(min,max))
df <- data.frame(as.Date(all_data$Date), all_data$`Scenario 1 Flow`, all_data$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=0.5) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=0.5)+
fixtheyscale+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(y = "Flow (cfs)")
ggsave("fig1.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
# Zoomed hydrograph in year of lowest 90-year flow -----
# Running scenario 1 calculations
f3_scenario1 <- zoo(all_data$`Scenario 1 Flow`, order.by = all_data$Date)
g2_scenario1 <- group2(f3_scenario1, year = 'water')
# Running scenairo 2 calculations
f3_scenario2 <- zoo(all_data$`Scenario 2 Flow`, order.by = all_data$Date)
g2_scenario2 <- group2(f3_scenario2, year = 'water')
yearly_scenario1_90DayMin <- g2_scenario1[,c(1,10)];
yearly_scenario2_90DayMin <- g2_scenario2[,c(1,10)];
low.year <- subset(yearly_scenario1_90DayMin, yearly_scenario1_90DayMin$`90 Day Min`==min(yearly_scenario1_90DayMin$`90 Day Min`));
low.year <- low.year$year;
low.year <- subset(all_data, year(all_data$Date)==low.year);
# Scaling using max/min
max <- max(c(max(low.year$`Scenario 1 Flow`), max(low.year$`Scenario 2 Flow`)));
min <- min(c(min(low.year$`Scenario 1 Flow`), min(low.year$`Scenario 2 Flow`)));
if (max > 10000){
max <- 100000
}else if (max > 1000){
max <- 10000
}else if (max > 100){
max <- 1000
}else if (max > 10){
max <- 100
}
if (min>100){
min<-100
}else if (min>10){
min<-10
}else
min<-1
if (min==100){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(100, 1000, 10000, 100000),
limits=c(min,max))
}else if (min==10){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(10, 100, 1000, 10000),
limits=c(min,max))
}else if (min==1){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(1, 10, 100, 1000, 10000),
limits=c(min,max))
}else
fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(),
labels=scaleFUN, limits=c(min,max))
df <- data.frame(as.Date(low.year$Date), low.year$`Scenario 1 Flow`, low.year$`Scenario 2 Flow`);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=0.7) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=0.7)+
fixtheyscale+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=11, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(y = "Flow (cfs)")
ggsave("fig2.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
#Flow exceedance plot -----
# Determining the "rank" (0-1) of the flow value
num_observations <- as.numeric(length(all_data$Date))
rank_vec <- as.numeric(c(1:num_observations))
# Calculating exceedance probability
prob_exceedance <- 100*((rank_vec) / (num_observations + 1))
exceed_scenario1 <- sort(all_data$`Scenario 1 Flow`, decreasing = TRUE)
exceed_scenario2 <- sort(all_data$`Scenario 2 Flow`, decreasing = TRUE)
scenario1_exceedance <- quantile(exceed_scenario1, probs = c(0.01, 0.05, 0.5, 0.95, 0.99))
scenario2_exceedance <- quantile(exceed_scenario2, probs = c(0.01, 0.05, 0.5, 0.95, 0.99))
# Determining max flow value for exceedance plot scale
max <- max(c(max(scenario1_exceedance), max(scenario2_exceedance)));
min <- min(c(min(scenario1_exceedance), min(scenario2_exceedance)));
if (max > 10000){
max <- 100000
}else if (max > 1000){
max <- 10000
}else if (max > 100){
max <- 1000
}else if (max > 10){
max <- 100
}
if (min>100){
min<-100
}else if (min>10){
min<-10
}else
min<-1
if (min==100){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(100, 1000, 10000, 100000),
limits=c(min,max))
}else if (min==10){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(10, 100, 1000, 10000),
limits=c(min,max))
}else if (min==1){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(1, 10, 100, 1000, 10000),
limits=c(min,max))
}else
fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(),
labels=scaleFUN, limits=c(min,max))
df <- data.frame(prob_exceedance, exceed_scenario1, exceed_scenario2);
colnames(df) <- c('Date', 'Scenario1', 'Scenario2')
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=0.8) +
geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=0.8)+
fixtheyscale+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(x= "Probability of Exceedance (%)", y = "Flow (cfs)")
ggsave("fig3.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
# Baseflow Indiviudal Graph -----
data1$year <- year(ymd(data1$date))
data1$month <- month(ymd(data1$date))
data1$day <- day(ymd(data1$date))
data2$year <- year(ymd(data2$date))
data2$month <- month(ymd(data2$date))
data2$day <- day(ymd(data2$date))
scenario1river <- createlfobj(data1, hyearstart = 10, baseflow = TRUE, meta = NULL)
scenario2river <- createlfobj(data2, hyearstart = 10, baseflow = TRUE, meta = NULL)
baseflowriver<- data.frame(scenario1river, scenario2river);
colnames(baseflowriver) <-c ('mday', 'mmonth', 'myear', 'mflow', 'mHyear', 'mBaseflow',
'gday', 'gmonth', 'gyear', 'gflow', 'gHyear', 'gBaseflow')
# removing NA values
baseflowriver<-baseflowriver[complete.cases(baseflowriver)==TRUE,]
scenario2river<- data.frame(baseflowriver$gday, baseflowriver$gmonth, baseflowriver$gyear,
baseflowriver$gflow, baseflowriver$gHyear, baseflowriver$gBaseflow);
scenario1river<- data.frame(baseflowriver$mday, baseflowriver$mmonth, baseflowriver$myear,
baseflowriver$mflow, baseflowriver$mHyear, baseflowriver$mBaseflow)
names(scenario1river) <- c('day', 'month', 'year', 'flow', 'hyear', 'baseflow')
names(scenario2river) <- c('day', 'month', 'year', 'flow', 'hyear', 'baseflow')
# Adding date vectors
scenario1river$date <- as.Date(paste0(scenario1river$year,"-",scenario1river$month,"-",scenario1river$day))
scenario2river$date <- as.Date(paste0(scenario2river$year,"-",scenario2river$month,"-",scenario2river$day))
# Determining max flow value for plot scale
max <- max(c(max(scenario1river$baseflow), max(scenario2river$baseflow), max(scenario1river$flow), max(scenario2river$flow)));
min <- min(c(min(scenario1river$baseflow), min(scenario2river$baseflow), min(scenario1river$flow), min(scenario2river$flow)));
if (max > 10000){
max <- 100000
}else if (max > 1000){
max <- 10000
}else if (max > 100){
max <- 1000
}else if (max > 10){
max <- 100
}
if (min>100){
min<-100
}else if (min>10){
min<-10
}else
min<-1
if (min==100){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(100, 1000, 10000, 100000),
limits=c(min,max))
}else if (min==10){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(10, 100, 1000, 10000),
limits=c(min,max))
}else if (min==1){
fixtheyscale<- scale_y_continuous(trans = log_trans(),
breaks = c(1, 10, 100, 1000, 10000),
limits=c(min,max))
}else
fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(),
labels=scaleFUN, limits=c(min,max))
df <- data.frame(as.Date(scenario1river$date), scenario1river$baseflow, scenario2river$baseflow, scenario1river$flow, scenario2river$flow);
colnames(df) <- c('Date', 'Scenario1Baseflow', 'Scenario2Baseflow','Scenario1Flow', 'Scenario2Flow')
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1Baseflow, color="VAHydro Scen. 1"), size=0.5) +
geom_line(aes(y=Scenario2Baseflow, color="VAHydro Scen. 2"), size=0.5)+
fixtheyscale+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(y = "Flow (cfs)")
ggsave("fig4.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
par(mfrow = c(1,1));
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_line(aes(y=Scenario1Flow, color=paste("Scen. 1 Flow")), size=1)+
geom_line(aes(y=Scenario2Flow, color=paste("Scen. 2 Flow")), size=0.5)+
geom_line(aes(y=Scenario1Baseflow, color=paste("Scen. 1 Baseflow")), size=0.5) +
geom_line(aes(y=Scenario2Baseflow, color=paste("Scen. 2 Baseflow")), size=1)+
fixtheyscale+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("black","red","grey", "light pink"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(y = "Flow (cfs)")
ggsave("fig5.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
# Setup for Residuals
data <- all_data[complete.cases(all_data),]
resid <- (data$`Scenario 2 Flow` - data$`Scenario 1 Flow`)
resid <- data.frame(data$Date, resid)
# Residuals plot for hydrograph
zeroline <- rep_len(0, length(data$Date))
quantresid <- data.frame(signif(quantile(resid$resid, na.rm = TRUE), digits=3))
min <- min(resid$resid)
max <- max(resid$resid)
names(quantresid) <- c('Percentiles')
df <- data.frame(as.Date(resid$data.Date), resid$resid, zeroline);
colnames(df) <- c('Date', 'Residual', 'Zeroline')
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_point(aes(y=Residual, color="Residual (CC - Base)"), size=1) +
geom_line(aes(y=Zeroline, color="Zeroline"), size=0.8)+
scale_y_continuous(limits=c(min,max))+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("dark green","black"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(y = "Flow Difference (cfs)")
ggsave("fig9A.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
# Setup for Residuals
data <- all_data[complete.cases(all_data),]
#data_weighted <- data/area
resid <- 1000000*((data$`Scenario 2 Flow` - data$`Scenario 1 Flow`)/area)
resid <- data.frame(data$Date, resid)
# Residuals plot for hydrograph
zeroline <- rep_len(0, length(data$Date))
quantresid <- data.frame(signif(quantile(resid$resid, na.rm =TRUE), digits=3))
min <- min(resid$resid)
max <- max(resid$resid)
names(quantresid) <- c('Percentiles')
df <- data.frame(as.Date(resid$data.Date), resid$resid, zeroline);
colnames(df) <- c('Date', 'Residual', 'Zeroline')
options(scipen=5, width = 1400, height = 950)
myplot <- ggplot(df, aes(x=Date)) +
geom_point(aes(y=Residual, color="Residual (CC - Base)"), size=1) +
geom_line(aes(y=Zeroline, color="Zeroline"), size=0.8)+
scale_y_continuous(limits=c(min,max))+
theme_bw()+
theme(legend.position="top",
legend.title=element_blank(),
legend.box = "horizontal",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="white"),
legend.text=element_text(size=12),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=14, colour="black"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.ticks = element_line(colour="black"),
panel.grid.major=element_line(colour = "light grey"),
panel.grid.minor=element_blank())+
scale_colour_manual(values=c("dark green","black"))+
guides(colour = guide_legend(override.aes = list(size=5)))+
labs(y = "Area Weighted Flow Difference*10^6 (ft/s)")
ggsave("fig9B.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
Fig. 1: Hydrograph
include_graphics("fig1.png")
Fig. 2: Zoomed Hydrograph
include_graphics("fig2.png")
Fig. 3: Flow Exceedance
include_graphics("fig3.png")
Fig. 4: Baseflow
include_graphics("fig4.png")
Fig. 5: Combined Baseflow
include_graphics("fig5.png")
Fig. 6: Largest Difference Segment
include_graphics("fig6.png")
Fig. 7: Second Largest Difference Segment
include_graphics("fig7.png")
Fig. 8: Third Largest Difference Segment
include_graphics("fig8.png")
Fig. 9A: Residuals Plot
include_graphics("fig9A.png")
Fig. 9B: Area Weighted Residuals Plot
include_graphics("fig9B.png")
dat <- vahydro_import_local.runoff.inflows_cfs(riv.seg, run.id1, token, site, start.date, end.date);
dat <- subset(dat, dat$date >= start.date & dat$date <= end.date);
boxplot(as.numeric(dat$flow.unit) ~ year(dat$date), outline = FALSE, ylab = 'Runit Flow (cfs)', xlab = 'Date')
dev.copy(png, 'fig10.png')
dev.off()
Fig. 10: VA Hydro Scen. 1 Runit Values (Outliers Excluded)
include_graphics("fig10.png")
Fig. 11: Smallest Difference Segment
include_graphics("fig11.png")
Fig. 12: Second Smallest Difference Segment
include_graphics("fig12.png")
Fig. 13: Third Smallest Difference Segment
include_graphics("fig13.png")
HARPgroup/cbp6 documentation built on Oct. 14, 2024, 4:19 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#token = params$token #riv.seg = params$riv.seg
# Loading Necessary Packages library(dataRetrieval) library(lubridate) library(plyr) library(zoo) library(knitr) library(ggplot2) library(stats) library(lfstat) library(rstudioapi) library(IHA) library(PearsonDS) library(lfstat) library(scales) library(readr) library(httr) library(dplyr) library(stringr) library(RCurl) library(rgeos) library(ggmap) library(ggsn) library(sp) library(rlist) library(tinytex) options(tinytex.verbose = TRUE) # # INPUTS riv.seg <- 'PU3_4450_4440' run.id1 <- '120' run.id2 <- '122' start.date <- '1984-01-01' end.date <- '2000-12-31' github_link <- "C:\\Users\\danie\\Documents\\HARP\\GitHub" cbp6_link <- paste0(github_link, "/cbp6/code"); site_url <- "http://deq2.bse.vt.edu/d.dh" site.or.server <- 'site' site <- site_url # obsolete inputs mod.phase1 <- 'p6/p6_gb604' #or "p532c-sova" (phase 5) mod.scenario1 <- 'CFBASE30Y20180615' #or "p532cal_062211" (phase 5) mod.phase2 <- 'p6/p6_gb604' #or "p532c-sova" (phase 5) mod.scenario2 <- 'CBASE1808L55CY55R45P50R45P50Y' #or "p532cal_062211" (phase 5) site_number <- '00000000' # GETTING MODEL DATA FROM VA HYDRO # SETUP # SETUP # Setting active directory # Setting working directory to the source file location
# Sourcing functions #source(paste0(cbp6_link,"/cbp6_functions.R")) #source(paste(github_link,"auth.private", sep = "/"));#load rest username and password, contained in #token <- rest_token(site, token, rest_uname, rest_pw); #options(timeout=1200); # set timeout to twice default level to avoid abort due to high traffic
source(paste0(cbp6_link,"/cbp6_functions.R")) source(paste(github_link,"auth.private", sep = "/"));#load rest username and password, contained in source(paste(cbp6_link, "/fn_vahydro-1.0.R", sep = '')) token <- rest_token(site, token, rest_uname, rest_pw); options(timeout=1200); # set timeout to twice default level to avoid abort due to high traffic # Should new or original data be used? new.or.original <- "new" # CARRYOVER IF MASTER IS BEING RUN ---------------------------------------- if (exists("container.master") == TRUE) { github_link <- container.master riv.seg <- riv.seg.master new.or.original <- new.or.original.master } # LOADING DATA ------------------------------------------------------------ if (site.or.server == 'site') { #data1 <- model_import_data_cfs(riv.seg, mod.phase1, mod.scenario1, start.date, end.date) data1 <- vahydro_import_data_cfs(riv.seg, run.id1, token, site = site_url, start.date, end.date) #data2 <- model_import_data_cfs(riv.seg, mod.phase2, mod.scenario2, start.date, end.date) data2 <- vahydro_import_data_cfs(riv.seg, run.id2, token, site = site_url, start.date, end.date) } else if (site.or.server == 'server') { data1 <- model_server_import_data_cfs(riv.seg, mod.phase1, mod.scenario1, start.date, end.date) data2 <- model_server_import_data_cfs(riv.seg, mod.phase2, mod.scenario2, start.date, end.date) } data1 <- water_year_trim(data1) data2 <- water_year_trim(data2)
title: "River Segment: r paste(riv.seg,"- VaHydro Run 120 (Base) vs. VAHydro Run 121 (CC)")"
header-includes:
- \usepackage{titling}
- \pretitle{\begin{flushleft}}
- \posttitle{\end{flushleft}}
output:
pdf_document
# #retrieve rest token # source(paste0(github_link, "/auth.private")); # # #load rest username and password, contained in auth.private file # # token <- rest_token(site_url, token, rest_uname, rest_pw); # options(timeout=120); # set timeout to twice default level to avoid abort due to high traffic
hydrocode = paste("vahydrosw_wshed_",riv.seg,sep=""); ftype = 'vahydro'; # nhd_huc8, nhd_huc10, vahydro inputs <- list ( hydrocode = hydrocode, bundle = 'watershed', ftype = 'vahydro' ) #property dataframe returned feature = FALSE; odata <- getFeature(inputs, token, site_url, feature); hydroid <- odata[1,"hydroid"]; fname <- as.character(odata[1,]$name ); print(paste("Retrieved hydroid",hydroid,"for", fname,riv.seg, sep=' '));
# Getting the local drainage area feature areainfo <- list( varkey = "wshed_drainage_area_sqmi", featureid = as.integer(as.character(hydroid)), entity_type = "dh_feature" ) model.area <- getPropertyALT(areainfo, site_url, model.area, token) area <- model.area$propvalue area <- area*27878400 #sq ft # Loading written gage description description <- try(read_file(paste0(cbp6_link, "/gage_descriptions/", site_number, ".txt"))) if (class(description) == "try-error") { description <- "" } # Generating gage location maps gis_img <- fn_gage_and_seg_mapperALT(riv.seg, site_number, site_url, cbp6_link, token) ggsave("gis.png", plot = gis_img, device = 'png', width = 8, height = 5.5, units = 'in')
# Creating Data Frame with calculated metrics metrics1 <- metrics_calc_all(data1) metrics2 <- metrics_calc_all(data2) percent_difference <- metrics_compare(metrics1, metrics2, riv.seg)
# initialize variables ---------------------------------------------------- # This section will create a hydrograph that will zoom in on 3 month segments where difference is high # It does so for the top three highest difference periods #all_data puts scenario flows and corresponding dates in one data frame all_data <- data.frame(data1$date, data1$flow, data2$flow) all_data$counter <- 1:length(all_data$data1.date) # counter fixes issues with row numbers later on in script colnames(all_data) <- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') # find the first date for which data is collected, (in date format) # and a date that is roughly one year and two months past the first date YearStart <- as.character(as.Date(all_data$Date[1])) fixer <- as.numeric(which(all_data$Date == paste0((year(YearStart)+1),"-11-30"))) YearEnd <- as.character(as.Date(all_data$Date[fixer])) # YearStart_Row and YearEnd_Row are the rows corresponding the the YearStart and YearEnd dates YearStart_Row <- as.numeric(which(all_data$Date==YearStart)) YearEnd_Row <- as.numeric(which(all_data$Date==YearEnd)) # initalize dataframes and counters, assign names for dataframe columns #AvgMonthlyDifference: used within nested for loop to create a 1x12 matrix that holds 1 year of 3 month difference segments #Timespan_Difference: used in large loop to store values from AvgMonthlyDifference; holds entire timespan of 3 month difference segments AvgMonthlyDifference <- data.frame(matrix(nrow=1,ncol=1)); names(AvgMonthlyDifference)<-'Difference' Timespan_Difference <- data.frame(matrix(nrow=1, ncol=1)); names(Timespan_Difference)<-'Difference' i <- 1; # used for first for loop to advance a year x <- 1 # x and y used to advance dataframes y <- 12 # start loops used for yearly and monthly data ------------------------------------------------------------- loop <- as.numeric(round(length(data1$date)/365, digits = 0))-1 for (i in 1:loop){ # run loop for an entire data series year <- all_data[YearStart_Row:YearEnd_Row,] # specify year: 10-01-year1 to 11-30-year2 m <- 1 # counter for nested loop MonthStart <- YearStart # first date for 3 month timespan doi <- as.Date(MonthStart) doi <- doi + seq(0,365,31) # doi= date of interest. dummy variable just to create the function next.month next.month <- function(doi) as.Date(as.yearmon(doi) + 1/12) + as.numeric(as.Date(doi)-(as.Date(as.yearmon(doi)))) #(re)initalize variables for the nested loop MonthEnd <-data.frame(next.month(doi)); # last date in 3 month timespan - used function to determine 3rd month MonthEnd <- MonthEnd[3,1]-2 # specifies end of month 3 as last date # (technically specifies 01 of month 4) # row numbers corresponding with start and end dates, as a number. See note below MonthStart_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthStart))) MonthEnd_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthEnd))) # Note: Counter column is used here to specify which row starts MonthStart and MonthEnd_Row. # When rows are pulled from year row numbers are also pulled, # so a counter must be used for proper row numbers. Start_new <- as.numeric(which(year$Counter==MonthStart_Row)) End_new <- as.numeric(which(year$Counter==MonthEnd_Row)) # begin nested loop for (m in 1:12){ month_time <- year[Start_new:End_new ,] #extract data for 3 month timespan within year of interest avgmonth_scenario1 <- mean(month_time$`Scenario 1 Flow`) # find average of scenario 1 flow for 3 months avgmonth_scenario2 <- mean(month_time$`Scenario 2 Flow`) # find average of scenario 2 flow for 3 months AvgMonthlyDifference[m,1] <- (avgmonth_scenario1 - avgmonth_scenario2)/ avgmonth_scenario1 * 100 # percent difference between scenarios MonthEnd<-as.Date(MonthEnd) MonthEnd<-MonthEnd+1 MonthEndyear <- year(MonthEnd) # Year associated with last month of extracted data MonthEndmonth <- month(MonthEnd) # Month associated with last month of extracted data # the next three lines are for the difference calculations -- stop on 1st of month 4 (31 of month 3) # Note: this DOES include the 1st of the next month in difference calculation # Put a control on what date the script advances by - if date is not 1st of month, reset it DateCheck <- as.Date(paste0(MonthEndyear,'-',MonthEndmonth,'-01')) if (MonthEnd != DateCheck) MonthEnd <- as.Date(paste0(MonthEndyear,'-', MonthEndmonth, '-01')) MonthEnd <- MonthEnd-1 stopdate <- as.Date(MonthEnd) AvgMonthlyDifference[m,2] <- stopdate # Advance to next month or count MonthStart <- next.month(MonthStart) MonthEnd <- MonthEnd+1 MonthEnd <- next.month(MonthEnd) MonthEnd <- MonthEnd-1 StartMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthStart)); StartMonth_Row <- as.numeric(StartMonth_Row) EndMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthEnd)); EndMonth_Row <- as.numeric(EndMonth_Row) Start_new <- which(year$Counter==StartMonth_Row) End_new <- which(year$Counter==EndMonth_Row) m <- m + 1 } Timespan_Difference[x:y, 1] <- AvgMonthlyDifference[,1] # save the difference entries from AvgMonthlyDifference Timespan_Difference[x:y, 2] <- AvgMonthlyDifference[,2] # save the dates # advance Timespan_Difference for next run x <- x + 12 y <- y + 12 YearStart <- as.Date(YearStart) + 365 # Advance 1 year YearEnd <- as.Date(YearEnd) + 365 # Advance 1 year & 2 months (from 10-01 to 11-30) # Put a control on what date the script advances by - if end date is not 11-30, reset it # - if begin date is not -10-01, reset it YearBeginyear <- year(YearStart) # pull year of beginning year YearBeginCheck <- as.Date(paste0(YearBeginyear,'-10-01')) if (YearBeginyear != YearBeginCheck) YearStart <- as.Date(paste0(YearBeginyear,'-10-01')) YearEndyear <- year(YearEnd) # pull year of ending year YearEndCheck <- as.Date(paste0(YearEndyear,'-11-30')) if (YearEnd != YearEndCheck) YearEnd <- as.Date(paste0(YearEndyear,'-11-30')) YearStart_Row <- which(as.Date(all_data$Date)== as.Date(YearStart)) YearEnd_Row <- which(as.Date(all_data$Date) == as.Date(YearEnd)) i <- i + 1 } # This section of code will plot timeframes with high difference. # count the number of 3 month periods over 20% difference, plot the highest 3 periods. Timespan_Difference$Logic <- Timespan_Difference$Difference>=20 | Timespan_Difference$Difference<= -20 over20 <- Timespan_Difference[Timespan_Difference$Logic=='TRUE',] HighDifference <- Timespan_Difference[order(abs(Timespan_Difference$Difference), decreasing = TRUE),] names(HighDifference)<-c('Difference', 'Date', 'Logic') # pull data for each of these 3 month segments. HighestDifferences <- HighDifference[1:3,] HighestDifferences$Date <- as.Date(HighestDifferences$Date) # initalize variables for loop differenceyear <- data.frame(matrix(nrow=1,ncol=6)) differencedates <- data.frame(matrix(nrow=1, ncol=2)) names(differenceyear)<- c('endyear', 'endmonth', 'enddate', 'startyear', 'startmonth', 'startdate') names(differencedates)<- c('start date row', 'end date row') storeplotdata1<- data.frame(matrix(nrow=1, ncol=4)) names(storeplotdata1)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') storeplotdata2<- data.frame(matrix(nrow=1, ncol=4)) names(storeplotdata2)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') storeplotdata3<- data.frame(matrix(nrow=1, ncol=4)) names(storeplotdata3)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') q <- 1 for (q in 1:length(HighestDifferences)){ differenceyear[q,1] <- year(HighestDifferences$Date[q]) # ending year differenceyear[q,2]<- month(HighestDifferences$Date[q]) + 1 # ending month differenceyear[q,4]<- year(HighestDifferences$Date[q]) #startyear differenceyear[q,5]<- month(HighestDifferences$Date[q])-2 #startmonth if (differenceyear[q,2] > 12) { # if end month is jan, must move year up differenceyear[q,4] <- differenceyear[q,1] differenceyear[q,1]<- differenceyear[q,1] + 1 # year for jan moves differenceyear[q,2] <- 1 }else if (differenceyear[q,5] == -1) { differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november differenceyear[q,5] <- 11 }else if (differenceyear[q,5] == 0) { differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november differenceyear[q,5] <- 12 } else{ differenceyear[q,1]<- differenceyear[q,1] #endyear differenceyear[q,2]<- differenceyear[q,2] #endmonth differenceyear[q,4]<- differenceyear[q,4] #startyear differenceyear[q,5]<- differenceyear[q,5] #startmonth } differenceyear[q,3]<- paste0(differenceyear[q,1], '-',differenceyear[q,2], '-01') #enddate differenceyear$enddate <- as.Date(differenceyear$enddate) differenceyear[q,6]<- as.Date(paste0(differenceyear[q,4], '-', differenceyear[q,5], '-01')) #startdate differenceyear$startdate <- as.Date(differenceyear$startdate) differencedates[q,1]<- as.character(differenceyear$startdate[q]) differencedates[q,2]<- as.character(differenceyear$enddate[q]-1) differencedates[q,3]<- which(as.Date(all_data$Date)==as.Date(differencedates$`start date row`[q])) differencedates[q,4]<- which(as.Date(all_data$Date)==as.Date(differencedates$`end date row`[q])) plot1<-all_data[differencedates$V3[q]:differencedates$V4[q],] if (q==1){ storeplotdata1<- plot1 }else if(q==2){ storeplotdata2<- plot1 }else if(q==3){ storeplotdata3<- plot1 } q <- q+1 } # # create and export 3 plots: \plot for info of row q difference1 <- signif(HighestDifferences$Difference[1], digits=3) #Create difference variable to display on graph difference2 <- signif(HighestDifferences$Difference[2], digits=3) difference3 <- signif(HighestDifferences$Difference[3], digits=3) # CREATES OUTPUT MATRIX ------------------------------------------------------- avg_scenario1 <- mean(data1$flow) avg_scenario2 <- mean(data2$flow) # also want to list the number of timespans that were over 20% difference. over20 <- signif(nrow(over20)/nrow(Timespan_Difference)*100, digits=3) OUTPUT_MATRIX <- matrix(c(avg_scenario1, avg_scenario2, over20), nrow=1, ncol=3) rownames(OUTPUT_MATRIX) = c("Flow") colnames(OUTPUT_MATRIX) = c('Scenario 1', 'Scenario 2', 'Difference>20 (%)') overall_difference <- signif((OUTPUT_MATRIX[1,1]-OUTPUT_MATRIX[1,2])/OUTPUT_MATRIX[1,1]*100, digits=3) OUTPUT_MATRIX <- matrix(c(over20, percent_difference[3,]), nrow=1, ncol=2) rownames(OUTPUT_MATRIX) = c("Percent") colnames(OUTPUT_MATRIX) = c('Difference > 20%', 'Overall Difference') OUTPUT_MATRIX <- signif(as.numeric(OUTPUT_MATRIX, digits = 2)) OUTPUT_MATRIXsaver <- OUTPUT_MATRIX OUTPUT_MATRIX <- kable(format(OUTPUT_MATRIX, digits = 3)) # plot for highest difference # Max/min for y axis scaling max <- max(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`))); min <- min(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`))); xpos1 <- min(storeplotdata1$Date)+20 xpos2 <- max(storeplotdata1$Date)-20 # Creating and exporting plot df <- data.frame(as.Date(storeplotdata1$Date), storeplotdata1$`Scenario 1 Flow`, storeplotdata1$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) difference1plot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=14), axis.text=element_text(size=14, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference1, '', '%'), size=3)+ annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '', min(storeplotdata1$Date),': ', max(storeplotdata1$Date)), size=3)+ labs(y = "Flow (cfs)") ggsave("fig6.png", plot = difference1plot, device = 'png', width = 8, height = 5.5, units = 'in') # plot for second highest difference # Max/min for y axis scaling max <- max(c(max(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`))); min <- min(c(max(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`))); xpos1 <- min(storeplotdata2$Date)+20 xpos2 <- max(storeplotdata2$Date)-20 # Creating and exporting plot df <- data.frame(as.Date(storeplotdata2$Date), storeplotdata2$`Scenario 1 Flow`, storeplotdata2$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) difference2plot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=14), axis.text=element_text(size=14, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference2, '', '%'), size=3)+ annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '', min(storeplotdata2$Date),': ', max(storeplotdata2$Date)), size=3)+ labs(y = "Flow (cfs)") ggsave("fig7.png", plot = difference2plot, device = 'png', width = 8, height = 5.5, units = 'in') # plot for third highest difference # Max/min for y axis scaling max <- max(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`))); min <- min(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`))); xpos1 <- min(storeplotdata3$Date)+20 xpos2 <- max(storeplotdata3$Date)-20 # Creating and exporting plot df <- data.frame(as.Date(storeplotdata3$Date), storeplotdata3$`Scenario 1 Flow`, storeplotdata3$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) difference3plot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=14), axis.text=element_text(size=14, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference3, '', '%'), size=3)+ annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '', min(storeplotdata3$Date),': ', max(storeplotdata3$Date)), size=3)+ labs(y = "Flow (cfs)") ggsave("fig8.png", plot = difference3plot, device = 'png', width = 8, height = 5.5, units = 'in')
# initialize variables ---------------------------------------------------- # This section will create a hydrograph that will zoom in on 3 month segments where difference is low # It does so for the top three lowest difference periods #all_data puts scenario flows and corresponding dates in one data frame all_data <- data.frame(data1$date, data1$flow, data2$flow) all_data$counter <- 1:length(all_data$data1.date) # counter fixes issues with row numbers later on in script colnames(all_data) <- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') # find the first date for which data is collected, (in date format) # and a date that is roughly one year and two months past the first date YearStart <- as.character(as.Date(all_data$Date[1])) fixer <- as.numeric(which(all_data$Date == paste0((year(YearStart)+1),"-11-30"))) YearEnd <- as.character(as.Date(all_data$Date[fixer])) # YearStart_Row and YearEnd_Row are the rows corresponding the the YearStart and YearEnd dates YearStart_Row <- as.numeric(which(all_data$Date==YearStart)) YearEnd_Row <- as.numeric(which(all_data$Date==YearEnd)) # initalize dataframes and counters, assign names for dataframe columns #AvgMonthlyDifference: used within nested for loop to create a 1x12 matrix that holds 1 year of 3 month difference segments #Timespan_Difference: used in large loop to store values from AvgMonthlyDifference; holds entire timespan of 3 month difference segments AvgMonthlyDifference <- data.frame(matrix(nrow=1,ncol=1)); names(AvgMonthlyDifference)<-'Difference' Timespan_Difference <- data.frame(matrix(nrow=1, ncol=1)); names(Timespan_Difference)<-'Difference' i <- 1; # used for first for loop to advance a year x <- 1 # x and y used to advance dataframes y <- 12 # start loops used for yearly and monthly data ------------------------------------------------------------- loop <- as.numeric(round(length(data1$date)/365, digits = 0))-1 for (i in 1:loop){ # run loop for an entire data series year <- all_data[YearStart_Row:YearEnd_Row,] # specify year: 10-01-year1 to 11-30-year2 m <- 1 # counter for nested loop MonthStart <- YearStart # first date for 3 month timespan doi <- as.Date(MonthStart) doi <- doi + seq(0,365,31) # doi= date of interest. dummy variable just to create the function next.month next.month <- function(doi) as.Date(as.yearmon(doi) + 1/12) + as.numeric(as.Date(doi)-(as.Date(as.yearmon(doi)))) #(re)initalize variables for the nested loop MonthEnd <-data.frame(next.month(doi)); # last date in 3 month timespan - used function to determine 3rd month MonthEnd <- MonthEnd[3,1]-2 # specifies end of month 3 as last date # (technically specifies 01 of month 4) # row numbers corresponding with start and end dates, as a number. See note below MonthStart_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthStart))) MonthEnd_Row <- as.numeric(which(as.Date(all_data$Date)==as.Date(MonthEnd))) # Note: Counter column is used here to specify which row starts MonthStart and MonthEnd_Row. # When rows are pulled from year row numbers are also pulled, # so a counter must be used for proper row numbers. Start_new <- as.numeric(which(year$Counter==MonthStart_Row)) End_new <- as.numeric(which(year$Counter==MonthEnd_Row)) # begin nested loop for (m in 1:12){ month_time <- year[Start_new:End_new ,] #extract data for 3 month timespan within year of interest avgmonth_scenario1 <- mean(month_time$`Scenario 1 Flow`) # find average of scenario 1 flow for 3 months avgmonth_scenario2 <- mean(month_time$`Scenario 2 Flow`) # find average of scenario 2 flow for 3 months AvgMonthlyDifference[m,1] <- (avgmonth_scenario1 - avgmonth_scenario2)/ avgmonth_scenario1 * 100 # percent difference between scenarios MonthEnd<-as.Date(MonthEnd) MonthEnd<-MonthEnd+1 MonthEndyear <- year(MonthEnd) # Year associated with last month of extracted data MonthEndmonth <- month(MonthEnd) # Month associated with last month of extracted data # the next three lines are for the difference calculations -- stop on 1st of month 4 (31 of month 3) # Note: this DOES include the 1st of the next month in difference calculation # Put a control on what date the script advances by - if date is not 1st of month, reset it DateCheck <- as.Date(paste0(MonthEndyear,'-',MonthEndmonth,'-01')) if (MonthEnd != DateCheck) MonthEnd <- as.Date(paste0(MonthEndyear,'-', MonthEndmonth, '-01')) MonthEnd <- MonthEnd-1 stopdate <- as.Date(MonthEnd) AvgMonthlyDifference[m,2] <- stopdate # Advance to next month or count MonthStart <- next.month(MonthStart) MonthEnd <- MonthEnd+1 MonthEnd <- next.month(MonthEnd) MonthEnd <- MonthEnd-1 StartMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthStart)); StartMonth_Row <- as.numeric(StartMonth_Row) EndMonth_Row <- which(as.Date(all_data$Date)==as.Date(MonthEnd)); EndMonth_Row <- as.numeric(EndMonth_Row) Start_new <- which(year$Counter==StartMonth_Row) End_new <- which(year$Counter==EndMonth_Row) m <- m + 1 } Timespan_Difference[x:y, 1] <- AvgMonthlyDifference[,1] # save the difference entries from AvgMonthlyDifference Timespan_Difference[x:y, 2] <- AvgMonthlyDifference[,2] # save the dates # advance Timespan_Difference for next run x <- x + 12 y <- y + 12 YearStart <- as.Date(YearStart) + 365 # Advance 1 year YearEnd <- as.Date(YearEnd) + 365 # Advance 1 year & 2 months (from 10-01 to 11-30) # Put a control on what date the script advances by - if end date is not 11-30, reset it # - if begin date is not -10-01, reset it YearBeginyear <- year(YearStart) # pull year of beginning year YearBeginCheck <- as.Date(paste0(YearBeginyear,'-10-01')) if (YearBeginyear != YearBeginCheck) YearStart <- as.Date(paste0(YearBeginyear,'-10-01')) YearEndyear <- year(YearEnd) # pull year of ending year YearEndCheck <- as.Date(paste0(YearEndyear,'-11-30')) if (YearEnd != YearEndCheck) YearEnd <- as.Date(paste0(YearEndyear,'-11-30')) YearStart_Row <- which(as.Date(all_data$Date)== as.Date(YearStart)) YearEnd_Row <- which(as.Date(all_data$Date) == as.Date(YearEnd)) i <- i + 1 } # This section of code will plot timeframes with high difference. # count the number of 3 month periods over 20% difference, plot the highest 3 periods. Timespan_Difference$Logic <- Timespan_Difference$Difference<=20 | Timespan_Difference$Difference>= -20 less20 <- Timespan_Difference[Timespan_Difference$Logic=='TRUE',] HighDifference <- Timespan_Difference[order(abs(Timespan_Difference$Difference), decreasing = FALSE),] names(HighDifference)<-c('Difference', 'Date', 'Logic') # pull data for each of these 3 month segments. HighestDifferences <- HighDifference[1:3,] HighestDifferences$Date <- as.Date(HighestDifferences$Date) # initalize variables for loop differenceyear <- data.frame(matrix(nrow=1,ncol=6)) differencedates <- data.frame(matrix(nrow=1, ncol=2)) names(differenceyear)<- c('endyear', 'endmonth', 'enddate', 'startyear', 'startmonth', 'startdate') names(differencedates)<- c('start date row', 'end date row') storeplotdata1<- data.frame(matrix(nrow=1, ncol=4)) names(storeplotdata1)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') storeplotdata2<- data.frame(matrix(nrow=1, ncol=4)) names(storeplotdata2)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') storeplotdata3<- data.frame(matrix(nrow=1, ncol=4)) names(storeplotdata3)<- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') q <- 1 for (q in 1:length(HighestDifferences)){ differenceyear[q,1] <- year(HighestDifferences$Date[q]) # ending year differenceyear[q,2]<- month(HighestDifferences$Date[q]) + 1 # ending month differenceyear[q,4]<- year(HighestDifferences$Date[q]) #startyear differenceyear[q,5]<- month(HighestDifferences$Date[q])-2 #startmonth if (differenceyear[q,2] > 12) { # if end month is jan, must move year up differenceyear[q,4] <- differenceyear[q,1] differenceyear[q,1]<- differenceyear[q,1] + 1 # year for jan moves differenceyear[q,2] <- 1 }else if (differenceyear[q,5] == -1) { differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november differenceyear[q,5] <- 11 }else if (differenceyear[q,5] == 0) { differenceyear[q,4] <- differenceyear[q,4] - 1 # if january, go back a year and start november differenceyear[q,5] <- 12 } else{ differenceyear[q,1]<- differenceyear[q,1] #endyear differenceyear[q,2]<- differenceyear[q,2] #endmonth differenceyear[q,4]<- differenceyear[q,4] #startyear differenceyear[q,5]<- differenceyear[q,5] #startmonth } differenceyear[q,3]<- paste0(differenceyear[q,1], '-',differenceyear[q,2], '-01') #enddate differenceyear$enddate <- as.Date(differenceyear$enddate) differenceyear[q,6]<- as.Date(paste0(differenceyear[q,4], '-', differenceyear[q,5], '-01')) #startdate differenceyear$startdate <- as.Date(differenceyear$startdate) differencedates[q,1]<- as.character(differenceyear$startdate[q]) differencedates[q,2]<- as.character(differenceyear$enddate[q]-1) differencedates[q,3]<- which(as.Date(all_data$Date)==as.Date(differencedates$`start date row`[q])) differencedates[q,4]<- which(as.Date(all_data$Date)==as.Date(differencedates$`end date row`[q])) plot1<-all_data[differencedates$V3[q]:differencedates$V4[q],] if (q==1){ storeplotdata1<- plot1 }else if(q==2){ storeplotdata2<- plot1 }else if(q==3){ storeplotdata3<- plot1 } q <- q+1 } # # create and export 3 plots: \plot for info of row q difference1 <- signif(HighestDifferences$Difference[1], digits=3) #Create difference variable to display on graph difference2 <- signif(HighestDifferences$Difference[2], digits=3) difference3 <- signif(HighestDifferences$Difference[3], digits=3) # CREATES OUTPUT MATRIX ------------------------------------------------------- avg_scenario1 <- mean(data1$flow) avg_scenario2 <- mean(data2$flow) # also want to list the number of timespans that were over 20% difference. less20 <- signif(nrow(less20)/nrow(Timespan_Difference)*100, digits=3) OUTPUT_MATRIX_LESS <- matrix(c(avg_scenario1, avg_scenario2, less20), nrow=1, ncol=3) rownames(OUTPUT_MATRIX_LESS) = c("Flow") colnames(OUTPUT_MATRIX_LESS) = c('Scenario 1', 'Scenario 2', 'Difference<20 (%)') overall_difference <- signif((OUTPUT_MATRIX_LESS[1,1]-OUTPUT_MATRIX_LESS[1,2])/OUTPUT_MATRIX_LESS[1,1]*100, digits=3) OUTPUT_MATRIX_LESS <- matrix(c(less20, percent_difference[3,]), nrow=1, ncol=2) rownames(OUTPUT_MATRIX_LESS) = c("Percent") colnames(OUTPUT_MATRIX_LESS) = c('Difference < 20%', 'Overall Difference') OUTPUT_MATRIX_LESS <- signif(as.numeric(OUTPUT_MATRIX_LESS, digits = 2)) OUTPUT_MATRIX_LESSsaver <- OUTPUT_MATRIX_LESS OUTPUT_MATRIX_LESS <- kable(format(OUTPUT_MATRIX_LESS, digits = 3)) # plot for highest difference # Max/min for y axis scaling max <- max(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`))); min <- min(c(max(storeplotdata1$`Scenario 1 Flow`), max(storeplotdata1$`Scenario 2 Flow`))); xpos1 <- min(storeplotdata1$Date)+20 xpos2 <- max(storeplotdata1$Date)-20 # Creating and exporting plot df <- data.frame(as.Date(storeplotdata1$Date), storeplotdata1$`Scenario 1 Flow`, storeplotdata1$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) difference1plot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=14), axis.text=element_text(size=14, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference1, '', '%'), size=3)+ annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '', min(storeplotdata1$Date),': ', max(storeplotdata1$Date)), size=3)+ labs(y = "Flow (cfs)") ggsave("fig11.png", plot = difference1plot, device = 'png', width = 8, height = 5.5, units = 'in') # plot for second highest difference # Max/min for y axis scaling max <- max(c(max(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`))); min <- max(c(min(storeplotdata2$`Scenario 1 Flow`), max(storeplotdata2$`Scenario 2 Flow`))); xpos1 <- min(storeplotdata2$Date)+20 xpos2 <- max(storeplotdata2$Date)-20 # Creating and exporting plot df <- data.frame(as.Date(storeplotdata2$Date), storeplotdata2$`Scenario 1 Flow`, storeplotdata2$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) difference2plot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=14), axis.text=element_text(size=14, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference2, '', '%'), size=3)+ annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '', min(storeplotdata2$Date),': ', max(storeplotdata2$Date)), size=3)+ labs(y = "Flow (cfs)") ggsave("fig12.png", plot = difference2plot, device = 'png', width = 8, height = 5.5, units = 'in') # plot for third highest difference # Max/min for y axis scaling max <- max(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`))); min <- min(c(max(storeplotdata3$`Scenario 1 Flow`), max(storeplotdata3$`Scenario 2 Flow`))); xpos1 <- min(storeplotdata3$Date)+20 xpos2 <- max(storeplotdata3$Date)-20 # Creating and exporting plot df <- data.frame(as.Date(storeplotdata3$Date), storeplotdata3$`Scenario 1 Flow`, storeplotdata3$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) difference3plot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=1)+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=14), axis.text=element_text(size=14, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ annotate("text", x=xpos1, y=1.05*max, label= paste0('Difference:', '',difference3, '', '%'), size=3)+ annotate("text", x=xpos2, y=1.05*max, label= paste0('Date Range: ', '', min(storeplotdata3$Date),': ', max(storeplotdata3$Date)), size=3)+ labs(y = "Flow (cfs)") ggsave("fig13.png", plot = difference3plot, device = 'png', width = 8, height = 5.5, units = 'in')
# OUTPUT TABLES ----- # Table 1: Monthly Average Flow Table1 <- matrix(c(percent_difference[1,1], percent_difference[1,14], percent_difference[1,15], percent_difference[1,16], percent_difference[1,17], percent_difference[1,18], percent_difference[1,19], percent_difference[1,20], percent_difference[1,21], percent_difference[1,22], percent_difference[1,23], percent_difference[1,24], percent_difference[1,25], percent_difference[2,1], percent_difference[2,14], percent_difference[2,15], percent_difference[2,16], percent_difference[2,17], percent_difference[2,18], percent_difference[2,19], percent_difference[2,20], percent_difference[2,21], percent_difference[2,22], percent_difference[2,23], percent_difference[2,24], percent_difference[2,25], percent_difference[3,1], percent_difference[3,14], percent_difference[3,15], percent_difference[3,16], percent_difference[3,17], percent_difference[3,18], percent_difference[3,19], percent_difference[3,20], percent_difference[3,21], percent_difference[3,22], percent_difference[3,23], percent_difference[3,24], percent_difference[3,25]), nrow = 13, ncol = 3); colnames(Table1) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference"); rownames(Table1) = c("Overall Mean Flow", "Jan. Mean Flow", "Feb. Mean Flow", "Mar. Mean Flow", "Apr. Mean Flow", "May Mean Flow", "Jun. Mean Flow", "Jul. Mean Flow", "Aug. Mean Flow", "Sep. Mean Flow", "Oct. Mean Flow", "Nov. Mean Flow", "Dec. Mean Flow"); Table1 <- signif(Table1, digits = 3) Table1 <- round(Table1, digits = 2) Table1 <- kable(format(Table1, digits = 3, drop0trailing = TRUE)) Table2 <- matrix(c(percent_difference[1,2], percent_difference[1,3], percent_difference[1,4], percent_difference[1,5], percent_difference[1,6], percent_difference[1,7], percent_difference[1,8], percent_difference[1,9], percent_difference[1,10], percent_difference[1,11], percent_difference[1,12], percent_difference[1,13], percent_difference[2,2], percent_difference[2,3], percent_difference[2,4], percent_difference[2,5], percent_difference[2,6], percent_difference[2,7], percent_difference[2,8], percent_difference[2,9], percent_difference[2,10], percent_difference[2,11], percent_difference[2,12], percent_difference[2,13], percent_difference[3,2], percent_difference[3,3], percent_difference[3,4], percent_difference[3,5], percent_difference[3,6], percent_difference[3,7], percent_difference[3,8], percent_difference[3,9], percent_difference[3,10], percent_difference[3,11], percent_difference[3,12], percent_difference[3,13]), nrow = 12, ncol = 3); colnames(Table2) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference"); rownames(Table2) = c("Jan. Low Flow", "Feb. Low Flow", "Mar. Low Flow", "Apr. Low Flow", "May Low Flow", "Jun. Low Flow", "Jul. Low Flow", "Aug. Low Flow", "Sep. Low Flow", "Oct. Low Flow", "Nov. Low Flow", "Dec. Low Flow"); Table2 <- signif(Table2, digits = 3) Table2 <- round(Table2, digits = 2) Table2 <- kable(format(Table2, digits = 3, drop0trailing = TRUE)) # Table 3: Monthly High Flow Table3 <- matrix(c(percent_difference[1,26], percent_difference[1,27], percent_difference[1,28],percent_difference[1,29], percent_difference[1,30], percent_difference[1,31],percent_difference[1,32], percent_difference[1,33], percent_difference[1,34],percent_difference[1,35], percent_difference[1,36], percent_difference[1,37],percent_difference[2,26],percent_difference[2,27], percent_difference[2,28], percent_difference[2,29],percent_difference[2,30], percent_difference[2,31], percent_difference[2,32],percent_difference[2,33], percent_difference[2,34], percent_difference[2,35],percent_difference[2,36], percent_difference[2,37],percent_difference[3,26], percent_difference[3,27], percent_difference[3,28],percent_difference[3,29], percent_difference[3,30], percent_difference[3,31],percent_difference[3,32], percent_difference[3,33], percent_difference[3,34],percent_difference[3,35], percent_difference[3,36], percent_difference[3,37]), nrow = 12, ncol = 3); colnames(Table3) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference"); rownames(Table3) = c("Jan. High Flow", "Feb. High Flow", "Mar. High Flow", "Apr. High Flow", "May High Flow", "Jun. High Flow", "Jul. High Flow", "Aug. High Flow", "Sep. High Flow", "Oct. High Flow", "Nov. High Flow", "Dec. High Flow"); Table3 <- signif(Table3, digits = 3) Table3 <- round(Table3, digits = 2) Table3 <- kable(format(Table3, digits = 3, drop0trailing = TRUE)) Table4 <- matrix(c(percent_difference[1,38], percent_difference[1,43], percent_difference[1,39], percent_difference[1,44], percent_difference[1,40], percent_difference[1,45], percent_difference[1,41], percent_difference[1,46], percent_difference[1,42], percent_difference[1,47], percent_difference[1,59], percent_difference[1,60], percent_difference[1,58], percent_difference[1,67], percent_difference[2,38], percent_difference[2,43], percent_difference[2,39], percent_difference[2,44], percent_difference[2,40], percent_difference[2,45], percent_difference[2,41], percent_difference[2,46], percent_difference[2,42], percent_difference[2,47], percent_difference[2,59], percent_difference[2,60], percent_difference[2,58], percent_difference[2,67], percent_difference[3,38], percent_difference[3,43], percent_difference[3,39], percent_difference[3,44], percent_difference[3,40], percent_difference[3,45], percent_difference[3,41], percent_difference[3,46], percent_difference[3,42], percent_difference[3,47], percent_difference[3,59], percent_difference[3,60], percent_difference[3,58], percent_difference[3,67]), nrow = 14, ncol = 3); colnames(Table4) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference"); rownames(Table4) = c("Min. 1 Day Min", "Med. 1 Day Min", "Min. 3 Day Min", "Med. 3 Day Min", "Min. 7 Day Min", "Med. 7 Day Min", "Min. 30 Day Min", "Med. 30 Day Min", "Min. 90 Day Min", "Med. 90 Day Min", "7Q10", "Year of 90-Day Min. Flow", "Drought Year Mean", "Mean Baseflow"); Table4 <- signif(Table4, digits = 3) Table4 <- round(Table4, digits = 2) Table4 <- kable(format(Table4, digits = 3, drop0trailing = TRUE)) # Table 5: Period High Flows Table5 <- matrix(c(percent_difference[1,48], percent_difference[1,53], percent_difference[1,49], percent_difference[1,54], percent_difference[1,50], percent_difference[1,55], percent_difference[1,51], percent_difference[1,56], percent_difference[1,52], percent_difference[1,57], percent_difference[2,48], percent_difference[2,53], percent_difference[2,49], percent_difference[2,54], percent_difference[2,50], percent_difference[2,55], percent_difference[2,51], percent_difference[2,56], percent_difference[2,52], percent_difference[2,57], percent_difference[3,48], percent_difference[3,53], percent_difference[3,49], percent_difference[3,54], percent_difference[3,50], percent_difference[3,55], percent_difference[3,51], percent_difference[3,56], percent_difference[3,52], percent_difference[3,57]), nrow = 10, ncol = 3); colnames(Table5) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference"); rownames(Table5) = c("Max. 1 Day Max", "Med. 1 Day Max", "Max. 3 Day Max", "Med. 3 Day Max", "Max. 7 Day Max", "Med. 7 Day Max", "Max. 30 Day Max", "Med. 30 Day Max", "Max. 90 Day Max", "Med. 90 Day Max"); Table5 <- signif(Table5, digits = 3) Table5 <- round(Table5, digits = 2) Table5 <- kable(format(Table5, digits = 3, drop0trailing = TRUE)) # Table 6: Non-Exceedance Flows Table6 <- matrix(c(percent_difference[1,62], percent_difference[1,63], percent_difference[1,64],percent_difference[1,65], percent_difference[1,66], percent_difference[1,61],percent_difference[2,62], percent_difference[2,63], percent_difference[2,64],percent_difference[2,65], percent_difference[2,66], percent_difference[2,61],percent_difference[3,62], percent_difference[3,63], percent_difference[3,64],percent_difference[3,65], percent_difference[3,66], percent_difference[3,61]), nrow = 6, ncol = 3); colnames(Table6) = c("VAHydro Scen. 1", "VAHydro Scen. 2", "Pct. Difference"); rownames(Table6) = c("1% Non-Exceedance", "5% Non-Exceedance", "50% Non-Exceedance", "95% Non-Exceedance", "99% Non-Exceedance", "Sept. 10% Non-Exceedance"); Table6 <- signif(Table6, digits = 3) Table6 <- round(Table6, digits = 2) Table6 <- kable(format(Table6, digits = 3, drop0trailing = TRUE)) # Creating names for plot legends name_scenario1 <- paste('Scenario 1', mod.scenario1); name_scenario2 <- paste('Scenario 2', mod.scenario2); # CREATE FUNCTIONS FOR PLOTTING ------------------------------------------ base_breaks <- function(n = 10){ function(x) { axisTicks(log10(range(x, na.rm = TRUE)), log = TRUE, n = n) } } scaleFUN <- function(x) sprintf("%.0f", x)
include_graphics("gis.png")
description.cont <- paste0(" The average daily discharge change between scenario 1 and scenario 2 for the 20 year timespan was ", OUTPUT_MATRIXsaver[2] ,"%, with ", OUTPUT_MATRIXsaver[1], "% of its rolling three month time spans above 20% difference.") cat(description) cat(description.cont)
\pagebreak
Table 1: Monthly Low Flows
Table2
Table 2: Monthly Average Flows
Table1
\pagebreak
Table 3: Monthly High Flows
Table3
Table 4: Period Low Flows
Table4
Table 5: Period High Flows
Table5
Table 6: Non-Exceedance Flows
Table6
# SETTING UP PLOTS # Basic hydrograph ----- # Max/min for y axis scaling max <- max(c(max(all_data$`Scenario 1 Flow`), max(all_data$`Scenario 2 Flow`))); min <- min(c(min(all_data$`Scenario 1 Flow`), min(all_data$`Scenario 2 Flow`))); if (max > 10000){ max <- 100000 }else if (max > 1000){ max <- 10000 }else if (max > 100){ max <- 1000 }else if (max > 10){ max <- 100 } if (min>100){ min<-100 }else if (min>10){ min<-10 }else min<-1 if (min==100){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(100, 1000, 10000, 100000), limits=c(min,max)) }else if (min==10){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(10, 100, 1000, 10000), limits=c(min,max)) }else if (min==1){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(1, 10, 100, 1000, 10000), limits=c(min,max)) }else fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(), labels=scaleFUN, limits=c(min,max)) df <- data.frame(as.Date(all_data$Date), all_data$`Scenario 1 Flow`, all_data$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=0.5) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=0.5)+ fixtheyscale+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=12, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(y = "Flow (cfs)") ggsave("fig1.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in') # Zoomed hydrograph in year of lowest 90-year flow ----- # Running scenario 1 calculations f3_scenario1 <- zoo(all_data$`Scenario 1 Flow`, order.by = all_data$Date) g2_scenario1 <- group2(f3_scenario1, year = 'water') # Running scenairo 2 calculations f3_scenario2 <- zoo(all_data$`Scenario 2 Flow`, order.by = all_data$Date) g2_scenario2 <- group2(f3_scenario2, year = 'water') yearly_scenario1_90DayMin <- g2_scenario1[,c(1,10)]; yearly_scenario2_90DayMin <- g2_scenario2[,c(1,10)]; low.year <- subset(yearly_scenario1_90DayMin, yearly_scenario1_90DayMin$`90 Day Min`==min(yearly_scenario1_90DayMin$`90 Day Min`)); low.year <- low.year$year; low.year <- subset(all_data, year(all_data$Date)==low.year); # Scaling using max/min max <- max(c(max(low.year$`Scenario 1 Flow`), max(low.year$`Scenario 2 Flow`))); min <- min(c(min(low.year$`Scenario 1 Flow`), min(low.year$`Scenario 2 Flow`))); if (max > 10000){ max <- 100000 }else if (max > 1000){ max <- 10000 }else if (max > 100){ max <- 1000 }else if (max > 10){ max <- 100 } if (min>100){ min<-100 }else if (min>10){ min<-10 }else min<-1 if (min==100){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(100, 1000, 10000, 100000), limits=c(min,max)) }else if (min==10){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(10, 100, 1000, 10000), limits=c(min,max)) }else if (min==1){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(1, 10, 100, 1000, 10000), limits=c(min,max)) }else fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(), labels=scaleFUN, limits=c(min,max)) df <- data.frame(as.Date(low.year$Date), low.year$`Scenario 1 Flow`, low.year$`Scenario 2 Flow`); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=0.7) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=0.7)+ fixtheyscale+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=11, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(y = "Flow (cfs)") ggsave("fig2.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in') #Flow exceedance plot ----- # Determining the "rank" (0-1) of the flow value num_observations <- as.numeric(length(all_data$Date)) rank_vec <- as.numeric(c(1:num_observations)) # Calculating exceedance probability prob_exceedance <- 100*((rank_vec) / (num_observations + 1)) exceed_scenario1 <- sort(all_data$`Scenario 1 Flow`, decreasing = TRUE) exceed_scenario2 <- sort(all_data$`Scenario 2 Flow`, decreasing = TRUE) scenario1_exceedance <- quantile(exceed_scenario1, probs = c(0.01, 0.05, 0.5, 0.95, 0.99)) scenario2_exceedance <- quantile(exceed_scenario2, probs = c(0.01, 0.05, 0.5, 0.95, 0.99)) # Determining max flow value for exceedance plot scale max <- max(c(max(scenario1_exceedance), max(scenario2_exceedance))); min <- min(c(min(scenario1_exceedance), min(scenario2_exceedance))); if (max > 10000){ max <- 100000 }else if (max > 1000){ max <- 10000 }else if (max > 100){ max <- 1000 }else if (max > 10){ max <- 100 } if (min>100){ min<-100 }else if (min>10){ min<-10 }else min<-1 if (min==100){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(100, 1000, 10000, 100000), limits=c(min,max)) }else if (min==10){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(10, 100, 1000, 10000), limits=c(min,max)) }else if (min==1){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(1, 10, 100, 1000, 10000), limits=c(min,max)) }else fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(), labels=scaleFUN, limits=c(min,max)) df <- data.frame(prob_exceedance, exceed_scenario1, exceed_scenario2); colnames(df) <- c('Date', 'Scenario1', 'Scenario2') options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1, color="VAHydro Scen. 1"), size=0.8) + geom_line(aes(y=Scenario2, color="VAHydro Scen. 2"), size=0.8)+ fixtheyscale+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=12, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(x= "Probability of Exceedance (%)", y = "Flow (cfs)") ggsave("fig3.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in') # Baseflow Indiviudal Graph ----- data1$year <- year(ymd(data1$date)) data1$month <- month(ymd(data1$date)) data1$day <- day(ymd(data1$date)) data2$year <- year(ymd(data2$date)) data2$month <- month(ymd(data2$date)) data2$day <- day(ymd(data2$date)) scenario1river <- createlfobj(data1, hyearstart = 10, baseflow = TRUE, meta = NULL) scenario2river <- createlfobj(data2, hyearstart = 10, baseflow = TRUE, meta = NULL) baseflowriver<- data.frame(scenario1river, scenario2river); colnames(baseflowriver) <-c ('mday', 'mmonth', 'myear', 'mflow', 'mHyear', 'mBaseflow', 'gday', 'gmonth', 'gyear', 'gflow', 'gHyear', 'gBaseflow') # removing NA values baseflowriver<-baseflowriver[complete.cases(baseflowriver)==TRUE,] scenario2river<- data.frame(baseflowriver$gday, baseflowriver$gmonth, baseflowriver$gyear, baseflowriver$gflow, baseflowriver$gHyear, baseflowriver$gBaseflow); scenario1river<- data.frame(baseflowriver$mday, baseflowriver$mmonth, baseflowriver$myear, baseflowriver$mflow, baseflowriver$mHyear, baseflowriver$mBaseflow) names(scenario1river) <- c('day', 'month', 'year', 'flow', 'hyear', 'baseflow') names(scenario2river) <- c('day', 'month', 'year', 'flow', 'hyear', 'baseflow') # Adding date vectors scenario1river$date <- as.Date(paste0(scenario1river$year,"-",scenario1river$month,"-",scenario1river$day)) scenario2river$date <- as.Date(paste0(scenario2river$year,"-",scenario2river$month,"-",scenario2river$day)) # Determining max flow value for plot scale max <- max(c(max(scenario1river$baseflow), max(scenario2river$baseflow), max(scenario1river$flow), max(scenario2river$flow))); min <- min(c(min(scenario1river$baseflow), min(scenario2river$baseflow), min(scenario1river$flow), min(scenario2river$flow))); if (max > 10000){ max <- 100000 }else if (max > 1000){ max <- 10000 }else if (max > 100){ max <- 1000 }else if (max > 10){ max <- 100 } if (min>100){ min<-100 }else if (min>10){ min<-10 }else min<-1 if (min==100){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(100, 1000, 10000, 100000), limits=c(min,max)) }else if (min==10){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(10, 100, 1000, 10000), limits=c(min,max)) }else if (min==1){ fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = c(1, 10, 100, 1000, 10000), limits=c(min,max)) }else fixtheyscale<- scale_y_continuous(trans = log_trans(), breaks = base_breaks(), labels=scaleFUN, limits=c(min,max)) df <- data.frame(as.Date(scenario1river$date), scenario1river$baseflow, scenario2river$baseflow, scenario1river$flow, scenario2river$flow); colnames(df) <- c('Date', 'Scenario1Baseflow', 'Scenario2Baseflow','Scenario1Flow', 'Scenario2Flow') options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1Baseflow, color="VAHydro Scen. 1"), size=0.5) + geom_line(aes(y=Scenario2Baseflow, color="VAHydro Scen. 2"), size=0.5)+ fixtheyscale+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=12, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(y = "Flow (cfs)") ggsave("fig4.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in') par(mfrow = c(1,1)); options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_line(aes(y=Scenario1Flow, color=paste("Scen. 1 Flow")), size=1)+ geom_line(aes(y=Scenario2Flow, color=paste("Scen. 2 Flow")), size=0.5)+ geom_line(aes(y=Scenario1Baseflow, color=paste("Scen. 1 Baseflow")), size=0.5) + geom_line(aes(y=Scenario2Baseflow, color=paste("Scen. 2 Baseflow")), size=1)+ fixtheyscale+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=12, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("black","red","grey", "light pink"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(y = "Flow (cfs)") ggsave("fig5.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in') # Setup for Residuals data <- all_data[complete.cases(all_data),] resid <- (data$`Scenario 2 Flow` - data$`Scenario 1 Flow`) resid <- data.frame(data$Date, resid) # Residuals plot for hydrograph zeroline <- rep_len(0, length(data$Date)) quantresid <- data.frame(signif(quantile(resid$resid, na.rm = TRUE), digits=3)) min <- min(resid$resid) max <- max(resid$resid) names(quantresid) <- c('Percentiles') df <- data.frame(as.Date(resid$data.Date), resid$resid, zeroline); colnames(df) <- c('Date', 'Residual', 'Zeroline') options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_point(aes(y=Residual, color="Residual (CC - Base)"), size=1) + geom_line(aes(y=Zeroline, color="Zeroline"), size=0.8)+ scale_y_continuous(limits=c(min,max))+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=12, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("dark green","black"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(y = "Flow Difference (cfs)") ggsave("fig9A.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in') # Setup for Residuals data <- all_data[complete.cases(all_data),] #data_weighted <- data/area resid <- 1000000*((data$`Scenario 2 Flow` - data$`Scenario 1 Flow`)/area) resid <- data.frame(data$Date, resid) # Residuals plot for hydrograph zeroline <- rep_len(0, length(data$Date)) quantresid <- data.frame(signif(quantile(resid$resid, na.rm =TRUE), digits=3)) min <- min(resid$resid) max <- max(resid$resid) names(quantresid) <- c('Percentiles') df <- data.frame(as.Date(resid$data.Date), resid$resid, zeroline); colnames(df) <- c('Date', 'Residual', 'Zeroline') options(scipen=5, width = 1400, height = 950) myplot <- ggplot(df, aes(x=Date)) + geom_point(aes(y=Residual, color="Residual (CC - Base)"), size=1) + geom_line(aes(y=Zeroline, color="Zeroline"), size=0.8)+ scale_y_continuous(limits=c(min,max))+ theme_bw()+ theme(legend.position="top", legend.title=element_blank(), legend.box = "horizontal", legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="white"), legend.text=element_text(size=12), axis.text=element_text(size=12, colour="black"), axis.title=element_text(size=14, colour="black"), axis.line = element_line(colour = "black", size = 0.5, linetype = "solid"), axis.ticks = element_line(colour="black"), panel.grid.major=element_line(colour = "light grey"), panel.grid.minor=element_blank())+ scale_colour_manual(values=c("dark green","black"))+ guides(colour = guide_legend(override.aes = list(size=5)))+ labs(y = "Area Weighted Flow Difference*10^6 (ft/s)") ggsave("fig9B.png", plot = myplot, device = 'png', width = 8, height = 5.5, units = 'in')
Fig. 1: Hydrograph
include_graphics("fig1.png")
Fig. 2: Zoomed Hydrograph
include_graphics("fig2.png")
Fig. 3: Flow Exceedance
include_graphics("fig3.png")
Fig. 4: Baseflow
include_graphics("fig4.png")
Fig. 5: Combined Baseflow
include_graphics("fig5.png")
Fig. 6: Largest Difference Segment
include_graphics("fig6.png")
Fig. 7: Second Largest Difference Segment
include_graphics("fig7.png")
Fig. 8: Third Largest Difference Segment
include_graphics("fig8.png")
Fig. 9A: Residuals Plot
include_graphics("fig9A.png")
Fig. 9B: Area Weighted Residuals Plot
include_graphics("fig9B.png")
dat <- vahydro_import_local.runoff.inflows_cfs(riv.seg, run.id1, token, site, start.date, end.date); dat <- subset(dat, dat$date >= start.date & dat$date <= end.date); boxplot(as.numeric(dat$flow.unit) ~ year(dat$date), outline = FALSE, ylab = 'Runit Flow (cfs)', xlab = 'Date') dev.copy(png, 'fig10.png') dev.off()
Fig. 10: VA Hydro Scen. 1 Runit Values (Outliers Excluded)
include_graphics("fig10.png")
Fig. 11: Smallest Difference Segment
include_graphics("fig11.png")
Fig. 12: Second Smallest Difference Segment
include_graphics("fig12.png")
Fig. 13: Third Smallest Difference Segment
include_graphics("fig13.png")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.