-
Notifications
You must be signed in to change notification settings - Fork 1
/
pfocal.R
204 lines (184 loc) · 7.15 KB
/
pfocal.R
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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
#' Fast, parallel implementation of grid data convolution
#'
#' A fast, parallel implementation of convolutions for grid-type
#' data (matrices, rasters and other grid based objects).
#'
#' @param data **\[matrix-type\]** Grid to compute onto.
#' @param kernel **\[matrix\]** Computation kernel (neighborhood).
#' @param edge_value **\[numeric\]** Numeric value, `NA` or `NaN`.The value to
#' give at the edge of the kernel matrix when the moving window overflows
#' the data grid.
#' @param transform_function **\[character\]** The function to apply to the
#' cell values covered by the kernel. For possible values, see
#' [pfocal_info_transform()]. Default to `"MULTIPLY"`.
#' @param reduce_function **\[character\]** The function to apply to the kernel
#' values after the function passed in `transform_function` has been applied
#' (the function that summarize the data). For possible values, see
#' [pfocal_info_reduce()]. Default to `"SUM"`.
#' @param mean_divider **\[character\]** Optional, allows to specify how the
#' final value at each cell is divided by a value that can be function of
#' the intermediate data resulting of applying `transform_function`. For
#' possible values, see [pfocal_info_mean_divisor()]. Default to "ONE" (for
#' no division).
#' @param variance **\[logical\]** Whether to return the "variance" of the
#' intermediate values at each point (for more details please see
#' [pfocal_info_variance()]). Default to `FALSE` (just returns the value
#' at each point).
#' @param na.rm **\[NA OR logical\]** The behavior to adopt for dealing with
#' missing values, default to `NA` (faster option). For possible values see
#' [pfocal_info_nan_policy()].
#' @param mp **\[logical\]** Whether to use the open_mp implementation,
#' default to `TRUE`.
#' @param debug_use_r_implementation **\[logical\]** Used for debugging purposes
#' whether to use the slow R implementation instead of the fast C++ code.
#' Default to `FALSE`.
#' @param ... None used at the moment
#' .
#'
#' @return
#' The updated, convoluted grid.
#'
#' @details
#' Note that the memory allocation for the output is of size
#' `sizeof(double) * ncol * nrow` and for the intermediate
#' values, `sizeof(double) * (ncol + kernel_ncol) * (nrow + kernel_nrow/2)`.
#'
#' @examples
#'
#' data <- matrix(nrow = 10, ncol = 10, data = runif(10 * 10))
#' kernel <- matrix(1 / 9, nrow = 3, ncol = 3)
#' pfocal(data = data, kernel = kernel)
#'
#' @export
pfocal <- function(data, kernel, edge_value = 0, transform_function = "MULTIPLY",
reduce_function = "SUM", mean_divider = "ONE", variance = FALSE,
na.rm = NA, mp = TRUE, debug_use_r_implementation = FALSE, ...) {
if (methods::is(data, "matrix")) {
return(pfocal.matrix(
data, kernel, edge_value, transform_function,
reduce_function, mean_divider, variance,
na.rm, mp, debug_use_r_implementation, ...
))
} else if (methods::is(data, "RasterLayer")){
data_matrix <- raster::as.matrix(data)
data_transformed <- pfocal.matrix(
data_matrix, kernel, edge_value, transform_function,
reduce_function, mean_divider, variance,
na.rm, mp, debug_use_r_implementation, ...
)
raster::values(data) <- data_transformed
return(data)
} else if (methods::is(data, "SpatRaster")){
data_matrix <- terra::as.matrix(data, wide = TRUE)
data_transformed <- pfocal.matrix(
data_matrix, kernel, edge_value, transform_function,
reduce_function, mean_divider, variance,
na.rm, mp, debug_use_r_implementation, ...
)
terra::values(data) <- data_transformed
return(data)
} else if (methods::is(data, "stars")) {
return(pfocal.stars(
data, kernel, edge_value, transform_function,
reduce_function, mean_divider, variance,
na.rm, mp, debug_use_r_implementation, ...
))
} else {
stop('unsupported type, x must be a "matrix", "RasterLayer", "SpatRaster", or "stars"')
}
}
# Matrix routine ----------------------------------------------------------
pfocal.matrix <- function(data, kernel, edge_value = NA, transform_function = "MULTIPLY",
reduce_function = "SUM", mean_divider = "ONE", variance = FALSE,
na.rm = NA, mp = TRUE, debug_use_r_implementation = FALSE, ...) {
.pfocal(
data, kernel, edge_value, transform_function, reduce_function, mean_divider,
variance, na.rm, mp, debug_use_r_implementation
)
}
# Stars object routine ----------------------------------------------------
pfocal.stars <- function(data, ...) {
# Code from github.com/michaeldorman/starsExtra R/focal2.R:focal2
template <- data
input <- starsExtra::layer_to_matrix(template, check = TRUE)
output <- pfocal.matrix(input, ...)
# Back to 'stars'
output <- t(output)
template[[1]] <- output
# Return
return(template)
}
# General routine ---------------------------------------------------------
.pfocal <- function(data, kernel, edge_value = NA, transform_function = "MULTIPLY",
reduce_function = "SUM", mean_divider = "ONE", variance = FALSE,
na.rm = NA, mp = TRUE, debug_use_r_implementation = FALSE) {
pfocal_f <- if (debug_use_r_implementation) {
.p_focal_r
} else {
.p_focal_cpp
}
if (is.na(na.rm)) {
na.rm <- "NA" # 'Fast' version
} else if (na.rm) {
na.rm <- "TRUE"
} else {
na.rm <- "FALSE"
}
t_index <- pfocal_narrow_transform(transform_function)
r_index <- pfocal_narrow_reduce(reduce_function)
n_index <- pfocal_narrow_nan_policy(na.rm)
m_index <- pfocal_narrow_mean_divisor(mean_divider)
v_index <- pfocal_narrow_variance(variance)
if (is.na(v_index)) {
stop("variance must be logical or 'TRUE' or 'FALSE'")
}
if (is.na(m_index)) {
stop(paste0(
"Unknown mean_divider. If you don't want to ",
"divide the output, leave this set to 'ONE'"
))
}
if (is.na(n_index)) {
stop(paste0(
"Unknown na.rm value. It must be logical, or NA",
" if you want more speed and don't care about",
" how NA is handeled"
))
}
if (is.na(t_index)) {
# special transform cases
if (toupper(transform_function) == "SUBTRACT") {
data <- data * -1
t_index <- pfocal_narrow_transform("ADD")
} else if (toupper(transform_function) == "L_DIVIDE") {
data <- 1 / data
t_index <- pfocal_narrow_transform("MULTIPLY")
} else if (toupper(transform_function) %in% c("DIVIDE", "R_DIVIDE")) {
kernel <- 1 / kernel
t_index <- pfocal_narrow_transform("MULTIPLY")
} else {
stop("Unknown transform_function")
}
}
if (is.na(r_index)) {
if (toupper(reduce_function) == "RANGE") {
return(pfocal_f(
data, kernel, edge_value, t_index, pfocal_narrow_reduce("MAX"),
n_index, m_index, v_index, mp
) -
pfocal_f(
data, kernel, edge_value, t_index, pfocal_narrow_reduce("MIN"),
n_index, m_index, v_index, mp
))
} else {
stop("Unknown reduce_function")
}
}
pf <- pfocal_f(
data, kernel, edge_value, t_index, r_index, n_index, m_index,
v_index, mp
)
pf[is.nan(pf)] <- NA
pf[is.infinite(pf)] <- Inf
return(pf)
}