R/season.R

In martindurocher/floodRFA: A R-package for regional frequency of floods

#' 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 mlabel Labels to print. By default they are month abbreviation.
#' @param ... Other parameter passed to \code{\link{points}}
#'
#' @export
#'
#' @examples
#'
#' ang <- runif(20, 0, 2*pi)
#' rr <- runif(20,0,1)
#' julianPlot(polar2xy(ang,rr))
#'
julianPlot <- function(x, y = NULL, rose.col = "gray40", rose.lwd = 1.5,
                       rose.cex = 1.5, rose.radius = seq(.25,1,.25),
                       mlabel = NULL, ...){

  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]
  }

  mlabel <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul",
              "Aug", "Sep", "Oct", "Nov", "Dec")


  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(mlabel, 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)

  points(x,y, ...)
}

#' @export
polar2xy <- function(ang,radius){
   x <- cos(ang) * radius
   y <- sin(ang) * radius
   return(list(x = x, y = y))
}

#' @export
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' ,...)
  }
}

#' @export
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],...)
  }
}


##################################################################
#' Seasonal statistics for peaks
#'
#'   Return the seasonal or circular statistics fo the occurences of peaks.
#'
#' @param x Dates
#'
#' @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
#'
#' xx <- canadaFlood$daily
#'
#' xa <- ExtractAmax(flow~date, xx)
#'
#' st <- SeasonStat(xa$date)
#' print(st)
#'
#' ## this indicate that the annual maximum is regular due to the
#' ## spring snowmelt and occur generally at the end of April.
#'
#' julianPlot(NULL)
#' points(st[1], st[2], col = 2, pch = 16)
#'
#'
SeasonStat <- function(x){

  deg <- 2 * pi * DecimalDay(x)

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

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

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


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

######################################################################
#' Convert the day of the year in decimal value
#'
#'  Return the decimal value of a day of the years.
#'  Take the leap year into account.
#'  The formula is \code{(i-0.5)/n} where i is the day of the year and
#'  n is the total number of days (365 or 366).
#'
#' @param x
#'
#' @export
#'
#' @examples
#'
#' DecimalDay(as.Date("2000-1-1")) ## leap year
#' DecimalDay(as.Date("2001-1-1"))
#' x <- as.Date("2001-5-13")
#' DecimalDay(x)
#' DecimalDay(x) * 365
#' format(x,'%j')
#'
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)

}