#' rr_trans_KNMI14.R
#' @description Function 'transforms' a specific reference-dataset with daily
#' precipitation sums (mm) to a dataset representative for a future climate
#' scenario.
#' @param obs data.frame or matrix: \cr
#' first column provides datestring YYYYMMDD \cr
#' other columns provide precipitation (mm) time series
#' (each column represents specific station)
#' @param deltas data.frame or matrix that contains deltas (=change factors
#' for the transformation) should contain following columns
#' indicated by following headers
#' HEADER\cr
#' "wdf" relative change (\%) in wet-day frequency
#' (wet day is defined as day with 0.1 mm or more
#' precipitation) \cr
#' "ave" relative change (\%) in mean precipitation \cr
#' "P99" relative change (\%) in the 99th percentile of
#' wet-day amounts
#' @keywords internal
rr_trans_KNMI14 <- function(obs, deltas) {
flog.debug("Running rr_trans_KNMI14")
# DEFINE CONSTANTS
th <- 0.1 # wet-day threshold
# PREPARE DATA
# explore observations
mm <- ObtainMonth(obs[, 1]) # the month of a day (1, 2, ..., 12)
climatology <- CalculateClimatology(obs[, -1, drop = FALSE], deltas, mm, th)
# future values (filled with NA)
fut <- obs
# TRANSFORMATION
# apply transformation per station / time series
fut[, -1] <- DryWetDays(obs, deltas$wdf, th, mm)
fut[, -1] <- WetDryDays(fut[, -1, drop = FALSE], deltas$wdf, th, mm)
fut[, -1] <- TransformWetDayAmounts(fut[, -1, drop = FALSE], climatology, mm, th)
return(fut)
}
DryWetDays <- function(obs, wdf, th, mm) {
# DRYING WET DAYS ##########################################################
if (sum(wdf < 0) > 0) {
# check if reduction in wet days is needed
flog.debug("Drying wet days")
# add very small number (based on datestring) to ensure that all numbers in
# time series are unique.
# necessary for the selection of 'days to dry'
makeUnique <- obs[, 1] * 1e-10
for (is in 2:ncol(obs)) {
Y <- obs[, is]
# select target values
target.values <- vector()
target.months <- vector()
# make Y unique
X <- ifelse(Y < th, Y, Y + makeUnique)
# loop all months for which a reduction of the wet day is projected
for (im in which(wdf < 0)) {
# sorted vector of all wet day amounts that fall in month <im>
Xw <- sort(X[which(X >= th & mm == im)])
# number of days 'to dry'
ndry <- round(-wdf[im] / 100 * length(Xw))
if (ndry > 0) {
# step size to step through wet day amount vector <Xw> (NOT AN INTEGER)
step <- length(Xw) / ndry
# determine target values for month <im> (homogeneously selected from <Xw>)
# and remember specific month <im> that belongs to target.values
target.values <- c(target.values, Xw[round(((1:ndry) - 0.5) * step)]) #nolint
target.months <- c(target.months, rep(im, ndry))
}
}
# assign rank order to all target.values (and belonging month-id)]
target.months <- target.months[order(target.values)]
target.values <- target.values[order(target.values)]
# selection of days to dry
toDry <- SelectDaysToDry(mm, target.values, target.months, X, th)
Y[toDry] <- 0 # actual drying of original time series
obs[, is] <- Y
} # END DRYING WET DAYS
}
return(obs[, -1])
}
SelectDaysToDry <- function(mm, target.values, target.months, X, th) {
nr <- length(mm)
# selection of days to dry
droogmaken <- vector() # vector containing the 'days to dry'
# step through all target values from small to large
for (idry in 1:length(target.values)) {
# select all days that are currently available for drying
# (during the drying procedure new wet days may become available for drying)
# daysInTargetMonth <-
available <- which(mm == target.months[idry] & # all days within same motarget value #nolint
X >= th & # all wet days
(c(0, X[-nr]) < th | # all days preceeded and/or succeeded by dry day #nolint
c(X[-1], 0) < th))
# determine which of all available days is closest
# to the target.value zit en put day(id) in vector
# containing days to dry
droogmaken <- c(droogmaken, available[
which(abs(X[available] - target.values[idry]) ==
min(abs(X[available] - target.values[idry])))])
X[droogmaken[idry]] <- 0 # dry specific day in vector of adjusted values
}
droogmaken
}
WetDryDays <- function(fut, wdf, th, mm) {
flog.debug("Wetting dry days")
for (is in 1:ncol(fut)) {
# WETTING DRY DAYS #######################
Y <- fut[, is]
X <- Y # time series after drying
nr <- length(mm)
X1 <- c(1, X[-nr]) # precipitation of preceding day (preceding day of first day is
# assigned 1 mm)
for (im in 1:12) {
# loop through 12 calendar months
if (wdf[im] > 0) {
# in case an increase of wdf is projected
rows <- which(mm == im) # identify all days in month <im>
Xm <- X[rows] # subset all days in month <im>
X1m <- X1[rows] # and all preceding days
Xw <- sort(Xm[which(Xm >= th)]) # sort all wet day values
dwet <- round((wdf[im] / 100) * length(Xw)) # number of 'dry days to wet' #nolint
if (dwet > 0) {
# select target values
# step size to step through sorted step
step <- length(Xw) / dwet
# determine target.values for month <im>
target.values <- Xw[round( ( (1:dwet) - 0.5) * step)]
# (homogeneously selected from sorted subset)
# select days to wet
# cumulative number of preceding wet days in month <im>
preceding.wet <- cumsum(Xm >= th) + step / 2
add <- vector() # vector with days that should be wetted
for (id in 1:dwet) {
# select 'first' 'dry' day that succeeds a wet' day,
# for which <preceding.wet> exceeds the <step> size
# and add this day(id) to vector <add>
add <- c(add,
which(Xm < th &
X1m >= th &
preceding.wet >= step)[1])
if (is.na(add[id])) {
add <- add[-id]
} else {
# and decrease vector <preceding.wet> with <step>
preceding.wet <- preceding.wet - step
preceding.wet[1:add[id]] <- 0
}
}
# Finally, target.values are assigned to selected days
# on the basis of the rank order of the precipitation amount of the
# preceding wet day
Y[rows[add]] <- target.values[rank(X1m[add], ties.method = "first")]
}
}
}
fut[, is] <- Y
}
return(fut)
}
CalculateClimatology <- function(obs, deltas, mm, th) {
flog.debug("Calculate climatology")
# national median of monthly ratios between wet-day 0.99-quantile and
# 0.90-quantile for 240 precipitation stations (to make qq1 more robust)
ratio <- c(2.22,
2.271,
2.366,
2.147,
2.346,
2.166,
2.276,
2.404,
2.476,
2.087,
2.336,
2.18)
# determine observed:
# wet-day frequency (wdf.obs),
# mean (mean.obs),
# # wet-day mean (mwet.obs),
# wet-day 99th percentile (q1.obs)
wdf.obs <- as.matrix(aggregate(obs, by = list(mm),
function(x) mean( x >= th )))[, -1, drop = FALSE]
mean.obs <- as.matrix(aggregate(obs, by = list(mm),
function(x) mean(x )))[, -1, drop = FALSE]
q2.obs <- as.matrix(aggregate(obs, by = list(mm),
function(x) quantile(x[x >= th], 0.90)))[, -1, drop = FALSE]
q1.obs <- q2.obs * ratio
# apply deltas to observed climatology to obtain future climatology
wdf.fut <- wdf.obs * (1 + deltas$wdf / 100)
mean.fut <- mean.obs * (1 + deltas$ave / 100)
mwet.fut <- mean.fut / wdf.fut
q1.fut <- q1.obs * (1 + deltas$P99 / 100)
list(#mwet.obs = RemoveDimNames(mwet.obs),
mwet.fut = RemoveDimNames(mwet.fut),
qobs = RemoveDimNames(q1.obs),
qfut = RemoveDimNames(q1.fut))
}
RemoveDimNames <- function(x) {
dims <- dim(x)
x <- as.numeric(x)
dim(x) <- dims
x
}
TransformWetDayAmounts <- function(fut, climatology, mm, th) {
flog.debug("Transform wet day amounts")
for (is in 1:ncol(fut)) {
Y <- fut[, is]
for (im in 1:12) {
# identify all wet days within calendar month <im>
wet.im <- which(im == mm & Y >= th)
Xm <- Y[wet.im] # select all wet day amounts
# get climatologies for reference and future period for the month at hand
qobs <- climatology$qobs[im, is]
mfut <- climatology$mwet.fut[im, is]
qfut <- climatology$qfut[im, is]
b <- floor(DeterminePowerLawExponentCpp(Xm, qfut, qobs, mfut) * 1000) / 1000
# straightforward estimation of coefficients a and c
a <- qfut / (qobs^b) #nolint
c <- a * qobs^b / qobs # factor for values larger than q99 # nolint
# actual transformation of wet-day amounts (application of transformation function)
Y[wet.im] <- ifelse(Xm < qobs, a * Xm^b, c * Xm) #nolint
# prevent days being dried by the wet-day transformation
Y[wet.im][which(Y[wet.im] < th)] <- th
}
# END TRANSFORMATION WET-DAY AMOUNTS
fut[, is] <- Y
}
return(fut)
}
# DeterminePowerLawExponent <- function(Xm, qfut, qobs, mfut) {
# # determine exponent 'b' of power-law correction function
#
# # function to minimise to find optimal value for coefficient 'b'
# f <- function(b) {
# qfut / mfut -
# (qobs^b) / mean(ifelse(Xm < qobs, Xm^b, Xm * (qobs^b) / qobs))
# }
#
# # root finding algorithm requires that both sides of search space are
# # of opposite sign
# if(f(0.1) * f(3) < 0) {
# rc <- uniroot(f, lower = 0.1, upper = 3, tol = 0.00001) # root finding
# return(rc$root)
# } else {
# # if root is non-existent, alternative estimation for b
# # value closest to zero is searched for
# bs <- (1:300) / 100 # determine search space for 'b'
# fs <- bs # fs = f(bs)
# for(ifs in 1:length(fs)) {
# fs[ifs] <- f(bs[ifs])
# }
# return(bs[which(abs(fs) == min(abs(fs)))]) # b for which f(b) is smallest is chosen
# }
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.