Skip to content

Commit

Permalink
Use line_length_linter(90)
Browse files Browse the repository at this point in the history
  • Loading branch information
MartinRoth committed Jun 21, 2017
1 parent 5c665b0 commit 66a1554
Show file tree
Hide file tree
Showing 12 changed files with 217 additions and 119 deletions.
5 changes: 3 additions & 2 deletions R/Helpers.R
Expand Up @@ -20,11 +20,12 @@ ReadInput <- function(var, ifile) {
}

flog.info("Reading reference data, file={%s}", ifile)
H.comments <- scan(ifile, character(0), sep = "\n", quiet = TRUE) # select lines with "#" from reference file
# select lines with "#" from reference file
H.comments <- scan(ifile, character(0), sep = "\n", quiet = TRUE)
flog.debug("Scanning of the reference data returned n={%i} lines.", length(H.comments))
H.comments <- H.comments[grep("#", H.comments)] # (only necessary for output file)

obs <- read.table(ifile, header = F) # read reference data (header wordt niet apart ingelezen)
obs <- read.table(ifile, header = F) # read reference data
header <- obs[which(obs[, 1] == 0), ] # header met stations meta-data etc.
header[, 1] <- paste0(rep(0, 8), collapse = "") # "00000000"

Expand Down
4 changes: 2 additions & 2 deletions R/MakkinkEvaporation.R
Expand Up @@ -4,8 +4,8 @@ makkink <- function(Tg, Q) {
rho <- 1000 # water mass density [kg/m3]
vps <- 6.107 * 10^((7.5*Tg)/(237.3+Tg)) # saturated vapor pressure [hPa] #nolint
delta <- ((7.5*237.3)/(237.3+Tg)^2) * log(10) * vps # vps gradient [hPa/K] #nolint
gamma <- PsychrometricConstant(Tg) #0.646 + 0.0006*Tg # psychrometric constant [hPa/K]
lambda <- WaterVaporizationEnthalpy(Tg) #1000 * (2501-2.38*Tg) # Enthalpy of vaporization [J/kg]
gamma <- PsychrometricConstant(Tg) # psychrometric constant [hPa/K]
lambda <- WaterVaporizationEnthalpy(Tg) # Enthalpy of vaporization [J/kg]
evmk <- 1000*Q * 1000*0.65*delta / ((delta+gamma)*rho*lambda) # [mm/day] #nolint
return(evmk)
}
Expand Down
51 changes: 33 additions & 18 deletions R/PrecipitationTransformation.R
Expand Up @@ -148,28 +148,36 @@ WetDryDays <- function(fut, wdf, th, mm) {
if (dwet > 0) {

# select target values
step <- length(Xw) / dwet # step size to step through sorted step
target.values <- Xw[round(((1:dwet) - 0.5) * step)] # determine target.values for month <im> #nolint
# 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
preceding.wet <- cumsum(Xm >= th) + step / 2 # cumulative number of preceding wet days in month <im>
add <- vector() # vector with days that should be wetted
# 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 & # select 'first' 'dry' day that succeeds a wet' day,
X1m >= th & # for which <preceding.wet> exceeds the <step> size
preceding.wet >= step)[1]) # and add this day(id) to vector <add>
which(Xm < th &
X1m >= th &
preceding.wet >= step)[1])
if (is.na(add[id])) {
add <- add[-id]
} else {
preceding.wet <- preceding.wet - step # and decrease vector <preceding.wet> with <step>
# 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
# 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")]

} # dfwet > 0
Expand All @@ -184,9 +192,12 @@ CalculateClimatology <- function(obs, deltas, mm, th) {

flog.debug("Calculate climatology")

# qq1 <- 0.99 # quantile of wet-day amounts that is used to estimate transformation coefficients
# qq2 <- 0.90 # quantile of wet-day amounts that is used to estimate qq1 (robustly)
# national median of monthly ratios between qq1 and qq2 for 240 precipitation stations
# qq1 <- 0.99 # quantile of wet-day amounts that is used to estimate
# transformation coefficients
# qq2 <- 0.90 # quantile of wet-day amounts that is used to estimate qq1
# (robustly)
# national median of monthly ratios between qq1 and qq2 for 240
# precipitation stations
ratio <- c(2.22,
2.271,
2.366,
Expand All @@ -205,10 +216,12 @@ CalculateClimatology <- function(obs, deltas, mm, th) {
# 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]
#mwet.obs <- as.matrix(aggregate(obs, by=list(mm), function(x) mean(x[x>=th] )))[, -1, drop = FALSE]
q2.obs <- as.matrix(aggregate(obs, by=list(mm), function(x) quantile(x[x>=th], 0.90)))[, -1, drop = FALSE]
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
Expand Down Expand Up @@ -239,7 +252,8 @@ TransformWetDayAmounts <- function(fut, climatology, mm, th) {
Y <- fut[, is]

for (im in 1:12) {
wet.im <- which(im == mm & Y >= th) # identify all wet days within calendar month <im>
# 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
Expand All @@ -257,7 +271,8 @@ TransformWetDayAmounts <- function(fut, climatology, mm, th) {
# actual transformation of wet-day amounts (application of transformation function)
Y[wet.im] <- ifelse(Xm < qobs, a * Xm^b, c*Xm)

Y[wet.im][which(Y[wet.im] < th)] <- th # prevent days being dried by the wet-day transformation
# 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
Expand Down
6 changes: 3 additions & 3 deletions R/RadiationTransformation.R
Expand Up @@ -38,14 +38,14 @@ rsds_trans_KNMI14 <- function(obs,
for (is in 1:ns) {
for (im in 1:12) {

days.im <- which(mm == im) # all days within in calendar month <im>
X <- obs[days.im, is+1] # select all obs in month <im> of station <is>
days.im <- which(mm == im) # all days within in calendar month <im>
X <- obs[days.im, is+1] # select all obs in month <im> of station <is>

# clear sky radiation for all days in month <im>
Rx <- 0.75 * ObtainAngotRadiation(obs[days.im, 1], lat[is])

# determine coefficient a for transformation function
delta <- deltas[im, 2]/100 # relative change of average in month <im>
delta <- deltas[im, 2]/100 # relative change of average in month <im>
if (delta > 0) {
a <- BoundedScaling(X, Rx, delta)
} else {
Expand Down
28 changes: 18 additions & 10 deletions R/TemperatureTransformation.R
Expand Up @@ -44,8 +44,10 @@ tm_trans_KNMI14 <- function(obs,
fut[, -1] <- NA # future values (filled with NA)

# region information
if (!("regio" %in% colnames(deltas))) deltas$regio <- "NLD" # add column <regio> to deltas if not provided
deltas <- deltas[c("regio", "maand", "P01", "P05", "P50", "P95", "P99")] # arrange deltas
# add column <regio> to deltas if not provided
if (!("regio" %in% colnames(deltas))) deltas$regio <- "NLD"
# arrange deltas
deltas <- deltas[c("regio", "maand", "P01", "P05", "P50", "P95", "P99")]

# TRANSFORMATION
# apply transformation per station <is> / time series
Expand All @@ -57,33 +59,39 @@ tm_trans_KNMI14 <- function(obs,
percentile.changes <- as.numeric(
deltas[which(deltas[, 1]==regio & deltas[, 2]==im), -1:-2] )

days.im <- which(mm == im) # all days within in calendar month <im>
X <- obs[days.im, is+1] # select all obs in month <im> of station <is>
days.im <- which(mm == im) # all days within in calendar month <im>
X <- obs[days.im, is+1] # select all obs in month <im> of station <is>
Y <- rep(NA, length(X))
X.percentiles <- as.numeric(quantile(X, c(1, 5, 50, 95, 99)/100, na.rm=T)) # observed percentiles
Y.percentiles <- X.percentiles + percentile.changes # estimation of future percentiles
# observed percentiles
X.percentiles <- as.numeric(quantile(X, c(1, 5, 50, 95, 99)/100, na.rm=T))
# estimation of future percentiles
Y.percentiles <- X.percentiles + percentile.changes

# linear transformation: for intervals X<qq5, qq5<X<qq50, qq50<X<qq95, qq95<X
# ip = percentile id

ip <- 2 # X < X.percentile[2] (5th percentile)
x.ids <- which(X < X.percentiles[ip]) # id's for all values smaller than second smallest
# id's for all values smaller than second smallest
x.ids <- which(X < X.percentiles[ip])
# estimate linear function for percentile range
a <- (Y.percentiles[ip] - Y.percentiles[ip-1]) /
(X.percentiles[ip] - X.percentiles[ip-1])
b <- Y.percentiles[ip] - a * X.percentiles[ip]
Y[x.ids] <- a * X[x.ids] + b # apply function to all temperatures below 5th percentile
# apply function to all temperatures below 5th percentile
Y[x.ids] <- a * X[x.ids] + b

for (ip in 3:(length(X.percentiles)-1)) {
x.ids <- which(X >= X.percentiles[ip-1] & X < X.percentiles[ip]) # id's all values within analysed percentile range
# id's all values within analysed percentile range
x.ids <- which(X >= X.percentiles[ip-1] & X < X.percentiles[ip])
a <- (Y.percentiles[ip] - Y.percentiles[ip-1]) /
(X.percentiles[ip] - X.percentiles[ip-1])
b <- Y.percentiles[ip] - a * X.percentiles[ip]
Y[x.ids] <- a * X[x.ids] + b
}

ip <- length(X.percentiles)
x.ids <- which(X >= X.percentiles[ip-1]) # id's for all values larger than second largest
# id's for all values larger than second largest
x.ids <- which(X >= X.percentiles[ip-1])
a <- (Y.percentiles[ip] - Y.percentiles[ip-1]) /
(X.percentiles[ip] - X.percentiles[ip-1])
b <- Y.percentiles[ip] - a * X.percentiles[ip]
Expand Down
2 changes: 1 addition & 1 deletion inst/.lintr
@@ -1,6 +1,6 @@
linters: with_defaults(
infix_spaces_linter = NULL, # 569
line_length_linter = NULL, # 189
line_length_linter(90), # 189
snake_case_linter = NULL, # should be incorporated
camel_case_linter = NULL, # 60
commented_code_linter = NULL, # 52
Expand Down
3 changes: 2 additions & 1 deletion tests/testthat/test-evaporation.R
Expand Up @@ -54,7 +54,8 @@ test_that("2030 decadal prediction", {
regions = regions,
rounding = TRUE)

expect_equal_to_reference(tmp, "regressionOutput/evaporation/KNMI14___2030_evmk_rounded.rds")
expect_equal_to_reference(tmp,
"regressionOutput/evaporation/KNMI14___2030_evmk_rounded.rds")
})

test_that("Scenario WL", {
Expand Down
38 changes: 24 additions & 14 deletions tests/testthat/test-precipitation_centr.R
Expand Up @@ -9,31 +9,33 @@ context("rr transformation (centr) - Entire station set")

subscenario <- "centr"
input <- ReadInput("rr", system.file("refdata",
"KNMI14____ref_rrcentr___19810101-20101231_v3.2.txt",
package="knmitransformer"))
"KNMI14____ref_rrcentr___19810101-20101231_v3.2.txt",
package="knmitransformer"))
ofile <- NA
rounding <- FALSE

test_that("2030 decadal prediction", {
scenario <- "GL"
horizon <- 2030
tmp <- TransformPrecip(input=system.file("refdata",
"KNMI14____ref_rrcentr___19810101-20101231_v3.2.txt",
package="knmitransformer"),
"KNMI14____ref_rrcentr___19810101-20101231_v3.2.txt",
package="knmitransformer"),
ofile="tmp.txt",
scenario=scenario,
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14___2030_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14___2030_rr_centr.rds")

tmp <- TransformPrecip(input=input,
ofile=ofile,
scenario=scenario,
horizon = horizon,
subscenario=subscenario,
rounding = TRUE)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14___2030_rr_centr_rounded.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14___2030_rr_centr_rounded.rds")
})

test_that("Scenario GL", {
Expand All @@ -45,7 +47,8 @@ test_that("Scenario GL", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_GL_2050_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_GL_2050_rr_centr.rds")

horizon <- 2085
tmp <- TransformPrecip(input=input,
Expand All @@ -54,7 +57,8 @@ test_that("Scenario GL", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_GL_2085_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_GL_2085_rr_centr.rds")
})

test_that("Scenario GH", {
Expand All @@ -66,7 +70,8 @@ test_that("Scenario GH", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_GH_2050_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_GH_2050_rr_centr.rds")

horizon <- 2085
tmp <- TransformPrecip(input=input,
Expand All @@ -75,7 +80,8 @@ test_that("Scenario GH", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_GH_2085_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_GH_2085_rr_centr.rds")
})

test_that("Scenario WH", {
Expand All @@ -87,7 +93,8 @@ test_that("Scenario WH", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_WH_2050_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_WH_2050_rr_centr.rds")

horizon <- 2085
tmp <- TransformPrecip(input=input,
Expand All @@ -96,7 +103,8 @@ test_that("Scenario WH", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_WH_2085_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_WH_2085_rr_centr.rds")
})

test_that("Scenario WL", {
Expand All @@ -108,7 +116,8 @@ test_that("Scenario WL", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_WL_2050_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_WL_2050_rr_centr.rds")

horizon <- 2085
tmp <- TransformPrecip(input=input,
Expand All @@ -117,5 +126,6 @@ test_that("Scenario WL", {
horizon = horizon,
subscenario=subscenario,
rounding = rounding)
expect_equal_to_reference(tmp, "regressionOutput/precipitation/KNMI14_WL_2085_rr_centr.rds")
expect_equal_to_reference(tmp,
"regressionOutput/precipitation/KNMI14_WL_2085_rr_centr.rds")
})

0 comments on commit 66a1554

Please sign in to comment.