Commit 5212b6cd authored by Jana Ulrich's avatar Jana Ulrich
Browse files

changed gev.d.diag legend

parent 10d6852a
(requireNamespace('IDF',quietly = TRUE))
?IDF.short
IDF.short(int.vec = test.data, durs = rep(c(1,2,3),each=100), n.y = 100)
install.packages("~/Downloads/IDF-covariates-6c1b09f1f9fc97a0b518b9ae025ba535095f63ea.tar.gz", repos = NULL, type = "source")
library(ismev)
library(IDF)
test.data <- matrix(NA,ncol=3,nrow=100)
for(i in 1:3) test.data[,i] <- rgev.d(n = 100,d=i)
test.data <- as.vector(test.data)
?rgev.d
for(i in 1:3) test.data[,i] <- rgev.d(n = 100,d=i)
?rgev.d
?rgev
?IDF.agg
IDF.agg(test.data,agg.lev=c(1,2,3))
detach("package:IDF", unload=TRUE)
remove.packages("IDF", lib="~/R/x86_64-pc-linux-gnu-library/3.3")
library(IDF)
install.packages("~/Downloads/IDF-covariates-6c1b09f1f9fc97a0b518b9ae025ba535095f63ea.tar.gz", repos = NULL, type = "source")
library(IDF)
IDF.agg(test.data,agg.lev=c(1,2,3))
detach("package:IDF", unload=TRUE)
remove.packages("IDF", lib="~/R/x86_64-pc-linux-gnu-library/3.3")
library(IDF)
install.packages("~/Downloads/IDF-covariates-6c1b09f1f9fc97a0b518b9ae025ba535095f63ea(1).tar.gz", repos = NULL, type = "source")
library(ismev)
library(IDF)
IDF.agg(test.data,agg.lev=c(1,2,3))
?IDF.short
detach("package:IDF", unload=TRUE)
remove.packages("IDF", lib="~/R/x86_64-pc-linux-gnu-library/3.3")
?IDF.short
?install.packages
devtools::install_git('https://gitlab.met.fu-berlin.de/Rpackages/IDF/tree/covariates')
devtools::install_git('https://gitlab.met.fu-berlin.de/Rpackages/IDF')
devtools::install_git('https://gitlab.met.fu-berlin.de/Rpackages/IDF/tree/covariates.git')
.Rproj.user
.Rhistory
.RData
.Ruserdata
.Rbuildignore
.Rproj
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: Sweave
LaTeX: pdfLaTeX
BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
...@@ -25,7 +25,7 @@ importFrom(graphics,plot) ...@@ -25,7 +25,7 @@ importFrom(graphics,plot)
importFrom(graphics,points) importFrom(graphics,points)
importFrom(graphics,title) importFrom(graphics,title)
importFrom(ismev,gev.fit) importFrom(ismev,gev.fit)
importFrom(pbapply,pblapply) importFrom(pbapply,pbsapply)
importFrom(stats,lm) importFrom(stats,lm)
importFrom(stats,optim) importFrom(stats,optim)
importFrom(zoo,rollapply) importFrom(zoo,rollapply)
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
#' standard date format. #' standard date format.
#' @param ds numeric vector of aggregation durations. #' @param ds numeric vector of aggregation durations.
#' (Must be multiples of time resolution at all stations.) #' (Must be multiples of time resolution at all stations.)
#' @param na.accept numeric giving maximum number of missing values for which annual max. should still be calculated
#' @param which.stations optional, subset of stations. Either numeric vector or character vector #' @param which.stations optional, subset of stations. Either numeric vector or character vector
#' containing names of elements in data. If not given, all elements in `data` will be used. #' containing names of elements in data. If not given, all elements in `data` will be used.
#' @param which.mon optional, subset of months of which to calculate the annual maxima from. #' @param which.mon optional, subset of months of which to calculate the annual maxima from.
...@@ -26,13 +27,13 @@ ...@@ -26,13 +27,13 @@
#' different durations, IDF.agg needs to be run seperately for the different groups of stations. #' different durations, IDF.agg needs to be run seperately for the different groups of stations.
#' Afterwards he results can be joint together using `rbind`. #' Afterwards he results can be joint together using `rbind`.
#' #'
#' @return data.frame containing the annual maxima in `$xdat`, the corresponding duration in `$ds` #' @return data.frame containing the annual intensity maxima [mm/h] in `$xdat`, the corresponding duration in `$ds`
#' and the station id or name in `$station`. #' and the station id or name in `$station`.
#' #'
#' @seealso \code{\link{pgev.d}} #' @seealso \code{\link{pgev.d}}
#' #'
#' @export #' @export
#' @importFrom pbapply pblapply #' @importFrom pbapply pbsapply
#' @importFrom zoo rollapply #' @importFrom zoo rollapply
#' #'
#' @examples #' @examples
...@@ -46,7 +47,9 @@ ...@@ -46,7 +47,9 @@
#'## 2 0.4112304 24 1 #'## 2 0.4112304 24 1
#'## 3 0.1650978 48 1 #'## 3 0.1650978 48 1
#'## 4 0.2356849 48 1 #'## 4 0.2356849 48 1
IDF.agg <- function(data,ds,which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){ IDF.agg <- function(data,ds,na.accept = 0,
which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){
if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data)) if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data))
# function 2: aggregate station data over durations and find annual maxima: # function 2: aggregate station data over durations and find annual maxima:
...@@ -65,14 +68,18 @@ IDF.agg <- function(data,ds,which.stations = NULL,which.mon = 0:11,names = c('da ...@@ -65,14 +68,18 @@ IDF.agg <- function(data,ds,which.stations = NULL,which.mon = 0:11,names = c('da
# function 1: aggregate over single durations and find annual maxima: # function 1: aggregate over single durations and find annual maxima:
agg.ts <- function(ds){ agg.ts <- function(ds){
runmean <- rollapply(data.s[,names[2]],ds/dtime,FUN=sum,fill =NA,align='right',na.rm= TRUE) runmean <- rollapply(data.s[,names[2]],ds/dtime,FUN=sum,fill =NA,align='right')
runmean <- runmean/ds #intensity per hour runmean <- runmean/ds #intensity per hour
subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon) subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
max <- tapply(runmean[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],max,na.rm=T) max <- tapply(runmean[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],
function(vec){
n.na <- sum(is.na(vec))
max <- ifelse(n.na <= na.accept,max(vec,na.rm = TRUE),NA)
return(max)})
return(max) # maxima for single durations return(max) # maxima for single durations
} }
# call function 1 in lapply to aggregate over all durations at single station # call function 1 in lapply to aggregate over all durations at single station
data.agg <- sapply(ds,agg.ts,simplify=TRUE) data.agg <- pbsapply(ds,agg.ts,simplify=TRUE,cl=cl)
df <- data.frame(xdat = as.vector(data.agg), ds = rep(ds,each=length(data.agg[,1]))) df <- data.frame(xdat = as.vector(data.agg), ds = rep(ds,each=length(data.agg[,1])))
df$station <- rep(station,length(df$xdat)) df$station <- rep(station,length(df$xdat))
...@@ -81,12 +88,11 @@ IDF.agg <- function(data,ds,which.stations = NULL,which.mon = 0:11,names = c('da ...@@ -81,12 +88,11 @@ IDF.agg <- function(data,ds,which.stations = NULL,which.mon = 0:11,names = c('da
# which stations should be used? # which stations should be used?
if(is.null(which.stations))which.stations <- 1:length(data) if(is.null(which.stations))which.stations <- 1:length(data)
# call function 2 in lapply to aggregate over all durations at all stations # call function 2 in lapply to aggregate over all durations at all stations
station.list <- pblapply(which.stations,agg.station,cl=cl) station.list <- lapply(which.stations,agg.station)
return(do.call('rbind',station.list)) return(do.call('rbind',station.list))
} }
#### IDF.plot #### #### IDF.plot ####
...@@ -163,6 +169,5 @@ IDF.plot <- function(data,fitparams,probs=c(0.5,0.9,0.99),calc.dur=NULL, ...@@ -163,6 +169,5 @@ IDF.plot <- function(data,fitparams,probs=c(0.5,0.9,0.99),calc.dur=NULL,
if(legend){## plot legend if(legend){## plot legend
legend(x="topright",title = st.name,legend=c(dt.name,paste(probs,"quantile",sep=" ")), 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))), 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)))) pch=c(pch,rep(NA,length(probs))),lwd=c(NA,rep(lwd,length(probs))))}
}
} }
...@@ -306,7 +306,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -306,7 +306,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
title=c('Residual Probability Plot','Residual Quantile Plot'), title=c('Residual Probability Plot','Residual Quantile Plot'),
emp.lab='Empirical',mod.lab='Model',...){ emp.lab='Empirical',mod.lab='Model',...){
# check parameter: # check parameter:
if(!is.element(which,c('both','pp','qq'))) stop("Parameter `which`= ",which, if(!is.element(which,c('both','pp','qq'))) stop("Parameter 'which'= ",which,
" but only 'both','pp' or 'qq' are allowed.") " but only 'both','pp' or 'qq' are allowed.")
# subset data # subset data
df <- data.frame(data=fit$data,ds=fit$ds) df <- data.frame(data=fit$data,ds=fit$ds)
...@@ -335,7 +335,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -335,7 +335,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
emp.lab, ylab = mod.lab,col=cols[df$cval],pch=pch,...) emp.lab, ylab = mod.lab,col=cols[df$cval],pch=pch,...)
abline(0, 1, col = 1,lwd=1) abline(0, 1, col = 1,lwd=1)
title(title[1]) title(title[1])
if(legend){legend('bottomright',legend = durs,pch=pch, if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
col = cols[1:length(durs)],title = 'Durations',ncol = 2)} col = cols[1:length(durs)],title = 'Durations',ncol = 2)}
} }
if(which=='both'|which=='qq'){ if(which=='both'|which=='qq'){
...@@ -344,7 +344,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -344,7 +344,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
emp.lab, xlab = mod.lab,col=cols[df$cval],pch=pch,...) emp.lab, xlab = mod.lab,col=cols[df$cval],pch=pch,...)
abline(0, 1, col = 1,lwd=1) abline(0, 1, col = 1,lwd=1)
title(title[2]) title(title[2])
if(legend){legend('bottomright',legend = durs,pch=pch, if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
col = cols[1:length(durs)],title = 'Durations',ncol = 2)} col = cols[1:length(durs)],title = 'Durations',ncol = 2)}
} }
if(which=='both') par(mfrow=c(1,1)) # reset par if(which=='both') par(mfrow=c(1,1)) # reset par
......
...@@ -4,8 +4,8 @@ ...@@ -4,8 +4,8 @@
\alias{IDF.agg} \alias{IDF.agg}
\title{Aggregation and annual maxima for choosen durations} \title{Aggregation and annual maxima for choosen durations}
\usage{ \usage{
IDF.agg(data, ds, which.stations = NULL, which.mon = 0:11, IDF.agg(data, ds, na.accept = 0, which.stations = NULL,
names = c("date", "RR"), cl = NULL) which.mon = 0:11, names = c("date", "RR"), cl = NULL)
} }
\arguments{ \arguments{
\item{data}{list of data.frames containing time series for every station. \item{data}{list of data.frames containing time series for every station.
...@@ -16,6 +16,8 @@ standard date format.} ...@@ -16,6 +16,8 @@ standard date format.}
\item{ds}{numeric vector of aggregation durations. \item{ds}{numeric vector of aggregation durations.
(Must be multiples of time resolution at all stations.)} (Must be multiples of time resolution at all stations.)}
\item{na.accept}{numeric giving maximum number of missing values for which annual max. should still be calculated}
\item{which.stations}{optional, subset of stations. Either numeric vector or character vector \item{which.stations}{optional, subset of stations. Either numeric vector or character vector
containing names of elements in data. If not given, all elements in `data` will be used.} containing names of elements in data. If not given, all elements in `data` will be used.}
...@@ -26,7 +28,7 @@ containing names of elements in data. If not given, all elements in `data` will ...@@ -26,7 +28,7 @@ containing names of elements in data. If not given, all elements in `data` will
\item{cl}{optional, number of cores to be used from \code{\link[pbapply]{pblapply}} for parallelization.} \item{cl}{optional, number of cores to be used from \code{\link[pbapply]{pblapply}} for parallelization.}
} }
\value{ \value{
data.frame containing the annual maxima in `$xdat`, the corresponding duration in `$ds` data.frame containing the annual intensity maxima [mm/h] in `$xdat`, the corresponding duration in `$ds`
and the station id or name in `$station`. and the station id or name in `$station`.
} }
\description{ \description{
......
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