Overview

This vignette generates synthetic spatial data over the state of North Carolina and walks through the core spatialkit workflow:

  1. Build four tessellation types (Voronoi, hex, square, Delaunay)
  2. Assign observation points to cells and compute cell-level summaries
  3. Produce choropleth maps showing cell-level mean response
  4. Fit a GWR model and visualise residuals across tessellations
  5. Run 5-fold cross-validation

Everything is self-contained — the boundary comes from the nc.shp demo shapefile bundled with the sf package, so no external files are needed.

1. Load Packages & Create Boundary

library(spatialkit)
library(sf)
library(dplyr)
library(ggplot2)

set.seed(42)

We load the North Carolina county boundaries shipped with sf, dissolve them into a single state outline, and project to NAD83 / North Carolina (ftUS) (EPSG:2264) for proper distance-based tessellations:

nc_counties <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc_boundary <- nc_counties |>
  st_union() |>
  st_transform(2264) |>
  st_as_sf()

2. Generate Synthetic Observation Data

We scatter 300 points inside the state boundary with two predictors and a spatially-varying response:

  • elevation — gradient increasing west → east, plus noise
  • pop_density — decays with distance from two fake “cities” (Charlotte, Raleigh)
  • y — response driven by predictors + a spatial sine/cosine trend
n_points <- 300

# Sample points inside the NC boundary
pts_raw <- st_sample(nc_boundary, size = n_points, type = "random")
pts_coords <- st_coordinates(pts_raw)
x_coords <- pts_coords[, 1]
y_coords <- pts_coords[, 2]

# Predictors
elevation <- scale(x_coords)[, 1] * 500 + rnorm(n_points, 3000, 400)

# Approximate projected coords for Charlotte & Raleigh in EPSG:2264
city1 <- c(1530000, 550000)   # Charlotte-ish
city2 <- c(2150000, 750000)   # Raleigh-ish
dist_to_city <- pmin(
  sqrt((x_coords - city1[1])^2 + (y_coords - city1[2])^2),
  sqrt((x_coords - city2[1])^2 + (y_coords - city2[2])^2)
)
pop_density <- exp(-dist_to_city / 400000) * 5000 + rnorm(n_points, 200, 100)
pop_density <- pmax(pop_density, 10)

# Response
y_response <- 50 +
  0.01  * elevation +
  0.005 * pop_density +
  2.0   * sin(x_coords / 300000) * cos(y_coords / 300000) +
  rnorm(n_points, 0, 5)

points_sf <- st_sf(
  y           = y_response,
  elevation   = elevation,
  pop_density = pop_density,
  geometry    = pts_raw
)

Quick sanity check — points over boundary:

ggplot() +
  geom_sf(data = nc_boundary, fill = "grey95", color = "black") +
  geom_sf(data = points_sf, aes(color = y), size = 1.2) +
  scale_color_viridis_c(name = "Response (y)") +
  theme_void() +
  ggtitle("Raw observation points — North Carolina")

3. Build Four Tessellation Types

3a. Voronoi

seeds <- get_voronoi_seeds(
  boundary      = nc_boundary,
  sample_points = points_sf,
  method        = "kmeans",
  n             = 40
)

tess_voronoi <- build_tessellation(
  points_sf, boundary = nc_boundary,
  method = "voronoi", clip = TRUE, quiet = TRUE
)

3b. Hexagonal Grid (~50 cells)

tess_hex <- build_tessellation(
  points_sf, boundary = nc_boundary,
  method = "hex", approx_n_cells = 50, clip = TRUE, quiet = TRUE
)

3c. Square Grid (~50 cells)

tess_square <- build_tessellation(
  points_sf, boundary = nc_boundary,
  method = "square", approx_n_cells = 50, clip = TRUE, quiet = TRUE
)

3d. Delaunay Triangles

tess_tri <- tryCatch(
  build_tessellation(
    points_sf, boundary = nc_boundary,
    method = "triangles", clip = TRUE, quiet = TRUE
  ),
  error = function(e) {
    message("Delaunay skipped: ", conditionMessage(e))
    NULL
  }
)

Cell counts

