MacroPotAnalysis

From mn/geo/geoit
Revision as of 10:54, 9 February 2015 by Annefou@uio.no (talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
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)
}