This function reduces the resolution of a raster by multiple aggregation methods and then quantifies the relative error induced by each method from the differences between the original values and the aggregated values. To implement the function, a raster (x) must be supplied as well as an aggregation factor (fact) and a named list of functions (stat) used to aggregate the raster. The raster is aggregated using each method (function) and mapped back onto the original resolution for calculation of the differences between the original raster and the aggregated raster(s). The function returns a visual statistical summary of the differences (if plot = TRUE) and a named list comprising the aggregated raster(s) and the re-sampled version(s) of those mapped back onto the original resolution.

process_surface(
  x,
  fact = 2L,
  stat = list(mean = mean),
  ...,
  plot = TRUE,
  cl = NULL,
  varlist = NULL,
  verbose = TRUE
)

Arguments

x

A raster to be processed. For implementations preceding a call to of one of flapper's particle filtering algorithms, x should be planar (i.e., Universal Transverse Mercator projection) with equal resolution in the x, y directions and identical units in the x, y and z directions (see pf).

fact

A positive integer that defines by how much x should be aggregated (see aggregate).

stat

A named list of functions used to aggregate x (see the fun argument of aggregate).

...

Additional arguments passed to aggregate to control aggregation.

plot

A logical input that defines whether or not to plot a summary of the differences between the original raster (x) and the aggregated raster(s). If specified, the minimum, median and maximum difference are shown for each statistic (stat).

cl, varlist

(optional) Parallelisation options. cl is (a) a cluster object from makeCluster or (b) an integer that defines the number of child processes. varlist is a character vector of variables for export (see cl_export). Exported variables must be located in the global environment. If a cluster is supplied, the connection to the cluster is closed within the function (see cl_stop). For further information, see cl_lapply and flapper-tips-parallel.

verbose

A logical input that defines whether or not to print messages to the console to relay function progress.

Value

