Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#' Calculate averages (mean) for defined window sizes
#'
#' Define window sizes (nyear_window) to be used to calculate averages (mean)
#' for each window (`dim(x)[3] / nyear_window`). Instead of discrete windows,
#' also moving averages can be computed as well as years inbetween interpolated.
#'
#' @param x LPJmL output array with `dim(x)=c(cell, month, year)`
#'
#' @param nyear_window integer, if supplied it defines the years for each window
#' to be averaged over in `dim(x)[3]`. If `nyear_window == 1` values are used
#' directly (instead of calculating an average). nyear_window has to be smaller
#' than `dim(x)[3]` and `dim(x)[3]` is ideally a multipe of nyear_window.
#' Defaults to `NULL`
#'
#' @param moving_average logical. If `TRUE` moving average is computed. start
#' and end are interpolated using spline interpolation.
#'
#' @param interpolate logical. If `TRUE` and nyear_window is defined (with
#' `moving_average == FALSE` years are interpolated (spline) to return array
#' with same dimensions as `x` (mainly`dim(x)[3]` -> year).
#'
#' @param nyear_reference integer, if supplied (default NULL), it defines a
#' time_span for ideally reference runs to be used as a baseline. E.g.
#' `nyear_reference = 30` to be used for preindustrial climate reference.
#'
#' @return array with same amount of cells and months as x. 3rd dimension is
#' defined by nyear_window, basically `dim(x)[3]/nyear_window` or equal to
#' dim(x)[3] if `moving_average == TRUE` or `interpolate == TRUE`
#'
#' @md
#' @importFrom magrittr %>%
#' @export
average_nyear_window <- function(x, # nolint
nyear_window = NULL,
moving_average = FALSE,
interpolate = FALSE,
nyear_reference = NULL) {
third_dim <- names(dim(x))[
!names(dim(x)) %in% c("cell", "year")
] %>% {
if (rlang::is_empty(.)) NULL else .
}
if (length(third_dim) > 1) {
stop(paste0("x has to have dimensions \"cell\", \"year\" and can have ",
"one third dimension (e.g. \"month\", \"year\""))
}
# check validity of x dimensions
if (!all(names(dim(x)) %in% c("cell", third_dim, "year"))) {
stop(paste0("x has to have dimensions \"cell\"",
ifelse(is.null(third_dim),
"",
paste0(", \"", third_dim, "\"")),
" and \"year\"."))
}
# moving average function - spline interpolation to fill NAs at start/end
moving_average_fun <- function(x, n) {
stats::filter(x, rep(1 / n, n), sides = 2) %>%
zoo::na.spline()
}
# utility function to interpolate inbetween averaging windows via spline
interpolate_spline <- function(x, y, nyear_window) {
rep(NA, dim(y)["year"]) %>%
`[<-`(seq(round(nyear_window / 2), dim(y)["year"], nyear_window),
value = x) %>%
zoo::na.spline()
}
# if nyear_window is supplied not all years are used for averaging
if (!is.null(nyear_window)) {
if (!is.null(nyear_reference)) {
orig_x <- x
x <- lpjmlkit::asub(x, year = 1:nyear_reference)
}
# only valid for nyear_window < years of x (dim(x)[3])
if (nyear_window > 1 & nyear_window <= dim(x)["year"]) {
# check if multiple (can also be left out)
# if (dim(x)[3] %% nyear_window == 0) {
if (moving_average) {
y <- aperm(apply(x,
c("cell", third_dim),
moving_average_fun,
nyear_window),
c("cell", third_dim, ""))
} else {
# calculate mean for discret windows/bins with size of nyear_window
y <- array(x,
# set correct dimensions (with names)
dim = c(dim(x)[c("cell", third_dim)],
year = nyear_window,
window = dim(x)[["year"]] / nyear_window),
# set correct dimensions names with nyear and windows
dimnames = append(dimnames(x)[c("cell", third_dim)],
list(year = seq_len(nyear_window),
window = dimnames(x)[["year"]][
seq(round(nyear_window / 2),
dim(x)[["year"]],
nyear_window)
]))) %>%
apply(c("cell", third_dim, "window"), mean)
if (interpolate) {
y <- aperm(apply(y,
c("cell", third_dim),
interpolate_spline,
x,
nyear_window),
c("cell", third_dim, ""))
}
}
# }
} else if (nyear_window == 1) {
y <- x
} else {
stop(paste0("Amount of nyear_window (", nyear_window, ") not supported."))
}
# recycle nyear_reference subset for original x (years)
if (!is.null(nyear_reference)) {
# multiple factor
nmultiple <- round(dim(orig_x)[["year"]] / nyear_reference)
replace_multiple_id <- nmultiple * dim(y)[[3]]
# if average window is returned as year dimension
if (!moving_average & !interpolate) {
z <- array(NA,
dim = c(dim(y)[c("cell", third_dim)],
window = replace_multiple_id),
dimnames = append(dimnames(y)[c("cell", third_dim)],
list(window = rep(dimnames(y)[[3]],
nmultiple))))
# return as original year dimension
} else {
# years vector also for non multiples (subset only partly recycled)
years <- rep(NA, dim(orig_x)[["year"]]) %>%
`[<-`(, value = dimnames(x)[["year"]]) %>%
suppressWarnings()
z <- array(NA,
dim = dim(orig_x),
dimnames = append(dimnames(y)[c("cell", third_dim)],
list(year = years)))
}
# recycle subset y for rest of original (x) years in z
z[, , seq_len(replace_multiple_id)] <- y
# check if not multiple - then only partly recylce array
if ((dim(z)[3] - replace_multiple_id) > 0) {
z[, , (replace_multiple_id + 1):dim(z)[3]] <- (
y[, , seq_len(dim(z)[3] - replace_multiple_id)]
)
}
return(z)
} else {
# rename dimnames of array
if (all(dim(y) == dim(x))) {
dim(y) <- dim(x)
dimnames(y) <- dimnames(x)
}
}
} else {
y <- apply(x, c("cell", third_dim), mean)
if (is.null(dim(y))) {
y <- array(
y,
dim = c(cell = dim(x)[["cell"]], 1),
dimnames = list(cell = dimnames(x)[["cell"]], 1)
)
}
}
return(y)
}