1 |
# Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
# All rights reserved. |
|
3 |
# |
|
4 |
# This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
# You may not use this file except in compliance with the License. You may |
|
6 |
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
# |
|
8 |
# Unless required by applicable law or agreed to in writing, software |
|
9 |
# distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
# See the License for the specific language governing permissions and |
|
12 |
# limitations under the License. |
|
13 | ||
14 |
if (!isGeneric("exact_extract")) { |
|
15 | 186x |
setGeneric("exact_extract", function(x, y, ...) |
16 | 186x |
standardGeneric("exact_extract")) |
17 |
} |
|
18 | ||
19 |
#' Extract or summarize values from Raster* objects |
|
20 |
#' |
|
21 |
#' Extracts the values of cells in a Raster* that are covered by a |
|
22 |
#' simple feature collection containing polygonal geometries, as well as the |
|
23 |
#' fraction of each cell that is covered by the polygon. Returns either |
|
24 |
#' the result of a summary operation or function applied to the values |
|
25 |
#' and coverage fractions (if \code{fun} is specified), or a data frame |
|
26 |
#' containing the values and coverage fractions themselves (if \code{fun} |
|
27 |
#' is \code{NULL}.) |
|
28 |
#' |
|
29 |
#' The value of \code{fun} may be set to a string (or vector of strings) |
|
30 |
#' representing summary operations supported by the exactextract library. |
|
31 |
#' If the input raster has a single layer and a single summary operation |
|
32 |
#' is specified, \code{exact_extract} will return a vector with the result |
|
33 |
#' of the summary operation for each feature in the input. If the input |
|
34 |
#' raster has multiple layers, or if multiple summary operations are specified, |
|
35 |
#' \code{exact_extract} will return a data frame with a row for each feature |
|
36 |
#' and a column for each summary operation / layer combination. (The |
|
37 |
#' \code{force_df} can be used to always return a data frame instead of a vector.) |
|
38 |
#' In all of the summary operations, \code{NA} values in the raster are ignored |
|
39 |
#' (i.e., \code{na.rm = TRUE}.) |
|
40 |
#' |
|
41 |
#' The following summary operations are supported: |
|
42 |
#' |
|
43 |
#' \itemize{ |
|
44 |
#' \item{\code{min} - the minimum defined value in any raster cell wholly or |
|
45 |
#' partially covered by the polygon} |
|
46 |
#' \item{\code{max} - the maximum defined value in any raster cell wholly or |
|
47 |
#' partially covered by the polygon} |
|
48 |
#' \item{\code{count} - the sum of fractions of raster cells with defined values |
|
49 |
#' covered by the polygon} |
|
50 |
#' \item{\code{sum} - the sum of defined raster cell values, multiplied by |
|
51 |
#' the fraction of the cell that is covered by the polygon} |
|
52 |
#' \item{\code{mean} - the mean cell value, weighted by the fraction of each cell |
|
53 |
#' that is covered by the polygon} |
|
54 |
#' \item{\code{median} - the median cell value, weighted by the fraction of each |
|
55 |
#' cell that is covered by the polygon} |
|
56 |
#' \item{\code{quantile} - arbitrary quantile(s) of cell values, specified in |
|
57 |
#' \code{quantiles}, weighted by the fraction of each |
|
58 |
#' cell that is covered by the polygon} |
|
59 |
#' \item{\code{mode} - the most common cell value, weighted by the fraction of |
|
60 |
#' each cell that is covered by the polygon. Where multiple |
|
61 |
#' values occupy the same maximum number of weighted cells, |
|
62 |
#' the largest value will be returned.} |
|
63 |
#' \item{\code{majority} - synonym for \code{mode}} |
|
64 |
#' \item{\code{minority} - the least common cell value, weighted by the fraction |
|
65 |
#' of each cell that is covered by the polygon. Where |
|
66 |
#' multiple values occupy the same minimum number of |
|
67 |
#' weighted cells, the smallest value will be returned.} |
|
68 |
#' \item{\code{variety} - the number of distinct values in cells that are wholly |
|
69 |
#' or partially covered by the polygon.} |
|
70 |
#' \item{\code{variance} - the population variance of cell values, weighted by the |
|
71 |
#' fraction of each cell that is covered by the polygon.} |
|
72 |
#' \item{\code{stdev} - the population standard deviation of cell values, weighted |
|
73 |
#' by the fraction of each cell that is covered by the polygon.} |
|
74 |
#' \item{\code{coefficient_of_variation} - the population coefficient of variation of |
|
75 |
#' cell values, weighted by the fraction of each cell that is |
|
76 |
#' covered by the polygon.} |
|
77 |
#' \item{\code{weighted_mean} - the mean cell value, weighted by the product of |
|
78 |
#' the fraction of each cell covered by the polygon |
|
79 |
#' and the value of a second weighting raster provided |
|
80 |
#' as \code{weights}} |
|
81 |
#' \item{\code{weighted_sum} - the sum of defined raster cell values, multiplied by |
|
82 |
#' the fraction of each cell that is covered by the polygon |
|
83 |
#' and the value of a second weighting raster provided |
|
84 |
#' as \code{weights}} |
|
85 |
#' } |
|
86 |
#' |
|
87 |
#' Alternatively, an R function may be provided as \code{fun}. The function will be |
|
88 |
#' called for each feature with with vectors of cell values and weights as arguments. |
|
89 |
#' \code{exact_extract} will then return a vector of the return values of \code{fun}. |
|
90 |
#' |
|
91 |
#' If \code{fun} is not specified, \code{exact_extract} will return a list with |
|
92 |
#' one data frame for each feature in the input feature collection. The data |
|
93 |
#' frame will contain a column with values from each layer in the input `Raster*`, |
|
94 |
#' and a final column indicating the fraction of the cell that is covered by the |
|
95 |
#' polygon. |
|
96 |
#' |
|
97 |
#' @param x a \code{RasterLayer}, \code{RasterStack}, or \code{RasterBrick} |
|
98 |
#' @param y a sf object with polygonal geometries |
|
99 |
#' @param fun an optional function or character vector, as described below |
|
100 |
#' @param weights a weighting raster to be used with the \code{weighted_mean} |
|
101 |
#' and \code{weighted_sum} summary operations. |
|
102 |
#' @param quantiles quantiles to be computed when \code{fun == 'quantile'} |
|
103 |
#' @param append_cols when \code{fun} is not \code{NULL}, an optional |
|
104 |
#' character vector of columns from \code{y} to be |
|
105 |
#' included in returned data frame. |
|
106 |
#' @param force_df always return a data frame instead of a vector, even if |
|
107 |
#' \code{x} has only one layer and \code{fun} has length 1 |
|
108 |
#' @param full_colnames include the names of \code{x} in the names of the |
|
109 |
#' returned data frame, even if \code{x} has only one |
|
110 |
#' layer. This is useful when the results of multiple |
|
111 |
#' calls to \code{exact_extract} are combined with |
|
112 |
#' \code{cbind}. |
|
113 |
#' @param include_cell if \code{TRUE}, and \code{fun} is \code{NULL}, augment |
|
114 |
#' the returned data frame for each feature with a column |
|
115 |
#' for the cell index (\code{cell}). If \code{TRUE} and |
|
116 |
#' \code{fun} is not \code{NULL}, add \code{cell} to the |
|
117 |
#' data frame passed to \code{fun} for each feature. |
|
118 |
#' @param include_cols an optional character vector of column names in |
|
119 |
#' \code{y} to be added to the data frame for each |
|
120 |
#' feature that is either returned (when \code{fun} is |
|
121 |
#' \code{NULL}) or passed to \code{fun}. |
|
122 |
#' @param include_xy if \code{TRUE}, and \code{fun} is \code{NULL}, augment |
|
123 |
#' the returned data frame for each feature with columns |
|
124 |
#' for cell center coordinates (\code{x} and \code{y}). If |
|
125 |
#' \code{TRUE} and \code{fun} is not \code{NULL}, add |
|
126 |
#' \code{x} and {y} to the data frame passed to \code{fun} |
|
127 |
#' for each feature. |
|
128 |
#' @param stack_apply if \code{TRUE}, apply \code{fun} to each layer of |
|
129 |
#' \code{x} independently. If \code{FALSE}, apply \code{fun} |
|
130 |
#' to all layers of \code{x} simultaneously. |
|
131 |
#' @param max_cells_in_memory the maximum number of raster cells to load at |
|
132 |
#' a given time when using a named summary operation |
|
133 |
#' for \code{fun} (as opposed to a function defined using |
|
134 |
#' R code). If a polygon covers more than \code{max_cells_in_memory} |
|
135 |
#' raster cells, it will be processed in multiple chunks. |
|
136 |
#' @param progress if \code{TRUE}, display a progress bar during processing |
|
137 |
#' @param ... additional arguments to pass to \code{fun} |
|
138 |
#' @return a vector or list of data frames, depending on the type of \code{x} and the |
|
139 |
#' value of \code{fun} (see Details) |
|
140 |
#' @examples |
|
141 |
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10) |
|
142 |
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))') |
|
143 |
#' |
|
144 |
#' # named summary operation on RasterLayer, returns vector |
|
145 |
#' exact_extract(rast, poly, 'mean') |
|
146 |
#' |
|
147 |
#' # two named summary operations on RasterLayer, returns data frame |
|
148 |
#' exact_extract(rast, poly, c('min', 'max')) |
|
149 |
#' |
|
150 |
#' # named summary operation on RasterStack, returns data frame |
|
151 |
#' stk <- raster::stack(list(a=rast, b=sqrt(rast))) |
|
152 |
#' exact_extract(stk, poly, 'mean') |
|
153 |
#' |
|
154 |
#' # named weighted summary operation, returns vector |
|
155 |
#' weights <- raster::raster(matrix(runif(100), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10) |
|
156 |
#' exact_extract(rast, poly, 'weighted_mean', weights=weights) |
|
157 |
#' |
|
158 |
#' # custom summary function, returns vector |
|
159 |
#' exact_extract(rast, poly, function(value, cov_frac) length(value[cov_frac > 0.9])) |
|
160 |
#' |
|
161 |
#' @name exact_extract |
|
162 |
NULL |
|
163 | ||
164 |
#' @import sf |
|
165 |
#' @import raster |
|
166 |
#' @useDynLib exactextractr |
|
167 |
#' @rdname exact_extract |
|
168 |
#' @export |
|
169 |
setMethod('exact_extract', signature(x='Raster', y='sf'), |
|
170 |
function(x, y, fun=NULL, ..., |
|
171 |
include_xy=FALSE, |
|
172 |
progress=TRUE, |
|
173 |
max_cells_in_memory=30000000, |
|
174 |
include_cell=FALSE, |
|
175 |
force_df=FALSE, |
|
176 |
full_colnames=FALSE, |
|
177 |
stack_apply=FALSE, |
|
178 |
append_cols=NULL, |
|
179 |
include_cols=NULL, |
|
180 |
quantiles=NULL) { |
|
181 | 13x |
.exact_extract(x, y, fun=fun, ..., |
182 | 13x |
include_xy=include_xy, |
183 | 13x |
progress=progress, |
184 | 13x |
max_cells_in_memory=max_cells_in_memory, |
185 | 13x |
include_cell=include_cell, |
186 | 13x |
force_df=force_df, |
187 | 13x |
full_colnames=full_colnames, |
188 | 13x |
stack_apply=stack_apply, |
189 | 13x |
append_cols=append_cols, |
190 | 13x |
include_cols=include_cols, |
191 | 13x |
quantiles=quantiles) |
192 |
}) |
|
193 | ||
194 |
# Return the number of standard (non-...) arguments in a supplied function that |
|
195 |
# do not have a default value. This is used to fail if the summary function |
|
196 |
# provided by the user cannot accept arguments of values and weights. |
|
197 |
.num_expected_args <- function(fun) { |
|
198 | 40x |
a <- formals(args(fun)) |
199 | 40x |
a <- a[names(a) != '...'] |
200 | 40x |
sum(sapply(a, nchar) == 0) |
201 |
} |
|
202 | ||
203 |
emptyVector <- function(rast) { |
|
204 | ! |
switch(substr(raster::dataType(rast), 1, 3), |
205 | ! |
LOG=logical(), |
206 | ! |
INT=integer(), |
207 | ! |
numeric()) |
208 |
} |
|
209 | ||
210 |
.exact_extract <- function(x, y, fun=NULL, ..., |
|
211 |
weights=NULL, |
|
212 |
include_xy=FALSE, |
|
213 |
progress=TRUE, |
|
214 |
max_cells_in_memory=30000000, |
|
215 |
include_cell=FALSE, |
|
216 |
force_df=FALSE, |
|
217 |
full_colnames=FALSE, |
|
218 |
stack_apply=FALSE, |
|
219 |
append_cols=NULL, |
|
220 |
include_cols=NULL, |
|
221 |
quantiles=NULL) { |
|
222 | 181x |
if(!is.null(append_cols)) { |
223 | 2x |
if (!inherits(y, 'sf')) { |
224 | ! |
stop(sprintf('append_cols only supported for sf arguments (received %s)', |
225 | ! |
paste(class(y), collapse = ' '))) |
226 |
} |
|
227 | ||
228 | 2x |
force_df <- TRUE |
229 |
} |
|
230 | ||
231 | 181x |
if(sf::st_geometry_type(y, by_geometry = FALSE) == 'GEOMETRY') { |
232 | 2x |
if (!all(sf::st_dimension(y) == 2)) { |
233 | 1x |
stop("Features in sfc_GEOMETRY must be polygonal") |
234 |
} |
|
235 | 1x |
y <- sf::st_cast(y, 'MULTIPOLYGON') |
236 |
} |
|
237 | ||
238 | 180x |
if(!is.null(weights)) { |
239 | 63x |
if (!startsWith(class(weights), 'Raster')) { |
240 | 1x |
stop("Weights must be a Raster object.") |
241 |
} |
|
242 | ||
243 | 62x |
if (is.character(fun) && !any(startsWith(fun, "weighted"))) { |
244 | 9x |
warning("Weights provided but no requested operations use them.") |
245 |
} |
|
246 | ||
247 | 62x |
if (!is.na(sf::st_crs(x))) { |
248 | 53x |
if (is.na(sf::st_crs(weights))) { |
249 | 1x |
warning("No CRS specified for weighting raster; assuming it has the same CRS as the value raster.") |
250 | 52x |
} else if (sf::st_crs(x) != sf::st_crs(weights)) { |
251 | 1x |
stop("Weighting raster does not have the same CRS as value raster.") |
252 |
} |
|
253 |
} |
|
254 |
} |
|
255 | ||
256 | 178x |
if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) { |
257 | 1x |
warning("No CRS specified for raster; assuming it has the same CRS as the polygons.") |
258 | 177x |
} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) { |
259 | 2x |
warning("No CRS specified for polygons; assuming they have the same CRS as the raster.") |
260 | 175x |
} else if(sf::st_crs(x) != sf::st_crs(y)) { |
261 | 1x |
y <- sf::st_transform(y, sf::st_crs(x)) |
262 | 1x |
warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")") |
263 |
} |
|
264 | ||
265 | 178x |
if (!is.null(fun) && !is.character(fun) && .num_expected_args(fun) < 2) { |
266 | 5x |
stop("exact_extract was called with a function that does not appear to ", |
267 | 5x |
"be of the form `function(values, coverage_fractions, ...)`") |
268 |
} |
|
269 | ||
270 | 173x |
if (is.character(fun)) { |
271 | 117x |
if (length(list(...)) > 0) { |
272 | 1x |
stop("exact_extract was called with a named summary operation that ", |
273 | 1x |
"does not accept additional arguments ...") |
274 |
} |
|
275 | ||
276 | 116x |
if (include_xy) { |
277 | 1x |
stop("include_xy must be FALSE for named summary operations") |
278 |
} |
|
279 | ||
280 | 115x |
if (include_cell) { |
281 | 1x |
stop("include_cell must be FALSE for named summary operations") |
282 |
} |
|
283 | ||
284 | 114x |
if (!is.null(include_cols)) { |
285 | 1x |
stop("include_cols not supported for named_summary operations (see argument append_cols)") |
286 |
} |
|
287 |
} |
|
288 | ||
289 | 169x |
geoms <- sf::st_geometry(y) |
290 | ||
291 | 169x |
if (progress && length(geoms) > 1) { |
292 | 5x |
n <- length(geoms) |
293 | 5x |
pb <- utils::txtProgressBar(min = 0, max = n, initial=0, style=3) |
294 | 5x |
update_progress <- function() { |
295 | 54x |
i <- 1 + utils::getTxtProgressBar(pb) |
296 | 54x |
utils::setTxtProgressBar(pb, i) |
297 | 54x |
if (i == n) { |
298 | 5x |
close(pb) |
299 |
} |
|
300 |
} |
|
301 |
} else { |
|
302 | 164x |
update_progress <- function() {} |
303 |
} |
|
304 | ||
305 | 169x |
tryCatch({ |
306 | 169x |
x <- raster::readStart(x) |
307 | 169x |
if (!is.null(weights)) { |
308 | 61x |
weights <- raster::readStart(weights) |
309 |
} |
|
310 | ||
311 | 169x |
if (is.character(fun)) { |
312 |
# Compute all stats in C++. |
|
313 |
# CPP_stats returns a matrix, which gets turned into a column by sapply |
|
314 |
# Results has one column per feature and one row per stat/raster layer |
|
315 | 113x |
if(is.null(weights) && any(.isWeighted(fun))) { |
316 | 2x |
stop("Weighted stat requested but no weights provided.") |
317 |
} |
|
318 | ||
319 | 111x |
results <- sapply(sf::st_as_binary(geoms, EWKB=TRUE), function(wkb) { |
320 | 146x |
ret <- CPP_stats(x, weights, wkb, fun, max_cells_in_memory, quantiles) |
321 | 137x |
update_progress() |
322 | 137x |
return(ret) |
323 |
}) |
|
324 | ||
325 | 102x |
if (!is.matrix(results) && !force_df) { |
326 |
# Single stat? Return a vector unless asked otherwise via force_df. |
|
327 | 84x |
return(results) |
328 |
} else { |
|
329 |
# Return a data frame with a column for each stat |
|
330 | 18x |
colnames <- .resultColNames(names(x), names(weights), fun, full_colnames, quantiles) |
331 | ||
332 | 18x |
if (is.matrix(results)) { |
333 | 16x |
results <- t(results) |
334 |
} else { |
|
335 | 2x |
results <- matrix(results, nrow=length(results)) |
336 |
} |
|
337 | ||
338 | 18x |
dimnames(results) <- list(NULL, colnames) |
339 | 18x |
ret <- as.data.frame(results) |
340 | ||
341 |
# drop duplicated columns (occurs when an unweighted stat is |
|
342 |
# requested alongside a weighted stat, with a stack of weights) |
|
343 | 18x |
ret <- ret[, unique(names(ret)), drop = FALSE] |
344 | ||
345 | 18x |
if (!is.null(append_cols)) { |
346 | 1x |
ret <- cbind(sf::st_drop_geometry(y[, append_cols]), ret) |
347 |
} |
|
348 | ||
349 | 18x |
return(ret) |
350 |
} |
|
351 |
} else { |
|
352 | 56x |
num_values <- raster::nlayers(x) |
353 | 56x |
num_weights <- ifelse(is.null(weights), 0, raster::nlayers(weights)) |
354 | 56x |
value_names <- names(x) |
355 | 56x |
weight_names <- names(weights) |
356 | ||
357 | 56x |
if (stack_apply || (num_values == 1 && num_weights <= 1)) { |
358 | 45x |
apply_layerwise <- TRUE |
359 | ||
360 | 45x |
if (num_values > 1 && num_weights > 1 && num_values != num_weights) { |
361 | 1x |
stop(sprintf("Can't apply function layerwise with stacks of %d value layers and %d layers", num_values, num_weights)) |
362 |
} |
|
363 | ||
364 | 44x |
result_names <- .resultColNames(value_names, weight_names, fun, full_colnames) |
365 | 44x |
num_results <- max(num_weights, num_values) |
366 | 44x |
ind <- .valueWeightIndexes(num_values, num_weights) |
367 |
} else { |
|
368 | 11x |
apply_layerwise <- FALSE |
369 |
} |
|
370 | ||
371 |
# For R summary functions and data frame output, we avoid using |
|
372 |
# input layer names if there is only one value/weight layer |
|
373 | 55x |
if (length(value_names) == 1) { |
374 | 41x |
value_names <- 'value' |
375 |
} |
|
376 | 55x |
if (length(weight_names) == 1) { |
377 | 11x |
weight_names <- 'weight' |
378 |
} |
|
379 | ||
380 | 55x |
ret <- lapply(seq_along(geoms), function(feature_num) { |
381 | 89x |
wkb <- sf::st_as_binary(geoms[[feature_num]], EWKB=TRUE) |
382 | ||
383 | 89x |
if (!is.null(include_cols)) { |
384 | 4x |
include_col_values <- sf::st_drop_geometry(y[feature_num, include_cols]) |
385 |
} else { |
|
386 | 85x |
include_col_values <- NULL |
387 |
} |
|
388 | ||
389 |
# only raise a disaggregation warning for the first feature |
|
390 | 89x |
warn_on_disaggregate <- feature_num == 1 |
391 | ||
392 | 89x |
col_list <- CPP_exact_extract(x, weights, wkb, include_xy, include_cell, include_col_values, value_names, weight_names, warn_on_disaggregate) |
393 | 89x |
if (!is.null(include_cols)) { |
394 |
# Replicate the include_cols vectors to be as long as the other columns, |
|
395 |
# so we can use quickDf |
|
396 | 4x |
nrow <- length(col_list$coverage_fraction) |
397 | 4x |
col_list[include_cols] <- lapply(col_list[include_cols], rep, nrow) |
398 |
} |
|
399 | 89x |
df <- .quickDf(col_list) |
400 | ||
401 | 89x |
update_progress() |
402 | ||
403 | 89x |
if (is.null(fun)) { |
404 | 22x |
return(df) |
405 |
} |
|
406 | ||
407 | 67x |
included_cols_df <- df[, !(names(df) %in% c(value_names, weight_names, 'coverage_fraction')), drop = FALSE] |
408 | 67x |
vals_df <- df[, value_names, drop = FALSE] |
409 | 67x |
weights_df <- df[, weight_names, drop = FALSE] |
410 | 67x |
cov_fracs <- df$coverage_fraction |
411 | ||
412 | 67x |
if (apply_layerwise) { |
413 | 61x |
result <- lapply(seq_len(num_results), function(i) { |
414 | 71x |
vx <- vals_df[, ind$values[i]] |
415 | 71x |
if (ncol(included_cols_df) > 0) { |
416 | 7x |
vx <- cbind(data.frame(value = vx), included_cols_df) |
417 |
} |
|
418 | 71x |
if (num_weights == 0) { |
419 | 59x |
fun(vx, cov_fracs, ...) |
420 |
} else { |
|
421 | 12x |
fun(vx, cov_fracs, weights_df[, ind$weights[i]], ...) |
422 |
} |
|
423 |
}) |
|
424 | ||
425 | 61x |
if (num_results == 1) { |
426 | 54x |
return(result[[1]]) |
427 |
} |
|
428 | ||
429 | 7x |
names(result) <- result_names |
430 | ||
431 | 7x |
.quickDf(result) |
432 |
} else { |
|
433 |
# Pass all layers to callback, to be handled together |
|
434 |
# Included columns (x/y/cell) are passed with the values. |
|
435 |
# Pass single-column data frames as vectors. |
|
436 | 6x |
vals_df <- .singleColumnToVector(cbind(vals_df, included_cols_df)) |
437 | 6x |
weights_df <- .singleColumnToVector(weights_df) |
438 | ||
439 | 6x |
if (num_weights == 0) { |
440 | 3x |
return(fun(vals_df, cov_fracs, ...)) |
441 |
} else { |
|
442 | 3x |
return(fun(vals_df, cov_fracs, weights_df, ...)) |
443 |
} |
|
444 |
} |
|
445 |
}) |
|
446 | ||
447 | 54x |
if (!is.null(fun)) { |
448 | 33x |
if (all(sapply(ret, is.data.frame))) { |
449 |
# function returned a data frame for each polygon? rbind them |
|
450 | 8x |
if (requireNamespace('dplyr', quietly = TRUE)) { |
451 | 8x |
ret <- dplyr::bind_rows(ret) # handle column name mismatches |
452 |
} else { |
|
453 | ! |
ret <- do.call(rbind, ret) |
454 |
} |
|
455 |
} else { |
|
456 |
# function returned something else; combine the somethings into |
|
457 |
# an array |
|
458 | 25x |
ret <- simplify2array(ret) |
459 | ||
460 | 25x |
if (force_df) { |
461 | 3x |
ret <- data.frame(result = ret) |
462 |
} |
|
463 |
} |
|
464 |
} |
|
465 | ||
466 | 54x |
if (!is.null(append_cols)) { |
467 | 1x |
ret <- cbind(sf::st_drop_geometry(y[, append_cols]), ret) |
468 |
} |
|
469 | ||
470 | 54x |
return(ret) |
471 |
} |
|
472 | 169x |
}, finally={ |
473 | 169x |
raster::readStop(x) |
474 | 169x |
if (!is.null(weights)) { |
475 | 61x |
raster::readStop(weights) |
476 |
} |
|
477 |
}) |
|
478 |
} |
|
479 | ||
480 |
# faster replacement for as.data.frame when input is a named list |
|
481 |
# with equal-length columns |
|
482 |
# from Advanced R, sec. 24.4.2 |
|
483 |
.quickDf <- function(lst) { |
|
484 | 96x |
class(lst) <- 'data.frame' |
485 | 96x |
attr(lst, 'row.names') <- .set_row_names(length(lst[[1]])) |
486 | 96x |
lst |
487 |
} |
|
488 | ||
489 |
.singleColumnToVector <- function(df) { |
|
490 | 12x |
if (ncol(df) == 1) { |
491 | 2x |
df[, 1] |
492 |
} else { |
|
493 | 10x |
df |
494 |
} |
|
495 |
} |
|
496 | ||
497 |
.appendXY <- function(vals_df, rast, first_row, nrow, first_col, ncol) { |
|
498 | ! |
if (nrow(vals_df) == 0) { |
499 | ! |
vals_df$x <- numeric() |
500 | ! |
vals_df$y <- numeric() |
501 |
} else { |
|
502 | ! |
x_coords <- raster::xFromCol(rast, col=seq(first_col, first_col + ncol - 1)) |
503 | ! |
y_coords <- raster::yFromRow(rast, row=seq(first_row, first_row + nrow - 1)) |
504 | ||
505 | ! |
vals_df$x <- rep(x_coords, times=nrow) |
506 | ! |
vals_df$y <- rep(y_coords, each=ncol) |
507 |
} |
|
508 | ||
509 | ! |
return(vals_df) |
510 |
} |
|
511 | ||
512 |
.appendCell <- function(vals_df, rast, first_row, nrow, first_col, ncol) { |
|
513 | ! |
if (nrow(vals_df) == 0) { |
514 | ! |
vals_df$cell <- numeric() |
515 |
} else { |
|
516 | ! |
rows <- rep(seq(first_row, first_row + nrow - 1), each = ncol) |
517 | ! |
cols <- rep(seq(first_col, first_col + ncol - 1), times = nrow) |
518 | ||
519 | ! |
vals_df$cell <- raster::cellFromRowCol(rast, row=rows, col=cols) |
520 |
} |
|
521 | ||
522 | ! |
return(vals_df) |
523 |
} |
|
524 | ||
525 |
#' @useDynLib exactextractr |
|
526 |
#' @rdname exact_extract |
|
527 |
#' @export |
|
528 |
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygonsDataFrame'), |
|
529 | 1x |
function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...)) |
530 | ||
531 |
#' @useDynLib exactextractr |
|
532 |
#' @rdname exact_extract |
|
533 |
#' @export |
|
534 |
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygons'), |
|
535 | 1x |
function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...)) |
536 | ||
537 |
#' @useDynLib exactextractr |
|
538 |
#' @rdname exact_extract |
|
539 |
#' @export |
|
540 |
setMethod('exact_extract', signature(x='Raster', y='sfc_MULTIPOLYGON'), .exact_extract) |
|
541 | ||
542 |
#' @useDynLib exactextractr |
|
543 |
#' @rdname exact_extract |
|
544 |
#' @export |
|
545 |
setMethod('exact_extract', signature(x='Raster', y='sfc_POLYGON'), .exact_extract) |
|
546 | ||
547 |
#' @useDynLib exactextractr |
|
548 |
#' @rdname exact_extract |
|
549 |
#' @export |
|
550 |
setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRY'), .exact_extract) |
1 |
# Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
# All rights reserved. |
|
3 |
# |
|
4 |
# This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
# You may not use this file except in compliance with the License. You may |
|
6 |
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
# |
|
8 |
# Unless required by applicable law or agreed to in writing, software |
|
9 |
# distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
# See the License for the specific language governing permissions and |
|
12 |
# limitations under the License. |
|
13 | ||
14 |
#' Return column names to be used for summary operations |
|
15 |
#' |
|
16 |
#' @param value_names names of value raster layers |
|
17 |
#' @param weight_names names of weighting raster layers |
|
18 |
#' @param fun functions or names of summary operations |
|
19 |
#' @param full_colnames return a complete column name even when there is no |
|
20 |
#' ambiguity? |
|
21 |
#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile} |
|
22 |
#' @return character vector of column names |
|
23 |
.resultColNames <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric()) { |
|
24 | 84x |
if (inherits(fun, 'standardGeneric')) { |
25 | 16x |
stat_names <- fun@generic[1] |
26 | 68x |
} else if (is.function(fun)) { |
27 | 24x |
stat_names <- 'fun' |
28 |
} else { |
|
29 | 44x |
stat_names <- fun |
30 |
} |
|
31 | ||
32 | 84x |
quantile_index = which(stat_names == 'quantile') |
33 | 84x |
if (length(quantile_index) != 0) { |
34 | 1x |
stat_names <- c(stat_names[seq_along(stat_names) < quantile_index], |
35 | 1x |
sprintf('q%02d', as.integer(100 * quantiles)), |
36 | 1x |
stat_names[seq_along(stat_names) > quantile_index]) |
37 |
} |
|
38 | ||
39 | 84x |
ind <- .valueWeightIndexes(length(value_names), length(weight_names)) |
40 | 84x |
vn <- value_names[ind$values] |
41 | 84x |
wn <- weight_names[ind$weights] |
42 | ||
43 | 84x |
if (length(vn) > 1 || full_colnames) { |
44 |
# determine all combinations of index and stat |
|
45 | 33x |
z <- expand.grid(index=seq_along(vn), |
46 | 33x |
stat=stat_names, stringsAsFactors=FALSE) |
47 | 33x |
z$value <- vn[z$index] |
48 | 33x |
if (is.null(wn)) { |
49 | 14x |
z$weights <- NA |
50 |
} else { |
|
51 | 19x |
z$weights <- wn[z$index] |
52 |
} |
|
53 | ||
54 |
# construct column names for each index, stat |
|
55 |
# add weight layer name only if layer is ambiguously weighted |
|
56 | 33x |
mapply(function(stat, value, weight) { |
57 | 123x |
ret <- stat |
58 | 123x |
if (full_colnames || length(value_names) > 1) { |
59 | 108x |
ret <- paste(ret, value, sep='.') |
60 |
} |
|
61 | 123x |
if (.includeWeightInColName(stat) && ((full_colnames & length(weight_names) > 0) |
62 | 123x |
|| length(weight_names) > 1)) { |
63 | 45x |
ret <- paste(ret, weight, sep='.') |
64 |
} |
|
65 | ||
66 | 123x |
return(ret) |
67 | 33x |
}, z$stat, z$value, z$weight, USE.NAMES = FALSE) |
68 |
} else { |
|
69 | 51x |
stat_names |
70 |
} |
|
71 |
} |
|
72 | ||
73 |
.includeWeightInColName <- function(fun) { |
|
74 | 123x |
.isWeighted(fun) | fun == 'fun' |
75 |
} |
|
76 | ||
77 |
.isWeighted <- function(stat_name) { |
|
78 | 193x |
stat_name %in% c('weighted_mean', 'weighted_sum') |
79 |
} |
|
80 | ||
81 |
#' Compute indexes for the value and weight layers that should be |
|
82 |
#' processed together |
|
83 |
#' |
|
84 |
#' @param num_values number of layers in value raster |
|
85 |
#' @param num_weights number of layers in weighting raster |
|
86 |
#' @return list with \code{values} and \code{weights} elements |
|
87 |
#' providing layer indexes |
|
88 |
.valueWeightIndexes <- function(num_values, num_weights) { |
|
89 | 128x |
if (num_weights == 0) { |
90 | 85x |
vi <- seq_len(num_values) |
91 | 85x |
wi <- NA |
92 | 43x |
} else if (num_values == num_weights) { |
93 |
# process in parallel |
|
94 | 28x |
vi <- seq_len(num_values) |
95 | 28x |
wi <- seq_len(num_weights) |
96 | 15x |
} else if (num_values == 1 && num_weights > 1) { |
97 |
# recycle values |
|
98 | 7x |
vi <- rep.int(1, num_weights) |
99 | 7x |
wi <- seq_len(num_weights) |
100 | 8x |
} else if (num_values > 1 && num_weights == 1) { |
101 |
# recycle weights |
|
102 | 8x |
vi <- seq_len(num_values) |
103 | 8x |
wi <- rep.int(1, num_values) |
104 |
} |
|
105 | ||
106 | 128x |
list(values = vi, weights = wi) |
107 |
} |
|
108 |
1 |
# Copyright (c) 2020 ISciences, LLC. |
|
2 |
# All rights reserved. |
|
3 |
# |
|
4 |
# This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
# You may not use this file except in compliance with the License. You may |
|
6 |
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
# |
|
8 |
# Unless required by applicable law or agreed to in writing, software |
|
9 |
# distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
# See the License for the specific language governing permissions and |
|
12 |
# limitations under the License. |
|
13 | ||
14 |
if (!isGeneric("exact_resample")) { |
|
15 | 7x |
setGeneric("exact_resample", function(x, y, ...) |
16 | 7x |
standardGeneric("exact_resample")) |
17 |
} |
|
18 | ||
19 |
#' Resample a raster to a new grid |
|
20 |
#' |
|
21 |
#' @param x a \code{RasterLayer} to be resampled |
|
22 |
#' @param y a \code{RasterLayer} with a grid definition to which \code{x} |
|
23 |
#' should be resampled |
|
24 |
#' @param fun a named summary operation to be used for the resampling |
|
25 |
#' @return a resampled version of \code{x} |
|
26 |
#' |
|
27 |
#' @name exact_resample |
|
28 |
NULL |
|
29 | ||
30 |
#' @import sf |
|
31 |
#' @import raster |
|
32 |
#' @useDynLib exactextractr |
|
33 |
#' @rdname exact_resample |
|
34 |
#' @export |
|
35 |
setMethod('exact_resample', |
|
36 |
signature(x='RasterLayer', y='RasterLayer'), |
|
37 |
function(x, y, fun) { |
|
38 | 7x |
if (sf::st_crs(x) != sf::st_crs(y)) { |
39 | 3x |
if (is.na(sf::st_crs(x))) { |
40 | 1x |
warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.") |
41 | 2x |
} else if (is.na(sf::st_crs(y))) { |
42 | 1x |
warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.") |
43 |
} else { |
|
44 | 1x |
stop('Destination raster must have same CRS as source.') |
45 |
} |
|
46 |
} |
|
47 | ||
48 | 6x |
x <- raster::readStart(x) |
49 | 6x |
tryCatch({ |
50 | 6x |
CPP_resample(x, y, fun) |
51 | 6x |
}, finally={ |
52 | 6x |
raster::readStop(x) |
53 |
}) |
|
54 |
}) |
1 |
# Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
# All rights reserved. |
|
3 |
# |
|
4 |
# This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
# You may not use this file except in compliance with the License. You may |
|
6 |
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
# |
|
8 |
# Unless required by applicable law or agreed to in writing, software |
|
9 |
# distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
# See the License for the specific language governing permissions and |
|
12 |
# limitations under the License. |
|
13 | ||
14 |
if (!isGeneric('coverage_fraction')) { |
|
15 | 12x |
setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...) |
16 | 12x |
standardGeneric('coverage_fraction')) |
17 |
} |
|
18 | ||
19 |
.coverage_fraction <- function(x, y, crop) { |
|
20 | 12x |
if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) { |
21 | 1x |
warning("No CRS specified for raster; assuming it has the same CRS as the polygons.") |
22 | 11x |
} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) { |
23 | 1x |
warning("No CRS specified for polygons; assuming they have the same CRS as the raster.") |
24 | 10x |
} else if(sf::st_crs(x) != sf::st_crs(y)) { |
25 | 1x |
y <- sf::st_transform(y, sf::st_crs(x)) |
26 | 1x |
warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")") |
27 |
} |
|
28 | ||
29 | 12x |
lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) { |
30 | 12x |
CPP_coverage_fraction(x, wkb, crop) |
31 |
}) |
|
32 |
} |
|
33 | ||
34 |
#' Compute the fraction of raster cells covered by a polygon |
|
35 |
#' |
|
36 |
#' @param x a (possibly empty) \code{RasterLayer} whose resolution and |
|
37 |
#' extent will be used for the generated \code{RasterLayer}. |
|
38 |
#' @param y a \code{sf} object with polygonal geometries |
|
39 |
#' @param crop if \code{TRUE}, each generated \code{RasterLayer} will be |
|
40 |
#' cropped to the extent of its associated feature. |
|
41 |
#' @return a list with a \code{RasterLayer} for each feature in \code{y}. |
|
42 |
#' Values of the raster represent the fraction of each |
|
43 |
#' cell in \code{x} that is covered by \code{y}. |
|
44 |
#' @examples |
|
45 |
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10) |
|
46 |
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))') |
|
47 |
#' |
|
48 |
#' cov_frac <- coverage_fraction(rast, poly)[[1]] |
|
49 |
#' @name coverage_fraction |
|
50 |
NULL |
|
51 | ||
52 |
#' @import sf |
|
53 |
#' @import raster |
|
54 |
#' @useDynLib exactextractr |
|
55 |
#' @rdname coverage_fraction |
|
56 |
#' @export |
|
57 |
setMethod('coverage_fraction', signature(x='RasterLayer', y='sf'), function(x, y, crop=FALSE) { |
|
58 | ! |
coverage_fraction(x, sf::st_geometry(y), crop) |
59 |
}) |
|
60 | ||
61 |
#' @rdname coverage_fraction |
|
62 |
#' @export |
|
63 |
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_MULTIPOLYGON'), .coverage_fraction) |
|
64 | ||
65 |
#' @rdname coverage_fraction |
|
66 |
#' @export |
|
67 |
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_POLYGON'), .coverage_fraction) |
|
68 |
1 |
// Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#include <Rcpp.h> |
|
15 | ||
16 |
#include "geos_r.h" |
|
17 | ||
18 |
#include "raster_utils.h" |
|
19 | ||
20 |
#include "exactextract/src/raster.h" |
|
21 |
#include "exactextract/src/raster_cell_intersection.h" |
|
22 | ||
23 |
using exactextract::RasterView; |
|
24 |
using exactextract::raster_cell_intersection; |
|
25 | ||
26 |
// [[Rcpp::export]] |
|
27 | 12x |
Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, const Rcpp::RawVector & wkb, bool crop) |
28 |
{ |
|
29 | 24x |
GEOSAutoHandle geos; |
30 | 36x |
Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster"); |
31 | 36x |
Rcpp::Function rasterFn = raster["raster"]; |
32 | 36x |
Rcpp::Function crsFn = raster["crs"]; |
33 | ||
34 | 12x |
auto grid = make_grid(rast); |
35 | 24x |
auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get()); |
36 | ||
37 | 12x |
if (crop) { |
38 | 2x |
grid = coverage_fraction.grid(); |
39 |
} |
|
40 | 24x |
RasterView<float> coverage_view(coverage_fraction, grid); |
41 | ||
42 | 12x |
Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()), |
43 | 12x |
static_cast<int>(grid.cols())}; |
44 | ||
45 | 629x |
for (size_t i = 0; i < grid.rows(); i++) { |
46 | 159018x |
for (size_t j = 0; j < grid.cols(); j++) { |
47 | 158401x |
weights(i, j) = coverage_view(i, j); |
48 | 158401x |
if (!crop && std::isnan(weights(i, j))) { |
49 | 139692x |
weights(i, j) = 0; |
50 |
} |
|
51 |
} |
|
52 |
} |
|
53 | ||
54 |
return rasterFn(weights, |
|
55 | 24x |
Rcpp::Named("xmn")=grid.xmin(), |
56 | 24x |
Rcpp::Named("xmx")=grid.xmax(), |
57 | 24x |
Rcpp::Named("ymn")=grid.ymin(), |
58 | 24x |
Rcpp::Named("ymx")=grid.ymax(), |
59 | 48x |
Rcpp::Named("crs")=crsFn(rast)); |
60 |
} |
1 |
// Copyright (c) 2018 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_GRID_H |
|
15 |
#define EXACTEXTRACT_GRID_H |
|
16 | ||
17 |
#include <numeric> |
|
18 |
#include <stdexcept> |
|
19 |
#include <vector> |
|
20 | ||
21 |
#include "box.h" |
|
22 | ||
23 |
namespace exactextract { |
|
24 |
struct infinite_extent { |
|
25 |
static const size_t padding = 1; |
|
26 |
}; |
|
27 | ||
28 |
struct bounded_extent { |
|
29 |
static const size_t padding = 0; |
|
30 |
}; |
|
31 | ||
32 | 239x |
static inline bool is_integral(double d, double tol) { |
33 | 239x |
return std::abs(d - std::round(d)) <= tol; |
34 |
} |
|
35 | ||
36 |
template<typename extent_tag> |
|
37 |
class Grid { |
|
38 | ||
39 |
public: |
|
40 | 1833890x |
Grid(const Box & extent, double dx, double dy) : |
41 |
m_extent{extent}, |
|
42 |
m_dx{dx}, |
|
43 |
m_dy{dy}, |
|
44 | 1833890x |
m_num_rows{2*extent_tag::padding + (extent.ymax > extent.ymin ? static_cast<size_t>(std::round((extent.ymax - extent.ymin) / dy)) : 0)}, |
45 | 3667780x |
m_num_cols{2*extent_tag::padding + (extent.xmax > extent.xmin ? static_cast<size_t>(std::round((extent.xmax - extent.xmin) / dx)) : 0)} |
46 |
{} |
|
47 | ||
48 | 118119x |
static Grid make_empty() { |
49 | 118119x |
return Grid({0, 0, 0, 0}, 0, 0); |
50 |
} |
|
51 | ||
52 | 4789762x |
size_t get_column(double x) const { |
53 |
if (extent_tag::padding) { |
|
54 | 6380700x |
if (x < m_extent.xmin) |
55 | 858x |
return 0; |
56 | 6379842x |
if (x > m_extent.xmax) |
57 | 1122x |
return m_num_cols - 1; |
58 | 6378720x |
if (x == m_extent.xmax) // special case, returning the cell for which xmax is the right |
59 | 3190154x |
return m_num_cols - 2; |
60 |
} else { |
|
61 | 1599412x |
if (x < m_extent.xmin || x > m_extent.xmax) |
62 | ! |
throw std::out_of_range("x"); |
63 | ||
64 | 1599412x |
if (x == m_extent.xmax) |
65 | 800356x |
return m_num_cols - 1; |
66 |
} |
|
67 | ||
68 |
// Since we've already range-checked our x value, make sure that |
|
69 |
// the computed column index is no greater than the column index |
|
70 |
// associated with xmax. |
|
71 | 7180017x |
return std::min( |
72 | 2393339x |
extent_tag::padding + static_cast<size_t>(std::floor((x - m_extent.xmin) / m_dx)), |
73 | 4786678x |
get_column(m_extent.xmax)); |
74 |
} |
|
75 | ||
76 | 4788530x |
size_t get_row(double y) const { |
77 |
if (extent_tag::padding) { |
|
78 | 6378008x |
if (y > m_extent.ymax) |
79 | 2206x |
return 0; |
80 | 6375802x |
if (y < m_extent.ymin) |
81 | 1848x |
return m_num_rows - 1; |
82 | 6373954x |
if (y == m_extent.ymin) // special case, returning the cell for which ymin is the bottom |
83 | 3188080x |
return m_num_rows - 2; |
84 |
} else { |
|
85 | 1599526x |
if (y < m_extent.ymin || y > m_extent.ymax) |
86 | ! |
throw std::out_of_range("y"); |
87 | ||
88 | 1599526x |
if (y == m_extent.ymin) |
89 | 800356x |
return m_num_rows - 1; |
90 |
} |
|
91 | ||
92 |
// Since we've already range-checked our y value, make sure that |
|
93 |
// the computed row index is no greater than the column index |
|
94 |
// associated with ymin. |
|
95 | 7176321x |
return std::min( |
96 | 2392107x |
extent_tag::padding + static_cast<size_t>(std::floor((m_extent.ymax - y) / m_dy)), |
97 | 4784214x |
get_row(m_extent.ymin)); |
98 |
} |
|
99 | ||
100 | 1035458x |
bool empty() const { return m_num_rows <= 2*extent_tag::padding && m_num_cols <= 2*extent_tag::padding; } |
101 | ||
102 | 10621282x |
size_t rows() const { return m_num_rows; } |
103 | ||
104 | 12042395x |
size_t cols() const { return m_num_cols; } |
105 | ||
106 | 131x |
size_t size() const { return rows()*cols(); } |
107 | ||
108 | 5314174x |
double xmin() const { return m_extent.xmin; } |
109 | ||
110 | 3808835x |
double xmax() const { return m_extent.xmax; } |
111 | ||
112 | 3667155x |
double ymin() const { return m_extent.ymin; } |
113 | ||
114 | 5549438x |
double ymax() const { return m_extent.ymax; } |
115 | ||
116 | 5830049x |
double dx() const { return m_dx; } |
117 | ||
118 | 6067104x |
double dy() const { return m_dy; } |
119 | ||
120 | 2233705x |
const Box& extent() const { return m_extent; } |
121 | ||
122 | 399922x |
size_t row_offset(const Grid & other) const { return static_cast<size_t>(std::round(std::abs(other.m_extent.ymax - m_extent.ymax) / m_dy)); } |
123 | ||
124 | 399922x |
size_t col_offset(const Grid & other) const { return static_cast<size_t>(std::round(std::abs(m_extent.xmin - other.m_extent.xmin) / m_dx)); } |
125 | ||
126 | 283x |
double x_for_col(size_t col) const { return m_extent.xmin + ((col - extent_tag::padding) + 0.5) * m_dx; } |
127 | ||
128 | 898x |
double y_for_row(size_t row) const { return m_extent.ymax - ((row - extent_tag::padding) + 0.5) * m_dy; } |
129 | ||
130 | 1192x |
Grid<extent_tag> crop(const Box & b) const { |
131 | 1192x |
if (extent().intersects(b)) { |
132 | 1109x |
return shrink_to_fit(b.intersection(extent())); |
133 |
} else { |
|
134 | 83x |
return make_empty(); |
135 |
} |
|
136 |
} |
|
137 | ||
138 | 400178x |
Grid<extent_tag> shrink_to_fit(const Box & b) const { |
139 | 400178x |
if (b.xmin < m_extent.xmin || b.ymin < m_extent.ymin || b.xmax > m_extent.xmax || b.ymax > m_extent.ymax) { |
140 | ! |
throw std::range_error("Cannot shrink extent to bounds larger than original."); |
141 |
} |
|
142 | ||
143 | 400178x |
size_t col0 = get_column(b.xmin); |
144 | 400178x |
size_t row1 = get_row(b.ymax); |
145 | ||
146 |
// Shrink xmin and ymax to fit the upper-left corner of the supplied extent |
|
147 | 400178x |
double snapped_xmin = m_extent.xmin + (col0 - extent_tag::padding) * m_dx; |
148 | 400178x |
double snapped_ymax = m_extent.ymax - (row1 - extent_tag::padding) * m_dy; |
149 | ||
150 |
// Make sure x0 and y1 are within the reduced extent. Because of |
|
151 |
// floating point round-off errors, this is not always the case. |
|
152 | 400178x |
if (b.xmin < snapped_xmin) { |
153 | ! |
snapped_xmin -= m_dx; |
154 | ! |
col0--; |
155 |
} |
|
156 | ||
157 | 400178x |
if (b.ymax > snapped_ymax) { |
158 | ! |
snapped_ymax += m_dy; |
159 | ! |
row1--; |
160 |
} |
|
161 | ||
162 | 400178x |
size_t col1 = get_column(b.xmax); |
163 | 400178x |
size_t row0 = get_row(b.ymin); |
164 | ||
165 | 400178x |
size_t num_rows = 1 + (row0 - row1); |
166 | 400178x |
size_t num_cols = 1 + (col1 - col0); |
167 | ||
168 |
// If xmax or ymin falls cleanly on a cell boundary, we don't |
|
169 |
// need as many rows or columns as we otherwise would, because |
|
170 |
// we assume that the rightmost cell of the grid is a closed |
|
171 |
// interval in x, and the lowermost cell of the grid is a |
|
172 |
// closed interval in y. |
|
173 | 400178x |
if (num_rows > 2 && (snapped_ymax - (num_rows-1)*m_dy <= b.ymin)) { |
174 | 370x |
num_rows--; |
175 |
} |
|
176 | 400178x |
if (num_cols > 2 && (snapped_xmin + (num_cols-1)*m_dx >= b.xmax)) { |
177 | 315x |
num_cols--; |
178 |
} |
|
179 | ||
180 |
// Perform offsets relative to the new xmin, ymax origin |
|
181 |
// points. If this is not done, then floating point roundoff |
|
182 |
// error can cause progressive shrink() calls with the same |
|
183 |
// inputs to produce different results. |
|
184 | 400178x |
Box reduced_box = { |
185 |
snapped_xmin, |
|
186 | 400178x |
std::min(snapped_ymax - num_rows * m_dy, b.ymin), |
187 | 400178x |
std::max(snapped_xmin + num_cols * m_dx, b.xmax), |
188 |
snapped_ymax |
|
189 |
}; |
|
190 | ||
191 |
// Fudge computed xmax and ymin, if needed, to prevent extent |
|
192 |
// from growing during a shrink operation. |
|
193 | 400178x |
if (reduced_box.xmax > m_extent.xmax) { |
194 | ! |
if (std::round((reduced_box.xmax - reduced_box.xmin)/m_dx) == |
195 | ! |
std::round((m_extent.xmax - reduced_box.xmin)/m_dx)) { |
196 | ! |
reduced_box.xmax = m_extent.xmax; |
197 |
} else { |
|
198 | ! |
throw std::runtime_error("Shrink operation failed."); |
199 |
} |
|
200 |
} |
|
201 | 400178x |
if (reduced_box.ymin < m_extent.ymin) { |
202 | ! |
if (std::round((reduced_box.ymax - reduced_box.ymin)/m_dy) == |
203 | ! |
std::round((reduced_box.ymax - m_extent.ymin)/m_dy)) { |
204 | ! |
reduced_box.ymin = m_extent.ymin; |
205 |
} else { |
|
206 | ! |
throw std::runtime_error("Shrink operation failed."); |
207 |
} |
|
208 |
} |
|
209 | ||
210 | 400178x |
Grid<extent_tag> reduced{reduced_box, m_dx, m_dy}; |
211 | ||
212 | ! |
if (b.xmin < reduced.xmin() || b.ymin < reduced.ymin() || b.xmax > reduced.xmax() || b.ymax > reduced.ymax()) { |
213 | ! |
throw std::runtime_error("Shrink operation failed."); |
214 |
} |
|
215 | ||
216 | 800356x |
return reduced; |
217 |
} |
|
218 | ||
219 |
template<typename extent_tag2> |
|
220 | 60x |
bool compatible_with(const Grid<extent_tag2> &b) const { |
221 |
// Define a tolerance for grid compatibility, to be used in the following ways: |
|
222 |
// Grid resolutions must be integer multiples of each other within a factor of 'compatibility_tol' |
|
223 |
// Grid origin points must differ by an integer multiple of the smaller of the two grid resolutions, |
|
224 |
// within a factor of 'compatibility_tol'. |
|
225 |
// Perhaps it's possible to remove this hardcoded tolerance and express something in terms of |
|
226 |
// std::numeric_limits<double>::epsilon(), but I wasn't able to make that work. |
|
227 | 60x |
constexpr double compatability_tol = 1e-6; |
228 | ||
229 |
if (empty() || b.empty()) { |
|
230 | ! |
return true; |
231 |
} |
|
232 | ||
233 |
// Check x-resolution compatibility |
|
234 | 60x |
if (!is_integral(std::max(m_dx, b.m_dx) / std::min(m_dx, b.m_dx), std::min(m_dx, b.m_dx)*compatability_tol)) { |
235 | ! |
return false; |
236 |
} |
|
237 | ||
238 |
// Check y-resolution compatibility |
|
239 | 60x |
if (!is_integral(std::max(m_dy, b.m_dy) / std::min(m_dy, b.m_dy), std::min(m_dy, b.m_dy)*compatability_tol)) { |
240 | ! |
return false; |
241 |
} |
|
242 | ||
243 |
// Check left-hand boundary compatibility |
|
244 | 60x |
if (!is_integral(std::abs(b.m_extent.xmin - m_extent.xmin) / std::min(m_dx, b.m_dx), std::min(m_dx, b.m_dx)*compatability_tol)) { |
245 | 1x |
return false; |
246 |
} |
|
247 | ||
248 |
// Check upper boundary compatibility |
|
249 | 59x |
if (!is_integral(std::abs(b.m_extent.ymax - m_extent.ymax) / std::min(m_dy, b.m_dy), std::min(m_dy, b.m_dy)*compatability_tol)) { |
250 | ! |
return false; |
251 |
} |
|
252 | ||
253 | 59x |
return true; |
254 |
} |
|
255 | ||
256 |
template<typename extent_tag2> |
|
257 | 60x |
Grid<extent_tag> common_grid(const Grid<extent_tag2> &b) const { |
258 | 60x |
if (!compatible_with(b)) { |
259 | 1x |
throw std::runtime_error("Incompatible extents."); |
260 |
} |
|
261 | ||
262 | 59x |
if (b.empty()) { |
263 | ! |
return *this; |
264 |
} |
|
265 | ||
266 | 59x |
const double common_dx = std::min(m_dx, b.m_dx); |
267 | 59x |
const double common_dy = std::min(m_dy, b.m_dy); |
268 | ||
269 | 59x |
const double common_xmin = std::min(m_extent.xmin, b.m_extent.xmin); |
270 | 59x |
const double common_ymax = std::max(m_extent.ymax, b.m_extent.ymax); |
271 | ||
272 | 59x |
double common_xmax = std::max(m_extent.xmax, b.m_extent.xmax); |
273 | 59x |
double common_ymin = std::min(m_extent.ymin, b.m_extent.ymin); |
274 | ||
275 | 59x |
const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx)); |
276 | 59x |
const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy)); |
277 | ||
278 | 59x |
common_xmax = std::max(common_xmax, common_xmin + nx*common_dx); |
279 | 59x |
common_ymin = std::min(common_ymin, common_ymax - ny*common_dy); |
280 | ||
281 | 59x |
return {{ common_xmin, common_ymin, common_xmax, common_ymax}, common_dx, common_dy }; |
282 |
} |
|
283 | ||
284 |
template<typename extent_tag2> |
|
285 |
Grid<extent_tag> overlapping_grid(const Grid<extent_tag2> &b) const { |
|
286 |
if (!compatible_with(b)) { |
|
287 |
throw std::runtime_error("Incompatible extents."); |
|
288 |
} |
|
289 | ||
290 |
if (empty() || b.empty()) { |
|
291 |
return make_empty(); |
|
292 |
} |
|
293 | ||
294 |
const double common_dx = std::min(m_dx, b.m_dx); |
|
295 |
const double common_dy = std::min(m_dy, b.m_dy); |
|
296 | ||
297 |
const double common_xmin = std::max(m_extent.xmin, b.m_extent.xmin); |
|
298 |
const double common_ymax = std::min(m_extent.ymax, b.m_extent.ymax); |
|
299 | ||
300 |
double common_xmax = std::min(m_extent.xmax, b.m_extent.xmax); |
|
301 |
double common_ymin = std::max(m_extent.ymin, b.m_extent.ymin); |
|
302 | ||
303 |
const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx)); |
|
304 |
const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy)); |
|
305 | ||
306 |
common_xmax = std::max(common_xmax, common_xmin + nx*common_dx); |
|
307 |
common_ymin = std::min(common_ymin, common_ymax - ny*common_dy); |
|
308 | ||
309 |
return {{ common_xmin, common_ymin, common_xmax, common_ymax}, common_dx, common_dy }; |
|
310 |
} |
|
311 | ||
312 |
bool operator==(const Grid<extent_tag> &b) const { |
|
313 |
return |
|
314 |
m_extent == b.m_extent && |
|
315 |
m_dx == b.m_dx && |
|
316 |
m_dy == b.m_dy; |
|
317 |
} |
|
318 | ||
319 |
bool operator!=(const Grid<extent_tag> &b) const { |
|
320 |
return !(*this == b); |
|
321 |
} |
|
322 | ||
323 |
private: |
|
324 |
Box m_extent; |
|
325 | ||
326 |
double m_dx; |
|
327 |
double m_dy; |
|
328 | ||
329 |
size_t m_num_rows; |
|
330 |
size_t m_num_cols; |
|
331 |
}; |
|
332 | ||
333 |
Box grid_cell(const Grid<bounded_extent> & grid, size_t row, size_t col); |
|
334 |
Box grid_cell(const Grid<infinite_extent> & grid, size_t row, size_t col); |
|
335 | ||
336 |
Grid<infinite_extent> make_infinite(const Grid<bounded_extent> & grid); |
|
337 |
Grid<bounded_extent> make_finite(const Grid<infinite_extent> & grid); |
|
338 | ||
339 |
std::vector<Grid<bounded_extent>> subdivide(const Grid<bounded_extent> & grid, size_t max_size); |
|
340 | ||
341 |
template<typename T> |
|
342 |
Grid<bounded_extent> common_grid(T begin, T end) { |
|
343 |
if (begin == end) { |
|
344 |
return Grid<bounded_extent>::make_empty(); |
|
345 |
} else if (std::next(begin) == end) { |
|
346 |
return begin->grid(); |
|
347 |
} |
|
348 |
return std::accumulate( |
|
349 |
std::next(begin), |
|
350 |
end, |
|
351 |
begin->grid(), |
|
352 |
[](auto& acc, auto& op) { |
|
353 |
return acc.common_grid(op.grid()); |
|
354 |
}); |
|
355 |
} |
|
356 | ||
357 |
} |
|
358 | ||
359 |
#endif //EXACTEXTRACT_INFINITEGRID_H |
1 |
// Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_MATRIX_H |
|
15 |
#define EXACTEXTRACT_MATRIX_H |
|
16 | ||
17 |
#include <iomanip> |
|
18 |
#include <iterator> |
|
19 |
#include <memory> |
|
20 |
#include <cstring> |
|
21 |
#include <vector> |
|
22 | ||
23 |
namespace exactextract { |
|
24 | ||
25 |
template<typename T> |
|
26 |
class Matrix { |
|
27 | ||
28 |
public: |
|
29 |
using value_type = T; |
|
30 | ||
31 | 1831570x |
Matrix(size_t rows, size_t cols) : |
32 |
m_rows{rows}, |
|
33 | 1831570x |
m_cols{cols} |
34 |
{ |
|
35 | 1831570x |
if (m_rows > 0 && m_cols > 0) { |
36 |
// new T[]() initializes to zero |
|
37 |
m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]()); |
|
38 |
} |
|
39 |
} |
|
40 | ||
41 | 213x |
Matrix(size_t rows, size_t cols, T value) : |
42 |
m_rows{rows}, |
|
43 | 213x |
m_cols{cols} |
44 |
{ |
|
45 | 213x |
if (m_rows > 0 && m_cols > 0) { |
46 |
// new T[] does not initialize |
|
47 |
m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]); |
|
48 |
} |
|
49 | ||
50 | 213x |
std::fill(m_data.get(), m_data.get() + m_rows*m_cols, value); |
51 |
} |
|
52 | ||
53 |
explicit Matrix(const std::vector<std::vector<T>> & data) : |
|
54 |
m_rows{data.size()}, |
|
55 |
m_cols{data[0].size()} |
|
56 |
{ |
|
57 |
m_data = std::unique_ptr<T[]>(new T[m_rows*m_cols]()); |
|
58 | ||
59 |
auto lastpos = m_data.get(); |
|
60 |
for (auto& row : data) { |
|
61 |
lastpos = std::copy(row.begin(), row.end(), lastpos); |
|
62 |
} |
|
63 |
} |
|
64 | ||
65 | 516717x |
Matrix(Matrix<T>&& m) noexcept : |
66 |
m_rows{m.rows()}, |
|
67 | 516717x |
m_cols{m.cols()} |
68 |
{ |
|
69 | 516717x |
m_data = std::move(m.m_data); |
70 |
} |
|
71 | ||
72 | 3647870x |
T& operator()(size_t row, size_t col) { |
73 | 3647870x |
check(row, col); |
74 | 3647870x |
return m_data[row*m_cols + col]; |
75 |
} |
|
76 | ||
77 | 1247057x |
T operator()(size_t row, size_t col) const { |
78 | 1247057x |
check(row, col); |
79 | 1247057x |
return m_data[row*m_cols + col]; |
80 |
} |
|
81 | ||
82 |
bool operator==(const Matrix<T> & other) const { |
|
83 |
if (m_rows != other.m_rows) { |
|
84 |
return false; |
|
85 |
} |
|
86 |
if (m_cols != other.m_cols) { |
|
87 |
return false; |
|
88 |
} |
|
89 | ||
90 |
return 0 == memcmp(m_data.get(), other.m_data.get(), m_rows*m_cols*sizeof(T)); |
|
91 |
} |
|
92 | ||
93 | 623378x |
void increment(size_t row, size_t col, const T & val) { |
94 | 623378x |
check(row, col); |
95 | 623378x |
m_data[row*m_cols + col] += val; |
96 |
} |
|
97 | ||
98 | 1436188x |
size_t rows() const { return m_rows; } |
99 | 1727862x |
size_t cols() const { return m_cols; } |
100 | ||
101 |
T* row(size_t row) { |
|
102 |
return &(m_data[row*m_cols]); |
|
103 |
} |
|
104 | ||
105 |
T* data() { |
|
106 |
return m_data.get(); |
|
107 |
} |
|
108 | ||
109 |
#ifdef MATRIX_CHECK_BOUNDS |
|
110 |
void check(size_t row, size_t col) const { |
|
111 |
if (row + 1 > m_rows) { |
|
112 |
throw std::out_of_range("Row " + std::to_string(row) + " is out of range."); |
|
113 |
} |
|
114 |
if (col + 1 > m_cols) { |
|
115 |
throw std::out_of_range("Col " + std::to_string(col) + " is out of range."); |
|
116 |
} |
|
117 |
} |
|
118 |
#else |
|
119 | 3694370x |
void check(size_t, size_t) const {} |
120 |
#endif |
|
121 | ||
122 |
private: |
|
123 |
std::unique_ptr<T[]> m_data; |
|
124 | ||
125 |
size_t m_rows; |
|
126 |
size_t m_cols; |
|
127 | ||
128 |
}; |
|
129 | ||
130 |
template<typename T> |
|
131 |
std::ostream& operator<<(std::ostream & os, const Matrix<T> & m) { |
|
132 |
for (size_t i = 0; i < m.rows(); i++) { |
|
133 |
for (size_t j = 0; j < m.cols(); j++) { |
|
134 |
if (m(i, j) != 0) { |
|
135 |
os << std::right << std::fixed << std::setw(10) << std::setprecision(6) << |
|
136 |
m(i, j) << " "; |
|
137 |
} else { |
|
138 |
os << " "; |
|
139 |
} |
|
140 |
} |
|
141 |
os << std::endl; |
|
142 |
} |
|
143 | ||
144 |
return os; |
|
145 |
} |
|
146 | ||
147 |
} |
|
148 | ||
149 |
#endif |
1 |
// Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_RASTER_H |
|
15 |
#define EXACTEXTRACT_RASTER_H |
|
16 | ||
17 |
#include <cmath> |
|
18 |
#include <limits> |
|
19 | ||
20 |
#include "grid.h" |
|
21 |
#include "matrix.h" |
|
22 | ||
23 |
namespace exactextract { |
|
24 |
template<typename T> |
|
25 |
class AbstractRaster { |
|
26 |
public: |
|
27 |
using value_type = T; |
|
28 | ||
29 | 1833604x |
explicit AbstractRaster(const Grid<bounded_extent> & ex) : |
30 |
m_grid{ex}, |
|
31 | 1833604x |
m_nodata{std::is_floating_point<T>::value ? std::numeric_limits<T>::quiet_NaN() : std::numeric_limits<T>::min()}, |
32 | 1833604x |
m_has_nodata{false} |
33 |
{} |
|
34 | ||
35 |
AbstractRaster(const Grid<bounded_extent> & ex, const T& nodata_val) : m_grid{ex}, m_nodata{nodata_val}, m_has_nodata{true} {} |
|
36 | ||
37 | 916802x |
virtual ~AbstractRaster() = default; |
38 | ||
39 | 3238180x |
size_t rows() const { |
40 | 3238180x |
return m_grid.rows(); |
41 |
} |
|
42 | ||
43 | 4735006x |
size_t cols() const { |
44 | 4735006x |
return m_grid.cols(); |
45 |
} |
|
46 | ||
47 | 798226x |
double xres() const { |
48 | 798226x |
return m_grid.dx(); |
49 |
} |
|
50 | ||
51 | 798226x |
double yres() const { |
52 | 798226x |
return m_grid.dy(); |
53 |
} |
|
54 | ||
55 | 798226x |
double xmin() const { |
56 | 798226x |
return m_grid.xmin(); |
57 |
} |
|
58 | ||
59 |
double ymin() const { |
|
60 |
return m_grid.ymin(); |
|
61 |
} |
|
62 | ||
63 |
double xmax() const { |
|
64 |
return m_grid.xmax(); |
|
65 |
} |
|
66 | ||
67 | 798226x |
double ymax() const { |
68 | 798226x |
return m_grid.ymax(); |
69 |
} |
|
70 | ||
71 | 915692x |
const Grid<bounded_extent>& grid() const { |
72 | 915692x |
return m_grid; |
73 |
} |
|
74 | ||
75 |
virtual T operator()(size_t row, size_t col) const = 0; |
|
76 | ||
77 | 798226x |
bool has_nodata() const { return m_has_nodata; } |
78 | ||
79 | 279412x |
T nodata() const { return m_nodata; } |
80 | ||
81 | 605057x |
bool get(size_t row, size_t col, T & val) { |
82 | 605057x |
val = operator()(row, col); |
83 | ||
84 | ! |
if (m_has_nodata && val == m_nodata) { |
85 | ! |
return false; |
86 |
} |
|
87 | 605057x |
if (std::is_floating_point<T>::value && std::isnan(val)) { |
88 | 23x |
return false; |
89 |
} |
|
90 | ||
91 | 605034x |
return true; |
92 |
} |
|
93 | ||
94 | ! |
void set_nodata(const T & val) { |
95 | ! |
m_has_nodata = true; |
96 | ! |
m_nodata = val; |
97 |
} |
|
98 | ||
99 |
bool operator==(const AbstractRaster<T> & other) const { |
|
100 |
if (rows() != other.rows()) |
|
101 |
return false; |
|
102 | ||
103 |
if (cols() != other.cols()) |
|
104 |
return false; |
|
105 | ||
106 |
if (xres() != other.xres()) |
|
107 |
return false; |
|
108 | ||
109 |
if (yres() != other.yres()) |
|
110 |
return false; |
|
111 | ||
112 |
if (xmin() != other.xmin()) |
|
113 |
return false; |
|
114 | ||
115 |
if (ymin() != other.ymin()) |
|
116 |
return false; |
|
117 | ||
118 |
// Do the rasters differ in their definition of NODATA? If so, we need to do some more detailed |
|
119 |
// checking below. |
|
120 |
bool nodata_differs = has_nodata() != other.has_nodata() || |
|
121 |
(has_nodata() && other.has_nodata() && nodata() != other.nodata()); |
|
122 | ||
123 |
for (size_t i = 0; i < rows(); i++) { |
|
124 |
for (size_t j = 0; j < cols(); j++) { |
|
125 |
if(operator()(i, j) != other(i, j)) { |
|
126 |
// Override default behavior of NAN != NAN |
|
127 |
if (!std::isnan(operator()(i, j)) || !std::isnan(other(i, j))) { |
|
128 |
return false; |
|
129 |
} |
|
130 |
} else if (nodata_differs && (operator()(i, j) == nodata() || other(i, j) == other.nodata())) { |
|
131 |
// For data types that do not have NAN, or even for floating point types where the user has |
|
132 |
// selected some other value to represent NODATA, we need to reverse a positive equality test |
|
133 |
// where the value is considered to be NODATA in one raster but not the other. |
|
134 |
return false; |
|
135 |
} |
|
136 |
} |
|
137 |
} |
|
138 | ||
139 |
return true; |
|
140 |
} |
|
141 | ||
142 |
bool operator!=(const AbstractRaster<T> & other) const { |
|
143 |
return !(operator==(other)); |
|
144 |
} |
|
145 |
private: |
|
146 |
Grid<bounded_extent> m_grid; |
|
147 |
T m_nodata; |
|
148 |
bool m_has_nodata; |
|
149 |
}; |
|
150 | ||
151 |
template<typename T> |
|
152 |
class Raster : public AbstractRaster<T> { |
|
153 |
public: |
|
154 |
Raster(Matrix<T>&& values, const Box & box) : AbstractRaster<T>( |
|
155 |
Grid<bounded_extent>( |
|
156 |
box, |
|
157 |
(box.xmax - box.xmin) / values.cols(), |
|
158 |
(box.ymax - box.ymin) / values.rows())), |
|
159 |
m_values{std::move(values)} {} |
|
160 | ||
161 | 516717x |
Raster(Matrix<T>&& values, const Grid<bounded_extent> & g) : |
162 |
AbstractRaster<T>(g), |
|
163 | 516717x |
m_values{std::move(values)} {} |
164 | ||
165 |
Raster(const Box & box, size_t nrow, size_t ncol) : |
|
166 |
AbstractRaster<T>(Grid<bounded_extent>(box, (box.xmax-box.xmin) / ncol, (box.ymax-box.ymin) / nrow)), |
|
167 |
m_values{nrow, ncol} |
|
168 |
{} |
|
169 | ||
170 |
explicit Raster(const Grid<bounded_extent> & ex) : |
|
171 |
AbstractRaster<T>(ex), |
|
172 |
m_values{ex.rows(), ex.cols()} |
|
173 |
{} |
|
174 | ||
175 |
Matrix<T>& data() { |
|
176 |
return m_values; |
|
177 |
} |
|
178 | ||
179 |
T& operator()(size_t row, size_t col) { |
|
180 |
return m_values(row, col); |
|
181 |
} |
|
182 | ||
183 | 623683x |
T operator()(size_t row, size_t col) const override { |
184 | 623683x |
return m_values(row, col); |
185 |
} |
|
186 | ||
187 |
private: |
|
188 |
Matrix<T> m_values; |
|
189 |
}; |
|
190 | ||
191 |
template<typename T> |
|
192 |
class RasterView : public AbstractRaster<T>{ |
|
193 |
public: |
|
194 |
// Construct a view of a raster r at an extent ex that is larger |
|
195 |
// and/or of finer resolution than r |
|
196 | 798226x |
RasterView(const AbstractRaster<T> & r, Grid<bounded_extent> ex) : AbstractRaster<T>(ex), m_raster{r} { |
197 | 798226x |
double disaggregation_factor_x = r.xres() / ex.dx(); |
198 | 798226x |
double disaggregation_factor_y = r.yres() / ex.dy(); |
199 | ||
200 | ! |
if (std::abs(disaggregation_factor_x - std::floor(disaggregation_factor_x)) > 1e-6 || |
201 | 798226x |
std::abs(disaggregation_factor_y - std::floor(disaggregation_factor_y)) > 1e-6) { |
202 | ! |
throw std::runtime_error("Must construct view at resolution that is an integer multiple of original."); |
203 |
} |
|
204 | ||
205 | 798226x |
if (disaggregation_factor_x < 0 || disaggregation_factor_y < 0) { |
206 | ! |
throw std::runtime_error("Must construct view at equal or higher resolution than original."); |
207 |
} |
|
208 | ||
209 | 798226x |
m_x_off = static_cast<long>(std::round((ex.xmin() - r.xmin()) / ex.dx())); |
210 | 798226x |
m_y_off = static_cast<long>(std::round((r.ymax() - ex.ymax()) / ex.dy())); |
211 | 798226x |
m_rx = static_cast<size_t>(disaggregation_factor_x); |
212 | 798226x |
m_ry = static_cast<size_t>(disaggregation_factor_y); |
213 | ||
214 | 798226x |
if (r.has_nodata()) { |
215 | ! |
this->set_nodata(r.nodata()); |
216 |
} |
|
217 |
} |
|
218 | ||
219 | 766555x |
T operator()(size_t row, size_t col) const override { |
220 | ! |
if (m_x_off < 0 && static_cast<size_t>(-m_x_off) > col) { |
221 | ! |
return this->nodata(); |
222 |
} |
|
223 | 732837x |
if (m_y_off < 0 && static_cast<size_t>(-m_y_off) > row) { |
224 | 19144x |
return this->nodata(); |
225 |
} |
|
226 | ||
227 | 713693x |
size_t i0 = (row + m_y_off) / m_ry; |
228 | 713693x |
size_t j0 = (col + m_x_off) / m_rx; |
229 | ||
230 | 713693x |
if (i0 > m_raster.rows() - 1 || j0 > m_raster.cols() - 1) { |
231 | 86844x |
return this->nodata(); |
232 |
} |
|
233 | ||
234 | 626849x |
return m_raster(i0, j0); |
235 |
} |
|
236 | ||
237 |
private: |
|
238 |
const AbstractRaster<T>& m_raster; |
|
239 | ||
240 |
long m_x_off; |
|
241 |
long m_y_off; |
|
242 |
size_t m_rx; |
|
243 |
size_t m_ry; |
|
244 |
}; |
|
245 | ||
246 | ||
247 |
template<typename T> |
|
248 |
std::ostream& operator<<(std::ostream & os, const AbstractRaster<T> & m) { |
|
249 |
for (size_t i = 0; i < m.rows(); i++) { |
|
250 |
for (size_t j = 0; j < m.cols(); j++) { |
|
251 |
if (m(i, j) != 0) { |
|
252 |
os << std::right << std::fixed << std::setw(10) << std::setprecision(6) << |
|
253 |
m(i, j) << " "; |
|
254 |
} else { |
|
255 |
os << " "; |
|
256 |
} |
|
257 |
} |
|
258 |
os << std::endl; |
|
259 |
} |
|
260 | ||
261 |
return os; |
|
262 |
} |
|
263 |
} |
|
264 | ||
265 |
#endif //EXACTEXTRACT_RASTER_H |
1 |
// Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
// [[Rcpp::plugins("cpp14")]] |
|
15 | ||
16 |
#pragma once |
|
17 | ||
18 |
#include <memory> |
|
19 | ||
20 |
#include <Rcpp.h> |
|
21 | ||
22 |
#include <geos_c.h> |
|
23 | ||
24 |
using geom_ptr = std::unique_ptr<GEOSGeometry, std::function<void(GEOSGeometry*)>>; |
|
25 |
using wkb_reader_ptr = std::unique_ptr<GEOSWKBReader, std::function<void(GEOSWKBReader*)>>; |
|
26 | ||
27 |
// GEOS warning handler |
|
28 | ! |
static void geos_warn(const char* fmt, ...) { |
29 | ! |
char buf[BUFSIZ] = { '\0' }; |
30 | ||
31 |
va_list msg; |
|
32 | ! |
va_start(msg, fmt); |
33 | ! |
vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg); |
34 | ! |
va_end(msg); |
35 | ||
36 | ! |
Rcpp::Function warning("warning"); |
37 | ! |
warning(buf); |
38 |
} |
|
39 | ||
40 |
// GEOS error handler |
|
41 | ! |
static void geos_error(const char* fmt, ...) { |
42 | ! |
char buf[BUFSIZ] = { '\0' }; |
43 | ||
44 |
va_list msg; |
|
45 | ! |
va_start(msg, fmt); |
46 | ! |
vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg); |
47 | ! |
va_end(msg); |
48 | ||
49 | ! |
Rcpp::stop(buf); |
50 |
} |
|
51 | ||
52 |
// GEOSContextHandle wrapper to ensure finishGEOS is called. |
|
53 |
struct GEOSAutoHandle { |
|
54 | 247x |
GEOSAutoHandle() { |
55 | 247x |
handle = initGEOS_r(geos_warn, geos_error); |
56 |
} |
|
57 | ||
58 | 494x |
~GEOSAutoHandle() { |
59 | 247x |
finishGEOS_r(handle); |
60 |
} |
|
61 | ||
62 |
GEOSContextHandle_t handle; |
|
63 |
}; |
|
64 | ||
65 |
// Return a smart pointer to a Geometry, given WKB input |
|
66 | 246x |
static inline geom_ptr read_wkb(const GEOSContextHandle_t & context, const Rcpp::RawVector & wkb) { |
67 | 738x |
wkb_reader_ptr wkb_reader{ GEOSWKBReader_create_r(context), [context](GEOSWKBReader* r) { GEOSWKBReader_destroy_r(context, r); } }; |
68 | ||
69 |
geom_ptr geom{ GEOSWKBReader_read_r(context, |
|
70 |
wkb_reader.get(), |
|
71 | 492x |
std::addressof(wkb[0]), |
72 | 246x |
wkb.size()), |
73 | 492x |
[context](GEOSGeometry* g) { GEOSGeom_destroy_r(context, g); } }; |
74 | ||
75 | 246x |
if (geom.get() == nullptr) { |
76 | ! |
Rcpp::stop("Failed to parse WKB geometry"); |
77 |
} |
|
78 | ||
79 | 492x |
return geom; |
80 |
} |
1 |
// Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
// [[Rcpp::plugins("cpp14")]] |
|
15 |
#include <memory> |
|
16 | ||
17 |
#include <Rcpp.h> |
|
18 | ||
19 |
#include "geos_r.h" |
|
20 |
#include "raster_utils.h" |
|
21 | ||
22 |
#include "s4_raster_source.h" |
|
23 | ||
24 |
#include "exactextract/src/geos_utils.h" |
|
25 |
#include "exactextract/src/grid.h" |
|
26 |
#include "exactextract/src/matrix.h" |
|
27 |
#include "exactextract/src/raster_cell_intersection.h" |
|
28 |
#include "exactextract/src/raster_source.h" |
|
29 |
#include "exactextract/src/raster_stats.h" |
|
30 | ||
31 |
using exactextract::Box; |
|
32 |
using exactextract::Matrix; |
|
33 |
using exactextract::Grid; |
|
34 |
using exactextract::subdivide; |
|
35 |
using exactextract::bounded_extent; |
|
36 |
using exactextract::Raster; |
|
37 |
using exactextract::RasterView; |
|
38 |
using exactextract::raster_cell_intersection; |
|
39 |
using exactextract::RasterStats; |
|
40 |
using exactextract::RasterSource; |
|
41 | ||
42 |
#include <vector> |
|
43 |
#include <string> |
|
44 | ||
45 |
// [[Rcpp::export]] |
|
46 | 89x |
Rcpp::List CPP_exact_extract(Rcpp::S4 & rast, |
47 |
Rcpp::Nullable<Rcpp::S4> & weights, |
|
48 |
const Rcpp::RawVector & wkb, |
|
49 |
bool include_xy, |
|
50 |
bool include_cell_number, |
|
51 |
Rcpp::Nullable<Rcpp::List> & include_cols, |
|
52 |
Rcpp::CharacterVector & src_names, |
|
53 |
Rcpp::Nullable<Rcpp::CharacterVector> & p_weights_names, |
|
54 |
bool warn_on_disaggregate) { |
|
55 | 178x |
GEOSAutoHandle geos; |
56 | 267x |
Rcpp::Function names("names"); |
57 | ||
58 | 89x |
auto grid = make_grid(rast); |
59 | 89x |
auto weights_grid = exactextract::Grid<bounded_extent>::make_empty(); |
60 | 89x |
auto common_grid = grid; |
61 | ||
62 | 178x |
S4RasterSource rsrc(rast); |
63 | 89x |
int src_nlayers = get_nlayers(rast); |
64 | ||
65 | 89x |
std::unique_ptr<S4RasterSource> rweights; |
66 | 89x |
int weights_nlayers = 0; |
67 | 178x |
Rcpp::CharacterVector weights_names; |
68 | 89x |
if (!weights.isNull()) { |
69 | 34x |
Rcpp::S4 weights_s4 = weights.get(); |
70 | 17x |
weights_nlayers = get_nlayers(weights_s4); |
71 | 17x |
weights_grid = make_grid(weights_s4); |
72 | 17x |
common_grid = grid.common_grid(weights_grid); |
73 | ||
74 | 17x |
rweights = std::make_unique<S4RasterSource>(weights_s4); |
75 | 17x |
weights_names = p_weights_names.get(); |
76 | ||
77 | 17x |
if (warn_on_disaggregate && (common_grid.dx() < grid.dx() || common_grid.dy() < grid.dy())) { |
78 | 3x |
Rcpp::warning("value raster implicitly disaggregated to match higher resolution of weights"); |
79 |
} |
|
80 |
} |
|
81 | ||
82 | 178x |
auto geom = read_wkb(geos.handle, wkb); |
83 | ||
84 | 89x |
auto bbox = exactextract::geos_get_box(geos.handle, geom.get()); |
85 | ||
86 | 89x |
common_grid = common_grid.crop(bbox); |
87 | ||
88 | 178x |
auto coverage_fractions = raster_cell_intersection(common_grid, geos.handle, geom.get()); |
89 | 89x |
auto& cov_grid = coverage_fractions.grid(); |
90 | ||
91 |
// We can't construct an Rcpp::DataFrame with a variable number of columns |
|
92 |
// because of a bug in Rcpp (see https://github.com/RcppCore/Rcpp/pull/1099) |
|
93 |
// Instead, we build up a list, and then create a data frame from the list. |
|
94 |
// Profiling shows that this ends up in as.data.frame, which is slow. |
|
95 |
// Once Rcpp 1.0.6 is released, this can be reworked to be simpler and faster. |
|
96 | 89x |
Rcpp::List cols; |
97 | ||
98 | 178x |
Rcpp::NumericVector coverage_fraction_vec = as_vector(coverage_fractions); |
99 | 178x |
Rcpp::LogicalVector covered = coverage_fraction_vec > 0; |
100 | ||
101 | 89x |
if (include_cols.isNotNull()) { |
102 | 8x |
Rcpp::List include_cols_list = include_cols.get(); |
103 | 8x |
Rcpp::CharacterVector include_names = include_cols_list.attr("names"); |
104 | ||
105 | 8x |
for (int i = 0; i < include_names.size(); i++) { |
106 | 4x |
std::string name(include_names[i]); |
107 | 4x |
cols[name] = include_cols_list[name]; |
108 |
} |
|
109 |
} |
|
110 | ||
111 | 202x |
for (int i = 0; i < src_nlayers; i++) { |
112 | 226x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
113 | 113x |
const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get()); |
114 | ||
115 |
// TODO Perhaps extend this to preserve types (integer, logical.) |
|
116 |
// A bit challenging, since we don't know the type of the raster |
|
117 |
// until we call getValuesBlock, and even then we can't be sure |
|
118 |
// that the returned type is correct (as in a RasterStack with mixed types.) |
|
119 |
// Since R integers are only 32-bit, we are not going to lose data by |
|
120 |
// converting everything to numeric, although we pay a storage penalty. |
|
121 | 113x |
Rcpp::NumericVector value_vec = r->vec(); |
122 | 211x |
if (grid.dx() != common_grid.dx() || grid.dy() != common_grid.dy() || |
123 | 98x |
value_vec.size() != covered.size()) { |
124 |
// Transform values to common grid |
|
125 | 20x |
RasterView<double> rt(*r, common_grid); |
126 | 20x |
value_vec = as_vector(rt); |
127 |
} |
|
128 | ||
129 | 113x |
value_vec = value_vec[covered]; |
130 | 113x |
cols[std::string(src_names[i])] = value_vec; |
131 |
} |
|
132 | ||
133 | 114x |
for (int i = 0; i < weights_nlayers; i++) { |
134 | 50x |
auto values = rweights->read_box(cov_grid.extent(), i); |
135 | 25x |
const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get()); |
136 | ||
137 | 50x |
Rcpp::NumericVector weight_vec = r->vec(); |
138 | ||
139 | 48x |
if (weights_grid.dx() != common_grid.dx() || weights_grid.dy() != common_grid.dy() || |
140 | 23x |
weight_vec.size() != covered.size()) { |
141 |
// Transform weights to common grid |
|
142 | 2x |
RasterView<double> rt (*r, common_grid); |
143 | ||
144 | 2x |
weight_vec = as_vector(rt); |
145 |
} |
|
146 | ||
147 | 25x |
weight_vec = weight_vec[covered]; |
148 | ||
149 | 25x |
std::string colname(weights_names[i]); |
150 | 25x |
if (cols.containsElementNamed(colname.c_str())) { |
151 |
// append ".1" to the column name, to match the behavior of |
|
152 |
// data.frame(). We're safe to just add ".1" (instead of incrementing |
|
153 |
// .1 to .2, for example) because duplicated names within the values |
|
154 |
// or weight stack will already have been made unique by raster::stack() |
|
155 | 2x |
colname = colname + ".1"; // append .1 |
156 |
} |
|
157 | 25x |
cols[colname] = weight_vec; |
158 |
} |
|
159 | ||
160 | 89x |
if (include_xy) { |
161 |
// Include xy values from whichever input raster has the higher resolution |
|
162 | 14x |
if (weights_nlayers > 0 && (weights_grid.dx() < grid.dx() || weights_grid.dy() < grid.dy())) { |
163 | 1x |
Rcpp::S4 weights_s4 = weights.get(); |
164 | 1x |
cols["x"] = get_x_values(weights_s4, cov_grid)[covered]; |
165 | 1x |
cols["y"] = get_y_values(weights_s4, cov_grid)[covered]; |
166 |
} else { |
|
167 | 13x |
cols["x"] = get_x_values(rast, cov_grid)[covered]; |
168 | 13x |
cols["y"] = get_y_values(rast, cov_grid)[covered]; |
169 |
} |
|
170 |
} |
|
171 | ||
172 | 89x |
if (include_cell_number) { |
173 | 6x |
cols["cell"] = get_cell_numbers(rast, cov_grid)[covered]; |
174 |
} |
|
175 | ||
176 | 89x |
cols["coverage_fraction"] = coverage_fraction_vec[covered]; |
177 | ||
178 | 178x |
return cols; |
179 |
} |
|
180 | ||
181 |
// Return a matrix with one row per stat and one row per raster layer |
|
182 |
// [[Rcpp::export]] |
|
183 | 146x |
Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast, |
184 |
Rcpp::Nullable<Rcpp::S4> weights, |
|
185 |
const Rcpp::RawVector & wkb, |
|
186 |
const Rcpp::StringVector & stats, |
|
187 |
int max_cells_in_memory, |
|
188 |
const Rcpp::Nullable<Rcpp::NumericVector> & quantiles) { |
|
189 |
try { |
|
190 | 292x |
GEOSAutoHandle geos; |
191 | ||
192 | 146x |
if (max_cells_in_memory < 1) { |
193 | 1x |
Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory); |
194 |
} |
|
195 | ||
196 | 145x |
int nlayers = get_nlayers(rast); |
197 | ||
198 | 290x |
S4RasterSource rsrc(rast); |
199 | ||
200 | 145x |
std::unique_ptr<S4RasterSource> rweights; |
201 | 145x |
bool weighted = false; |
202 | 145x |
int nweights = 0; |
203 | 145x |
if (weights.isNotNull()) { |
204 | 43x |
Rcpp::S4 weights_s4 = weights.get(); |
205 | 43x |
nweights = get_nlayers(weights_s4); |
206 | ||
207 | 43x |
if (nlayers > 1 && nweights > 1 && nlayers != nweights) { |
208 | ! |
Rcpp::stop("Incompatible number of layers in value and weighting rasters"); |
209 |
} |
|
210 | ||
211 | 43x |
rweights = std::make_unique<S4RasterSource>(weights_s4); |
212 | 43x |
weighted = true; |
213 |
} |
|
214 | ||
215 | 290x |
auto geom = read_wkb(geos.handle, wkb); |
216 | 145x |
auto bbox = exactextract::geos_get_box(geos.handle, geom.get()); |
217 | ||
218 | 145x |
auto grid = weighted ? rsrc.grid().common_grid(rweights->grid()) : rsrc.grid(); |
219 | ||
220 | 144x |
bool disaggregated = (grid.dx() < rsrc.grid().dx() || grid.dy() < rsrc.grid().dy()); |
221 | ||
222 | 144x |
bool store_values = false; |
223 | 144x |
int stat_result_size = 0; |
224 | 305x |
for (const auto & stat : stats) { |
225 | 165x |
store_values = store_values || requires_stored_values(stat); |
226 | ||
227 | 338x |
if (disaggregated && (stat == std::string("count") || |
228 | 173x |
stat == std::string("sum"))) { |
229 | 2x |
Rcpp::stop("Cannot compute 'count' or 'sum' when value raster is disaggregated to resolution of weights."); |
230 |
} |
|
231 | ||
232 | 163x |
if (stat == std::string("quantile")) { |
233 | 8x |
int num_quantiles = 0; |
234 | ||
235 | 8x |
if (quantiles.isNotNull()) { |
236 | 7x |
Rcpp::NumericVector qvec = quantiles.get(); |
237 | 7x |
num_quantiles = qvec.size(); |
238 |
} |
|
239 | ||
240 | 8x |
if (num_quantiles == 0) { |
241 | 2x |
Rcpp::stop("Quantiles not specified."); |
242 |
} |
|
243 | ||
244 | 6x |
stat_result_size += num_quantiles; |
245 |
} else { |
|
246 | 155x |
stat_result_size += 1; |
247 |
} |
|
248 |
} |
|
249 | ||
250 | 140x |
int nresults = std::max(nlayers, nweights); |
251 | ||
252 | 280x |
std::vector<RasterStats<double>> raster_stats; |
253 | 140x |
raster_stats.reserve(nresults); |
254 | 291x |
for (int i = 0; i < nresults; i++) { |
255 | 151x |
raster_stats.emplace_back(store_values); |
256 |
} |
|
257 | ||
258 | 280x |
Rcpp::NumericMatrix stat_results = Rcpp::no_init(nresults, stat_result_size); |
259 | ||
260 | 140x |
if (bbox.intersects(grid.extent())) { |
261 | 131x |
auto cropped_grid = grid.crop(bbox); |
262 | ||
263 | 297x |
for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) { |
264 | 332x |
auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get()); |
265 | 166x |
auto& cov_grid = coverage_fraction.grid(); |
266 | ||
267 | 166x |
if (!cov_grid.empty()) { |
268 | 166x |
if (weighted) { |
269 | 40x |
if (nlayers > nweights) { |
270 |
// recycle weights |
|
271 | 4x |
auto weights = rweights->read_box(cov_grid.extent(), 0); |
272 | ||
273 | 7x |
for (int i = 0; i < nlayers; i++) { |
274 | 10x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
275 | 5x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
276 |
} |
|
277 | 38x |
} else if (nweights > nlayers) { |
278 |
// recycle values |
|
279 | 2x |
auto values = rsrc.read_box(cov_grid.extent(), 0); |
280 | ||
281 | 4x |
for (int i = 0; i < nweights; i++) { |
282 | 6x |
auto weights = rweights->read_box(cov_grid.extent(), i); |
283 | 3x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
284 |
} |
|
285 |
} else { |
|
286 |
// process values and weights in parallel |
|
287 | 76x |
for (int i = 0; i < nlayers; i++) { |
288 | 78x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
289 | 78x |
auto weights = rweights->read_box(cov_grid.extent(), i); |
290 | ||
291 | 39x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
292 |
} |
|
293 |
} |
|
294 |
} else { |
|
295 | 256x |
for (int i = 0; i < nlayers; i++) { |
296 | 260x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
297 | 130x |
raster_stats[i].process(coverage_fraction, *values); |
298 |
} |
|
299 |
} |
|
300 |
} |
|
301 |
} |
|
302 |
} |
|
303 | ||
304 | 288x |
for (int j = 0; j < nresults; j++) { |
305 | 151x |
const auto& rs = raster_stats[j]; |
306 | ||
307 | 151x |
int i = 0; |
308 | 326x |
for(const auto& stat : stats) { |
309 | 178x |
if (stat == std::string("mean")) stat_results(j, i++) = rs.mean(); |
310 | ||
311 | 150x |
else if (stat == std::string("sum")) stat_results(j, i++) = rs.sum(); |
312 | 93x |
else if (stat == std::string("count")) stat_results(j, i++) = rs.count(); |
313 | ||
314 | 86x |
else if (stat == std::string("min")) stat_results(j, i++) = rs.min().value_or(NA_REAL); |
315 | 80x |
else if (stat == std::string("max")) stat_results(j, i++) = rs.max().value_or(NA_REAL); |
316 | ||
317 | 74x |
else if (stat == std::string("median")) stat_results(j, i++) = rs.quantile(0.5).value_or(NA_REAL); |
318 | ||
319 | 73x |
else if (stat == std::string("mode")) stat_results(j, i++) = rs.mode().value_or(NA_REAL); |
320 | 68x |
else if (stat == std::string("majority")) stat_results(j, i++) = rs.mode().value_or(NA_REAL); |
321 | 65x |
else if (stat == std::string("minority")) stat_results(j, i++) = rs.minority().value_or(NA_REAL); |
322 | ||
323 | 61x |
else if (stat == std::string("variety")) stat_results(j, i++) = rs.variety(); |
324 | 49x |
else if (stat == std::string("weighted_mean")) stat_results(j, i++) = rs.weighted_mean(); |
325 | 15x |
else if (stat == std::string("weighted_sum")) stat_results(j, i++) = rs.weighted_sum(); |
326 | ||
327 | 10x |
else if (stat == std::string("variance")) stat_results(j, i++) = rs.variance(); |
328 | 9x |
else if (stat == std::string("stdev")) stat_results(j, i++) = rs.stdev(); |
329 | 8x |
else if (stat == std::string("coefficient_of_variation")) stat_results(j, i++) = rs.coefficient_of_variation(); |
330 | ||
331 | 7x |
else if (stat == std::string("quantile")) { |
332 | 12x |
Rcpp::NumericVector qvec = quantiles.get(); |
333 | 15x |
for (double q : qvec) { |
334 | 11x |
stat_results(j, i++) = rs.quantile(q).value_or(NA_REAL); |
335 |
} |
|
336 |
} |
|
337 | ||
338 | 4x |
else Rcpp::stop("Unknown stat: " + stat); |
339 |
} |
|
340 |
} |
|
341 | ||
342 | 274x |
return stat_results; |
343 | 18x |
} catch (std::exception & e) { |
344 |
// throw predictable exception class |
|
345 | 9x |
Rcpp::stop(e.what()); |
346 |
} |
|
347 |
} |
1 |
// Copyright (c) 2018-2019 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_BOX_H |
|
15 |
#define EXACTEXTRACT_BOX_H |
|
16 | ||
17 |
#include "coordinate.h" |
|
18 |
#include "crossing.h" |
|
19 |
#include "side.h" |
|
20 | ||
21 |
#include <limits> |
|
22 | ||
23 |
namespace exactextract { |
|
24 | ||
25 |
struct Box { |
|
26 |
double xmin; |
|
27 |
double ymin; |
|
28 |
double xmax; |
|
29 |
double ymax; |
|
30 | ||
31 | 8721020x |
Box(double p_xmin, double p_ymin, double p_xmax, double p_ymax) : |
32 |
xmin{p_xmin}, |
|
33 |
ymin{p_ymin}, |
|
34 |
xmax{p_xmax}, |
|
35 | 8721020x |
ymax{p_ymax} {} |
36 | ||
37 |
static Box maximum_finite() { |
|
38 |
return { |
|
39 |
std::numeric_limits<double>::lowest(), |
|
40 |
std::numeric_limits<double>::lowest(), |
|
41 |
std::numeric_limits<double>::max(), |
|
42 |
std::numeric_limits<double>::max() |
|
43 |
}; |
|
44 |
} |
|
45 | ||
46 | 267x |
static Box make_empty() { |
47 | 267x |
return {0, 0, 0, 0}; |
48 |
} |
|
49 | ||
50 | 6372234x |
double width() const { |
51 | 6372234x |
return xmax - xmin; |
52 |
} |
|
53 | ||
54 | 6372234x |
double height() const { |
55 | 6372234x |
return ymax - ymin; |
56 |
} |
|
57 | ||
58 | 6366768x |
double area() const { |
59 | 6366768x |
return width() * height(); |
60 |
} |
|
61 | ||
62 | 2733x |
double perimeter() const { |
63 | 2733x |
return 2 * width() + 2 * height(); |
64 |
} |
|
65 | ||
66 | 400715x |
bool intersects(const Box & other) const { |
67 | 400715x |
if (other.ymin > ymax) |
68 | 42x |
return false; |
69 | 400673x |
if (other.ymax < ymin) |
70 | 50x |
return false; |
71 | 400623x |
if (other.xmin > xmax) |
72 | ! |
return false; |
73 | 400623x |
if (other.xmax < xmin) |
74 | ! |
return false; |
75 | ||
76 | 400623x |
return true; |
77 |
} |
|
78 | ||
79 | 4098885x |
Box intersection(const Box & other) const { |
80 |
return { |
|
81 | 4098885x |
std::max(xmin, other.xmin), |
82 | 4098885x |
std::max(ymin, other.ymin), |
83 | 4098885x |
std::min(xmax, other.xmax), |
84 | 4098885x |
std::min(ymax, other.ymax) |
85 |
}; |
|
86 |
} |
|
87 | ||
88 |
Box translate(double dx, double dy) const { |
|
89 |
return { |
|
90 |
xmin + dx, |
|
91 |
ymin + dy, |
|
92 |
xmax + dx, |
|
93 |
ymax + dy |
|
94 |
}; |
|
95 |
} |
|
96 | ||
97 |
Coordinate upper_left() const { |
|
98 |
return Coordinate{xmin, ymax}; |
|
99 |
} |
|
100 | ||
101 |
Coordinate upper_right() const { |
|
102 |
return Coordinate{xmax, ymax}; |
|
103 |
} |
|
104 | ||
105 |
Coordinate lower_left() const { |
|
106 |
return Coordinate{xmin, ymin}; |
|
107 |
} |
|
108 | ||
109 |
Coordinate lower_right() const { |
|
110 |
return Coordinate{xmax, ymin}; |
|
111 |
} |
|
112 | ||
113 |
Side side(const Coordinate &c) const; |
|
114 | ||
115 |
Crossing crossing(const Coordinate &c1, const Coordinate &c2) const; |
|
116 | ||
117 | 516987x |
bool empty() const { |
118 | 516987x |
return xmin >= xmax || ymin >= ymax; |
119 |
} |
|
120 | ||
121 | 3x |
Box expand_to_include(const Box & other) const { |
122 | 3x |
if (empty()) { |
123 | ! |
return other; |
124 |
} |
|
125 | ||
126 | 3x |
if (other.empty()) { |
127 | ! |
return *this; |
128 |
} |
|
129 | ||
130 | 3x |
return { std::min(xmin, other.xmin), |
131 | 3x |
std::min(ymin, other.ymin), |
132 | 3x |
std::max(xmax, other.xmax), |
133 | 3x |
std::max(ymax, other.ymax) }; |
134 |
} |
|
135 | ||
136 |
bool contains(const Box &b) const; |
|
137 | ||
138 |
bool contains(const Coordinate &c) const; |
|
139 | ||
140 |
bool strictly_contains(const Coordinate &c) const; |
|
141 | ||
142 | 1242x |
bool operator==(const Box& other) const { |
143 | 1242x |
return xmin == other.xmin && xmax == other.xmax && ymin == other.ymin && ymax == other.ymax; |
144 |
} |
|
145 | ||
146 |
friend std::ostream &operator<<(std::ostream &os, const Box &c); |
|
147 |
}; |
|
148 | ||
149 |
std::ostream &operator<<(std::ostream &os, const Box &c); |
|
150 | ||
151 |
} |
|
152 | ||
153 |
#endif |
1 |
// Copyright (c) 2018-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_RASTER_STATS_H |
|
15 |
#define EXACTEXTRACT_RASTER_STATS_H |
|
16 | ||
17 |
#include <algorithm> |
|
18 |
#include <iostream> |
|
19 |
#include <limits> |
|
20 |
#include <unordered_map> |
|
21 | ||
22 |
#include "raster_cell_intersection.h" |
|
23 |
#include "weighted_quantiles.h" |
|
24 |
#include "variance.h" |
|
25 | ||
26 |
#include "../vend/optional.hpp" |
|
27 | ||
28 |
namespace exactextract { |
|
29 | ||
30 |
template<typename T> |
|
31 |
class RasterStats { |
|
32 | ||
33 |
public: |
|
34 |
/** |
|
35 |
* Compute raster statistics from a Raster representing intersection percentages, |
|
36 |
* a Raster representing data values, and (optionally) a Raster representing weights. |
|
37 |
* and a set of raster values. |
|
38 |
*/ |
|
39 | 516601x |
explicit RasterStats(bool store_values = false) : |
40 | 516601x |
m_min{std::numeric_limits<T>::max()}, |
41 | 516601x |
m_max{std::numeric_limits<T>::lowest()}, |
42 |
m_sum_ciwi{0}, |
|
43 |
m_sum_ci{0}, |
|
44 |
m_sum_xici{0}, |
|
45 |
m_sum_xiciwi{0}, |
|
46 | 516601x |
m_store_values{store_values} {} |
47 | ||
48 | 398938x |
void process(const Raster<float> & intersection_percentages, const AbstractRaster<T> & rast) { |
49 | 797876x |
RasterView<T> rv{rast, intersection_percentages.grid()}; |
50 | ||
51 | 904507x |
for (size_t i = 0; i < rv.rows(); i++) { |
52 | 1106375x |
for (size_t j = 0; j < rv.cols(); j++) { |
53 | 600806x |
float pct_cov = intersection_percentages(i, j); |
54 |
T val; |
|
55 | 600806x |
if (pct_cov > 0 && rv.get(i, j, val)) { |
56 | 600197x |
process_value(val, pct_cov, 1.0); |
57 |
} |
|
58 |
} |
|
59 |
} |
|
60 |
} |
|
61 | ||
62 | 47x |
void process(const Raster<float> & intersection_percentages, const AbstractRaster<T> & rast, const AbstractRaster<T> & weights) { |
63 |
// Process the entire intersection_percentages grid, even though it may |
|
64 |
// be outside the extent of the weighting raster. Although we've been |
|
65 |
// provided a weighting raster, we still need to calculate correct values |
|
66 |
// for unweighted stats. |
|
67 | 47x |
auto& common = intersection_percentages.grid(); |
68 | ||
69 | 47x |
if (common.empty()) |
70 | ! |
return; |
71 | ||
72 | 94x |
RasterView<float> iv{intersection_percentages, common}; |
73 | 94x |
RasterView<T> rv{rast, common}; |
74 | 94x |
RasterView<T> wv{weights, common}; |
75 | ||
76 | 307x |
for (size_t i = 0; i < rv.rows(); i++) { |
77 | 3013x |
for (size_t j = 0; j < rv.cols(); j++) { |
78 | 2753x |
float pct_cov = iv(i, j); |
79 |
T weight; |
|
80 |
T val; |
|
81 | ||
82 | 2753x |
if (pct_cov > 0 && rv.get(i, j, val)) { |
83 | 2420x |
if (wv.get(i, j, weight)) { |
84 | 2417x |
process_value(val, pct_cov, weight); |
85 |
} else { |
|
86 |
// Weight is NODATA, convert to NAN |
|
87 | 3x |
process_value(val, pct_cov, std::numeric_limits<double>::quiet_NaN()); |
88 |
} |
|
89 |
} |
|
90 |
} |
|
91 |
} |
|
92 |
} |
|
93 | ||
94 |
/** |
|
95 |
* The mean value of cells covered by this polygon, weighted |
|
96 |
* by the percent of the cell that is covered. |
|
97 |
*/ |
|
98 | 28x |
float mean() const { |
99 | 28x |
return sum() / count(); |
100 |
} |
|
101 | ||
102 |
/** |
|
103 |
* The mean value of cells covered by this polygon, weighted |
|
104 |
* by the percent of the cell that is covered and a secondary |
|
105 |
* weighting raster. |
|
106 |
* |
|
107 |
* If any weights are undefined, will return NAN. If this is undesirable, |
|
108 |
* caller should replace undefined weights with a suitable default |
|
109 |
* before computing statistics. |
|
110 |
*/ |
|
111 | 34x |
float weighted_mean() const { |
112 | 34x |
return weighted_sum() / weighted_count(); |
113 |
} |
|
114 | ||
115 |
/** The fraction of weighted cells to unweighted cells. |
|
116 |
* Meaningful only when the values of the weighting |
|
117 |
* raster are between 0 and 1. |
|
118 |
*/ |
|
119 |
float weighted_fraction() const { |
|
120 |
return weighted_sum() / sum(); |
|
121 |
} |
|
122 | ||
123 |
/** |
|
124 |
* The raster value occupying the greatest number of cells |
|
125 |
* or partial cells within the polygon. When multiple values |
|
126 |
* cover the same number of cells, the greatest value will |
|
127 |
* be returned. Weights are not taken into account. |
|
128 |
*/ |
|
129 | 16208x |
nonstd::optional<T> mode() const { |
130 | 16208x |
if (variety() == 0) { |
131 | 91x |
return nonstd::nullopt; |
132 |
} |
|
133 | ||
134 | 32234x |
return std::max_element(m_freq.cbegin(), |
135 |
m_freq.cend(), |
|
136 | 27262x |
[](const auto &a, const auto &b) { |
137 | 27262x |
return a.second < b.second || (a.second == b.second && a.first < b.first); |
138 | 16117x |
})->first; |
139 |
} |
|
140 | ||
141 |
/** |
|
142 |
* The minimum value in any raster cell wholly or partially covered |
|
143 |
* by the polygon. Weights are not taken into account. |
|
144 |
*/ |
|
145 | 6x |
nonstd::optional<T> min() const { |
146 | 6x |
if (m_sum_ci == 0) { |
147 | 1x |
return nonstd::nullopt; |
148 |
} |
|
149 | 5x |
return m_min; |
150 |
} |
|
151 | ||
152 |
/** |
|
153 |
* The maximum value in any raster cell wholly or partially covered |
|
154 |
* by the polygon. Weights are not taken into account. |
|
155 |
*/ |
|
156 | 6x |
nonstd::optional<T> max() const { |
157 | 6x |
if (m_sum_ci == 0) { |
158 | 1x |
return nonstd::nullopt; |
159 |
} |
|
160 | 5x |
return m_max; |
161 |
} |
|
162 | ||
163 |
/** |
|
164 |
* The given quantile (0-1) of raster cell values. Coverage fractions |
|
165 |
* are taken into account but weights are not. |
|
166 |
*/ |
|
167 | 12x |
nonstd::optional<T> quantile(double q) const { |
168 | 12x |
if (m_sum_ci == 0) { |
169 | ! |
return nonstd::nullopt; |
170 |
} |
|
171 | ||
172 |
// The weighted quantile computation is not processed incrementally. |
|
173 |
// Create it on demand and retain it in case we want multiple quantiles. |
|
174 | 12x |
if (!m_quantiles) { |
175 | 7x |
m_quantiles = std::make_unique<WeightedQuantiles>(); |
176 | ||
177 | 62x |
for (const auto& entry : m_freq) { |
178 | 55x |
m_quantiles->process(entry.first, entry.second); |
179 |
} |
|
180 |
} |
|
181 | ||
182 | 22x |
return m_quantiles->quantile(q); |
183 |
} |
|
184 | ||
185 |
/** |
|
186 |
* The sum of raster cells covered by the polygon, with each raster |
|
187 |
* value weighted by its coverage fraction. |
|
188 |
*/ |
|
189 | 500335x |
float sum() const { |
190 | 500335x |
return (float) m_sum_xici; |
191 |
} |
|
192 | ||
193 |
/** |
|
194 |
* The sum of raster cells covered by the polygon, with each raster |
|
195 |
* value weighted by its coverage fraction and weighting raster value. |
|
196 |
* |
|
197 |
* If any weights are undefined, will return NAN. If this is undesirable, |
|
198 |
* caller should replace undefined weights with a suitable default |
|
199 |
* before computing statistics. |
|
200 |
*/ |
|
201 | 39x |
float weighted_sum() const { |
202 | 39x |
return (float) m_sum_xiciwi; |
203 |
} |
|
204 | ||
205 |
/** |
|
206 |
* The number of raster cells with a defined value |
|
207 |
* covered by the polygon. Weights are not taken |
|
208 |
* into account. |
|
209 |
*/ |
|
210 | 35x |
float count() const { |
211 | 35x |
return (float) m_sum_ci; |
212 |
} |
|
213 | ||
214 |
/** |
|
215 |
* The population variance of raster cells touched |
|
216 |
* by the polygon. Cell coverage fractions are taken |
|
217 |
* into account; values of a weighting raster are not. |
|
218 |
*/ |
|
219 | 1x |
float variance() const { |
220 | 1x |
return static_cast<float>(m_variance.variance()); |
221 |
} |
|
222 | ||
223 |
/** |
|
224 |
* The population standard deviation of raster cells |
|
225 |
* touched by the polygon. Cell coverage fractions |
|
226 |
* are taken into account; values of a weighting |
|
227 |
* raster are not. |
|
228 |
*/ |
|
229 | 1x |
float stdev() const { |
230 | 1x |
return static_cast<float>(m_variance.stdev()); |
231 |
} |
|
232 | ||
233 |
/** |
|
234 |
* The population coefficient of variation of raster |
|
235 |
* cells touched by the polygon. Cell coverage fractions |
|
236 |
* are taken into account; values of a weighting |
|
237 |
* raster are not. |
|
238 |
*/ |
|
239 | 1x |
float coefficient_of_variation() const { |
240 | 1x |
return static_cast<float>(m_variance.coefficent_of_variation()); |
241 |
} |
|
242 | ||
243 |
/** |
|
244 |
* The sum of weights for each cell covered by the |
|
245 |
* polygon, with each weight multiplied by the coverage |
|
246 |
* coverage fraction of each cell. |
|
247 |
* |
|
248 |
* If any weights are undefined, will return NAN. If this is undesirable, |
|
249 |
* caller should replace undefined weights with a suitable default |
|
250 |
* before computing statistics. |
|
251 |
*/ |
|
252 | 34x |
float weighted_count() const { |
253 | 34x |
return (float) m_sum_ciwi; |
254 |
} |
|
255 | ||
256 |
/** |
|
257 |
* The raster value occupying the least number of cells |
|
258 |
* or partial cells within the polygon. When multiple values |
|
259 |
* cover the same number of cells, the lowest value will |
|
260 |
* be returned. |
|
261 |
* |
|
262 |
* Cell weights are not taken into account. |
|
263 |
*/ |
|
264 | 4x |
nonstd::optional<T> minority() const { |
265 | 4x |
if (variety() == 0) { |
266 | 2x |
return nonstd::nullopt; |
267 |
} |
|
268 | ||
269 | 4x |
return std::min_element(m_freq.cbegin(), |
270 |
m_freq.cend(), |
|
271 | 11x |
[](const auto &a, const auto &b) { |
272 | 11x |
return a.second < b.second || (a.second == b.second && a.first < b.first); |
273 | 2x |
})->first; |
274 |
} |
|
275 | ||
276 |
/** |
|
277 |
* The number of distinct defined raster values in cells wholly |
|
278 |
* or partially covered by the polygon. |
|
279 |
*/ |
|
280 | 16224x |
size_t variety() const { |
281 | 16224x |
return m_freq.size(); |
282 |
} |
|
283 | ||
284 |
bool stores_values() const { |
|
285 |
return m_store_values; |
|
286 |
} |
|
287 | ||
288 |
private: |
|
289 |
T m_min; |
|
290 |
T m_max; |
|
291 | ||
292 |
// ci: coverage fraction of pixel i |
|
293 |
// wi: weight of pixel i |
|
294 |
// xi: value of pixel i |
|
295 |
double m_sum_ciwi; |
|
296 |
double m_sum_ci; |
|
297 |
double m_sum_xici; |
|
298 |
double m_sum_xiciwi; |
|
299 |
WestVariance m_variance; |
|
300 | ||
301 |
mutable std::unique_ptr<WeightedQuantiles> m_quantiles; |
|
302 | ||
303 |
std::unordered_map<T, float> m_freq; |
|
304 | ||
305 |
bool m_store_values; |
|
306 | ||
307 | 602617x |
void process_value(const T& val, float coverage, double weight) { |
308 | 602617x |
m_sum_ci += static_cast<double>(coverage); |
309 | 602617x |
m_sum_xici += val*static_cast<double>(coverage); |
310 | ||
311 | 602617x |
m_variance.process(val, coverage); |
312 | ||
313 | 602617x |
double ciwi = static_cast<double>(coverage)*weight; |
314 | 602617x |
m_sum_ciwi += ciwi; |
315 | 602617x |
m_sum_xiciwi += val * ciwi; |
316 | ||
317 | 602617x |
if (val < m_min) { |
318 | 477392x |
m_min = val; |
319 |
} |
|
320 | ||
321 | 602617x |
if (val > m_max) { |
322 | 478491x |
m_max = val; |
323 |
} |
|
324 | ||
325 |
// TODO should weights factor in here? |
|
326 | 602617x |
if (m_store_values) { |
327 | 64127x |
m_freq[val] += coverage; |
328 | 64127x |
m_quantiles.reset(); |
329 |
} |
|
330 |
} |
|
331 |
}; |
|
332 | ||
333 |
template<typename T> |
|
334 |
std::ostream& operator<<(std::ostream& os, const RasterStats<T> & stats) { |
|
335 |
os << "{" << std::endl; |
|
336 |
os << " \"count\" : " << stats.count() << "," << std::endl; |
|
337 | ||
338 |
os << " \"min\" : "; |
|
339 |
if (stats.min().has_value()) { |
|
340 |
os << stats.min().value(); |
|
341 |
} else { |
|
342 |
os << "null"; |
|
343 |
} |
|
344 |
os << "," << std::endl; |
|
345 | ||
346 |
os << " \"max\" : "; |
|
347 |
if (stats.max().has_value()) { |
|
348 |
os << stats.max().value(); |
|
349 |
} else { |
|
350 |
os << "null"; |
|
351 |
} |
|
352 |
os << "," << std::endl; |
|
353 | ||
354 |
os << " \"mean\" : " << stats.mean() << "," << std::endl; |
|
355 |
os << " \"sum\" : " << stats.sum() << "," << std::endl; |
|
356 |
os << " \"weighted_mean\" : " << stats.weighted_mean() << "," << std::endl; |
|
357 |
os << " \"weighted_sum\" : " << stats.weighted_sum(); |
|
358 |
if (stats.stores_values()) { |
|
359 |
os << "," << std::endl; |
|
360 |
os << " \"mode\" : "; |
|
361 |
if (stats.mode().has_value()) { |
|
362 |
os << stats.mode().value(); |
|
363 |
} else { |
|
364 |
os << "null"; |
|
365 |
} |
|
366 |
os << "," << std::endl; |
|
367 | ||
368 |
os << " \"minority\" : "; |
|
369 |
if (stats.minority().has_value()) { |
|
370 |
os << stats.minority().value(); |
|
371 |
} else { |
|
372 |
os << "null"; |
|
373 |
} |
|
374 |
os << "," << std::endl; |
|
375 | ||
376 |
os << " \"variety\" : " << stats.variety() << std::endl; |
|
377 |
} else { |
|
378 |
os << std::endl; |
|
379 |
} |
|
380 |
os << "}" << std::endl; |
|
381 |
return os; |
|
382 |
} |
|
383 | ||
384 | ||
385 |
} |
|
386 | ||
387 |
#endif |
1 |
// Copyright (c) 2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_VARIANCE_H |
|
15 |
#define EXACTEXTRACT_VARIANCE_H |
|
16 | ||
17 |
#include <cmath> |
|
18 | ||
19 |
namespace exactextract { |
|
20 | ||
21 |
class WestVariance { |
|
22 |
/** \brief Implements an incremental algorithm for weighted standard |
|
23 |
* deviation, variance, and coefficient of variation, as described in |
|
24 |
* formula WV2 of West, D.H.D. (1979) "Updating Mean and Variance |
|
25 |
* Estimates: An Improved Method". Communications of the ACM 22(9). |
|
26 |
*/ |
|
27 | ||
28 |
private: |
|
29 |
double sum_w = 0; |
|
30 |
double mean = 0; |
|
31 |
double t = 0; |
|
32 | ||
33 |
public: |
|
34 |
/** \brief Update variance estimate with another value |
|
35 |
* |
|
36 |
* @param x value to add |
|
37 |
* @param w weight of `x` |
|
38 |
*/ |
|
39 | 602617x |
void process(double x, double w) { |
40 | 602617x |
double mean_old = mean; |
41 | ||
42 | 602617x |
sum_w += w; |
43 | 602617x |
mean += (w / sum_w) * (x - mean_old); |
44 | 602617x |
t += w * (x - mean_old) * (x - mean); |
45 |
} |
|
46 | ||
47 |
/** \brief Return the population variance. |
|
48 |
*/ |
|
49 | 3x |
constexpr double variance() const { |
50 | 3x |
return t / sum_w; |
51 |
} |
|
52 | ||
53 |
/** \brief Return the population standard deviation |
|
54 |
*/ |
|
55 | 2x |
double stdev() const { |
56 | 2x |
return std::sqrt(variance()); |
57 |
} |
|
58 | ||
59 |
/** \brief Return the population coefficient of variation |
|
60 |
*/ |
|
61 | 1x |
double coefficent_of_variation() const { |
62 | 1x |
return stdev() / mean; |
63 |
} |
|
64 | ||
65 |
}; |
|
66 | ||
67 |
} |
|
68 | ||
69 |
#endif //EXACTEXTRACT_VARIANCE_H |
1 |
// Copyright (c) 2019-2020 ISciences, LLC. |
|
2 |
// All rights reserved. |
|
3 |
// |
|
4 |
// This software is licensed under the Apache License, Version 2.0 (the "License"). |
|
5 |
// You may not use this file except in compliance with the License. You may |
|
6 |
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. |
|
7 |
// |
|
8 |
// Unless required by applicable law or agreed to in writing, software |
|
9 |
// distributed under the License is distributed on an "AS IS" BASIS, |
|
10 |
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
11 |
// See the License for the specific language governing permissions and |
|
12 |
// limitations under the License. |
|
13 | ||
14 |
#ifndef EXACTEXTRACT_WEIGHTED_QUANTILES_H |
|
15 |
#define EXACTEXTRACT_WEIGHTED_QUANTILES_H |
|
16 | ||
17 |
#include <algorithm> |
|
18 |
#include <cmath> |
|
19 |
#include <stdexcept> |
|
20 |
#include <vector> |
|
21 | ||
22 |
namespace exactextract { |
|
23 | ||
24 |
class WeightedQuantiles { |
|
25 |
// Compute weighted quantiles based on https://stats.stackexchange.com/a/13223 |
|
26 | ||
27 |
public: |
|
28 | ||
29 | 55x |
void process(double x, double w) { |
30 | 55x |
if (w < 0) { |
31 | ! |
throw std::runtime_error("Weighted quantile calculation does not support negative weights."); |
32 |
} |
|
33 | ||
34 | 55x |
if (!std::isfinite(w)) { |
35 | ! |
throw std::runtime_error("Weighted quantile does not support non-finite weights."); |
36 |
} |
|
37 | ||
38 | 55x |
m_ready_to_query = false; |
39 | ||
40 | 55x |
m_elems.emplace_back(x, w); |
41 |
} |
|
42 | ||
43 |
double quantile(double q) const; |
|
44 | ||
45 |
private: |
|
46 |
struct elem_t { |
|
47 | 65x |
elem_t(double _x, double _w) : x(_x), w(_w), cumsum(0), s(0) {} |
48 | ||
49 |
double x; |
|
50 |
double w; |
|
51 |
double cumsum; |
|
52 |
double s; |
|
53 |
}; |
|
54 | ||
55 |
void prepare() const; |
|
56 | ||
57 |
mutable std::vector<elem_t> m_elems; |
|
58 |
mutable double m_sum_w; |
|
59 |
mutable bool m_ready_to_query; |
|
60 |
}; |
|
61 | ||
62 |
} |
|
63 | ||
64 |
#endif //EXACTEXTRACT_WEIGHTED_QUANTILES_H |
1 |
// |
|
2 |
// Copyright (c) 2014-2018 Martin Moene |
|
3 |
// |
|
4 |
// https://github.com/martinmoene/optional-lite |
|
5 |
// |
|
6 |
// Distributed under the Boost Software License, Version 1.0. |
|
7 |
// (See accompanying file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
|
8 | ||
9 |
#pragma once |
|
10 | ||
11 |
#ifndef NONSTD_OPTIONAL_LITE_HPP |
|
12 |
#define NONSTD_OPTIONAL_LITE_HPP |
|
13 | ||
14 |
#define optional_lite_MAJOR 3 |
|
15 |
#define optional_lite_MINOR 2 |
|
16 |
#define optional_lite_PATCH 0 |
|
17 | ||
18 |
#define optional_lite_VERSION optional_STRINGIFY(optional_lite_MAJOR) "." optional_STRINGIFY(optional_lite_MINOR) "." optional_STRINGIFY(optional_lite_PATCH) |
|
19 | ||
20 |
#define optional_STRINGIFY( x ) optional_STRINGIFY_( x ) |
|
21 |
#define optional_STRINGIFY_( x ) #x |
|
22 | ||
23 |
// optional-lite configuration: |
|
24 | ||
25 |
#define optional_OPTIONAL_DEFAULT 0 |
|
26 |
#define optional_OPTIONAL_NONSTD 1 |
|
27 |
#define optional_OPTIONAL_STD 2 |
|
28 | ||
29 |
#if !defined( optional_CONFIG_SELECT_OPTIONAL ) |
|
30 |
# define optional_CONFIG_SELECT_OPTIONAL ( optional_HAVE_STD_OPTIONAL ? optional_OPTIONAL_STD : optional_OPTIONAL_NONSTD ) |
|
31 |
#endif |
|
32 | ||
33 |
// Control presence of exception handling (try and auto discover): |
|
34 | ||
35 |
#ifndef optional_CONFIG_NO_EXCEPTIONS |
|
36 |
# if defined(__cpp_exceptions) || defined(__EXCEPTIONS) || defined(_CPPUNWIND) |
|
37 |
# define optional_CONFIG_NO_EXCEPTIONS 0 |
|
38 |
# else |
|
39 |
# define optional_CONFIG_NO_EXCEPTIONS 1 |
|
40 |
# endif |
|
41 |
#endif |
|
42 | ||
43 |
// C++ language version detection (C++20 is speculative): |
|
44 |
// Note: VC14.0/1900 (VS2015) lacks too much from C++14. |
|
45 | ||
46 |
#ifndef optional_CPLUSPLUS |
|
47 |
# if defined(_MSVC_LANG ) && !defined(__clang__) |
|
48 |
# define optional_CPLUSPLUS (_MSC_VER == 1900 ? 201103L : _MSVC_LANG ) |
|
49 |
# else |
|
50 |
# define optional_CPLUSPLUS __cplusplus |
|
51 |
# endif |
|
52 |
#endif |
|
53 | ||
54 |
#define optional_CPP98_OR_GREATER ( optional_CPLUSPLUS >= 199711L ) |
|
55 |
#define optional_CPP11_OR_GREATER ( optional_CPLUSPLUS >= 201103L ) |
|
56 |
#define optional_CPP11_OR_GREATER_ ( optional_CPLUSPLUS >= 201103L ) |
|
57 |
#define optional_CPP14_OR_GREATER ( optional_CPLUSPLUS >= 201402L ) |
|
58 |
#define optional_CPP17_OR_GREATER ( optional_CPLUSPLUS >= 201703L ) |
|
59 |
#define optional_CPP20_OR_GREATER ( optional_CPLUSPLUS >= 202000L ) |
|
60 | ||
61 |
// C++ language version (represent 98 as 3): |
|
62 | ||
63 |
#define optional_CPLUSPLUS_V ( optional_CPLUSPLUS / 100 - (optional_CPLUSPLUS > 200000 ? 2000 : 1994) ) |
|
64 | ||
65 |
// Use C++17 std::optional if available and requested: |
|
66 | ||
67 |
#if optional_CPP17_OR_GREATER && defined(__has_include ) |
|
68 |
# if __has_include( <optional> ) |
|
69 |
# define optional_HAVE_STD_OPTIONAL 1 |
|
70 |
# else |
|
71 |
# define optional_HAVE_STD_OPTIONAL 0 |
|
72 |
# endif |
|
73 |
#else |
|
74 |
# define optional_HAVE_STD_OPTIONAL 0 |
|
75 |
#endif |
|
76 | ||
77 |
#define optional_USES_STD_OPTIONAL ( (optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_STD) || ((optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_DEFAULT) && optional_HAVE_STD_OPTIONAL) ) |
|
78 | ||
79 |
// |
|
80 |
// in_place: code duplicated in any-lite, expected-lite, optional-lite, value-ptr-lite, variant-lite: |
|
81 |
// |
|
82 | ||
83 |
#ifndef nonstd_lite_HAVE_IN_PLACE_TYPES |
|
84 |
#define nonstd_lite_HAVE_IN_PLACE_TYPES 1 |
|
85 | ||
86 |
// C++17 std::in_place in <utility>: |
|
87 | ||
88 |
#if optional_CPP17_OR_GREATER |
|
89 | ||
90 |
#include <utility> |
|
91 | ||
92 |
namespace nonstd { |
|
93 | ||
94 |
using std::in_place; |
|
95 |
using std::in_place_type; |
|
96 |
using std::in_place_index; |
|
97 |
using std::in_place_t; |
|
98 |
using std::in_place_type_t; |
|
99 |
using std::in_place_index_t; |
|
100 | ||
101 |
#define nonstd_lite_in_place_t( T) std::in_place_t |
|
102 |
#define nonstd_lite_in_place_type_t( T) std::in_place_type_t<T> |
|
103 |
#define nonstd_lite_in_place_index_t(K) std::in_place_index_t<K> |
|
104 | ||
105 |
#define nonstd_lite_in_place( T) std::in_place_t{} |
|
106 |
#define nonstd_lite_in_place_type( T) std::in_place_type_t<T>{} |
|
107 |
#define nonstd_lite_in_place_index(K) std::in_place_index_t<K>{} |
|
108 | ||
109 |
} // namespace nonstd |
|
110 | ||
111 |
#else // optional_CPP17_OR_GREATER |
|
112 | ||
113 |
#include <cstddef> |
|
114 | ||
115 |
namespace nonstd { |
|
116 |
namespace detail { |
|
117 | ||
118 |
template< class T > |
|
119 |
struct in_place_type_tag {}; |
|
120 | ||
121 |
template< std::size_t K > |
|
122 |
struct in_place_index_tag {}; |
|
123 | ||
124 |
} // namespace detail |
|
125 | ||
126 |
struct in_place_t {}; |
|
127 | ||
128 |
template< class T > |
|
129 |
inline in_place_t in_place( detail::in_place_type_tag<T> /*unused*/ = detail::in_place_type_tag<T>() ) |
|
130 |
{ |
|
131 |
return in_place_t(); |
|
132 |
} |
|
133 | ||
134 |
template< std::size_t K > |
|
135 |
inline in_place_t in_place( detail::in_place_index_tag<K> /*unused*/ = detail::in_place_index_tag<K>() ) |
|
136 |
{ |
|
137 |
return in_place_t(); |
|
138 |
} |
|
139 | ||
140 |
template< class T > |
|
141 |
inline in_place_t in_place_type( detail::in_place_type_tag<T> /*unused*/ = detail::in_place_type_tag<T>() ) |
|
142 |
{ |
|
143 |
return in_place_t(); |
|
144 |
} |
|
145 | ||
146 |
template< std::size_t K > |
|
147 |
inline in_place_t in_place_index( detail::in_place_index_tag<K> /*unused*/ = detail::in_place_index_tag<K>() ) |
|
148 |
{ |
|
149 |
return in_place_t(); |
|
150 |
} |
|
151 | ||
152 |
// mimic templated typedef: |
|
153 | ||
154 |
#define nonstd_lite_in_place_t( T) nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T> ) |
|
155 |
#define nonstd_lite_in_place_type_t( T) nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T> ) |
|
156 |
#define nonstd_lite_in_place_index_t(K) nonstd::in_place_t(&)( nonstd::detail::in_place_index_tag<K> ) |
|
157 | ||
158 |
#define nonstd_lite_in_place( T) nonstd::in_place_type<T> |
|
159 |
#define nonstd_lite_in_place_type( T) nonstd::in_place_type<T> |
|
160 |
#define nonstd_lite_in_place_index(K) nonstd::in_place_index<K> |
|
161 | ||
162 |
} // namespace nonstd |
|
163 | ||
164 |
#endif // optional_CPP17_OR_GREATER |
|
165 |
#endif // nonstd_lite_HAVE_IN_PLACE_TYPES |
|
166 | ||
167 |
// |
|
168 |
// Using std::optional: |
|
169 |
// |
|
170 | ||
171 |
#if optional_USES_STD_OPTIONAL |
|
172 | ||
173 |
#include <optional> |
|
174 | ||
175 |
namespace nonstd { |
|
176 | ||
177 |
using std::optional; |
|
178 |
using std::bad_optional_access; |
|
179 |
using std::hash; |
|
180 | ||
181 |
using std::nullopt; |
|
182 |
using std::nullopt_t; |
|
183 | ||
184 |
using std::operator==; |
|
185 |
using std::operator!=; |
|
186 |
using std::operator<; |
|
187 |
using std::operator<=; |
|
188 |
using std::operator>; |
|
189 |
using std::operator>=; |
|
190 |
using std::make_optional; |
|
191 |
using std::swap; |
|
192 |
} |
|
193 | ||
194 |
#else // optional_USES_STD_OPTIONAL |
|
195 | ||
196 |
#include <cassert> |
|
197 |
#include <utility> |
|
198 | ||
199 |
// optional-lite alignment configuration: |
|
200 | ||
201 |
#ifndef optional_CONFIG_MAX_ALIGN_HACK |
|
202 |
# define optional_CONFIG_MAX_ALIGN_HACK 0 |
|
203 |
#endif |
|
204 | ||
205 |
#ifndef optional_CONFIG_ALIGN_AS |
|
206 |
// no default, used in #if defined() |
|
207 |
#endif |
|
208 | ||
209 |
#ifndef optional_CONFIG_ALIGN_AS_FALLBACK |
|
210 |
# define optional_CONFIG_ALIGN_AS_FALLBACK double |
|
211 |
#endif |
|
212 | ||
213 |
// Compiler warning suppression: |
|
214 | ||
215 |
#if defined(__clang__) |
|
216 |
# pragma clang diagnostic push |
|
217 |
# pragma clang diagnostic ignored "-Wundef" |
|
218 |
#elif defined(__GNUC__) |
|
219 |
# pragma GCC diagnostic push |
|
220 |
# pragma GCC diagnostic ignored "-Wundef" |
|
221 |
#elif defined(_MSC_VER ) |
|
222 |
# pragma warning( push ) |
|
223 |
#endif |
|
224 | ||
225 |
// half-open range [lo..hi): |
|
226 |
#define optional_BETWEEN( v, lo, hi ) ( (lo) <= (v) && (v) < (hi) ) |
|
227 | ||
228 |
// Compiler versions: |
|
229 |
// |
|
230 |
// MSVC++ 6.0 _MSC_VER == 1200 (Visual Studio 6.0) |
|
231 |
// MSVC++ 7.0 _MSC_VER == 1300 (Visual Studio .NET 2002) |
|
232 |
// MSVC++ 7.1 _MSC_VER == 1310 (Visual Studio .NET 2003) |
|
233 |
// MSVC++ 8.0 _MSC_VER == 1400 (Visual Studio 2005) |
|
234 |
// MSVC++ 9.0 _MSC_VER == 1500 (Visual Studio 2008) |
|
235 |
// MSVC++ 10.0 _MSC_VER == 1600 (Visual Studio 2010) |
|
236 |
// MSVC++ 11.0 _MSC_VER == 1700 (Visual Studio 2012) |
|
237 |
// MSVC++ 12.0 _MSC_VER == 1800 (Visual Studio 2013) |
|
238 |
// MSVC++ 14.0 _MSC_VER == 1900 (Visual Studio 2015) |
|
239 |
// MSVC++ 14.1 _MSC_VER >= 1910 (Visual Studio 2017) |
|
240 | ||
241 |
#if defined(_MSC_VER ) && !defined(__clang__) |
|
242 |
# define optional_COMPILER_MSVC_VER (_MSC_VER ) |
|
243 |
# define optional_COMPILER_MSVC_VERSION (_MSC_VER / 10 - 10 * ( 5 + (_MSC_VER < 1900 ) ) ) |
|
244 |
#else |
|
245 |
# define optional_COMPILER_MSVC_VER 0 |
|
246 |
# define optional_COMPILER_MSVC_VERSION 0 |
|
247 |
#endif |
|
248 | ||
249 |
#define optional_COMPILER_VERSION( major, minor, patch ) ( 10 * (10 * (major) + (minor) ) + (patch) ) |
|
250 | ||
251 |
#if defined(__GNUC__) && !defined(__clang__) |
|
252 |
# define optional_COMPILER_GNUC_VERSION optional_COMPILER_VERSION(__GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__) |
|
253 |
#else |
|
254 |
# define optional_COMPILER_GNUC_VERSION 0 |
|
255 |
#endif |
|
256 | ||
257 |
#if defined(__clang__) |
|
258 |
# define optional_COMPILER_CLANG_VERSION optional_COMPILER_VERSION(__clang_major__, __clang_minor__, __clang_patchlevel__) |
|
259 |
#else |
|
260 |
# define optional_COMPILER_CLANG_VERSION 0 |
|
261 |
#endif |
|
262 | ||
263 |
#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 140 ) |
|
264 |
# pragma warning( disable: 4345 ) // initialization behavior changed |
|
265 |
#endif |
|
266 | ||
267 |
#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 150 ) |
|
268 |
# pragma warning( disable: 4814 ) // in C++14 'constexpr' will not imply 'const' |
|
269 |
#endif |
|
270 | ||
271 |
// Presence of language and library features: |
|
272 | ||
273 |
#define optional_HAVE(FEATURE) ( optional_HAVE_##FEATURE ) |
|
274 | ||
275 |
#ifdef _HAS_CPP0X |
|
276 |
# define optional_HAS_CPP0X _HAS_CPP0X |
|
277 |
#else |
|
278 |
# define optional_HAS_CPP0X 0 |
|
279 |
#endif |
|
280 | ||
281 |
// Unless defined otherwise below, consider VC14 as C++11 for optional-lite: |
|
282 | ||
283 |
#if optional_COMPILER_MSVC_VER >= 1900 |
|
284 |
# undef optional_CPP11_OR_GREATER |
|
285 |
# define optional_CPP11_OR_GREATER 1 |
|
286 |
#endif |
|
287 | ||
288 |
#define optional_CPP11_90 (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1500) |
|
289 |
#define optional_CPP11_100 (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1600) |
|
290 |
#define optional_CPP11_110 (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1700) |
|
291 |
#define optional_CPP11_120 (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1800) |
|
292 |
#define optional_CPP11_140 (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1900) |
|
293 |
#define optional_CPP11_141 (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1910) |
|
294 | ||
295 |
#define optional_CPP14_000 (optional_CPP14_OR_GREATER) |
|
296 |
#define optional_CPP17_000 (optional_CPP17_OR_GREATER) |
|
297 | ||
298 |
// Presence of C++11 language features: |
|
299 | ||
300 |
#define optional_HAVE_CONSTEXPR_11 optional_CPP11_140 |
|
301 |
#define optional_HAVE_IS_DEFAULT optional_CPP11_140 |
|
302 |
#define optional_HAVE_NOEXCEPT optional_CPP11_140 |
|
303 |
#define optional_HAVE_NULLPTR optional_CPP11_100 |
|
304 |
#define optional_HAVE_REF_QUALIFIER optional_CPP11_140 |
|
305 | ||
306 |
// Presence of C++14 language features: |
|
307 | ||
308 |
#define optional_HAVE_CONSTEXPR_14 optional_CPP14_000 |
|
309 | ||
310 |
// Presence of C++17 language features: |
|
311 | ||
312 |
#define optional_HAVE_NODISCARD optional_CPP17_000 |
|
313 | ||
314 |
// Presence of C++ library features: |
|
315 | ||
316 |
#define optional_HAVE_CONDITIONAL optional_CPP11_120 |
|
317 |
#define optional_HAVE_REMOVE_CV optional_CPP11_120 |
|
318 |
#define optional_HAVE_TYPE_TRAITS optional_CPP11_90 |
|
319 | ||
320 |
#define optional_HAVE_TR1_TYPE_TRAITS (!! optional_COMPILER_GNUC_VERSION ) |
|
321 |
#define optional_HAVE_TR1_ADD_POINTER (!! optional_COMPILER_GNUC_VERSION ) |
|
322 | ||
323 |
// C++ feature usage: |
|
324 | ||
325 |
#if optional_HAVE( CONSTEXPR_11 ) |
|
326 |
# define optional_constexpr constexpr |
|
327 |
#else |
|
328 |
# define optional_constexpr /*constexpr*/ |
|
329 |
#endif |
|
330 | ||
331 |
#if optional_HAVE( IS_DEFAULT ) |
|
332 |
# define optional_is_default = default; |
|
333 |
#else |
|
334 |
# define optional_is_default {} |
|
335 |
#endif |
|
336 | ||
337 |
#if optional_HAVE( CONSTEXPR_14 ) |
|
338 |
# define optional_constexpr14 constexpr |
|
339 |
#else |
|
340 |
# define optional_constexpr14 /*constexpr*/ |
|
341 |
#endif |
|
342 | ||
343 |
#if optional_HAVE( NODISCARD ) |
|
344 |
# define optional_nodiscard [[nodiscard]] |
|
345 |
#else |
|
346 |
# define optional_nodiscard /*[[nodiscard]]*/ |
|
347 |
#endif |
|
348 | ||
349 |
#if optional_HAVE( NOEXCEPT ) |
|
350 |
# define optional_noexcept noexcept |
|
351 |
#else |
|
352 |
# define optional_noexcept /*noexcept*/ |
|
353 |
#endif |
|
354 | ||
355 |
#if optional_HAVE( NULLPTR ) |
|
356 |
# define optional_nullptr nullptr |
|
357 |
#else |
|
358 |
# define optional_nullptr NULL |
|
359 |
#endif |
|
360 | ||
361 |
#if optional_HAVE( REF_QUALIFIER ) |
|
362 |
// NOLINTNEXTLINE( bugprone-macro-parentheses ) |
|
363 |
# define optional_ref_qual & |
|
364 |
# define optional_refref_qual && |
|
365 |
#else |
|
366 |
# define optional_ref_qual /*&*/ |
|
367 |
# define optional_refref_qual /*&&*/ |
|
368 |
#endif |
|
369 | ||
370 |
// additional includes: |
|
371 | ||
372 |
#if optional_CONFIG_NO_EXCEPTIONS |
|
373 |
// already included: <cassert> |
|
374 |
#else |
|
375 |
# include <stdexcept> |
|
376 |
#endif |
|
377 | ||
378 |
#if optional_CPP11_OR_GREATER |
|
379 |
# include <functional> |
|
380 |
#endif |
|
381 | ||
382 |
#if optional_HAVE( INITIALIZER_LIST ) |
|
383 |
# include <initializer_list> |
|
384 |
#endif |
|
385 | ||
386 |
#if optional_HAVE( TYPE_TRAITS ) |
|
387 |
# include <type_traits> |
|
388 |
#elif optional_HAVE( TR1_TYPE_TRAITS ) |
|
389 |
# include <tr1/type_traits> |
|
390 |
#endif |
|
391 | ||
392 |
// Method enabling |
|
393 | ||
394 |
#if optional_CPP11_OR_GREATER |
|
395 | ||
396 |
#define optional_REQUIRES_0(...) \ |
|
397 |
template< bool B = (__VA_ARGS__), typename std::enable_if<B, int>::type = 0 > |
|
398 | ||
399 |
#define optional_REQUIRES_T(...) \ |
|
400 |
, typename = typename std::enable_if< (__VA_ARGS__), nonstd::optional_lite::detail::enabler >::type |
|
401 | ||
402 |
#define optional_REQUIRES_R(R, ...) \ |
|
403 |
typename std::enable_if< (__VA_ARGS__), R>::type |
|
404 | ||
405 |
#define optional_REQUIRES_A(...) \ |
|
406 |
, typename std::enable_if< (__VA_ARGS__), void*>::type = nullptr |
|
407 | ||
408 |
#endif |
|
409 | ||
410 |
// |
|
411 |
// optional: |
|
412 |
// |
|
413 | ||
414 |
namespace nonstd { namespace optional_lite { |
|
415 | ||
416 |
namespace std11 { |
|
417 | ||
418 |
#if optional_CPP11_OR_GREATER |
|
419 |
using std::move; |
|
420 |
#else |
|
421 |
template< typename T > T & move( T & t ) { return t; } |
|
422 |
#endif |
|
423 | ||
424 |
#if optional_HAVE( CONDITIONAL ) |
|
425 |
using std::conditional; |
|
426 |
#else |
|
427 |
template< bool B, typename T, typename F > struct conditional { typedef T type; }; |
|
428 |
template< typename T, typename F > struct conditional<false, T, F> { typedef F type; }; |
|
429 |
#endif // optional_HAVE_CONDITIONAL |
|
430 | ||
431 |
} // namespace std11 |
|
432 | ||
433 |
#if optional_CPP11_OR_GREATER |
|
434 | ||
435 |
/// type traits C++17: |
|
436 | ||
437 |
namespace std17 { |
|
438 | ||
439 |
#if optional_CPP17_OR_GREATER |
|
440 | ||
441 |
using std::is_swappable; |
|
442 |
using std::is_nothrow_swappable; |
|
443 | ||
444 |
#elif optional_CPP11_OR_GREATER |
|
445 | ||
446 |
namespace detail { |
|
447 | ||
448 |
using std::swap; |
|
449 | ||
450 |
struct is_swappable |
|
451 |
{ |
|
452 |
template< typename T, typename = decltype( swap( std::declval<T&>(), std::declval<T&>() ) ) > |
|
453 |
static std::true_type test( int /*unused*/ ); |
|
454 | ||
455 |
template< typename > |
|
456 |
static std::false_type test(...); |
|
457 |
}; |
|
458 | ||
459 |
struct is_nothrow_swappable |
|
460 |
{ |
|
461 |
// wrap noexcept(expr) in separate function as work-around for VC140 (VS2015): |
|
462 | ||
463 |
template< typename T > |
|
464 |
static constexpr bool satisfies() |
|
465 |
{ |
|
466 |
return noexcept( swap( std::declval<T&>(), std::declval<T&>() ) ); |
|
467 |
} |
|
468 | ||
469 |
template< typename T > |
|
470 |
static auto test( int /*unused*/ ) -> std::integral_constant<bool, satisfies<T>()>{} |
|
471 | ||
472 |
template< typename > |
|
473 |
static auto test(...) -> std::false_type; |
|
474 |
}; |
|
475 | ||
476 |
} // namespace detail |
|
477 | ||
478 |
// is [nothow] swappable: |
|
479 | ||
480 |
template< typename T > |
|
481 |
struct is_swappable : decltype( detail::is_swappable::test<T>(0) ){}; |
|
482 | ||
483 |
template< typename T > |
|
484 |
struct is_nothrow_swappable : decltype( detail::is_nothrow_swappable::test<T>(0) ){}; |
|
485 | ||
486 |
#endif // optional_CPP17_OR_GREATER |
|
487 | ||
488 |
} // namespace std17 |
|
489 | ||
490 |
/// type traits C++20: |
|
491 | ||
492 |
namespace std20 { |
|
493 | ||
494 |
template< typename T > |
|
495 |
struct remove_cvref |
|
496 |
{ |
|
497 |
typedef typename std::remove_cv< typename std::remove_reference<T>::type >::type type; |
|
498 |
}; |
|
499 | ||
500 |
} // namespace std20 |
|
501 | ||
502 |
#endif // optional_CPP11_OR_GREATER |
|
503 | ||
504 |
/// class optional |
|
505 | ||
506 |
template< typename T > |
|
507 |
class optional; |
|
508 | ||
509 |
namespace detail { |
|
510 | ||
511 |
// for optional_REQUIRES_T |
|
512 | ||
513 |
#if optional_CPP11_OR_GREATER |
|
514 |
enum class enabler{}; |
|
515 |
#endif |
|
516 | ||
517 |
// C++11 emulation: |
|
518 | ||
519 |
struct nulltype{}; |
|
520 | ||
521 |
template< typename Head, typename Tail > |
|
522 |
struct typelist |
|
523 |
{ |
|
524 |
typedef Head head; |
|
525 |
typedef Tail tail; |
|
526 |
}; |
|
527 | ||
528 |
#if optional_CONFIG_MAX_ALIGN_HACK |
|
529 | ||
530 |
// Max align, use most restricted type for alignment: |
|
531 | ||
532 |
#define optional_UNIQUE( name ) optional_UNIQUE2( name, __LINE__ ) |
|
533 |
#define optional_UNIQUE2( name, line ) optional_UNIQUE3( name, line ) |
|
534 |
#define optional_UNIQUE3( name, line ) name ## line |
|
535 | ||
536 |
#define optional_ALIGN_TYPE( type ) \ |
|
537 |
type optional_UNIQUE( _t ); struct_t< type > optional_UNIQUE( _st ) |
|
538 | ||
539 |
template< typename T > |
|
540 |
struct struct_t { T _; }; |
|
541 | ||
542 |
union max_align_t |
|
543 |
{ |
|
544 |
optional_ALIGN_TYPE( char ); |
|
545 |
optional_ALIGN_TYPE( short int ); |
|
546 |
optional_ALIGN_TYPE( int ); |
|
547 |
optional_ALIGN_TYPE( long int ); |
|
548 |
optional_ALIGN_TYPE( float ); |
|
549 |
optional_ALIGN_TYPE( double ); |
|
550 |
optional_ALIGN_TYPE( long double ); |
|
551 |
optional_ALIGN_TYPE( char * ); |
|
552 |
optional_ALIGN_TYPE( short int * ); |
|
553 |
optional_ALIGN_TYPE( int * ); |
|
554 |
optional_ALIGN_TYPE( long int * ); |
|
555 |
optional_ALIGN_TYPE( float * ); |
|
556 |
optional_ALIGN_TYPE( double * ); |
|
557 |
optional_ALIGN_TYPE( long double * ); |
|
558 |
optional_ALIGN_TYPE( void * ); |
|
559 | ||
560 |
#ifdef HAVE_LONG_LONG |
|
561 |
optional_ALIGN_TYPE( long long ); |
|
562 |
#endif |
|
563 | ||
564 |
struct Unknown; |
|
565 | ||
566 |
Unknown ( * optional_UNIQUE(_) )( Unknown ); |
|
567 |
Unknown * Unknown::* optional_UNIQUE(_); |
|
568 |
Unknown ( Unknown::* optional_UNIQUE(_) )( Unknown ); |
|
569 | ||
570 |
struct_t< Unknown ( * )( Unknown) > optional_UNIQUE(_); |
|
571 |
struct_t< Unknown * Unknown::* > optional_UNIQUE(_); |
|
572 |
struct_t< Unknown ( Unknown::* )(Unknown) > optional_UNIQUE(_); |
|
573 |
}; |
|
574 | ||
575 |
#undef optional_UNIQUE |
|
576 |
#undef optional_UNIQUE2 |
|
577 |
#undef optional_UNIQUE3 |
|
578 | ||
579 |
#undef optional_ALIGN_TYPE |
|
580 | ||
581 |
#elif defined( optional_CONFIG_ALIGN_AS ) // optional_CONFIG_MAX_ALIGN_HACK |
|
582 | ||
583 |
// Use user-specified type for alignment: |
|
584 | ||
585 |
#define optional_ALIGN_AS( unused ) \ |
|
586 |
optional_CONFIG_ALIGN_AS |
|
587 | ||
588 |
#else // optional_CONFIG_MAX_ALIGN_HACK |
|
589 | ||
590 |
// Determine POD type to use for alignment: |
|
591 | ||
592 |
#define optional_ALIGN_AS( to_align ) \ |
|
593 |
typename type_of_size< alignment_types, alignment_of< to_align >::value >::type |
|
594 | ||
595 |
template< typename T > |
|
596 |
struct alignment_of; |
|
597 | ||
598 |
template< typename T > |
|
599 |
struct alignment_of_hack |
|
600 |
{ |
|
601 |
char c; |
|
602 |
T t; |
|
603 |
alignment_of_hack(); |
|
604 |
}; |
|
605 | ||
606 |
template< size_t A, size_t S > |
|
607 |
struct alignment_logic |
|
608 |
{ |
|
609 |
enum { value = A < S ? A : S }; |
|
610 |
}; |
|
611 | ||
612 |
template< typename T > |
|
613 |
struct alignment_of |