R/hotspot-workflow.R

In spatialrisk: Spatial Concentration and Radius-Based Risk Calculations

#' Prepare fixed-radius concentration hotspot analysis
#'
#' @description
#' `prepare_spatialrisk()`, `select_candidates()`, and `optimize_hotspot()`
#' expose the main steps used by \code{\link{concentration_hotspot}}. They are useful
#' when the intermediate search state needs to be inspected or when the same
#' prepared portfolio is used in more than one hotspot search strategy.
#'
#' @param data A data.frame containing point-level exposures. Must include
#'   longitude and latitude in EPSG:4326 and the value of interest. Coordinates
#'   are projected internally to `crs_metric`.
#' @param value A string giving the numeric column in `data` to aggregate
#'   within each radius.
#' @param radius Numeric. Radius of the circle in meters.
#' @param cell_size Numeric. Size of the raster cells used for the initial
#'   screening raster.
#' @param lon A string giving the longitude column in `data`.
#' @param lat A string giving the latitude column in `data`.
#' @param crs_metric Numeric. EPSG code for a projected CRS with meter units.
#'   The default `3035` is ETRS89 / LAEA Europe.
#' @param x A prepared spatial-risk workflow object returned by
#'   `prepare_spatialrisk()` or `select_candidates()`.
#' @param grid_spacing Numeric. Spacing between candidate grid centres in the
#'   units of `crs_metric`; for the default metric CRS these units are meters.
#'   Used for grid-based refinement.
#' @param max_refinement_points Positive integer. Maximum number of local points
#'   used for pair-intersection refinement before falling back to grid
#'   refinement in a search state produced by `select_candidates()`. Direct
#'   optimisation of a prepared object never falls back to a grid; values above
#'   this limit produce a computational-cost warning instead.
#' @param method Hotspot search strategy. `"continuous"` is the default and
#'   searches for centres that may lie between observed points. `"observed"`
#'   searches only observed point locations. `"grid"` uses the grid-refinement
#'   workflow.
#' @param threshold Optional numeric lower bound for candidate focal cells.
#'   If `NULL`, the lower bound is estimated using the same preliminary
#'   refinement step as `concentration_hotspot()`.
#' @param n_hotspots Positive integer. Number of non-overlapping hotspots to
#'   return.
#' @param top_n Deprecated. Use `n_hotspots` instead.
#' @param grid_precision Deprecated. Use `grid_spacing` instead.
#' @param progress Logical. Whether to print progress messages.
#' @param type Plot type. `"auto"` shows the prepared raster before candidate
#'   selection and selected focal candidate cells afterwards.
#' @param basemaps Character vector of basemap providers for `plot()`.
#'   The first provider is shown initially; the others are available in the
#'   layer control. Defaults to `"Esri.WorldGrayCanvas"`, `"OpenStreetMap"`,
#'   and `"Esri.WorldImagery"`. Passed to `mapview::mapview()`.
#' @param ... Additional arguments passed to `mapview::mapview()`.
#'
#' @return
#' `prepare_spatialrisk()` and `select_candidates()` return an object of class
#' `spatialrisk_hotspot_workflow`. `optimize_hotspot()` returns the same
#' `hotspot` object structure as \code{\link{concentration_hotspot}}.
#'
#' @details
#' The three-step interface decomposes the hotspot workflow without replacing
#' `concentration_hotspot()`. The wrapper remains the simplest public function
#' for normal use, while the decomposed functions make the intermediate
#' candidate selection visible.
#'
#' `optimize_hotspot()` optimises over the candidate search state supplied to
#' it. Called directly on the output of `prepare_spatialrisk()`, it performs a
#' full geometric search: candidate centres are the active observed locations
#' and the valid radius-circle intersections generated by every active point
#' pair no farther than twice the radius apart. This route does not use raster
#' screening or grid fallback and is intended mainly for small portfolios,
#' diagnostics, validation, and methodological benchmarks. Pair generation can
#' be computationally expensive.
#'
#' In `select_candidates()`, `threshold = NULL` estimates a lower bound by
#' taking the highest focal raster cells, refining those cells on a small local
#' grid, evaluating those centres in the projected metric coordinates, and using
#' the best feasible value as the candidate-cell threshold. The
#' focal cells first survive when their moving-window sum is at least this
#' lower bound. For continuous search with non-negative values and an automatic
#' threshold, a second bound sums only active points whose minimum distance to
#' the closed centre-cell rectangle is within the radius, allowing for the
#' scoring tolerance. It does not count the entire value of an extra raster
#' cell. The five highest surviving point bounds also seed feasible trial
#' centres; only their actual portfolio sums can raise the lower bound.
#' Cells with a point bound strictly below that lower bound are removed.
#' The focal window includes the requested radius plus a raster-cell
#' diagonal. With non-negative values, its sum is therefore an upper bound for
#' every exact centre located in that focal cell. Under the default automatic
#' threshold, observed and pair-intersection centres are mapped back to the
#' raster and exact evaluation is restricted to centres that lie in a selected
#' cell. For a point pair, its two possible circle centres are screened
#' separately. For a single hotspot, a streaming Rcpp angular sweep maintains
#' conservative interval totals over all active neighbours, including distance
#' slack. Competitive candidates are confirmed with the projected radius scorer;
#' this avoids materialising and separately querying every centre. With
#' complete pair-intersection refinement and exact
#' scoring, and provided no grid fallback occurs, this removes no cell that can
#' contain a strictly improving centre and therefore preserves a global optimum
#' of the single-disk problem. The guarantee does not apply to a
#' user-supplied threshold, negative values, or grid fallback.
#'
#' Candidate screening never defines the evaluation portfolio. Every retained
#' centre is scored against all records in the complete active portfolio, so a
#' point outside the candidate-generation subset still contributes when it lies
#' within the radius. Point-to-raster-cell membership is stored during
#' preparation and reused to retrieve nearby records; candidate regions within
#' one hotspot iteration share the same active-portfolio evaluation state. When
#' `optimize_hotspot(n_hotspots > 1)` or
#' `concentration_hotspot(n_hotspots > 1)` is used,
#' the points in the selected hotspot are removed and the candidate-selection
#' logic is run again for the next hotspot. Therefore the number of candidate
#' cells shown by `select_candidates()` for the first iteration does not limit
#' the number of hotspots returned by `n_hotspots`.
#'
#' @examples
#' portfolio <- Groningen[1:200, c("lon", "lat", "amount")]
#'
#' model <- prepare_spatialrisk(portfolio, value = "amount", radius = 200,
#'                              cell_size = 100)
#'
#' # Full geometric reference search
#' full <- optimize_hotspot(model, progress = FALSE)
#'
#' # Screened continuous search
#' screened <- model |>
#'   select_candidates(progress = FALSE) |>
#'   optimize_hotspot(progress = FALSE)
#'
#' full$hotspots
#' screened$hotspots
#'
#' @author Martin Haringa
#'
#' @export
prepare_spatialrisk <- function(data, value, radius = 200, cell_size = 100,
                                lon = "lon", lat = "lat",
                                crs_metric = 3035) {
  value <- validate_hotspot_value(value)
  data <- as.data.frame(data)
  check_hotspot_columns(data, value, lon, lat)
  check_hotspot_crs_metric(crs_metric)
  check_hotspot_radius(radius)
  check_hotspot_cell_size(cell_size)

  original <- data
  original$ix <- seq_len(nrow(original))
  # Metric coordinates are stored once so all later distance calculations use
  # Euclidean geometry in the chosen projected CRS.
  metric <- convert_crs_df(original, 4326, crs_metric, lon, lat, "x", "y")
  # The prepared object keeps both the rasterised portfolio and the metric
  # points, making the paper workflow reproducible step by step.
  state <- initialise_prepare_hotspot_state(original, value, radius, cell_size,
                                            lon, lat, crs_metric)

  out <- list(
    data = original,
    metric = metric,
    params = list(
      value = value,
      radius = radius,
      lon = lon,
      lat = lat,
      crs_metric = crs_metric,
      cell_size = cell_size
    ),
    candidates = NULL,
    state = state
  )
  class(out) <- c("spatialrisk_hotspot_workflow", class(out))
  out
}

