Package {spatialkit}

spatialkit: Spatial Tessellation, Modeling, and Cross-Validation Toolkit

Description

A modular toolkit for spatial analysis workflows including coordinate reference system management, Voronoi/Delaunay/grid tessellation, feature-to-polygon assignment, geographically weighted regression (GWR, via 'GWmodel'), Bayesian spatial Gaussian process regression (via 'brms'), spatial cross-validation with block and buffered strategies, and model comparison. Provides an S3 class system ('spatial_fit') with consistent predict, fitted, and residuals methods across model backends.

Author(s)

Maintainer: Justin Chase jchase.msu@gmail.com [copyright holder]

Authors:

• Justin Chase jchase.msu@gmail.com [copyright holder]

See Also

Useful links:

• https://github.com/elkronos/gis_modeling_toolkit

• Report bugs at https://github.com/elkronos/gis_modeling_toolkit/issues

Fast center points of per-feature bounding boxes

Description

Fast center points of per-feature bounding boxes

Usage

.bbox_center_sfc(x)

Arguments

Value

An sfc (POINT) vector.

Assign features to polygons and attach a polygon ID

Description

Joins an sf layer of input features to a polygon layer via spatial join.

Usage

assign_features_to_polygons(
  features_sf,
  polygons_sf,
  polygon_id_col = "poly_id",
  keep_unassigned = FALSE,
  predicate = sf::st_intersects,
  largest = TRUE,
  tie_break = c("smallest_area", "first")
)

Arguments

Value

An sf object with polygon_id_col attached.

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = runif(20, 0, 100), y = runif(20, 0, 100), val = rnorm(20)),
  coords = c("x", "y"), crs = 32632
)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
grid <- create_grid_polygons(bnd, target_cells = 9, type = "square")
assigned <- assign_features_to_polygons(pts, grid)
table(assigned$poly_id)

Build a tessellation (Voronoi, Delaunay triangles, hex grid, or square grid)

Description

Build a tessellation (Voronoi, Delaunay triangles, hex grid, or square grid)

Usage

build_tessellation(
  points_sf,
  boundary = NULL,
  method = c("voronoi", "triangles", "hex", "square"),
  approx_n_cells = NULL,
  cellsize = NULL,
  expand = 0,
  clip = TRUE,
  keep_duplicates = FALSE,
  crs = NULL,
  quiet = FALSE
)

Arguments

Value

A list with cells, index, boundary, method, params.

See Also

Other tessellation: create_grid_polygons(), create_voronoi_polygons(), get_voronoi_seeds(), plot_tessellation_map()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = runif(20, 0, 100), y = runif(20, 0, 100)),
  coords = c("x", "y"), crs = 32632
)
tess <- build_tessellation(pts, method = "voronoi", quiet = TRUE)
tess$cells

Clear cached fitted values for a Bayesian spatial model

Description

Removes the lazily-cached fitted() result so that the next call recomputes from the posterior. This is only necessary if the underlying brmsfit engine or training data has been manually mutated after fitting — normal usage never requires it.

Usage

clear_fitted_cache(object)

Arguments

Value

object, invisibly (called for side effect).

Clear the in-session grid cache

Description

Removes all memoized grid results from the internal cache environment.

Usage

clear_grid_cache(cache_env = .gmt_cache)

Arguments

Value

Invisibly, the number of entries removed.

Build a polygonal clip target from points and/or a boundary

Description

Build a polygonal clip target from points and/or a boundary

Usage

clip_target_for(points_sf, boundary = NULL, expand = 0, quiet = FALSE)

Arguments

Value

An sf polygon layer representing the clip target.

Extract Bayesian model fixed-effect summaries

Description

Extract Bayesian model fixed-effect summaries

Usage

## S3 method for class 'bayesian_fit'
coef(object, ...)

Arguments

Value

A data.frame of fixed-effect posterior summaries.

Extract GWR local coefficients

Description

Extract GWR local coefficients

Usage

## S3 method for class 'gwr_fit'
coef(object, ...)

Arguments

Value

A data.frame of local coefficient estimates (one row per obs).

Coerce arbitrary geometries to representative points

Description

Converts the geometry column of an sf object to POINTs using one of several strategies.

Usage

coerce_to_points(
  x,
  mode = c("auto", "centroid", "point_on_surface", "surface", "line_midpoint",
    "bbox_center"),
  tmp_project = TRUE
)

Arguments

Value

An sf object with geometry coerced to POINTs.

