riv.seg <- params$riv.seg
dat.source1 <- params$dat.source1
dat.source2 <- params$dat.source2
start.date <- params$start.date
end.date <- params$end.date
github_location <- params$github_location
site <- params$site
site.or.server <- params$site.or.server
run.id1 <- params$run.id1
run.id2 <- params$run.id2
gage_number <- params$gage_number
mod.phase1 <- params$mod.phase1
mod.scenario1 <- params$mod.scenario1
mod.phase2 <- params$mod.phase2
mod.scenario2 <- params$mod.scenario2
cn1 <- params$cn1
cn2 <- params$cn2
gage.timespan.trimmed <- params$gage.timespan.trimmed
export_path_custom <- params$export_path_custom
basepath='/var/www/R';
source(paste(basepath,'config.R',sep='/'))
library("knitr")
library("hydrotools")
ds <- RomDataSource$new("http://deq1.bse.vt.edu/d.dh", rest_uname)
ds$get_token(rest_pw)

cbp6_link <- paste0(github_location, "/cbp6/code");
if (!is.logical(export_path_custom)) {
  # allow override from command line, otherwise use value from config files
  export_path <- export_path_custom
}
if (gage_number == FALSE | is.na(gage_number) == TRUE) {
  gage_number <- '00000000'
}
mrun_name1 <- paste0('runid_', run.id1)
mrun_name2 <- paste0('runid_', run.id2)
# LOADING DATA ------------------------------------------------------------
# Load feature
hydrocode = paste0("vahydrosw_wshed_", riv.seg);
message(paste("Calling om_get_feature(",site, ", ", hydrocode, ', watershed', ', vahydro'))
feature <- RomFeature$new(ds,list(hydrocode=hydrocode, bundle='watershed', ftype='vahydro'),TRUE)
hydroid = feature$hydroid
if (dat.source1 == 'gage') {
  id1 <- gage_number
} else {
  id1 <- riv.seg
  if (gage.timespan.trimmed == TRUE) {
    mrun_name1 <- paste0(mrun_name1,'_gage_timespan')
  }
}
data1 <- om_get_watershed_model_data(dat.source1, id1, run.id1, 
  mod.phase1, mod.scenario1, start.date, end.date
)
message("Getting 1st Model")
m1 <- RomProperty$new(
  ds, list(featureid = hydroid, entity_type = 'dh_feature', propcode = mod.scenario1)
)
message("Getting 1st Scenario element")
model_run1 <- RomProperty$new(
  ds, list(featureid = m1$pid, entity_type = 'dh_properties', propname = mrun_name1)
)
scenprop1 <- model_run1$pid

# Now get data set to compare against
if (dat.source2 == 'gage') {
  id2 <- gage_number
} else {
  id2 <- riv.seg
  if (gage.timespan.trimmed == TRUE) {
    mrun_name2 <- paste0(mrun_name2,'_gage_weighted')
  }
}

message("Getting 2nd Model")
m2 <- RomProperty$new(
  ds, 
  list(
    featureid = hydroid, entity_type = 'dh_feature', propcode = mod.scenario2
  )
)
message("Getting 2nd Scenario element")
model_run2 <- RomProperty$new(
  ds, list(featureid = m2$pid, entity_type = 'dh_properties', propname = mrun_name1)
)
scenprop2 <- model_run2$pid

if (dat.source2 == 'vahydro') {
  elid2 = om_get_model_elementid(site, m2$pid) 
  data2 <- om_get_rundata(elid = elid2, runid = run.id2, site = omsite, FALSE, TRUE)
  data2$thisdate <- as.character(as.Date(index(data2)))
  dat.date <- as.Date(as.character(data2$thisdate))
  dat.flow <- as.numeric(data2$Qout)
  data2 <- data.frame(dat.date, dat.flow, row.names = NULL)
  colnames(data2) <- c('date','flow')
} else {
  # we still use the old plumbing for this
  data2 <- om_get_watershed_model_data(
    dat.source2, id2, run.id2, 
    mod.phase2, mod.scenario2, 
    start.date, end.date, 
    token = token, site = omsite
  )
}
# TRIMMING DATA TO PROPER WATER YEAR
data1 <- water_year_trim(data1)
data2 <- water_year_trim(data2)
if (dat.source1 == 'vahydro') {
  title.part.1 <- paste('VA Hydro Run ', run.id1, sep = '')
  name.part1 <- paste('runid_', run.id1, sep = '')
} else if (dat.source1 == 'gage') {
  title.part.1 <- paste('USGS Gage ', gage_number, sep = '')
  name.part1 <- paste('gage_', gage_number, sep = '')
} else if (dat.source1 == 'cbp_model') {
  title.part.1 <- paste('Scenario ', mod.scenario1, sep = '')
  name.part1 <- paste('scen_', mod.scenario1, sep = '')
}
if (dat.source2 == 'vahydro') {
  title.part.2 <- paste('VA Hydro Run ', run.id2, sep = '')
  name.part2 <- paste('runid_', run.id2, sep = '')
} else if (dat.source2 == 'gage') {
  title.part.2 <- paste('USGS Gage ', gage_number, sep = '')
  name.part2 <- paste('gage_', gage_number, sep = '')
} else if (dat.source2 == 'cbp_model') {
  title.part.2 <- paste('Scenario ', mod.scenario2, sep = '')
  name.part2 <- paste('scen_', mod.scenario2, sep = '')
}
dashboard_title <- paste('River Segment ', riv.seg, ': ', title.part.1, ' vs. ', title.part.2, sep = '')