#' @rdname prepare_spatialrisk
#' @export
select_candidates <- function(x, grid_spacing = 1,
                              max_refinement_points = 1500,
                              method = c("continuous", "observed", "grid"),
                              threshold = NULL, progress = TRUE,
                              grid_precision = lifecycle::deprecated()) {
  UseMethod("select_candidates")
}

#' @export
select_candidates.spatialrisk_hotspot_workflow <- function(
    x,
    grid_spacing = 1,
    max_refinement_points = 1500,
    method = c("continuous", "observed", "grid"),
    threshold = NULL,
    progress = TRUE,
    grid_precision = lifecycle::deprecated()
) {
  method <- match.arg(method)
  args <- resolve_hotspot_deprecated_args(
    n_hotspots = 1,
    top_n = lifecycle::deprecated(),
    top_n_supplied = FALSE,
    grid_spacing = grid_spacing,
    grid_precision = grid_precision,
    grid_precision_supplied = lifecycle::is_present(grid_precision),
    caller = "select_candidates"
  )
  grid_spacing <- args$grid_spacing
  p <- complete_hotspot_workflow_params(x$params)
  check_hotspot_threshold(threshold)
  validate_hotspot_search_settings(
    data = x$data,
    value = p$value,
    top_n = 1,
    radius = p$radius,
    cell_size = p$cell_size,
    grid_precision = grid_spacing,
    max_refinement_points = max_refinement_points,
    lon = p$lon,
    lat = p$lat,
    crs_metric = p$crs_metric,
    progress = progress
  )

  x$params$grid_spacing <- grid_spacing
  x$params$grid_precision <- grid_spacing
  x$params$max_refinement_points <- max_refinement_points
  x$params$method <- method
  x$params$threshold <- threshold

  if (method == "observed") {
    # Observed-point candidates form a fast deterministic baseline; they do not
    # guarantee the continuous optimum when the best centre lies between points.
    hotspot_progress(progress, "Selecting observed point candidates.")
    x$candidates <- data.frame(
      data_row = x$data$ix,
      x = x$metric$x,
      y = x$metric$y
    )
    x$state <- NULL
    return(x)
  }

  hotspot_progress(progress, "Selecting terra focal candidate area.")
  state <- x$state
  if (is.null(state)) {
    state <- initialise_prepare_hotspot_state(
      x$data, p$value, p$radius, p$cell_size, p$lon, p$lat, p$crs_metric
    )
  }
  if (is.null(state$focal)) {
    stop("Candidate selection requires a focal raster. Try a smaller ",
         "`cell_size`, a larger dataset extent, or `method = \"observed\"`.",
         call. = FALSE)
  }
  if (is.null(threshold)) {
    threshold <- estimate_hotspot_candidate_threshold(state$focal, x$data, p)
  }
  # Candidate cells are an interpretable inspection layer: their focal value is
  # high enough to warrant local refinement in the next optimisation step.
  candidate_cells <- cells_above_threshold_with_values(state$focal, threshold)
  if (method == "continuous" && is.null(x$params$threshold) &&
      all(x$metric[[p$value]] >= 0)) {
    screened <- tighten_hotspot_candidate_cells(
      candidate_cells, state, x$metric, p$value, p$radius, threshold
    )
    candidate_cells <- screened$cells
    threshold <- screened$threshold
  }
  approximate <- terra_screening_center(state$focal)
  approximate_cell <- terra::cellFromXY(
    state$raster,
    matrix(c(approximate$x[1], approximate$y[1]), ncol = 2)
  )
  local_ix <- local_pair_refine_subset(
    x$metric, approximate$x[1], approximate$y[1], p$radius, p$cell_size,
    state = state, cell = approximate_cell
  )

  x$state <- state
  x$candidates <- list(
    approximate_center = approximate,
    threshold = threshold,
    cells = candidate_cells,
    local_data_rows = x$metric$ix[local_ix],
    local_points = length(local_ix),
    refinement = if (method == "continuous" &&
                     length(local_ix) <= max_refinement_points) {
      "pair_intersections"
    } else if (method == "continuous") {
      "grid"
    } else {
      "grid"
    }
  )
  x
}

