#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.