cat(sprintf(
  "Voronoi: %d | Hex: %d | Square: %d | Triangles: %s\n",
  nrow(tess_voronoi$cells),
  nrow(tess_hex$cells),
  nrow(tess_square$cells),
  if (!is.null(tess_tri)) nrow(tess_tri$cells) else "skipped"
))
## Voronoi: 300 | Hex: 49 | Square: 45 | Triangles: 10

4. Choropleth Maps — Cell-Level Mean Response

For each tessellation we assign observation points to cells, compute the mean response y per cell, and render a filled choropleth with a clean look.

#' Assign points → cells, compute mean, and produce a clean choropleth
make_choropleth <- function(tess, boundary, points, fill_var = "y",
                            palette = "viridis", title = NULL,
                            legend_title = "Mean Response (y)") {
  cells <- tess$cells

  # Identify the id column
  id_col <- if ("cell_id" %in% names(cells)) "cell_id" else "poly_id"
  if (!id_col %in% names(cells)) {
    cells$cell_id <- seq_len(nrow(cells))
    id_col <- "cell_id"
  }

  # Assign points to cells and compute cell-level mean
  assigned <- assign_features_to_polygons(points, cells, polygon_id_col = id_col)

  cell_summary <- assigned |>
    st_drop_geometry() |>
    group_by(.data[[id_col]]) |>
    summarise(
      fill_value = mean(.data[[fill_var]], na.rm = TRUE),
      n_obs      = n(),
      .groups    = "drop"
    )

  cells <- left_join(cells, cell_summary, by = id_col)

  # Build the choropleth via plot_tessellation_map
  plot_tessellation_map(
    tessellation_sf = cells,
    boundary        = boundary,
    fill_col        = "fill_value",
    palette         = palette,
    tile_alpha      = 0.9,
    outline_col     = "white",
    outline_size    = 0.3,
    boundary_col    = "grey20",
    boundary_size   = 0.8,
    legend_title    = legend_title,
    title           = title,
    subtitle        = sprintf("%d cells  |  %d observations", nrow(cells), nrow(points))
  )
}

4a. Voronoi Choropleth

make_choropleth(tess_voronoi, nc_boundary, points_sf,
                title = "Voronoi Tessellation — Mean Response")
Voronoi choropleth — mean response per cell

Voronoi choropleth — mean response per cell

4b. Hexagonal Grid Choropleth

make_choropleth(tess_hex, nc_boundary, points_sf,
                title = "Hexagonal Grid — Mean Response")
Hex grid choropleth — mean response per cell

Hex grid choropleth — mean response per cell

4c. Square Grid Choropleth

make_choropleth(tess_square, nc_boundary, points_sf,
                title = "Square Grid — Mean Response")
Square grid choropleth — mean response per cell

Square grid choropleth — mean response per cell

4d. Delaunay Choropleth

make_choropleth(tess_tri, nc_boundary, points_sf,
                title = "Delaunay Triangulation — Mean Response")
Delaunay choropleth — mean response per cell

Delaunay choropleth — mean response per cell

5. Side-by-Side Comparison

if (requireNamespace("patchwork", quietly = TRUE)) {
  library(patchwork)

  p1 <- make_choropleth(tess_voronoi, nc_boundary, points_sf,
                         title = "Voronoi")
  p2 <- make_choropleth(tess_hex,     nc_boundary, points_sf,
                         title = "Hex Grid")
  p3 <- make_choropleth(tess_square,  nc_boundary, points_sf,
                         title = "Square Grid")

  (p1 | p2 | p3) +
    plot_annotation(
      title    = "Tessellation Comparison — Cell-Level Mean Response",
      subtitle = sprintf("%d observations, North Carolina", n_points),
      theme    = theme(
        plot.title    = element_text(size = 16, face = "bold"),
        plot.subtitle = element_text(size = 11, color = "grey40")
      )
    )
} else {
  cat("Install 'patchwork' for the side-by-side panel: install.packages('patchwork')")
}
All tessellations at a glance

All tessellations at a glance

6. GWR Model Fitting

Fit a geographically weighted regression: y ~ elevation + pop_density.

response_var   <- "y"
predictor_vars <- c("elevation", "pop_density")

