R/get_seasonality.R

In USGS-R/EflowStats: Hydrologic Indicator and Alteration Stats

#' Function to compute the seasonal factors (amplitude and phase) for a given data series
#' 
#' This function accepts a data frame containing daily streamflow data, then computes get_seasonality 
#' variables by first standardizing flows, the fitting relation 
#' A*cos(2*pi*t) + B*sin(2*pi*t)1) Get decimal yearand returns the amplitude and phase
#' 
#' @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 yearType A charcter of either "water" or "calendar" indicating whether to use water years or calendar years, respectively.
#' @param wyMonth A numeric. The month of the year in which the water year starts 
#' (1=January, 12=December). The water year begins on the first day of wyMonth.
#' @importFrom stats .lm.fit coef
#' @return get_seasonality vector of seasonal factors (amplitude and phase)
#' @examples 
#' x <- sampleData[c("date","discharge")]
#' get_seasonality(x=x)
#' @export
get_seasonality <- function(x,yearType = "water",wyMonth=10L) {
        
        #Check data inputs
        x <- validate_data(x,yearType=yearType,wyMonth=wyMonth)
        
        if(isFALSE(x)) stop("input data not valid")
        
        #calculate some stuff for use later
        x$month_val <- lubridate::month(x$date)
        
        rawdates<-x$date
        dateaschar<-as.character(rawdates)
        jday<-lubridate::yday(x$date)
        decimal_year<-as.numeric(x$year_val)+(jday/365.25)
        #2) Standardize flows
        std_flows<-scale(x$discharge, center = TRUE, scale = TRUE)
        #3) Use linear model to fit 
        x_mat = cbind(1, sin(2*pi*decimal_year), cos(2*pi*decimal_year))
        seasonfit<-.lm.fit(x_mat, std_flows)
        b1<-as.vector(coef(seasonfit)[2])
        b2<-as.vector(coef(seasonfit)[3]) 
        #Now compute the amplitude and phase of the seasonal signal
        amplitude<-round(sqrt((b2^2)+(b1^2)),digits=2)
        #phase<-round(atan((-seasonB)/seasonA),digits=2)
        MaxDay <- function(b1,b2){
                version1 <- 365.25 * ((pi/2)-atan(b2/b1))/(2*pi)
                version2 <- 365.25 * ((pi/2)-pi-atan(b2/b1))/(2*pi)
                MaxDay <- if(b1 > 0) version1 else 365.25 + version2
                MaxDay <- if(b1 == 0 & b2 > 0) 365.25 else MaxDay
                MaxDay <- if(b1 == 0 & b2 < 0) 365.25/2 else MaxDay
                MaxDay <- if(b1 == 0 & b2 == 0) NA else MaxDay
                return(MaxDay)
        }
        phase <- MaxDay(b1,b2)
        get_seasonalityv <- cbind(amplitude,phase)
        return(get_seasonalityv)
}