# install and load packages, you only need to do this once. 

rprofile_path = file.path(Sys.getenv("HOME"), ".Rprofile")
write('\noptions(repos=c(getOption(\'repos\'),
          CRAN=\'https://cloud.r-project.org\',
          USGS=\'https://owi.usgs.gov/R\'))\n',
      rprofile_path, 
      append =  TRUE)
install.packages('smwrBase')
install.packages('rloadest')
install.packages('car')
install.packages('unitted')
install.packages('dplyr')
packageurl <- "https://github.com/USGS-R/loadflex/archive/v1.0.1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")
install.packages('dplyr')

options(warn=-1)

library('loadflex')
library('rloadest')

QS_chem <- readRDS("data/QS_chem_UNH.rds")
QSDischargeShort <- readRDS("data/QS_Discharge_USGS.rds")
QSWatershedArea <- 261.5888


QS_chemNa <- QS_chem[complete.cases(QS_chem["Na_mg_L"]),]
QS_chemNa <- QS_chemNa[complete.cases(QS_chemNa["CFS"]),]


meta <- metadata(constituent="Na_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1")


lr <- loadReg2(loadReg(Na_mg_L ~ model(1), data=QS_chemNa,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))


lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemNa, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

# preds_lrNa <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

preds_lcNa <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcNa <- aggregateSolute(preds_lcNa,meta, format="flux total", se.preds=preds_lcNa$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemNa, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lcNaConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcNaConc <- aggregateSolute(preds_lcNaConc,meta, format="conc", se.preds=preds_lcNaConc$se.pred, agg.by="calendar year")

aggs_lcNa$Na_Kg_Ha_yr <- aggs_lcNa$Flux_Total / QSWatershedArea
aggs_lcNa$Na_mg_l <- aggs_lcNaConc$Conc
FluxTotalsdf <- aggs_lcNa[, c('Na_mg_l','Na_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdf$Calendar_Year
FluxTotalsdf$Calendar_Year <- NULL
FluxTotalsdf <- round(FluxTotalsdf, 2)
FluxTotalsdf$Calendar_Year <- years


QS_chemCa <- QS_chem[complete.cases(QS_chem["Ca_mg_L"]),]
QS_chemCa <- QS_chemCa[complete.cases(QS_chemCa["CFS"]),]

meta <- metadata(constituent="Ca_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="Rio Piedras, PR")


lr <- loadReg2(loadReg(Ca_mg_L ~ model(1), data=QS_chemCa,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemCa, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lrCa <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrCa <- aggregateSolute(preds_lrCa,meta, format="flux total", se.preds=preds_lrCa$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemCa, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lcCaConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcCaConc <- aggregateSolute(preds_lcCaConc,meta, format="conc", se.preds=preds_lcCaConc$se.pred, agg.by="calendar year")

aggs_lrCa$Ca_Kg_Ha_yr <- aggs_lrCa$Flux_Total / QSWatershedArea
aggs_lrCa$Ca_mg_l <- aggs_lcCaConc$Conc
FluxTotalsdfCa <- aggs_lrCa[, c('Ca_mg_l','Ca_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfCa$Calendar_Year
FluxTotalsdfCa$Calendar_Year <- NULL
FluxTotalsdfCa <- round(FluxTotalsdfCa, 2)
FluxTotalsdfCa$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfCa,FluxTotalsdf, by = "Calendar_Year")


# FluxTtotalsdf <- merge(FluxTtotalsdfCa,FluxTtotalsdf, by.x='Calendar_Year', by.y='Calendar_Year')
FluxTotals


#Mg
QS_chemMg <- QS_chem[complete.cases(QS_chem["Mg_mg_L"]),]
QS_chemMg <- QS_chemMg[complete.cases(QS_chemMg["CFS"]),]

meta <- metadata(constituent="Mg_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(Mg_mg_L ~ model(1), data=QS_chemMg,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemMg, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lrMg <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrMg <- aggregateSolute(preds_lrMg,meta, format="flux total", se.preds=preds_lrMg$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemMg, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lcMgConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcMgConc <- aggregateSolute(preds_lcMgConc,meta, format="conc", se.preds=preds_lcMgConc$se.pred, agg.by="calendar year")

aggs_lrMg$Mg_Kg_Ha_yr <- aggs_lrMg$Flux_Total / QSWatershedArea
aggs_lrMg$Mg_mg_l <- aggs_lcMgConc$Conc
FluxTotalsdfMg <- aggs_lrMg[, c('Mg_mg_l','Mg_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfMg$Calendar_Year
FluxTotalsdfMg$Calendar_Year <- NULL
FluxTotalsdfMg <- round(FluxTotalsdfMg, 2)
FluxTotalsdfMg$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfMg,FluxTotals, by = "Calendar_Year")
FluxTotals

QS_chemNO3 <- QS_chem[complete.cases(QS_chem["NO3.N_ug_N_L"]),]
QS_chemNO3 <- QS_chemNO3[complete.cases(QS_chemNO3["CFS"]),]

meta <- metadata(constituent="NO3.N_ug_N_L", flow="CFS", 
                 dates="date", conc.units="ug/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(NO3.N_ug_N_L ~ model(1), data=QS_chemNO3,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="ug/L", load.units="kg"))
lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemNO3, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcNO3Conc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcNO3Conc <- aggregateSolute(preds_lcNO3Conc,meta, format="conc", se.preds=preds_lcNO3Conc$se.pred, agg.by="calendar year")

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemNO3, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrNO3 <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrNO3 <- aggregateSolute(preds_lrNO3,meta, format="flux total", se.preds=preds_lrNO3$se.pred, agg.by="calendar year")
aggs_lrNO3$NO3_Kg_Ha_yr <- aggs_lrNO3$Flux_Total / QSWatershedArea
aggs_lrNO3$NO3_mg_l <- aggs_lcNO3Conc$Conc / 1000.0
FluxTotalsdfNO3 <- aggs_lrNO3[, c('NO3_mg_l','NO3_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfNO3$Calendar_Year
FluxTotalsdfNO3$Calendar_Year <- NULL
FluxTotalsdfNO3 <- round(FluxTotalsdfNO3, 3)
FluxTotalsdfNO3$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfNO3,FluxTotals, by = "Calendar_Year")
FluxTotals

#SO4
QS_chemSO4 <- QS_chem[complete.cases(QS_chem["SO4_S_mg_L"]),]
QS_chemSO4 <- QS_chemSO4[complete.cases(QS_chemSO4["CFS"]),]

meta <- metadata(constituent="SO4_S_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(SO4_S_mg_L ~ model(1), data=QS_chemSO4,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))


lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemSO4, abs.or.rel.resids="absolute", interp.function=triangularInterpolation)

preds_lcSO4Conc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcSO4Conc <- aggregateSolute(preds_lcSO4Conc,meta, format="conc", se.preds=preds_lcSO4Conc$se.pred, agg.by="calendar year")

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemSO4, abs.or.rel.resids="relative", interp.function=triangularInterpolation)

preds_lrSO4 <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrSO4 <- aggregateSolute(preds_lrSO4,meta, format="flux total", se.preds=preds_lrSO4$se.pred, agg.by="calendar year")


aggs_lrSO4$SO4_Kg_Ha_yr <- aggs_lrSO4$Flux_Total / QSWatershedArea
aggs_lrSO4$SO4_mg_l <- aggs_lcSO4Conc$Conc
FluxTotalsdfSO4 <- aggs_lrSO4[, c('SO4_mg_l','SO4_Kg_Ha_yr','Calendar_Year')]


#round values 
years <- FluxTotalsdfSO4$Calendar_Year
FluxTotalsdfSO4$Calendar_Year <- NULL
FluxTotalsdfSO4 <- round(FluxTotalsdfSO4, 3)
FluxTotalsdfSO4$Calendar_Year <- years


FluxTotals <- merge(FluxTotalsdfSO4,FluxTotals, by = "Calendar_Year")
FluxTotals

#F
QS_chemF <- QS_chem[complete.cases(QS_chem["Fluoride_mg_L"]),]
QS_chemF <- QS_chemF[complete.cases(QS_chemF["CFS"]),]
QSDischargeShort2<- QSDischargeShort[QSDischargeShort$date >= as.POSIXct("2013-01-01") 
                                                     & QSDischargeShort$date <= as.POSIXct("2014-12-31"),]

meta <- metadata(constituent="Fluoride_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(Fluoride_mg_L ~ model(1), data=QS_chemF,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))


lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemF, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcFConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort2, se.pred=TRUE, date=TRUE)

aggs_lcFConc <- aggregateSolute(preds_lcFConc,meta, format="conc", se.preds=preds_lcFConc$se.pred, agg.by="calendar year")

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemF, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrF <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort2, se.pred=TRUE, date=TRUE)

#get some really high wierd values just remove them. 
#preds_lrF <-preds_lrF[preds_lrF$fit < 50,]

aggs_lrF <- aggregateSolute(preds_lrF,meta, format="flux total", se.preds=preds_lrF$se.pred, agg.by="calendar year")
aggs_lcFConc <- aggregateSolute(preds_lcFConc,meta, format="conc", se.preds=preds_lcFConc$se.pred, agg.by="calendar year")
aggs_lrF$F_Kg_Ha_yr <- aggs_lrF$Flux_Total / QSWatershedArea
aggs_lrF$F_mg_l <- aggs_lcFConc$Conc
FluxTotalsdfF <- aggs_lrF[, c('F_mg_l','F_Kg_Ha_yr','Calendar_Year')]

# round values
years <- FluxTotalsdfF$Calendar_Year
FluxTotalsdfF$Calendar_Year <- NULL
FluxTotalsdfF <- round(FluxTotalsdfF, 3)
FluxTotalsdfF$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfF,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals

QS_chemSiO2 <- QS_chem[complete.cases(QS_chem["SiO2_mg_L"]),]
QS_chemSiO2 <- QS_chemSiO2[complete.cases(QS_chemSiO2["CFS"]),]

meta <- metadata(constituent="SiO2_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(SiO2_mg_L ~ model(1), data=QS_chemSiO2,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemSiO2, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lcSiO2Conc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcSiO2Conc <- aggregateSolute(preds_lcSiO2Conc,meta, format="conc", se.preds=preds_lcSiO2Conc$se.pred, agg.by="calendar year")




lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemSiO2, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lrSiO2 <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrSiO2 <- aggregateSolute(preds_lrSiO2,meta, format="flux total", se.preds=preds_lrSiO2$se.pred, agg.by="calendar year")

aggs_lrSiO2$SiO2_Kg_Ha_yr <- aggs_lrSiO2$Flux_Total / QSWatershedArea
aggs_lrSiO2$SiO2_mg_l <- aggs_lcSiO2Conc$Conc
FluxTotalsdfSiO2 <- aggs_lrSiO2[, c('SiO2_mg_l','SiO2_Kg_Ha_yr','Calendar_Year')]


#round values 
years <- FluxTotalsdfSiO2$Calendar_Year
FluxTotalsdfSiO2$Calendar_Year <- NULL
FluxTotalsdfSiO2 <- round(FluxTotalsdfSiO2, 2)
FluxTotalsdfSiO2$Calendar_Year <- years


FluxTotals <- merge(FluxTotalsdfSiO2,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals


#PO4
QS_chemPO4 <- QS_chem[complete.cases(QS_chem["PO4.P_ug_P_L"]),]
QS_chemPO4 <- QS_chemPO4[complete.cases(QS_chemPO4["CFS"]),]

meta <- metadata(constituent="PO4.P_ug_P_L", flow="CFS", 
                 dates="date", conc.units="ug/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(PO4.P_ug_P_L ~ model(1), data=QS_chemPO4,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="ug/L", load.units="kg"))
lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemPO4, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcPO4Conc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcPO4Conc <- aggregateSolute(preds_lcPO4Conc,meta, format="conc", se.preds=preds_lcPO4Conc$se.pred, agg.by="calendar year")

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemPO4, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrPO4 <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrPO4 <- aggregateSolute(preds_lrPO4,meta, format="flux total", se.preds=preds_lrPO4$se.pred, agg.by="calendar year")

aggs_lrPO4$PO4_Kg_Ha_yr <- aggs_lrPO4$Flux_Total / QSWatershedArea
aggs_lrPO4$PO4_mg_l <- aggs_lcPO4Conc$Conc / 1000.0
FluxTotalsdfPO4 <- aggs_lrPO4[, c('PO4_mg_l','PO4_Kg_Ha_yr','Calendar_Year')]


#round values 
years <- FluxTotalsdfPO4$Calendar_Year
FluxTotalsdfPO4$Calendar_Year <- NULL
FluxTotalsdfPO4$PO4_mg_l <- round(FluxTotalsdfPO4$PO4_mg_l, 3)
FluxTotalsdfPO4$PO4_Kg_Ha_yr <- round(FluxTotalsdfPO4$PO4_Kg_Ha_yr, 2)

FluxTotalsdfPO4$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfPO4,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals

#NH4 

QS_chemNH4 <- QS_chem[complete.cases(QS_chem["NH4.N_ug_L"]),]
QS_chemNH4 <- QS_chemNH4[complete.cases(QS_chemNH4["CFS"]),]

meta <- metadata(constituent="NH4.N_ug_L", flow="CFS", 
                 dates="date", conc.units="ug/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(NH4.N_ug_L ~ model(1), data=QS_chemNH4,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="ug/L", load.units="kg"))
lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemNH4, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcNH4Conc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcNH4Conc <- aggregateSolute(preds_lcNH4Conc,meta, format="conc", se.preds=preds_lcNH4Conc$se.pred, agg.by="calendar year")

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemNH4, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrNH4 <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrNH4 <- aggregateSolute(preds_lrNH4,meta, format="flux total", se.preds=preds_lrNH4$se.pred, agg.by="calendar year")

aggs_lrNH4$NH4_Kg_Ha_yr <- aggs_lrNH4$Flux_Total / QSWatershedArea
aggs_lrNH4$NH4_mg_l <- aggs_lcNH4Conc$Conc
FluxTotalsdfNH4 <- aggs_lrNH4[, c('NH4_mg_l','NH4_Kg_Ha_yr','Calendar_Year')]


#round values 
years <- FluxTotalsdfNH4$Calendar_Year
FluxTotalsdfNH4$Calendar_Year <- NULL
FluxTotalsdfNH4$NH4_mg_l <- round(FluxTotalsdfNH4$NH4_mg_l, 3)
FluxTotalsdfNH4$NH4_Kg_Ha_yr <- round(FluxTotalsdfNH4$NH4_Kg_Ha_yr, 2)
FluxTotalsdfNH4$Calendar_Year <- years


FluxTotals <- merge(FluxTotalsdfNH4,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals

QS_chemTDN <- QS_chem[complete.cases(QS_chem["TDN_mg_L"]),]
QS_chemTDN <- QS_chemTDN[complete.cases(QS_chemTDN["CFS"]),]

meta <- metadata(constituent="TDN_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(TDN_mg_L ~ model(1), data=QS_chemTDN,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemTDN, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrTDN <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrTDN <- aggregateSolute(preds_lrTDN,meta, format="flux total", se.preds=preds_lrTDN$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemTDN, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcTDNConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcTDNConc <- aggregateSolute(preds_lcTDNConc,meta, format="conc", se.preds=preds_lcTDNConc$se.pred, agg.by="calendar year")

aggs_lrTDN$TDN_Kg_Ha_yr <- aggs_lrTDN$Flux_Total / QSWatershedArea
aggs_lrTDN$TDN_mg_l <- aggs_lcTDNConc$Conc
FluxTotalsdfTDN <- aggs_lrTDN[, c('TDN_mg_l','TDN_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfTDN$Calendar_Year
FluxTotalsdfTDN$Calendar_Year <- NULL
FluxTotalsdfTDN <- round(FluxTotalsdfTDN, 2)
FluxTotalsdfTDN$Calendar_Year <- years


FluxTotals <- merge(FluxTotalsdfTDN,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals


#DON
QS_chemDON <- QS_chem[complete.cases(QS_chem["DON_mg_L"]),]
QS_chemDON <- QS_chemDON[complete.cases(QS_chemDON["CFS"]),]

meta <- metadata(constituent="DON_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(DON_mg_L ~ model(1), data=QS_chemDON,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemDON, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrDON <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrDON <- aggregateSolute(preds_lrDON,meta, format="flux total", se.preds=preds_lrDON$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemDON, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcDONConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcDONConc <- aggregateSolute(preds_lcDONConc,meta, format="conc", se.preds=preds_lcDONConc$se.pred, agg.by="calendar year")


aggs_lrDON$DON_Kg_Ha_yr <- aggs_lrDON$Flux_Total / QSWatershedArea
aggs_lrDON$DON_mg_l <- aggs_lcDONConc$Conc
FluxTotalsdfDON <- aggs_lrDON[, c('DON_mg_l','DON_Kg_Ha_yr','Calendar_Year')]


#round values 
years <- FluxTotalsdfDON$Calendar_Year
FluxTotalsdfDON$Calendar_Year <- NULL
FluxTotalsdfDON <- round(FluxTotalsdfDON, 2)
FluxTotalsdfDON$Calendar_Year <- years


FluxTotals <- merge(FluxTotalsdfDON,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals


#DOC
QS_chemDOC <- QS_chem[complete.cases(QS_chem["DOC_mg_L"]),]
QS_chemDOC <- QS_chemDOC[complete.cases(QS_chemDOC["CFS"]),]

meta <- metadata(constituent="DOC_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(DOC_mg_L ~ model(1), data=QS_chemDOC,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemDOC, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lrDOC <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrDOC <- aggregateSolute(preds_lrDOC,meta, format="flux total", se.preds=preds_lrDOC$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemDOC, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lcDOCConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcDOCConc <- aggregateSolute(preds_lcDOCConc,meta, format="conc", se.preds=preds_lcDOCConc$se.pred, agg.by="calendar year")

aggs_lrDOC$DOC_Kg_Ha_yr <- aggs_lrDOC$Flux_Total / QSWatershedArea
aggs_lrDOC$DOC_mg_l <- aggs_lcDOCConc$Conc
FluxTotalsdfDOC <- aggs_lrDOC[, c('DOC_mg_l','DOC_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfDOC$Calendar_Year
FluxTotalsdfDOC$Calendar_Year <- NULL
FluxTotalsdfDOC <- round(FluxTotalsdfDOC, 2)
FluxTotalsdfDOC$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfDOC,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals


#Cl
QS_chemCl <- QS_chem[complete.cases(QS_chem["Cl_mg_L"]),]
QS_chemCl <- QS_chemCl[complete.cases(QS_chemCl["CFS"]),]

meta <- metadata(constituent="Cl_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(Cl_mg_L ~ model(1), data=QS_chemCl,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemCl, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lrCl <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrCl <- aggregateSolute(preds_lrCl,meta, format="flux total", se.preds=preds_lrCl$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemCl, abs.or.rel.resids="relative", interp.function=rectangularInterpolation)

preds_lcClConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lcClConc <- aggregateSolute(preds_lcClConc,meta, format="conc", se.preds=preds_lcClConc$se.pred, agg.by="calendar year")


aggs_lrCl$Cl_Kg_Ha_yr <- aggs_lrCl$Flux_Total / QSWatershedArea
aggs_lrCl$Cl_mg_l <- aggs_lcClConc$Conc
FluxTotalsdfCl <- aggs_lrCl[, c('Cl_mg_l','Cl_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfCl$Calendar_Year
FluxTotalsdfCl$Calendar_Year <- NULL
FluxTotalsdfCl <- round(FluxTotalsdfCl, 2)
FluxTotalsdfCl$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfCl,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals

QS_chemK <- QS_chem[complete.cases(QS_chem["K_mg_L"]),]
QS_chemK <- QS_chemK[complete.cases(QS_chemK["CFS"]),]

meta <- metadata(constituent="K_mg_L", flow="CFS", 
                 dates="date", conc.units="mg/L", flow.units="cfs", load.units="kg", 
                 load.rate.units="kg d^-1", station="QS, PR")


lr <- loadReg2(loadReg(K_mg_L ~ model(1), data=QS_chemK,
                       flow="CFS", dates="date", time.step="instantaneous", 
                       flow.units="cfs", conc.units="mg/L", load.units="kg"))

lc <- loadComp(reg.model=lr, interp.format="flux", 
               interp.data=QS_chemK, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lrK <- predictSolute(lc, "flux",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)

aggs_lrK <- aggregateSolute(preds_lrK,meta, format="flux total", se.preds=preds_lrK$se.pred, agg.by="calendar year")

lcConc <- loadComp(reg.model=lr, interp.format="conc", 
                   interp.data=QS_chemK, abs.or.rel.resids="absolute", interp.function=rectangularInterpolation)

preds_lcKConc <- predictSolute(lcConc, "conc",interval="prediction", QSDischargeShort, se.pred=TRUE, date=TRUE)


aggs_lrK$K_Kg_Ha_yr <- aggs_lrK$Flux_Total / QSWatershedArea
aggs_lcKConc <- aggregateSolute(preds_lcKConc,meta, format="conc", se.preds=preds_lcKConc$se.pred, agg.by="calendar year")
aggs_lrK$K_mg_l <- aggs_lcKConc$Conc
FluxTotalsdfK <- aggs_lrK[, c('K_mg_l','K_Kg_Ha_yr','Calendar_Year')]

#round values 
years <- FluxTotalsdfK$Calendar_Year
FluxTotalsdfK$Calendar_Year <- NULL
FluxTotalsdfK <- round(FluxTotalsdfK, 2)
FluxTotalsdfK$Calendar_Year <- years

FluxTotals <- merge(FluxTotalsdfK,FluxTotals, by = "Calendar_Year", all.y=TRUE)
FluxTotals

rownames(FluxTotals) <- FluxTotals$Calendar_Year
FluxTotals2 <- as.data.frame(t(FluxTotals))
FluxTotals2$Calendar_Year <- NULL
FluxTotals2