gwr_fit <- tryCatch({
  fit_gwr_model(
    data_sf        = points_sf,
    response_var   = response_var,
    predictor_vars = predictor_vars,
    adaptive       = TRUE,
    kernel         = "bisquare"
  )
}, error = function(e) {
  message("GWR skipped: ", conditionMessage(e))
  NULL
})
## Adaptive bandwidth (number of nearest neighbours): 193 AICc value: 1805.596 
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 1806.03 
## Adaptive bandwidth (number of nearest neighbours): 234 AICc value: 1806.232 
## Adaptive bandwidth (number of nearest neighbours): 167 AICc value: 1805.02 
## Adaptive bandwidth (number of nearest neighbours): 152 AICc value: 1805.208 
## Adaptive bandwidth (number of nearest neighbours): 177 AICc value: 1805.013 
## Adaptive bandwidth (number of nearest neighbours): 183 AICc value: 1805.186 
## Adaptive bandwidth (number of nearest neighbours): 173 AICc value: 1805.065 
## Adaptive bandwidth (number of nearest neighbours): 179 AICc value: 1805.174 
## Adaptive bandwidth (number of nearest neighbours): 175 AICc value: 1805.042 
## Adaptive bandwidth (number of nearest neighbours): 177 AICc value: 1805.013
cat(sprintf("Bandwidth: %.1f  |  R²: %.3f  |  RMSE: %.3f\n",
            gwr_fit$info$bandwidth,
            gwr_fit$metrics$r_squared,
            gwr_fit$metrics$rmse))

GWR Residuals — Choropleth per Tessellation

Map the mean absolute residual onto each tessellation type, producing clean filled maps so you can compare how tessellation geometry aggregates model error.

pts_with_gwr <- points_sf
pts_with_gwr$gwr_fitted   <- as.numeric(fitted(gwr_fit))
pts_with_gwr$gwr_residual <- as.numeric(residuals(gwr_fit))
pts_with_gwr$abs_error    <- abs(pts_with_gwr$gwr_residual)

tess_list <- list(
  list(tess = tess_voronoi, label = "Voronoi"),
  list(tess = tess_hex,     label = "Hex Grid"),
  list(tess = tess_square,  label = "Square Grid")
)

for (info in tess_list) {
  cells <- info$tess$cells
  id_col <- if ("cell_id" %in% names(cells)) "cell_id" else "poly_id"
  if (!id_col %in% names(cells)) {
    cells$cell_id <- seq_len(nrow(cells))
    id_col <- "cell_id"
  }

  asgn <- assign_features_to_polygons(pts_with_gwr, cells, polygon_id_col = id_col)

  cell_err <- asgn |>
    st_drop_geometry() |>
    group_by(.data[[id_col]]) |>
    summarise(mean_abs_error = mean(abs_error, na.rm = TRUE), .groups = "drop")

  cells <- left_join(cells, cell_err, by = id_col)

  p <- plot_tessellation_map(
    tessellation_sf = cells,
    boundary        = nc_boundary,
    fill_col        = "mean_abs_error",
    palette         = "magma",
    tile_alpha      = 0.9,
    outline_col     = "white",
    outline_size    = 0.3,
    boundary_col    = "grey20",
    boundary_size   = 0.8,
    legend_title    = "Mean |Residual|",
    title           = sprintf("GWR Residuals — %s", info$label),
    subtitle        = sprintf("Abs. residual aggregated to %d cells", nrow(cells))
  )
  print(p)
  cat("\n\n")
}

7. Cross-Validation (5-Fold GWR)

