Skip to content

Commit

Permalink
Use linters infix spaces and snake case
Browse files Browse the repository at this point in the history
  • Loading branch information
MartinRoth committed Jun 21, 2017
1 parent 66a1554 commit 7cf953b
Show file tree
Hide file tree
Showing 19 changed files with 533 additions and 500 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
@@ -1,8 +1,8 @@
Package: knmitransformer
Type: Package
Title: Transformation program
Version: 0.2.7
Date: 2017-06-20
Version: 0.2.8
Date: 2017-06-21
Authors@R: as.person(c(
"Alexander Bakker [aut]",
"Martin Roth <roth@knmi.nl> [aut, cre]",
Expand Down
2 changes: 1 addition & 1 deletion R/Helpers.R
Expand Up @@ -303,5 +303,5 @@ CheckRegions <- function(regions, nStations) {
#' @param day integer in format yyyymmdd
#' @keywords internal
ObtainMonth <- function(day) {
(day%/%100)%%100 # the month that a day belongs to (1, 2, ..., 12)
(day %/% 100) %% 100 # the month that a day belongs to (1, 2, ..., 12)
}
45 changes: 20 additions & 25 deletions R/PrecipitationTransformation.R
Expand Up @@ -140,8 +140,8 @@ WetDryDays <- function(fut, wdf, th, mm) {
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>
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
Expand Down Expand Up @@ -180,9 +180,9 @@ WetDryDays <- function(fut, wdf, th, mm) {
# preceding wet day
Y[rows[add]] <- target.values[rank(X1m[add], ties.method = "first")]

} # dfwet > 0
} # days need to be added
} # calander month
}
}
}
fut[, is] <- Y
}
return(fut)
Expand All @@ -192,12 +192,8 @@ 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
# 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,
Expand All @@ -216,19 +212,19 @@ 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]
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
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)
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)
q1.fut <- q1.obs * (1 + deltas$P99 / 100)

list(#mwet.obs = RemoveDimNames(mwet.obs),
mwet.fut = RemoveDimNames(mwet.fut),
Expand Down Expand Up @@ -257,19 +253,18 @@ TransformWetDayAmounts <- function(fut, climatology, mm, th) {
Xm <- Y[wet.im] # select all wet day amounts

# get climatologies for reference and future period for the month at hand
#mobs <- climatology$mwet.obs[im,is]
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)
c <- a*qobs^b / qobs # multiplication factor for values larger than q99
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)
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
Expand Down
48 changes: 27 additions & 21 deletions R/RadiationTransformation.R
Expand Up @@ -39,20 +39,20 @@ rsds_trans_KNMI14 <- function(obs,
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>
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 {
a <- 1 + delta
}
Y <- TransformRad(a, X, Rx)
fut[days.im, is+1] <- Y
fut[days.im, is + 1] <- Y
}
} # END TRANSFORMATION LOOP

