# spatialkit

> Spatial tessellation, modeling, and cross-validation toolkit for R

[![R >= 4.1](https://img.shields.io/badge/R-%3E%3D%204.1-blue)](https://www.r-project.org/)
[![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](LICENSE)
[![Version](https://img.shields.io/badge/version-1.0.0-orange)](DESCRIPTION)

`spatialkit` is a modular R package for end-to-end spatial analysis workflows built on [`sf`](https://r-spatial.github.io/sf/): CRS management, Voronoi/Delaunay/grid tessellation, feature-to-polygon assignment and aggregation, geographically weighted regression (GWR via `GWmodel`), Bayesian spatial Gaussian-process regression (via `brms`), and — the part most spatial ML pipelines get wrong — **spatially aware cross-validation** with block and buffered leave-one-out strategies that respect the autocorrelation structure of the data.

All model backends return a common `spatial_fit` S3 object with consistent `predict()`, `fitted()`, `residuals()`, `coef()`, and `summary()` methods, so models can be fitted, cross-validated, and compared through one interface.

**The idea in three steps:**

1. **Cut geography into cells.** Grow Voronoi regions from k-means seeds, lay down hex or square grids, or triangulate — clipped to your study area, with stable reproducible cell IDs.
2. **Model what varies across space.** Fit geographically weighted regression or a Bayesian Gaussian process through one common interface, and aggregate observations per cell with honest, autocorrelation-aware standard errors.
3. **Validate without fooling yourself.** Build cross-validation folds that hold out whole regions sized to the data's autocorrelation range, instead of leaking spatial signal between train and test.

Three ways to cut the same geography. The first panel shows raw observations of a smooth spatial field over North Carolina (high in the west, with an eastern hotspot). The other panels show the same points aggregated into Voronoi regions grown from k-means seeds, a hex grid, and a square grid — each cell coloured by the mean of the points that fall inside it (computed with `assign_features_to_polygons()` + `summarize_by_cell()`). All panels share one colour scale, so each tessellation should look like a mosaic version of the raw data:

![Raw spatial field and three tessellation methods aggregating it over North Carolina](https://raw.githubusercontent.com/elkronos/gis_modeling_toolkit/main/man/figures/readme-tessellations.png)

### How many cells?

Resolution is a modeling decision, not a cosmetic one: too few cells smooth the signal away, too many leave each cell with a handful of noisy observations. `determine_optimal_levels()` chooses candidate cell counts for you by combining two criteria — a geometric elbow on the within-cluster sum of squares from a k-means sweep, and a model-aware check that computes Moran's I on OLS residuals at each candidate resolution. The Moran's I profile measures how much spatial structure in the response remains *unexplained* at that scale, so the combined ranking picks the resolution that balances parsimony against residual spatial independence:

```r
k <- determine_optimal_levels(pts, response_var = "price",
                              predictor_vars = "elev", max_levels = 15)
```

The same field cut at three resolutions, next to the raw observations — too coarse blurs the hotspot, too fine chases noise with near-empty cells, and the selected k preserves the trend without overfitting geography:

![Raw observations and the same spatial field tessellated at three resolutions, including the automatically selected one](https://raw.githubusercontent.com/elkronos/gis_modeling_toolkit/main/man/figures/readme-resolution.png)

## Why spatial CV?

Random k-fold CV leaks information between spatially correlated train and test observations, producing optimistic performance estimates. `spatialkit` builds folds from spatial blocks sized to the data's autocorrelation range (estimated from directional variograms via `estimate_sac_range()`), or from distance-buffered leave-one-out splits, so that reported metrics reflect true out-of-sample spatial generalization.

The difference is easy to see. On the same North Carolina sites, random folds scatter test points among their spatially correlated training neighbours, while block folds hold out contiguous regions:

![Random k-fold versus spatial block k-fold assignments across North Carolina](https://raw.githubusercontent.com/elkronos/gis_modeling_toolkit/main/man/figures/readme-spatial-cv.png)

## Installation

```r
# From GitHub
# install.packages("remotes")
remotes::install_github("elkronos/gis_modeling_toolkit")

# From a local clone
devtools::install("path/to/spatialkit")
```

Hard dependencies (`sf`, `dplyr`, `logger`, `digest`) install automatically. Optional backends are declared in `Suggests` and checked at runtime:

| Feature | Requires |
|---|---|
| GWR modeling | `GWmodel`, `sp` |
| Bayesian GP modeling | `brms` (with `rstan`, or the non-CRAN `cmdstanr` backend if installed) |
| Delaunay triangulation | `geometry` |
| Variogram-based autocorrelation range | `gstat` |
| PSIS-LOO information criterion | `loo` |
| Plotting | `ggplot2` |
| Fast k-NN weights (recommended for n > 5,000) | `FNN`, `Matrix` |

## Quick start

A self-contained example with synthetic data:

```r
library(spatialkit)
library(sf)

# --- Synthetic point data in a projected CRS -------------------------------
set.seed(1)
n  <- 200
df <- data.frame(x = runif(n, 0, 5000), y = runif(n, 0, 5000))
df$elev  <- rnorm(n)
df$price <- 50 + 0.004 * df$x + 3 * df$elev + rnorm(n)
pts <- st_as_sf(df, coords = c("x", "y"), crs = 32632)

# --- Tessellate the study area ---------------------------------------------
tess <- build_tessellation(pts, method = "voronoi")
plot_tessellation_map(tess$cells, features_sf = pts)

# --- Aggregate observations per cell ---------------------------------------
assigned <- assign_features_to_polygons(pts, tess$cells, polygon_id_col = "cell_id")
cells    <- summarize_by_cell(assigned, response_var = "price",
                              predictor_vars = "elev", deff = "kish")

# --- Spatially blocked cross-validation folds ------------------------------
folds <- make_folds(pts, k = 5, method = "block_kfold",
                    auto_range = TRUE, response_var = "price",
                    predictor_vars = "elev", seed = 42)

# --- Fit and cross-validate a GWR model ------------------------------------
fit <- fit_gwr_model(pts, response_var = "price", predictor_vars = "elev")
cv  <- cv_gwr(pts, "price", "elev", folds = folds)
cv$overall        # pooled out-of-sample RMSE, MAE, R2, ...
cv$fold_metrics   # per-fold breakdown

# --- Residual spatial autocorrelation diagnostic ---------------------------
residual_morans_i(fit)
```

Compare backends head-to-head (each is fitted and cross-validated on identical folds):

```r
comparison <- compare_models_cv(pts, "price", "elev",
                                models = c("GWR", "Bayesian"), k = 5)
comparison$overall
```

## Function overview

| Area | Key functions |
|---|---|
| CRS & geometry | `ensure_projected()`, `harmonize_crs()`, `coerce_to_points()`, `prep_model_data()` |
| Tessellation | `build_tessellation()`, `create_voronoi_polygons()`, `create_grid_polygons()`, `create_grid_polygons_cached()`, `clip_target_for()`, `ensure_stable_poly_id()` |
| Seeding & resolution | `get_voronoi_seeds()`, `voronoi_seeds_kmeans()`, `voronoi_seeds_random()`, `determine_optimal_levels()` |
| Assignment & aggregation | `assign_features_to_polygons()`, `summarize_by_cell()` |
| Modeling | `fit_gwr_model()`, `fit_bayesian_spatial_model()`, `gp_lengthscale_bounds()`; S3: `predict()`, `fitted()`, `residuals()`, `coef()`, `summary()`, `model_metrics()` |
| Cross-validation | `make_folds()`, `estimate_sac_range()`, `cv_gwr()`, `cv_bayes()`, `cv_spatial()` |
| Comparison & diagnostics | `compare_models()`, `compare_models_cv()`, `evaluate_insample()`, `residual_morans_i()` |
| Plotting | `plot_tessellation_map()` |

`cv_spatial()` is the extensibility point: pass any `fit_fn(train_sf)` that returns a `spatial_fit` object and it plugs into the same fold infrastructure, metrics, and comparison tooling.

## Prediction

Both backends support true out-of-sample prediction — `newdata` needs only the predictor columns and geometry, not the response:

```r
new_sites <- st_sf(
  elev = c(0.5, -1.2),
  geometry = st_sfc(st_point(c(1000, 2000)), st_point(c(3000, 4000)), crs = 32632)
)
preds <- predict(fit, newdata = new_sites)
```

`newdata` is automatically transformed to the CRS used during fitting, rows with missing or non-finite predictors return `NA` (output length always matches `nrow(newdata)`), and non-point geometries are coerced to representative points. The Bayesian method additionally supports `type = "predict"` for full posterior predictive draws and `draws = TRUE` for the raw draw matrix.

`model_metrics()` *does* require the response in `newdata`, since it computes error metrics against observed values.

## Parallel cross-validation

All three CV functions accept a `parallel` argument for fold-level parallelism via `parallel::mclapply()` (macOS/Linux; falls back to sequential on Windows with a message). This matters most for `cv_bayes()`, where every fold is a full MCMC run:

```r
cv <- cv_bayes(pts, "price", "elev", k = 5, parallel = TRUE)  # auto-detect cores
cv <- cv_gwr(pts, "price", "elev", k = 5, parallel = 4L)      # explicit count
```

## Diagnostics & statistical notes

**Residual Moran's I.** `residual_morans_i()` computes Moran's I on model residuals with the Cliff & Ord randomisation variance, using row-standardised k-NN weights by default (sparse via `FNN` + `Matrix` when available) or a user-supplied weight matrix (base or sparse `Matrix`). `compare_models()` runs it automatically and warns when residual spatial structure remains.

**Aggregation standard errors.** The `..se_*` columns from `summarize_by_cell()` are IID standard errors by default, which are anticonservative under within-cell spatial correlation. Pass `deff = "kish"` to apply Kish's design-effect correction from estimated intra-class correlations (separate ICCs for response and predictors), or a fixed numeric design effect. Inspect what was applied via `attr(result, "deff_applied")`.

**GWR collinearity.** `fit_gwr_model()` checks the global condition number of the predictor matrix *and* spot-checks local condition numbers within bandwidth windows at sampled locations, since spatially clustered subsets can be collinear even when the global matrix is not.

**GWR bandwidth fallback.** If automatic bandwidth selection fails, a heuristic fallback is used, a `warning()` is raised, and `fit$info$bandwidth_is_fallback = TRUE` is set so downstream comparisons can flag the result. Supply an explicit `bandwidth` if you see this.

**Bayesian GP anisotropy.** `fit_bayesian_spatial_model()` standardizes X and Y coordinates independently before the GP term. This stabilises fitting but makes the kernel anisotropic in the original CRS whenever `sd(X) ≠ sd(Y)`; the strategy is recorded in `fit$info$coord_scaling$scaling_type`. A data-informed GP length-scale prior is derived automatically from the inter-point distance distribution (see `gp_lengthscale_bounds()`).

## Logging

Detailed diagnostics are logged to a session temp file, and warnings are echoed to the console. Logging is scoped to the `"spatialkit"` namespace and never touches your global logger configuration. To customize:

```r
logger::log_appender(logger::appender_file("my_analysis.log"), namespace = "spatialkit")
logger::log_threshold(logger::WARN, namespace = "spatialkit")
```

## Documentation

- Every exported function is documented: see `?fit_gwr_model`, `?make_folds`, etc.
- A worked end-to-end demo on the classic North Carolina dataset ships as a vignette (`vignettes/spatialkit_nc_demo.Rmd`); build vignettes at install time with `remotes::install_github(..., build_vignettes = TRUE)`, then `vignette("spatialkit_nc_demo")`.
- A runnable version of the demo is installed with the package: `system.file("scripts", "example_nc_demo.R", package = "spatialkit")`.

## Development

```r
devtools::load_all()   # interactive development
devtools::test()       # run the test suite (~180 tests)
devtools::document()   # regenerate NAMESPACE + man/ from roxygen2
devtools::check()      # full R CMD check (vignette build requires pandoc)
```

The README figures are generated from actual package output; regenerate them with `Rscript dev/make_readme_figures.R`.

The test suite covers the core geometry/tessellation pipeline plus targeted regression tests for the statistical internals (Moran's I variance, CV fold/row-ID alignment, CRPS, hex-grid sizing, CRS alignment in prediction, and more). Tests that need optional backends (`GWmodel`, `brms`) skip automatically when those packages are not installed.

Contributions are welcome — please open an issue describing the bug or proposed change, and include a regression test with any fix.

## License

MIT © Justin Chase. See [LICENSE](LICENSE).

## Disclaimer

This is a personal project. It is not affiliated with, endorsed by, or connected to any organization. It uses public data sources only and was developed independently on personal time. No confidential, proprietary, or non-public information is included.
