#' @title Seasonality trend 1960:2010
#' @description Seasonality trend 1960:2010
#' @return Vector with RP, threshold, lm coefficients both for
#' all threshold exceedances and annual peaks of those
#' @author Berry Boessenkool, \email{berry-b@@gmx.de}, Sep 2017
#' @seealso \code{\link{rfsApp}("trend")}
#' @keywords aplot
#' @importFrom extremeStat distLextreme
#' @importFrom berryFunctions almost.equal seasonality between linReg colPointsLegend addAlpha owa
#' @importFrom stats coef lm
#' @importFrom graphics plot.new
#' @export
#' @examples
#' load(seasFolder("data/dismeta.Rdata"))
#' seasTrend("Koeln")
#' seasTrend("Koeln", seasargs=list(xaxt="n", tlab="Empty"))
#' seasTrend("Koeln", legargs=list(y1=0.7, y2=0.95))
#'
#' @param n Character: Name of gauge to be analyzed, see \code{\link{gnames}}
#' @param disdf Dataframe with columns "date" and \code{n}. DEFAULT: dis
#' @param RP Return Period for threshold. DEFAULT: 1.2
#' @param plot Logical: plot result? DEFAULT: TRUE
#' @param map Logical: draw \code{\link{minimap}}? DEFAULT: FALSE
#' @param trex Logical: Plot trend for threshold excesses? DEFAULT: TRUE
#' @param peak Logical: Plot trend for annual peaks of threshold excesses? DEFAULT: TRUE
#' @param startmonth Integer between 1 and 12. To select which months
#' are used in the trend analysis. DEFAULT: 1
#' @param nmonths Integer between 1 and 12. Hom many months are used? DEFAULT: 12
#' @param legpos Position of trend info \code{\link{legend}}. DEFAULT: "left"
#' @param trendargs List of further arguments for trend info. DEFAULT: NULL
#' @param legargs List of arguments for Streamflow legend via
#' \code{\link{colPointsLegend}}. DEFAULT: NULL
#' @param shift Shift passed to \code{\link{seasonality}}. DEFAULT: 61
#' @param seasargs List of arguments passed to \code{\link{seasonality}},
#' like e.g. seasargs=list(xaxt="n"). DEFAULT: NULL
#' @param \dots Further arguments passed to \code{\link{linReg}}
#'
seasTrend <- function(
n,
disdf=get("dis"),
RP=1.2,
plot=TRUE,
map=FALSE,
trex=TRUE,
peak=TRUE,
startmonth=1,
nmonths=12,
legpos="left",
trendargs=NULL,
legargs=NULL,
shift=61,
seasargs=NULL,
...
)
{
# Get threshold from RP:
threshold <- if(RP>=1) thresfuns[[n]](RP) else 0
if(RP>10)
{
plot.new()
text(0.5,0.5,paste0("RP is ", RP, ", but may not be larger than 10."))
return()
}
# Seasonality for all values > threshold:
if(!startmonth %in% 1:12) stop("startmonth must be between 1 and 12, not ", startmonth)
large <- disdf[,n]>=threshold
if(startmonth!=1 | nmonths!=12)
large <- large & between(as.numeric(format(disdf$date, "%m")),
startmonth, startmonth+nmonths-1)
large <- which(large)
meta <- get("meta")
col <- seqPal(100)
if(!trex) col <- NA
seasdef <- list(dates=disdf[large,"date"], values=disdf[large,n], shift=shift,
nmax=1, pch=15, main=paste0(n, ", ", meta[n,"river"]),
returnall=TRUE, col=col, maxargs=list(col=if(peak) "purple" else NA),
adj=0.4, drange=1930:2012, legend=FALSE, plot=ifelse(plot,1,0),
hlines=TRUE)
s <- do.call(seasonality, owa(seasdef, seasargs))
# Trend line:
s_trex <- s$data [between(s$data$year, 1960, 2010),]
s_peak <- s$annmax[between(s$annmax$year, 1960, 2010),]
c_trex <- as.vector(coef(lm(doy~year, data=s_trex)))
c_peak <- as.vector(coef(lm(doy~year, data=s_peak)))
l_trex <- paste0( "All threshold exceedances\nn = ",nrow(s_trex),
", trend slope = ", round(-10*c_trex[2],1), " days/decade\n")
l_peak <- paste0("Annual peaks above threshold\nn = ",nrow(s_peak),
", trend slope = ", round(-10*c_peak[2],1), " days/decade")
if(plot)
{
abline(v=seq(1930, 2010, by=5), col=c(1,8)); box()
if(peak) linReg(doy~year, data=s_peak, add=TRUE, plotrange=1960:2011, col="purple",
colband=addAlpha("purple", alpha=0.08), lwd=2, pos1=NA, ...)
if(trex) linReg(doy~year, data=s_trex, add=TRUE, plotrange=1960:2011, col="orange", pos1=NA, lwd=2, ...)
trenddef <- list(x=legpos,
legend=c(ifelse(trex, l_trex, " \n \n "), ifelse(peak, l_peak, " \n ")),
lty=0, text.col=c("orange","purple"), bg=addAlpha("white"), box.lty=0)
do.call(legend, owa(trenddef, trendargs))
# Legend and map:
legdef <- list(z=disdf[large,n], nlab=4, title=paste0("Streamflow > ",
round(min(disdf[large,n], na.rm=TRUE))," [m\U{00B3}/s]"),
y1=0.89, y2=1, colors=col)
do.call(colPointsLegend, owa(legdef, legargs))
if(map) minimap(n)
}
# Output
out <- c(RP=RP, threshold=threshold, n=length(large),
slope_trex=c_trex[2], slope_peak=c_peak[2],
inter_trex=c_trex[1], inter_peak=c_peak[1])
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.