## ----setup, include=FALSE-----------------------------------------------------
# ggplot2 is a Suggests package, and every chunk below either draws something
# or feeds something that does.  R CMD build runs this file, so the whole
# vignette is gated on ggplot2 being installed rather than calling library()
# unguarded and failing the build on a machine without it.
has_ggplot  <- requireNamespace("ggplot2",  quietly = TRUE)
has_geom    <- requireNamespace("geometry", quietly = TRUE)
has_ranger  <- requireNamespace("ranger",   quietly = TRUE)
has_gstat   <- requireNamespace("gstat",    quietly = TRUE)
has_gwmodel <- requireNamespace("GWmodel",  quietly = TRUE) &&
               requireNamespace("sp",       quietly = TRUE)

knitr::opts_chunk$set(
  echo       = TRUE,
  message    = FALSE,
  warning    = FALSE,
  fig.align  = "center",
  out.width  = "100%",
  eval       = has_ggplot
)

## ----no-ggplot, echo=FALSE, eval=!has_ggplot, results='asis'------------------
# cat("**ggplot2 is not installed, so the code below is shown but not run.**",
#     "Install it with `install.packages(\"ggplot2\")` and rebuild.\n")

## ----backend-report, echo=FALSE-----------------------------------------------
data.frame(
  package   = c("ggplot2", "geometry", "ranger", "gstat", "GWmodel + sp"),
  available = c(has_ggplot, has_geom, has_ranger, has_gstat, has_gwmodel),
  used_for  = c("all plots", "true Delaunay triangulation",
                "random forest backend", "variogram / autocorrelation range",
                "GWR backend")
)

## ----load-packages------------------------------------------------------------
library(spatialkit)
library(sf)
library(dplyr)
library(ggplot2)

set.seed(42)

## ----boundary-----------------------------------------------------------------
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()

## ----fake-data----------------------------------------------------------------
n_points <- 300

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]

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

city1 <- c(1530000, 550000)   # Charlotte-ish in EPSG:2264
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 <- pmax(exp(-dist_to_city / 400000) * 5000 + rnorm(n_points, 200, 100), 10)

spatial_field <- 20 * sin(x_coords / 250000) * cos(y_coords / 250000)

y_response <- 50 +
  0.01  * elevation +
  0.005 * pop_density +
  spatial_field +
  rnorm(n_points, 0, 5)

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

## ----quick-peek, fig.height=3.6-----------------------------------------------
ggplot() +
  geom_sf(data = nc_boundary, fill = "grey95", colour = "black") +
  geom_sf(data = points_sf, aes(colour = y), size = 1.1) +
  scale_colour_viridis_c(name = "Response (y)") +
  theme_void() +
  ggtitle("Raw observation points, North Carolina")

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

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

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

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

## ----cell-id-check------------------------------------------------------------
names(tess_hex$cells)

## ----tess-tri, eval=has_ggplot && has_geom------------------------------------
tess_tri <- build_tessellation(
  points_sf, boundary = nc_boundary,
  method = "triangles", clip = TRUE, quiet = TRUE
)
nrow(tess_tri$cells)

## ----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 (has_geom) nrow(tess_tri$cells) else "skipped (install 'geometry')"
))

## ----summarize----------------------------------------------------------------
assigned <- assign_features_to_polygons(points_sf, tess_voronoi$cells,
                                        polygon_id_col = "cell_id")

cell_stats <- summarize_by_cell(
  assigned,
  response_var   = "y",
  predictor_vars = c("elevation", "pop_density"),
  id_col         = "cell_id"
)

head(as.data.frame(cell_stats)[, c("cell_id", "n", "resp_mean_y",
                                   "..sd_resp_y", "..se_resp_y")])

## ----summarize-deff-----------------------------------------------------------
cell_kish <- summarize_by_cell(
  assigned,
  response_var   = "y",
  predictor_vars = c("elevation", "pop_density"),
  id_col         = "cell_id",
  deff           = "kish"
)

deff <- attr(cell_kish, "deff_applied")
cat(sprintf("method = %s | ICC(response) = %.3f | median design effect = %.2f\n",
            deff$method, deff$icc_resp, stats::median(deff$deff, na.rm = TRUE)))