Expand All @@ -62,35 +62,41 @@ rsds_trans_KNMI14 <- function(obs,

# X will not exceed Xmax
TransformRad <- function(a, X, Xmax) {
apply(cbind(a*X, Xmax), 1, min)
apply(cbind(a * X, Xmax), 1, min)
}

# iterative estimation of coefficient a in function <tf>
BoundedScaling <- function(X, Xmax, delta) {
f <- function(a) mean(TransformRad(a, X, Xmax)) - (1+delta)*mean(X)
rc <- uniroot(f, lower=0, upper=4, tol = 0.01)
f <- function(a) mean(TransformRad(a, X, Xmax)) - (1 + delta) * mean(X)
rc <- uniroot(f, lower = 0, upper = 4, tol = 0.01)
a <- rc$root
return(a)
}

# "Angot" or "Top Of Atmoshphere" radiation
ObtainAngotRadiation <- function(datestring_YYYYMMDD, lat) {
#' Obtain Top of Atmossphere radiation
#' @inheritParams daynumber
#' @param lat latitude
#' @keywords internal
ObtainAngotRadiation <- function(datestring, lat) {
Gsc <- 0.0820 # solar constant [MJ/m2/min]
phi <- pi*lat/180
J <- daynumber(datestring_YYYYMMDD)
dr <- 1 + 0.033*cos(2*pi*J/365)
delta <- 0.409*sin(2*pi*J/365-1.39)
omega <- acos(-tan(phi)*tan(delta))
Ra <- (24*60/pi) * Gsc * dr*(omega *sin(phi)*sin(delta) + sin(omega)*cos(phi)*cos(delta)) #nolint
1000*Ra # unit is kJ/m2
phi <- pi * lat / 180
J <- daynumber(datestring)
dr <- 1 + 0.033 * cos(2 * pi * J / 365)
delta <- 0.409 * sin(2 * pi * J / 365 - 1.39)
omega <- acos(-tan(phi) * tan(delta))
Ra <- (24 * 60 / pi) * Gsc * dr * (omega * sin(phi) * sin(delta) +
sin(omega) * cos(phi) * cos(delta))
1000 * Ra # unit is kJ/m2
}

# Derive daynumber within year (1-366)
daynumber <- function(datestring_YYYYMMDD) {
#' Derive daynumber within year (1-366)
#' @param datestring integer in format YYYYMMDD
#' @keywords internal
daynumber <- function(datestring) {
dpm <- c(0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334)
id <- floor( datestring_YYYYMMDD %% 100)
im <- floor((datestring_YYYYMMDD %% 10000) / 100) #nolint
iy <- floor( datestring_YYYYMMDD / 10000) %% 4
dnr <- dpm[im] + id + (iy==0 & im >2)
id <- floor( datestring %% 100)
im <- floor( (datestring %% 10000) / 100)
iy <- floor( datestring / 10000) %% 4
dnr <- dpm[im] + id + (iy == 0 & im > 2)
return(dnr)
}
26 changes: 13 additions & 13 deletions R/TemperatureTransformation.R
Expand Up @@ -57,13 +57,13 @@ tm_trans_KNMI14 <- function(obs,
for (im in 1:12) {
# get change factors for specific station and month
percentile.changes <- as.numeric(
deltas[which(deltas[, 1]==regio & deltas[, 2]==im), -1:-2] )
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>
X <- obs[days.im, is + 1] # select obs in month <im> of station <is>
Y <- rep(NA, length(X))
# observed percentiles
X.percentiles <- as.numeric(quantile(X, c(1, 5, 50, 95, 99)/100, na.rm=T))
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

Expand All @@ -74,30 +74,30 @@ tm_trans_KNMI14 <- function(obs,
# 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])
a <- (Y.percentiles[ip] - Y.percentiles[ip - 1]) /
(X.percentiles[ip] - X.percentiles[ip - 1])
b <- Y.percentiles[ip] - a * X.percentiles[ip]
# apply function to all temperatures below 5th percentile
Y[x.ids] <- a * X[x.ids] + b

for (ip in 3:(length(X.percentiles)-1)) {
for (ip in 3:(length(X.percentiles) - 1)) {
# 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])
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)
# 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])
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]
Y[x.ids] <- a * X[x.ids] + b

fut[days.im, is+1] <- Y
fut[days.im, is + 1] <- Y
} # END MONTHLY LOOP

} # END OF TRANSFORMATION LOOP
Expand Down
39 changes: 20 additions & 19 deletions R/knmitransformer.R
Expand Up @@ -166,10 +166,10 @@ TransformPrecip <- function(input,
#' global radiation and calculates the Makkink evaporation for 'future time
#' series' that match a certain climate
#' @inheritParams TransformTemp
#' @param input_tg Name of the input file for temperature
#' @param input_rsds Name of the input file for radiation
#' @param inputTemp Name of the input file for temperature
#' @param inputRad Name of the input file for radiation
#' @export
TransformEvap <- function(input_tg, input_rsds,
TransformEvap <- function(inputTemp, inputRad,
scenario = "GL",
horizon = 2030,
regions = "NLD",
Expand All @@ -180,32 +180,33 @@ TransformEvap <- function(input_tg, input_rsds,

CheckPeriod(horizon)

userProvided <- CheckIfUserProvided(input_tg) |
CheckIfUserProvided(input_rsds)
userProvided <- CheckIfUserProvided(inputTemp) |
CheckIfUserProvided(inputRad)


rsds_input <- TransformRadiation(input = input_rsds, scenario=scenario,
horizon = horizon,
rounding = FALSE)
tg_input <- TransformTemp(input = input_tg, var="tg", scenario=scenario,
horizon = horizon, regions = regions,
rounding = FALSE)
inputRad <- TransformRadiation(input = inputRad, scenario = scenario,
horizon = horizon,
rounding = FALSE)
inputTemp <- TransformTemp(input = inputTemp, var = "tg", scenario = scenario,
horizon = horizon, regions = regions,
rounding = FALSE)

rsds <- rsds_input[-c(1:5), ]
tg <- tg_input[-c(1:5), ]
if (!all(rsds_input[1:5, ] == tg_input[1:5, ])) {
flog.error("Same stations should be used for temperature and radiation")
stop("Same stations should be used for temperature and radiation")
rsds <- inputRad[-c(1:5), ]
tg <- inputTemp[-c(1:5), ]
if (!all(inputRad[1:5, ] == inputTemp[1:5, ])) {
err <- "Same stations should be used for temperature and radiation"
flog.error(err)
stop(err)
}

fut <- rsds
fut[, 2:ncol(fut)] <- NA
fut[, 2:ncol(fut)] <- makkink(tg[, 2:ncol(fut), with=FALSE],
rsds[, 2:ncol(fut), with=FALSE])
fut[, 2:ncol(fut)] <- makkink(tg[, 2:ncol(fut), with = FALSE],
rsds[, 2:ncol(fut), with = FALSE])

# Have to add a test to make sure that the header here is the same as the
# header in the regressionInput files
header <- rsds_input[1:5, ]
header <- inputRad[1:5, ]
H.comments <- "# Makkink Evaporation [mm] as derived from transformed tg & rsds "

# OUTPUT #####################################################################
Expand Down
5 changes: 2 additions & 3 deletions inst/.lintr
@@ -1,7 +1,6 @@
linters: with_defaults(
infix_spaces_linter = NULL, # 569
line_length_linter(90), # 189
snake_case_linter = NULL, # should be incorporated
line_length_linter(90), #slightly longer lines
snake_case_linter, # object snake case
camel_case_linter = NULL, # 60
commented_code_linter = NULL, # 52
NULL
Expand Down
17 changes: 17 additions & 0 deletions man/ObtainAngotRadiation.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

6 changes: 3 additions & 3 deletions man/TransformEvap.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

15 changes: 15 additions & 0 deletions man/daynumber.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

0 comments on commit 7cf953b

Please sign in to comment.