cv_results <- tryCatch({
  cv_gwr(
    data_sf        = points_sf,
    response_var   = response_var,
    predictor_vars = predictor_vars,
    k              = 5,
    adaptive       = TRUE
  )
}, error = function(e) {
  message("CV skipped: ", conditionMessage(e))
  NULL
})
## Adaptive bandwidth (number of nearest neighbours): 155 AICc value: 1453.163 
## Adaptive bandwidth (number of nearest neighbours): 104 AICc value: 1455.129 
## Adaptive bandwidth (number of nearest neighbours): 188 AICc value: 1452.97 
## Adaptive bandwidth (number of nearest neighbours): 207 AICc value: 1452.754 
## Adaptive bandwidth (number of nearest neighbours): 220 AICc value: 1452.48 
## Adaptive bandwidth (number of nearest neighbours): 227 AICc value: 1452.542 
## Adaptive bandwidth (number of nearest neighbours): 214 AICc value: 1452.598 
## Adaptive bandwidth (number of nearest neighbours): 222 AICc value: 1452.461 
## Adaptive bandwidth (number of nearest neighbours): 225 AICc value: 1452.429 
## Adaptive bandwidth (number of nearest neighbours): 225 AICc value: 1452.429 
## Adaptive bandwidth (number of nearest neighbours): 157 AICc value: 1463.314 
## Adaptive bandwidth (number of nearest neighbours): 105 AICc value: 1462.364 
## Adaptive bandwidth (number of nearest neighbours): 71 AICc value: 1465.435 
## Adaptive bandwidth (number of nearest neighbours): 124 AICc value: 1462.743 
## Adaptive bandwidth (number of nearest neighbours): 91 AICc value: 1463.106 
## Adaptive bandwidth (number of nearest neighbours): 111 AICc value: 1462.68 
## Adaptive bandwidth (number of nearest neighbours): 98 AICc value: 1462.531 
## Adaptive bandwidth (number of nearest neighbours): 106 AICc value: 1462.447 
## Adaptive bandwidth (number of nearest neighbours): 101 AICc value: 1462.379 
## Adaptive bandwidth (number of nearest neighbours): 104 AICc value: 1462.363 
## Adaptive bandwidth (number of nearest neighbours): 107 AICc value: 1462.543 
## Adaptive bandwidth (number of nearest neighbours): 105 AICc value: 1462.364 
## Adaptive bandwidth (number of nearest neighbours): 106 AICc value: 1462.447 
## Adaptive bandwidth (number of nearest neighbours): 105 AICc value: 1462.364 
## Adaptive bandwidth (number of nearest neighbours): 105 AICc value: 1462.364 
## Adaptive bandwidth (number of nearest neighbours): 104 AICc value: 1462.363 
## Adaptive bandwidth (number of nearest neighbours): 155 AICc value: 1446.256 
## Adaptive bandwidth (number of nearest neighbours): 104 AICc value: 1450.545 
## Adaptive bandwidth (number of nearest neighbours): 188 AICc value: 1444.664 
## Adaptive bandwidth (number of nearest neighbours): 207 AICc value: 1444.455 
## Adaptive bandwidth (number of nearest neighbours): 220 AICc value: 1444.289 
## Adaptive bandwidth (number of nearest neighbours): 227 AICc value: 1444.225 
## Adaptive bandwidth (number of nearest neighbours): 232 AICc value: 1444.144 
## Adaptive bandwidth (number of nearest neighbours): 235 AICc value: 1444.098 
## Adaptive bandwidth (number of nearest neighbours): 237 AICc value: 1444.069 
## Adaptive bandwidth (number of nearest neighbours): 238 AICc value: 1444.062 
## Adaptive bandwidth (number of nearest neighbours): 239 AICc value: 1444.055 
## Adaptive bandwidth (number of nearest neighbours): 239 AICc value: 1444.055 
## Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 1413.537 
## Adaptive bandwidth (number of nearest neighbours): 103 AICc value: 1416.397 
## Adaptive bandwidth (number of nearest neighbours): 186 AICc value: 1413.782 
## Adaptive bandwidth (number of nearest neighbours): 134 AICc value: 1412.956 
## Adaptive bandwidth (number of nearest neighbours): 122 AICc value: 1414.002 
## Adaptive bandwidth (number of nearest neighbours): 142 AICc value: 1413.02 
## Adaptive bandwidth (number of nearest neighbours): 129 AICc value: 1413.528 
## Adaptive bandwidth (number of nearest neighbours): 136 AICc value: 1413.033 
## Adaptive bandwidth (number of nearest neighbours): 131 AICc value: 1413.269 
## Adaptive bandwidth (number of nearest neighbours): 134 AICc value: 1412.956 
## Adaptive bandwidth (number of nearest neighbours): 156 AICc value: 1458.753 
## Adaptive bandwidth (number of nearest neighbours): 104 AICc value: 1460.506 
## Adaptive bandwidth (number of nearest neighbours): 188 AICc value: 1459.461 
## Adaptive bandwidth (number of nearest neighbours): 135 AICc value: 1459.209 
## Adaptive bandwidth (number of nearest neighbours): 167 AICc value: 1458.956 
## Adaptive bandwidth (number of nearest neighbours): 147 AICc value: 1458.678 
## Adaptive bandwidth (number of nearest neighbours): 143 AICc value: 1458.782 
## Adaptive bandwidth (number of nearest neighbours): 151 AICc value: 1458.647 
## Adaptive bandwidth (number of nearest neighbours): 152 AICc value: 1458.62 
## Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 1458.726 
## Adaptive bandwidth (number of nearest neighbours): 152 AICc value: 1458.62
cat(sprintf("CV RMSE: %.3f  |  CV R²: %.3f  |  CV MAE: %.3f\n",
            cv_results$summary$rmse,
            cv_results$summary$r_squared,
            cv_results$summary$mae))