The function returns a plot of the differences between the original and aggregated raster(s), if plot = TRUE, and a named list of (a) the aggregated raster(s) (`agg_by_stat'), (b) the aggregated, resampled raster(s) (`agg_by_stat_rs') and (c) the summary statistics plotted.

Details

This function was motivated by the particle filtering algorithms in flapper (e.g., pf). For these algorithms, it is computationally beneficial to reduce raster resolution, where possible, by aggregation. To facilitate this process, this function quantifies the relative error induced by different aggregation functions. If appropriate, the particle filtering algorithm(s) can then be implemented using the aggregated raster that minimises the error, with the magnitude of that error propagated via the depth_error parameter.

See also

Author

Edward Lavender

Examples

# Define the raster for which to implement the function
x <- dat_gebco
blank <- raster::raster(raster::extent(x), crs = raster::crs(x), resolution = 250)
x <- raster::resample(x, blank, method = "bilinear")
#> Warning: aggregation factor is larger than the number of columns
#> Warning: aggregation factor is larger than the number of rows
# Implement function using a list of statistics
out <- process_surface(x, fact = 2, stat = list(min = min, mean = mean, median = median, max = max))
#> flapper::process_surface() called (@ 2023-08-29 15:46:08)... 
#> ... Aggregating raster... 
#> Warning: aggregation factor is larger than the number of columns
#> Warning: aggregation factor is larger than the number of rows
#> Warning: aggregation factor is larger than the number of columns
#> Warning: aggregation factor is larger than the number of rows
#> Warning: aggregation factor is larger than the number of columns
#> Warning: aggregation factor is larger than the number of rows
#> Warning: aggregation factor is larger than the number of columns
#> Warning: aggregation factor is larger than the number of rows
#> ... Resampling aggregated raster(s) back onto the original resolution... 
#> ... Computing differences between the original and aggregated raster(s)... 
#> ... Summarising the differences between rasters across statistic(s)... 
#> Warning: no non-missing arguments to min; returning Inf
#> Warning: no non-missing arguments to min; returning Inf
#> Warning: no non-missing arguments to min; returning Inf
#> Warning: no non-missing arguments to min; returning Inf
#> Warning: no non-missing arguments to max; returning -Inf
#> Warning: no non-missing arguments to max; returning -Inf
#> Warning: no non-missing arguments to max; returning -Inf
#> Warning: no non-missing arguments to max; returning -Inf
#> 4 observation pair(s) in x are NA; these are removed.
#> Error in (function (side = 1:4, x = list(), lim = list(), pretty = list(n = 5),     units = list(), axis = list(), pi_notation = NULL, control_axis = list(las = TRUE),     control_sci_notation = list(), control_digits = NULL, control_factor_lim = 0.5,     axis_ls = NULL, add = FALSE, return_list = TRUE, ...) {    if (is.null(axis_ls)) {        if (length(side) == 3)             stop("Three axis sides are not currently supported.")        if (length(side) > 1) {            if (!(side[1] %in% c(1, 3)))                 stop("side[1] should refer to an x axis (i.e., side[1] = 1 or 3).")            if (!(side[2] %in% c(2, 4)))                 stop("side[2] should refer to a y axis (i.e., side[2] = 2 or 4).")        }        if (length(x) == 0) {            if (length(lim) == 0) {                stop("'x' or 'lim' must be supplied.")            }            else {                message("pretty_axis() 'x' values taken from 'lim'.")                x <- lim            }        }        check_input_class(arg = "x", input = x, if_class = NULL,             to_class = "list", type = "stop")        lim <- empty_list_to_list_null("lim", lim)        pretty <- empty_list_to_list_null("pretty", pretty)        units <- empty_list_to_list_null("units", units)        axis <- empty_list_to_list_null("axis", axis)        if (is.null(pi_notation))             pi_notation <- lapply(1:length(side), function(x) -1)        pi_notation <- empty_list_to_list_null("pi_notation",             pi_notation)        mapply(list(x, lim, units), list("x", "lim", "units"),             FUN = function(elm, elm_name) {                if (length(elm) > 0 & length(elm) > length(side)) {                  stop(paste0(length(side), " side(s) supplied but argument '",                     elm_name, "' contains ", length(elm), " elements."))                }            })        x <- mapply(x, 1:length(x), FUN = function(elm, i) {            arg <- paste0("x[[", i, "]]")            elm <- check_input_class(arg = arg, input = elm,                 if_class = "character", to_class = "factor",                 type = "warning", coerce_input = factor)            elm <- check_tz(arg, elm)            return(elm)        }, SIMPLIFY = FALSE)        lim <- mapply(lim, 1:length(lim), FUN = function(elm,             i) {            arg <- paste0("lim[[", i, "]]")            elm <- check_tz(arg, elm)            return(elm)        }, SIMPLIFY = FALSE)        lx <- sapply(x, length)        if (length(unique(lx)) != 1) {            message("'x' contains elements with different numbers of observations; collapsing each element (i) in 'x' to range(i).")            x <- lapply(x, function(e) {                if (inherits(e, "factor")) {                  rng <- range_factor(e)                }                else {                  rng <- range(e, na.rm = TRUE)                }                return(rng)            })            lx <- c(2, 2)        }        dat <- data.frame(do.call(cbind, x))        dat <- dat[stats::complete.cases(dat), , drop = FALSE]        nrw <- nrow(dat)        if (nrw != lx[1]) {            lna <- lx[1] - nrw            message(paste(lna, "observation pair(s) in x are NA; these are removed."))            x <- lapply(1:ncol(dat), function(i) return(dat[,                 i]))            if (length(x[[1]]) < 1)                 stop("No non NA observations left in x.")        }        lim <- list_adjust(l = lim, f = length, side = side)        pretty <- list_adjust(l = pretty, f = list_depth, side = side)        units <- list_adjust(l = units, f = length, side = side)        axis <- list_adjust(l = axis, f = list_depth, side = side)        pi_notation <- list_adjust(l = pi_notation, f = list_depth,             side = side)        dots <- list(...)        dnms <- names(dots)        if (!is.null(dnms)) {            for (param in c("cex.axis", "cex.lab", "col.axis",                 "col.lab", "font.axis", "font.lab", "las")) {                if (param %in% dnms)                   control_axis[[param]] <- dots[[param]]            }        }        axis_ls <- mapply(side, x, lim, pretty, units, axis,             pi_notation, SIMPLIFY = FALSE, FUN = function(iside,                 ix, ilim, ipretty, iunits, iaxis, ipi) {                testing <- FALSE                if (testing) {                  iside <- side[1]                  ix <- x[[1]]                  ilim <- lim[[1]]                  ipretty <- pretty[[1]]                  iunits <- units[[1]]                  iaxis <- axis[[1]]                  ipi <- pi_notation[[1]]                }                if (inherits(ilim, "list"))                   ilim <- Reduce(c, ilim)                ipretty <- compact(ipretty)                iunits <- unlist(compact(iunits))                iaxis <- compact(iaxis)                ipi <- compact(ipi)                ipi <- compact(ipi)                if (all(sapply(1:length(ipi), function(i) isTRUE(unlist(ipi)[i] ==                   -1))))                   ipi <- NULL                if (is_number(ix)) {                  if (!is.null(ipi)) {                    ix <- ix/pi                    if (!is.null(ilim))                       ilim <- ilim/pi                  }                }                else {                  if (!is.null(ipi))                     warning("pi_notation argument(s) ignored for non-numeric variable(s).")                }                iunits_warn <- ifelse(length(iunits) != 0, TRUE,                   FALSE)                iunits <- units_x(iunits, x = ix)                ilim <- define_lim_init(x = ix, at = iaxis$at,                   lim = ilim)                if (is.null(iaxis$at)) {                  if (length(ipretty) > 0) {                    if (iunits_warn)                       warning("Both pretty and units specified for an axis. pretty arguments implemented.")                    pretty_seq_ls <- pretty_seq(x = ix, lim = ilim,                       pretty_args = ipretty)                    iaxis$at <- pretty_seq_ls$at                    ilim <- pretty_seq_ls$lim                  }                  else {                    iaxis$at <- seq_x(ilim[1], ilim[2], units = iunits)                  }                }                if (is.factor(ix)) {                  ilim <- c(min(ilim) - control_factor_lim, max(ilim) +                     control_factor_lim)                  attributes(ilim)$user <- c(FALSE, FALSE)                }                if (is.null(iaxis$labels))                   iaxis$labels <- pretty_labels(x = ix, at = iaxis$at,                     n = control_digits, pi_notation_args = ipi,                     sci_notation_args = control_sci_notation)                if (is_number(ix) & !is.null(ipi)) {                  ilim <- ilim * pi                  iaxis$at <- iaxis$at * pi                }                iaxis$side <- iside                if (length(control_axis) > 0)                   iaxis <- rlist::list.merge(iaxis, control_axis)                out <- list(axis = iaxis, lim = ilim)                return(out)            })        names(axis_ls) <- side        s1s <- c("1", "2", "3", "4")        s2s <- list(c("2", "4"), c("1", "3"), c("2", "4"), c("1",             "3"))        sfs <- list(min, min, max, max)        for (s in 1:4) {            s1 <- s1s[s]            s2 <- s2s[[s]]            s21 <- s2[1]            s22 <- s2[2]            sf <- sfs[[s]]            if (s1 %in% side) {                if (is.null(axis_ls[[s1]]$axis$pos) & !is.null(axis_ls[[s21]]$axis) |                   !is.null(axis_ls[[s22]]$axis)) {                  s2select <- s2[which(c(!is.null(axis_ls[[s21]]$axis),                     !is.null(axis_ls[[s22]]$axis)))]                  s2select <- s2select[1]                  axis_ls[[s1]]$axis$pos <- sf(axis_ls[[s2select]]$lim)                }            }        }    }    if (any(is.na(names(axis_ls)))) {        stop("names(axis_ls) contains NA(s). The number of sides and the number of elements in argument lists is not aligned.")    }    if (add) {        lapply(axis_ls, function(elem) {            tmp_axis_args <- elem$axis            tmp_axis_args$at <- elem$lim            tmp_axis_args$lwd.ticks <- c(0, 0)            tmp_axis_args$labels <- c("", "")            at <- elem$axis$at            add_axis <- choose_foo_axis(at)            do.call(add_axis, tmp_axis_args)            elem$axis$col.ticks <- elem$axis$col            elem$axis$col <- NA            if (is.null(elem$axis$col.ticks))                 elem$axis$col.ticks <- "black"            do.call(add_axis, elem$axis)        })    }    if (return_list)         return(axis_ls)})(add = FALSE, type = "n", side = 1:2, pretty = list(n = 5),     x = list(structure(1:4, levels = c("min", "mean", "median",     "max"), class = "factor"), c(min = NaN, mean = NaN, median = NaN,     max = NaN)), lim = list(x = NULL, y = c(-Inf, Inf))): No non NA observations left in x.
summary(out)
#> Error in summary(out): object 'out' not found