Commit ee93972a authored by Jana Ulrich's avatar Jana Ulrich
Browse files

IDF.plot without data points, because it makes more sense to add them individually

parent b6c9c384
......@@ -9,3 +9,4 @@
/*.Rproj
.Rbuildignore
.Rproj.user/
.Rproj.user
......@@ -99,20 +99,15 @@ IDF.agg <- function(data,ds,na.accept = 0,
#' Plotting of IDF curves at a chosen station
#'
#' @param data matrix or dataframe containing: first column maxima,
#' second column coresponding durations
#' @param durations vector of durations for which to calculate the quantiles.
#' @param fitparams vector containing parameters mut, sigma0, xi, theta, eta
#' (modified location, scale, shape, duration offset, duration exponent) for chosen station
#' as obtained from \code{\link{gev.d.fit}}.
#' as obtained from \code{\link{gev.d.fit}}
#' (or \code{\link{gev.d.params}} for model with covariates).
#' @param probs vector of exeedance probabilities for which to plot IDF curves (p = 1-1/ReturnPeriod)
#' @param calc.dur vector of durations for which to calculate IDF curves. If `NULL` (the default),
#' durations from data are taken
#' @param cols vector of colors for IDF curves. Should have same length as \code{probs}
#' @param add logical indicating if plot should be added to existing plot
#' @param xlim,ylim vectors of x- / y-plot-range
#' @param legend logical indicating if legend should be plotted
#' @param st.name string containing legend title (for example station name)
#' @param dt.name string containing name for data points in legend
#' @param ... additional parameters passed on to the \code{plot} function
#'
#' @export
......@@ -120,54 +115,54 @@ IDF.agg <- function(data,ds,na.accept = 0,
#' @importFrom graphics axis box lines plot points
#' @examples
#' data('example',package = 'IDF')
#' # fit d-gev
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
#' ,mul = c(1,2),sigl = 1)
#' # get parameters for cov1 = 1, cov2 = 1
#' par <- gev.d.params(fit = fit, ydat = matrix(1,1,2))
#' IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par),
#' calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example')
IDF.plot <- function(data,fitparams,probs=c(0.5,0.9,0.99),calc.dur=NULL,
cols=4:2,add=FALSE,ylim=NULL,xlim=NULL,
legend=TRUE, st.name='Station',dt.name="observations",...){
#' # plot quantiles
#' IDF.plot(durations = seq(0.5,35,0.2),fitparams = par)
#' # add data points
#' points(example[example$cov1==1,]$d,example[example$cov1==1,]$dat)
IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
cols=4:2,add=FALSE,
legend=TRUE,...){
# if cols is to short, make longer
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
# if no calc.dur given, take durations from data
if(is.null(calc.dur))calc.dur <- data[,2]
## calculate IDF values for given probability and durations
idf.array <- simplify2array(qgev.d(probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3],
theta=fitparams[4],eta=fitparams[5],d=calc.dur)) # array[probs,durs]
theta=fitparams[4],eta=fitparams[5],d=durations)) # array[probs,durs]
if(!add){ #new plot
## initialize plot window
# check if limits where passed
if(is.null(ylim)){ylim <- range(idf.array,data[,1],na.rm=TRUE)}
if(is.null(xlim)){xlim <- range(data[,2],calc.dur,na.rm=TRUE)}
# empty plot
plot(NA,axes=F,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy",...)
box('plot')
axis(1,at=data[,2],labels=round(data[,2],digits = 2))
axis(2)
# check if `pch` was passed in ...
if(is.element('pch',names(list(...))))pch <- list(...)[['pch']] else pch <- 16
## plot data points
points(data[,2],data[,1],pch=pch,col=rgb(0,0,0,0.5))
}else points(data[,2],data[,1],pch=16,col=rgb(0,0,0,0.5)) # add to existing plot
# check if `lty,ltd` were passed in ...
if(is.element('lty',names(list(...))))lty <- list(...)[['lty']] else lty <- 1
if(is.element('lwd',names(list(...))))lwd <- list(...)[['lwd']] else lwd <- 1
# check if limits were passed
if(is.element('ylim',names(list(...)))){
ylim <- list(...)[['ylim']]
}else{ylim <- range(idf.array,na.rm=TRUE)}
if(is.element('ylim',names(list(...)))){
xlim <- list(...)[['xlim']]
}else{xlim <- range(durations,na.rm=TRUE)}
# empty plot
plot(NA,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy")
}
## plot IDF curves
for(i in 1:length(probs))
lines(calc.dur,idf.array[i,],col=cols[i],lty=lty,lwd=lwd)
for(i in 1:length(probs)){
lines(durations,idf.array[i,],col=cols[i],...)
}
if(legend){## plot legend
legend(x="topright",title = st.name,legend=c(dt.name,paste(probs,"quantile",sep=" ")),
col=c(rgb(0,0,0,0.5),cols),lty=c(NA,rep(lty,length(probs))),
pch=c(pch,rep(NA,length(probs))),lwd=c(NA,rep(lwd,length(probs))))}
}
# check if lwd, lty were passed
if(is.element('lwd',names(list(...)))){
lwd <- list(...)[['lwd']]
}else{lwd <- 1}
if(is.element('lty',names(list(...)))){
lty <- list(...)[['lty']]
}else{lty <- 1}
legend(x="topright",title = 'p-quantile',legend=probs,
col=cols,lty=lty,lwd=lwd)
}
}
\ No newline at end of file
......@@ -153,7 +153,8 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
qfun <- function(d){
if(d+theta<=0){return(rep(NA,length(p)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(qgev(p,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
return(qgev(p,loc=as.numeric(mut*sigma.d)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
}
# calculate for each duration in d
durations <- as.list(d)
......
......@@ -4,35 +4,25 @@
\alias{IDF.plot}
\title{Plotting of IDF curves at a chosen station}
\usage{
IDF.plot(data, fitparams, probs = c(0.5, 0.9, 0.99), calc.dur = NULL,
cols = 4:2, add = FALSE, ylim = NULL, xlim = NULL,
legend = TRUE, st.name = "Station", dt.name = "observations", ...)
IDF.plot(durations, fitparams, probs = c(0.5, 0.9, 0.99), cols = 4:2,
add = FALSE, legend = TRUE, ...)
}
\arguments{
\item{data}{matrix or dataframe containing: first column maxima,
second column coresponding durations}
\item{durations}{vector of durations for which to calculate the quantiles.}
\item{fitparams}{vector containing parameters mut, sigma0, xi, theta, eta
(modified location, scale, shape, duration offset, duration exponent) for chosen station
as obtained from \code{\link{gev.d.fit}}.}
as obtained from \code{\link{gev.d.fit}}
(or \code{\link{gev.d.params}} for model with covariates).}
\item{probs}{vector of exeedance probabilities for which to plot IDF curves (p = 1-1/ReturnPeriod)}
\item{calc.dur}{vector of durations for which to calculate IDF curves. If `NULL` (the default),
durations from data are taken}
\item{cols}{vector of colors for IDF curves. Should have same length as \code{probs}}
\item{add}{logical indicating if plot should be added to existing plot}
\item{xlim, ylim}{vectors of x- / y-plot-range}
\item{legend}{logical indicating if legend should be plotted}
\item{st.name}{string containing legend title (for example station name)}
\item{dt.name}{string containing name for data points in legend}
\item{...}{additional parameters passed on to the \code{plot} function}
}
\description{
......@@ -40,9 +30,13 @@ Plotting of IDF curves at a chosen station
}
\examples{
data('example',package = 'IDF')
# fit d-gev
fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
,mul = c(1,2),sigl = 1)
# get parameters for cov1 = 1, cov2 = 1
par <- gev.d.params(fit = fit, ydat = matrix(1,1,2))
IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par),
calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example')
# plot quantiles
IDF.plot(durations = seq(0.5,35,0.2),fitparams = par)
# add data points
points(example[example$cov1==1,]$d,example[example$cov1==1,]$dat)
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment