library(regthresh) library(copula) library(rgdal) library(extRemes) library(potreg) library(ggmap) library(data.table) library(doParallel) registerDoParallel()
sites <- 16 years <- 50 blockLength <- 92 days <- years * blockLength probs <- seq(.9, .995, by = 0.001) retPeriods <- c(5, 50, 100)
xpar <- list() xpar$S <- sites xpar$T <- days xpar$fpar <- function(p, xpar) { loc <- xpar$loc scale <- p[1] * loc shape <- matrix(p[2], xpar$T, xpar$S) list(loc = loc, scale = scale, shape = shape) } xpar$nOfP <- 2 xpar$start <- c(0.5, 0)
tau <- 0.95 epsilon <- 0.25 gamma <- 0.5 xi <- 0.15 kappasV <- rbeta(xpar$S, 2, 5) + 0.5 kappasP <- rep(calculateCorrespondingWeibullShape(0.15, tau), xpar$S) betas <- runif(16, 2, 4) mParamsV <- setMarginalParameters(tau = tau, epsilon = epsilon, betas = betas ,kappas = kappasV, gamma = gamma, xi = xi, xpar = xpar) mParamsP <- setMarginalParameters(tau = tau, epsilon = epsilon, betas = betas ,kappas = kappasP, gamma = gamma, xi = xi, xpar = xpar)
x <- seq(0, 50, by = 0.1) plot(x, pDW(x, mParamsP[[1]]), type = "l") for(i in 2 : 5) lines(x, pDW(x, mParamsV[[i]]), lty = i) abline(h = 0.95, col = 2, lty = 3) x <- seq(12, 16, by = 0.1) plot(x, dDW(x, mParamsP[[1]]), type = "l", ylim = c(0, 0.025)) lines(x, dDW(x, mParamsV[[1]]), col = 2) for(i in 2 : 5) lines(x, dDW(x, mParamsV[[i]]), lty = i) x <- seq(9, 20, by = 0.1) plot(x, hazDW(x, mParamsP[[1]]), type = "l", ylim = c(0, .3)) lines(x, hazDW(x, mParamsV[[1]]), col = 2) for(i in 2 : 5) lines(x, hazDW(x, mParamsV[[i]]), lty = i)
copulaInd <- normalCopula(0, sites) dataVInd <- generateDWData(mParamsV, copulaInd, xpar) dataPInd <- generateDWData(mParamsP, copulaInd, xpar)
Plot area of Vallei en Veluwe waterboard.
waterschap.rg <- readOGR("./inst/extData/Waterschapsgrenzen_28992/","AU_AdministrativeUnitPolygon") waterschap.sp <- spTransform(waterschap.rg, CRS("+proj=longlat")) ValleiEnVeluwe.sp <- subset(waterschap.sp, waterschap.sp$name == "Vallei & Veluwe") ValleiEnVeluwe.dat <- fortify(ValleiEnVeluwe.sp) basemap <- ggmap(get_map(location = c(lon = 5, lat = 52.5), color = "color", source = "stamen", maptype = "watercolor", zoom = 7)) focusBasemap <- ggmap(get_map(location = c(lon = 5.75, lat = 52.25), color = "color", source = "stamen", maptype = "watercolor", zoom = 9))
Load rain data
prDat <- knmiR::HomogenPrecip(ValleiEnVeluwe.sp, period="1951/2015") prDat[, year := year(date)] prDat[, month := month(date)] prDat <- prDat[month %in% c(6,7,8)]
overview <- basemap + geom_path(aes(x = long, y = lat), col = "red", size = 1, data = ValleiEnVeluwe.dat) print(overview) setkey(prDat, stationId) focusview <- focusBasemap + geom_path(aes(x = long, y = lat), col = "red", size = 1, data = ValleiEnVeluwe.dat) + geom_point(aes(x = lon, y = lat), col = "black", shape = 3, data = knmiR::stationMetaData[.(prDat[, unique(stationId)]), ]) + geom_point(aes(x = lon, y = lat), col = "blue", shape = 4, data = knmiR::stationMetaData[stationName=="Putten", ]) print(focusview)
summary(prDat)
prDat[, jdays := julian(date, origin=as.Date("1951-01-01"))] setkey(prDat, stationId) #prDat[, prDcl2 := NA_real_] #system.time(prDat[, prDcl := potreg::declusterData(jdays, pr, 1), by = stationId]) #system.time({ prDcl <- foreach(stnId=unique(prDat[["stationId"]]), .combine = "c") %dopar% { prDat[.(stnId), declusterData(jdays, pr, 1)] } prDat[, prDcl := prDcl] #})
ObtainTSPlot <- function(x, tau) { tmp <- foreach(q = iter(tau), .combine = "rbind") %dopar% { data.table(tau = q, shape = fevd(x, threshold = quantile(x, q), type="GP")$results$par[2]) } return(tmp) } setQuantiles <- function(dat, tau) { return(matrix(apply(dat, 2, quantile, tau), nrow=nrow(dat), ncol=ncol(dat), byrow=TRUE)) } getExceedances <- function(data, u) { index <- data > u data[!index] <- NA return(data) } getExcesses <- function(data, u) { excesses <- getExceedances(data - u, 0) return(excesses) }
mTS1 <- prDat[stationId == knmiR::stationMetaData[stationName=="Putten", stationId], ObtainTSPlot(prDcl, probs)] #plot(probs[-(93:96)], mTS1[-(93:96)], xlab = "", ylab = "") #ggplot(mTS1, aes(x = tau, y = shape)) + geom_line() + xlim(c(0.9, 0.985)) + ylim(c(0, 0.3))
mSIVInd <- prDat[, ObtainTSPlot(prDcl, probs), by = stationId] mSIVInd <- mSIVInd[, .(type = "individual", shape = mean(shape)), by = tau] #ggplot(mSIVInd, aes(x = tau, y = shape)) + geom_line() + xlim(c(0.9, 0.985)) + ylim(c(0.025, 0.2)) tmpData <- as.matrix(dcast(prDat, date ~ stationId, value.var = 'prDcl')[, -1, with=FALSE]) tmpXpar <- list() tmpXpar$S <- ncol(tmpData) tmpXpar$T <- nrow(tmpData) fpar <- function(p, xpar) { loc = xpar$loc scale = p[1]*loc shape = matrix(p[2], nrow = tmpXpar$T, ncol = tmpXpar$S) list(loc = loc, scale = scale, shape = shape) } mTSCV <- foreach(q = iter(probs), .combine = "rbind") %dopar% { tmpXpar$loc <- setQuantiles(tmpData, q) excess <- getExcesses(tmpData, tmpXpar$loc) tmp <- fitGpdFlex(excess, tmpXpar, fpar, hessian = TRUE, numberOfParameters = 2, start = c(0.5, 0), method = "BFGS", control = list(maxit = 5000, ndeps = c(1e-03, 1e-05))) return(data.table(type = "joint", tau = q, shape = tmp$estimate[2])) } ggplot(rbind(mSIVInd, mTSCV), aes(x = tau, y = shape, col = type)) + geom_line() + xlim(c(0.9, 0.985)) + ylim(c(0.025, 0.2))
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.