Merge branch 'master' into master

DOI-USGS · Aug 9, 2017 · b2f03c0 · b2f03c0
2 parents f48141d + f18046b
commit b2f03c0
Show file tree

Hide file tree

Showing 34 changed files with 356 additions and 869 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -18,3 +18,4 @@ test_stats_huc12_modeled_performance.R
 appveyor.yml
 EflowStats.bbl
 ^appveyor\.yml$
+CONDUCT.md
diff --git a/CONDUCT.md b/CONDUCT.md
@@ -0,0 +1,25 @@
+# Contributor Code of Conduct
+
+As contributors and maintainers of this project, we pledge to respect all people who 
+contribute through reporting issues, posting feature requests, updating documentation,
+submitting pull requests or patches, and other activities.
+
+We are committed to making participation in this project a harassment-free experience for
+everyone, regardless of level of experience, gender, gender identity and expression,
+sexual orientation, disability, personal appearance, body size, race, ethnicity, age, or religion.
+
+Examples of unacceptable behavior by participants include the use of sexual language or
+imagery, derogatory comments or personal attacks, trolling, public or private harassment,
+insults, or other unprofessional conduct.
+
+Project maintainers have the right and responsibility to remove, edit, or reject comments,
+commits, code, wiki edits, issues, and other contributions that are not aligned to this 
+Code of Conduct. Project maintainers who do not follow the Code of Conduct may be removed 
+from the project team.
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by 
+opening an issue or contacting one or more of the project maintainers.
+
+This Code of Conduct is adapted from the Contributor Covenant 
+(http:contributor-covenant.org), version 1.0.0, available at 
+http://contributor-covenant.org/version/1/0/0/
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,11 +1,12 @@
 Package: EflowStats
 Type: Package