title: "r dashboard_title" header-includes: - \usepackage{titling} - \pretitle{\begin{flushleft}} - \posttitle{\end{flushleft}} output: pdf_document


# POSTING METRICS TO SCENARIO PROPERTIES ON VA HYDRO
metrics1 <- vahydro_import_all_metrics_from_scenprop(scenprop1, site, token)
metrics2 <- vahydro_import_all_metrics_from_scenprop(scenprop2, site, token)
percent_difference <- metrics_compare(metrics1, metrics2, riv.seg)
all_data <- all_data_maker(data1, data2)
# linking land-river segments
cbp6.landunits <- link.cbp6.lrseg.hydrocodes(riv.seg, psk, site, token)
if (dat.source1 != 'vahydro') {
  run.id <- run.id2
} else {
  run.id <- run.id1
}
# Downloading local runoff inflow data
lri.dat <- vahydro_import_local.runoff.inflows_cfs(riv.seg, run.id, token, site, start.date, end.date);
if (!is.logical(lri.dat)) {
  lri.dat <- subset(lri.dat, lri.dat$date >= start.date & lri.dat$date <= end.date);
}  
# Creating all tables
Table1 <- tab1.monthly.low.flows(percent_difference, cn1, cn2)
Table2 <- tab2.monthly.average.flows(percent_difference, cn1, cn2)
Table3 <- tab3.monthly.high.flows(percent_difference, cn1, cn2)
Table4 <- tab4.period.low.flows(percent_difference, cn1, cn2)
Table5 <- tab5.period.high.flows(percent_difference, cn1, cn2)
Table6 <- tab6.nonexceedence.flows(percent_difference, cn1, cn2)
if (is.data.frame(cbp6.landunits)) {

  for (i in 1:length(cbp6.landunits)) {
    tmp.data <- vahydro_import_lrseg_all_flows(cbp6.landunits[i], run.id, token, site, start.date, end.date)
    tmp.data <- tmp.data[complete.cases(tmp.data),]

    namer <- paste0('tab.', cbp6.landunits[i], '.means.by.flow')
    tmp.tab <- tab.means.by.lrseg.flow(tmp.data)
    assign(namer, tmp.tab)

    namer <- paste0('tab.', cbp6.landunits[i], '.zero.day.ratios.by.flow')
    tmp.tab <- tab.zero.day.ratios.by.lrseg.flow(tmp.data)
    assign(namer, tmp.tab)

    namer <- paste0('tab.SURO.', cbp6.landunits[i], '.iqr.by.lrseg.flow.annual')
    tmp.tab <- tab.iqr.by.lrseg.flow.annual(tmp.data, flow.abbreviation = 'suro')
    assign(namer, tmp.tab)

    namer <- paste0('tab.IFWO.', cbp6.landunits[i], '.iqr.by.lrseg.flow.annual')
    tmp.tab <- tab.iqr.by.lrseg.flow.annual(tmp.data, flow.abbreviation = 'ifwo')
    assign(namer, tmp.tab)

    namer <- paste0('tab.AGWO.', cbp6.landunits[i], '.iqr.by.lrseg.flow.annual')
    tmp.tab <- tab.iqr.by.lrseg.flow.annual(tmp.data, flow.abbreviation = 'agwo')
    assign(namer, tmp.tab)

    namer <- paste0('tab.', cbp6.landunits[i], '.means.by.land.use')
    tmp.tab <- tab.means.by.lrseg.land.use(tmp.data)
    assign(namer, tmp.tab)

    namer <- paste0('tab.', cbp6.landunits[i], '.zero.day.ratios.by.land.use')
    tmp.tab <- tab.zero.day.ratios.by.lrseg.land.use(tmp.data)
    assign(namer, tmp.tab)
  }
}
wd <- getwd()
if (!is.logical(lri.dat)) {
  tab.iqr.by.lrseg.lri.annual <- tab.iqr.by.lrseg.lri.annual(lri.dat)
}
# Creating all plots
fig1.hydrograph(all_data, cn1, cn2, export_path)
fig2.zoomed.hydrograph(all_data, cn1, cn2, export_path)
fig3.flow.exceedance(all_data, cn1, cn2, export_path)
fig4.baseflow.hydrograph(all_data, cn1, cn2, export_path)
fig5.combined.hydrograph(all_data, export_path)
OUTPUT_MATRIXsaver <- figs6to8.largest.diff.periods(all_data, cn1, cn2, export_path)
fig9a.residual.plot(all_data, cn1, cn2, export_path)
# commented out due to retrieving data inside the function which breaks.
#fig9b.area.weighted.residual.plot(all_data, riv.seg, token, site, cn1, cn2, export_path)
if (!is.logical(lri.dat)) {
  fig10.runit.boxplot(lri.dat, export_path)
}
figs11to13.smallest.diff.periods(all_data, cn1, cn2, export_path)
fig.gis(riv.seg, site_number = gage_number, site_url = site, cbp6_link, token, export_path)
description <- try(read_file(paste0(cbp6_link, "/gage_descriptions/", gage_number, ".txt")))
if (class(description) == "try-error") {
  description <- ""
}
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.")
if (dat.source1 == 'gage' | dat.source2 == 'gage') {
  if (dat.source1 == 'gage') {
    gage.dat <- data1
    scen.dat <- data2
    nse.val <- nse(gage.dat, scen.dat)
  } else if (dat.source2 == 'gage') {
    scen.dat <- data1
    gage.dat <- data2
    nse.val <- nse(gage.dat, scen.dat)
  }
  description.cont2 <- paste0(' The Nash-Sutcliffe Efficiency of the model, calculated between the gage and scenario data, was found to be ', nse.val, '.  Gage data was available for ', signif(length(data1$flow)*100/length(data2$flow), 4), '% of the model timespan.')
} else {
  description.cont2 <- ''
}
include_graphics(paste0(export_path,"gis.png"))
# Loading written gage description
cat(description)
cat(description.cont)
cat(description.cont2)

