Updates for instantaneous vignette

Author: dlorenz-usgs

R/predConc.R

@@ -9,7 +9,9 @@
 #' If \code{allow.incomplete} is \code{TRUE}, then concentrations will be 
 #'computed based on all nonmissing values, otherwise missing values
 #'\code{NAs} will be returned. For this application, missing values 
-#'includes \code{NAs} and incomplete days.
+#'includes \code{NAs} and incomplete days. For prediction by "day" when 
+#'there are variable number of unit values per day, \code{allow.incomplete}
+#'must be set to \code{TRUE}.
 #'
 #' The term confidence interval is used here as in the original 
 #'documentation for LOADEST, but the values that are reported are 
@@ -211,6 +213,9 @@ predConc <- function(fit, newdata, by="day",
       Kdy <- as.integer(KDate)
       KDate <- unique(KDate)
       Kdy <- Kdy - Kdy[1L] + 1L # make relative to first day (Index)
+      if(length(unique(table(Kdy))) > 1L && !allow.incomplete) {
+        warning("Variable observations per day, either set the allow.incomplete argument to TRUE or use the resampleUVdata function to construct a uniform series")
+      }
       KinAll <- unique(Kdy)
       ## Make it daily flow, Flow0 indicates a partial 0 flow
       Flow0 <- tapply(Flow, Kdy, min) 

R/predLoad.R

@@ -14,7 +14,9 @@
 #'\code{NAs} will be returned. For this application, missing values 
 #'includes \code{NAs} and gaps in the record, except for \code{by} 
 #'set to "total" or user defined groups where missing values only
-#'includes \code{NAs}.
+#'includes \code{NAs}. For prediction by "day" when there are variable
+#'number of unit values per day, \code{allow.incomplete} must be
+#'set to \code{TRUE}.
 #'
 #' The term confidence interval is used here as in the original 
 #'documentation for LOADEST, but the values that are reported are 
@@ -52,7 +54,6 @@
 #' @useDynLib rloadest estlday
 #' @useDynLib rloadest estltot
 #' @export
-
 predLoad <- function(fit, newdata, load.units=fit$load.units, by="total", 
                      seopt="exact", allow.incomplete=FALSE, 
                      conf.int=0.95, print=FALSE) {
@@ -239,6 +240,9 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
       Kdy <- as.integer(KDate)
       KDate <- unique(KDate)
       Kdy <- Kdy - Kdy[1L] + 1L # make relative to first day (Index)
+      if(length(unique(table(Kdy))) > 1L && !allow.incomplete) {
+        warning("Variable observations per day, either set the allow.incomplete argument to TRUE or use the resampleUVdata function to construct a uniform series")
+      }
       KinAll <- unique(Kdy)
       ## Make it daily flow, Flow0 indicates a partial 0 flow
       Flow0 <- tapply(Flow, Kdy, min) 
@@ -470,14 +474,19 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
     cat("Streamflow Summary Statistics\n",
         "-------------------------------------------\n\n", sep="")
     Qsum <- rbind(Cal.=fit$Sum.flow, Est.=summary(Flow))
-    if(diff(range(Qsum[, ncol(Qsum)])) > 0)
+    if(Qsum[2L, ncol(Qsum)] > Qsum[1L, ncol(Qsum)]) {
       cat("WARNING: The maximum estimation data set steamflow exceeds the maximum\n",
           "calibration data set streamflow. Load estimates require extrapolation.\n\n",
           sep="")
-    else
+    } else if(Qsum[2L, 1L] < Qsum[1L, 1L]) {
+	   cat("WARNING: The minimum estimation data set steamflow exceeds the minimum\n",
+          "calibration data set streamflow. Load estimates require extrapolation.\n\n",
+          sep="")
+	} else {
       cat("The maximum estimation data set streamflow does not exceed the maximum\n",
           "calibration data set streamflow. No extrapolation is required.\n\n",
           sep="")
+	}
     if(!(by %in% c("day", "unit"))) {
       cat("\n-------------------------------------------------------------\n",
           "Constituent Output File Part IIb: Estimation (Load Estimates)\n",

R/print.loadReg.R

@@ -116,7 +116,7 @@ print.loadReg <- function(x, digits=4, brief=TRUE, load.only=brief, ...) {
   else
     pval <- format(round(pval, 4), scientific=5)
   ## Compute the PPCC
-  Res <- residuals(x$lfit, type="working")
+  Res <- residuals(x$lfit, type="working", suppress.na.action=TRUE)
   ppcc <- censPPCC.test(as.lcens(Res, censor.codes=x$lfit$CENSFLAG))
   cat("\n", x$method, " Regression Statistics\n",
       "Residual variance: ", 

R/rloadest-package.R

@@ -27,17 +27,11 @@
 #'Furhtermore, the model building capability in the rloadest functions make easier
 #'to explore other forms of rating-curve models than LOADEST.
 #'
-#' @name LOADEST-package
+#' @name rloadest-package
 #' @docType package
 #' @author Dave Lorenz \email{lorenz@@usgs.gov}
 #' @keywords load estimation
 NULL
-
-#' Example Atrazine data included in LOADEST package
-#'
-#' Example data representing atrazine
-#'
-#' @name Atrazine
-#' @docType data
-#' @keywords water quality data
-NULL
+.onAttach <- function(libname, pkgname) {
+  packageStartupMessage("Although this software program has been used by the U.S. Geological Survey (USGS), no warranty, expressed or implied, is made by the USGS or the U.S. Government as to the accuracy and functioning of the program and related program material nor shall the fact of distribution constitute any such warranty, and no responsibility is assumed by the USGS in connection therewith.")
+}

demo/00Index

@@ -1 +1,2 @@
 BatchLoad	A script for batch data retrieval and load estimation. 
+ChangeLoadGraph	Create a surface plot showing the change in loading between two years.

demo/ChangeLoadGraph.r

@@ -0,0 +1,56 @@
+# Use app1 (model7) from the vignette
+library(rloadest)
+data(app1.calib)
+# construct the model, see the vignette for app1 for details
+app1.lr7 <- loadReg(Phosphorus ~ model(7), data = app1.calib, flow = "FLOW",
+  dates = "DATES", conc.units="mg/L",
+  station="Illinois River at Marseilles, Ill.")
+# construct the grids
+FLOW <- seq(3000, 60000, length.out=21)
+DATES75 <- seq(as.Date("1975-01-01"), as.Date("1976-01-01"), length.out=24)
+Tmp75 <- expand.grid(DATES=DATES75, FLOW=FLOW)
+Tmp75$C <- predConc(app1.lr7, Tmp75)$Conc
+Tmp75$L <- predLoad(app1.lr7, Tmp75, by="day")$Flux
+Z.conc75 <- matrix(Tmp75$C, nrow=24)
+Z.load75 <- matrix(Tmp75$L, nrow=24)
+
+DATES84 <- seq(as.Date("1984-01-01"), as.Date("1985-01-01"), length.out=24)
+Tmp84 <- expand.grid(DATES=DATES84, FLOW=FLOW)
+Tmp84$C <- predConc(app1.lr7, Tmp84)$Conc
+Tmp84$L <- predLoad(app1.lr7, Tmp84, by="day")$Flux
+Z.conc84 <- matrix(Tmp84$C, nrow=24)
+Z.load84 <- matrix(Tmp84$L, nrow=24)
+
+# The range in concentration should be from .24 to .92
+AA.lev <- seq(.24, .92, by=.02)
+# the maximum load is 84,000+, so set range to 100,000
+# Construct the 1975 load/concentration surface
+# The shape of the surface indicates the potential loading given
+# The date and flow. The color of the surface indicates the 
+# potential concentration given the date and flow.
+preSurface(DATES75, FLOW, Z.load75, zaxis.range = c(0, 100000),
+ batch="I") -> AA75.pre
+# I selected
+# Construct the 1984 load/concentration surface
+preSurface(DATES84, FLOW, Z.load84, zaxis.range = c(0, 100000),
+ batch="I") -> AA84.pre
+# I selected
+# Proceed:
+#set the page to landscape
+setPDF("land", basename="Change_Load")
+setLayout(width=c(4.5,4.5), height=5.5, explanation=list(right=1.2)) -> AA.lo
+setGraph(1, AA.lo)
+surfacePlot(pre=AA75.pre, z.color=Z.conc75, 
+						Surface=list(name="Concentration", levels=AA.lev),
+						xtitle="1975", ytitle="Streamflow", ztitle="Load") -> AA75.pl
+addCaption("Change in Phosphorus Loading in the Illinois River at Marseilles, Ill. from 1975 to 1984")
+# Construct the 1984 load/concentration surface
+setGraph(2, AA.lo)
+surfacePlot(pre=AA84.pre, z.color=Z.conc84, 
+						Surface=list(name="Concentration", levels=AA.lev),
+						xtitle="1984", ytitle="Streamflow", ztitle="Load") -> AA84.pl
+# The 
+setGraph("explanation", AA.lo)
+addExplanation(AA84.pl)
+dev.off()
+