#' @rdname prepare_spatialrisk
#' @export
optimize_hotspot <- function(x, n_hotspots = 1, progress = TRUE,
                             top_n = lifecycle::deprecated()) {
  UseMethod("optimize_hotspot")
}

#' @export
optimize_hotspot.spatialrisk_hotspot_workflow <- function(
    x,
    n_hotspots = 1,
    progress = TRUE,
    top_n = lifecycle::deprecated()
) {
  args <- resolve_hotspot_deprecated_args(
    n_hotspots = n_hotspots,
    top_n = top_n,
    top_n_supplied = lifecycle::is_present(top_n),
    grid_spacing = 1,
    grid_precision = lifecycle::deprecated(),
    grid_precision_supplied = FALSE,
    caller = "optimize_hotspot"
  )
  n_hotspots <- args$n_hotspots
  p <- complete_hotspot_workflow_params(x$params)
  validate_hotspot_search_settings(
    data = x$data,
    value = p$value,
    top_n = n_hotspots,
    radius = p$radius,
    cell_size = p$cell_size,
    grid_precision = p$grid_spacing,
    max_refinement_points = p$max_refinement_points,
    lon = p$lon,
    lat = p$lat,
    crs_metric = p$crs_metric,
    progress = progress
  )

  data <- x$data
  data$ix <- NULL
  full_search <- is.null(x$candidates)

  if (full_search) {
    if (p$method != "continuous") {
      stop("Direct optimisation of a prepared object is defined for ",
           "`method = \"continuous\"`. Run `select_candidates()` to choose ",
           "another search method.", call. = FALSE)
    }
    if (any(x$data[[p$value]] < 0)) {
      stop("Full geometric optimisation requires non-negative `value` ",
           "weights.", call. = FALSE)
    }
    if (nrow(x$data) > p$max_refinement_points) {
      warning(
        "Full geometric optimisation will generate pair candidates from ",
        nrow(x$data), " active points. This can be expensive; ",
        "`max_refinement_points` does not trigger grid fallback for a ",
        "prepared object.",
        call. = FALSE
      )
    }
    hotspot_progress(
      progress,
      "Using full geometric candidate search over the active portfolio."
    )
    return(concentration_hotspot_pair_intersections(
      data = data,
      value = p$value,
      top_n = n_hotspots,
      radius = p$radius,
      lon = p$lon,
      lat = p$lat,
      crs_metric = p$crs_metric,
      progress = progress
    ))
  }

  if (p$method == "observed") {
    hotspot_progress(progress, "Using observed-points hotspot search.")
    out <- concentration_hotspot_indexed(
      data = data,
      value = p$value,
      top_n = n_hotspots,
      radius = p$radius,
      lon = p$lon,
      lat = p$lat,
      crs_metric = p$crs_metric,
      print_progress = progress,
      cell_size = p$radius,
      candidate_rows = x$candidates$data_row
    )
    attr(out, "candidate_search") <- "screened"
    return(out)
  }

  if (p$method == "grid") {
    hotspot_progress(progress, "Using grid-refinement hotspot search.")
    out <- concentration_hotspot_terra(
      data = data,
      value = p$value,
      top_n = n_hotspots,
      radius = p$radius,
      cell_size = p$cell_size,
      grid_precision = p$grid_spacing,
      lon = p$lon,
      lat = p$lat,
      crs_metric = p$crs_metric,
      progress = progress
    )
    attr(out, "candidate_search") <- "screened"
    return(out)
  }

  out <- concentration_hotspot_pair_refine(
    data = data,
    value = p$value,
    top_n = n_hotspots,
    radius = p$radius,
    cell_size = p$cell_size,
    grid_precision = p$grid_spacing,
    max_refinement_points = p$max_refinement_points,
    threshold = p$threshold,
    lon = p$lon,
    lat = p$lat,
    crs_metric = p$crs_metric,
    progress = progress,
    initial_candidate_cells = x$candidates$cells,
    initial_state = x$state
  )
  attr(out, "candidate_search") <- "screened"
  out
}