\pagebreak

Table 1: Monthly Low Flows

Table1

Table 2: Monthly Average Flows

Table2

Table 3: Monthly High Flows

Table3

Table 4: Period Low Flows

Table4

Table 5: Period High Flows

Table5

Table 6: Non-Exceedance Flows

Table6

\pagebreak

Fig. 1: Hydrograph

include_graphics(paste0(export_path,"fig1.png"))

Fig. 2: Zoomed Hydrograph

include_graphics(paste0(export_path,"fig2.png"))

Fig. 3: Flow Exceedance

include_graphics(paste0(export_path,"fig3.png"))

Fig. 4: Baseflow

include_graphics(paste0(export_path,"fig4.png"))

Fig. 5: Combined Baseflow

include_graphics(paste0(export_path,"fig5.png"))

Fig. 6: Largest Difference Period

include_graphics(paste0(export_path,"fig6.png"))

Fig. 7: Second Largest Difference Period

include_graphics(paste0(export_path,"fig7.png"))

Fig. 8: Third Largest Difference Period

include_graphics(paste0(export_path,"fig8.png"))

Fig. 9A: Residuals Plot

include_graphics(paste0(export_path,"fig9A.png"))

Fig. 9B: Area Weighted Residuals Plot

#include_graphics(paste0(export_path,"fig9B.png"))
cat("Area-weighed residual plot disabled.")

Fig. 10: VA Hydro Scen. 1 Runit Values (Outliers Excluded)

cat('\n')
if (!is.logical(lri.dat)) {
include_graphics(paste0(export_path,"fig10.png"))
  tab.iqr.by.lrseg.lri.annual
  cat(paste0('## Tab: Annual IQR of Local Runoff Inflows'))
} else {
  cat(paste0('## Tab: Annual IQR of Local Runoff Inflows'))
  cat("Land Use Runoff data not available")
}

Fig. 11: Smallest Difference Period

include_graphics(paste0(export_path,"fig11.png"))

Fig. 12: Second Smallest Difference Period

