floodStat: FloodNet Regional Frequency Analysis

###########################################################
#' Extracts the Annual maxima of a daily time series
#'
#' Returns a dataset containing the annual maximums,
#' the date and the number of observed days during the year.
#'
#' @param form Formula of the form \code{value ~ date} that specifies the variable
#'   from which the annual maximums are extracted and a date variable.
#'
#' @param x Data containing the date and variable.
#'  If a \code{data.frame} is passed directly the first column must be the
#'  value and the second the date.
#'
#' @param tol Filter the return years a verify that the maxima is computed from
#'   at least \code{tol} days.

#' @export
#'
#' @examples
#'
#' ## This function extracts the annual maximum and keep only the years with
#' ## at least 350 days
#'
#' ExtractAmax(flow~date, flowStJohn, tol = 350)
#'

ExtractAmax <- function(x, ...) UseMethod('ExtractAmax',x)

#' @export
#' @rdname ExtractAmax
ExtractAmax.formula <- function(form, x, tol = 0){

  ## reformat dataset according to formula
  x <- model.frame(form,x)

  if(ncol(x) == 2){
    ## Case of one site
    ans <- ExtractAmax.data.frame(x, tol = tol)


  } else {
    ## case multiple sites

    ## split the site
    xlst <- split(x[,c(1,3)], x[,2])
    site.value <- sapply(split(x[,2], x[,2]), '[',1)

    ## extract all annual maximums
    ans <- lapply(xlst, ExtractAmax, tol = tol)

    ## merge the results in one dataset
    cname <- c(colnames(ans[[1]]), colnames(x)[2])
    rname <- unlist(lapply(ans, rownames))

    for(ii in seq_along(site.value))
      ans[[ii]] <- data.frame(ans[[ii]], site.value[ii])

    ans <- do.call('rbind', ans)
    colnames(ans) <- cname
    rownames(ans) <- rname
  }

  ans
}

#' @export
ExtractAmax.data.frame <- function(x, tol = 0){

	## Split data by years
	yy <- format(x[,2],'%Y')
	x <- cbind(x,seq(nrow(x)))
	lx <- split(x,yy)

	## Identify the annual maxima
	nx <- sapply(lx, nrow)
	mx <- sapply(lx, function(z) z[which.max(z[,1]),3])

	ans <- x[mx,-3]
	ans <- cbind(ans, n = nx, yy = as.numeric(format(ans[,2],'%Y')))

	ans <- ans[ans$n>=tol,]

	ans
}

##################################################################
#' Seasonal statistics for flood peaks
#'
#'   Return the circular (Seasonal)statistics for the occurences of flood peaks.
#' The angle represent the average timing of the floods and the radius its
#' regularity. For instance a radius of one represent perfect regularity.
#'
#' @param x Data in Date format
#'
#' @param form Formula
#'
#' @section References:
#'
#' Burn, D.H. (1997). Catchment similarity for regional flood frequency analysis
#'   using seasonality measures. Journal of Hydrology 202, 212–230.
#'   https://doi.org/10.1016/S0022-1694(97)00068-1
#'
#' @export
#'
#' @examples
#'
#' dt <- ExtractAmax(flow~date, flowStJohn)$date
#'
#' st <- SeasonStat(dt)
#' print(st)
#'
#' JulianPlot()
#' points(st['x'], st['y'], col = 2, pch = 16)
#'
#'
SeasonStat <- function(x, ...) UseMethod('SeasonStat', x)

#' @export
SeasonStat.default <- function(x){

  as.Date(x)

  deg <- 2 * pi * DecimalDay(x)

  ## compute average
  xbar <- mean(cos(deg))
  ybar <- mean(sin(deg))

  ## compute the circular statistic
  cs <- Xy2polar(xbar, ybar)

  ans <- c(xbar, ybar, cs$angle, cs$radius)
  names(ans) <- c('x','y','angle','radius')
  ans
}

#' @export
SeasonStat.data.frame <- function(x){
   ans <- lapply(split(x[,1], x[,2]), SeasonStat)
   do.call('rbind', ans)
}

#' @export
SeasonStat.formula <- function(form, x){
  x <- as.data.frame(x)
  x <- model.frame(form,x)
  SeasonStat(x)
}

## Convert the day of the year into a decimal value
## take leap year into account
DecimalDay <- function(x){

  ## Extract julian day
  x <- as.Date(x)
  yy <- as.numeric(format(x, "%Y"))
	dd <- as.numeric(format(x, "%j"))

	## verify for leap years
	isLeap <- is.leapyear(yy)

	## Convert in decimal
	(dd - .5) / (365 + isLeap)

}

## Logical is y a leap year
is.leapyear <- function(y)
  ((y %% 4 == 0) & (y %% 100 != 0)) | (y %% 400 == 0)


## Convert cartesian to polar coordinates
Xy2polar <- function(x,y){

  ## compute polar coordinates
 	ang <- atan(abs(y/x))
 	r <- sqrt(x^2 + y^2)

 	## correcting angle for
  if(sign(x) < 0 & sign(y) >= 0) # second quadrant
    ang <- pi - ang
 	else if(sign(x) < 0 & sign(y) < 0) # third quadrant
 	  ang <- pi + ang
 	else if(sign(x) >= 0 & sign(y) < 0) # fourth quadrant
 	 ang <- 2*pi - ang

 	list(radius = r, angle = ang)
}