#' @export
print.spatialrisk_hotspot_workflow <- function(x, ...) {
  p <- complete_hotspot_workflow_params(x$params)
  cat("<spatialrisk_hotspot_workflow>\n")
  cat("Rows:", nrow(x$data), "\n")
  cat("Value:", p$value, "\n")
  cat("Radius:", p$radius, "meters\n")
  cat("Method:", p$method, "\n")
  cat("Candidate search:",
      if (is.null(x$candidates)) "full active portfolio" else "selected state",
      "\n")
  if (!is.null(x$candidates)) {
    if (is.list(x$candidates)) {
      cat("Candidate cells:", nrow(x$candidates$cells), "\n")
      cat("Candidate threshold:", x$candidates$threshold, "\n")
      cat("Candidate local points:", x$candidates$local_points, "\n")
      cat("Refinement:", x$candidates$refinement, "\n")
    } else {
      cat("Candidate centres:", nrow(x$candidates), "\n")
    }
  }
  invisible(x)
}

#' @rdname prepare_spatialrisk
#' @method plot spatialrisk_hotspot_workflow
#' @export
plot.spatialrisk_hotspot_workflow <- function(x, type = c("auto", "raster",
                                                          "candidates"),
                                              ...,
                                              basemaps = c("Esri.WorldGrayCanvas",
                                                           "OpenStreetMap",
                                                           "Esri.WorldImagery")) {
  type <- match.arg(type)
  if (type == "auto") {
    type <- if (is.null(x$candidates)) "raster" else "candidates"
  }

  if (type == "candidates" && is.null(x$candidates)) {
    stop("Run `select_candidates()` before plotting candidate cells.",
         call. = FALSE)
  }

  if (!requireNamespace("mapview", quietly = TRUE)) {
    stop("mapview is needed for this function to work. Install it via ",
         "install.packages(\"mapview\")", call. = FALSE)
  }
  p <- complete_hotspot_workflow_params(x$params)

  if (type == "raster") {
    if (is.null(x$state) || is.null(x$state$rasterized)) {
      stop("`x` does not contain a prepared raster.", call. = FALSE)
    }
    rasterized <- x$state$rasterized
    names(rasterized) <- paste0(p$value, "_sum")
    return(mapview::mapview(
      rasterized,
      layer.name = paste0(p$value, " sum per raster cell"),
      map.types = basemaps,
      ...
    ))
  }

  if (is.data.frame(x$candidates)) {
    candidates <- convert_crs_df(x$candidates, p$crs_metric, 4326,
                                 lon_from = "x", lat_from = "y",
                                 lon_to = p$lon, lat_to = p$lat)
    candidates <- convert_df_to_sf(candidates, p$lon, p$lat, 4326,
                                   p$crs_metric)
    return(mapview::mapview(candidates, layer.name = "Observed candidates",
                            map.types = basemaps, ...))
  }

  candidate_polygons <- candidate_cells_polygons(x$state$focal,
                                                 x$candidates$cells)
  mapview::mapview(candidate_polygons,
                   zcol = "focal_value",
                   layer.name = "Candidate focal cells",
                   map.types = basemaps,
                   ...)
}