include_graphics(paste0(export_path,"fig12.png"))

Fig. 13: Third Smallest Difference Period

include_graphics(paste0(export_path,"fig13.png"))
cat("\n\n\\pagebreak\n")
if (is.data.frame(cbp6.landunits)) {
  for (i in 1:length(cbp6.landunits)) {
    namer <- paste0('data_', cbp6.landunits[i])
    tmp.data <- vahydro_import_lrseg_all_flows(cbp6.landunits[i], run.id, token, site, start.date, end.date)
    assign(namer, tmp.data)
  }
}

Additional Tables: Land-River Segment Flow Metrics

if (is.data.frame(cbp6.landunits)) {
  for (i in 1:length(cbp6.landunits)) {
    namer <- paste0('tab.', cbp6.landunits[i], '.means.by.flow')
    tmp.tab <- get(namer)
    cat('\n')
    cat(paste0('## Tab: Mean Flows by Flow Type: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat('\n')

    namer <- paste0('tab.', cbp6.landunits[i], '.zero.day.ratios.by.flow')
    tmp.tab <- get(namer)
    cat(paste0('## Tab: Ratio of Zero-Flow Days by Flow Type: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat("\n")

    namer <- paste0('tab.SURO.', cbp6.landunits[i], '.iqr.by.lrseg.flow.annual')
    tmp.tab <- get(namer)
    cat(paste0('## Tab: IQR for SURface Outflow: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat("\n\n\\pagebreak\n")

    namer <- paste0('tab.IFWO.', cbp6.landunits[i], '.iqr.by.lrseg.flow.annual')
    tmp.tab <- get(namer)
    cat('\n')
    cat(paste0('## Tab: IQR for InterFloW Outflow: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat("\n")

    namer <- paste0('tab.AGWO.', cbp6.landunits[i], '.iqr.by.lrseg.flow.annual')
    tmp.tab <- get(namer)
    cat(paste0('## Tab: IQR for Active GroundWater Outflow: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat("\n\n\\pagebreak\n")

    namer <- paste0('tab.', cbp6.landunits[i], '.means.by.land.use')
    tmp.tab <- get(namer)
    cat('\n')
    cat(paste0('## Tab: Mean Flows by Land Use: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat("\n\n\\pagebreak\n")

    namer <- paste0('tab.', cbp6.landunits[i], '.zero.day.ratios.by.land.use')
    tmp.tab <- get(namer)
    cat('\n')
    cat(paste0('## Tab: Ratio of Zero-Flow Days by Land Use: LR-Seg ', cbp6.landunits[i]))
    print(tmp.tab)
    cat("\n\n\\pagebreak\n")
  }
}

Additional Figures: Land-River Segment Flow Boxplots

if (is.data.frame(cbp6.landunits)) {
  for (i in 1:length(cbp6.landunits)) {
    namer <- paste0('data_', cbp6.landunits[i])
    tmp.data <- get(namer)

    cat('\n')
    cat(paste0('## Fig: Annual SURO Flows for LR-seg ', cbp6.landunits[i]))
    plot <- fig.boxplot.by.flow.for.dash(tmp.data, flow.abbreviation = 'suro')
    plot
    cat("\n\n\\pagebreak\n") 

    cat('\n')
    cat(paste0('## Fig: Annual IFWO Flows for LR-seg ', cbp6.landunits[i]))
    plot <- fig.boxplot.by.flow.for.dash(tmp.data, flow.abbreviation = 'ifwo')
    plot
    cat("\n\n\\pagebreak\n") 

    cat('\n')
    cat(paste0('## Fig: Annual AGWO Flows for LR-seg ', cbp6.landunits[i]))
    plot <- fig.boxplot.by.flow.for.dash(tmp.data, flow.abbreviation = 'agwo')
    plot
    cat("\n\n\\pagebreak\n") 
  }
}
vahydro.propcode <- get.overall.vahydro.prop(riv.seg, site, token)
dashboard.name = paste0("dashboard_file: ", name.part1, '_vs_', name.part2)
dashboard.url = paste0('http://deq2.bse.vt.edu/p6/p6_gb604/out/dashboards/run_', run.id1, '_vs_run_', run.id2, '/', riv.seg, '.pdf')
vahydro_post_metric_to_scenprop(scenprop.pid = vahydro.propcode, met.varkey = 'external_file',
                                met.propcode = dashboard.url, met.name = dashboard.name, 
                                met.value = 0, site = site, token = token)


HARPgroup/cbp6 documentation built on Oct. 14, 2024, 4:19 p.m.