-Title: Returns Hydrologic Indicator stats for a given set of data
-Version: 1.0
-Date: 2017-04-04
-Authors@R: c( person("Joe", "Mills", role = c("aut","cre"),
+
+Title: Hydrologic Indicator and Alteration Stats
+Version: 5.0.0
+Date: 2017-08-04
+Authors@R: c( person("Joe", "Mills", role = c("aut"),
     email = "jmills@ucar.edu"),
-              person("David", "Blodgett", role = c("aut"),
+              person("David", "Blodgett", role = c("aut","cre"),
     email = "dblodgett@usgs.gov"))
 Depends:
     R (>= 3.3)
@@ -23,7 +24,7 @@ Suggests:
     DT,
     dataRetrieval,
     readr
-Description: EflowStats is a reimplementation of the Hydrologic Index Tool 
+Description: A reimplementation of the Hydrologic Index Tool 
         (HIT; Henriksen et al, 2006) for calculating 171 hydrologic indices 
         for stream classification analysis. EflowStats also calculates 7 
         additional statistics used for streamflow classification reffered 

diff --git a/EflowStats.Rproj b/EflowStats.Rproj
@@ -6,7 +6,7 @@ AlwaysSaveHistory: Default
 
 EnableCodeIndexing: Yes
 UseSpacesForTab: Yes
-NumSpacesForTab: 8
+NumSpacesForTab: 2
 Encoding: UTF-8
 
 RnwWeave: Sweave
@@ -15,4 +15,4 @@ LaTeX: pdfLaTeX
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --with-keep.source
-PackageRoxygenize: rd,collate,namespace,vignette
+PackageRoxygenize: rd,collate,namespace
diff --git a/R/00_EflowStats.R b/R/00_EflowStats.R
@@ -21,10 +21,18 @@ month_val <- discharge <- leapYear <- day <-
         roll30Mean <- roll90Mean <- events <- 
         . <- highFlow <- CV <- block <- 
         peakQ <- flow <- maxFlow <-  
-        minFlow <- event <- ".dplyr.var"
+        minFlow <- event <- calendar_day <- 
+        ".dplyr.var"
 
 check_preference <- function(pref) {
         if(!pref %in% c("mean", "median")){ 
                 stop("Preference must be either mean or median") 
         }
+}
+
+.onAttach <- function(libname, pkgname) {
+        packageStartupMessage(paste(strwrap(
+                'USGS Research Package: 
+                https://owi.usgs.gov/R/packages.html#research'),
+                collapse='\n'))
 }
diff --git a/R/calc_timingHigh.R b/R/calc_timingHigh.R
@@ -19,6 +19,8 @@
 #' flow is less than the 1.67-year flood threshold and also occurs in all years. Accumulate nonflood days that 
 #' span all years. TH3 is maximum length of those flood-free periods divided by 365.  
 #' }
+#' Note: In these definitions, "Julian date" should be interpreted as the count of days starting with 1 on January 
+#' first of a given year, ending at 365 or 366 on December 31st.
 #' @return A data.frame of flow statistics
 #' @importFrom lubridate year
 #' @importFrom lubridate month
@@ -36,30 +38,33 @@ calc_timingHigh <- function(x,yearType = "water",digits=3,pref="mean",floodThres
 
         #calculate some stuff for use later
         x$month_val <- lubridate::month(x$date)
-        #x$jul_day <- lubridate::yday(x$date)
+        x$calendar_day <- lubridate::yday(x$date)
 
+        # Calculate max flow value and calendat day of max flow for every year.
         flowSum_year <- dplyr::summarize(dplyr::group_by(x,year_val),
                                             maxFlow = max(discharge),
-                                         maxFlowJulDay = min(day[discharge==max(discharge)])
+                                         maxFlowJulDay = min(calendar_day[discharge==max(discharge)])
                                          )
 
-        #th1
+        #th1 Convert the max flow julian day to psuedo-radians and take cos (x) and sin (y). 
         flowSum_year$np <- cos(flowSum_year$maxFlowJulDay*2*pi/365.25)
         flowSum_year$mdata <- sin(flowSum_year$maxFlowJulDay*2*pi/365.25)
 
+        # Average accross all years.
         xbar <- mean(flowSum_year$np)
         ybar <- mean(flowSum_year$mdata)
-        if (xbar>0) {
+        if (xbar>0) { # if x component is greater than 0, just return the angle.
                 th1_temp <- atan(ybar/xbar)*180/pi
-        } else if (xbar<0) {
-                th1_temp <- (atan(ybar/xbar)*180/pi)+180
+        } else if (xbar<0) {# if x component is less than 0 (angles pi/2 to 3pi/2), arctan returns -90 -> 0 -> 90 
+                th1_temp <- (atan(ybar/xbar)*180/pi)+180 # but we should return 90 -> 180 -> 270 so add 180 to output.
         } else if (xbar==0 && ybar>0) {
                 th1_temp <- 90
         } else if (xbar==0 && ybar<0) {
                 th1_temp <- 270
         }
 
-        th1_temp <- ifelse(th1_temp<0,th1_temp+360,th1_temp)
+        # if th1_temp is negative from above (for 270 or 3pi/2 to 360 or 2pi degrees) just add 360.
+        th1_temp <- ifelse(th1_temp<0,th1_temp+360,th1_temp) 
         th1 <- th1_temp*365.25/360
 
         #th2

diff --git a/R/calc_timingLow.R b/R/calc_timingLow.R
@@ -23,6 +23,8 @@
 #' threshold as the ratio of number of days to 365 or 366 (leap year) for each year. TL4 is the maximum of the yearly 
 #' ratios.
 #' }
+#' Note: In these definitions, "Julian date" should be interpreted as the count of days starting with 1 on January 
+#' first of a given year, ending at 365 or 366 on December 31st.
 #' @return A data.frame of flow statistics
 #' @importFrom lubridate year
 #' @importFrom lubridate month
@@ -40,30 +42,36 @@ calc_timingLow <- function(x,yearType = "water",digits=3,pref="mean",floodThresh
 
         #calculate some stuff for use later
         x$month_val <- lubridate::month(x$date)
+        x$calendar_day <- lubridate::yday(x$date)
 
+        # Calculate min flow value and julian date of min flow for every year.
         flowSum_year <- dplyr::summarize(dplyr::group_by(x,year_val),
                                          minFlow = min(discharge),
-                                         minFlowJulDay = min(day[discharge==min(discharge)])
+                                         minFlowJulDay = min(calendar_day[discharge==min(discharge)])
         )
 
-        #tl1
+        #tl1 Convert the min flow julian day to psuedo-radians and take cos (x) and sin (y). 
         flowSum_year$np <- cos(flowSum_year$minFlowJulDay*2*pi/365.25)
         flowSum_year$mdata <- sin(flowSum_year$minFlowJulDay*2*pi/365.25)
 
+        # Average accross all years.
         xbar <- mean(flowSum_year$np)
         ybar <- mean(flowSum_year$mdata)
-        if (xbar>0) {
+        if (xbar>0) { # if x component is greater than 0, just return the angle.
                 tl1_temp <- atan(ybar/xbar)*180/pi
-        } else if (xbar<0) {
-                tl1_temp <- (atan(ybar/xbar)*180/pi)+180
-        } else if (xbar==0 && ybar>0) {
+        } else if (xbar<0) { # if x component is less than 0 (angles pi/2 to 3pi/2), arctan returns -90 -> 0 -> 90 
+                tl1_temp <- (atan(ybar/xbar)*180/pi)+180 # but we should return 90 -> 180 -> 270 so add 180 to output.
+        } else if (xbar==0 && ybar>0) { # discontinuity
                 tl1_temp <- 90
-        } else if (xbar==0 && ybar<0) {
+        } else if (xbar==0 && ybar<0) { # discontinuity
                 tl1_temp <- 270
         }
 
+        # if tl1_temp is negative from above (for 270 or 3pi/2 to 360 or 2pi degrees) just add 360.
         tl1_temp <- ifelse(tl1_temp<0,tl1_temp+360,tl1_temp)
-        tl1 <- tl1_temp*365.25/360
+        tl1 <- round(tl1_temp*365.25/360)
+
+        if(tl1 == 0) tl <- 365.25 
 
         #tl2
         tl2_a <- sqrt((xbar*xbar)+(ybar*ybar))

diff --git a/R/get_peakThreshold.R b/R/get_peakThreshold.R
@@ -2,17 +2,24 @@
 #' 
 #' This function calculates flood thresholds for specified recurrence intervals. 
 #' 
-#' @param x A dataframe containing a vector of date values in the first column and vector of numeric flow values in the second column.
-#' @param peakValues A dataframe containing a vector of date values in the first column and vector of numeric annual peak flow values in the second column.
-#' @param yearType A charcter of either "water" or "calendar" indicating whether to use water years or calendar years, respectively.
+#' @param x A dataframe containing a vector of date values in the first column 
+#' and vector of numeric flow values in the second column.
+#' @param peakValues A dataframe containing a vector of date values in the first 
+#' column and vector of numeric annual peak flow values in the second column.
+#' @param yearType A charcter of either "water" or "calendar" indicating whether 
+#' to use water years or calendar years, respectively.
 #' @param perc value containing the desired percentile to be calculated
-#' @return thresh numeric containing the flood threshold for recurence interval specified by \code{perc}
+#' @return thresh numeric containing the flood threshold for recurence interval 
+#' specified by \code{perc}
 #' @details Compute the log10 of the daily 
-#' flows for the peak annual flow days. Calculate the coefficients for a linear regression equation for logs of peak 
-#' annual flow versus logs of average daily flow for peak days. Using the log peak flow for the 1.67-year recurrence 
-#' interval (60th percentile, \code{perc=0.6}) as input to the regression equation, predict the log10 of the average daily flow. The 
-#' threshold is 10 to the log10 (average daily flow) power (cubic feet per second). for the 5-year recurrence interval 
-#' (80th percentile, \code{perc=0.8}), used by indices TL3 and TL4, follow the same process, inputing a different 'perc' value.
+#' flows for the peak annual flow days. Calculate the coefficients for a linear
+#' regression equation for logs of peak annual flow versus logs of average daily
+#' flow for peak days. Using the log peak flow for the 1.67-year recurrence 
+#' interval (60th percentile, \code{perc=0.6}) as input to the regression
+#' equation, predict the log10 of the average daily flow. The threshold is 10 to
+#' the log10 (average daily flow) power (cubic feet per second). for the 5-year
+#' recurrence interval (80th percentile, \code{perc=0.8}), used by indices TL3
+#' and TL4, follow the same process, inputing a different 'perc' value.
 #' @importFrom stats na.omit quantile
 #' @export
 #' @examples
@@ -25,73 +32,88 @@
 #' get_peakThreshold(x,peakValues,.6,yearType="water")
 #' }
 get_peakThreshold <- function(x,peakValues,perc=0.6,yearType = "water") {
-
-        #Check x data input
-        x <- validate_data(x,yearType)
-
-        #Check peakValues data input
-        if (nrow(peakValues)<=1) {
-                stop("peakValues must have a minimum of two annual values")
-        }
-
-        col1_class = class(peakValues[,1])
-        col2_class = class(peakValues[,2])
-
-        if(col1_class != "Date" && col2_class != "numeric")
-        {
-                stop("First column of peakValues must contain a vector of class date.\nSecond column of peakValues must contain a vector of class numeric.") 
-        } else if (col1_class != "Date")
-        {
-                stop("First column of peakValues must contain a vector of class date.") 
-        } else if (col2_class != "numeric" & col2_class != "integer")
-        {
-                stop("Second column of peakValues must contain a vector of class numeric.") 
-        }
-
-        if(anyNA(peakValues))
-        {
-                stop("dataframe peakValues cannot contain NA values")
-        }
-
-
-        names(peakValues) <- c("date","peakQ")
-
-
-        if(yearType == "water")
-        {
-                peakValues$year_val <- get_waterYear(peakValues$date)
-        } else {
-                peakValues$year_val <- lubridate::year(peakValues$date)
-        }
-
-        if(any(duplicated(peakValues$year_val))) {
-                warning("peakValues data frame contains multiple peak values for one or more years. Only the maximum annual value will be retained.")
-                temp <- dplyr::summarize(dplyr::group_by(peakValues,year_val),
-                                         date = date[peakQ == max(peakQ)],
-                                         peakQ = max(peakQ))
-                peakValues <- dplyr::left_join(peakValues["date"],temp, by="date")
-                peakValues <- na.omit(peakValues)
-        }
-
-        peakDaily <- dplyr::inner_join(x,peakValues[c("date","peakQ")],by="date")
-        peakDaily$logQ <- log10(peakDaily$discharge)
-        peakDaily$logPeakQ <- log10(peakDaily$peakQ)
-
-        dailyMean <- mean(peakDaily$logQ)
-        instMean <- mean(peakDaily$logPeakQ)
-
-        num <- sum((peakDaily$logQ-dailyMean)*(peakDaily$logPeakQ-instMean))
-        denom <- sum((peakDaily$logPeakQ-instMean)*(peakDaily$logPeakQ-instMean))
-
-        b <- num/denom
-        a <- dailyMean - (b*instMean)
-
-        lfcrit <- quantile(peakDaily$logPeakQ,
-                           probs = perc,
-                           type=6,
-                           names = F)
-
-        lq167 <- a + (b*as.numeric(lfcrit))
-        thresh <- 10^(lq167)
-        return(thresh)
+
+  #Check x data input
+  x <- validate_data(x,yearType)
+
+  #Check peakValues data input
+  if (nrow(peakValues)<=1) {
+    stop("peakValues must have a minimum of two annual values")
+  }
+
+  col1_class = class(peakValues[,1])
+  col2_class = class(peakValues[,2])
+
+  if(col1_class != "Date" && col2_class != "numeric") {
+    stop(paste("First column of peakValues must contain a vector of class",
+               "date.\nSecond column of peakValues must contain a vector of", 
+               "class numeric."))
+  } else if (col1_class != "Date") {
+    stop("First column of peakValues must contain a vector of class date.") 
+  } else if (col2_class != "numeric" & col2_class != "integer") {
+    stop("Second column of peakValues must contain a vector of class numeric.") 
+  }
+
+  if(anyNA(peakValues)) stop("dataframe peakValues cannot contain NA values")
+
+  names(peakValues) <- c("date","peakQ")
+
+  if(yearType == "water") {
+    peakValues$year_val <- get_waterYear(peakValues$date)
+  } else {
+    peakValues$year_val <- lubridate::year(peakValues$date)
+  }
+
+  if(any(duplicated(peakValues$year_val))) { 
+    warning(paste("peakValues data frame contains multiple peak values", 
+                  "for one or more years. Only the maximum annual value", 
+                  "will be retained."))
+    temp <- dplyr::summarize(dplyr::group_by(peakValues,year_val),
+                             date = date[peakQ == max(peakQ)],
+                             peakQ = max(peakQ))
+    # temp contains the peakQ and date of peakQ for each year.
+    # The following eliminates ones that are less than max.
+    peakValues <- dplyr::left_join(peakValues["date"],temp, by="date")
+    peakValues <- na.omit(peakValues)
+  }
+
+  # Get the peak daily for the peak flow days. (peakQ is instantaneous!)
+  peakDaily <- dplyr::inner_join(x,peakValues[c("date","peakQ")],by="date")
+
+  #Get logs of both instantaneous and daily peak flows.
+  peakDaily$logQ <- log10(peakDaily$discharge)
+  peakDaily$logPeakQ <- log10(peakDaily$peakQ)
+
+  # Get mean of logs of flows.
+  dailyMean <- mean(peakDaily$logQ)
+  instMean <- mean(peakDaily$logPeakQ)
+
+  # Calculate the coefficients for a linear regression equation for logs 
+  # of peak annual flow versus logs of average daily flow for peak days.
+
+  # ... in other words
+  # numerator is the sum of the difference between the annual series values 
+  # and POR mean peaks for daily * instantaneous peaks
+  # Denominator is the sum of the difference between the annual series values 
+  # and POR mean peaks for instantaneous^2 peaks
+  num <- sum((peakDaily$logQ-dailyMean)*(peakDaily$logPeakQ-instMean))
+  denom <- sum((peakDaily$logPeakQ-instMean)*(peakDaily$logPeakQ-instMean))
+
+  # Calculates coeficients of a linear regression based on the above.
+  b <- num/denom
+  a <- dailyMean - (b*instMean)
+
+  # Using the log peak flow for the 1.67-year recurrence interval 
+  # (60th percentile) as input to the regression equation, 
+
+  lfcrit <- quantile(peakDaily$logPeakQ,
+                     probs = perc,
+                     type=6,
+                     names = F)
+
+  # predict the log10 of the average daily flow.
+  lq167 <- a + (b*as.numeric(lfcrit))
+  thresh <- 10^(lq167)
+
+  return(thresh)
 }