check_hotspot_radius <- function(radius) {
  if (!is.numeric(radius) || length(radius) != 1L ||
      is.na(radius) || !is.finite(radius) || radius <= 0) {
    stop("`radius` must be a single finite positive number.", call. = FALSE)
  }
  invisible(NULL)
}

check_hotspot_cell_size <- function(cell_size) {
  if (!is.numeric(cell_size) || length(cell_size) != 1L ||
      is.na(cell_size) || !is.finite(cell_size) || cell_size <= 0) {
    stop("`cell_size` must be a single finite positive number.",
         call. = FALSE)
  }
  invisible(NULL)
}

check_hotspot_threshold <- function(threshold) {
  if (!is.null(threshold) &&
      (!is.numeric(threshold) || length(threshold) != 1L ||
       is.na(threshold) || !is.finite(threshold))) {
    stop("`threshold` must be NULL or a single finite numeric value.",
         call. = FALSE)
  }
  invisible(NULL)
}

validate_hotspot_search_settings <- function(data, value, top_n, radius,
                                             cell_size, grid_precision,
                                             max_refinement_points, lon, lat,
                                             crs_metric, progress) {
  validate_pair_refine_input(data, value, top_n, radius, cell_size,
                             grid_precision, max_refinement_points, lon, lat,
                             crs_metric, progress)
}

complete_hotspot_workflow_params <- function(params) {
  defaults <- list(
    cell_size = 100,
    grid_spacing = 1,
    max_refinement_points = 1500,
    method = "continuous",
    threshold = NULL
  )
  out <- utils::modifyList(defaults, params)
  if (is.null(params$grid_spacing) && !is.null(params$grid_precision)) {
    out$grid_spacing <- params$grid_precision
  }
  out$grid_precision <- out$grid_spacing
  out
}

HOTSPOT_LOWER_BOUND_SEED_CELLS <- 5L
HOTSPOT_LOWER_BOUND_TRIALS_PER_AXIS <- 10L

estimate_hotspot_candidate_threshold <- function(focal, data, params) {
  top_focals <- top_n_focals(focal, n = HOTSPOT_LOWER_BOUND_SEED_CELLS)
  metric_hotspot_lower_bound(top_focals, data, params)$concentration[1]
}