# Inflation of the response standard error, cell by cell.
summary(cell_kish$`..se_resp_y` / cell_stats$`..se_resp_y`)

## ----choropleth-helper--------------------------------------------------------
#' Assign points to cells, summarise, and draw a choropleth.
make_choropleth <- function(tess, boundary, points, title = NULL) {
  cells    <- tess$cells
  assigned <- assign_features_to_polygons(points, cells, polygon_id_col = "cell_id")

  stats_df <- summarize_by_cell(assigned, response_var = "y",
                                id_col = "cell_id")
  cells    <- left_join(cells, as.data.frame(stats_df)[, c("cell_id", "resp_mean_y")],
                        by = "cell_id")

  plot_tessellation_map(
    tessellation_sf = cells,
    boundary        = boundary,
    fill_col        = "resp_mean_y",
    palette         = "viridis",
    tile_alpha      = 0.9,
    outline_col     = "white",
    outline_size    = 0.3,
    boundary_col    = "grey20",
    boundary_size   = 0.8,
    legend_title    = "Mean y",
    title           = title,
    subtitle        = sprintf("%d cells | %d observations",
                              nrow(cells), nrow(points))
  )
}

## ----choro-voronoi, fig.cap="Voronoi choropleth"------------------------------
make_choropleth(tess_voronoi, nc_boundary, points_sf, "Voronoi tessellation")

## ----choro-hex, fig.cap="Hex grid choropleth"---------------------------------
make_choropleth(tess_hex, nc_boundary, points_sf, "Hexagonal grid")

## ----choro-square, fig.cap="Square grid choropleth"---------------------------
make_choropleth(tess_square, nc_boundary, points_sf, "Square grid")

## ----choro-tri, fig.cap="Delaunay choropleth", eval=has_ggplot && has_geom----
make_choropleth(tess_tri, nc_boundary, points_sf, "Delaunay triangulation")

## ----comparison-panel, fig.width=7, fig.height=2.8, fig.cap="All three at a glance", eval=has_ggplot && requireNamespace("patchwork", quietly = TRUE)----
library(patchwork)

(make_choropleth(tess_voronoi, nc_boundary, points_sf, "Voronoi") |
 make_choropleth(tess_hex,     nc_boundary, points_sf, "Hex") |
 make_choropleth(tess_square,  nc_boundary, points_sf, "Square")) +
  plot_annotation(title = "Tessellation comparison, cell-level mean response")

## ----sac-range, eval=has_ggplot && has_gstat----------------------------------
sac <- estimate_sac_range(points_sf, response_var = "y",
                          predictor_vars = c("elevation", "pop_density"))
sac

if (is.na(sac)) {
  cat("no identified range:", attr(sac, "rejected_reason"), "\n")
} else {
  cat(sprintf("range = %.0f ft; anisotropy ratio %.2f; used %s\n",
              as.numeric(sac), attr(sac, "anisotropy"),
              attr(sac, "anisotropy_used")))
}

## ----folds--------------------------------------------------------------------
folds_random  <- make_folds(points_sf, k = 5, method = "random_kfold", seed = 42)
folds_blocked <- make_folds(points_sf, k = 5, method = "block_kfold",  seed = 42)

c(random = folds_random$k, blocked = folds_blocked$k)

## ----plot-folds, fig.width=7, fig.height=3, eval=has_ggplot && requireNamespace("patchwork", quietly = TRUE)----
library(patchwork)
plot_folds(folds_random,  points_sf, boundary = nc_boundary) +
  plot_folds(folds_blocked, points_sf, boundary = nc_boundary)

## ----cv-contrast, eval=has_ggplot && has_ranger-------------------------------
rf_args <- list(include_coords = TRUE, num_trees = 300, seed = 1)

cv_random  <- do.call(cv_rf, c(list(points_sf, "y", c("elevation", "pop_density"),
                                    folds = folds_random),  rf_args))
cv_blocked <- do.call(cv_rf, c(list(points_sf, "y", c("elevation", "pop_density"),
                                    folds = folds_blocked), rf_args))

data.frame(
  folds = c("random_kfold", "block_kfold"),
  R2    = c(cv_random$overall$R2,   cv_blocked$overall$R2),
  RMSE  = c(cv_random$overall$RMSE, cv_blocked$overall$RMSE)
)