Examples

library(sf)
poly <- st_sf(
  id = 1,
  geometry = st_sfc(st_polygon(list(rbind(
    c(0, 0), c(2, 0), c(2, 2), c(0, 2), c(0, 0)
  ))), crs = 32632)
)
coerce_to_points(poly, "auto")  # interior representative point

Side-by-side comparison of fitted spatial models

Description

Takes a named list of already-fit spatial_fit objects and produces a tidy comparison table including in-sample metrics and model-specific information criteria (AICc, LOOIC).

Usage

compare_models(fits, newdata = NULL, ...)

Arguments

Value

A data.frame comparing all models.

See Also

Other model evaluation: compare_models_cv(), evaluate_insample(), residual_morans_i()

Cross-validated comparison of spatial models

Description

Fits and cross-validates one or more model types, returning a unified comparison table. Unlike compare_models(), this function does perform fitting (inside CV folds), because CV inherently requires repeated fitting.

Usage

compare_models_cv(
  data_sf,
  response_var,
  predictor_vars,
  models = c("GWR", "Bayesian"),
  k = 5,
  seed = 123,
  folds = NULL,
  boundary = NULL,
  pointize = "auto",
  gwr_args = list(),
  bayes_args = list(),
  summary = c("mean", "median"),
  quiet = FALSE
)

Arguments

Value

A list with overall, by_fold, and per-model cv_results.

See Also

Other model evaluation: compare_models(), evaluate_insample(), residual_morans_i()

Create square or hexagonal grid polygons over a boundary

Description

Create square or hexagonal grid polygons over a boundary

Usage

create_grid_polygons(
  boundary,
  target_cells = NULL,
  type = c("square", "hex"),
  cellsize = NULL,
  n = NULL,
  clip = TRUE,
  crs = NULL,
  quiet = FALSE
)

Arguments

Value

An sf polygon layer with poly_id column.

See Also

Other tessellation: build_tessellation(), create_voronoi_polygons(), get_voronoi_seeds(), plot_tessellation_map()

Examples

library(sf)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
grid_sq  <- create_grid_polygons(bnd, target_cells = 25, type = "square")
grid_hex <- create_grid_polygons(bnd, target_cells = 25, type = "hex")
nrow(grid_sq)
nrow(grid_hex)  # near-target thanks to hex packing correction

Create and cache grid polygons over a boundary

Description

Builds a grid via create_grid_polygons() and memoizes the result so repeated calls with the same inputs return instantly.

Usage

create_grid_polygons_cached(
  boundary,
  target_cells,
  type = c("hex", "square"),
  ...,
  cache_env = .gmt_cache
)

Arguments

Value

An sf data frame with a stable poly_id column.

Create Voronoi polygons from points with robust CRS and optional clipping

Description

Create Voronoi polygons from points with robust CRS and optional clipping

Usage

create_voronoi_polygons(
  points_sf,
  boundary = NULL,
  expand = 0,
  clip = TRUE,
  keep_duplicates = FALSE,
  crs = NULL,
  quiet = FALSE
)

Arguments

Value

A list with cells, index, boundary, method, params.

See Also

Other tessellation: build_tessellation(), create_grid_polygons(), get_voronoi_seeds(), plot_tessellation_map()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = runif(15, 0, 100), y = runif(15, 0, 100)),
  coords = c("x", "y"), crs = 32632
)
res <- create_voronoi_polygons(pts, quiet = TRUE)
res$cells   # one polygon per unique point, with stable cell_id
res$index   # cell_id assignment for each input point

K-fold cross-validation for the Bayesian spatial model

Description

K-fold cross-validation for the Bayesian spatial model

Usage

cv_bayes(
  data_sf,
  response_var,
  predictor_vars,
  folds = NULL,
  k = 5,
  seed = 123,
  boundary = NULL,
  pointize = "auto",
  fit_args = list(),
  summary = c("mean", "median"),
  compute_pred_intervals = TRUE,
  coverage_levels = c(0.5, 0.8, 0.95),
  block_size = NULL,
  auto_range = FALSE,
  parallel = FALSE
)

Arguments

Value

A list with overall, fold_metrics, predictions, folds, formula, and predictive_coverage.

See Also

Other cross-validation: cv_gwr(), cv_spatial(), estimate_sac_range(), make_folds()

Examples