metric_hotspot_lower_bound <- function(candidate_cells, data, params,
                                       points =
                                         HOTSPOT_LOWER_BOUND_TRIALS_PER_AXIS,
                                       metric = NULL) {
  offsets <- seq(-params$cell_size / 2, params$cell_size / 2,
                 length.out = points)
  grid <- expand.grid(dx = offsets, dy = offsets)
  x_candidates <- rep(candidate_cells$x, each = nrow(grid)) +
    rep(grid$dx, times = nrow(candidate_cells))
  y_candidates <- rep(candidate_cells$y, each = nrow(grid)) +
    rep(grid$dy, times = nrow(candidate_cells))

  if (is.null(metric)) metric <- convert_crs_df(
    data,
    crs_from = 4326,
    crs_to = params$crs_metric,
    lon_from = params$lon,
    lat_from = params$lat,
    lon_to = "x",
    lat_to = "y"
  )
  ix <- if ("ix" %in% names(metric)) metric$ix else seq_len(nrow(metric))

  # The screening threshold must be a feasible objective value in exactly the
  # same projected Euclidean geometry used by pair-intersection refinement.
  indexed_concentration_best_cpp(
    x_candidates = x_candidates,
    y_candidates = y_candidates,
    x_ref = metric$x,
    y_ref = metric$y,
    value_ref = metric[[params$value]],
    ix_ref = as.integer(ix),
    radius = params$radius,
    cell_width = params$radius
  )
}

initialise_prepare_hotspot_state <- function(data, value, radius, cell_size,
                                             lon, lat, crs_metric) {
  tryCatch(
    initialise_terra_hotspot_state(data, value, radius, cell_size,
                                   lon, lat, crs_metric),
    error = function(e) {
      metric_sf <- convert_crs_df(data, 4326, crs_metric, lon, lat, "x", "y")
      terra_crs <- paste0("EPSG:", crs_metric)
      spatvctr <- terra::vect(metric_sf, geom = c("x", "y"), crs = terra_crs)
      raster <- hotspot_raster_template(spatvctr, cell_size, terra_crs, radius)
      rasterized <- terra::rasterize(spatvctr, raster, field = value, fun = sum)
      point_cells <- terra::cellFromXY(
        raster,
        as.matrix(metric_sf[, c("x", "y"), drop = FALSE])
      )
      list(
        spatvctr = spatvctr,
        raster = raster,
        rasterized = rasterized,
        moving_window = NULL,
        focal = NULL,
        point_cells = data.frame(ix = metric_sf$ix, cell = point_cells),
        points_by_cell = split(metric_sf$ix, as.character(point_cells)),
        raster_geometry = hotspot_raster_geometry(raster)
      )
    }
  )
}

cells_above_threshold_with_values <- function(focal, threshold) {
  cells <- cells_above_threshold(focal, threshold)
  if (is.null(cells) || nrow(cells) == 0L) {
    return(data.frame(cell = integer(), x = numeric(), y = numeric(),
                      focal_value = numeric()))
  }
  values <- terra::values(focal, mat = FALSE)
  cells$focal_value <- values[cells$cell]
  cells
}

tighten_hotspot_candidate_cells <- function(cells, state, metric, value,
                                            radius, threshold) {
  if (!nrow(cells)) return(list(cells = cells, threshold = threshold))
  geometry <- state$raster_geometry
  # The mapping survives greedy removal; pass only the currently active rows.
  point_cells <- state$point_cells$cell[
    match(metric$ix, state$point_cells$ix)
  ]
  cells$point_upper_bound <- cell_point_upper_bounds_cpp(
    as.integer(cells$cell), cells$x, cells$y, as.integer(point_cells),
    metric$x, metric$y, metric[[value]], radius,
    geometry$xres, geometry$yres, geometry$nrow, geometry$ncol
  )
  cells <- cells[cells$point_upper_bound >= threshold, , drop = FALSE]
  if (!nrow(cells)) return(list(cells = cells, threshold = threshold))
  seeds <- head(order(-cells$point_upper_bound, cells$cell),
                HOTSPOT_LOWER_BOUND_SEED_CELLS)
  # Bounds rank seed regions; only an actual full-portfolio disk value can
  # raise the lower bound. The bound itself is never used as a hotspot value.
  feasible <- metric_hotspot_lower_bound(
    cells[seeds, , drop = FALSE], data = NULL,
    params = list(cell_size = min(geometry$xres, geometry$yres),
                  radius = radius, value = value), metric = metric
  )$concentration[1]
  threshold <- max(threshold, feasible)
  cells <- cells[cells$point_upper_bound >= threshold, , drop = FALSE]
  list(cells = cells, threshold = threshold)
}

candidate_cells_polygons <- function(focal, candidate_cells) {
  out <- focal
  vals <- rep(NA_real_, length(terra::values(focal, mat = FALSE)))
  vals[candidate_cells$cell] <- candidate_cells$focal_value
  names(out) <- "focal_value"
  terra::values(out) <- vals
  terra::as.polygons(out, na.rm = TRUE)
}