## ----cv-contrast-text, echo=FALSE, results='asis', eval=has_ggplot && has_ranger----
cat(sprintf(
  "Random folds report R<sup>2</sup> = %.3f; blocked folds report %.3f on the same fitted model — a drop of %.0f%% of the reported skill.\n",
  cv_random$overall$R2, cv_blocked$overall$R2,
  100 * (cv_random$overall$R2 - cv_blocked$overall$R2) / cv_random$overall$R2))

## ----rf-fit, eval=has_ggplot && has_ranger------------------------------------
rf_fit <- fit_rf_model(points_sf, response_var = "y",
                       predictor_vars = c("elevation", "pop_density"))
rf_fit

## ----rf-summary, eval=has_ggplot && has_ranger--------------------------------
summary(rf_fit)
rf_fit$info$importance     # coef() on an rf_fit errors: a forest has no coefficients

## ----rf-resid, fig.height=3.6, eval=has_ggplot && has_ranger------------------
plot(rf_fit, type = "residuals")

## ----rf-moran, eval=has_ggplot && has_ranger----------------------------------
mi <- residual_morans_i(rf_fit)
cat(sprintf("Moran's I = %.4f (z = %.2f, p = %.3g)\n",
            mi$observed, mi$z, mi$p_value))

## ----gwr-fit, eval=has_ggplot && has_gwmodel----------------------------------
gwr_fit <- fit_gwr_model(
  data_sf        = points_sf,
  response_var   = "y",
  predictor_vars = c("elevation", "pop_density"),
  adaptive       = TRUE,
  kernel         = "bisquare"
)

gwr_met <- model_metrics(gwr_fit)
cat(sprintf("Bandwidth: %.1f | in-sample R2: %.3f | RMSE: %.3f\n",
            gwr_fit$info$bandwidth, gwr_met$R2, gwr_met$RMSE))

## ----compare, eval=has_ggplot && has_ranger-----------------------------------
cmp <- compare_models_cv(
  points_sf, "y", c("elevation", "pop_density"),
  models   = c("RF", "GWR"),
  folds    = folds_blocked,
  rf_args  = list(num_trees = 300)
)
cmp$overall

## ----surface, fig.height=3.6, eval=has_ggplot && has_ranger-------------------
surf <- predict_surface(rf_fit, n_cells = 3000,
                        covariates = points_sf, boundary = nc_boundary)

ggplot() +
  geom_sf(data = surf, aes(colour = .pred), size = 0.6) +
  geom_sf(data = nc_boundary, fill = NA, colour = "grey20") +
  scale_colour_viridis_c(name = "Predicted y") +
  theme_void() +
  ggtitle("Predicted surface")

## ----aoa, fig.height=3.6, eval=has_ggplot && has_ranger-----------------------
aoa <- area_of_applicability(surf, model = rf_fit, folds = folds_blocked)

cat(sprintf("inside: %d | outside: %d | undetermined: %d | DI threshold %.3f\n",
            aoa$n_inside, aoa$n_outside, aoa$n_na, aoa$threshold))

# AOA is NA wherever a predictor was missing or non-finite, and `!NA` is NA,
# which R silently skips in a subscripted assignment. Test for TRUE and let
# anything else count as outside.
inside <- aoa$aoa$AOA %in% TRUE
surf$.pred_masked <- ifelse(inside, surf$.pred, NA_real_)

ggplot() +
  geom_sf(data = surf, aes(colour = .pred_masked), size = 0.6) +
  geom_sf(data = nc_boundary, fill = NA, colour = "grey20") +
  scale_colour_viridis_c(name = "Predicted y", na.value = "grey85") +
  theme_void() +
  ggtitle("Predicted surface, extrapolations blanked out")

## ----aoa-outside, eval=has_ggplot && has_ranger-------------------------------
extreme <- st_sf(
  elevation   = c(mean(points_sf$elevation), max(points_sf$elevation) * 4),
  pop_density = c(mean(points_sf$pop_density), max(points_sf$pop_density) * 4),
  geometry    = st_geometry(points_sf)[1:2]
)
area_of_applicability(extreme, model = rf_fit, folds = folds_blocked)$aoa[, c("DI", "AOA")]

## ----session-info, echo=FALSE, eval=TRUE--------------------------------------
sessionInfo()