1 |
# Copyright (c) 2018-2022 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 |
.resultColNames <- function(...) { |
|
15 | 108x |
.resultColumns(...)$colname |
16 |
} |
|
17 | ||
18 |
#' Return column names to be used for summary operations |
|
19 |
#' |
|
20 |
#' @param value_names names of value raster layers |
|
21 |
#' @param weight_names names of weighting raster layers |
|
22 |
#' @param fun functions or names of summary operations |
|
23 |
#' @param full_colnames return a complete column name even when there is no |
|
24 |
#' ambiguity? |
|
25 |
#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile} |
|
26 |
#' @return character vector of column names |
|
27 |
#' @keywords internal |
|
28 |
.resultColumns <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric(), unique_values=numeric(), colname_fun = NULL) { |
|
29 | 246x |
if (inherits(fun, 'standardGeneric')) { |
30 | 22x |
stat_names <- fun@generic[1] |
31 | 224x |
} else if (is.function(fun)) { |
32 | 32x |
stat_names <- 'fun' |
33 | 192x |
} else if (is.null(fun)) { |
34 | 42x |
stat_names <- '' |
35 |
} else { |
|
36 | 150x |
stat_names <- fun |
37 |
} |
|
38 | ||
39 | 246x |
quantile_index = which(stat_names == 'quantile') |
40 | 246x |
if (length(quantile_index) != 0) { |
41 | 6x |
stat_names <- splice(stat_names, quantile_index, rep('quantile', length(quantiles))) |
42 |
} |
|
43 | ||
44 | 246x |
for (stat in c('frac', 'weighted_frac')) { |
45 | 492x |
frac_index = which(stat_names == stat) |
46 | 492x |
if (length(frac_index) != 0) { |
47 | 6x |
stat_names <- splice(stat_names, frac_index, rep(stat, length(unique_values))) |
48 |
} |
|
49 |
} |
|
50 | ||
51 | 246x |
ind <- .valueWeightIndexes(length(value_names), length(weight_names)) |
52 | 246x |
vn <- value_names[ind$values] |
53 | 246x |
wn <- weight_names[ind$weights] |
54 | ||
55 |
# determine all combinations of index and stat |
|
56 | 246x |
z <- expand.grid(index=seq_along(vn), |
57 | 246x |
stat_name=stat_names, stringsAsFactors=FALSE) |
58 | ||
59 | 246x |
z$values <- vn[z$index] |
60 | 246x |
z$base_value <- NA_real_ |
61 | ||
62 | 246x |
for (stat in c('frac', 'weighted_frac')) { |
63 | 492x |
ifrac <- which(z$stat_name == stat) |
64 | 492x |
z$base_value[ifrac] <- rep(unique_values, each = length(ifrac) / length(unique_values)) |
65 |
} |
|
66 | ||
67 | 246x |
iquantile <- which(z$stat_name == 'quantile') |
68 | 246x |
z$base_value[iquantile] <- rep(quantiles, each = length(iquantile) / length(quantiles)) |
69 | ||
70 | 246x |
if (is.null(wn)) { |
71 | 160x |
z$weights <- NA |
72 |
} else { |
|
73 | 86x |
z$weights <- wn[z$index] |
74 |
} |
|
75 | 246x |
z$weights[!.includeWeightInColName(z$stat_name)] <- NA |
76 | ||
77 | 246x |
if (is.null(colname_fun)) { |
78 | 245x |
colname_fun <- function(...) { |
79 | 387x |
.makeColname(full_colnames = full_colnames, ...) |
80 |
} |
|
81 |
} |
|
82 | ||
83 | 246x |
z$colname <- mapply(colname_fun, |
84 | 246x |
fun_name = z$stat_name, |
85 | 246x |
values = z$values, |
86 | 246x |
weights = z$weights, |
87 | 246x |
fun_value = z$base_value, |
88 | 246x |
MoreArgs = list(nvalues = length(value_names), |
89 | 246x |
nweights = length(weight_names)), |
90 | 246x |
USE.NAMES = FALSE) |
91 | ||
92 | 246x |
return(z) |
93 |
} |
|
94 | ||
95 |
.makeColname <- function(fun_name, values, weights, fun_value, full_colnames, nvalues, nweights) { |
|
96 |
# construct column names for each index, stat |
|
97 |
# add weight layer name only if layer is ambiguously weighted |
|
98 | 387x |
if (fun_name == 'quantile') { |
99 | 13x |
fun_component <- sprintf('q%02d', as.integer(100 * fun_value)) |
100 | 374x |
} else if (fun_name %in% c('frac', 'weighted_frac')) { |
101 | 26x |
fun_component <- sprintf('%s_%s', fun_name, fun_value) |
102 |
} else { |
|
103 | 348x |
fun_component <- fun_name |
104 |
} |
|
105 | ||
106 | 387x |
ret <- fun_component |
107 | 387x |
if (full_colnames || nvalues > 1) { |
108 | 137x |
ret <- paste(ret, values, sep='.') |
109 |
} |
|
110 | 387x |
if ((!is.na(weights)) && ((full_colnames & nweights > 0) || nweights > 1)) { |
111 | 54x |
ret <- paste(ret, weights, sep='.') |
112 |
} |
|
113 | ||
114 | 387x |
return(ret) |
115 |
} |
|
116 | ||
117 |
.includeWeightInColName <- function(fun) { |
|
118 | 246x |
.isWeighted(fun) | fun == 'fun' |
119 |
} |
|
120 | ||
121 |
.isWeighted <- function(stat_name) { |
|
122 | 334x |
stat_name %in% c('weighted_mean', 'weighted_sum', 'weighted_frac') |
123 |
} |
|
124 | ||
125 |
#' Compute indexes for the value and weight layers that should be |
|
126 |
#' processed together |
|
127 |
#' |
|
128 |
#' @param num_values number of layers in value raster |
|
129 |
#' @param num_weights number of layers in weighting raster |
|
130 |
#' @return list with \code{values} and \code{weights} elements |
|
131 |
#' providing layer indexes |
|
132 |
#' @keywords internal |
|
133 |
.valueWeightIndexes <- function(num_values, num_weights) { |
|
134 | 330x |
if (num_weights == 0) { |
135 | 223x |
vi <- seq_len(num_values) |
136 | 223x |
wi <- NA |
137 | 107x |
} else if (num_values == num_weights) { |
138 |
# process in parallel |
|
139 | 90x |
vi <- seq_len(num_values) |
140 | 90x |
wi <- seq_len(num_weights) |
141 | 17x |
} else if (num_values == 1 && num_weights > 1) { |
142 |
# recycle values |
|
143 | 8x |
vi <- rep.int(1, num_weights) |
144 | 8x |
wi <- seq_len(num_weights) |
145 | 9x |
} else if (num_values > 1 && num_weights == 1) { |
146 |
# recycle weights |
|
147 | 9x |
vi <- seq_len(num_values) |
148 | 9x |
wi <- rep.int(1, num_values) |
149 |
} |
|
150 | ||
151 | 330x |
list(values = vi, weights = wi) |
152 |
} |
|
153 | ||
154 |
.areaMethod <- function(crs_obj) { |
|
155 | 39x |
if (!(is.na(crs_obj)) && crs_obj$units_gdal == 'degree') { |
156 | 14x |
return('spherical') |
157 |
} else { |
|
158 | 25x |
return('cartesian') |
159 |
} |
|
160 |
} |
|
161 | ||
162 |
.validateFlag <- function(value, name) { |
|
163 | 2762x |
if(!(is.logical(value) && length(value) == 1 && !is.na(value))) { |
164 | 27x |
stop(name, ' must be TRUE or FALSE') |
165 |
} |
|
166 |
} |
|
167 | ||
168 |
.validateNumericScalar <- function(value, name) { |
|
169 | 576x |
if (!(is.numeric(value) && length(value) == 1 && !is.na(value))) { |
170 | 4x |
stop(name, ' must be a single numeric value') |
171 |
} |
|
172 |
} |
|
173 | ||
174 |
.validateNumericScalarOrNA <- function(value, name) { |
|
175 | 568x |
if (!(is.numeric(value) && length(value) == 1)) { |
176 | 4x |
stop(name, ' must be a single numeric value') |
177 |
} |
|
178 |
} |
|
179 | ||
180 |
.validateUniqueNames <- function(x) { |
|
181 | 563x |
nm <- names(x) |
182 | 563x |
if (!is.null(nm)) { |
183 | 370x |
if (length(nm) != length(unique(nm))) { |
184 | 1x |
stop('names of input rasters must be unique') |
185 |
} |
|
186 |
} |
|
187 |
} |
|
188 | ||
189 |
# faster replacement for as.data.frame when input is a named list |
|
190 |
# with equal-length columns |
|
191 |
# from Advanced R, sec. 24.4.2 |
|
192 |
.quickDf <- function(lst) { |
|
193 | 162x |
class(lst) <- 'data.frame' |
194 | 162x |
attr(lst, 'row.names') <- .set_row_names(length(lst[[1]])) |
195 | 162x |
lst |
196 |
} |
|
197 | ||
198 |
.singleColumnToVector <- function(df) { |
|
199 | 14x |
if (ncol(df) == 1) { |
200 | 3x |
df[, 1] |
201 |
} else { |
|
202 | 11x |
df |
203 |
} |
|
204 |
} |
|
205 | ||
206 |
# Return the number of standard (non-...) arguments in a supplied function that |
|
207 |
# do not have a default value. This is used to fail if the summary function |
|
208 |
# provided by the user cannot accept arguments of values and weights. |
|
209 |
.num_expected_args <- function(fun) { |
|
210 | 70x |
a <- formals(args(fun)) |
211 | 70x |
a <- a[names(a) != '...'] |
212 | 70x |
sum(sapply(a, nchar) == 0) |
213 |
} |
|
214 | ||
215 |
.startReading <- function(r) { |
|
216 | 521x |
if(inherits(r, 'BasicRaster')) { |
217 | 317x |
return(raster::readStart(r)) |
218 | 204x |
} else if (inherits(r, 'SpatRaster')) { |
219 | 35x |
terra::readStart(r) |
220 |
} |
|
221 | ||
222 | 204x |
return(r) |
223 |
} |
|
224 | ||
225 |
.stopReading <- function(r) { |
|
226 | 521x |
if(inherits(r, 'BasicRaster')) { |
227 | 317x |
return(raster::readStop(r)) |
228 | 204x |
} else if (inherits(r, 'SpatRaster')) { |
229 | 35x |
terra::readStop(r) |
230 |
} |
|
231 | ||
232 | 204x |
return(r) |
233 |
} |
|
234 | ||
235 |
.crs <- function(r) { |
|
236 | 12x |
if(inherits(r, 'BasicRaster')) { |
237 | 12x |
if (utils::packageVersion('raster') < numeric_version('3.5')) { |
238 | ! |
return(raster::crs(r)) |
239 |
} else { |
|
240 | 12x |
return(terra::crs(r)) |
241 |
} |
|
242 | ! |
} else if (inherits(r, 'SpatRaster')) { |
243 | ! |
return(terra::crs(r)) |
244 |
} else { |
|
245 | ! |
stop('Unknown type: ', class(r)) |
246 |
} |
|
247 |
} |
|
248 | ||
249 |
.setValues <- function(r, x) { |
|
250 | 13x |
if(inherits(r, 'BasicRaster')) { |
251 | 13x |
if (utils::packageVersion('raster') < numeric_version('3.5')) { |
252 | ! |
raster::values(r) <- x |
253 |
} else { |
|
254 | 13x |
terra::values(r) <- x |
255 |
} |
|
256 | ! |
} else if (inherits(r, 'SpatRaster')) { |
257 | ! |
raster::values(r) <- x |
258 |
} else { |
|
259 | ! |
stop('Unknown type: ', class(r)) |
260 |
} |
|
261 | ||
262 | 13x |
return(r) |
263 |
} |
|
264 | ||
265 |
.numLayers <- function(r) { |
|
266 | 640x |
if(inherits(r, 'BasicRaster')) { |
267 | 520x |
return(raster::nlayers(r)) |
268 | 120x |
} else if (inherits(r, 'SpatRaster')) { |
269 | 120x |
return(terra::nlyr(r)) |
270 |
} else { |
|
271 | ! |
stop('Unknown type: ', class(r)) |
272 |
} |
|
273 |
} |
|
274 | ||
275 |
.isRaster <- function(r) { |
|
276 | 191x |
inherits(r, 'BasicRaster') | inherits(r, 'SpatRaster') |
277 |
} |
|
278 | ||
279 |
.xFromCol <- function(r, col) { |
|
280 | 30x |
if (inherits(r, 'BasicRaster')) { |
281 | 17x |
raster::xFromCol(r, col) |
282 | 13x |
} else if (inherits(r, 'SpatRaster')) { |
283 | 13x |
terra::xFromCol(r, col) |
284 |
} else { |
|
285 | ! |
stop('Unknown type: ', class(r)) |
286 |
} |
|
287 |
} |
|
288 | ||
289 |
.colFromX <- function(r, x) { |
|
290 | 78x |
if (inherits(r, 'BasicRaster')) { |
291 | 30x |
raster::colFromX(r, x) |
292 | 48x |
} else if (inherits(r, 'SpatRaster')) { |
293 | 48x |
terra::colFromX(r, x) |
294 |
} else { |
|
295 | ! |
stop('Unknown type: ', class(r)) |
296 |
} |
|
297 |
} |
|
298 | ||
299 |
.yFromRow <- function(r, row) { |
|
300 | 30x |
if (inherits(r, 'BasicRaster')) { |
301 | 17x |
raster::yFromRow(r, row) |
302 | 13x |
} else if (inherits(r, 'SpatRaster')) { |
303 | 13x |
terra::yFromRow(r, row) |
304 |
} else { |
|
305 | ! |
stop('Unknown type: ', class(r)) |
306 |
} |
|
307 |
} |
|
308 | ||
309 |
.rowFromY <- function(r, y) { |
|
310 | 78x |
if (inherits(r, 'BasicRaster')) { |
311 | 30x |
raster::rowFromY(r, y) |
312 | 48x |
} else if (inherits(r, 'SpatRaster')) { |
313 | 48x |
terra::rowFromY(r, y) |
314 |
} else { |
|
315 | ! |
stop('Unknown type: ', class(r)) |
316 |
} |
|
317 |
} |
|
318 | ||
319 |
.cellFromRowCol <- function(r, row, col) { |
|
320 | 24x |
if (inherits(r, 'BasicRaster')) { |
321 | 11x |
raster::cellFromRowCol(r, row, col) |
322 | 13x |
} else if (inherits(r, 'SpatRaster')) { |
323 | 13x |
terra::cellFromRowCol(r, row, col) |
324 |
} else { |
|
325 | ! |
stop('Unknown type: ', class(r)) |
326 |
} |
|
327 |
} |
|
328 | ||
329 |
.extent <- function(r) { |
|
330 | 724x |
if (inherits(r, 'BasicRaster')) { |
331 | 572x |
ex <- r@extent |
332 | 572x |
c(ex@xmin, ex@ymin, ex@xmax, ex@ymax) |
333 | 152x |
} else if (inherits(r, 'SpatRaster')) { |
334 | 152x |
ex <- terra::ext(r) |
335 | 152x |
c(ex$xmin, ex$ymin, ex$xmax, ex$ymax) |
336 |
} else { |
|
337 | ! |
stop('Unknown type: ', class(r)) |
338 |
} |
|
339 |
} |
|
340 | ||
341 |
.res <- function(r) { |
|
342 | 724x |
if (inherits(r, 'BasicRaster')) { |
343 | 572x |
raster::res(r) |
344 | 152x |
} else if (inherits(r, 'SpatRaster')) { |
345 | 152x |
terra::res(r) |
346 |
} else { |
|
347 | ! |
stop('Unknown type: ', class(r)) |
348 |
} |
|
349 |
} |
|
350 | ||
351 |
.getValuesBlock <- function(r, row, nrows, col, ncols) { |
|
352 | 1054x |
if (inherits(r, 'BasicRaster')) { |
353 | 951x |
raster::getValuesBlock(r, row, nrows, col, ncols, format = 'm') |
354 | 103x |
} else if (inherits(r, 'SpatRaster')) { |
355 | 103x |
terra::readValues(r, row, nrows, col, ncols, mat = TRUE) |
356 |
} else { |
|
357 | ! |
stop('Unknown type: ', class(r)) |
358 |
} |
|
359 |
} |
|
360 | ||
361 |
.createProgress <- function(progress, n) { |
|
362 | 252x |
if (progress && n > 1) { |
363 | 8x |
pb <- utils::txtProgressBar(min = 0, max = n, initial=0, style=3) |
364 | 8x |
update_progress <- function() { |
365 | 76x |
i <- 1 + utils::getTxtProgressBar(pb) |
366 | 76x |
utils::setTxtProgressBar(pb, i) |
367 | 76x |
if (i == n) { |
368 | 8x |
close(pb) |
369 |
} |
|
370 |
} |
|
371 |
} else { |
|
372 | 244x |
update_progress <- function() {} |
373 |
} |
|
374 | ||
375 | 252x |
return(update_progress) |
376 |
} |
|
377 | ||
378 |
.isInMemory <- function(r) { |
|
379 | 297x |
if (inherits(r, 'BasicRaster')) { |
380 | 272x |
return(raster::inMemory(r)) |
381 | 25x |
} else if (inherits(r, 'SpatRaster')) { |
382 | 25x |
return(terra::inMemory(r[[1]])) |
383 |
} else { |
|
384 | ! |
stop('Unknown type: ', class(r)) |
385 |
} |
|
386 |
} |
|
387 | ||
388 |
.netCDFBlockSize <- function(fname, varname) { |
|
389 | 1x |
nc <- NULL |
390 | 1x |
sz <- NA |
391 | ||
392 | 1x |
tryCatch({ |
393 | 1x |
nc <- ncdf4::nc_open(fname) |
394 | 1x |
sz <- nc$var[[varname]]$chunksizes |
395 | ||
396 | 1x |
dim_index <- nc$var[[varname]]$dimids + 1L |
397 | 1x |
dim_names <- sapply(dim_index, function(i) nc$dim[[i]]$name) |
398 | ||
399 | 1x |
if (all(is.na(sz))) { |
400 |
# file is not compressed |
|
401 | ! |
sz <- rep.int(1, length(dim_names)) |
402 |
} |
|
403 | ||
404 | 1x |
names(sz) <- dim_names |
405 | 1x |
}, finally = { |
406 | 1x |
if (!is.null(nc)) { |
407 | 1x |
ncdf4::nc_close(nc) |
408 |
} |
|
409 |
}) |
|
410 | ||
411 |
# flip dimensions 1 and 2 so we return row/col |
|
412 | 1x |
return(sz[c(2, 1, seq_along(sz)[-(1:2)])]) |
413 |
} |
|
414 | ||
415 |
.blockSize <- function(r) { |
|
416 |
# set default return value in case file is uncompressed and has |
|
417 |
# has no block size, or we simply can't figure it out |
|
418 | 8x |
ret <- c(1, 1) |
419 | ||
420 | 8x |
if (inherits(r, 'BasicRaster')) { |
421 | 2x |
if (r[[1]]@file@driver == 'netcdf') { |
422 | 1x |
ret <- .netCDFBlockSize(r[[1]]@file@name, attr(r[[1]]@data, 'zvar')) |
423 | 1x |
} else if (r@file@driver == 'gdal') { |
424 | 1x |
ret <- c(r@file@blockrows, r@file@blockcols) |
425 |
} |
|
426 | 6x |
} else if (inherits(r, 'SpatRaster')) { |
427 | 6x |
ret <- terra::fileBlocksize(r)[1, ] |
428 |
} |
|
429 | ||
430 | 8x |
unname(ret) |
431 |
} |
|
432 | ||
433 |
.eagerLoad <- function(r, geoms, max_cells_in_memory, message_on_fail) { |
|
434 | 17x |
if (is.null(r)) { |
435 | 5x |
return(NULL) |
436 |
} |
|
437 | ||
438 | 12x |
cells_required <- .numCells(r, geoms) |
439 | ||
440 | 12x |
if (cells_required <= max_cells_in_memory) { |
441 | 6x |
box <- sf::st_bbox(geoms) |
442 | 6x |
geom_ext <- terra::ext(box[c('xmin', 'xmax', 'ymin', 'ymax')]) |
443 | ||
444 | 6x |
if (!inherits(r, 'SpatRaster')) { |
445 |
# current CRAN version of terra (1.4-22) does not preserve |
|
446 |
# names on conversion (https://github.com/rspatial/terra/issues/430) |
|
447 | ! |
nm <- names(r) |
448 | ! |
r <- terra::rast(r) |
449 | ! |
names(r) <- nm |
450 |
} |
|
451 | ||
452 | 6x |
overlap_ext <- terra::intersect(terra::ext(r), geom_ext) |
453 | 6x |
if (is.null(overlap_ext)) { |
454 |
# Extents do not overlap, and terra::crop will throw an error |
|
455 |
# if we try to crop. Return the input raster as-is; nothing will be |
|
456 |
# read from it anyway. |
|
457 | ! |
return(r) |
458 |
} |
|
459 | ||
460 | 6x |
r <- terra::crop(r, geom_ext, snap = 'out') |
461 | 6x |
} else if (message_on_fail) { |
462 | 3x |
message('Cannot preload entire working area of ', cells_required, |
463 | 3x |
' cells with max_cells_in_memory = ', max_cells_in_memory, '.', |
464 | 3x |
' Raster values will be read for each feature individually.') |
465 |
} |
|
466 | ||
467 | 12x |
return(r) |
468 |
} |
|
469 | ||
470 |
.numCells <- function(r, g) { |
|
471 | 12x |
if (is.null(r)) { |
472 | ! |
return(0) |
473 |
} |
|
474 | ||
475 | 12x |
box <- sf::st_bbox(g) |
476 | ||
477 | 12x |
top <- .rowFromY(r, box['ymax']) |
478 | 12x |
bottom <- .rowFromY(r, box['ymin']) |
479 | 12x |
left <- .colFromX(r, box['xmin']) |
480 | 12x |
right <- .colFromX(r, box['xmax']) |
481 | ||
482 | 12x |
if (is.na(top) && is.na(bottom)) { |
483 | ! |
return(0L) |
484 |
} |
|
485 | 12x |
if (is.na(left) && is.na(right)) { |
486 | ! |
return(0L) |
487 |
} |
|
488 | 12x |
if (is.na(top)) { |
489 | 1x |
top <- 1 |
490 |
} |
|
491 | 12x |
if (is.na(bottom)) { |
492 | ! |
bottom <- nrow(r) |
493 |
} |
|
494 | 12x |
if (is.na(left)) { |
495 | ! |
left <- 1 |
496 |
} |
|
497 | 12x |
if (is.na(right)) { |
498 | 1x |
right <- ncol(r) |
499 |
} |
|
500 | ||
501 | 12x |
return( (bottom - top + 1) * (right - left + 1) * .numLayers(r) ) |
502 |
} |
|
503 | ||
504 |
splice <- function(x, i, replacement) { |
|
505 | 12x |
c(x[seq_along(x) < i], |
506 | 12x |
replacement, |
507 | 12x |
x[seq_along(x) > i]) |
508 |
} |
1 |
# Copyright (c) 2018-2022 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 | 321x |
setGeneric("exact_extract", function(x, y, ...) |
15 | 321x |
standardGeneric("exact_extract")) |
16 | ||
17 |
#' Extract or summarize values from rasters |
|
18 |
#' |
|
19 |
#' Extracts the values of cells in a raster (`RasterLayer`, `RasterStack` |
|
20 |
#' `RasterBrick`, or `SpatRaster`) that are covered by polygons in a |
|
21 |
#' simple feature collection (`sf` or `sfc`) or `SpatialPolygonsDataFrame`. |
|
22 |
#' Returns either a summary of the extracted values or the extracted values |
|
23 |
#' themselves. |
|
24 |
#' |
|
25 |
#' @details |
|
26 |
#' `exact_extract` extracts the values of cells in a raster that are covered |
|
27 |
#' by polygonal features in a simple feature collection (`sf` or `sfc`) or |
|
28 |
#' `SpatialPolygonDataFrame`, as well as the fraction or area of each cell that |
|
29 |
#' is covered by the feature. Pixels covered by all parts of the polygon are |
|
30 |
#' considered. If an (invalid) multipart polygon covers the same pixels more |
|
31 |
#' than once, the pixel may have a coverage fraction greater than one. |
|
32 |
#' |
|
33 |
#' The function can either return pixel values directly to the caller, or can |
|
34 |
#' return the result of a predefined summary operation or user-defined R |
|
35 |
#' function applied to the values. These three approaches are described in the |
|
36 |
#' subsections below. |
|
37 |
#' |
|
38 |
#' ## Returning extracted values directly |
|
39 |
#' |
|
40 |
#' If `fun` is not specified, `exact_extract` will return a list with |
|
41 |
#' one data frame for each feature in the input feature collection. The data |
|
42 |
#' frame will contain a column with cell values from each layer in the input |
|
43 |
#' raster (and optional weighting raster) and a column indicating |
|
44 |
#' the fraction or area of the cell that is covered by the polygon. |
|
45 |
#' |
|
46 |
#' If the input rasters have only one layer, the value and weight columns in the |
|
47 |
#' data frame will be named `values` or `weights`. When the input rasters have |
|
48 |
#' more than one layer, the columns will be named according to `names(x)` and |
|
49 |
#' `names(weights)`. The column containing pixel coverage will be called |
|
50 |
#' `coverage_fraction` when `coverage_area = FALSE`, or `coverage_area` when |
|
51 |
#' `coverage_area = TRUE`. Additional columns can be added to the returned data |
|
52 |
#' frames with the `include_area`, `include_cell`, and `include_xy` arguments. |
|
53 |
#' |
|
54 |
#' If the output data frames for multiple features are to be combined (e.g., |
|
55 |
#' with `rbind`), it may be useful to include identifying column(s) from the |
|
56 |
#' input features in the returned data frames using `include_cols`. |
|
57 |
#' |
|
58 |
#' ## Predefined summary operations |
|
59 |
#' |
|
60 |
#' Often the individual pixel values are not needed; only one or more summary |
|
61 |
#' statistics (e.g., mean, sum) is required for each feature. Common summary |
|
62 |
#' statistics can be calculated by `exact_extract` directly using a predefined |
|
63 |
#' summary operation. Where possible, this approach is advantageous because it |
|
64 |
#' allows the package to calculate the statistics incrementally, avoiding the |
|
65 |
#' need to store all pixel values in memory at the same time. This allows the |
|
66 |
#' package to process arbitrarily large data with a small amount of memory. (The |
|
67 |
#' `max_pixels_in_memory` argument can be used to fine-tune the amount of memory |
|
68 |
#' made available to `exact_extract`.) |
|
69 |
#' |
|
70 |
#' To summarize pixel values using a predefined summary option, `fun` should be |
|
71 |
#' set to a character vector of one or more operation names. If the input raster |
|
72 |
#' has a single layer and a single summary operation is specified, |
|
73 |
#' `exact_extract` will return a vector with the result of the summary operation |
|
74 |
#' for each feature in the input. If the input raster has multiple layers, or if |
|
75 |
#' multiple summary operations are specified, `exact_extract` will return a data |
|
76 |
#' frame with a row for each feature and a column for each summary operation / |
|
77 |
#' layer combination. (The `force_df` option can be used to always return a data |
|
78 |
#' frame instead of a vector.) |
|
79 |
#' |
|
80 |
#' The following summary operations are supported: |
|
81 |
#' |
|
82 |
#' * `min` - the minimum non-`NA` value in any raster cell wholly or |
|
83 |
#' partially covered by the polygon |
|
84 |
#' * `max` - the maximum non-`NA` value in any raster cell wholly or |
|
85 |
#' partially covered by the polygon |
|
86 |
#' * `count` - the sum of fractions of raster cells with non-`NA` |
|
87 |
#' values covered by the polygon |
|
88 |
#' * `sum` - the sum of non-`NA` raster cell values, multiplied by |
|
89 |
#' the fraction of the cell that is covered by the polygon |
|
90 |
#' * `mean` - the mean cell value, weighted by the fraction of each cell |
|
91 |
#' that is covered by the polygon |
|
92 |
#' * `median` - the median cell value, weighted by the fraction of each cell |
|
93 |
#' that is covered by the polygon |
|
94 |
#' * `quantile` - arbitrary quantile(s) of cell values, specified in |
|
95 |
#' `quantiles`, weighted by the fraction of each cell that is |
|
96 |
#' covered by the polygon |
|
97 |
#' * `mode` - the most common cell value, weighted by the fraction of |
|
98 |
#' each cell that is covered by the polygon. Where multiple |
|
99 |
#' values occupy the same maximum number of weighted cells, |
|
100 |
#' the largest value will be returned. |
|
101 |
#' * `majority` - synonym for `mode` |
|
102 |
#' * `minority` - the least common cell value, weighted by the fraction |
|
103 |
#' of each cell that is covered by the polygon. Where |
|
104 |
#' multiple values occupy the same minimum number of |
|
105 |
#' weighted cells, the smallest value will be returned. |
|
106 |
#' * `variety` - the number of distinct values in cells that are wholly or |
|
107 |
#' partially covered by the polygon. |
|
108 |
#' * `variance` - the population variance of cell values, weighted by the |
|
109 |
#' fraction of each cell that is covered by the polygon. |
|
110 |
#' * `stdev` - the population standard deviation of cell values, weighted by |
|
111 |
#' the fraction of each cell that is covered by the polygon. |
|
112 |
#' * `coefficient_of_variation` - the population coefficient of variation of |
|
113 |
#' cell values, weighted by the fraction of each |
|
114 |
#' cell that is covered by the polygon. |
|
115 |
#' * `weighted_mean` - the mean cell value, weighted by the product of |
|
116 |
#' the fraction of each cell covered by the polygon |
|
117 |
#' and the value of a second weighting raster provided |
|
118 |
#' as `weights` |
|
119 |
#' * `weighted_sum` - the sum of defined raster cell values, multiplied by |
|
120 |
#' the fraction of each cell that is covered by the polygon |
|
121 |
#' and the value of a second weighting raster provided |
|
122 |
#' as `weights` |
|
123 |
#' * `frac` - returns one column for each possible value of `x`, with the |
|
124 |
#' the fraction of defined raster cells that are equal to that |
|
125 |
#' value. |
|
126 |
#' * `weighted_frac` - returns one column for each possible value of `x`, |
|
127 |
#' with the fraction of defined cells that are equal |
|
128 |
#' to that value, weighted by `weights.` |
|
129 |
#' |
|
130 |
#' In all of the summary operations, `NA` values in the the primary raster (`x`) |
|
131 |
#' raster are ignored (i.e., `na.rm = TRUE`.) If `NA` values occur in the |
|
132 |
#' weighting raster, the result of the weighted operation will be `NA`. `NA` |
|
133 |
#' values in both `x` and `weights` can be replaced on-the-fly using the |
|
134 |
#' `default_value` and `default_weight` arguments. |
|
135 |
#' |
|
136 |
#' ## User-defined summary functions |
|
137 |
#' |
|
138 |
#' If no predefined summary operation is suitable, a user-defined R function may |
|
139 |
#' be provided as `fun`. The function will be called once for each feature and |
|
140 |
#' must return either a single value or a data frame. The results of the |
|
141 |
#' function for each feature will be combined and returned by `exact_extract`. |
|
142 |
#' |
|
143 |
#' The simplest way to write a summary function is to set |
|
144 |
#' argument `summarize_df = TRUE`. (For backwards compatibility, this is not the |
|
145 |
#' default.) In this mode, the summary function takes the signature |
|
146 |
#' `function(df, ...)` where `df` is the same data frame that would be returned |
|
147 |
#' by `exact_extract` with `fun = NULL`. |
|
148 |
#' |
|
149 |
#' With `summarize_df = FALSE`, the function must have the signature |
|
150 |
#' `function(values, coverage_fractions, ...)` when weights are not used, and |
|
151 |
#' `function(values, coverage_fractions, weights, ...)` when weights are used. |
|
152 |
#' If the value and weight rasters each have a single layer, the function arguments |
|
153 |
#' will be vectors; if either has multiple layers, the function arguments will |
|
154 |
#' be data frames, with column names taken from the names of the value/weight |
|
155 |
#' rasters. Values brought in through the `include_xy`, `include_area`, |
|
156 |
#' `include_cell`, and `include_cols` arguments will be added to the `values` |
|
157 |
#' data frame. For most applications, it is simpler to set `summarize_df = TRUE` |
|
158 |
#' and work with all inputs in a single data frame. |
|
159 |
#' |
|
160 |
#' @param x a `RasterLayer`, `RasterStack`, `RasterBrick`, or `SpatRaster` |
|
161 |
#' @param y a `sf`, `sfc`, `SpatialPolygonsDataFrame`, or `SpatialPolygons` |
|
162 |
#' object with polygonal geometries |
|
163 |
#' @param fun an optional function or character vector, as described below |
|
164 |
#' @param weights a weighting raster to be used with the `weighted_mean` |
|
165 |
#' and `weighted_sum` summary operations or a user-defined |
|
166 |
#' summary function. When `weights` is set to `'area'`, the |
|
167 |
#' cell areas of `x` will be calculated and used as weights. |
|
168 |
#' @param coverage_area if `TRUE`, output pixel `coverage_area` |
|
169 |
#' instead of `coverage_fraction` |
|
170 |
#' @param default_value an optional value to use instead of `NA` in `x` |
|
171 |
#' @param default_weight an optional value to use instead of `NA` in `weights` |
|
172 |
#' @param quantiles quantiles to be computed when `fun = 'quantile'` |
|
173 |
#' @param append_cols when `fun` is not `NULL`, an optional character vector |
|
174 |
#' of columns from `y` to be included in returned data frame. |
|
175 |
#' @param force_df always return a data frame instead of a vector, even if |
|
176 |
#' `x` has only one layer and `fun` has length 1 |
|
177 |
#' @param summarize_df pass values, coverage fraction/area, and weights to |
|
178 |
#' `fun` as a single data frame instead of |
|
179 |
#' separate arguments. |
|
180 |
#' @param full_colnames include the names of `x` and `weights` in |
|
181 |
#' the names of the data frame for each feature, even if |
|
182 |
#' `x` or `weights` has only one layer. |
|
183 |
#' This is useful when the results of multiple |
|
184 |
#' calls to `exact_extract` are combined with |
|
185 |
#' `cbind`. |
|
186 |
#' @param colname_fun an optional function used to construct column names. |
|
187 |
#' Should accept arguments `values` (name of value layer), |
|
188 |
#' `weights` (name of weight layer), `fun_name` (value of |
|
189 |
#' `fun`), `fun_value` (value associated with `fun`, for |
|
190 |
#' `fun %in% c('quantile', 'frac', 'weighted_frac)` |
|
191 |
#' `nvalues` (number of value layers), `weights` |
|
192 |
#' (number of weight layers) |
|
193 |
#' @param include_area if `TRUE`, and `fun` is `NULL`, augment |
|
194 |
#' the data frame for each feature with a column |
|
195 |
#' for the cell area. If the units of the raster CRS are |
|
196 |
#' degrees, the area in square meters will be calculated |
|
197 |
#' based on a spherical approximation of Earth. Otherwise, |
|
198 |
#' a Cartesian area will be calculated (and will be the |
|
199 |
#' same for all pixels.) If `TRUE` and `fun` is |
|
200 |
#' not `NULL`, add `area` to the data frame passed |
|
201 |
#' to `fun` for each feature. |
|
202 |
#' @param include_cell if `TRUE`, and `fun` is `NULL`, augment |
|
203 |
#' the data frame for each feature with a column |
|
204 |
#' for the cell index (`cell`). If `TRUE` and |
|
205 |
#' `fun` is not `NULL`, add `cell` to the |
|
206 |
#' data frame passed to `fun` for each feature. |
|
207 |
#' @param include_cols an optional character vector of column names in |
|
208 |
#' `y` to be added to the data frame for each |
|
209 |
#' feature that is either returned (when `fun` is |
|
210 |
#' `NULL`) or passed to `fun`. |
|
211 |
#' @param include_xy if `TRUE`, and `fun` is `NULL`, augment |
|
212 |
#' the returned data frame for each feature with columns |
|
213 |
#' for cell center coordinates (`x` and `y`). If |
|
214 |
#' `TRUE` and `fun` is not `NULL`, add |
|
215 |
#' `x` and `y` to the data frame passed to `fun` |
|
216 |
#' for each feature. |
|
217 |
#' @param stack_apply if `TRUE`, apply `fun` independently to |
|
218 |
#' each layer or `x` (and its corresponding layer |
|
219 |
#' of `weights`, if provided.) The number of |
|
220 |
#' layers in `x` and `weights` must equal |
|
221 |
#' each other or `1`, in which case the |
|
222 |
#' single layer raster will be recycled. |
|
223 |
#' If `FALSE`, apply `fun` to all layers of |
|
224 |
#' `x` (and `weights`) simultaneously. |
|
225 |
#' @param max_cells_in_memory the maximum number of raster cells to load at |
|
226 |
#' a given time when using a named summary operation |
|
227 |
#' for `fun` (as opposed to a function defined using |
|
228 |
#' R code). If a polygon covers more than `max_cells_in_memory` |
|
229 |
#' raster cells, it will be processed in multiple chunks. |
|
230 |
#' @param grid_compat_tol require value and weight grids to align within |
|
231 |
#' `grid_compat_tol` times the smaller of the two |
|
232 |
#' grid resolutions. |
|
233 |
#' @param progress if `TRUE`, display a progress bar during processing |
|
234 |
#' @param ... additional arguments to pass to `fun` |
|
235 |
#' @return a vector, data frame, or list of data frames, depending on the type |
|
236 |
#' of `x` and the value of `fun` (see Details) |
|
237 |
#' @examples |
|
238 |
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10) |
|
239 |
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))') |
|
240 |
#' |
|
241 |
#' # named summary operation on RasterLayer, returns vector |
|
242 |
#' exact_extract(rast, poly, 'mean') |
|
243 |
#' |
|
244 |
#' # two named summary operations on RasterLayer, returns data frame |
|
245 |
#' exact_extract(rast, poly, c('min', 'max')) |
|
246 |
#' |
|
247 |
#' # named summary operation on RasterStack, returns data frame |
|
248 |
#' stk <- raster::stack(list(a=rast, b=sqrt(rast))) |
|
249 |
#' exact_extract(stk, poly, 'mean') |
|
250 |
#' |
|
251 |
#' # named weighted summary operation, returns vector |
|
252 |
#' weights <- raster::raster(matrix(runif(100), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10) |
|
253 |
#' exact_extract(rast, poly, 'weighted_mean', weights=weights) |
|
254 |
#' |
|
255 |
#' # custom summary function, returns vector |
|
256 |
#' exact_extract(rast, poly, function(value, cov_frac) length(value[cov_frac > 0.9])) |
|
257 |
#' |
|
258 |
#' @name exact_extract |
|
259 |
#' @md |
|
260 |
NULL |
|
261 | ||
262 |
.exact_extract <- function(x, y, fun=NULL, ..., |
|
263 |
weights=NULL, |
|
264 |
append_cols=NULL, |
|
265 |
coverage_area=FALSE, |
|
266 |
default_value=NA_real_, |
|
267 |
default_weight=NA_real_, |
|
268 |
include_area=FALSE, |
|
269 |
include_cell=FALSE, |
|
270 |
include_cols=NULL, |
|
271 |
include_xy=FALSE, |
|
272 |
force_df=FALSE, |
|
273 |
full_colnames=FALSE, |
|
274 |
stack_apply=FALSE, |
|
275 |
summarize_df=FALSE, |
|
276 |
quantiles=NULL, |
|
277 |
progress=TRUE, |
|
278 |
max_cells_in_memory=30000000, |
|
279 |
grid_compat_tol=1e-3, |
|
280 |
colname_fun=NULL |
|
281 |
) { |
|
282 | 317x |
area_weights <- is.character(weights) && length(weights) == 1 && weights == 'area' |
283 | 317x |
if (area_weights) { |
284 | 9x |
weights <- NULL |
285 |
} |
|
286 | ||
287 | 317x |
.validateFlag(coverage_area, 'coverage_area') |
288 | 314x |
.validateFlag(force_df, 'force_df') |
289 | 311x |
.validateFlag(full_colnames, 'full_colnames') |
290 | 308x |
.validateFlag(include_area, 'include_area') |
291 | 305x |
.validateFlag(include_cell, 'include_cell') |
292 | 302x |
.validateFlag(include_xy, 'include_xy') |
293 | 299x |
.validateFlag(progress, 'progress') |
294 | 296x |
.validateFlag(stack_apply, 'stack_apply') |
295 | 293x |
.validateFlag(summarize_df, 'summarize_df') |
296 | 290x |
.validateNumericScalar(max_cells_in_memory, 'max_cells_in_memory') |
297 | 286x |
.validateNumericScalar(grid_compat_tol, 'grid_compat_tol') |
298 | 286x |
.validateNumericScalarOrNA(default_value, 'default_value') |
299 | 282x |
.validateNumericScalarOrNA(default_weight, 'default_weight') |
300 | ||
301 | 282x |
.validateUniqueNames(x) |
302 | 281x |
.validateUniqueNames(weights) |
303 | ||
304 | 281x |
if(!is.null(append_cols)) { |
305 | 10x |
if (!inherits(y, 'sf')) { |
306 | 2x |
stop(sprintf('append_cols only supported for sf arguments (received %s)', |
307 | 2x |
paste(class(y), collapse = ' '))) |
308 |
} |
|
309 | ||
310 | 8x |
force_df <- TRUE |
311 |
} |
|
312 | ||
313 | 279x |
if(!is.null(include_cols)) { |
314 | 7x |
if (!inherits(y, 'sf')) { |
315 | 1x |
stop(sprintf('include_cols only supported for sf arguments (received %s)', |
316 | 1x |
paste(class(y), collapse = ' '))) |
317 |
} |
|
318 |
} |
|
319 | ||
320 | 278x |
if(sf::st_geometry_type(y, by_geometry = FALSE) == 'GEOMETRY') { |
321 | 2x |
if (!all(sf::st_dimension(y) == 2)) { |
322 | 1x |
stop("Features in sfc_GEOMETRY must be polygonal") |
323 |
} |
|
324 |
} |
|
325 | ||
326 | 277x |
if(!is.null(weights)) { |
327 | 89x |
if (!(.isRaster(weights) || weights == 'area')) { |
328 | 1x |
stop("Weights must be a Raster object or \"area\".") |
329 |
} |
|
330 | ||
331 | 88x |
if (is.character(fun) && !any(startsWith(fun, "weighted"))) { |
332 | 9x |
warning("Weights provided but no requested operations use them.") |
333 |
} |
|
334 | ||
335 | 88x |
if (inherits(weights, 'BasicRaster') && !is.na(sf::st_crs(x))) { |
336 | 72x |
if (is.na(sf::st_crs(weights))) { |
337 | 1x |
warning("No CRS specified for weighting raster; assuming it has the same CRS as the value raster.") |
338 | 71x |
} else if (sf::st_crs(x) != sf::st_crs(weights)) { |
339 | 1x |
stop("Weighting raster does not have the same CRS as value raster.") |
340 |
} |
|
341 |
} |
|
342 |
} |
|
343 | ||
344 | 275x |
analysis_crs <- sf::st_crs(x) |
345 | 275x |
if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) { |
346 | 1x |
warning("No CRS specified for raster; assuming it has the same CRS as the polygons.") |
347 | 1x |
analysis_crs <- sf::st_crs(y) |
348 | 274x |
} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) { |
349 | 2x |
warning("No CRS specified for polygons; assuming they have the same CRS as the raster.") |
350 | 272x |
} else if(sf::st_crs(x) != sf::st_crs(y)) { |
351 | 1x |
y <- sf::st_transform(y, sf::st_crs(x)) |
352 | 1x |
warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")") |
353 |
} |
|
354 | 275x |
if(area_weights || include_area || coverage_area) { |
355 | 22x |
area_method <- .areaMethod(analysis_crs) |
356 |
} else { |
|
357 | 253x |
area_method <- NULL |
358 |
} |
|
359 | ||
360 | 275x |
if (is.character(fun)) { |
361 | 157x |
if (length(fun) == 0) { |
362 | 1x |
stop("No summary operations provided.") |
363 |
} |
|
364 | ||
365 | 156x |
if (length(list(...)) > 0) { |
366 | 1x |
stop("exact_extract was called with a named summary operation that ", |
367 | 1x |
"does not accept additional arguments ...") |
368 |
} |
|
369 | ||
370 | 155x |
if (include_xy) { |
371 | 1x |
stop("include_xy must be FALSE for named summary operations") |
372 |
} |
|
373 | ||
374 | 154x |
if (include_cell) { |
375 | 1x |
stop("include_cell must be FALSE for named summary operations") |
376 |
} |
|
377 | ||
378 | 153x |
if (include_area) { |
379 | 1x |
stop("include_area must be FALSE for named summary operations") |
380 |
} |
|
381 | ||
382 | 152x |
if (!is.null(include_cols)) { |
383 | 1x |
stop("include_cols not supported for named_summary operations (see argument append_cols)") |
384 |
} |
|
385 | ||
386 | 151x |
if (summarize_df) { |
387 | 1x |
stop("summarize_df can only be used when `fun` is an R function") |
388 |
} |
|
389 | 118x |
} else if (is.function(fun)) { |
390 | 63x |
if (summarize_df) { |
391 | 8x |
if (.num_expected_args(fun) < 1) { |
392 | 2x |
stop("exact_extract was called with a function that does not appear to ", |
393 | 2x |
"be of the form `function(df, ...`).") |
394 |
} |
|
395 | 55x |
} else if (is.null(weights)) { |
396 | 40x |
if (.num_expected_args(fun) < 2) { |
397 | 5x |
stop("exact_extract was called with a function that does not appear to ", |
398 | 5x |
"be of the form `function(values, coverage_fractions, ...)`. If ", |
399 | 5x |
"the summary function should accept a single data frame argument, ", |
400 | 5x |
"set `summarize_df = TRUE`.") |
401 |
} |
|
402 | 15x |
} else if (.num_expected_args(fun) < 3) { |
403 | 3x |
stop("exact_extract was called with a function that does not appear to ", |
404 | 3x |
"be of the form `function(values, coverage_fractions, weights, ...)`.", |
405 | 3x |
"If the summary function should accept a single data frame argument, ", |
406 | 3x |
"set `summarize_df = TRUE`.") |
407 |
} |
|
408 | 55x |
} else if (is.null(fun)) { |
409 | 53x |
if (length(list(...)) > 0) { |
410 | 1x |
stop("Unexpected arguments: ", paste(names(list(...)), collapse = ',')) |
411 |
} |
|
412 | ||
413 | 52x |
if (summarize_df) { |
414 | 1x |
stop("summarize_df can only be used when `fun` is an R function") |
415 |
} |
|
416 | 51x |
if (stack_apply) { |
417 | 1x |
stop("stack_apply can only be used when `fun` is a summary operation or function") |
418 |
} |
|
419 | 50x |
if (!is.null(append_cols)) { |
420 | 1x |
stop("append_cols can only be used when `fun` is a summary operation or function. See `include_cols`.") |
421 |
} |
|
422 |
} else { |
|
423 | 2x |
stop("fun must be a character vector, function, or NULL") |
424 |
} |
|
425 | ||
426 | 252x |
geoms <- sf::st_geometry(y) |
427 | ||
428 | 252x |
if (requireNamespace('terra', quietly = TRUE)) { |
429 | 252x |
strategies <- c('check_block_size', 'eager_load') |
430 |
} else { |
|
431 | ! |
strategies = c() |
432 |
} |
|
433 | ||
434 | 252x |
x_orig <- NULL |
435 | 252x |
if ('eager_load' %in% strategies && length(geoms) > 1 && (!.isInMemory(x))) { |
436 |
# Eagerly load the entire area to be processed into memory. If the raster |
|
437 |
# block sizes are large, this potentially allows us to cache data in memory |
|
438 |
# that would not fit within the GDAL block cache (because we can cache an |
|
439 |
# area smaller than a single block). We only do this when the terra package |
|
440 |
# is available, because raster::crop throws an error when we set snap = |
|
441 |
# 'out'. |
|
442 | 8x |
if (include_cell) { |
443 |
# retain a reference to the original so that we can use it for include_cell |
|
444 | 2x |
x_orig <- x |
445 |
} |
|
446 | 8x |
x <- .eagerLoad(x, geoms, max_cells_in_memory, message_on_fail = progress) |
447 | 8x |
weights <- .eagerLoad(weights, geoms, max_cells_in_memory, message_on_fail = progress) |
448 |
} |
|
449 | ||
450 |
# At this point, if the data have not been preloaded into memory, either: |
|
451 |
# - we don't have the terra package available |
|
452 |
# - we can't fit the whole processing area in memory |
|
453 |
# Whether this is a problem depends on how the raster is chunked. If the |
|
454 |
# raster is poorly chunked, we can't do anything to improve the perormance, |
|
455 |
# but we try to detect the situation so we can alert the user. |
|
456 | 252x |
if (!.isInMemory(x)) { |
457 |
# - Do we have a RasterStack? Suggest using terra instead. |
|
458 | 5x |
if (inherits(x, 'RasterStack')) { |
459 | 1x |
message('exact_extract may perform poorly when extracting from a RasterStack. ', |
460 | 1x |
'It is recommended to use a SpatRaster from the terra package instead. ', |
461 | 1x |
'Convert using terra::rast(x)') |
462 |
} |
|
463 | ||
464 | 5x |
if (('check_block_size' %in% strategies) && inherits(x, 'SpatRaster')) { |
465 |
# Can we fit at least one block from every layer in GDAL's block cache? |
|
466 |
# If not, warn the user of expected poor performance. |
|
467 | 4x |
block_dim <- .blockSize(x) |
468 | 4x |
bytes_per_pixel <- 8L # terra has no datatype function? |
469 | 4x |
block_sz_mb <- prod(block_dim[1:2])* bytes_per_pixel / 1024 / 1024 |
470 | 4x |
cache_sz_mb <- terra::gdalCache() |
471 | ||
472 | 4x |
min_cache_needed_mb <- ceiling(block_sz_mb) * .numLayers(x) |
473 | ||
474 | ||
475 | 4x |
if (min_cache_needed_mb >= cache_sz_mb) { |
476 | 1x |
message('Loading one raster block of data per layer requires approximately ', |
477 | 1x |
min_cache_needed_mb, ' MB but the GDAL block size cache is only ', |
478 | 1x |
cache_sz_mb, ' MB. This is likely to result in very slow performance. ', |
479 | 1x |
'It is recommended to increase the cache size using ', |
480 | 1x |
'terra::gdalCache(new_size), reduce the number of layers in the ', |
481 | 1x |
'input raster, or rewrite the input raster using a smaller chunk size.') |
482 |
} |
|
483 |
} |
|
484 |
} |
|
485 | ||
486 | 252x |
update_progress <- .createProgress(progress, length(geoms)) |
487 | 252x |
tryCatch({ |
488 | 252x |
x <- .startReading(x) |
489 | 252x |
weights <- .startReading(weights) |
490 | ||
491 | 252x |
if (is.character(fun)) { |
492 |
# Compute all stats in C++. |
|
493 |
# CPP_stats returns a matrix, which gets turned into a column by sapply |
|
494 |
# Results has one column per feature and one row per stat/raster layer |
|
495 | 150x |
if((is.null(weights) && !area_weights) && any(.isWeighted(fun))) { |
496 | 2x |
stop("Weighted stat requested but no weights provided.") |
497 |
} |
|
498 | ||
499 | 148x |
results <- lapply(sf::st_as_binary(geoms, EWKB=TRUE), function(wkb) { |
500 | 214x |
ret <- CPP_stats(x, weights, wkb, default_value, default_weight, coverage_area, area_method, fun, max_cells_in_memory, grid_compat_tol, quantiles) |
501 | 204x |
update_progress() |
502 | 204x |
return(ret) |
503 |
}) |
|
504 | ||
505 | 138x |
if ('frac' %in% fun || 'weighted_frac' %in% fun) { |
506 | 5x |
unique_values <- unique(Reduce(c, lapply(results, function(r) attr(r, 'unique_values')))) |
507 |
} else { |
|
508 | 133x |
unique_values <- numeric() |
509 |
} |
|
510 | ||
511 | 138x |
result_cols <- .resultColumns(names(x), names(weights), fun, full_colnames, quantiles, unique_values, colname_fun = colname_fun) |
512 | 138x |
consistent_cols <- !(result_cols$stat_name %in% c('frac', 'weighted_frac')) |
513 | ||
514 |
# Construct an empty matrix that we can populate with stats returned |
|
515 |
# for each geometry, with 0 as a default where a geometry does not |
|
516 |
# return a stat (e.g., 'frac' for a specific value). Then convert the |
|
517 |
# filled matrix to a data frame. (This is faster than populating the |
|
518 |
# data frame directly, because the data frame assignment operator is |
|
519 |
# slow.) |
|
520 | 138x |
result_mat <- matrix(0.0, nrow = length(geoms), ncol = nrow(result_cols)) |
521 | 138x |
for (i in seq_along(results)) { |
522 | 204x |
cols <- consistent_cols | result_cols$base_value %in% attr(results[[i]], 'unique_values') |
523 | 204x |
result_mat[i, cols] <- results[[i]] |
524 |
} |
|
525 | 138x |
colnames(result_mat) <- result_cols$colname |
526 | ||
527 | 138x |
result_df <- as.data.frame(result_mat) |
528 | ||
529 | 138x |
if (ncol(result_df) == 1 && !force_df) { |
530 | 111x |
return(result_df[[1]]) |
531 |
} |
|
532 | ||
533 |
# drop duplicated columns (occurs when an unweighted stat is |
|
534 |
# requested alongside a weighted stat, with a stack of weights) |
|
535 | 27x |
result_df <- result_df[, unique(names(result_df)), drop = FALSE] |
536 | ||
537 | 27x |
if (!is.null(append_cols)) { |
538 | 2x |
result_df <- cbind(sf::st_drop_geometry(y[, append_cols]), result_df) |
539 |
} |
|
540 | ||
541 | 26x |
return(result_df) |
542 |
} else { |
|
543 | 102x |
num_values <- .numLayers(x) |
544 | 102x |
value_names <- names(x) |
545 | 102x |
if (.isRaster(weights)) { |
546 | 29x |
num_weights <- .numLayers(weights) |
547 | 29x |
weight_names <- names(weights) |
548 | 73x |
} else if (area_weights) { |
549 | 1x |
num_weights <- 1 |
550 | 1x |
weight_names <- 'area' |
551 | 1x |
weights <- NULL |
552 |
} else { |
|
553 | 72x |
num_weights <- 0 |
554 | 72x |
weight_names <- NULL |
555 |
} |
|
556 | ||
557 | 102x |
if (stack_apply || (num_values == 1 && num_weights <= 1)) { |
558 | 85x |
apply_layerwise <- TRUE |
559 | ||
560 | 85x |
if (num_values > 1 && num_weights > 1 && num_values != num_weights) { |
561 | 1x |
stop(sprintf("Can't apply function layerwise with stacks of %d value layers and %d layers", num_values, num_weights)) |
562 |
} |
|
563 | ||
564 | 84x |
result_names <- .resultColNames(value_names, weight_names, fun, full_colnames, colname_fun = colname_fun) |
565 | 84x |
num_results <- max(num_weights, num_values) |
566 | 84x |
ind <- .valueWeightIndexes(num_values, num_weights) |
567 |
} else { |
|
568 | 17x |
apply_layerwise <- FALSE |
569 |
} |
|
570 | ||
571 |
# For R summary functions and data frame output, we avoid using |
|
572 |
# input layer names if there is only one value/weight layer |
|
573 | 101x |
if (length(value_names) == 1) { |
574 | 80x |
value_names <- 'value' |
575 |
} |
|
576 | 101x |
if (length(weight_names) == 1) { |
577 | 19x |
weight_names <- 'weight' |
578 |
} |
|
579 | ||
580 | 101x |
ret <- lapply(seq_along(geoms), function(feature_num) { |
581 | 154x |
wkb <- sf::st_as_binary(geoms[[feature_num]], EWKB=TRUE) |
582 | ||
583 | 154x |
if (!is.null(include_cols)) { |
584 | 6x |
include_col_values <- sf::st_drop_geometry(y[feature_num, include_cols]) |
585 |
} else { |
|
586 | 148x |
include_col_values <- NULL |
587 |
} |
|
588 | ||
589 |
# only raise a disaggregation warning for the first feature |
|
590 | 153x |
warn_on_disaggregate <- feature_num == 1 |
591 | ||
592 | 153x |
col_list <- CPP_exact_extract(x, |
593 | 153x |
x_orig, |
594 | 153x |
weights, |
595 | 153x |
wkb, |
596 | 153x |
default_value, |
597 | 153x |
default_weight, |
598 | 153x |
include_xy, |
599 | 153x |
include_cell, |
600 | 153x |
include_area, |
601 | 153x |
area_weights, |
602 | 153x |
coverage_area, |
603 | 153x |
area_method, |
604 | 153x |
include_col_values, |
605 | 153x |
value_names, |
606 | 153x |
weight_names, |
607 | 153x |
warn_on_disaggregate, |
608 | 153x |
grid_compat_tol) |
609 | 153x |
if (coverage_area) { |
610 | 4x |
coverage_col <- 'coverage_area' |
611 |
} else { |
|
612 | 149x |
coverage_col <- 'coverage_fraction' |
613 |
} |
|
614 | ||
615 | 153x |
if (!is.null(include_cols)) { |
616 |
# Replicate the include_cols vectors to be as long as the other columns, |
|
617 |
# so we can use quickDf |
|
618 | 5x |
nrow <- length(col_list[[coverage_col]]) |
619 | 5x |
col_list[include_cols] <- lapply(col_list[include_cols], rep, nrow) |
620 |
} |
|
621 | 153x |
df <- .quickDf(col_list) |
622 | ||
623 | 153x |
update_progress() |
624 | ||
625 |
# No summary function? Return the whole data frame. |
|
626 | 153x |
if (is.null(fun)) { |
627 | 64x |
return(df) |
628 |
} |
|
629 | ||
630 |
# Summary function accepts a data frame? Pass it along. |
|
631 | 89x |
if (summarize_df && !apply_layerwise) { |
632 | 3x |
return(fun(df, ...)) |
633 |
} |
|
634 | ||
635 |
# Break the data frame into components that we can pass |
|
636 |
# separately to the summary functions. |
|
637 | 86x |
included_cols_df <- df[, !(names(df) %in% c(value_names, weight_names, coverage_col)), drop = FALSE] |
638 | 86x |
vals_df <- df[, value_names, drop = FALSE] |
639 | 86x |
weights_df <- df[, weight_names, drop = FALSE] |
640 | 86x |
cov_fracs <- df[[coverage_col]] |
641 | ||
642 | 86x |
if (apply_layerwise) { |
643 | 79x |
result <- lapply(seq_len(num_results), function(i) { |
644 | 92x |
if (summarize_df) { |
645 |
# Pack everything into a single data frame and pass it to `fun` |
|
646 | 5x |
arg_df <- cbind(value = vals_df[[ind$values[i]]], |
647 | 5x |
included_cols_df) |
648 | 5x |
if (num_weights > 0) { |
649 | 3x |
arg_df$weight <- weights_df[[ind$weights[i]]] |
650 |
} |
|
651 | 5x |
arg_df[[coverage_col]] <- df[[coverage_col]] |
652 | ||
653 | 5x |
return(fun(arg_df, ...)) |
654 |
} else { |
|
655 |
# Pull values and weights out into vectors, unless we have |
|
656 |
# included columns (x/y/cell, etc.) in which case values |
|
657 |
# remain a data frame. Retained for backwards compat. |
|
658 | ||
659 | 87x |
vx <- vals_df[, ind$values[i]] |
660 | 87x |
if (ncol(included_cols_df) > 0) { |
661 | 7x |
vx <- cbind(data.frame(value = vx), included_cols_df) |
662 |
} |
|
663 | 87x |
if (num_weights == 0) { |
664 | 73x |
return(fun(vx, cov_fracs, ...)) |
665 |
} else { |
|
666 | 14x |
return(fun(vx, cov_fracs, weights_df[, ind$weights[i]], ...)) |
667 |
} |
|
668 |
} |
|
669 |
}) |
|
670 | ||
671 | 79x |
if (num_results == 1) { |
672 | 70x |
return(result[[1]]) |
673 |
} |
|
674 | ||
675 | 9x |
names(result) <- result_names |
676 | ||
677 | 9x |
return(.quickDf(result)) |
678 |
} else { |
|
679 |
# Pass all layers to callback, to be handled together |
|
680 |
# Included columns (x/y/cell) are passed with the values. |
|
681 |
# Pass single-column data frames as vectors. |
|
682 | 7x |
vals_df <- .singleColumnToVector(cbind(vals_df, included_cols_df)) |
683 | 7x |
weights_df <- .singleColumnToVector(weights_df) |
684 | ||
685 | 7x |
if (num_weights == 0) { |
686 | 3x |
return(fun(vals_df, cov_fracs, ...)) |
687 |
} else { |
|
688 | 4x |
return(fun(vals_df, cov_fracs, weights_df, ...)) |
689 |
} |
|
690 |
} |
|
691 |
}) |
|
692 | ||
693 |
# if we have columns to append, iterate over the results and |
|
694 |
# append the columns to each, coercing into a data frame if |
|
695 |
# necessary. We need to iterate over the results before binding |
|
696 |
# them together because we don't know the legnth of each result |
|
697 |
# or even if their lengths are constant. |
|
698 | 99x |
if (!is.null(append_cols)) { |
699 | 5x |
for (i in seq_along(ret)) { |
700 | 9x |
if (!is.data.frame(ret[[i]])) { |
701 | 5x |
ret[[i]] = data.frame(result = ret[[i]]) |
702 |
} |
|
703 | ||
704 | 9x |
if (nrow(ret[[i]]) >= 1) { |
705 | 7x |
append_row <- i |
706 |
} else { |
|
707 |
# cbinding row 0 cleanly brings in zero-length columns |
|
708 |
# of the correct names and types, in the correct positions |
|
709 | 2x |
append_row <- 0 |
710 |
} |
|
711 | ||
712 | 9x |
ret[[i]] <- cbind(sf::st_drop_geometry(y[append_row, append_cols]), |
713 | 9x |
ret[[i]], |
714 | 9x |
row.names = NULL) |
715 |
} |
|
716 |
} |
|
717 | ||
718 | 98x |
if (!is.null(fun)) { |
719 | 50x |
if (all(sapply(ret, is.data.frame))) { |
720 |
# function returned a data frame for each polygon? rbind them |
|
721 | 14x |
if (requireNamespace('dplyr', quietly = TRUE)) { |
722 | 14x |
ret <- dplyr::bind_rows(ret) # handle column name mismatches |
723 |
} else { |
|
724 | ! |
ret <- do.call(rbind, ret) |
725 |
} |
|
726 |
} else { |
|
727 |
# function returned something else; combine the somethings into |
|
728 |
# an array |
|
729 | 36x |
ret <- simplify2array(ret) |
730 | ||
731 | 36x |
if (force_df) { |
732 | 2x |
ret <- data.frame(result = ret) |
733 |
} |
|
734 |
} |
|
735 |
} |
|
736 | ||
737 | 98x |
return(ret) |
738 |
} |
|
739 | 252x |
}, finally={ |
740 | 252x |
.stopReading(x) |
741 | 252x |
.stopReading(weights) |
742 |
}) |
|
743 |
} |
|
744 | ||
745 |
# faster replacement for as.data.frame when input is a named list |
|
746 |
# with equal-length columns |
|
747 |
# from Advanced R, sec. 24.4.2 |
|
748 |
.quickDf <- function(lst) { |
|
749 |
class(lst) <- 'data.frame' |
|
750 |
attr(lst, 'row.names') <- .set_row_names(length(lst[[1]])) |
|
751 |
lst |
|
752 |
} |
|
753 | ||
754 |
.singleColumnToVector <- function(df) { |
|
755 |
if (ncol(df) == 1) { |
|
756 |
df[, 1] |
|
757 |
} else { |
|
758 |
df |
|
759 |
} |
|
760 |
} |
|
761 | ||
762 |
#' @import sf |
|
763 |
#' @import raster |
|
764 |
#' @useDynLib exactextractr |
|
765 |
#' @rdname exact_extract |
|
766 |
#' @export |
|
767 |
setMethod('exact_extract', signature(x='Raster', y='sf'), |
|
768 |
.exact_extract) |
|
769 | ||
770 |
#' @useDynLib exactextractr |
|
771 |
#' @rdname exact_extract |
|
772 |
#' @export |
|
773 |
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygonsDataFrame'), |
|
774 | 1x |
function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...)) |
775 | ||
776 |
#' @useDynLib exactextractr |
|
777 |
#' @rdname exact_extract |
|
778 |
#' @export |
|
779 |
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygons'), |
|
780 | 1x |
function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...)) |
781 | ||
782 |
#' @useDynLib exactextractr |
|
783 |
#' @rdname exact_extract |
|
784 |
#' @export |
|
785 |
setMethod('exact_extract', signature(x='Raster', y='sfc_MULTIPOLYGON'), .exact_extract) |
|
786 | ||
787 |
#' @useDynLib exactextractr |
|
788 |
#' @rdname exact_extract |
|
789 |
#' @export |
|
790 |
setMethod('exact_extract', signature(x='Raster', y='sfc_POLYGON'), .exact_extract) |
|
791 | ||
792 |
#' @useDynLib exactextractr |
|
793 |
#' @rdname exact_extract |
|
794 |
#' @export |
|
795 |
setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRY'), .exact_extract) |
|
796 | ||
797 |
# CRAN version of sf does not explicitly declare this class. If we do not do it |
|
798 |
# ourselves, documentation is not generated for the `sfc_GEOMETRYCOLLECTION` |
|
799 |
# overload. |
|
800 |
setOldClass('sfc_GEOMETRYCOLLECTION') |
|
801 | ||
802 |
#' @useDynLib exactextractr |
|
803 |
#' @rdname exact_extract |
|
804 |
#' @export |
|
805 |
setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRYCOLLECTION'), .exact_extract) |
|
806 | ||
807 |
#' @import sf |
|
808 |
#' @useDynLib exactextractr |
|
809 |
#' @rdname exact_extract |
|
810 |
#' @export |
|
811 |
setMethod('exact_extract', signature(x='SpatRaster', y='sf'), |
|
812 |
.exact_extract) |
|
813 | ||
814 |
#' @useDynLib exactextractr |
|
815 |
#' @rdname exact_extract |
|
816 |
#' @export |
|
817 |
setMethod('exact_extract', signature(x='SpatRaster', y='SpatialPolygonsDataFrame'), |
|
818 | ! |
function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...)) |
819 | ||
820 |
#' @useDynLib exactextractr |
|
821 |
#' @rdname exact_extract |
|
822 |
#' @export |
|
823 |
setMethod('exact_extract', signature(x='SpatRaster', y='SpatialPolygons'), |
|
824 | ! |
function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...)) |
825 | ||
826 |
#' @useDynLib exactextractr |
|
827 |
#' @rdname exact_extract |
|
828 |
#' @export |
|
829 |
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_MULTIPOLYGON'), .exact_extract) |
|
830 | ||
831 |
#' @useDynLib exactextractr |
|
832 |
#' @rdname exact_extract |
|
833 |
#' @export |
|
834 |
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_POLYGON'), .exact_extract) |
|
835 | ||
836 |
#' @useDynLib exactextractr |
|
837 |
#' @rdname exact_extract |
|
838 |
#' @export |
|
839 |
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_GEOMETRY'), .exact_extract) |
|
840 | ||
841 |
#' @useDynLib exactextractr |
|
842 |
#' @rdname exact_extract |
|
843 |
#' @export |
|
844 |
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_GEOMETRYCOLLECTION'), .exact_extract) |
1 |
# Copyright (c) 2020-2022 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 | 23x |
setGeneric("exact_resample", function(x, y, ...) |
15 | 23x |
standardGeneric("exact_resample")) |
16 | ||
17 |
#' Resample a raster to a new grid |
|
18 |
#' |
|
19 |
#' @param x a \code{RasterLayer} or \code{SpatRaster} to be resampled |
|
20 |
#' @param y a raster of the same class as \code{x} with a grid definition to |
|
21 |
#' which \code{x} should be resampled |
|
22 |
#' @param fun a named summary operation or R function to be used for the resampling |
|
23 |
#' @param coverage_area use cell coverage areas instead of coverage fractions |
|
24 |
#' in \code{fun} |
|
25 |
#' @return a resampled version of \code{x}, returned as a \code{RasterLayer} or |
|
26 |
#' \code{SpatRaster}, depending on the values of \code{x} and \code{y} |
|
27 |
#' |
|
28 |
#' @name exact_resample |
|
29 |
NULL |
|
30 | ||
31 |
.exact_resample <- function(x, y, fun, coverage_area = FALSE) { |
|
32 | 23x |
analysis_crs <- sf::st_crs(x) |
33 | 23x |
if (sf::st_crs(x) != sf::st_crs(y)) { |
34 | 3x |
if (is.na(sf::st_crs(x))) { |
35 | 1x |
warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.") |
36 | 1x |
analysis_crs <- sf::st_crs(y) |
37 | 2x |
} else if (is.na(sf::st_crs(y))) { |
38 | 1x |
warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.") |
39 |
} else { |
|
40 | 1x |
stop('Destination raster must have same CRS as source.') |
41 |
} |
|
42 |
} |
|
43 | ||
44 | 22x |
if (is.character(fun)) { |
45 | 15x |
if (length(fun) != 1) { |
46 | 2x |
stop("Only a single operation may be used for resampling.") |
47 |
} |
|
48 | ||
49 | 13x |
if (startsWith(fun, 'weighted')) { |
50 | 1x |
stop("Weighted operations cannot be used for resampling.") |
51 |
} |
|
52 | ||
53 | 12x |
if (.numLayers(x) != 1) { |
54 | 1x |
stop("Raster to be resampled must have a single layer.") |
55 |
} |
|
56 | ||
57 | 11x |
summary_stat <- fun |
58 | 11x |
summary_fun <- NULL |
59 |
} else { |
|
60 | 7x |
if (!is.function(fun)) { |
61 | ! |
stop("fun must be a named summary operation or an R function") |
62 |
} |
|
63 | ||
64 | 7x |
if (.num_expected_args(fun) != 2) { |
65 | 1x |
stop("exact_extract was called with a function that does not appear to ", |
66 | 1x |
"be of the form `function(values, coverage_fractions)`.") |
67 |
} |
|
68 | ||
69 | 6x |
summary_stat <- NULL |
70 | 6x |
summary_fun <- fun |
71 |
} |
|
72 | ||
73 | 17x |
.validateFlag(coverage_area, 'coverage_area') |
74 | ||
75 | 17x |
area_method <- .areaMethod(analysis_crs) |
76 | ||
77 | 17x |
x <- .startReading(x) |
78 | 17x |
tryCatch({ |
79 | 17x |
ret <- CPP_resample(x, |
80 | 17x |
y, |
81 | 17x |
summary_stat, |
82 | 17x |
summary_fun, |
83 | 17x |
coverage_area, |
84 | 17x |
area_method) |
85 | 13x |
if (inherits(x, 'SpatRaster')) { |
86 | 7x |
terra::rast(ret) |
87 |
} else { |
|
88 | 6x |
ret |
89 |
} |
|
90 | 17x |
}, finally={ |
91 | 17x |
.stopReading(x) |
92 |
}) |
|
93 |
} |
|
94 | ||
95 |
#' @import sf |
|
96 |
#' @import raster |
|
97 |
#' @useDynLib exactextractr |
|
98 |
#' @rdname exact_resample |
|
99 |
#' @export |
|
100 |
setMethod('exact_resample', signature(x='RasterLayer', y='RasterLayer'), .exact_resample) |
|
101 | ||
102 |
#' @useDynLib exactextractr |
|
103 |
#' @rdname exact_resample |
|
104 |
#' @export |
|
105 |
setMethod('exact_resample', signature(x='SpatRaster', y='SpatRaster'), .exact_resample) |
1 |
# Copyright (c) 2018-2021 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 | 13x |
setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...) |
15 | 13x |
standardGeneric('coverage_fraction')) |
16 | ||
17 |
.coverage_fraction <- function(x, y, crop) { |
|
18 | 13x |
if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) { |
19 | 1x |
warning("No CRS specified for raster; assuming it has the same CRS as the polygons.") |
20 | 12x |
} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) { |
21 | 1x |
warning("No CRS specified for polygons; assuming they have the same CRS as the raster.") |
22 | 11x |
} else if(sf::st_crs(x) != sf::st_crs(y)) { |
23 | 1x |
y <- sf::st_transform(y, sf::st_crs(x)) |
24 | 1x |
warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")") |
25 |
} |
|
26 | ||
27 | 13x |
lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) { |
28 | 13x |
CPP_coverage_fraction(x, wkb, crop) |
29 |
}) |
|
30 |
} |
|
31 | ||
32 |
#' Compute the fraction of raster cells covered by a polygon |
|
33 |
#' |
|
34 |
#' @param x a (possibly empty) \code{RasterLayer} whose resolution and |
|
35 |
#' extent will be used for the generated \code{RasterLayer}. |
|
36 |
#' @param y a \code{sf} object with polygonal geometries |
|
37 |
#' @param crop if \code{TRUE}, each generated \code{RasterLayer} will be |
|
38 |
#' cropped to the extent of its associated feature. |
|
39 |
#' @return a list with a \code{RasterLayer} for each feature in \code{y}. |
|
40 |
#' Values of the raster represent the fraction of each |
|
41 |
#' cell in \code{x} that is covered by \code{y}. |
|
42 |
#' @examples |
|
43 |
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10) |
|
44 |
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))') |
|
45 |
#' |
|
46 |
#' cov_frac <- coverage_fraction(rast, poly)[[1]] |
|
47 |
#' @name coverage_fraction |
|
48 |
NULL |
|
49 | ||
50 |
#' @import sf |
|
51 |
#' @import raster |
|
52 |
#' @useDynLib exactextractr |
|
53 |
#' @rdname coverage_fraction |
|
54 |
#' @export |
|
55 |
setMethod('coverage_fraction', signature(x='RasterLayer', y='sf'), function(x, y, crop=FALSE) { |
|
56 | ! |
coverage_fraction(x, sf::st_geometry(y), crop) |
57 |
}) |
|
58 | ||
59 |
#' @rdname coverage_fraction |
|
60 |
#' @export |
|
61 |
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_MULTIPOLYGON'), .coverage_fraction) |
|
62 | ||
63 |
#' @rdname coverage_fraction |
|
64 |
#' @export |
|
65 |
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_POLYGON'), .coverage_fraction) |
|
66 | ||
67 |
#' @rdname coverage_fraction |
|
68 |
#' @export |
|
69 |
setMethod('coverage_fraction', signature(x='SpatRaster', y='sf'), function(x, y, crop=FALSE) { |
|
70 | ! |
coverage_fraction(x, sf::st_geometry(y), crop) |
71 |
}) |
|
72 | ||
73 |
#' @rdname coverage_fraction |
|
74 |
#' @export |
|
75 |
setMethod('coverage_fraction', signature(x='SpatRaster', y='sfc_MULTIPOLYGON'), .coverage_fraction) |
|
76 | ||
77 |
#' @rdname coverage_fraction |
|
78 |
#' @export |
|
79 |
setMethod('coverage_fraction', signature(x='SpatRaster', y='sfc_POLYGON'), .coverage_fraction) |
1 |
// Copyright (c) 2018-2021 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 | 13x |
Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, |
28 |
const Rcpp::RawVector & wkb, |
|
29 |
bool crop) |
|
30 |
{ |
|
31 |
try { |
|
32 | 26x |
GEOSAutoHandle geos; |
33 | 39x |
Rcpp::Environment xx = Rcpp::Environment::namespace_env("exactextractr"); |
34 | ||
35 | 13x |
auto grid = make_grid(rast); |
36 | 26x |
auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get()); |
37 | ||
38 | 13x |
if (crop) { |
39 | 2x |
grid = coverage_fraction.grid(); |
40 |
} |
|
41 | 26x |
RasterView<float> coverage_view(coverage_fraction, grid); |
42 | ||
43 | 13x |
Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()), |
44 | 26x |
static_cast<int>(grid.cols())}; |
45 | ||
46 | 904x |
for (size_t i = 0; i < grid.rows(); i++) { |
47 | 238204x |
for (size_t j = 0; j < grid.cols(); j++) { |
48 | 237313x |
weights(i, j) = coverage_view(i, j); |
49 | 237313x |
if (!crop && std::isnan(weights(i, j))) { |
50 | 209304x |
weights(i, j) = 0; |
51 |
} |
|
52 |
} |
|
53 |
} |
|
54 | ||
55 | 13x |
if (rast.inherits("SpatRaster")) { |
56 | 3x |
Rcpp::Environment terra = Rcpp::Environment::namespace_env("terra"); |
57 | 3x |
Rcpp::Function rastFn = terra["rast"]; |
58 | 3x |
Rcpp::Function extFn = terra["ext"]; |
59 | 3x |
Rcpp::Function crsFn = terra["crs"]; |
60 | ||
61 | 1x |
Rcpp::S4 ext = extFn(grid.xmin(), grid.xmax(), grid.ymin(), grid.ymax()); |
62 | ||
63 |
return rastFn(weights, |
|
64 | 2x |
Rcpp::Named("ext") = ext, |
65 | 3x |
Rcpp::Named("crs") = crsFn(rast)); |
66 |
} else { |
|
67 | 36x |
Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster"); |
68 | 36x |
Rcpp::Function rasterFn = raster["raster"]; |
69 | 24x |
Rcpp::Function crsFn = xx[".crs"]; |
70 | ||
71 |
return rasterFn(weights, |
|
72 | 24x |
Rcpp::Named("xmn")=grid.xmin(), |
73 | 24x |
Rcpp::Named("xmx")=grid.xmax(), |
74 | 24x |
Rcpp::Named("ymn")=grid.ymin(), |
75 | 24x |
Rcpp::Named("ymx")=grid.ymax(), |
76 | 36x |
Rcpp::Named("crs")=crsFn(rast)); |
77 |
} |
|
78 |
} catch (std::exception & e) { |
|
79 |
// throw predictable exception class |
|
80 |
#ifdef __SUNPRO_CC |
|
81 |
// Rcpp::stop crashes CRAN Solaris build |
|
82 |
// https://github.com/RcppCore/Rcpp/issues/1159 |
|
83 |
Rf_error(e.what()); |
|
84 |
return R_NilValue; |
|
85 |
#else |
|
86 |
Rcpp::stop(e.what()); |
|
87 |
#endif |
|
88 |
} |
|
89 |
} |
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 | 391x |
static inline bool is_integral(double d, double tol) { |
33 | 391x |
return std::abs(d - std::round(d)) <= tol; |
34 |
} |
|
35 | ||
36 |
template<typename extent_tag> |
|
37 |
class Grid { |
|
38 | ||
39 |
public: |
|
40 | 1836196x |
Grid(const Box & extent, double dx, double dy) : |
41 |
m_extent{extent}, |
|
42 |
m_dx{dx}, |
|
43 |
m_dy{dy}, |
|
44 | 1836196x |
m_num_rows{2*extent_tag::padding + (extent.ymax > extent.ymin ? static_cast<size_t>(std::round((extent.ymax - extent.ymin) / dy)) : 0)}, |
45 | 3672392x |
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 | 118372x |
static Grid make_empty() { |
49 | 118372x |
return Grid({0, 0, 0, 0}, 0, 0); |
50 |
} |
|
51 | ||
52 | 4794273x |
size_t get_column(double x) const { |
53 |
if (extent_tag::padding) { |
|
54 | 6384784x |
if (x < m_extent.xmin) |
55 | 878x |
return 0; |
56 | 6383906x |
if (x > m_extent.xmax) |
57 | 1180x |
return m_num_cols - 1; |
58 | 6382726x |
if (x == m_extent.xmax) // special case, returning the cell for which xmax is the right |
59 | 3192320x |
return m_num_cols - 2; |
60 |
} else { |
|
61 | 1601881x |
if (x < m_extent.xmin || x > m_extent.xmax) |
62 | ! |
throw std::out_of_range("x"); |
63 | ||
64 | 1601881x |
if (x == m_extent.xmax) |
65 | 801682x |
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 | 7186206x |
return std::min( |
72 | 2395402x |
extent_tag::padding + static_cast<size_t>(std::floor((x - m_extent.xmin) / m_dx)), |
73 | 4790804x |
get_column(m_extent.xmax)); |
74 |
} |
|
75 | ||
76 | 4793002x |
size_t get_row(double y) const { |
77 |
if (extent_tag::padding) { |
|
78 | 6382114x |
if (y > m_extent.ymax) |
79 | 2246x |
return 0; |
80 | 6379868x |
if (y < m_extent.ymin) |
81 | 1896x |
return m_num_rows - 1; |
82 | 6377972x |
if (y == m_extent.ymin) // special case, returning the cell for which ymin is the bottom |
83 | 3190236x |
return m_num_rows - 2; |
84 |
} else { |
|
85 | 1601945x |
if (y < m_extent.ymin || y > m_extent.ymax) |
86 | ! |
throw std::out_of_range("y"); |
87 | ||
88 | 1601945x |
if (y == m_extent.ymin) |
89 | 801682x |
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 | 7182393x |
return std::min( |
96 | 2394131x |
extent_tag::padding + static_cast<size_t>(std::floor((m_extent.ymax - y) / m_dy)), |
97 | 4788262x |
get_row(m_extent.ymin)); |
98 |
} |
|
99 | ||
100 | 29924212x |
bool empty() const { return m_num_rows <= 2*extent_tag::padding && m_num_cols <= 2*extent_tag::padding; } |
101 | ||
102 | 38787878x |
size_t rows() const { return m_num_rows; } |
103 | ||
104 | 115537109x |
size_t cols() const { return m_num_cols; } |
105 | ||
106 | 199x |
size_t size() const { return rows()*cols(); } |
107 | ||
108 | 5731163x |
double xmin() const { return m_extent.xmin; } |
109 | ||
110 | 3811117x |
double xmax() const { return m_extent.xmax; } |
111 | ||
112 | 3671541x |
double ymin() const { return m_extent.ymin; } |
113 | ||
114 | 5957169x |
double ymax() const { return m_extent.ymax; } |
115 | ||
116 | 6247212x |
double dx() const { return m_dx; } |
117 | ||
118 | 6475059x |
double dy() const { return m_dy; } |
119 | ||
120 | 2236374x |
const Box& extent() const { return m_extent; } |
121 | ||
122 | 400437x |
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 | 400437x |
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 | 26314x |
double x_for_col(size_t col) const { return m_extent.xmin + ((col - extent_tag::padding) + 0.5) * m_dx; } |
127 | ||
128 | 13356x |
double y_for_row(size_t row) const { return m_extent.ymax - ((row - extent_tag::padding) + 0.5) * m_dy; } |
129 | ||
130 | 1555x |
Grid<extent_tag> crop(const Box & b) const { |
131 | 1555x |
if (extent().intersects(b)) { |
132 | 1471x |
return shrink_to_fit(b.intersection(extent())); |
133 |
} else { |
|
134 | 84x |
return make_empty(); |
135 |
} |
|
136 |
} |
|
137 | ||
138 | 400841x |
Grid<extent_tag> shrink_to_fit(const Box & b) const { |
139 | 400841x |
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 | 400841x |
size_t col0 = get_column(b.xmin); |
144 | 400841x |
size_t row1 = get_row(b.ymax); |
145 | ||
146 |
// Shrink xmin and ymax to fit the upper-left corner of the supplied extent |
|
147 | 400841x |
double snapped_xmin = m_extent.xmin + (col0 - extent_tag::padding) * m_dx; |
148 | 400841x |
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 | 400841x |
if (b.xmin < snapped_xmin) { |
153 | ! |
snapped_xmin -= m_dx; |
154 | ! |
col0--; |
155 |
} |
|
156 | ||
157 | 400841x |
if (b.ymax > snapped_ymax) { |
158 | ! |
snapped_ymax += m_dy; |
159 | ! |
row1--; |
160 |
} |
|
161 | ||
162 | 400841x |
size_t col1 = get_column(b.xmax); |
163 | 400841x |
size_t row0 = get_row(b.ymin); |
164 | ||
165 | 400841x |
size_t num_rows = 1 + (row0 - row1); |
166 | 400841x |
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 | 400841x |
if (num_rows > 2 && (snapped_ymax - (num_rows-1)*m_dy <= b.ymin)) { |
174 | 465x |
num_rows--; |
175 |
} |
|
176 | 400841x |
if (num_cols > 2 && (snapped_xmin + (num_cols-1)*m_dx >= b.xmax)) { |
177 | 453x |
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 | 400841x |
Box reduced_box = { |
185 |
snapped_xmin, |
|
186 | 400841x |
std::min(snapped_ymax - num_rows * m_dy, b.ymin), |
187 | 400841x |
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 | 400841x |
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 | 400841x |
if (reduced_box.ymin < m_extent.ymin) { |
202 | 3x |
if (std::round((reduced_box.ymax - reduced_box.ymin)/m_dy) == |
203 | 3x |
std::round((reduced_box.ymax - m_extent.ymin)/m_dy)) { |
204 | 3x |
reduced_box.ymin = m_extent.ymin; |
205 |
} else { |
|
206 | ! |
throw std::runtime_error("Shrink operation failed."); |
207 |
} |
|
208 |
} |
|
209 | ||
210 | 400841x |
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 | 801682x |
return reduced; |
217 |
} |
|
218 | ||
219 |
template<typename extent_tag2> |
|
220 | 98x |
bool compatible_with(const Grid<extent_tag2> &b, double compatability_tol) 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 | ||
226 |
if (empty() || b.empty()) { |
|
227 | ! |
return true; |
228 |
} |
|
229 | ||
230 | 98x |
double xtol = std::min(m_dx, b.m_dx) * compatability_tol; |
231 | 98x |
double ytol = std::min(m_dy, b.m_dy) * compatability_tol; |
232 | ||
233 |
// Check x-resolution compatibility |
|
234 | 98x |
if (!is_integral(std::max(m_dx, b.m_dx) / std::min(m_dx, b.m_dx), xtol)) { |
235 | ! |
return false; |
236 |
} |
|
237 | ||
238 |
// Check y-resolution compatibility |
|
239 | 98x |
if (!is_integral(std::max(m_dy, b.m_dy) / std::min(m_dy, b.m_dy), ytol)) { |
240 | ! |
return false; |
241 |
} |
|
242 | ||
243 |
// Check left-hand boundary compatibility |
|
244 | 98x |
if (!is_integral(std::abs(b.m_extent.xmin - m_extent.xmin) / std::min(m_dx, b.m_dx), xtol)) { |
245 | 1x |
return false; |
246 |
} |
|
247 | ||
248 |
// Check upper boundary compatibility |
|
249 | 97x |
if (!is_integral(std::abs(b.m_extent.ymax - m_extent.ymax) / std::min(m_dy, b.m_dy), ytol)) { |
250 | ! |
return false; |
251 |
} |
|
252 | ||
253 | 97x |
return true; |
254 |
} |
|
255 | ||
256 |
template<typename extent_tag2> |
|
257 | 98x |
Grid<extent_tag> common_grid(const Grid<extent_tag2> &b, double tol = 1e-6) const { |
258 | 98x |
if (!compatible_with(b, tol)) { |
259 | 1x |
throw std::runtime_error("Incompatible extents."); |
260 |
} |
|
261 | ||
262 | 97x |
if (b.empty()) { |
263 | ! |
return *this; |
264 |
} |
|
265 | ||
266 | 97x |
const double common_dx = std::min(m_dx, b.m_dx); |
267 | 97x |
const double common_dy = std::min(m_dy, b.m_dy); |
268 | ||
269 | 97x |
const double common_xmin = std::min(m_extent.xmin, b.m_extent.xmin); |
270 | 97x |
const double common_ymax = std::max(m_extent.ymax, b.m_extent.ymax); |
271 | ||
272 | 97x |
double common_xmax = std::max(m_extent.xmax, b.m_extent.xmax); |
273 | 97x |
double common_ymin = std::min(m_extent.ymin, b.m_extent.ymin); |
274 | ||
275 | 97x |
const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx)); |
276 | 97x |
const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy)); |
277 | ||
278 | 97x |
common_xmax = std::max(common_xmax, common_xmin + nx*common_dx); |
279 | 97x |
common_ymin = std::min(common_ymin, common_ymax - ny*common_dy); |
280 | ||
281 | 97x |
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, double tol = 1e-6) const { |
|
286 |
if (!compatible_with(b, tol)) { |
|
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 | 399262x |
bool operator==(const Grid<extent_tag> &b) const { |
313 |
return |
|
314 | 399262x |
m_extent == b.m_extent && |
315 | 399509x |
m_dx == b.m_dx && |
316 | 399509x |
m_dy == b.m_dy; |
317 |
} |
|
318 | ||
319 | 399262x |
bool operator!=(const Grid<extent_tag> &b) const { |
320 | 399262x |
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 | 1832812x |
Matrix(size_t rows, size_t cols) : |
32 |
m_rows{rows}, |
|
33 | 1832812x |
m_cols{cols} |
34 |
{ |
|
35 | 1832812x |
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 | 341x |
Matrix(size_t rows, size_t cols, T value) : |
42 |
m_rows{rows}, |
|
43 | 341x |
m_cols{cols} |
44 |
{ |
|
45 | 341x |
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 | 341x |
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 | 517024x |
Matrix(Matrix<T>&& m) noexcept : |
66 |
m_rows{m.rows()}, |
|
67 | 517024x |
m_cols{m.cols()} |
68 |
{ |
|
69 | 517024x |
m_data = std::move(m.m_data); |
70 |
} |
|
71 | ||
72 | 237684786x |
T& operator()(size_t row, size_t col) { |
73 | 237684786x |
check(row, col); |
74 | 237684786x |
return m_data[row*m_cols + col]; |
75 |
} |
|
76 | ||
77 | 54194793x |
T operator()(size_t row, size_t col) const { |
78 | 54194793x |
check(row, col); |
79 | 54194793x |
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 | 26984725x |
void increment(size_t row, size_t col, const T & val) { |
94 | 26984725x |
check(row, col); |
95 | 26984725x |
m_data[row*m_cols + col] += val; |
96 |
} |
|
97 | ||
98 | 10022068x |
size_t rows() const { return m_rows; } |
99 | 115488381x |
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 | 318864304x |
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-2021 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 <array> |
|
18 |
#include <cmath> |
|
19 |
#include <limits> |
|
20 | ||
21 |
#include "grid.h" |
|
22 |
#include "matrix.h" |
|
23 | ||
24 |
namespace exactextract { |
|
25 |
template<typename T> |
|
26 |
class AbstractRaster { |
|
27 |
public: |
|
28 |
using value_type = T; |
|
29 | ||
30 | 917408x |
explicit AbstractRaster(const Grid<bounded_extent> & ex) : |
31 |
m_grid{ex}, |
|
32 | 917408x |
m_nodata{std::is_floating_point<T>::value ? std::numeric_limits<T>::quiet_NaN() : std::numeric_limits<T>::min()}, |
33 | 917408x |
m_has_nodata{false} |
34 |
{} |
|
35 | ||
36 |
AbstractRaster(const Grid<bounded_extent> & ex, const T& nodata_val) : m_grid{ex}, m_nodata{nodata_val}, m_has_nodata{true} {} |
|
37 | ||
38 | 917408x |
virtual ~AbstractRaster() = default; |
39 | ||
40 | 58719556x |
size_t rows() const { |
41 | 58719556x |
return m_grid.rows(); |
42 |
} |
|
43 | ||
44 | 210677100x |
size_t cols() const { |
45 | 210677100x |
return m_grid.cols(); |
46 |
} |
|
47 | ||
48 |
size_t size() const { |
|
49 |
return rows() * cols(); |
|
50 |
} |
|
51 | ||
52 | 399125x |
double xres() const { |
53 | 399125x |
return m_grid.dx(); |
54 |
} |
|
55 | ||
56 | 399125x |
double yres() const { |
57 | 399125x |
return m_grid.dy(); |
58 |
} |
|
59 | ||
60 | 399125x |
double xmin() const { |
61 | 399125x |
return m_grid.xmin(); |
62 |
} |
|
63 | ||
64 |
double ymin() const { |
|
65 |
return m_grid.ymin(); |
|
66 |
} |
|
67 | ||
68 |
double xmax() const { |
|
69 |
return m_grid.xmax(); |
|
70 |
} |
|
71 | ||
72 | 399125x |
double ymax() const { |
73 | 399125x |
return m_grid.ymax(); |
74 |
} |
|
75 | ||
76 | 61204366x |
const Grid<bounded_extent>& grid() const { |
77 | 61204366x |
return m_grid; |
78 |
} |
|
79 | ||
80 |
virtual T operator()(size_t row, size_t col) const = 0; |
|
81 | ||
82 | 399142x |
bool has_nodata() const { return m_has_nodata; } |
83 | ||
84 | 209318x |
T nodata() const { return m_nodata; } |
85 | ||
86 | 9478416x |
bool get(size_t row, size_t col, T & val) const { |
87 | 9478416x |
val = operator()(row, col); |
88 | ||
89 | ! |
if (m_has_nodata && val == m_nodata) { |
90 | ! |
return false; |
91 |
} |
|
92 | 9478416x |
if (std::is_floating_point<T>::value && std::isnan(val)) { |
93 | 23x |
return false; |
94 |
} |
|
95 | ||
96 | 9478393x |
return true; |
97 |
} |
|
98 | ||
99 | ! |
void set_nodata(const T & val) { |
100 | ! |
m_has_nodata = true; |
101 | ! |
m_nodata = val; |
102 |
} |
|
103 | ||
104 |
bool operator==(const AbstractRaster<T> & other) const { |
|
105 |
if (rows() != other.rows()) |
|
106 |
return false; |
|
107 | ||
108 |
if (cols() != other.cols()) |
|
109 |
return false; |
|
110 | ||
111 |
if (xres() != other.xres()) |
|
112 |
return false; |
|
113 | ||
114 |
if (yres() != other.yres()) |
|
115 |
return false; |
|
116 | ||
117 |
if (xmin() != other.xmin()) |
|
118 |
return false; |
|
119 | ||
120 |
if (ymin() != other.ymin()) |
|
121 |
return false; |
|
122 | ||
123 |
// Do the rasters differ in their definition of NODATA? If so, we need to do some more detailed |
|
124 |
// checking below. |
|
125 |
bool nodata_differs = has_nodata() != other.has_nodata() || |
|
126 |
(has_nodata() && other.has_nodata() && nodata() != other.nodata()); |
|
127 | ||
128 |
for (size_t i = 0; i < rows(); i++) { |
|
129 |
for (size_t j = 0; j < cols(); j++) { |
|
130 |
if(operator()(i, j) != other(i, j)) { |
|
131 |
// Override default behavior of NAN != NAN |
|
132 |
if (!std::isnan(operator()(i, j)) || !std::isnan(other(i, j))) { |
|
133 |
return false; |
|
134 |
} |
|
135 |
} else if (nodata_differs && (operator()(i, j) == nodata() || other(i, j) == other.nodata())) { |
|
136 |
// For data types that do not have NAN, or even for floating point types where the user has |
|
137 |
// selected some other value to represent NODATA, we need to reverse a positive equality test |
|
138 |
// where the value is considered to be NODATA in one raster but not the other. |
|
139 |
return false; |
|
140 |
} |
|
141 |
} |
|
142 |
} |
|
143 | ||
144 |
return true; |
|
145 |
} |
|
146 | ||
147 |
bool operator!=(const AbstractRaster<T> & other) const { |
|
148 |
return !(operator==(other)); |
|
149 |
} |
|
150 | ||
151 |
class Iterator : public std::iterator<std::forward_iterator_tag, T> { |
|
152 |
public: |
|
153 |
Iterator(const AbstractRaster<T>* r, size_t i, size_t j) : |
|
154 |
m_rast(r), m_row(i), m_col(j) {} |
|
155 | ||
156 |
const T& operator*() const { |
|
157 |
m_val = m_rast->operator()(m_row, m_col); |
|
158 |
return m_val; |
|
159 |
} |
|
160 | ||
161 |
Iterator& operator++() { |
|
162 |
m_col++; |
|
163 |
if (m_col == m_rast->cols()) { |
|
164 |
m_col = 0; |
|
165 |
m_row++; |
|
166 |
} |
|
167 |
return *this; |
|
168 |
} |
|
169 | ||
170 |
friend bool operator==(const Iterator& a, |
|
171 |
const Iterator& b) { |
|
172 |
return a.m_rast == b.m_rast && a.m_row == b.m_row && a.m_col == b.m_col; |
|
173 |
} |
|
174 | ||
175 |
friend bool operator!=(const Iterator& a, |
|
176 |
const Iterator& b) { |
|
177 |
return a.m_row != b.m_row || a.m_col != b.m_col || a.m_rast != b.m_rast; |
|
178 |
} |
|
179 | ||
180 |
private: |
|
181 |
const AbstractRaster<T>* m_rast; |
|
182 |
mutable T m_val; |
|
183 |
size_t m_row; |
|
184 |
size_t m_col; |
|
185 |
}; |
|
186 | ||
187 |
Iterator begin() const { |
|
188 |
return Iterator(this, 0, 0); |
|
189 |
} |
|
190 | ||
191 |
Iterator end() const { |
|
192 |
return Iterator(this, rows(), 0); |
|
193 |
} |
|
194 |
private: |
|
195 |
Grid<bounded_extent> m_grid; |
|
196 |
T m_nodata; |
|
197 |
bool m_has_nodata; |
|
198 |
}; |
|
199 | ||
200 |
template<typename T> |
|
201 |
class Raster : public AbstractRaster<T> { |
|
202 |
public: |
|
203 |
Raster(Matrix<T>&& values, const Box & box) : AbstractRaster<T>( |
|
204 |
Grid<bounded_extent>( |
|
205 |
box, |
|
206 |
(box.xmax - box.xmin) / values.cols(), |
|
207 |
(box.ymax - box.ymin) / values.rows())), |
|
208 |
m_values{std::move(values)} {} |
|
209 | ||
210 | 517024x |
Raster(Matrix<T>&& values, const Grid<bounded_extent> & g) : |
211 |
AbstractRaster<T>(g), |
|
212 | 517024x |
m_values{std::move(values)} {} |
213 | ||
214 |
Raster(const Box & box, size_t nrow, size_t ncol) : |
|
215 |
AbstractRaster<T>(Grid<bounded_extent>(box, (box.xmax-box.xmin) / ncol, (box.ymax-box.ymin) / nrow)), |
|
216 |
m_values{nrow, ncol} |
|
217 |
{} |
|
218 | ||
219 |
explicit Raster(const Grid<bounded_extent> & ex) : |
|
220 |
AbstractRaster<T>(ex), |
|
221 |
m_values{ex.rows(), ex.cols()} |
|
222 |
{} |
|
223 | ||
224 |
Matrix<T>& data() { |
|
225 |
return m_values; |
|
226 |
} |
|
227 | ||
228 | 504x |
T& operator()(size_t row, size_t col) { |
229 | 504x |
return m_values(row, col); |
230 |
} |
|
231 | ||
232 | 27210069x |
T operator()(size_t row, size_t col) const override { |
233 | 27210069x |
return m_values(row, col); |
234 |
} |
|
235 | ||
236 |
private: |
|
237 |
Matrix<T> m_values; |
|
238 |
}; |
|
239 | ||
240 |
template<typename T> |
|
241 |
class RasterView : public AbstractRaster<T>{ |
|
242 |
public: |
|
243 |
// Construct a view of a raster r at an extent ex that is larger |
|
244 |
// and/or of finer resolution than r |
|
245 | 399142x |
RasterView(const AbstractRaster<T> & r, Grid<bounded_extent> ex) : |
246 |
AbstractRaster<T>(ex), |
|
247 |
m_raster{r}, |
|
248 |
m_x_off{0}, |
|
249 |
m_y_off{0}, |
|
250 |
m_rx{1}, |
|
251 | 399142x |
m_ry{1} { |
252 | 399142x |
if (!this->grid().empty()) { |
253 | 399125x |
double disaggregation_factor_x = r.xres() / ex.dx(); |
254 | 399125x |
double disaggregation_factor_y = r.yres() / ex.dy(); |
255 | ||
256 | ! |
if (std::abs(disaggregation_factor_x - std::round(disaggregation_factor_x)) > 1e-6 || |
257 | 399125x |
std::abs(disaggregation_factor_y - std::round(disaggregation_factor_y)) > 1e-6) { |
258 | ! |
throw std::runtime_error("Must construct view at resolution that is an integer multiple of original."); |
259 |
} |
|
260 | ||
261 | 399125x |
if (disaggregation_factor_x < 0 || disaggregation_factor_y < 0) { |
262 | ! |
throw std::runtime_error("Must construct view at equal or higher resolution than original."); |
263 |
} |
|
264 | ||
265 | 399125x |
m_x_off = static_cast<long>(std::round((ex.xmin() - r.xmin()) / ex.dx())); |
266 | 399125x |
m_y_off = static_cast<long>(std::round((r.ymax() - ex.ymax()) / ex.dy())); |
267 | 399125x |
m_rx = static_cast<size_t>(std::round(disaggregation_factor_x)); |
268 | 399125x |
m_ry = static_cast<size_t>(std::round(disaggregation_factor_y)); |
269 |
} |
|
270 | ||
271 | 399142x |
if (r.has_nodata()) { |
272 | ! |
this->set_nodata(r.nodata()); |
273 |
} |
|
274 |
} |
|
275 | ||
276 | 28488528x |
T operator()(size_t row, size_t col) const override { |
277 | 28488528x |
if (m_raster.grid().empty()) { |
278 | 8x |
return this->nodata(); |
279 |
} |
|
280 | ||
281 | 237313x |
if (m_x_off < 0 && static_cast<size_t>(-m_x_off) > col) { |
282 | 50432x |
return this->nodata(); |
283 |
} |
|
284 | 186881x |
if (m_y_off < 0 && static_cast<size_t>(-m_y_off) > row) { |
285 | 28670x |
return this->nodata(); |
286 |
} |
|
287 | ||
288 | 28409418x |
size_t i0 = (row + m_y_off) / m_ry; |
289 | 28409418x |
size_t j0 = (col + m_x_off) / m_rx; |
290 | ||
291 | 28409418x |
if (i0 > m_raster.rows() - 1 || j0 > m_raster.cols() - 1) { |
292 | 130208x |
return this->nodata(); |
293 |
} |
|
294 | ||
295 | 28279210x |
return m_raster(i0, j0); |
296 |
} |
|
297 | ||
298 |
private: |
|
299 |
const AbstractRaster<T>& m_raster; |
|
300 | ||
301 |
long m_x_off; |
|
302 |
long m_y_off; |
|
303 |
size_t m_rx; |
|
304 |
size_t m_ry; |
|
305 |
}; |
|
306 | ||
307 | ||
308 |
template<typename T> |
|
309 |
std::ostream& operator<<(std::ostream & os, const AbstractRaster<T> & m) { |
|
310 |
for (size_t i = 0; i < m.rows(); i++) { |
|
311 |
for (size_t j = 0; j < m.cols(); j++) { |
|
312 |
if (m(i, j) != 0) { |
|
313 |
os << std::right << std::fixed << std::setw(10) << std::setprecision(6) << |
|
314 |
m(i, j) << " "; |
|
315 |
} else { |
|
316 |
os << " "; |
|
317 |
} |
|
318 |
} |
|
319 |
os << std::endl; |
|
320 |
} |
|
321 | ||
322 |
return os; |
|
323 |
} |
|
324 |
} |
|
325 | ||
326 |
#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 | 380x |
GEOSAutoHandle() { |
55 | 380x |
handle = initGEOS_r(geos_warn, geos_error); |
56 |
} |
|
57 | ||
58 | 760x |
~GEOSAutoHandle() { |
59 | 380x |
finishGEOS_r(handle); |
60 |
} |
|
61 | ||
62 |
GEOSContextHandle_t handle; |
|
63 |
}; |
|
64 | ||
65 |
// Return a smart pointer to a Geometry, given WKB input |
|
66 | 379x |
static inline geom_ptr read_wkb(const GEOSContextHandle_t & context, const Rcpp::RawVector & wkb) { |
67 | 1137x |
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 | 758x |
std::addressof(wkb[0]), |
72 | 379x |
wkb.size()), |
73 | 758x |
[context](GEOSGeometry* g) { GEOSGeom_destroy_r(context, g); } }; |
74 | ||
75 | 379x |
if (geom.get() == nullptr) { |
76 | ! |
Rcpp::stop("Failed to parse WKB geometry"); |
77 |
} |
|
78 | ||
79 | 758x |
return geom; |
80 |
} |
1 |
// Copyright (c) 2018-2022 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 <Rcpp.h> |
|
16 | ||
17 |
#include <memory> |
|
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 | 153x |
Rcpp::List CPP_exact_extract(Rcpp::S4 & rast, |
47 |
Rcpp::Nullable<Rcpp::S4> & rast_uncropped, |
|
48 |
Rcpp::Nullable<Rcpp::S4> & weights, |
|
49 |
const Rcpp::RawVector & wkb, |
|
50 |
double default_value, |
|
51 |
double default_weight, |
|
52 |
bool include_xy, |
|
53 |
bool include_cell_number, |
|
54 |
bool include_area, |
|
55 |
bool area_weights, |
|
56 |
bool coverage_areas, |
|
57 |
Rcpp::Nullable<Rcpp::CharacterVector> & p_area_method, |
|
58 |
Rcpp::Nullable<Rcpp::List> & include_cols, |
|
59 |
Rcpp::CharacterVector & src_names, |
|
60 |
Rcpp::Nullable<Rcpp::CharacterVector> & p_weights_names, |
|
61 |
bool warn_on_disaggregate, |
|
62 |
double grid_compat_tol) { |
|
63 |
try { |
|
64 | 306x |
GEOSAutoHandle geos; |
65 | 459x |
Rcpp::Function names("names"); |
66 | ||
67 | 153x |
auto grid = make_grid(rast); |
68 | 153x |
auto weights_grid = exactextract::Grid<bounded_extent>::make_empty(); |
69 | 153x |
auto common_grid = grid; |
70 | ||
71 | 306x |
S4RasterSource rsrc(rast, default_value); |
72 | 153x |
int src_nlayers = get_nlayers(rast); |
73 | ||
74 | 153x |
std::unique_ptr<S4RasterSource> rweights; |
75 | 153x |
int weights_nlayers = 0; |
76 | 306x |
Rcpp::CharacterVector weights_names; |
77 | ||
78 | 306x |
std::string area_method; |
79 | 153x |
if (p_area_method.isNotNull()) { |
80 | 11x |
Rcpp::CharacterVector amethod = p_area_method.get(); |
81 | 11x |
area_method = amethod[0]; |
82 |
} |
|
83 | ||
84 | 153x |
if (weights.isNotNull()) { |
85 | 86x |
Rcpp::S4 weights_s4 = weights.get(); |
86 | 43x |
weights_nlayers = get_nlayers(weights_s4); |
87 | 43x |
weights_grid = make_grid(weights_s4); |
88 | ||
89 | 43x |
common_grid = grid.common_grid(weights_grid, grid_compat_tol); |
90 | ||
91 | 43x |
rweights = std::make_unique<S4RasterSource>(weights_s4, default_weight); |
92 | 43x |
weights_names = p_weights_names.get(); |
93 | ||
94 | 43x |
if (warn_on_disaggregate && (common_grid.dx() < grid.dx() || common_grid.dy() < grid.dy())) { |
95 | 3x |
Rcpp::warning("value raster implicitly disaggregated to match higher resolution of weights"); |
96 |
} |
|
97 | 110x |
} else if (area_weights) { |
98 | 1x |
weights_nlayers = 1; |
99 | 1x |
weights_names = p_weights_names.get(); |
100 |
} |
|
101 | ||
102 | 306x |
auto geom = read_wkb(geos.handle, wkb); |
103 | ||
104 | 153x |
auto bbox = exactextract::geos_get_box(geos.handle, geom.get()); |
105 | ||
106 | 153x |
common_grid = common_grid.crop(bbox); |
107 | ||
108 | 306x |
auto coverage_fractions = raster_cell_intersection(common_grid, geos.handle, geom.get()); |
109 | 153x |
auto& cov_grid = coverage_fractions.grid(); |
110 | ||
111 |
// We can't construct an Rcpp::DataFrame with a variable number of columns |
|
112 |
// because of a bug in Rcpp (see https://github.com/RcppCore/Rcpp/pull/1099) |
|
113 |
// Instead, we build up a list, and then create a data frame from the list. |
|
114 |
// Profiling shows that this ends up in as.data.frame, which is slow. |
|
115 |
// Once Rcpp 1.0.6 is released, this can be reworked to be simpler and faster. |
|
116 | 306x |
Rcpp::List cols; |
117 | ||
118 | 306x |
Rcpp::NumericVector coverage_vec = as_vector(coverage_fractions); |
119 | 153x |
if (coverage_areas) { |
120 | 8x |
auto areas = get_area_raster(area_method, cov_grid); |
121 | 4x |
Rcpp::NumericVector area_vec = as_vector(*areas); |
122 | 4x |
coverage_vec = coverage_vec * area_vec; |
123 |
} |
|
124 | 306x |
Rcpp::LogicalVector covered = coverage_vec > 0; |
125 | ||
126 | 153x |
if (include_cols.isNotNull()) { |
127 | 10x |
Rcpp::List include_cols_list = include_cols.get(); |
128 | 10x |
Rcpp::CharacterVector include_names = include_cols_list.attr("names"); |
129 | ||
130 | 10x |
for (int i = 0; i < include_names.size(); i++) { |
131 | 5x |
std::string name(include_names[i]); |
132 | 5x |
cols[name] = include_cols_list[name]; |
133 |
} |
|
134 |
} |
|
135 | ||
136 | 341x |
for (int i = 0; i < src_nlayers; i++) { |
137 | 376x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
138 | 188x |
const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get()); |
139 | ||
140 |
// TODO Perhaps extend this to preserve types (integer, logical.) |
|
141 |
// A bit challenging, since we don't know the type of the raster |
|
142 |
// until we call getValuesBlock, and even then we can't be sure |
|
143 |
// that the returned type is correct (as in a RasterStack with mixed types.) |
|
144 |
// Since R integers are only 32-bit, we are not going to lose data by |
|
145 |
// converting everything to numeric, although we pay a storage penalty. |
|
146 | 188x |
Rcpp::NumericVector value_vec = r->vec(); |
147 | 356x |
if (grid.dx() != cov_grid.dx() || grid.dy() != cov_grid.dy() || |
148 | 168x |
value_vec.size() != covered.size()) { |
149 |
// Transform values to same grid as coverage fractions |
|
150 | 40x |
RasterView<double> rt(*r, cov_grid); |
151 | 40x |
value_vec = as_vector(rt); |
152 |
} |
|
153 | ||
154 | 188x |
value_vec = value_vec[covered]; |
155 | 188x |
cols[std::string(src_names[i])] = value_vec; |
156 |
} |
|
157 | ||
158 | 212x |
for (int i = 0; i < weights_nlayers; i++) { |
159 | 118x |
Rcpp::NumericVector weight_vec; |
160 | 59x |
if (area_weights) { |
161 | 1x |
auto weights = get_area_raster(area_method, cov_grid); |
162 | 1x |
weight_vec = as_vector(*weights); |
163 |
} else { |
|
164 | 116x |
auto values = rweights->read_box(cov_grid.extent(), i); |
165 | 58x |
const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get()); |
166 | 58x |
weight_vec = r->vec(); |
167 | ||
168 | 96x |
if (weights_grid.dx() != cov_grid.dx() || weights_grid.dy() != cov_grid.dy() || |
169 | 38x |
weight_vec.size() != covered.size()) { |
170 |
// Transform weights to same grid as coverage fractions |
|
171 | 22x |
RasterView<double> rt (*r, cov_grid); |
172 | ||
173 | 22x |
weight_vec = as_vector(rt); |
174 |
} |
|
175 |
} |
|
176 | ||
177 | 59x |
weight_vec = weight_vec[covered]; |
178 | ||
179 | 59x |
std::string colname(weights_names[i]); |
180 | 59x |
if (cols.containsElementNamed(colname.c_str())) { |
181 |
// append ".1" to the column name, to match the behavior of |
|
182 |
// data.frame(). We're safe to just add ".1" (instead of incrementing |
|
183 |
// .1 to .2, for example) because duplicated names within the values |
|
184 |
// or weight stack will already have been made unique by raster::stack() |
|
185 | 2x |
colname = colname + ".1"; // append .1 |
186 |
} |
|
187 | 59x |
cols[colname] = weight_vec; |
188 |
} |
|
189 | ||
190 | 153x |
if (include_xy) { |
191 |
// Include xy values from whichever input raster has the higher resolution |
|
192 | 30x |
if (weights_nlayers > 0 && (weights_grid.dx() < grid.dx() || weights_grid.dy() < grid.dy())) { |
193 | 1x |
Rcpp::S4 weights_s4 = weights.get(); |
194 | 1x |
cols["x"] = get_x_values(weights_s4, cov_grid)[covered]; |
195 | 1x |
cols["y"] = get_y_values(weights_s4, cov_grid)[covered]; |
196 |
} else { |
|
197 | 29x |
cols["x"] = get_x_values(rast, cov_grid)[covered]; |
198 | 29x |
cols["y"] = get_y_values(rast, cov_grid)[covered]; |
199 |
} |
|
200 |
} |
|
201 | ||
202 | 153x |
if (include_cell_number) { |
203 | 24x |
if (rast_uncropped.isNull()) { |
204 | 12x |
cols["cell"] = get_cell_numbers(rast, cov_grid)[covered]; |
205 |
} else { |
|
206 | 12x |
Rcpp::S4 tmp = rast_uncropped.get(); |
207 | 12x |
cols["cell"] = get_cell_numbers(tmp, cov_grid)[covered]; |
208 |
} |
|
209 |
} |
|
210 | ||
211 | 153x |
if (include_area) { |
212 | 12x |
auto area_rast = get_area_raster(area_method, cov_grid); |
213 | 6x |
Rcpp::NumericVector area_vec = as_vector(*area_rast); |
214 | 6x |
cols["area"] = area_vec[covered]; |
215 |
} |
|
216 | ||
217 | 153x |
if (coverage_areas) { |
218 | 4x |
cols["coverage_area"] = coverage_vec[covered]; |
219 |
} else { |
|
220 | 149x |
cols["coverage_fraction"] = coverage_vec[covered]; |
221 |
} |
|
222 | ||
223 | 306x |
return cols; |
224 |
} catch (std::exception & e) { |
|
225 |
// throw predictable exception class |
|
226 |
#ifdef __SUNPRO_CC |
|
227 |
// Rcpp::stop crashes CRAN Solaris build |
|
228 |
// https://github.com/RcppCore/Rcpp/issues/1159 |
|
229 |
Rf_error(e.what()); |
|
230 |
return R_NilValue; |
|
231 |
#else |
|
232 |
Rcpp::stop(e.what()); |
|
233 |
#endif |
|
234 |
} |
|
235 |
} |
|
236 | ||
237 |
enum class WeightingMethod { |
|
238 |
NONE, |
|
239 |
RASTER, |
|
240 |
AREA |
|
241 |
}; |
|
242 | ||
243 | ||
244 | 208x |
static int get_num_stats(const Rcpp::StringVector & stats, |
245 |
const std::size_t num_quantiles, |
|
246 |
const std::size_t num_unique_values) { |
|
247 | 208x |
int num_stats = 0; |
248 | 452x |
for (const auto & stat : stats) { |
249 | 724x |
if (stat == std::string("frac") || |
250 | 480x |
stat == std::string("weighted_frac")) { |
251 | 10x |
num_stats += static_cast<int>(num_unique_values); |
252 | 234x |
} else if (stat == std::string("quantile")) { |
253 | 12x |
num_stats += static_cast<int>(num_quantiles); |
254 |
} else { |
|
255 | 222x |
num_stats += 1; |
256 |
} |
|
257 |
} |
|
258 | ||
259 | 208x |
return num_stats; |
260 |
} |
|
261 | ||
262 |
// Return a matrix with one row per stat and one row per raster layer |
|
263 |
// [[Rcpp::export]] |
|
264 | 214x |
Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast, |
265 |
Rcpp::Nullable<Rcpp::S4> weights, |
|
266 |
const Rcpp::RawVector & wkb, |
|
267 |
double default_value, |
|
268 |
double default_weight, |
|
269 |
bool coverage_areas, |
|
270 |
Rcpp::Nullable<Rcpp::CharacterVector> & p_area_method, |
|
271 |
const Rcpp::StringVector & stats, |
|
272 |
int max_cells_in_memory, |
|
273 |
double grid_compat_tol, |
|
274 |
const Rcpp::Nullable<Rcpp::NumericVector> & quantiles) { |
|
275 |
try { |
|
276 | 428x |
GEOSAutoHandle geos; |
277 | ||
278 | 214x |
if (max_cells_in_memory < 1) { |
279 | 1x |
Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory); |
280 |
} |
|
281 | ||
282 | 213x |
int nlayers = get_nlayers(rast); |
283 | ||
284 | 426x |
S4RasterSource rsrc(rast, default_value); |
285 | ||
286 | 213x |
std::unique_ptr<S4RasterSource> rweights; |
287 | 426x |
std::string area_method; |
288 | 213x |
if (p_area_method.isNotNull()) { |
289 | 10x |
area_method = ((Rcpp::CharacterVector) p_area_method.get())[0]; |
290 |
} |
|
291 | ||
292 | 213x |
WeightingMethod weighting = WeightingMethod::NONE; |
293 | 213x |
int nweights = 0; |
294 | 213x |
if (weights.isNotNull()) { |
295 | 55x |
weighting = WeightingMethod::RASTER; |
296 | 55x |
Rcpp::S4 weights_s4 = weights.get(); |
297 | 55x |
nweights = get_nlayers(weights_s4); |
298 | ||
299 | 55x |
if (nlayers > 1 && nweights > 1 && nlayers != nweights) { |
300 | ! |
Rcpp::stop("Incompatible number of layers in value and weighting rasters"); |
301 |
} |
|
302 | ||
303 | 55x |
rweights = std::make_unique<S4RasterSource>(weights_s4, default_weight); |
304 | 158x |
} else if (p_area_method.isNotNull()) { |
305 | 9x |
weighting = WeightingMethod::AREA; |
306 | 9x |
nweights = 1; |
307 |
} |
|
308 | ||
309 | 426x |
auto geom = read_wkb(geos.handle, wkb); |
310 | 213x |
auto bbox = exactextract::geos_get_box(geos.handle, geom.get()); |
311 | ||
312 |
auto grid = weighting == WeightingMethod::RASTER ? |
|
313 | 213x |
rsrc.grid().common_grid(rweights->grid(), grid_compat_tol) : rsrc.grid(); |
314 | ||
315 | 212x |
bool disaggregated = (grid.dx() < rsrc.grid().dx() || grid.dy() < rsrc.grid().dy()); |
316 | ||
317 | 212x |
bool store_values = false; |
318 | 212x |
bool calc_value_set = false; |
319 | 212x |
std::size_t num_unique_values = 0; |
320 | 212x |
std::size_t num_quantiles = 0; |
321 | ||
322 | 458x |
for (const auto & stat : stats) { |
323 | 250x |
store_values = store_values || requires_stored_values(stat); |
324 | ||
325 | 508x |
if (disaggregated && (stat == std::string("count") || |
326 | 258x |
stat == std::string("sum"))) { |
327 | 2x |
Rcpp::stop("Cannot compute 'count' or 'sum' when value raster is disaggregated to resolution of weights."); |
328 |
} |
|
329 | ||
330 | 736x |
if (stat == std::string("frac") || |
331 | 488x |
stat == std::string("weighted_frac")) { |
332 | 10x |
calc_value_set = true; |
333 |
} |
|
334 | ||
335 | 248x |
if (stat == std::string("quantile")) { |
336 | 14x |
if (quantiles.isNotNull()) { |
337 | 13x |
Rcpp::NumericVector qvec = quantiles.get(); |
338 | 13x |
num_quantiles = qvec.size(); |
339 |
} |
|
340 | ||
341 | 14x |
if (num_quantiles == 0) { |
342 | 2x |
Rcpp::stop("Quantiles not specified."); |
343 |
} |
|
344 |
} |
|
345 |
} |
|
346 | ||
347 | 208x |
int nresults = std::max(nlayers, nweights); |
348 | ||
349 | 416x |
std::vector<RasterStats<double>> raster_stats; |
350 | 208x |
raster_stats.reserve(nresults); |
351 | 431x |
for (int i = 0; i < nresults; i++) { |
352 | 223x |
raster_stats.emplace_back(store_values); |
353 |
} |
|
354 | ||
355 | 208x |
if (bbox.intersects(grid.extent())) { |
356 | 199x |
auto cropped_grid = grid.crop(bbox); |
357 | ||
358 | 437x |
for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) { |
359 | 476x |
auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get()); |
360 | 238x |
if (coverage_areas) { |
361 | 4x |
auto areas = get_area_raster(area_method, subgrid); |
362 | 6x |
for (size_t i = 0; i < coverage_fraction.rows(); i++) { |
363 | 12x |
for (size_t j = 0; j < coverage_fraction.cols(); j++) { |
364 | 8x |
coverage_fraction(i, j) = coverage_fraction(i, j) * (*areas)(i, j); |
365 |
} |
|
366 |
} |
|
367 |
} |
|
368 | ||
369 | 238x |
auto& cov_grid = coverage_fraction.grid(); |
370 | ||
371 | 238x |
if (!cov_grid.empty()) { |
372 | 237x |
if (weighting == WeightingMethod::NONE) { |
373 | 360x |
for (int i = 0; i < nlayers; i++) { |
374 | 368x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
375 | 184x |
raster_stats[i].process(coverage_fraction, *values); |
376 |
} |
|
377 | 61x |
} else if (weighting == WeightingMethod::AREA) { |
378 | 18x |
auto weights = get_area_raster(area_method, cov_grid); |
379 | ||
380 | 18x |
for (int i = 0; i < nlayers; i++) { |
381 | 18x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
382 | 9x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
383 |
} |
|
384 |
} else { |
|
385 | 52x |
if (nlayers > nweights) { |
386 |
// recycle weights |
|
387 | 4x |
auto weights = rweights->read_box(cov_grid.extent(), 0); |
388 | ||
389 | 7x |
for (int i = 0; i < nlayers; i++) { |
390 | 10x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
391 | 5x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
392 |
} |
|
393 | 50x |
} else if (nweights > nlayers) { |
394 |
// recycle values |
|
395 | 2x |
auto values = rsrc.read_box(cov_grid.extent(), 0); |
396 | ||
397 | 4x |
for (int i = 0; i < nweights; i++) { |
398 | 6x |
auto weights = rweights->read_box(cov_grid.extent(), i); |
399 | 3x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
400 |
} |
|
401 |
} else { |
|
402 |
// process values and weights in parallel |
|
403 | 100x |
for (int i = 0; i < nlayers; i++) { |
404 | 102x |
auto values = rsrc.read_box(cov_grid.extent(), i); |
405 | 102x |
auto weights = rweights->read_box(cov_grid.extent(), i); |
406 | ||
407 | 51x |
raster_stats[i].process(coverage_fraction, *values, *weights); |
408 |
} |
|
409 |
} |
|
410 |
} |
|
411 |
} |
|
412 |
} |
|
413 |
} |
|
414 | ||
415 | ||
416 | 416x |
std::set<double> value_set; |
417 | 208x |
if (calc_value_set) { |
418 | 22x |
for (const auto& rs : raster_stats) { |
419 | 36x |
for (const auto& value : rs) { |
420 | 24x |
value_set.insert(value); |
421 |
} |
|
422 |
} |
|
423 | 10x |
num_unique_values = value_set.size(); |
424 |
} |
|
425 | ||
426 | 208x |
auto stat_result_size = get_num_stats(stats, num_quantiles, num_unique_values); |
427 | 416x |
Rcpp::NumericMatrix stat_results = Rcpp::no_init(nresults, stat_result_size); |
428 | ||
429 | 208x |
if (calc_value_set) { |
430 | 10x |
Rcpp::NumericVector unique_values(value_set.begin(), value_set.end()); |
431 | 10x |
stat_results.attr("unique_values") = unique_values; |
432 |
} |
|
433 | ||
434 |
// rows (j) represent layers |
|
435 |
// cols (i) represent stats |
|
436 | 427x |
for (int j = 0; j < nresults; j++) { |
437 | 223x |
const auto& rs = raster_stats[j]; |
438 | ||
439 | 223x |
int i = 0; |
440 | 490x |
for(const auto& stat : stats) { |
441 | 271x |
if (stat == std::string("mean")) stat_results(j, i++) = rs.mean(); |
442 | ||
443 | 233x |
else if (stat == std::string("sum")) stat_results(j, i++) = rs.sum(); |
444 | 172x |
else if (stat == std::string("count")) stat_results(j, i++) = rs.count(); |
445 | ||
446 | 155x |
else if (stat == std::string("min")) stat_results(j, i++) = rs.min().value_or(NA_REAL); |
447 | 149x |
else if (stat == std::string("max")) stat_results(j, i++) = rs.max().value_or(NA_REAL); |
448 | ||
449 | 143x |
else if (stat == std::string("median")) stat_results(j, i++) = rs.quantile(0.5).value_or(NA_REAL); |
450 | ||
451 | 142x |
else if (stat == std::string("mode")) stat_results(j, i++) = rs.mode().value_or(NA_REAL); |
452 | 107x |
else if (stat == std::string("majority")) stat_results(j, i++) = rs.mode().value_or(NA_REAL); |
453 | 104x |
else if (stat == std::string("minority")) stat_results(j, i++) = rs.minority().value_or(NA_REAL); |
454 | ||
455 | 100x |
else if (stat == std::string("variety")) stat_results(j, i++) = rs.variety(); |
456 | 88x |
else if (stat == std::string("weighted_mean")) stat_results(j, i++) = rs.weighted_mean(); |
457 | 41x |
else if (stat == std::string("weighted_sum")) stat_results(j, i++) = rs.weighted_sum(); |
458 | ||
459 | 31x |
else if (stat == std::string("variance")) stat_results(j, i++) = rs.variance(); |
460 | 30x |
else if (stat == std::string("stdev")) stat_results(j, i++) = rs.stdev(); |
461 | 29x |
else if (stat == std::string("coefficient_of_variation")) stat_results(j, i++) = rs.coefficient_of_variation(); |
462 | ||
463 | 28x |
else if (stat == std::string("quantile")) { |
464 | 28x |
Rcpp::NumericVector qvec = quantiles.get(); |
465 | 39x |
for (double q : qvec) { |
466 | 27x |
stat_results(j, i++) = rs.quantile(q).value_or(NA_REAL); |
467 |
} |
|
468 |
} |
|
469 | ||
470 | 14x |
else if (stat == std::string("frac")) { |
471 | 34x |
for (double v : value_set) { |
472 | 24x |
stat_results(j, i++) = rs.frac(v).value_or(0); |
473 |
} |
|
474 |
} |
|
475 | ||
476 | 4x |
else if (stat == std::string("weighted_frac")) { |
477 | 6x |
for (double v : value_set) { |
478 | 4x |
stat_results(j, i++) = rs.weighted_frac(v).value_or(0); |
479 |
} |
|
480 |
} |
|
481 | ||
482 | 8x |
else Rcpp::stop("Unknown stat: " + stat); |
483 |
} |
|
484 |
} |
|
485 | ||
486 |
// FIXME set dimnames here |
|
487 | ||
488 | 408x |
return stat_results; |
489 | 20x |
} catch (std::exception & e) { |
490 |
// throw predictable exception class |
|
491 |
#ifdef __SUNPRO_CC |
|
492 |
// Rcpp::stop crashes CRAN Solaris build |
|
493 |
// https://github.com/RcppCore/Rcpp/issues/1159 |
|
494 |
Rf_error(e.what()); |
|
495 |
return R_NilValue; |
|
496 |
#else |
|
497 | 10x |
Rcpp::stop(e.what()); |
498 |
#endif |
|
499 |
} |
|
500 |
} |
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 | 8934600x |
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 | 8934600x |
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 | 404x |
static Box make_empty() { |
47 | 404x |
return {0, 0, 0, 0}; |
48 |
} |
|
49 | ||
50 | 6973418x |
double width() const { |
51 | 6973418x |
return xmax - xmin; |
52 |
} |
|
53 | ||
54 | 6973418x |
double height() const { |
55 | 6973418x |
return ymax - ymin; |
56 |
} |
|
57 | ||
58 | 6568936x |
double area() const { |
59 | 6568936x |
return width() * height(); |
60 |
} |
|
61 | ||
62 | 202241x |
double perimeter() const { |
63 | 202241x |
return 2 * width() + 2 * height(); |
64 |
} |
|
65 | ||
66 | 401621x |
bool intersects(const Box & other) const { |
67 | 401621x |
if (other.ymin > ymax) |
68 | 42x |
return false; |
69 | 401579x |
if (other.ymax < ymin) |
70 | 51x |
return false; |
71 | 401528x |
if (other.xmin > xmax) |
72 | ! |
return false; |
73 | 401528x |
if (other.xmax < xmin) |
74 | ! |
return false; |
75 | ||
76 | 401528x |
return true; |
77 |
} |
|
78 | ||
79 | 4101192x |
Box intersection(const Box & other) const { |
80 |
return { |
|
81 | 4101192x |
std::max(xmin, other.xmin), |
82 | 4101192x |
std::max(ymin, other.ymin), |
83 | 4101192x |
std::min(xmax, other.xmax), |
84 | 4101192x |
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 | 517462x |
bool empty() const { |
118 | 517462x |
return xmin >= xmax || ymin >= ymax; |
119 |
} |
|
120 | ||
121 | 13x |
Box expand_to_include(const Box & other) const { |
122 | 13x |
if (empty()) { |
123 | ! |
return other; |
124 |
} |
|
125 | ||
126 | 13x |
if (other.empty()) { |
127 | ! |
return *this; |
128 |
} |
|
129 | ||
130 | 13x |
return { std::min(xmin, other.xmin), |
131 | 13x |
std::min(ymin, other.ymin), |
132 | 13x |
std::max(xmax, other.xmax), |
133 | 13x |
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 | 400883x |
bool operator==(const Box& other) const { |
143 | 400883x |
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) 2021 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_AREA_H |
|
15 |
#define EXACTEXTRACT_RASTER_AREA_H |
|
16 | ||
17 |
#include "raster.h" |
|
18 | ||
19 |
namespace exactextract { |
|
20 |
template<typename T> |
|
21 |
class CartesianAreaRaster : public AbstractRaster<T> { |
|
22 |
public: |
|
23 | 14x |
explicit CartesianAreaRaster(const Grid<bounded_extent> & ex) : |
24 |
AbstractRaster<T>(ex), |
|
25 | 14x |
m_area(grid_cell(AbstractRaster<T>::grid(), 0, 0).area()) |
26 |
{} |
|
27 | ||
28 | 7116x |
T operator()(size_t row, size_t column) const override { |
29 |
(void) row; |
|
30 |
(void) column; |
|
31 | 7116x |
return m_area; |
32 |
} |
|
33 |
private: |
|
34 |
double m_area; |
|
35 |
}; |
|
36 | ||
37 |
template<typename T> |
|
38 |
class SphericalAreaRaster : public AbstractRaster<T> { |
|
39 |
public: |
|
40 | 25x |
explicit SphericalAreaRaster(const Grid<bounded_extent> &ex) : |
41 |
AbstractRaster<T>(ex), |
|
42 | 25x |
m_areas(ex.rows()) |
43 |
{ |
|
44 | 25x |
const auto& g = AbstractRaster<T>::grid(); |
45 | ||
46 | 25x |
double dlon = g.dx(); |
47 | 25x |
double dlat = g.dy(); |
48 | ||
49 | 99x |
for (size_t i = 0; i < g.rows(); i++) { |
50 | 74x |
double y = g.y_for_row(i); |
51 | 74x |
double ymin = y - 0.5 * dlat; |
52 | 74x |
double ymax = y + 0.5 * dlat; |
53 | ||
54 | 74x |
m_areas[i] = EARTH_RADIUS_SQ * PI_180 * std::abs(std::sin(ymin * PI_180) - std::sin(ymax * PI_180)) * dlon; |
55 |
} |
|
56 |
} |
|
57 | ||
58 | 276x |
T operator()(size_t row, size_t column) const override { |
59 |
(void) column; |
|
60 | 276x |
return m_areas[row]; |
61 |
} |
|
62 | ||
63 |
private: |
|
64 |
static constexpr double EARTH_RADIUS = 6378137; |
|
65 |
static constexpr double EARTH_RADIUS_SQ = EARTH_RADIUS * EARTH_RADIUS; |
|
66 |
static constexpr double PI_180 = 3.141592653589793238462643383279502884197169399375105 / 180; |
|
67 | ||
68 |
std::vector<double> m_areas; |
|
69 |
}; |
|
70 |
} |
|
71 | ||
72 |
#endif //EXACTEXTRACT_RASTER_AREA_H |
1 |
// Copyright (c) 2018-2022 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: |