## Not run: 
if (requireNamespace("brms", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = runif(n, 0, 1000), y = runif(n, 0, 1000), elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * st_coordinates(dat)[, 1] + 2 * dat$elev + rnorm(n)
  cv <- cv_bayes(dat, "price", "elev", k = 2,
                 fit_args = list(chains = 2, iter = 500))
  cv$overall
  cv$predictive_coverage  # coverage at 50/80/95% plus mean CRPS
}

## End(Not run)

K-fold cross-validation for GWR

Description

K-fold cross-validation for GWR

Usage

cv_gwr(
  data_sf,
  response_var,
  predictor_vars,
  folds = NULL,
  k = 5,
  seed = 123,
  adaptive = TRUE,
  bandwidth = NULL,
  kernel = c("bisquare", "gaussian", "tricube", "boxcar", "exponential"),
  boundary = NULL,
  pointize = "auto",
  block_size = NULL,
  auto_range = FALSE,
  parallel = FALSE
)

Arguments

Value

A list with overall, fold_metrics, predictions, folds, formula, adaptive.

See Also

Other cross-validation: cv_bayes(), cv_spatial(), estimate_sac_range(), make_folds()

Examples

if (requireNamespace("GWmodel", quietly = TRUE) &&
    requireNamespace("sp", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = runif(n, 0, 1000), y = runif(n, 0, 1000), elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * st_coordinates(dat)[, 1] + 2 * dat$elev + rnorm(n)
  cv <- cv_gwr(dat, "price", "elev", k = 3, bandwidth = 30)
  cv$overall
  cv$fold_metrics
}

Model-agnostic spatial cross-validation

Description

Run K-fold CV for any model that returns a spatial_fit object. This is the extensibility point: to plug in a new model type, supply a fit_fn(train_sf) that returns a spatial_fit.

Usage

cv_spatial(
  data_sf,
  response_var,
  predictor_vars,
  fit_fn,
  folds = NULL,
  k = 5,
  seed = 123,
  boundary = NULL,
  pointize = "auto",
  predict_args = list(),
  fold_info_fn = NULL,
  p = NULL,
  block_size = NULL,
  auto_range = FALSE,
  parallel = FALSE
)

Arguments

Value

A list with overall, fold_metrics, predictions, folds.

See Also

Other cross-validation: cv_bayes(), cv_gwr(), estimate_sac_range(), make_folds()

Determine an optimal number of spatial levels via an elbow heuristic

Description

Computes a WSS curve over k=1..K_max using k-means on projected feature coordinates and selects candidate k values around the elbow.

Usage

determine_optimal_levels(
  data_sf,
  max_levels = 12L,
  top_n = 3L,
  sample_n = 1500L,
  set_seed = 123L,
  response_var = NULL,
  predictor_vars = NULL,
  criterion = c("geometric", "morans_i", "combined")
)

Arguments

Details

When response_var and predictor_vars are provided, the geometric WSS elbow is supplemented with Moran's I computed on OLS residuals at each candidate k. The Moran's I profile measures how much spatial autocorrelation in the response remains unexplained at a given tessellation resolution — a direct reflection of the spatial process being modeled, rather than mere geometric compactness of coordinates. The combined criterion selects the k that best balances geometric parsimony and residual spatial independence.

To keep memory use and runtime bounded for large max_levels, the initial k-means sweep records only within-cluster sum-of-squares (WSS) without retaining cluster assignments. Moran's I is then evaluated lazily: k-means is re-run only for a focused neighbourhood around the elbow (±4 by default, or ±top_n if larger), so that only the most promising candidate k values incur the cost of the full Moran's I computation.

Value

An integer vector of candidate level counts. When criterion != "geometric", an attribute "diagnostics" is attached with per-k Moran's I values.

Examples

library(sf)
set.seed(1)
# Two clearly separated clusters: the elbow should sit near k = 2
pts <- st_as_sf(
  data.frame(x = c(runif(25, 0, 10), runif(25, 90, 100)),
             y = c(runif(25, 0, 10), runif(25, 90, 100))),
  coords = c("x", "y"), crs = 32632
)
determine_optimal_levels(pts, max_levels = 6)

Ensure an object has a projected CRS (with sensible defaults)

Description

Coerces spatial objects to a projected coordinate reference system suitable for distance/area calculations.

Usage

ensure_projected(x, target_crs = NULL)

Arguments

Value

x, potentially with a new projected CRS.

Examples

library(sf)
pts_ll <- st_as_sf(
  data.frame(lon = c(9.1, 9.2), lat = c(48.7, 48.8)),
  coords = c("lon", "lat"), crs = 4326
)
st_crs(ensure_projected(pts_ll))$epsg  # auto-selected UTM zone (32632)

Create deterministic, stable polygon IDs based on spatial sort keys

Description

Ensures that a polygon layer has a reproducible, deterministic identifier column by sorting features using representative point coordinates (and secondary tie-breakers) and then assigning sequential IDs.

Usage

ensure_stable_poly_id(
  polygons_sf,
  id_col = "poly_id",
  method = c("centroid", "surface_point", "bbox_center"),
  make_valid = TRUE,
  transform_for_sort = 4326
)

Arguments

Value

An sf polygon layer re-ordered with sequential IDs in id_col.

Estimate the spatial autocorrelation range from data

Description

Fits exponential (or spherical) variogram models and returns the effective range — the distance at which the semivariance reaches ~95 \

Usage

estimate_sac_range(
  points_sf,
  response_var,
  predictor_vars = NULL,
  n_max = 5000L,
  cutoff = 0.5,
  seed = NULL
)

Arguments

Details

To guard against anisotropy, the function first estimates directional variograms at 0° (N–S) and 90° (E–W) azimuths (tolerance 22.5°, which avoids double-counting point pairs near the 45° diagonal but requires denser point clouds for stable estimates). When both fits succeed the maximum of the two directional ranges is returned, which is the conservative choice for spatial block CV — blocks must be at least as large as the longest autocorrelation range to avoid information leakage.

If either directional fit fails (e.g., too few point pairs in a direction), the function falls back to an omnidirectional (isotropic) variogram.

A log warning is emitted when notable anisotropy is detected (ratio of directional ranges > 1.5).

The returned range is in the coordinate units of the (projected) data and can be passed directly to make_folds(block_size = ...) to ensure that CV blocks are at least as wide as the autocorrelation range.

Value

A single positive numeric value (the effective range in projected coordinate units), or NA_real_ if estimation fails.

See Also

Other cross-validation: cv_bayes(), cv_gwr(), cv_spatial(), make_folds()

Examples

if (requireNamespace("gstat", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 80
  xy <- data.frame(x = runif(n, 0, 1000), y = runif(n, 0, 1000))
  xy$z <- sin(xy$x / 200) + rnorm(n, sd = 0.2)
  pts <- st_as_sf(xy, coords = c("x", "y"), crs = 32632)
  estimate_sac_range(pts, response_var = "z")
}

Compute in-sample (or out-of-sample) metrics for fitted spatial models

Description

Accepts a single spatial_fit object or a named list of them. Does NOT refit — uses fitted() for in-sample and predict() for new data.

Usage

evaluate_insample(fits, newdata = NULL, ...)

Arguments

Value

A data.frame with one row per model and columns for model name and all regression metrics.

See Also

Other model evaluation: compare_models(), compare_models_cv(), residual_morans_i()

Evaluate spatial models (legacy interface)

Description

Thin wrapper that preserves the original evaluate_models() call signature. New code should use compare_models() (for already-fit objects) or compare_models_cv() (for CV) instead.

Usage

evaluate_models(
  data_sf,
  response_var,
  predictor_vars,
  do_cv = TRUE,
  folds = NULL,
  k = 5,
  seed = 123,
  boundary = NULL,
  pointize = "auto",
  gwr_args = list(),
  bayes_args = list(),
  summary = c("mean", "median"),
  models = c("GWR", "Bayesian"),
  quiet = FALSE
)

Arguments

Value

A list with CV or in-sample results.

Cross-validated comparison with optional tessellation (legacy interface)

Description

Thin wrapper preserving the original evaluate_models_cv() call signature. New code should use compare_models_cv() directly.

Usage

evaluate_models_cv(
  data_sf,
  response_var,
  predictor_vars,
  k = 5,
  seed = 123,
  folds = NULL,
  boundary = NULL,
  pointize = "auto",
  tess_method = c("grid", "hex", "square", "voronoi", "triangles"),
  tess_args = list(),
  summary = c("mean", "median"),
  models = c("GWR", "Bayesian"),
  gwr_args = list(),
  bayes_args = list(),
  quiet = FALSE
)

Arguments

Value

A list with overall, by_fold, tessellation.

Fit a Bayesian spatial regression with a 2D Gaussian Process (via brms)

Description

Fit a Bayesian spatial regression with a 2D Gaussian Process (via brms)

Usage

fit_bayesian_spatial_model(
  data_sf,
  response_var,
  predictor_vars,
  family = NULL,
  gp_k = NULL,
  gp_c = 1.5,
  prior = NULL,
  chains = 4,
  iter = 2000,
  warmup = floor(iter/2),
  cores = max(1L, parallel::detectCores() - 1L),
  seed = 123,
  backend = c("auto", "cmdstanr", "rstan"),
  control = list(adapt_delta = 0.9, max_treedepth = 12),
  compute_loo = TRUE,
  standardize_predictors = FALSE,
  check_convergence = TRUE,
  pointize = "auto",
  boundary = NULL,
  .already_prepped = FALSE
)

Arguments

Details

Coordinate scaling and anisotropy. Before fitting the GP, X and Y coordinates are each centered and divided by their own standard deviation (lines 110–111). Because the two axes are scaled independently, an isotropic squared-exponential kernel in the scaled space corresponds to an anisotropic kernel in the original CRS: the effective length-scale in the X direction (in CRS units) differs from the Y direction whenever sd(X) != sd(Y).

This per-axis standardization is deliberate — it stabilises the GP numerically when the coordinate extents differ dramatically (common in projected CRSs where easting and northing span very different ranges) — but users who expect the GP to be isotropic in geographic distance should be aware of this behaviour.

If true isotropy in the original CRS is desired, one could use a single scaling factor such as max(sd(X), sd(Y)) for both axes. The stored $info$coord_scaling list includes a scaling_type element ("anisotropic") so downstream code can detect which strategy was used.

Value

A bayesian_fit object (inherits from spatial_fit). Supports predict(), fitted(), residuals(), coef(), summary(), and model_metrics(). Model-specific metadata lives in $info (coord_scaling, predictor_scaling, gp_k, loo, looic, convergence_ok, convergence_diagnostics). The raw brmsfit is in $engine.

See Also

Other model fitting: fit_gwr_model(), prep_model_data()

Examples

## Not run: 
if (requireNamespace("brms", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = runif(n, 0, 1000), y = runif(n, 0, 1000), elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * st_coordinates(dat)[, 1] + 2 * dat$elev + rnorm(n)
  fit <- fit_bayesian_spatial_model(dat, "price", "elev",
                                    chains = 2, iter = 500,
                                    compute_loo = FALSE)
  summary(fit)
  head(predict(fit, newdata = dat))
}

## End(Not run)

Fit a Geographically Weighted Regression (GWR) via GWmodel

Description

Fits a GWR using GWmodel on an sf dataset with either adaptive or fixed bandwidth.

Usage

fit_gwr_model(
  data_sf,
  response_var,
  predictor_vars,
  adaptive = TRUE,
  bandwidth = NULL,
  kernel = c("bisquare", "gaussian", "tricube", "boxcar", "exponential"),
  .already_prepped = FALSE
)

Arguments

Value

A gwr_fit object (inherits from spatial_fit). Supports predict(), fitted(), residuals(), coef(), summary(), and model_metrics(). Model-specific metadata lives in $info (bandwidth, adaptive, kernel, AICc). The raw GWmodel result is in $engine.

Collinearity diagnostics

The function checks the condition number of the predictor matrix and warns when it exceeds a threshold. A global condition number is computed on the full predictor matrix. In addition, a local spot-check is performed at a small random sample of locations: for each sampled point, the nearest neighbours within the bandwidth window are selected and the condition number of that local (weighted) design sub-matrix is evaluated. If the fraction of sampled locations with an extreme local condition number (> 1e6) exceeds 25\ separate warning is issued.

Because the local spot-check examines only a subset of locations (up to 30 by default), it may not detect every problematic neighbourhood. Users working with highly clustered data or near-collinear predictors should consider a full local-collinearity audit as a post-fit diagnostic.

See Also

Other model fitting: fit_bayesian_spatial_model(), prep_model_data()

Examples

if (requireNamespace("GWmodel", quietly = TRUE) &&
    requireNamespace("sp", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = runif(n, 0, 1000), y = runif(n, 0, 1000), elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * st_coordinates(dat)[, 1] + 2 * dat$elev + rnorm(n)
  fit <- fit_gwr_model(dat, "price", "elev", bandwidth = 30)
  summary(fit)
  head(predict(fit, newdata = dat))   # newdata is re-projected if needed
}

Generate seed points for Voronoi tessellation

Description

Creates an sf POINT layer of "seed" locations. Multiple strategies are supported: user-provided points, uniform random sampling within a boundary, or k-means clustering of a sampling cloud.

Usage

get_voronoi_seeds(
  boundary = NULL,
  method = c("kmeans", "random", "provided"),
  n = NULL,
  seeds = NULL,
  sample_points = NULL,
  kmeans_nstart = 10,
  kmeans_iter = 100,
  set_seed = NULL
)

Arguments

Value

An sf POINT object with seed_id and method columns.

See Also

Other tessellation: build_tessellation(), create_grid_polygons(), create_voronoi_polygons(), plot_tessellation_map()

Examples

library(sf)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
get_voronoi_seeds(bnd, method = "random", n = 5, set_seed = 1)

Heuristic length-scale bounds for a squared-exponential GP

Description

Computes sensible prior bounds for the GP length-scale parameter \ell of a squared-exponential (exponentiated-quadratic) kernel, k(h) = \exp(-h^2 / (2\ell^2)). The "effective range" where correlation drops to ~5\ \ell \sqrt{2 \ln 20} \approx 2.45\,\ell.

Usage

gp_lengthscale_bounds(coords_xy, q_small = 0.25, max_n = 1000L)

phi_prior_bounds(coords_xy, q_small = 0.25, max_n = 1000L)

Arguments

Details

Subsamples large datasets to avoid O(n^2) memory and time cost.

phi_prior_bounds is a deprecated alias retained for backward compatibility. Previous versions computed bounds calibrated for an exponential covariance (3 / d); the current implementation delegates to gp_lengthscale_bounds which is calibrated for the squared-exponential kernel used by brms::gp().

Value

Named numeric vector c(lower, upper) on the length-scale.

Examples

set.seed(1)
xy <- cbind(runif(50), runif(50))
gp_lengthscale_bounds(xy)

Harmonize CRS between two spatial objects

Description

Aligns two sf objects to a common CRS.

Usage

harmonize_crs(
  a,
  b,
  prefer = c("a", "b"),
  target_crs = NULL,
  on_transform_error = c("stop", "set_crs")
)

Arguments

Value

A named list with components a and b.

Create spatial cross-validation folds

Description

Builds train/test splits using random K-fold, spatial block K-fold, or buffered leave-one-out strategies.

Usage

make_folds(
  points_sf,
  k,
  method = c("random_kfold", "block_kfold", "buffered_loo"),
  seed = NULL,
  block_nx = NULL,
  block_ny = NULL,
  block_multiplier = 3,
  block_size = NULL,
  auto_range = FALSE,
  response_var = NULL,
  predictor_vars = NULL,
  boundary = NULL,
  buffer = NULL,
  drop_empty_blocks = TRUE
)

Arguments

Details

For block_kfold, the default grid sizing is purely geometric and unrelated to the autocorrelation range of the data. When blocks are smaller than the autocorrelation range, spatially correlated observations leak across folds and CV metrics become optimistic. Use block_size to set a minimum block edge length (in CRS units), or set auto_range = TRUE to estimate the range from an empirical variogram and enforce it automatically.

Value

A list with method, k, folds, assignment, params. The train/test elements of each fold contain ..row_id values (equal to row positions when the input has no pre-existing ..row_id column), consistent with the assignment tibble.

See Also

Other cross-validation: cv_bayes(), cv_gwr(), cv_spatial(), estimate_sac_range()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = runif(30, 0, 1000), y = runif(30, 0, 1000)),
  coords = c("x", "y"), crs = 32632
)
folds <- make_folds(pts, k = 3, method = "block_kfold", seed = 42)
folds$assignment          # fold membership per row
lengths(folds$folds[[1]]) # train/test row-ID splits

# Buffered leave-one-out: neighbours within 100 units excluded from training
loo <- make_folds(pts, k = 1, method = "buffered_loo", buffer = 100)

Compute goodness-of-fit metrics for a spatial model

Description

Compute goodness-of-fit metrics for a spatial model

Usage

model_metrics(object, ...)

## S3 method for class 'spatial_fit'
model_metrics(object, newdata = NULL, ...)

Arguments

Value

A data.frame with n, RMSE, MAE, MAPE, SMAPE, R2, Adj_R2.

Build a spatial_fit S3 object

Description

Low-level constructor used by fit_gwr_model() and fit_bayesian_spatial_model(). Users should not call this directly.

Usage

new_spatial_fit(
  subclass,
  engine,
  formula,
  response_var,
  predictor_vars,
  data_sf,
  info = list()
)

Arguments

Value

An object of class c(subclass, "spatial_fit").

Plot a tessellation map with optional boundary, seeds, and features

Description

Builds a layered ggplot2 map of polygon tessellations and optional overlays for a study boundary, seed points, and additional features.

Usage

plot_tessellation_map(
  tessellation_sf,
  boundary = NULL,
  seeds_sf = NULL,
  features_sf = NULL,
  fill_col = NULL,
  palette = "viridis",
  na_fill = "grey90",
  tile_alpha = 0.9,
  outline_col = "white",
  outline_size = 0.2,
  features_col = "#333333",
  features_size = 0.5,
  seeds_col = "#1f77b4",
  seeds_size = 1.5,
  boundary_col = "#111111",
  boundary_size = 0.6,
  labels = FALSE,
  label_col = "grid_id",
  label_size = 2.7,
  legend = TRUE,
  legend_title = NULL,
  theme = ggplot2::theme_void(),
  target_crs = NULL,
  title = NULL,
  subtitle = NULL,
  caption = NULL,
  xlim = NULL,
  ylim = NULL,
  expand = TRUE
)

Arguments

Value

A ggplot2 object.

See Also

Other tessellation: build_tessellation(), create_grid_polygons(), create_voronoi_polygons(), get_voronoi_seeds()

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  pts <- st_as_sf(
    data.frame(x = runif(20, 0, 100), y = runif(20, 0, 100)),
    coords = c("x", "y"), crs = 32632
  )
  tess <- build_tessellation(pts, method = "voronoi", quiet = TRUE)
  p <- plot_tessellation_map(tess$cells, features_sf = pts,
                             fill_col = "cell_id", legend = FALSE)
  p
}

Predict from a Bayesian spatial GP model

Description

Applies the same newdata preparation pipeline as predict.gwr_fit(): non-point geometries are coerced to points, the data is projected to the CRS used during fitting (via ensure_projected()), and rows with missing or non-finite values are dropped. Coordinate scaling and predictor standardisation stored at fit time are then applied before delegating to brms::posterior_epred() or brms::posterior_predict().

Usage

## S3 method for class 'bayesian_fit'
predict(
  object,
  newdata = NULL,
  summary = c("mean", "median"),
  type = c("epred", "predict"),
  draws = FALSE,
  ...
)

Arguments

Value

Numeric vector (or matrix when draws=TRUE).

Predict from a GWR spatial model

Description

When newdata is NULL, returns the in-sample fitted values. Otherwise uses GWmodel::gwr.predict() on the new locations. newdata is first transformed to the CRS used during fitting (via ensure_projected()), so predictions are computed in a single coordinate system regardless of the CRS newdata arrives in.

Usage

## S3 method for class 'gwr_fit'
predict(object, newdata = NULL, ...)

Arguments

Value

Numeric vector of predictions.

Prepare and sanitize an sf dataset for spatial modeling

Description

Ensures point geometry, projected CRS, and removes rows with missing or non-finite values in modeling columns. All non-POINT geometries (including MULTIPOINT) are coerced to representative points via coerce_to_points(), so downstream coordinate extraction always aligns one row per observation.

Usage

prep_model_data(
  data_sf,
  response_var,
  predictor_vars,
  boundary = NULL,
  pointize = c("auto", "surface", "centroid", "line_midpoint", "bbox_center"),
  require_response = TRUE
)

Arguments

Value

An sf object (points) in a projected CRS, cleaned.

See Also

Other model fitting: fit_bayesian_spatial_model(), fit_gwr_model()

Examples

library(sf)
dat <- st_as_sf(
  data.frame(x = 1:5, y = 5:1,
             resp = c(1, 2, NA, 4, 5),
             pred = c(1, 2, 3, 4, Inf)),
  coords = c("x", "y"), crs = 32632
)
prep_model_data(dat, "resp", "pred")  # drops rows 3 (NA) and 5 (Inf)

Compute Moran's I on the residuals of a fitted spatial model

Description

Given a spatial_fit object (GWR or Bayesian), extracts the residuals and the observation coordinates, builds a spatial weight matrix, and computes the Moran's I statistic together with its analytical expectation and variance under the randomisation assumption (Cliff & Ord). A z-score and two-sided p-value are provided so the caller can assess whether statistically significant spatial autocorrelation remains after fitting.

Usage

residual_morans_i(
  fit,
  alternative = c("two.sided", "greater", "less"),
  weights = NULL,
  k = 8L
)

Arguments

Details

By default, weights are constructed as a k-nearest-neighbour (k = 8) binary matrix, row-standardised. Users may supply their own weight matrix via the weights argument.

Value

A list with components:

observed

Numeric scalar, Moran's I statistic.

expected

Expected I under the null of no spatial autocorrelation, -1/(n-1).

sd

Standard deviation of I under the randomisation assumption.

z

Standardised z-score, (I - E[I]) / sd(I).

p_value

Two-sided (or one-sided) p-value from the normal approximation.

n

Number of observations used.

Returns NULL with a warning if computation fails (e.g. fewer than 4 valid residuals).

See Also

Other model evaluation: compare_models(), compare_models_cv(), evaluate_insample()

Examples

if (requireNamespace("GWmodel", quietly = TRUE) &&
    requireNamespace("sp", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = runif(n, 0, 1000), y = runif(n, 0, 1000), elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * st_coordinates(dat)[, 1] + 2 * dat$elev + rnorm(n)
  fit <- fit_gwr_model(dat, "price", "elev", bandwidth = 30)
  residual_morans_i(fit)  # z near 0 / p large = no residual structure
}

Summarize features by polygon/cell ID

Description

Aggregates an sf point dataset into one row per cell. By default computes counts and means, but the aggregation function is configurable.

Usage

summarize_by_cell(
  assigned_points_sf,
  response_var = NULL,
  predictor_vars = NULL,
  id_col = "poly_id",
  agg_funs = list(mean = function(x) mean(x, na.rm = TRUE)),
  cells_sf = NULL,
  deff = 1,
  quiet = TRUE
)

Arguments

Details

In addition to user-specified aggregation functions, this function always computes within-cell standard deviation (⁠..sd_<var>⁠) and standard error (⁠..se_<var>⁠) for every numeric response/predictor column, plus a cell_weight column equal to the observation count. These columns let downstream models account for the fact that a cell with 2 observations carries more aggregation uncertainty than one with 200.

Value

A tibble/data.frame (or sf if cells_sf given) with per-cell summaries including n, cell_weight, and ⁠..sd_*⁠ / ⁠..se_*⁠ columns. When deff != 1, an attribute "deff_applied" is attached to the result recording the design effect(s) used.

Spatial autocorrelation and standard-error bias

Important: By default (deff = 1), the ⁠..se_*⁠ columns are computed as sd / sqrt(n), which assumes observations within each cell are independent. When data are spatially autocorrelated — the common case for the spatial workflows this package supports — within-cell observations are typically positively correlated, so the effective sample size is smaller than n. The naive SE is therefore anticonservative (too small), and downstream weighted regressions using cell_weight or ⁠..se_*⁠ columns will produce overconfident standard errors for cells with strong intra-cell correlation.

Setting deff = "kish" applies an approximate correction using Kish's design effect. Separate intra-class correlations (ICCs) are estimated for response and predictor variables via a one-way random-effects decomposition across all cells. Each variable type's ICC is used for its own SE adjustment, and each cell's effective sample size is reduced to n_i / (1 + (n_i - 1) * rho). This is a first-order correction that does not require a full spatial covariance model but does require enough cells and observations for a stable ICC estimate. You may also pass a fixed numeric design effect (e.g. deff = 2) to uniformly inflate standard errors.

Even with the Kish correction, the adjusted SE is an approximation. For rigorous inference under spatial dependence, consider fitting an explicit spatial covariance model (e.g. via fit_bayesian_spatial_model).

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = runif(40, 0, 100), y = runif(40, 0, 100), val = rnorm(40)),
  coords = c("x", "y"), crs = 32632
)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
grid <- create_grid_polygons(bnd, target_cells = 9, type = "square")
assigned <- assign_features_to_polygons(pts, grid)

# IID standard errors (default) vs Kish design-effect adjustment
cells <- summarize_by_cell(assigned, response_var = "val", deff = "kish")
cells
attr(cells, "deff_applied")

K-means seed generation from point coordinates

Description

K-means seed generation from point coordinates

Usage

voronoi_seeds_kmeans(points_sf, k, set_seed = 456)

Arguments

Value

An sf object of k cluster center POINTs.

Random seed generation within a polygonal boundary

Description

Random seed generation within a polygonal boundary

Usage

voronoi_seeds_random(boundary, k, set_seed = 456)

Arguments

Value

An sf object of k random POINTs.