##############################################################################
#' Distance in seasonal space
#'
#' Return a square matrix of the distance between points in the seasonal space.
#' It is equivalent Euclidean distance applied to regularity (radius)
#' and timing (angle) as cartesian coordinates.
#'
#' @param r,a Coordinates in the seasonal space: radius `r` and angle `a`.
#'
#' @param form,x Formula and dataset providing the coordinates of the
#'   seasonal space. Must be of the form `radius ~ angle`.
#'
#' @param w Scale factor to standardize difference between angles.
#'
#' @section References:
#'
#' Durocher, M., Burn, D. H., & Ashkar, F. (2019). Comparison of estimation
#'   methods for a nonstationary index-flood model in flood frequency
#'   analysis using peaks over threshold. https://doi.org/10.31223/osf.io/rnepc
#'
#' @export
#'
#' @examples
#'
#' scoord <- data.frame(radius = runif(5), angle = runif(5,0,2*pi))
#' DistSeason(radius ~ angle , scoord)
#'
#'
DistSeason <- function(x,...)  UseMethod('DistSeason', x)

#' @export
DistSeason.matrix <- function(x, w = pi)
  DistSeason(x[,1], x[,2], w)

#' @export
DistSeason.data.frame <- function(x, w = pi)
  DistSeason(x[,1], x[,2], w)

#' @export
DistSeason.formula <- function(form, x, w = pi){
  x <- as.data.frame(x)
  x <- model.frame(form, x) ## form = r ~ a
  DistSeason(x[,1], x[,2], w)
}

#' @export
DistSeason.numeric <- function(r, a, w = pi){

  ## Extract the pairs or every angles
  n <- length(a)

  if(length(r) != n)
    stop('Coordinates must be of the same length')

  id <- expand.grid(1:n, 1:n)
  aii <- a[id[,1]]
  ajj <- a[id[,2]]

  ## Compute the standardized absolute differences between angles
  mn <- pmin(aii,ajj)
  d <- pmax(aii-mn, ajj-mn)
  d <- pmin(2*pi-d, d)/w
  a.mat <- matrix(d, nrow = n)

  ## Compute the absolute differences between radius
  r.mat <- as.matrix(dist(r, method = 'man'))

  ## squared distances
  sqrt(r.mat^2 + a.mat^2)
}

#' Scatter Plot on the unit polar space
#'
#' Plot circular statistics for seasonal data insde the unitary circle.
#'
#' @param x,y Sample to plot, could be a matrix or a list. If \code{x} is a
#'   formula then \code{y} must be a data.frame.
#'
#' @param rose.col,rose.lwd,rose.cex Property of the polar axes.
#' @param ... Other parameter passed to \code{\link{points}}
#'
#' @export
#'
#' @examples
#'
#' JulianPlot()
#'
JulianPlot <- function(rose.col = "gray40", rose.lwd = 1.5,
                       rose.cex = 1.5, rose.radius = seq(.25,1,.25), ...){

  plot(1, pch = '', ylim = c(-1,1)*1.2, xlim = c(-1,1)*1.2, axes = FALSE,
       ylab = "", xlab = "")

  drawCircle(0,0, radius = rose.radius, col = rose.col, lwd = rose.lwd)

  textCircle(month.abb, radius = 1.1, col = rose.col, cex = rose.cex)

  segments(0,-1,0,1, col = rose.col, lwd = rose.lwd)
  segments(-1,0,1,0, col = rose.col, lwd = rose.lwd)

}


drawCircle <- function(x = 0, y = NULL, radius = 1, res = 500, ...){

  if(class(x) %in% c('data.frame','matrix')){
    y <- x[,2]
    x <- x[,1]
  } else if(class(x) == 'list'){
    y <- x$y
    x <- x$x
  } else if(class(x) == 'formula'){
    xd <- model.frame(x,y)
    y <- x[,2]
    x <- x[,1]
  }

  if(is.null(y)) stop('Locations not correctly specified')

  ## Series of angles
  tt <- c(seq(0,2*pi, len = res),0)

  if(length(x) == length(y) & length(y) == length(radius)){
    for(ii in seq_along(radius))
        lines(radius[ii]*cos(tt) + x[ii],
              radius[ii]*sin(tt) + y[ii], type = 'l' ,...)

  } else {
    for(ii in seq_along(radius))
      lines(radius[ii]*cos(tt) + x[1],
            radius[ii]*sin(tt) + y[1], type = 'l' ,...)
  }
}


textCircle <- function(label, x = 0, y = 0, radius = 1, ...){

  pp <- 1/(length(label))

  for(ii in seq_along(label)){
    ang <- 2*pi*pp*(ii-1)
    xii <- radius*cos(ang) + x
    yii <- radius*sin(ang) + y
    text(xii, yii,labels = label[ii],...)
  }
}