Summary

Tessellation Cells Notes
Voronoi 300 Adapts to point density via k-means seeds
Hex grid 49 Uniform hexagons, good for regular sampling
Square grid 45 Simplest regular grid
Delaunay 10 One triangle per point triplet, finest resolution

The choropleth maps show how each tessellation aggregates the response variable spatially across North Carolina. The GWR residual maps let you compare how each tessellation captures model error — hexes and squares give a uniform view, while Voronoi cells highlight where the observation network is dense or sparse.

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Chicago
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] patchwork_1.3.2  ggplot2_4.0.3    dplyr_1.2.1      sf_1.1-2        
## [5] spatialkit_1.0.0
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6           xfun_0.59              bslib_0.11.0          
##  [4] lattice_0.22-9         LearnBayes_2.15.2      vctrs_0.7.3           
##  [7] tools_4.6.1            generics_0.1.4         sandwich_3.1-2        
## [10] spdep_1.4-2            parallel_4.6.1         tibble_3.3.1          
## [13] proxy_0.4-29           spacetime_1.3-3        DEoptimR_1.2-0        
## [16] xts_0.14.2             pkgconfig_2.0.3        Matrix_1.7-6          
## [19] KernSmooth_2.23-26     data.table_1.18.4      RColorBrewer_1.1-3    
## [22] S7_0.2.2               lifecycle_1.0.5        compiler_4.6.1        
## [25] farver_2.1.2           deldir_2.0-4           FNN_1.1.4.1           
## [28] codetools_0.2-20       marginaleffects_0.32.0 htmltools_0.5.9       
## [31] class_7.3-23           sass_0.4.10            yaml_2.3.12           
## [34] pillar_1.11.1          jquerylib_0.1.4        MASS_7.3-65           
## [37] classInt_0.4-11        cachem_1.1.0           wk_0.9.5              
## [40] spatialreg_1.4-3       multcomp_1.4-31        boot_1.3-32           
## [43] nlme_3.1-169           robustbase_0.99-7      tidyselect_1.2.1      
## [46] digest_0.6.39          mvtnorm_1.4-2          splines_4.6.1         
## [49] labeling_0.4.3         GWmodel_2.4-1          fastmap_1.2.0         
## [52] grid_4.6.1             cli_3.6.6              logger_0.4.2          
## [55] magrittr_2.0.5         survival_3.8-6         TH.data_1.1-5         
## [58] e1071_1.7-17           withr_3.0.3            backports_1.5.1       
## [61] scales_1.4.0           sp_2.2-3               spData_2.3.5          
## [64] rmarkdown_2.31         otel_0.2.0             zoo_1.8-15            
## [67] coda_0.19-4.1          evaluate_1.0.5         knitr_1.51            
## [70] viridisLite_0.4.3      s2_1.1.11              rlang_1.3.0           
## [73] Rcpp_1.1.1-1.1         glue_1.8.1             DBI_1.3.0             
## [76] rstudioapi_0.19.0      jsonlite_2.0.0         R6_2.6.1              
## [79] intervals_0.15.5       units_1.0-1