#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 # Set defaults for river segment and run ID riv.seg <- 'YP3_6330_6700' run.id <- '11' # load command line args argst <- commandArgs(trailingOnly=T) # Use command line args for riverseg and runid if supplied, i.e., this script is run as: # Rscript Working_Gage_Vs_Dashboard_2019.Rmd RU5_6030_0001 121 if (length(argst) >= 1) { riv.seg=argst[1] } if (length(argst) >= 2) { run.id=argst[2] } 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) start.date <- '1984-01-01' end.date <- '2014-12-31' site_number <- '01671020' 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 <- "http://deq2.bse.vt.edu/d.dh" # 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 <- gage_import_data_cfs(site_number, 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.id, token, site = site_url, start.date = start.date, end.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) } #all_data puts scenario flows and corresponding dates in one data frame all_data <- merge(data1, data2, by.x = 'date', by.y = 'date') all_data$counter <- 1:length(all_data$date) # counter fixes issues with row numbers later on in script colnames(all_data) <- c('Date', 'Scenario 1 Flow', 'Scenario 2 Flow', 'Counter') data1 <- data1[data1$date %in% all_data$Date,] data2 <- data2[data2$date %in% all_data$Date,] data1 <- water_year_trim(data1) data2 <- water_year_trim(data2)
title: "River Segment: r paste(riv.seg,'- Scenario : ',mod.scenario1,": Gage ", site_number, " vs. VAHydro")
"
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)
# 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="USGS Gage"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro"), 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)") # 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="USGS Gage"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro"), 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)") # 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="USGS Gage"), size=1) + geom_line(aes(y=Scenario2, color="VAHydro"), 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)")
# 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("USGS Gage", "VAHydro", "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 <- 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("USGS Gage", "VAHydro", "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 <- 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("USGS Gage", "VAHydro", "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 <- 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("USGS Gage", "VAHydro", "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 <- 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("USGS Gage", "VAHydro", "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 <- 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("USGS Gage", "VAHydro", "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 <- 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)
gis_img
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
Table2
Table1
\pagebreak
Table3
Table4
Table5
Table6
# 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="USGS Gage"), size=0.5) + geom_line(aes(y=Scenario2, color="VAHydro"), 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)") myplot
# 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="USGS Gage"), size=0.7) + geom_line(aes(y=Scenario2, color="VAHydro"), 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)") myplot
#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="USGS Gage"), size=0.8) + geom_line(aes(y=Scenario2, color="VAHydro"), 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)") myplot
# 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="USGS Gage"), size=0.5) + geom_line(aes(y=Scenario2Baseflow, color="VAHydro"), 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)") myplot
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("USGS Gage Flow")), size=1)+ geom_line(aes(y=Scenario2Flow, color=paste("VAHydro Flow")), size=0.5)+ geom_line(aes(y=Scenario1Baseflow, color=paste("USGS Gage Baseflow")), size=0.5) + geom_line(aes(y=Scenario2Baseflow, color=paste("VAHydro 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)") myplot
difference1plot
difference2plot
difference3plot
# 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)") myplot
``` {r, echo=FALSE, warning=FALSE, message=FALSE}
# Setup for Residuals data <- all_data[complete.cases(all_data),]
resid <- 1000000*((data$Scenario 2 Flow
- data$Scenario 1 Flow
)/area)
resid <- data.frame(data$Date, resid)
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)") myplot
```
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.