MacroPotAnalysis
From mn/geo/geoit
select_POT <- function(data, threshold, run){ pot_x <- data[data[,2] > threshold, ] pot_x_decl <- declusterGEO4350(pot_x, run) pot_x_decl } declusterGEO4350 <- function (data, spacing){ # data: data frame with 2 columns; time and runoff # threshold: threshold above which events are selected # spacing: number of days that should be between to extreme events >= 0 olddata <- data # transform date values olddata[,1] <- as.POSIXct(strptime(olddata[,1], format="%Y%m%d")) # First run before loop is started to initializ values max <- max(olddata[,2]) # maximum of the runoff values posmax <- match(max, olddata[,2]) # position of the maximum newdata <- olddata[posmax,] olddata <- olddata[-posmax, ] # Loop # The maximim of olddata is taken out and checked, if its date is far enough from the dates of newdata. # If so, the value is inserted, otherwise it is neglected for (i in 1:length(olddata[,1])){ max <- max(olddata[,2]) # maximum of the runoff values posmax <- match(max, olddata[,2]) # position of the maximum datemax <- olddata[posmax,1] # date of the maximum runoff value # check, that the date of the max is not within the spacing distance # of data points in newdata; otherwise it is not further considered date_distance <- abs(as.numeric(difftime(datemax, newdata[,1], units="days"))) if (all(date_distance > spacing)) { # determine insertion position in the new data series n <- length(newdata[,1]) # current length of the new, filtered data series newdata # define bolean vector 'laterbol' with # FALSE for dates of 'newdata' that are smaller than the date of the new maximum and # TRUE for dates of 'newdata' that are bigger than the date of the new maximum laterbol <- datemax < newdata[,1] if (any(laterbol)){ insertposition <- match(1, as.numeric(laterbol)) - 1 } else{ insertposition <- n } # insert the max data point into the new data series if (insertposition == 0){ newdata <- rbind(olddata[posmax,], newdata) } if (insertposition == n){ newdata <- rbind(newdata, olddata[posmax,]) } if (insertposition != 0 & insertposition != n){ newdata <- rbind(newdata[1:insertposition,], olddata[posmax,], newdata[(insertposition+1):n,]) } } # remove the max data point from olddata olddata <- olddata[-posmax, ] } newdata[,1] <- as.integer(format(newdata[,1], format="%Y%m%d")) cat("Data reduced from",length(data[,1]),"to", length(newdata[,1]), "values\n") newdata } qq.gpd <- function(data, xi, beta, threshold){ empirical <- sort(data) n <- length(data) plottingpos <- (1:n)/(n+1) theoretical <- qgpd(plottingpos, xi, threshold, beta) min <- 0.99*min(empirical,theoretical) max <- 1.01*max(empirical, theoretical) plot(empirical, theoretical, xlim=range(min,max), ylim=range(min,max), xlab="Empirical quantiles", ylab="Theoretical quantiles", main="Quantile plot") abline(0,1) } pp.gpd <- function(data, xi, beta, threshold){ n <- length(data) plottingpos <- (1:n)/(n+1) empirical <- plottingpos theoretical <- pgpd(sort(data), xi, threshold, beta) min <- 0 max <- 1 plot(empirical, theoretical, xlim=range(min,max), ylim=range(min,max), xlab="Empirical probability", ylab="Theoretical probability", main="Probability plot: non-exceedance probabilites") abline(0,1) } gevpot.q <- function(k,a,u,l,p){ if(k != 0){ u + (a/k) *(1-(-log(1-p)/l)^k) } else u - a *log(-log(1-p)/l) }