                     Changes in version 2.0.0.9000                      

New features

  - fold_separation() measures what blocked cross-validation is for: the
    distance from every held-out point to its nearest training point,
    summarised per fold. Until now the only evidence a fold scheme
    worked was a proxy --- the block size compared against the estimated
    autocorrelation range, which make_folds() warns about. That is a
    statement about the design; this is the result, and the two can
    disagree in the direction that matters. Measured on 300 points over
    a 1000-unit square with a fitted range of 292.5: blocks of 343 units
    are wider than the range, so make_folds() raises no warning at all,
    and yet 76 percent of the held-out points still sit closer to a
    training point than the range, the nearest of them 23.3 units away.
    Blocks wider than the range leak wherever a test point sits near a
    block edge with training data just across it, which in a fine grid
    is most of them. The returned table carries n_train, n_test,
    n_blocks, min_dist, median_dist and within_range (the share inside
    the range), and its print() method closes with what that share means
    for the score. Nothing is estimated: the distances come from the
    geometry and the range is the one the folds already carry or the one
    you pass. A range that records its CRS (the folds' own
    params$sac_range and params$crs, or an estimate_sac_range() result)
    is compared with distances measured in that CRS, so a copy of the
    layer in other units gives the same answer (a US-foot copy used to
    report 19 percent of the hold-out inside a metre range where the
    metre layer reported 89); sac and the distances are in the units of
    the CRS the crs attribute names, metres for lon/lat input. fold is
    the fold_id a cv_*() result's $folds carries, so a join on it pairs
    each fold's error with its own distances (it was the list position,
    which after a dropped fold paired one fold's error with another's
    distances). A units object for sac is refused. As in cv_*(), a
    make_folds() result whose recorded rows sit at other locations in
    data_sf is refused. Folds built on prep_model_data()'s 200 points
    and measured on the 195 that assign_features_to_polygons() kept used
    to be matched by position, and reported 100 percent of the hold-out
    inside the range against 23--48 percent on the right layer. The
    location check is skipped when one layer is POINT and the other is
    not, so folds built on polygons still measure their pointized copy.
    The print() verdict now depends on the fold scheme: "widen the
    blocks" is said only of blocked folds, buffered leave-one-out is
    told to widen the buffer, and random, leave-location-out and
    hand-made splits are told to use blocked or buffered folds. NNDM
    folds are no longer called optimistic: they are built to reproduce
    the prediction-to-data distances, so the share describes the
    prediction task. The range line gives the CRS's unit ("in metres"),
    where it used to print the CRS code as if it were a unit ("in
    EPSG:32617 units").

  - summary() on a resolution_profile() puts every criterion's pick in
    one table: the level each prefers, the flat region around it,
    whether a ladder bound is doing the choosing, and the levels that
    lie in every band. Reading the criteria off a profile took one
    select_resolution() call per criterion, and the comparison had to be
    assembled by hand. The intersection of the bands is often empty,
    which is a result rather than a failure: it says the field has no
    single resolution that satisfies every way of asking. Nothing in the
    table chooses, and the help page says so. Each region comes back in
    full on the "bands" attribute.

  - Every cv_*() result now says what became of each fold. fold_status
    is a data.frame with one row per fold supplied --- fold, status,
    message --- where status is "ok", "error" (the fit or its predict()
    threw; message is the error text), "skipped" (nothing scorable: too
    few matched rows, a prediction of the wrong length, no finite
    observed/predicted pair), "dropped" (an empty test set or fewer than
    two training rows once incomplete rows were removed, so the fold
    never reached the fitter) or "worker_error" (a parallel worker
    died). The per-fold error text was already collected and thrown away
    except when every fold failed, so a partial failure --- "3 of 5
    folds produced predictions" --- left its causes only in console
    scrollback, which a script, a callr job or a knitted document does
    not keep. This is worth most where a run is expensive: a cv_bayes()
    fold whose sampler failed now names the reason in the returned
    object. Beside it, orphan_rows holds the row IDs no fold names (they
    enter no training set and are never scored --- non-empty only when
    the folds were built on a different or subsetted layer),
    n_unknown_ids counts fold entries naming rows the data does not
    have, and n_dropped the rows prep_model_data() removed before any
    fold was fitted. The four together account for every row and every
    fold, so n_folds_attempted - n_folds_succeeded never has to be
    explained from the log. A parallel worker that was killed outright
    (for lack of memory, say) is a "worker_error" too and enters the
    first-error text; it used to show as "skipped" ("no result
    returned") and was left out of that text. An "ok" row carries a
    message when the fold's fold_info_fn failed.

  - prep_model_data() records what it removed. attr(x, "dropped") is a
    list with n, n_geometry, which (positions in the input), row_id
    (when the layer carries ..row_id) and reason, one of "geometry",
    "missing" or "non_finite" per dropped row. The three masks behind
    that decision were already computed and collapsed into a log line;
    the row identities reached nothing, and through the eight-plus
    internal call sites even the count was invisible, so a fit's $n was
    the post-cleaning row count with nothing saying how many rows were
    lost or why. Every fit now carries the count as $info$n_dropped, and
    every cv_*() result as n_dropped. Silently losing a third of the
    rows is a classic cause of a suspiciously good score. The layer
    carries the class "spatialkit_rows" after "sf" (c("sf",
    "spatialkit_rows", "data.frame")): [ returns a plain layer without
    the record, and dplyr::bind_rows(), vctrs::vec_rbind() and
    dplyr::union_all() return a plain sf (with the class ahead of "sf"
    they failed on two layers with the same record, 'attr(obj,
    "sf_column") does not point to a geometry column'). dplyr::filter(),
    slice(), arrange() and distinct() now drop the record as [ does
    (they kept it, so filter() down to 3 of 5 rows still reported the
    parent's counts). Each record carries n_rows, the number of rows it
    was made for. sf::st_drop_geometry() keeps the rows and the record;
    binding such data frames keeps the first one's record, whose n_rows
    then no longer matches, and the package's readers ignore it.

  - make_folds(method = "block_kfold") returns the block design it built
    the folds from. assignment gains a third column, block_id; params
    gains blocks (an sf layer of the block polygons in the CRS the folds
    were built in, numbered to match), block_sizes (points per block,
    indexed by block_id, so a block kept empty by drop_empty_blocks =
    FALSE shows as a zero) and fold_blocks (which blocks were packed
    into each fold). blocks$source_row is the row each block came from
    in the layer it originated in --- a cell's index in the full grid_nx
    by grid_ny grid, or the row of the blocks argument --- because
    dropping the empty blocks renumbers the rest: nine supplied blocks
    of which three are empty come back as six rows numbered 1 to 6, and
    a join by row position would mis-attribute every block after the
    first gap. blocks[params$blocks$source_row, ] recovers them with
    their own columns and in their own order. The folds account for
    every block exactly once, empty ones included, so a fold's territory
    on the map is all of its blocks rather than only those holding
    points; and blocks_used is the number of blocks the design has,
    which equals nrow(params$blocks) and length(params$block_sizes),
    with sum(params$block_sizes > 0) giving how many of them hold
    points. The grid is the design of a blocked cross-validation:
    without it a user could not draw the blocks over their data, see
    that 40 of 64 blocks were empty, or tell whether a fold is one
    contiguous region or several. plot_folds() now draws those outlines
    under the points when the folds carry them, and takes blocks = FALSE
    to suppress them. Its subtitle states the parameter that decides
    whether the scheme leaks (the block size, the buffer, the number of
    location groups or the median NNDM exclusion), with the units of the
    folds' CRS; folds built from points without a CRS get no units
    rather than "(NA units)".

  - estimate_sac_range()'s returns are uniform. The four directional
    variograms and their fits run unconditionally on every call, and the
    rejected-range paths used to discard them --- precisely where a user
    most needs to know whether the field is anisotropic. Every classed
    return now carries directional, anisotropy and anisotropy_used, and
    three new attributes report the sweep rather than collapsing it:
    directional_status (per azimuth, why that direction is NA in
    directional --- "ok", "over_cutoff", "not_converged" or "no_fit",
    which were indistinguishable before), directional_fitted (the range
    each direction's fit reported whether or not it was usable, so a
    refused directional range --- the most informative number in an
    anisotropic failure --- stays recoverable) and, under
    keep_directional_fits = TRUE, directional_fits (each direction's
    empirical variogram and fitted model). Those four objects are off by
    default because they dominate the result when present --- 42.2 KB of
    a 58.8 KB object at n = 400, against 16.6 KB without them --- and
    make_folds(auto_range = TRUE) calls this on every build. print()
    names the reason and the refused value for a direction it cannot
    use.

  - A roster of quantities the package already computed and dropped are
    now returned. summarize_by_cell() attaches the Kish ICCs it
    estimated as attr(, "icc") whether or not either was large enough to
    apply --- the case with no "deff_applied" is exactly the one where a
    user wants to know what the ICC came out as --- and the "variogram"
    path adds deff_rows, the per-cell design effect at the cell's row
    count that the log line reduced to a median and a max. A
    summarize_by_cell() call that requests a design effect (any deff
    other than 1) also returns a logical deff_applied column, TRUE on
    every row when the correction was applied and FALSE when it fell
    back to the uncorrected standard errors; unlike the "deff_applied"
    attribute it survives rbind() and dplyr::bind_rows() of many
    results. The default frame is unchanged.
    assign_features_to_polygons() reports the features that matched more
    than one polygon and had the tie_break rule decide for them, as
    attr(, "ties") and a log line; a tie-break firing on a third of the
    features means the polygon layer overlaps and every cell count built
    from it is suspect. The layer carries the "spatialkit_rows" class
    after "sf", as prep_model_data()'s does, so binding such layers
    returns a plain sf. The "ties" record is dropped by the same dplyr
    row verbs as prep_model_data()'s, and with largest = TRUE it also
    counts features that overlap two or more polygons by exactly the
    same area. ensure_projected() attaches crs_choice, the projections
    it considered with each one's measured worst-case distance error.
    Every path that picks a local projection reports what it picked, the
    two that compare nothing included: a UTM zone on a local extent, and
    the equal-area projection chosen for a layer straddling the
    antimeridian. residual_morans_i() returns the residual kurtosis the
    randomisation variance conditions on, the design rank p behind the
    residual moments, exact, whether those moments are exact for these
    residuals, and weights_summary, a description of the weight matrix
    it used: n, storage (the matrix class), neighbours (the smallest and
    largest number of neighbours any row has, NA for a dense matrix,
    where counting them would allocate a second one), kept and desc, the
    line print() shows. The matrix itself comes back as weights under
    keep_weights = TRUE and is NULL otherwise, because it is n by n: at
    n = 500 it is 50.3 KB of a 53.5 KB object in its sparse form, 112.7
    KB at n = 120 when the dense fallback is taken (going sparse needs
    both FNN and Matrix, so a no-Suggests install always falls back),
    and 191 MB at the n = 5000 that fallback is capped at --- against
    the 3.2 KB everything else occupies --- and scoring a list of fits
    would hold one matrix per fit. The result is classed "morans_i" and
    prints through a print() method that shows the statistic, its null
    and the weights line; [ drops the class, and $, [[ and unlist() read
    the result exactly as for a plain list. fit_gwr_model() keeps
    info$nonfinite_coef, the per-row, per-term mask behind
    n_local_singular, and fit_bayesian_spatial_model() keeps rhat_failed
    / neff_failed, the parameters that failed each convergence check by
    name --- "max R-hat 1.09" is not actionable where "sdgp_gp..x..y has
    R-hat 1.09" is. determine_optimal_levels()'s model-aware diagnostics
    gain knee_k and failed_k (whose interpolated WSS entries are not
    measurements). build_tessellation() records the points whose cell
    assignment was repaired by nearest-cell snapping, and how far
    outside each sat, as params$snapped --- a comment had long said it
    should. area_of_applicability() returns the scaling (per-predictor
    training centre and SD) the dissimilarity index is computed in,
    without which a location's DI cannot be traced to the predictor that
    put it outside, and n_outliers, the training DI values the
    threshold's fence set aside. get_voronoi_seeds(method = "kmeans")
    returns the clustering as attr(, "kmeans") --- which cloud points
    fed which seed, the cluster sizes and the within-cluster sums of
    squares. gwr_model_selection() reports criterion_by_name,
    criterion_column and criterion_verified, so a script can gate on the
    case its log calls "unverified" instead of reading the label.

  - estimate_sac_range()'s result gains a plot() method.
    plot(estimate_sac_range(pts, "z")) draws the empirical variogram,
    with the fitted model and the effective range overlaid where a range
    was identified, and a subtitle saying why not where it was not: the
    variogram never reached a sill, both model fits were singular, the
    optimiser halted, or a range_frac below 1 refused a range inside the
    fitted lags. Nothing is recomputed --- the plot reads the attributes
    the estimate already carries --- and the same drawing routine now
    serves plot.spatial_fit(type = "variogram"), so the two pictures
    agree. The "attached for inspection" messages estimate_sac_range()
    logs when it returns NA used to point at plot(type = "variogram"),
    which is the method for a fitted model and could not take the
    estimate; they now point at plot() on the returned value.

  - summarize_by_cell() gains conf_level. With conf_level = 0.95, every
    numeric response and predictor column gets four more columns beside
    its ..sd_* and ..se_*: ..neff_*, the column's effective sample size
    in the cell (its non-missing count over its design effect --- the
    per-column version of cell_weight); ..df_*, the degrees of freedom
    the interval uses; and ..ci_lo_* / ..ci_hi_*, a t interval for the
    cell mean as an estimate of the grand mean, built on the
    design-effect-corrected standard error. The default conf_level =
    NULL returns exactly the frame it always did. The degrees of freedom
    are n - 1 at deff = 1, for a numeric deff and for deff = "kish",
    because the interval's spread is estimated from the within-cell
    variance, whose n - 1 df survive exchangeable correlation whatever
    the design effect (measured 95% coverage on the Kish path at an ICC
    of 0.2 / 0.6 / 0.9: 0.954 / 0.953 / 0.952; an effective-sample-size
    df of neff - 1 gives 0.992 / 1.000 / 1.000 and is not used). Under
    deff = "variogram" the df are the Satterthwaite (1946)
    moment-matched df of the within-cell variance under the fitted
    correlation, a fraction of n - 1 that shrinks with the range: 0.960
    and 0.958 coverage at exponential ranges of 150 and 400 on
    a 1000-unit domain, against 0.931 and 0.918 with n - 1. The help
    page's "Confidence intervals" section has the reasoning and the
    numbers. ..neff_*, like ..df_* and the interval, is NA for a column
    with a single non-missing value in the cell, where cell_weight still
    counts it.

  - compare_models_cv() gains block_size and auto_range, and now hands
    response_var and predictor_vars to make_folds() when it builds the
    shared fold set. With the defaults the folds are the same geometric
    blocks as before, so an existing comparison does not move; what
    changes is that the fold-leakage diagnostic (above) can now fire for
    the one function that compares models, which until now was the one
    whose folds could never be checked against the range.
    select_features_forward() gains auto_range for its inner folds for
    the same reason, and records it in $params. A shared fold set that
    cannot be built is an error. It used to be a log line and a
    fall-back to each backend's default blocks, which were given neither
    argument: block_size = 1e6 (one block, which cv_rf() refuses)
    returned a five-fold comparison on geometric blocks, and auto_range
    = TRUE the small blocks it exists to prevent, with no R condition.

  - compare_models_cv()$overall carries the Bayesian backend's
    calibration when a Bayesian model ran: coverage_50, coverage_80,
    coverage_95 and mean_CRPS, the same fold-weighted summary cv_bayes()
    returns as predictive_coverage, with NA on the GWR and RF rows. A
    model that predicts well on average and covers badly (Heaton et
    al. 2019) is now visible in the table a user picks from, not only in
    $bayes_cv. model stays the last column.

  - fit_rf_model() gains replace and sample_fraction, which reach
    ranger::ranger() as replace and sample.fraction; both were already
    accepted through ..., but are now recorded in $info and printed with
    the fit ("Sampling: bootstrap, with replacement (100.0% of rows per
    tree)"). The defaults are ranger's, so no forest changes. The help
    page carries Strobl et al.'s (2007) case for replace = FALSE.
    Passing the ranger spellings through ... is now refused like the
    other arguments the wrapper sets. replace = FALSE with
    sample_fraction = 1 grows every tree on every row, so no row is out
    of bag. The fit now warns that fitted(), the OOB error and the
    permutation importance are all NaN (they were NaN in silence), and
    print() on such a fit says the OOB error and the permutation
    importance are undefined, where it dropped the OOB line and printed
    the importance line empty. area_of_applicability(weights =
    pmax(imp, 0)) failed on that importance with a hint to use the very
    pmax() it had been given, because pmax() keeps NaN; it now names the
    predictors whose weight is NaN, says why (no row is out of bag), and
    says what to do instead: refit with out-of-bag rows, or pass weights
    = NULL. With permutation importance, the fit's own warning now says
    as much: area_of_applicability() cannot be weighted by that
    importance, so pass weights = NULL.

  - area_of_applicability() records the method of the folds its
    threshold came from as $params$folds_method ("block_kfold",
    "random_kfold", ... from a make_folds() result; "labels" or "splits"
    when the input cannot say) and prints it, because the threshold is a
    statistic of a cross-validated hold-out and pairs with the CV error
    from the same kind of hold-out. An AOA built on random_kfold folds
    logs a caution saying it does not belong beside a blocked cv_*()
    result.

  - estimate_sac_range() gains detrend = c("ols", "reml") and
    reml_max_n. A variogram fitted to least-squares residuals
    underestimates the range, because the trend fit absorbs part of the
    long-wavelength variation (Lark, Cullis and Welham 2006). Measured
    for this estimator on simulated fields (n = 300, true effective
    range 300), as the median ratio to the estimate from the true field:
    a white-noise covariate 1.00; a spatially smooth covariate 0.97; a
    linear trend in the coordinates 0.92; a quadratic one 0.75. detrend
    = "reml" fits the trend and an exponential-plus-nugget covariance
    together by REML with nlme::gls() and returns the REML range
    --- 0.95 to 1.06 on the same designs --- with the empirical
    variogram of the REML residuals attached for inspection. It is cubic
    in n, so it runs on at most reml_max_n (400) points; a fit that does
    not converge falls back to OLS with a warning. The default stays
    "ols", so nothing changes unless asked; the help page's new section
    carries the numbers. Iterating GLS trend fits against variogram
    refits (Neuman and Jacobson 1984) was measured too and recovers only
    part of the bias (0.80 in the quadratic case), so it was not added.
    nlme joins Suggests. A fitted range too short for enough pairs of
    points to lie inside it is refused, with rejected_reason = "fitted
    range is below the shortest lag fitted" and the bound it fell short
    of as the range_floor attribute. For the least-squares fits the
    bound is the shortest lag the empirical variogram resolves (the mean
    separation in its first bin); it did not fire on white noise or on
    fields with a 60 m range. The REML range is fitted to the point
    pairs, not to those bins, so its bound is the distance within which
    30 pairs of the points the REML fit used lie, or that first lag when
    it is shorter. On white noise (n = 300, 30 draws) the REML fit
    returned ranges of 0.18--23.6 m in 19 draws (one of 0.27 m sized
    a 3642 x 3676 block grid, refused as a unit mistake); the 16 of them
    up to 12.5 m are refused and the estimate is finite in 13 of 30.
    REML estimates of a true 30 m range that the first lag alone would
    have refused (14.5--28.8 m, in 10 of 20 draws at n = 400 and 6 of 15
    at n = 300) are returned, and at n = 400 the same holds at cutoff
    = 0.5 and 0.1. The bound is about identification, not a test for
    spatial structure. A refused REML result keeps its reml list.

  - sac_nugget() returns the nugget variance behind an
    estimate_sac_range() result, and every classed result --- identified
    or rejected --- now carries it as a nugget attribute (NA when no
    model could be fitted). It was reachable before only by reading the
    Nug row of the attached gstat model. Results also record
    detrend_method ("ols", "reml" or NA), and with detrend = "reml" a
    reml list (n_used, subsampled, nugget_prop, sigma2).

  - estimate_sac_range() refuses a range fitted through an empirical
    variogram that decreases with distance over its shorter lags: a net
    fall of more than 15% of the mean semivariance there, weighted by
    pairs. A model that rises to a sill has nothing to identify on such
    a curve, and the result is NA with rejected_reason = "empirical
    variogram decreases with distance", the refused value in
    rejected_range, and the variogram attached; plot() captions it. The
    shape is what a periodic (hole-effect) structure or a variance that
    differs between a dense cluster and the rest of the layer produces
    --- measured on 60 draws each: 98% of fields with a periodic
    component and 100% of the clustered case are flagged, against 0% of
    an exponential field with a range of 100 or more on a 1000-unit
    extent, 2--3% at very short ranges, 2% of white noise --- and not
    what a trend produces, which is a variogram that rises without a
    sill and is refused as before. The one other path that returned a
    bare NA from inside a completed fit (a non-positive fitted range)
    now returns the classed, inspectable shape too.

  - resolution_profile() and select_resolution(): the number of cells,
    scored on every criterion at once. resolution_profile() runs a
    log-spaced ladder of level counts from a floor the autocorrelation
    range implies (ceiling(area / range^2)) to a ceiling the support
    implies (floor(n / min_cell_n), where n is every point in the layer,
    not the sample_n subsample the k-means runs on), fits each level as
    the best of 25 k-means++ restarts, and returns a data.frame with the
    WSS elbow statistic (read on log-log axes, and NA when the points
    have no cluster structure to bend the curve), Mallows' C_p of the
    piecewise-constant approximation of the response (or of its
    residuals on the predictors, fitted on the rows with a complete
    response and predictors, with the nugget from the fitted variogram
    as the noise variance and the penalty set for the whole layer), the
    standardised residual Moran's z of the cell means, and an analytic
    reliability of the cell means --- the share of their spread that is
    between-cell signal rather than sampling noise, from the variogram
    alone via Krige's additivity relation (Cressie 1996), the shrinkage
    factor of Fay and Herriot (1979) --- plus the cell-support and
    cell-diameter columns and the between-restart spread. A floor above
    the ceiling is reported as a finding rather than resolved silently.
    select_resolution() reads a level off one criterion together with
    its flat region and says when a bound, not the criterion, is
    choosing. The flat region is a set, not an interval: the criterion
    curves are not monotone, so a region can skip a rung of the ladder,
    and it is printed as the runs the criterion accepts ("19
    to 21, 26, 31 to 33") rather than a range that would quietly include
    the levels it rejected. Two things measured before this shipped,
    both on the help page: on smooth fields with a small nugget C_p
    descends to the support ceiling (every replicate at effective
    ranges 90--900 with nugget 0.3 on a unit sill; interior only at
    nugget 2), and the reliability optimum agrees with the empirical one
    from true block means on simulated fields but is broad --- flat to
    within 2 percent over a factor of 3--6 in the number of cells. Read
    the flat region. determine_optimal_levels() is unchanged in shape
    and keeps its integer-vector interface. A supplied sac whose range
    was refused (an NA with a rejected_reason) is not read as accepted:
    cp keeps that fit's nugget and warns, naming the reason, reliability
    is NA (it needs the range), and a fit that did not converge, or
    whose range is below the shortest lag fitted, gives neither its
    nugget nor its range (the whole refused fit used to be used, and a
    correlation function whose range could be many times the extent
    pinned the reliability optimum to the first level, with no R
    warning). A range below the shortest lag means the structure cannot
    be told from a nugget, so the nugget is not identified either: on
    white noise detrended by REML it was 6e-7 on a sill of 0.99, and
    C_p, with no penalty, ran to the support ceiling (33 cells). A
    nugget of 0 (under 1e-4 of the sill, since a REML fit stops short of
    its bound: 6e-7 passed a test for exactly 0), on which C_p has no
    penalty and falls to the ceiling, is warned about, and so is a sac
    whose detrended flag does not match the variable scored (a residual
    variogram on the raw response moved the C_p pick from 2--4 cells to
    the ceiling of 44 in five of five simulated fields); attr(x,
    "variogram") records detrended. When no sac is passed, these
    warnings name the variogram the profile estimated (kept in attr(x,
    "sac")), not a sac argument the caller never gave. A supplied sac is
    read in its own CRS: the points are transformed to attr(sac, "crs")
    first, as summarize_by_cell() does (a range in US feet put the floor
    of a metre layer at 2 where the same range in metres put it at 8). A
    sac given as a units object (set_units(1.5, "km")) or a character
    string is refused by name; a units object used to be read as a
    number in the CRS units, so 1.5 km became a range of 1.5 m and a
    floor of 41 million cells. A plain number is taken as the range
    alone: it sets the floor, and reliability is NA. Wherever the
    variogram gives no nugget (none was fitted, its fit did not
    converge, its range is below the shortest lag, or sac is a range
    alone), cp takes Mallows' own noise variance, the residual mean
    square RSS / (m - L) of the finest level whose cells hold two scored
    rows each on average, and the profile warns; attr(x, "cp_noise")
    records which variance was used, and the print says so. It counts
    the structure within those cells as noise too, so on average it is
    no smaller than the nugget and errs towards fewer cells. cp used to
    be NA at every level in all four cases, so a workflow that reads C_p
    by default (build_tessellation(approx_n_cells = <profile>),
    summary()) fell through to another criterion with a log line at
    most. Reliability's domain term is taken over the convex hull the
    area is measured on, not the bounding box: on a 3000 x 120 strip the
    reliability pick is 6 cells whether the strip lies axis-aligned or
    rotated by 45 degrees (it was 8 and 2), and 10 either way on a
    square (it was 10 and 8), and reliability values shift slightly on
    every profile. Row order does not change the profile (see the
    determine_optimal_levels() item under Bug fixes); over permutations
    of one 2000-point layer, WSS used to move by up to 2.6 percent and
    the C_p pick across its flat region (222, 173 and 135 cells). The
    ceiling reaches the number of distinct locations when locations
    repeat (and explicit levels up to it are kept) instead of stopping
    one short; with no location repeated, the print says the ceiling is
    one short of the points rather than crediting the distinct
    locations. Points with empty or non-finite coordinates are dropped
    with an R warning, as in determine_optimal_levels(); it was a log
    line alone. range_floor = FALSE starts the ladder at 2 whatever the
    range and only reports the floor, so profiles whose range estimates
    differ (cross-validation folds) have comparable ladders; with the
    floor applying on one fold and not the next, reliability and the
    elbow moved by a factor of 10--15 between folds. The default,
    range_floor = TRUE, keeps the floor, and the help page describes
    that regime switch. Under select_on = "split" a supplied sac is
    flagged in the log, since it must be fitted on the selection half,
    and the help page shows the two-call workflow.

  - select_on = c("all", "split") on determine_optimal_levels(),
    resolution_profile() and select_features_forward(). Whenever a
    selection reads the response --- a level count chosen with
    response_var and predictor_vars, or a predictor set chosen by a
    forward sweep --- what is estimated afterwards on the result is
    post-selection, and its standard errors are descriptive rather than
    at nominal coverage (Gao, Bien and Witten 2022). select_on = "split"
    is sample splitting: the layer is cut into two spatially blocked
    halves (make_folds(k = 2, method = "block_kfold")), the selection
    reads the response on the first only, and the row positions of both
    come back (as a "split" attribute on the first two functions, as
    $split on the third) so the estimation can be done on the half the
    selection never saw. The positions index the layer as passed, for
    all three functions (select_features_forward()'s were positions
    after its completeness filter, so with seven incomplete rows
    pts[fs$split$estimation, ] held 47 rows of the selection half and 4
    of the dropped ones). When the default six-block split leaves fewer
    than 10 points in a half (a small layer, or a small group far from
    the rest), it is retried on grids of about 16, 36 and 100 blocks,
    with a warning, instead of stopping, and the split records grid,
    n_blocks, balance (a 2:1 split on clustered layers used to pass
    silently) and each half's extent. The halves are not independent:
    blocking reduces the dependence across their shared border but does
    not remove it (85 percent of the estimation points of the package's
    test layer lie within the fitted range of a selection point), and
    the seed decides only which side selects, not where the cut falls.
    The level-count functions still draw their cells on every point, so
    the count they return is a count for the whole layer it will be
    applied to. select_features_forward() also returns score_holdout:
    the selected set fitted on the selection half and scored on the
    other, the honest number its selection-internal score is not. The
    cost is precision --- half the points estimate, and a contiguous
    spatial half is less efficient than an exchangeable one (García
    Rasines and Young 2023). Data thinning (Neufeld et al. 2024) and
    data fission (Leiner et al. 2023) keep the whole sample and are
    noted on the help page, not implemented. The help page also now says
    plainly that supplying both response_var and predictor_vars upgrades
    the level-selection criterion to "combined", so the selection
    depends on the response without that having been asked for.

  - make_folds() gains blocks: a polygon layer (sf or sfc) to use as the
    blocks of method = "block_kfold" in place of the grid it would
    otherwise build --- the $cells of a build_tessellation() result,
    hexagons, watersheds, administrative units, the $blocks of
    blockCV::cv_spatial(). Each point takes the block that contains it
    and the blocks are assigned to folds exactly as grid cells are, so
    the fold builder can now consume every shape the tessellation half
    of the package produces. The grid-sizing arguments and boundary are
    ignored with a log line, auto_range compares the estimated range
    against the blocks instead of resizing them, and the leakage warning
    uses the median over blocks of the side of the square with the
    block's area (params$block_scale). Points inside no block are
    assigned to the nearest one with a warning that counts them, except
    points within a millionth of the extent of a block --- an edge that
    reprojection moved by a rounding error; points inside more than one
    block take the first, with a warning when the blocks concerned
    overlap in area rather than share an edge. params gains n_blocks
    (before empties were dropped; for a grid, grid_nx * grid_ny, the
    cells a boundary clips away included), blocks_supplied and
    block_scale on every block_kfold result; grid_nx/grid_ny are NA for
    supplied blocks. A blocks layer without a CRS is brought into the
    points' CRS with an R warning naming it, as boundary is.

  - make_folds() gains balance_tol, the largest-to-smallest fold size
    ratio above which block_kfold reports its folds as imbalanced. The
    check was always there at a hard-coded 3:1, and it was a log line
    only; it is now an R warning (Inf disables it), the ratio achieved
    is returned as params$balance_ratio, and the help page says which
    methods balance what: only block_kfold balances point counts, by
    packing blocks largest first into the fold with the fewest points so
    far. No search over packings was added, because the packing is not
    where the imbalance comes from. Measured against the optimum by
    enumeration (two folds, up to ten blocks, heavy-tailed sizes), the
    greedy packing is optimal in 72 percent of cases and within two
    points of optimal on average; a local search with 30 random restarts
    moved the ratio by 0.004 on average over 480 block-size vectors and
    never brought one of the 91 above 3:1 below it. The remedy for an
    imbalance past the tolerance is the block design --- the warning now
    says so --- and blocks is how to supply one: on clustered layouts
    where the geometric grid exceeded 3:1 in 22 percent of draws
    (median 1.8, worst 6.3), Voronoi cells around 15
    get_voronoi_seeds(method = "kmeans") seeds never exceeded 1.3
    (median 1.09). With the default tolerance the folds of every
    existing call are unchanged.

  - metrics on cv_spatial(), cv_gwr(), cv_bayes(), cv_rf() and
    compare_models_cv(): a scoring function of your own, function(y,
    yhat) returning a named numeric vector (a named list of scalars or a
    one-row data frame also serve), applied the way the built-in metrics
    are --- once per fold, so each name becomes a column of
    fold_metrics, and once to the pooled out-of-sample predictions, so
    each name becomes a column of overall --- on the same finite pairs
    RMSE uses. This is the way to score what the Gaussian set cannot: a
    Poisson deviance, a log score on a probability, a weighted loss. The
    contract is strict where it should be (every element named, names
    unique and not a built-in column, one number per name --- anything
    else is an error, because a scoring function of the wrong shape is a
    mistake to surface) and forgiving where it should be (a function
    that throws on a fold is logged and its columns are NA there; a fold
    is never dropped for it). Nor may a name reuse a per-fold extra
    (CRPS, coverage_*, gp_k, gp_n_basis, n_draws, bandwidth, or a name
    the fold_info_fn returns) or mean_CRPS: each used to replace the
    package's value silently, and cv_bayes()$predictive_coverage then
    reported the user's number. The empty frames of a run where every
    fold failed carry the columns, typed, when the function can be
    called on zero-length input. compare_models_cv() hands one function
    to every backend and protects it like the fold arguments, so the
    columns of its overall are comparable across rows. fold_info_fn is
    documented as the per-fold half of the same mechanism, with access
    to the fitted object and the held-out layer.

  - ensure_projected() gains purpose = c("distance", "area"). The
    default is what it always did: for lon/lat input, the candidate that
    distorts distances least. With "area" --- densities or rates per
    cell are going to be computed --- the choice is made among
    equal-area projections only (a Lambert azimuthal centred on the
    data, or an Albers conic where its parallels do not degenerate,
    whichever distorts distances less), which a UTM zone never enters,
    and global coverage gets Equal Earth rather than Web Mercator.
    Already-projected input is still returned untouched, but its area
    distortion over the extent is now measured --- the spread of
    planar-to-geodesic area ratios over probe polygons --- and logged as
    a warning above 1 percent. Measured: a UTM zone edge to edge 0.25
    percent, a 2.5-degree extent inside one 0.04 percent, the
    conterminous United States forced into one zone 14 percent, Web
    Mercator over 2.5 degrees of latitude at 48N 4 percent, an
    equal-area projection a few tenths of a percent (the sphere the
    geodesic areas are computed on against the ellipsoid). The
    measurement works whether sf_use_s2() is on or off (with it off and
    no lwgeom it returned NA, so summarize_by_cell(area = TRUE) refused
    every grid and purpose = "area" skipped its warning), and its probe
    polygons are densified before their geodesic area is taken, so an
    equal-area grid over a near-global extent no longer measures 10
    percent. With purpose = "area", a single polygon's candidates are
    scored rather than falling back to the Lambert azimuthal unmeasured.

  - summarize_by_cell() gains area = TRUE: with cells_sf, the result
    carries cell_area (planar, in the squared units of the cells' CRS;
    for lon/lat cells the geodesic area in square metres, not a planar
    area in squared degrees) and n_per_area, a point density; a rate of
    anything else is its agg_funs sum over cell_area. The request is
    refused with an error --- not answered with a number --- when the
    cells' CRS distorts areas across them by more than 1 percent by the
    measurement above, because a density is a comparison between cells
    and means nothing where the map scale differs from one cell to the
    next; the message names the CRS, the figure and the remedy. A cell
    with no observations gets NA, not zero. Lon/lat cells pass the
    distortion check by construction; with s2 switched off their area
    needs lwgeom, and the request is refused without it. The measured
    spread is attached as attr(, "area_error"). A cells_sf with no ID
    column the summaries can be joined on is an error under area = TRUE,
    since the area columns cannot be produced.

  - build_tessellation(approx_n_cells = ) and get_voronoi_seeds(n = )
    accept what the level-selection step returned: the integer vector of
    ranked candidates from determine_optimal_levels() (its first element
    is used), a select_resolution() result (its $best), or a
    resolution_profile() (read with select_resolution() at its default
    criterion). Both were hard errors before, so nothing that worked
    changes; the count used and where it came from are recorded as
    params$approx_n_cells / params$approx_n_cells_from and as
    attr(seeds, "n_from"). The two functions the pipeline documents as a
    pair are now connected: get_voronoi_seeds(n =
    determine_optimal_levels(pts)) needs no number carried between the
    calls by hand. A hex or square lattice records params$cells_occupied
    and params$cells_empty, and a count taken from a profile or a
    selection warns when fewer than three quarters of it end up
    occupied: the count is of k-means cells, all of them occupied, and
    on clustered points a lattice leaves about half its cells empty.

  - Five diagnostic plots that show the curve behind a chosen point, the
    folds behind a pooled number, or the distribution behind a count.
    None recomputes anything; each draws what the result already
    carries.
    
      - plot_cv_metrics(cv, metric): one point per fold, sized by the
        held-out rows it contributed, with the pooled value from overall
        as a dashed line; a compare_models_cv() result gets one panel
        per model on a shared scale. Any column of fold_metrics can be
        drawn, including backend extras and columns a metrics function
        added; a column that is NA in every fold is refused with the
        reason (Adj_R2 without p, coverage without draws) rather than
        drawn empty, and a per-fold extra with no pooled counterpart
        draws without the line and says so. So does a count (n_pred,
        n_MAPE, n_SMAPE), whose overall value is the total over the
        folds; it was drawn as the pooled line, at 150 against folds
        of 30. A model with no finite per-fold value gets no panel and
        no pooled line, and the caption names it, whether or not overall
        has a value for it (RF's bandwidth was dropped without a
        mention). A compare_models_cv() result keeps its model strip
        when only one model is left to draw.
      - plot.aoa(): the dissimilarity index of the prediction locations
        against the training DI --- cross-validated over the folds
        passed, or else each point's distance to its nearest other
        training point; the legend says which --- as ECDFs, or a
        histogram with the training curve, threshold marked, with the
        share outside and how close the inside ones run to the edge in
        the subtitle, and whether the threshold came from
        cross-validated folds in the caption. Prediction locations
        outside on a predictor dropped for having no training variance
        (DI = Inf) count in the prediction curve, which then tops out
        below 1, and the caption says how many are off the axis. It
        counts the DI = NA rows (a missing predictor) as well. The curve
        used to leave the Inf rows out: it read 0.97 inside at the
        threshold while the subtitle counted 11 of 40 outside. A result
        with every row at DI = Inf was refused as "every row had a
        missing or non-finite predictor", and the error now gives the
        true reason.
      - plot.spatial_fit(type = "variogram") overlays the response's own
        variogram (hollow points, dashed fit) on the residual variogram,
        on the same points and lags, so the structure the model absorbed
        is the gap between the two curves. The response curve is
        whichever variogram estimate_sac_range() returns for the
        response, which is its widest single direction when the
        all-pairs fit is unusable (a response with a trend, whose
        residuals are fine). A single-direction response curve is
        labelled with its azimuth, and the caption compares the sills
        only when both ranges were identified and the two curves cover
        the same point pairs (such a curve used to be labelled plainly
        as the response, its sill set against the all-pairs residual
        curve's: "Residual sill is 45% of the response sill" from a
        quarter of the pairs). response = FALSE restores the residual
        curve alone.
      - plot_calibration(cv): observed against nominal coverage of
        cv_bayes()'s posterior predictive intervals, pooled (blue) and
        per fold (grey), with the diagonal and a one-line verdict. Each
        level is drawn at the nominal value cv_bayes() records in
        coverage_levels (0.995 was drawn at 1.00), so coverage_levels =
        seq(0.1, 0.9, by = 0.1) gives a full curve; the default three
        levels are unchanged. Its error for all-NA coverage names
        compute_pred_intervals = FALSE as well as failed draws.
      - One sweep drawer behind three methods: plot.resolution_profile()
        (a panel per criterion, the level each selects marked, its flat
        region shaded, a note when a bound rather than the criterion is
        choosing); plot.feature_selection() (the accepted variable's
        score at each step as the path, every other candidate faint, the
        stop in red, the best candidate at the step after the stop,
        which was scored and not added, hollow and labelled "not added"
        rather than drawn like the accepted variables, the hold-out
        score as a separate mark when select_on = "split" computed one
        --- and a caption saying whether the intercept-only model was
        scored, since for the RF and GWR backends it usually is not, so
        the path starts at the first variable);
        plot.gwr_model_selection() (every model's AICc against its size,
        the best of each size joined, the winner marked, its lead over
        the runner-up in the subtitle). select_features_forward()'s
        result now carries class "feature_selection" so plot() finds the
        method; it is the same list otherwise.

  - cv_block_size_sweep(): the same cross-validation at a ladder of
    block sizes, with random folds as the leaky reference, returned as a
    table with the fold-to-fold spread at each size and the estimated
    autocorrelation range alongside; plot() draws the curve with the
    range marked. Blocks smaller than the range leak, so the curve rises
    from the random-fold value towards the range and plateaus beyond it,
    and the height of the rise is what the random-fold number overstated
    --- measured on simulated fields with a 100--130-unit range and a
    random forest with coordinates, the plateau begins at one to two
    times the estimated range. Each size is a full cv_spatial(), so a
    fit budget (max_fits, default 60) refuses to start rather than run
    past it. The default ladder runs up to the largest size, at most
    half the shorter side, whose grid still holds k cells. At the
    default k = 5 half the side is a 2 x 2 grid of four cells on any
    extent less than 1.5 times as long as it is wide, so that rung used
    to be dropped every time. n_sizes = 6 then ran five
    cross-validations, the last at 0.30 of the side, and missed ranges
    up to a third of it that a 3 x 3 grid reaches. The top is now about
    a third of the side on such an extent. Sizes the caller passes in
    block_sizes whose grid holds fewer than k cells are not run, and a
    warning names them and the largest size that gives k blocks; they
    were dropped with only a log line. A units object for block_sizes is
    refused by name. On clustered data make_folds() can still lower k at
    the top sizes, which the k column shows. The default ladder also
    skips sizes whose grid would exceed the 1,000,000 blocks
    make_folds() builds (a 6 km x 2 m transect used to abort with an
    error about block_size), runs along the line for points on one
    axis-parallel line (they were refused as having "no extent"), and
    warns when every rung is below the estimated range while longer
    blocks would still fit k times along the longer side (a 10 km x 100
    m corridor: rungs of 4 to 50 m against a 1.7 km range, a flat curve
    that read as no leakage). On a roughly square extent at k = 5 that
    warning also fired, as in the block-size tour script on a 996 m
    square, with two wrong numbers. It called the top rung it had run
    (300) "half the shorter side" (498). It said blocks only up to the
    longer side over k (199, smaller than rungs already run) still gave
    k blocks, when blocks up to 332 did. It now fires only where a
    longer block would still give k blocks, and it names the top rung
    run and the largest size that gives k blocks. User-supplied
    block_sizes over the grid cap are refused before any fit, and a sac
    estimated in another CRS is converted to the sweep's units with a
    warning (a metre range on a US-foot axis was drawn 3.3 times too
    short). The plot's caption reads every coverage_* column as "closer
    to the nominal level is better" (only coverage_50, coverage_80 and
    coverage_95 were known, as higher-is-better, so coverage_97.5 from
    cv_bayes()'s full-precision names was captioned "lower is better").

  - fit_gwr_model() keeps its local collinearity survey. Every fitting
    window --- not a sample of 30 --- has scaled condition indices of
    its kernel-weighted local design computed, the way Wheeler and
    Tiefelsdorf diagnose GWR collinearity, and the fit carries them as
    info$local_collinearity (one row per observation: coordinates,
    window size, and two condition indices: cn, Belsley's uncentred
    index with the intercept, and cn_slopes, the predictors centred in
    the window and scaled by their study-area standard deviation), with
    info$n_local_collinear, info$n_local_singular and the global
    info$condition_index (on the centred predictors) beside AICc.
    n_local_collinear and the warning count the windows whose slopes are
    collinear (cn_slopes above 30, or cn above 1e6), so a predictor's
    origin (degrees C or kelvin) does not change the verdict. The
    warning is now the exact fraction of locations rather than a sampled
    one, and its thresholds (a quarter of the locations, or any) are
    unchanged; above a quarter it now says that an exactly singular
    window stops the fit, instead of promising non-finite coefficients.
    The weighted survey sees what the unweighted spot-check could not: a
    bisquare window's edge points contribute almost nothing to the fit,
    so they contribute almost nothing to its conditioning. The survey
    also runs with a single numeric predictor, and its kernel weights
    equal GWmodel::gw.weight() exactly: a boxcar keeps a point at the
    kernel's edge, and a zero-width adaptive kernel gives NaN, counted
    as singular, instead of weight 1 at the co-located points. (A grid
    with a boxcar bandwidth equal to its spacing used to be reported
    collinear at every location while GWmodel fitted windows of 3 to 5
    points.)

  - plot(fit, type = "coefficients") for a GWR fit maps one local
    coefficient (term) at the training locations, which is the reason to
    fit GWR at all --- and masks the locations where it is not to be
    believed: a collinear local design (for a slope, the slope condition
    index above 30 or singular; for the intercept, the condition index
    with the intercept above 30) or a non-finite coefficient is drawn
    hollow and grey, counted in the subtitle, because the smooth surface
    a naive map draws over them is the picture of an unstable estimate.
    mask = FALSE draws them anyway and says how many it is drawing. A
    fit with a duplicated predictor name is drawn rather than refused as
    "every location is masked". A slope map in kelvin is drawn like the
    one in degrees C; before, every location was masked and the map
    refused.

  - kriging_adequacy(): what a block-kriging aggregator would deliver on
    a set of cells, computed beside the plain means and changing none of
    them. Per cell, from a fitted variogram (estimate_sac_range()'s, or
    estimated here): the block-kriging estimate and variance, that
    variance as a share of the variance the cell's mean would have with
    no data at all (kr_ratio, the coverage score --- near 1 the data
    tell the cell nothing; a share of the point sill never came near 1
    for cells larger than the range), whether it exceeds the
    design-based s^2/n of the plain mean (kr_exceeds_design), and the
    kriged-minus-plain shift in standard errors (kr_shift); plus the
    variance of the standardised errors from blocked cross-validation
    (attr(, "cv")$zscore_var), which is about 1 when the kriging
    variance is right. Measured on simulated exponential
    fields: 0.93--1.07 with the true variogram, 0.85--1.01 with the
    estimated one; under blocked folds it checks the sill and range
    rather than the nugget (0.95--1.24 with the nugget understated
    tenfold), and random folds are the instrument for the nugget. The
    comparison the function exists for: under uniform sampling kriged
    and plain means differed by more than one standard error in 11--24
    percent of cells; under clustered sampling in 34--63 percent,
    with 3--27 of 16--64 cells empty and kriged anyway. Repeat visits to
    one location are kriged from their mean, with the part of the nugget
    that varies between visits divided by their count; kriging them as
    separate rows made every kriging system singular and returned NA
    everywhere. Any cell or held-out location gstat still cannot solve
    is counted in a warning. This is the first kriging path in the
    package (gstat::krige()); its model families are the ones the
    package interprets elsewhere, and any other is refused by name. The
    cross-validation runs each make_folds() split on its own training
    set, so "buffered_loo" and "nndm" keep their exclusion zones
    (reduced to fold labels they ran as plain leave-one-out: with a 250
    m buffer, a standardised-error variance of 1.02 and RMSE 0.757,
    against 0.873 and 0.944 with the buffer kept), and print() names the
    scheme instead of calling every one "blocked CV". A cell holding
    locations that the nmax nearest its centre leave out is kriged from
    all of its own locations plus the nmax nearest outside it (a cell
    of 1,500 points was otherwise estimated from its middle 50: 0.43 off
    against 0.07), and the column kr_n_used says how many locations each
    cell was kriged from. max_neighbours (default 2000) leaves out, with
    a warning, a cell that would need a larger kriging system, and
    max_box_ratio (default 1000) a cell whose bounding box exceeds its
    area that many times, because gstat discretises over the whole box
    (+592 MB for one strip at 7,072); both are counted in attr(,
    "cells_left_out") and by print(). The points are put in the CRS the
    variogram was fitted in (attr(sac, "crs")), as summarize_by_cell()
    does: a sac fitted in metres used on points in km gave a CV
    statistic of 3.05 against 0.85. A sac whose variogram is of
    residuals (detrended = TRUE) is used with a warning that the
    response is kriged without its predictors and the variances come out
    too small (4.3--5.2 against 0.67--1.53). attr(, "rejected_reason")
    records why a sac's range was refused, and the warning and print()
    say it instead of "sill never reached" for every refusal. The cell
    ID is found and matched as summarize_by_cell() finds it (cells keyed
    by id or grid_id too, and a double ID of 1e5 matches an
    integer 100000, where it used to leave that cell with n = 0), a
    point whose ID matches no cell is counted in a warning, and a layer
    with no CRS is taken to be in the other's. Folds whose recorded rows
    sit at other locations in assigned_points_sf are refused, as in
    cv_*(). Folds built before assign_features_to_polygons() dropped
    points used to be applied by position, holding out the wrong points
    (cross-validation RMSE 1.82 against 2.11 with folds built on the
    assigned layer), with nothing said.

  - MAPE and SMAPE now say how many rows they were averaged over. Every
    metrics frame --- model_metrics(), summary(), evaluate_insample(),
    compare_models(), and the overall and fold_metrics of cv_gwr(),
    cv_bayes(), cv_spatial(), cv_rf() and compare_models_cv() --- gains
    two trailing integer columns, n_MAPE and n_SMAPE: the rows each
    percentage error actually used once those where its denominator is
    zero were dropped (y for MAPE, |y| + |yhat| for SMAPE, zero meaning
    no larger than 100 machine epsilons times the data's own magnitude).
    They equal n (n_pred in the CV frames) when nothing was dropped, and
    are 0 in an empty frame. The values themselves are unchanged, apart
    from what now counts as zero (see Bug fixes): a MAPE over 58 of 120
    rows is the same number 2.0.0 reported, but it now arrives labelled,
    where before nothing in the frame recorded that it was a subset
    average. print(summary(fit)) appends "(over k of n rows)" to its
    SMAPE line when the two differ. The columns sit after Adj_R2 so code
    addressing the seven metric columns by position is unaffected; code
    pinning the exact column set needs the two names added.

  - voronoi_seeds_kmeans() gains nstart, the number of stats::kmeans()
    starts (default 10, as before).

  - cv_bayes() and compare_models() say when a Bayesian fit did not
    converge. Both scored such a fit like any other, and only the fit's
    WARN log lines (which tryCatch() and knitr never see, and
    spatialkit_quiet() hides) said that its posterior was not to be
    trusted. cv_bayes()'s fold_metrics gains convergence_ok: TRUE or
    FALSE as fit_bayesian_spatial_model() judged that fold's sampler
    (R-hat, effective sample size, divergences), and NA when fit_args
    sets check_convergence = FALSE. A run with any FALSE raises one
    warning naming those folds. compare_models() gains the same column
    and warns once for each model that did not converge. Code that pins
    the exact column set of either table needs the new name added.

  - select_resolution() returns $seeds, the centres of the partition the
    profile scored at the chosen level, and get_voronoi_seeds(method =
    "kmeans") and voronoi_seeds_kmeans() return those centres when given
    the selection or the profile, instead of running a k-means of their
    own. A k-means partition is the Voronoi partition of its centres, so
    these seeds rebuild the very cells the criteria judged. The fresh
    k-means the seeding functions ran (Hartigan-Wong, 10 starts, on
    whatever cloud they were given) was a different partition: on 300
    points it put 5--16 percent of them in a different cell from the one
    they were scored in, and on 2000 points (profiled on a 1500-point
    subsample) 26--29 percent, so the count was defended on cells nobody
    built. The profile keeps every level's centres in attr(x, "centres")
    (with attr(x, "centres_crs")), the seeds come back in the CRS of the
    boundary or points they are for, and attr(seeds, "kmeans") assigns
    sample_points to them, each point to its nearest seed. A count
    passed as a number (n = sel$best) still runs a fresh k-means, as
    does a selection made before seeds were kept. voronoi_seeds_kmeans()
    did not accept a selection or profile as k before, and the advice in
    build_tessellation()'s occupancy warning now points at the selection
    rather than the bare count.

Bug fixes

  - ensure_stable_poly_id() stopped on a cell that is valid in its own
    CRS and crosses itself in lon/lat, and took
    build_tessellation(method = "voronoi") with it. The sort key is
    measured on a lon/lat copy with s2. Only the vertices are
    transformed, and s2 joins them with great-circle arcs, so a long
    straight edge that passes about a metre from another vertex of the
    same ring can come out on the other side of that vertex. s2 refuses
    such a ring, the repair of the lon/lat copy does not split crossing
    edges, and st_centroid() and st_area() stopped with "Loop 1 is not
    valid: Edge 36 crosses edge 52". Met on 1 of 291 Voronoi cells of
    Texas clipped to the state outline in EPSG:5070, where a 36.9 km
    cell edge passes 1.47 m from a vertex of the coastline and the great
    circle through its ends passes 0.37 m beyond that vertex: no cell of
    the state got an ID. A feature s2 still refuses after the repair now
    gets a second repair of its own sort copy with the crossing edges
    split, and if s2 refuses that too its centroid and area are measured
    as plane geometry in the layer's own CRS and the centroid is
    transformed to the sort CRS. The function warns once, saying how
    many features took each route, and create_voronoi_polygons() and
    build_tessellation(method = "voronoi") pass the warning on. Those
    features' keys are close to the spherical ones and not equal to
    them: the second repair moved the Texas cell's area by about 40
    square metres in 2,949 square kilometres, and its centroid taken in
    the plane of EPSG:5070 lay 20 m from the one the second repair gave
    on the sphere. Every other feature keeps exactly the key it had, so
    the other 290 Texas cells are numbered in the order they always
    would have been, the geometry returned is still the caller's own,
    and the same 291 cells get the same IDs in EPSG:5070, EPSG:3083, UTM
    zone 14N, Web Mercator and lon/lat. A feature keyed in the plane is
    measured in the CRS the layer arrived in (the Texas cell's plane
    centroid lies 9 to 73 m from the spherical one depending on that
    CRS), so it can take another place in the order when the same layer
    arrives in another projection, and the IDs between the two places
    move with it; the warning says so whenever that route is taken. With
    make_valid = FALSE no repair is tried, as the caller asked, and a
    feature s2 refuses goes straight to the plane: such a call used to
    stop with s2's error and now returns IDs with the warning. All three
    methods are covered: for "surface_point" and "bbox_center" the point
    never went to s2, but the area, the third key of each, did, and both
    stopped there. An error that is not s2 refusing a feature is raised
    as it was.

  - resolution_profile() took the noise variance of C_p from the nugget
    of a variogram fit that had not converged, whenever the range was
    refused for another reason. Its documentation says a model that did
    not converge gives cp no nugget, the nugget and sill being where the
    optimiser stopped, and the profile decided that by looking for
    "converge" in rejected_reason. estimate_sac_range() gives a refused
    range one reason, and a falling variogram and a range past the
    largest lag come before the optimiser in it, so a fit that stopped
    at gstat's iteration limit with a range past the fitted lags is
    recorded as "fitted range exceeds the largest lag fitted". That is
    most fits that do not converge: over 200 simulated layers (trends,
    ranges longer than the extent, 40--60 points, white noise) 96 of
    the 102 all-pairs fits that did not converge carried that reason
    and 6 "variogram model did not converge". With an east-west trend on
    an exponential field (300 points, nugget 0.1, 25 draws) every
    all-pairs fit stopped at the limit, every one was refused as past
    the lags, and C_p was given nuggets of 1.18--1.62; in the first draw
    the range reported was 29,180 against a largest lag of 647. Whether
    the model converged is now read from the model (attr(attr(sac,
    "variogram_model"), "converged"), which estimate_sac_range() already
    set); a sac that carries no such flag (made by hand, or built from
    REML parameters) is judged by its reason, as before, and a reason
    that says the fit did not converge is believed whatever the flag
    says. For a refused range whose model did not converge, cp takes its
    noise variance from the finest level whose cells hold at least two
    scored rows on average, attr(x, "cp_noise")$source reads
    "finest-level residual mean square", attr(x, "variogram") is NULL,
    and the warning names both facts: why the range was refused, and
    that the model did not converge. The existing notice that names the
    level and the variance cp was given follows it, so such a profile
    raises two warnings where it raised one. The finest-level value can
    be larger or smaller than the nugget it replaces (larger in 41 of 52
    simulated layers, smaller in 11), so the cp pick can move either
    way: it stayed in 49, moved to fewer cells in 2 and to more in 1.
    When no level holds two scored rows a cell, cp is NA at every level
    for such a sac, where it was finite; the warning's last clause says
    so, and the notice is not raised. A fit that converged to a range
    past the lags, or on a falling variogram, still gives cp its nugget.
    estimate_sac_range() returns what it did; ?estimate_sac_range now
    says that the count of "variogram model did not converge" reasons
    understates how many fits did not converge, and where to read it.

  - A clipped cell that also touched the boundary from outside was
    dropped, and its points left with no cell. st_intersection() returns
    such a cell as a GEOMETRYCOLLECTION of its area and a line or point
    (on an L-shaped boundary, the cell over the inner corner), and the
    Voronoi, grid and Delaunay builders kept only POLYGON and
    MULTIPOLYGON rows, so the whole cell went, area and all. On a 4 x 4
    L with square cells of side 2, ten of thirty points had no cell and
    the cells covered 8 of the boundary's 10 square units. Each
    collection is now reduced to its polygonal part, one row per cell
    still; a piece with no area is dropped as before.

  - coerce_to_points() crashed R on an empty line feature. sf's
    st_cast() turns an empty MULTILINESTRING into one empty LINESTRING,
    not zero parts, and st_line_sample() on it segfaulted and took the
    session with it. A null geometry in a line layer loads exactly like
    this from a GeoPackage or a shapefile, and prep_model_data(), every
    cv_*() function and fold_separation() go through this path. GEOS's
    st_point_on_surface() crashed the same way on a line feature holding
    an empty part beside real ones, which ensure_projected() reached on
    ordinary lon/lat input. Empty parts are now removed before either
    call, and an empty feature becomes an empty POINT in its own row,
    which prep_model_data() and make_folds() drop like any empty
    geometry. An empty LINESTRING used to raise an error; it now gives
    an empty POINT like every other geometry type.

  - cv_gwr(parallel = n) hung forever once a GWR had been fitted in the
    session. GWmodel is built with OpenMP, and GNU libgomp is not
    fork-safe: after fit_gwr_model(), a sequential cv_gwr() or a bare
    GWmodel::bw.gwr(), the forked mclapply() workers blocked on a futex
    and never returned, with no timeout. That is the ordinary
    fit-then-validate order on Linux. cv_gwr() now runs its folds one
    after another whenever parallel would fork, and says so in a
    warning. This gives up the speed-up parallel folds had in a fresh
    session; the other cv_*() functions still fork, since ranger does
    not use libgomp. A cv_spatial() fit_fn that calls GWmodel can still
    hang with parallel, which its help page now says.

  - When some folds failed, overall quietly left them out. A partial
    failure was only logged, so overall pooled the surviving folds with
    no R condition, and the folds that fail are usually the hardest to
    predict (a region or a factor level no training fold covers).
    cv_gwr(), cv_bayes(), cv_spatial() and cv_rf() now warn, naming each
    failed fold and how many rows overall covers. Folds dropped before
    fitting are already warned about and are not counted twice.

  - select_features_forward() could pick a variable for making folds
    fail. Each candidate was scored on whatever rows its CV run
    predicted, so a set whose fit failed on a fold (a factor level found
    in one block, an ordinary case under block folds) was scored on
    fewer, easier rows. A pure-noise factor beat the true driver:
    RMSE 2.35 on 192 rows against 2.63 on 250, although the driver
    scored 1.73 on those same 192 rows. Every candidate is now scored on
    one fixed row set (the rows the null model predicts, or, with no
    null model, the rows the step-1 sets predict), a set that leaves any
    of them unpredicted scores NA with a warning, and history gains
    n_pred.

  - compare_models_cv() ranked models scored on different rows. Shared
    folds guarantee the same splits, not the same scored rows: when one
    backend lost a fold or returned NA predictions its overall row
    pooled a subset, and the help page called the columns comparable. A
    fixed-bandwidth GWR scored on 158 of 200 rows ranked above a random
    forest scored on all 200, although the forest was 45 percent better
    on the 158 rows both predicted. When the predicted row sets differ,
    every model is now rescored on the rows all of them predicted, with
    a warning, overall carries n_pred, and the all-rows numbers stay in
    attr(overall, "all_rows"). The per-fold table and the Bayesian
    coverage columns are not rescored.

  - determine_optimal_levels() read an elbow into points that have none.
    The elbow was the level furthest below the chord of the WSS curve on
    linear axes. For points with no cluster structure WSS falls like c /
    k, and the furthest point below that chord is exactly sqrt(a * b)
    for a ladder from a to b, so the answer was set by the ladder's
    ends: 1,500 uniform points gave 4, 4, 7, 9 and 13 for max_levels
    of 12, 20, 40, 80 and 160. At the default max_levels = 12 it also
    missed well-separated clusters (four clusters came back as 3). The
    elbow is now read on log-log axes, where c / k is a straight line,
    and counts only when the curve sags clearly below it. Two, three,
    four, five, eight and ten well-separated clusters are now recovered
    at the default on every seed tried (six came back as five on one
    seed in five). With no elbow the function still returns its
    linear-axis answer, but warns that the ladder chose it;
    build_tessellation() and get_voronoi_seeds() refuse a geometry-only
    resolution_profile() with no elbow instead of drawing a count from
    it. A ladder of two levels (max_levels of 1 or 2, or three points)
    has no line to test. The warning now says the ladder is too short to
    read an elbow from and names what ended it. It used to say the curve
    fell in a straight line "as it does for points with no cluster
    structure", which two clusters 90 m apart at max_levels = 2 were
    told.

  - One invalid polygon changed the assignment rule for the whole layer.
    assign_features_to_polygons() wrapped st_join(largest = TRUE) in a
    retry without largest on any error. The comment blamed a predicate
    that cannot take largest, but sf never calls the predicate on that
    path; the retry fired when GEOS threw on an invalid ring, and every
    straddling feature was then assigned by tie_break instead: 95 of 200
    buffered parcels changed cell, and a polygon 91 percent inside one
    cell went to its neighbour, with no warning. Invalid features and
    cells are now repaired with sf::st_make_valid() for the join only,
    with a warning counting them, and a join that still fails stops and
    names largest = FALSE.

  - Lon/lat polygons were assigned to projected cells on bent cell
    edges. assign_features_to_polygons() moved the cells into the
    features' CRS, so lon/lat features pulled the package's own
    projected cells into lon/lat, and with s2 the largest-overlap join
    failed and fell into the silent retry above: 55 of North
    Carolina's 100 counties went to a cell other than their largest
    overlap on a 36-cell grid. The join now runs in the cells' CRS
    whenever it is projected, and the features come back with the
    geometry they arrived with. Lon/lat points near a cell edge can
    change cell as a result; they now agree with a join done in the
    projected CRS.

  - build_tessellation(method = "triangles") dropped most points at UTM
    coordinates. qhull lifts each point onto x^2 + y^2, and at projected
    magnitudes (a northing near 5e6) the lift had no precision left to
    separate points a few metres apart, so they never became
    vertices: 200 points over 100 m gave 26 triangles instead of 386,
    and the help page's own example lost 8 of its 20 points. The points
    are now centred before triangulation. Triangles that were already
    right are the same triangles, but qhull returns them in a different
    order, so triangle cell_id values change.

  - Data around a pole were projected to Web Mercator. A layer spanning
    more than 180 degrees of longitude with no gap was treated as global
    coverage. Antarctic stations came out with worst-case distance
    errors near 20,000 percent (the South Pole at y = -2.4e8 m), and
    Voronoi cells put 9.5 percent of Arctic locations in a station's
    cell that was not their nearest. A layer that lies wholly on one
    side of the equator now gets a Lambert azimuthal equal-area
    projection centred on its pole whenever that measures a smaller
    distance error than the global fallback: about 2 percent on the same
    stations.

  - GWR mixed elevation into its distances. prep_model_data() kept the Z
    (and M) coordinate of POINT Z input, which GPS layers and
    st_as_sf(coords = c("x", "y", "z")) produce. predict.gwr_fit() then
    handed GWmodel three coordinate columns, which it reshaped into two,
    scrambling the prediction locations with no warning;
    gwr_model_selection() ranked models on 3-D distances; and
    fit_gwr_model() and cv_gwr() failed. Z and M are now dropped in
    prep_model_data() and again where the data are handed to GWmodel.

  - predict.gwr_fit() lost every prediction to one bad location. It went
    through GWmodel::gwr.predict(), which returns nothing for any row if
    one location's window is empty or singular, so fixed-bandwidth block
    CV failed every fold. It also never assigns its distance matrix once
    training and new rows together exceed 10,000 (a 100 by 100
    prediction grid came back all NA), and it built the full training
    hat matrix for a variance it then discarded, which is cubic in the
    training size. Predictions now come from gwr.basic(regression.points
    = ) in chunks, a failing chunk is redone location by location, and
    only a truly singular location is NA, with a warning counting them.
    Where the old path worked, the values are identical.

  - On brms 2.17 to 2.22, one far-off row moved every Bayesian
    prediction in the call. brms rebuilt the Hilbert-space GP's boundary
    from the rows being predicted, and the package's padding rows could
    only widen it, so a single row outside the training envelope changed
    the basis for every row, interior ones included: a 40 by 40 grid
    padded 15 percent moved the in-bbox cells by a mean of 22 percent of
    the surface's standard deviation, and predict_surface() depended on
    chunk_size. On those versions the GP term's boundary factor is now
    rescaled per call so the boundary stays at its fitted value, and a
    row beyond the boundary, where the basis means nothing, is NA with a
    warning. brms 2.23 stores the boundary itself; there only the NA
    rule applies. The help page no longer says the boundary "has to
    grow".

  - A bayesian_fit's cached fitted values could come from another model.
    The cache lives in an environment, so it is shared by every copy of
    a fit, and the entry was stamped with the row count and a digest of
    the training data only. Two fits over the same data hash identically
    however different their engines are, so refit <- fit; refit$engine
    <- <re-estimated> made fitted(), residuals(), model_metrics() and
    summary() on either object return the other's numbers --- and
    clear_fitted_cache(), which the help page offers for exactly this
    case, could not fix it, because clearing through one copy cleared
    the one shared entry and the next call re-wrote it. The entry now
    carries the environment rstan and brms create for each sampling run
    (@.MISC on the stanfit), which tells engines apart as well as the
    engine itself does and is written once, and is used only for the
    engine it came from (identical(), which settles the common case by
    pointer, so nothing is slower and no memory is held that the fit did
    not already hold, in a session or on disk). A stale entry is no
    longer deleted on a miss either: the fit that wrote it still wants
    it. new_spatial_fit() now always builds a fresh cache rather than
    adopting one that arrived in info, and summary() no longer carries
    .cache at all --- a summary was not a value snapshot (its contents
    changed when anyone later called fitted() on the fit),
    clear_fitted_cache() on a summary emptied the fit's cache, and
    saveRDS() on one serialised an environment holding the full
    n-vector.

  - Every CRPS was NA above 46,340 posterior draws. .crps_energy()
    formed its weights with m * m, where m is nrow(draws) and therefore
    an integer, so the product overflowed and took CRPS, mean_CRPS and
    the whole predictive_coverage summary to NA behind one "NAs produced
    by integer overflow" warning. 48,000 draws is an ordinary
    cv_bayes(fit_args = list(chains = 4, iter = 13000)) run. The
    arithmetic is now done in double precision; the values it produces
    were, and remain, exact against the closed-form Gaussian CRPS.

  - Fold numbers depended on the machine's collation. Character fold
    labels are turned into a factor, whose levels as.factor() sorts
    under LC_COLLATE, and a fold's number is its level's position.
    Labels differing only in case or punctuation --- "north" and "North"
    --- therefore landed in a different order under C than under en_US,
    so fold_metrics$fold, predictions$fold and fold_status$fold named
    different groups on different machines from the same data and the
    same seed. The partition was never affected, so pooled scores were
    right. Character levels are now sorted with method = "radix", which
    is always C collation, making the numbering a property of the labels
    alone. The C order applies to character labels only: numeric labels
    are numbered in numeric order, as in 2.0.0, and a factor by its own
    level order (an earlier development build sorted numbers as strings
    too, so with ten or more numeric labels fold 2 was the user's
    label 10). area_of_applicability() given a vector of fold labels now
    numbers them by the same rule; it sorted character labels with
    as.factor() under the session's collation, so a message could name a
    different fold from the one cv_*() named for the same label. Its
    partition and threshold were never affected.

  - residual_morans_i(k = ) silently answered a different question. k
    reached the weight builder unvalidated, where min() collapses a
    vector: k = c(4, 8) built the k = 4 matrix and returned a statistic
    for neighbours the caller never asked for, with no condition raised,
    and k = NA aborted on "missing value where TRUE/FALSE needed". k is
    now checked. Two neighbouring gaps are closed with it: a
    user-supplied weights matrix holding any NA, NaN or Inf is refused
    by name instead of aborting inside a guard after a log line reading
    "row sums range from NA to NA", and a row whose geometry is empty is
    now dropped with a count, as make_folds() and estimate_sac_range()
    already do --- its residual is perfectly finite, so it used to
    survive into FNN::get.knn() and abort with "Data include NAs".

  - A scalar argument that reaches as.integer() is validated. Seven
    exported functions took a count or a distance and passed it straight
    into as.integer() or an if () test, where NA, Inf, a length-2 vector
    or any value above .Machine$integer.max aborts with "missing value
    where TRUE/FALSE needed" or "'length = 2' in coercion to
    'logical(1)'" --- errors that name nothing the caller passed.
    fit_gwr_model(bandwidth = ) under adaptive = TRUE,
    predict_surface(chunk_size = ) (Inf being the natural way to ask for
    one chunk), voronoi_seeds_kmeans(k = ), voronoi_seeds_random(k = )
    and the count resolved by build_tessellation() / get_voronoi_seeds()
    now report the argument and the bound. voronoi_seeds_kmeans(k = 0)
    and a negative k used to return one seed in silence, which "at most
    k" does not describe.

  - expand was ignored rather than refused. clip_target_for(expand =
    c(0.05, 0.05)) returned a clip target byte-identical to expand = 0
    with no condition raised, and create_voronoi_polygons() did the same
    for a non-numeric expand, because both tested is.numeric(expand) &&
    expand > 0 and short-circuited to FALSE. A malformed expand is now
    an error naming the argument; the same values on a degenerate
    bounding box used to abort inside the expansion helper instead.

  - summarize_by_cell(deff = "variogram") aborted when given no value
    column. With neither response_var nor predictor_vars the internal
    primary column is NULL by design, and df[[NULL]] raised "attempt to
    select less than one element in get1index". The variogram design
    effect is a function of the cell's coordinates and the fitted
    correlation, not of any column's values, so every point in the cell
    now counts and the call returns its counts and cell_weight as
    documented.

  - compare_models() aborted on a list holding no spatial_fit.
    evaluate_insample() warns and skips a non-fit, and returned NULL
    when every element was skipped (it is now an error, below); the NULL
    then became a bare list and seq_len(nrow(NULL)) raised "argument
    must be coercible to non-negative integer". It now says which
    argument is wrong and what belongs there.

  - fit_rf_model(include_coords = TRUE)'s caveat said "once per session"
    and was not. .log_warn_once() records the key in a package-level
    environment, and under cv_rf(parallel = ) that write happens inside
    a forked worker and dies with it: the paragraph printed once per
    worker, the parent's registry stayed empty, and the next sequential
    fit printed it again. cv_rf() now raises it in the parent before
    dispatching any fold, so each worker inherits the already-warned
    flag and stays quiet. Measured on two workers: three occurrences
    before, one after.

  - The grid cache's order registry outlived the environments it
    described. Insertion order is kept in a package-level environment
    keyed by the cache environment's printed address, and
    clear_grid_cache() removes only the entry for the environment it is
    handed. Passing a fresh cache_env per call therefore left one
    permanent character vector per call, for environments that had since
    been garbage-collected and could no longer be named. A finalizer now
    removes an entry with its environment.

  - plot(fit, type = "variogram") no longer runs its subtitle off the
    edge of the figure. ggplot2 clips a label that is wider than the
    plot instead of wrapping it, and the sentence saying why no range
    was identified is up to 150 characters, so on a six-inch figure it
    was cut mid-word. Labels built from a fit's own numbers are now
    wrapped at draw time.

  - ensure_stable_poly_id() could not give IDs to a tessellation this
    package had just built. It repaired the geometry in the layer's own
    CRS and then transformed it to the sort CRS, but validity is a
    property of the geometry in the CRS it is measured in: two vertices
    a centimetre apart in a projected CRS can land on one longitude and
    latitude, and s2 calls the ring degenerate. On a clipped hex
    tessellation of North Carolina, 2 of 18 cells that are valid
    projected are invalid once transformed, and st_centroid() on one of
    them aborted the call with "Loop 0 is not valid: Edge 1 is
    degenerate". The sort copy is now repaired after the transform as
    well. The geometry returned is still the caller's own, and the same
    cell gets the same ID whether the layer arrives projected, in
    lon/lat or in Web Mercator, except where two fine cells' centres lie
    within the sort key's rounding step of the same longitude (36
    of 2,500 100 m cells straddling a UTM central meridian changed ID
    via EPSG:3035); join such layers on geometry.

  - select_features_forward() now says when fit_fn is ignoring the
    variables it is handed. The learner it takes is a function of
    (train_sf, predictor_vars), but the one cv_spatial() takes is a
    function of train_sf alone, and a learner written for that ---
    function(train_sf, ...) --- swallows the second argument and fits
    the same model every time. Nothing errored, because the training
    layer still carries every column: each candidate scored exactly what
    the intercept-only model scored, no candidate improved on it, and
    the result was an empty selected and an NA score with no word about
    why. Two different predictor sets do not produce the same
    cross-validated metric to the last digit, so that pattern is now
    recognised at the first step and reported as a warning naming the
    fix. The result is still returned.

  - create_grid_polygons() warns when both cellsize and target_cells are
    supplied. It already did for cellsize and n, and the documentation
    says to supply exactly one of the three, but target_cells was
    dropped in silence when cellsize was present. Through
    build_tessellation() that meant method = "hex", approx_n_cells = 25,
    cellsize = 10 returned however many cells a 10-unit lattice holds
    and said nothing about the 25. cellsize still wins; the override is
    now logged like its sibling.

  - build_tessellation() warns when approx_n_cells or cellsize is
    supplied with method = "voronoi" or "triangles". Both arguments size
    the hex and square lattices and nothing else --- Voronoi grows one
    cell per input point and Delaunay one triangle per neighbouring
    triple --- and both used to be dropped in silence. So
    build_tessellation(pts, method = "voronoi", approx_n_cells = 25)
    returned one cell per observation: the degenerate nearest-neighbour
    case, where every cell holds a single point, there is no within-cell
    variation and every standard error is NA. The warning names the
    argument and points at get_voronoi_seeds(), which is where a Voronoi
    cell count is actually set. It adds that params does not record the
    request either, so a saved result carries no sign of it: the Voronoi
    branch returns create_voronoi_polygons()'s own list, which has no
    slot for the argument, and the triangles branch no longer echoes
    approx_n_cells back (below). The warning fires under quiet = TRUE,
    which gates this function's message()s and is documented not to
    silence R warnings.

  - determine_optimal_levels() fits each k as the best of 25 k-means++
    restarts (Arthur and Vassilvitskii 2007; Fränti and Sieranoja 2019;
    Steinley 2003) instead of stats::kmeans(nstart = 5). The WSS curve
    is read for its shape, and with a handful of random restarts it
    carried optimisation noise: on eight-cluster layouts a sweep over k
    = 1..30 rose at one or two steps in three of five draws, and an
    earlier form of the elbow rule once selected such a bump. With the
    new budget the same sweeps rose at no step. A curve that still rises
    is now logged as a warning naming the number of rising steps, and
    the model-aware diagnostics carry it as wss_bumps beside wss_spread
    (the relative spread of WSS across restarts at each k) and nstart. A
    selection made on a curve that had a bump can differ from before ---
    those were the cases that were wrong; a clean curve gives the same
    answer. Under a model-aware criterion the function also now warns
    before the sweep when max_levels leaves no k above the nine-cell
    floor, rather than fitting every k first and falling back
    afterwards. The seeding draws each centre by inverting the
    cumulative squared distance rather than with sample.int(prob = ):
    the same law, and a resolution_profile() of 3000 points takes 10 s
    instead of 35 s (determine_optimal_levels() on 3000 points, 4 s
    instead of 13 s), but it gives different centres for the same seed
    than earlier development builds did.

  - make_folds(method = "block_kfold") can now raise its "block
    dimension < autocorrelation range" warning. The comparison was
    always there, but the range it compared against was estimated only
    under auto_range = TRUE, the one setting in which the blocks had
    already been sized from that range and the warning could never fire;
    on every default call the diagnostic was dead code. With auto_range
    = FALSE (the default) and a response_var to hand --- always the case
    when a cv_*() function builds the folds --- a range is now estimated
    for the diagnostic alone. It sizes nothing: the blocks are the same
    geometric blocks as before and the folds do not change. A hand-set
    block_size below the range raises the same warning. The estimate's
    own log lines stay off the console, and the check is skipped (with
    an INFO log line saying so) when gstat is not installed or there are
    fewer than 30 points. Only a grid dimension that is split is
    compared with the range, because a single row (or column) of blocks
    borders no other block across its width. On points along a line the
    comparison was with 0, so the warning fired on every call that had a
    response: a 10 km line had 15 blocks 667 m long against a 499 m
    range. Its advice, block_size = 499, made the blocks shorter (19
    of 523 m). A 10 km x 100 m corridor was compared with its 100 m
    width in the same way. A response_var that names no column of
    points_sf is now an error for block_kfold, whether or not auto_range
    is set. With auto_range off, a misspelt name used to switch the
    check off silently.

  - fit_rf_model(include_coords = TRUE) logs its caution once per
    session rather than once per fit. Inside a five-fold cv_rf() or a
    twenty-fit cv_block_size_sweep() the same paragraph printed on every
    fit, which reads as twenty problems rather than one decision; the
    message now says it will not repeat.

  - fit_rf_model() reports what ranger actually objected to. ranger
    diagnoses a bad argument in its C++ layer, writes the diagnosis
    straight to stderr and then throws "User interrupt or internal
    error." --- so mtry = 99 on a two-predictor forest printed "mtry can
    not be larger than number of variables in data. Ranger will EXIT
    now." to the console and raised an error naming neither the argument
    nor the problem. That line is not an R condition, so
    suppressMessages(), withCallingHandlers() and tryCatch() all missed
    it: it escaped every handler to the console (and into CI logs, where
    it reads as an error from a test that is passing), while cv_*()
    recorded the placeholder as the fold's cause. The message stream is
    now diverted for the duration of the call, so the diagnosis becomes
    the reported reason --- fold_status$message included --- and nothing
    is printed behind the caller's back. Output from a call that
    succeeds is passed through unchanged, and when the stream is already
    diverted (under testthat, knitr or capture.output(type = "message"),
    where only one sink is permitted) the call runs exactly as before.
    So does a call made after the session temp directory has been
    deleted, when there is nowhere to divert the stream to; it used to
    fail with "cannot open the connection".

  - One empty point made determine_optimal_levels() return 1. The
    function reduced every feature to a point and projected it but never
    dropped an empty or non-finite one, so a single POINT EMPTY made the
    k = 1 WSS NA, k-means failed at every k, and the failure handler
    shrank the sweep to nothing: two well-separated clusters that gave
    2 1 3 gave 1 once one empty row was added, with only a log line
    about interpolating the WSS. Such rows are now dropped with a
    warning that gives their number, and under select_on = "split" the
    returned positions still index the layer as passed.

  - A few missing predictor values took Moran's I away from
    determine_optimal_levels(). A cell mean over a row with one missing
    predictor was NA and dropped the whole cell, and how many cells that
    removed depended on the points per cell, so z was computed on a
    different subset of cells at each k. Three missing values in 400
    rows made every z in the evaluated window NA, and the call fell back
    to the geometric ranking with a warning that did not mention missing
    values. Rows with a missing or non-finite response or predictor now
    stay in the WSS sweep and the cells and are left out of Moran's I,
    with a logged count, so every cell mean uses the same rows.

  - A misspelt response_var or predictor_vars in
    determine_optimal_levels() silently changed the criterion. A column
    that was not there counted as "no model variables": with the default
    criterion a typo kept "geometric" where supplying both variables
    upgrades to "combined" (on one 400-point layer 7 6 8 instead of
    11 10 7, with no warning and no diagnostics), and under "morans_i"
    or "combined" the one log line said the variables were required
    although both had been supplied. A named column that is not in the
    layer is now an error, as it is in resolution_profile(), and a
    model-aware criterion given no predictor_vars (or no response_var)
    falls back with an R warning that says which is missing, rather than
    a log line.

  - determine_optimal_levels() stopped one level short when locations
    repeat. The sweep was capped one short of the number of distinct
    locations, the bound stats::kmeans() needs only when no location
    repeats: five stations visited thirty times each could not reach k
    = 5. The cap is now the number of distinct locations, and still one
    short of the number of points. A k whose WSS is 0 (to within 1e-12
    of the total: a cell on every location) is left out of the log-log
    elbow line, and when the rest of the curve has no elbow, the fall to
    zero is the elbow. Leaving the level out and reading the rest made
    the call warn that the five stations had no cluster structure and
    return 3 2 4 (one metre of jitter gave 5 4 6); it now returns 5 4.
    Two stations visited thirty times each give 2 1, not 1 2. Zero is
    relative because k-means leaves floating-point residue: two groups
    of ten stations visited ten times each have a WSS of 7.8e-17 at k
    = 20, which dragged the whole line down and, at max_levels = 30,
    reported no cluster structure (still answering 2).
    resolution_profile() reads its elbow column the same way. The
    no-elbow warning names the bound that ended the ladder (the distinct
    locations, the points or max_levels); it named max_levels whichever
    bound it was.

  - The same layer with its rows in another order gave
    determine_optimal_levels() another answer. The subsample and every
    k-means start index rows, so a permutation of the input moved the
    WSS curve and could move the count. The rows are now put in
    coordinate order (response and predictors breaking ties) before
    either, so any permutation of a layer gives the same result. Results
    for a given seed differ from 2.0.0's for that reason too.

  - determine_optimal_levels() blamed the wrong cause when its
    model-aware criteria had nothing to score, and its help page
    understated the fix. The model-aware pass scores only the elbow's
    neighbourhood, and on points with no cluster structure the elbow
    sits near sqrt(max_levels), so every candidate stayed at or below
    the nine-cell floor until max_levels was about 40 (measured on 1000
    uniform points: 12, 20 and 30 all fell back, and 40 scored k = 10
    and 11 alone). The help page said "above roughly 10", and the
    fallback was logged as "Moran's I could not be computed". The logged
    warning now names the window and the floor and points to
    resolution_profile(), and the help page gives the measured numbers.
    Which candidates are scored is unchanged.

  - build_tessellation() could not tessellate CRS-less planar points
    inside a CRS-less boundary. ensure_projected() marks such points
    crs_assumed = "none" ("planar, leave alone"), and
    build_tessellation() read that mark as a CRS name: st_crs("none")
    failed with "invalid crs: none" for every method, including the
    documented boundary = clip_target_for(pts), so CRS-less planar data
    could not be gridded at all. Only a real assumption (EPSG:4326) is
    now given to the boundary, which is refused if its coordinates
    cannot be degrees (below); with no assumption both stay in the same
    unnamed space.

  - When only one of the points and the boundary had a CRS, the
    tessellation builders stopped on sf's bare "st_crs(x) == st_crs(y)
    is not TRUE". UTM points read from a CSV with a UTM boundary failed
    in all four methods of build_tessellation() and in
    create_voronoi_polygons(); projected points with a boundary that had
    lost its .prj failed in create_voronoi_polygons() and method =
    "triangles", and clip_target_for() returned a target with no CRS ---
    for lon/lat points with a CRS-less lon/lat boundary, in degrees,
    with expand = 20 buffering by 20 degrees. The side without a CRS is
    now interpreted in the other's, as harmonize_crs() does, with a
    warning: lon/lat-looking coordinates are reprojected from EPSG:4326,
    others are stamped. CRS-less points that do not look like lon/lat
    are refused, with a message saying what to do, when the boundary is
    geographic, because stamping degrees on them would be wrong; a
    CRS-less boundary beside geographic points is read as lon/lat or
    refused (below).

  - A boundary without a CRS got a log line in make_folds(), cv_*() and
    predict_surface(), where every other function raises an R warning,
    and cv_rf() warned twice about one lon/lat boundary. With projected
    points, build_tessellation(), clip_target_for(), plot_folds() and
    the rest said "boundary has no CRS ... stamping" as an R warning,
    while make_folds(), cv_*() and predict_surface() stamped it with a
    "WARN ensure_projected(): input has no CRS" log line that tryCatch()
    and knitr never see and spatialkit_quiet() hides, naming neither the
    function nor the argument. A boundary whose coordinates looked like
    lon/lat got two R warnings from one cv_rf() call, one from
    prep_model_data() and one from make_folds(), both naming
    ensure_projected(). All three now warn as the others do, naming
    themselves and boundary (make_folds() also prediction_points,
    predict_surface() also grid and covariates), and a cv_*() call warns
    once, naming the cv_*() function. Which CRS the layer ends up in is
    unchanged. For predict_surface() this matters most on grid and
    covariates: a layer stamped with the wrong CRS puts every covariate
    lookup in the wrong place. The stamping warning of every function
    now names the CRS it stamps ("stamping the target CRS ('EPSG:32632')
    WITHOUT reprojection") instead of "the supplied crs", an argument
    most of them do not have.

  - A geographic crs made every tessellation method work in degrees.
    build_tessellation(), create_voronoi_polygons() and
    create_grid_polygons() took crs = 4326 as the CRS to compute in, so
    Voronoi cells stopped being a nearest-point partition (150 points at
    53-57N: 20 percent of sampled locations lay in another point's
    cell), grid cells were neither square nor equal-area, and clipping
    them under s2 stopped with "Edge 0 is degenerate" (hex) or left a
    point inside the boundary with an NA index (square). A geographic
    crs is now the CRS the result is returned in: the cells are built
    and indexed in the local projected CRS ensure_projected() picks,
    then transformed with long edges densified. A hex or square grid
    sized by an explicit cellsize is still laid in degrees, since that
    is the unit cellsize is in.

  - ensure_projected() stopped on lon/lat polygons that s2 rejects, and
    chose its UTM zone differently when sf_use_s2() was off. The centre
    that places the zone was st_centroid(st_union()): with s2 on, a
    polygon with a repeated vertex (valid for GEOS, common in
    shapefiles) stopped ensure_projected(), create_grid_polygons() and
    prep_model_data(boundary =) with "Edge 1 is degenerate (duplicate
    vertex)"; with s2 off it was planar in degrees, so two clusters
    at 0N and 60N near -78 got UTM zone 17 in one session and zone 18 in
    another, and sf's warning and message about it got past quiet =
    TRUE. The centre is now always taken on the sphere, a geometry s2
    rejects is repaired first, and a mean of unit vectors is the last
    resort.

  - A single study-area polygon was never scored, so ensure_projected()
    kept a UTM zone at any extent. A polygon layer is scored on one
    point per feature, and one point makes no pair: every candidate
    scored NA and the selector fell back to the zone. A CONUS outline
    stayed in UTM zone 15 (12.8 percent worst-case distance error) where
    a Lambert azimuthal scores 2.1 percent, and prep_model_data(boundary
    =) moved the whole analysis into the zone with it. Layers with fewer
    than 40 features are now scored on their outline's vertices too.

  - create_grid_polygons() laid a grid over a near-global lon/lat
    boundary in Web Mercator. The cells were equal on the map and not on
    the ground: the true areas of whole cells differed nearly five-fold.
    When the CRS picked for distances distorts areas across the boundary
    by more than 1 percent, the grid is now laid in the equal-area CRS
    ensure_projected(purpose = "area") picks (Equal Earth here), with a
    logged warning; a local extent keeps its UTM zone.
    create_grid_polygons_cached() makes the same choice.

  - With sf_use_s2(FALSE) the package needed lwgeom, which it does not
    depend on. The stable-ID sort key is measured in lon/lat, so every
    Voronoi tessellation --- projected ones included, and the examples
    of create_voronoi_polygons() and build_tessellation() --- failed
    with "package lwgeom required", as did ensure_stable_poly_id(),
    create_grid_polygons_cached() on a cache miss, and random and
    k-means seeding on a lon/lat boundary. These measurements now run on
    the sphere (s2) whatever the session's setting, which is restored
    afterwards; the IDs are the ones s2-on sessions always got.

  - Hex cell size depended on which way the boundary lay. st_make_grid()
    builds hexagons from cellsize[1] alone, and that was the box's width
    over a rounded column count: a 1 x 1000 strip at target_cells = 9
    got 1734 hexagons where the same strip lying flat got 89. The size
    is now counted along the longer side, which leaves every boundary at
    least as wide as it is tall with exactly the grid it had.

  - Voronoi with expand > 0 returned the boundary before it was grown.
    The cells are clipped to the grown boundary, so they covered 1.93
    km^2 against a returned boundary of 1 km^2, and points up to expand
    outside it were indexed. create_voronoi_polygons() and
    build_tessellation(method = "voronoi") now return the grown
    boundary, and the documentation says expand grows the study area
    too.

  - Collinear points gave an empty Delaunay triangulation.
    build_tessellation(method = "triangles") on a transect returned no
    cells and an index of NAs, after logging that delaunayn() had
    failed, which it had not. It now stops with a message that says the
    points are collinear and points to method = "voronoi"; the
    fallback's log line names the reason that applies.

  - k-means seeding clustered CRS-less lon/lat points as if degrees were
    metres. voronoi_seeds_kmeans() and get_voronoi_seeds(method =
    "kmeans") projected only when a CRS said lon/lat, so 800 CRS-less
    points 89 km wide and 111 km tall were split east-west where the
    same points tagged EPSG:4326 were split north-south. They now apply
    the lon/lat heuristic ensure_projected() applies (with its warning)
    and return the seeds in the input's own coordinates.

  - voronoi_seeds_random() returned the same seeding on every call. Its
    default set_seed = 456 reset the random-number stream inside the
    call, so five calls under set.seed(1) to set.seed(5) gave one draw,
    and the sensitivity comparison its help page recommends compared a
    seeding with itself. The default is now NULL, as in
    get_voronoi_seeds(): the draw comes from the session's stream. Pass
    set_seed for a fixed seeding. This changes the default.

  - harmonize_crs() refused a layer as target_crs. A multi-row sf failed
    with "the condition has length > 1" and a one-row one with "cannot
    create a crs from an object of class sf", where ensure_projected()
    accepts both. An sf or sfc target now means its CRS.

  - summarize_by_cell() failed on deff = NA and mis-recorded deff = Inf.
    The check on a numeric deff came out NA for NA_real_ and NaN, so the
    call stopped with "missing value where TRUE/FALSE needed" instead of
    the documented warning and fallback to 1. Inf passed the check: the
    standard errors were the uncorrected ones, yet cell_weight was 0 in
    every cell and deff_applied recorded deff = Inf. Any deff that is
    not a single finite number of at least 1 now falls back to 1 with
    the warning.

  - summarize_by_cell(deff = "variogram") fell back to uncorrected
    standard errors with no R warning when no model could be fitted.
    With predictors only and no sac, without gstat, or when
    estimate_sac_range() returned no fit (fewer than 30 points), the
    only signal was a log line: no tryCatch(), withCallingHandlers() or
    warnings() saw it, and spatialkit_quiet() hid it, while the standard
    errors came out about 5 times smaller than the corrected ones in one
    check. It is now a warning that names the reason. Every
    design-effect fallback (a refused deff, a rejected or unsupported
    variogram, no model) is raised with class
    "spatialkit_deff_fallback", so a loop over many summaries can catch
    exactly that case. A rejected sac that a variogram estimated from
    response_var replaces is not a fallback: it gets a plain warning,
    and the classed warning is raised only when nothing replaces it,
    once per call, naming every reason. A pure-nugget model (no
    structured component) implies that distinct observations are
    uncorrelated, so it is applied as a design effect of 1 in every
    cell, with deff_applied = TRUE; it was reported as "the supplied
    model could not be read", with the fallback warning.

  - An empty point switched off summarize_by_cell(deff = "variogram")
    for its cell. Its missing coordinates made the cell's mean
    correlation NA, and the cell silently got a design effect of 1: the
    standard error of that cell was a quarter of the corrected one in
    one check. (In the development version the same input stopped the
    call with "missing value where TRUE/FALSE needed".) Points with
    empty or non-finite coordinates now count towards their cell's
    values but not towards its correlation, with a warning.

  - summarize_by_cell(deff = "variogram") read an anisotropic variogram
    as isotropic. Only model, psill and range were read, so vgm(0.8,
    "Exp", 300, 0.2, anis = c(0, 0.2)) gave a correlation of 0.677 at 50
    m east-west where gstat's is 0.348, and a median design effect of
    14.1 against 7.7: standard errors too wide by a factor of 1.7. A 2-D
    geometric anisotropy (ang1, anis1) is now applied as gstat applies
    it. (resolution_profile(), which uses the same correlation function
    on distances alone, still reads the major range.)

  - A misspelt or non-numeric response_var was dropped in silence.
    summarize_by_cell() reported it through a progress message, which
    the default quiet = TRUE suppresses, and returned a frame with no
    resp_* column and cell_weight equal to n. A missing or non-numeric
    response_var or predictor is now a warning, and a response_var of
    more than one name is an error instead of "the condition has length
    > 1".

  - A double cell ID of 100000 lost its cell in
    summarize_by_cell(cells_sf = ). When the points' and the cells' ID
    columns had different classes (an integer poly_id against a double
    from a GeoPackage Integer64 field or a CSV), both were converted
    with as.character(), which writes 1e+05 for the double and 100000
    for the integer. That cell came back with NA summaries and its
    points left the result: 11 of 30 points in one check, with no R
    warning. Whole numbers are now written out in full, and a summarised
    ID that matches no cell is reported with a warning.

  - summarize_by_cell(cells_sf = ) did not read the ID columns
    assign_features_to_polygons() writes from. Cells keyed by id or
    grid_id (common in shapefiles) were assigned cleanly and then
    summarised to a plain table with no geometry, cell_area or
    n_per_area, and a layer with id and a differently numbered cell_id
    was joined on cell_id, putting 89 percent of the summaries on the
    wrong polygons in one check. The cells are now searched in the order
    the assignment used, and a cells_sf that cannot be joined is a
    warning (an error with area = TRUE) rather than a log line. agg_funs
    = median or "median" is honoured instead of being replaced by the
    mean. A single function is named after the expression passed, so
    stats::median gives resp_median_* as median does (it gave
    resp_agg1_*).

  - assign_features_to_polygons(largest = TRUE) assigned polygon
    features that only touch the cells. sf keeps the largest
    intersection piece without checking its area, so under GEOS a
    feature sharing only an edge or a corner with the cell layer went to
    that cell with zero overlap, while under s2 (lon/lat) the same
    feature was unassigned: sf's nc counties against 50 of them as cells
    gave 70 rows projected and 50 in lon/lat. A feature with no overlap
    area is now unassigned in every CRS.

  - tie_break = "smallest_area" depended on the row order after all.
    Candidates of equal area --- the cells of any regular grid, for a
    point on a shared edge --- fell through to the first row, so
    reversing the cells' rows moved every edge point to the neighbouring
    cell (points at x = 100 went to cells 1, 4 and 7, or to 2, 5 and 8),
    and on a cached grid, whose IDs do not run row by row, one cell
    could take both of its edges. Equal areas (to 9 digits) are now
    decided by the lowest, then leftmost, bounding-box centre. On a
    create_grid_polygons() square grid that is the cell row order
    already picked (1,681 of 1,681 lattice points unchanged); on a hex
    grid 3 of 56 shared vertices move.

  - create_grid_polygons_cached() could return another site's grid. The
    cache key used the CRS's input name, which is "unknown" for any
    custom CRS read back from a GeoPackage or shapefile, so two
    site-centred CRSs with the same local boundary coordinates shared an
    entry, and the second site received the first one's grid 11,000 km
    away. The key now hashes the CRS's WKT; the only cost is a rebuild
    when one CRS arrives written two ways. target_cells now defaults to
    NULL, as in create_grid_polygons(), so cellsize = or n = work
    without it.

  - ensure_stable_poly_id() could number cells differently with s2 off.
    The sort key's centroid was taken planar in degrees when
    sf::sf_use_s2(FALSE), which moves it by far more than the key's
    rounding step, so near-tied cells swapped IDs between s2-on and
    s2-off sessions (4 of 2,000 Voronoi cells), and without lwgeom the
    s2-off call stopped at the area. The key is now taken on the sphere
    for the sort copy only, whatever the session setting.

  - summarize_by_cell(deff = "variogram") built each cell's correlation
    matrix up to five times per column. The same call with two partly
    missing columns and conf_level now builds 15 where it built 45.

  - estimate_sac_range() returned the widest direction that reached a
    sill when the all-pairs variogram had run past the fitted lags.
    Under an unremoved trend the pooled variogram rises without a sill,
    which the help page said the fitted-lag bound catches; but when two
    of the four directional variograms (those across the slope) did
    reach one, their maximum came back as the range with anisotropy_used
    = TRUE and nothing on the console. On an exponential field of
    range 150 with an east-west trend, 12 of 30 draws returned 131--596
    this way and 16 others NA, so make_folds(auto_range = TRUE) switched
    between a 2 x 2 grid of 430 m and geometric blocks from one draw to
    the next. The directions that reach a sill are the shorter ones, so
    their maximum is a lower bound, not an estimate. A converged
    all-pairs fit past the fitted lags is now refused whatever the
    directions found (rejected_reason = "fitted range exceeds the
    largest lag fitted", the directional ranges still attached); the
    directional maximum stands in only for an all-pairs fit that is
    singular or did not converge. On stationary fields whose range is
    close to the cutoff this also turns a few draws from a directional
    maximum into NA (3 of 30 at an effective range of 570 on a 1000 m
    square).

  - print() on an estimate_sac_range() result showed a bare number. For
    lon/lat input the number is in metres of a CRS the estimate picked,
    and with predictor_vars it is the range of the residuals rather than
    of the response; neither showed, so a mismatch with the layer it was
    about to be used on could not be seen. A last line now names the
    unit and the CRS (in metres of EPSG:32617) and whether the variogram
    is of the response or of its residuals (and by which detrend
    method). For a layer with no CRS the line says the range is in that
    layer's own coordinate units (in the coordinate units of a layer
    with no CRS), still with what was modelled; it used to be left out.

  - estimate_sac_range(predictor_vars = ) fitted its variogram to the
    raw response, trend and all, when a single predictor or response
    value was infinite. lm()'s na.exclude drops NA but not Inf, so one
    Inf among 200 rows stopped the detrending ("NA/NaN/Inf in 'x'"), and
    the function fell back to the raw response with a warning. The range
    came out at 2950 instead of 2168 (the answer with that row removed),
    and make_folds(auto_range = TRUE, predictor_vars = ) built blocks 37
    percent larger. resolution_profile(), which hides that warning, then
    warned that the variogram it had estimated itself was of the raw
    response. With detrend = "reml", a -Inf predictor (a log(0)
    covariate) first gave a false warning that the REML fit "did not
    converge", and then the OLS fallback failed the same way. Rows with
    a missing or non-finite response or predictor are now left out of
    the detrending fit and the variogram, with a logged count, as
    prep_model_data() and resolution_profile() already do. One Inf now
    gives the same range as that row set to NA (2167.5 detrended on the
    example; 1598.4 under REML).

  - make_folds(drop_empty_blocks = FALSE) could return folds with no
    test points. k was lowered only when the highest block id holding a
    point was below it, and with empty blocks kept that id says nothing
    about how many blocks hold points: two clusters on a 4 x 4 grid gave
    id 16, so k = 5 was kept for 2 occupied blocks and three folds came
    back empty, with no warning (the imbalance check skipped an empty
    fold). The cv_*() functions then ran on the two folds that had
    points and area_of_applicability(folds =) failed. k is now lowered
    to the number of blocks that hold points, as ?make_folds already
    said, the log line names the numbers, and every point in a single
    one of several blocks is the single-block error it always was with
    the default.

  - The automatic block grid gave an east-west corridor more than twice
    as many blocks as the same corridor running north-south. Only the
    row count was capped, so a layer more than about block_multiplier *
    k times as wide as it is tall got round(sqrt(15 * w/h)) columns at k
    = 5: 39 x 1 on a 10 km x 100 m corridor against 1 x 15 turned on its
    side, blocks less than half as long, and a scheme drifting towards
    random k-fold (1-NN CV RMSE 0.77 against 0.94 on the same values;
    lower in 16 of 20 seeds). Points on one horizontal line were treated
    as a square and got a 4 x 4 grid that collapsed onto the line,
    lowering k from 5 to 4. The column count is now capped at
    block_multiplier * k as the row count was, so both orientations and
    both lines get 15 blocks at k = 5. Folds change only for extents
    more than about block_multiplier * k + 1 times as wide as they are
    tall.

  - make_folds() ignored block_nx or block_ny given alone, and accepted
    invalid ones. Giving one dimension sent the call to the automatic
    grid without a word (block_nx = 10 alone gave a 3 x 4 grid); 0, a
    negative, NA or a vector failed inside sf or base R, and 2.7 was
    truncated to 2. The dimension given is now used and the other
    derived from the extent's aspect ratio (roughly square blocks), and
    each must be a single whole number >= 1.

  - With a non-rectangular boundary, make_folds(block_kfold) kept
    zero-area slivers as blocks. Where the boundary only touches a grid
    cell at a corner or along an edge, the clipped cell is a POINT or a
    LINESTRING, and it was kept as a block: packed into a fold under
    drop_empty_blocks = FALSE (3 of 13 blocks under a triangular
    boundary), and able to catch a data point lying exactly on the
    boundary as a one-point block of its own. Only the areal part of the
    grid is kept now.

  - make_folds() failed with R's own errors on a missing k or an invalid
    buffer. k = NULL (or no k) reached if (k < 2) and failed with
    "argument is of length zero"; a buffer of NA, numeric(0) or length 2
    failed with "missing value where TRUE/FALSE needed" or a length
    error, and an NA is exactly what estimate_sac_range() returns when
    no range is identified. Both are now refused by name (the
    leave-one-out methods still need no k), and so is a units object
    passed as buffer, block_size or block_nx/block_ny, which used to
    fail inside the units package without naming the argument.

  - make_folds(method = "buffered_loo") said nothing when the buffer
    excluded no neighbour. The buffer is in the units of the CRS the
    folds are built in, which for lon/lat input is metres, so a buffer
    in degrees (0.1) excluded nothing and the scheme was plain
    leave-one-out: 79 of 79 training points in every fold of an 80-point
    layer, no condition raised. It now warns when no fold excludes any
    neighbour, and ?make_folds says what unit buffer is in.

  - NNDM folds could be far more optimistic than their target with only
    a log-file line to say so. When min_train stops the matching ---
    samples clustered well inside the prediction domain, the layout NNDM
    is meant for --- the realised distances stay short: one cluster
    predicted onto a 20 km grid kept a median of 1171 m against a target
    of 6704 m, 96 of 100 folds held at the floor, while ?make_folds said
    the result is never more optimistic than the target. make_folds()
    now warns when the floor leaves more than one point's worth of
    excess below phi, records params$n_at_min_train, and the
    documentation states the guarantee only where neither phi nor
    min_train binds. The warning names the statistic that fires it, the
    largest excess of the realised nearest-neighbour ECDF over the
    target at distances up to phi, and gives the two medians only as
    context: they summarise all the distances, and the realised median
    can sit above the target's while the short distances are
    over-represented.

  - make_folds(auto_range = TRUE) fell back to geometric blocks with
    only a log line. When no range was identified (an unremoved trend, a
    range past the fitted lags, fewer than 30 points, gstat missing),
    the blocks the caller asked to be sized from the data were not, and
    under knitr, spatialkit_quiet or tryCatch() nothing showed it. This
    is now a warning that gives the rejection reason.
    estimate_sac_range() can give up before it fits anything: fewer
    than 30 points, fewer than 30 finite values, a constant response (or
    residuals, when the predictors explain the response exactly), points
    with no extent, or gstat missing. It used to return a bare NA then,
    with the reason only in a log line, and the warning said only
    "estimate_sac_range() returned NA". That NA now carries a
    rejected_reason attribute saying which (it is still unclassed, with
    no other attribute). The warning quotes it, as do
    kriging_adequacy()'s no-model error and summarize_by_cell()'s
    fallback warning.

  - NNDM fold construction releases FNN's copy of the neighbour tables
    as soon as it has them, so a second n x n/2 pair is no longer held
    through the sweep and the construction of the folds. The peak inside
    get.knn() is unchanged.

  - cv_bayes() rounded coverage levels to a whole percent, so close
    levels overwrote each other. Coverage columns were named
    sprintf("coverage_%.0f", 100 * level): coverage_levels =
    c(0.5, 0.975, 0.985, 0.995) gave three columns for four levels,
    coverage_98 holding the 0.985 value and the 0.975 value lost, with
    no condition raised. Levels given as percentages, c(50, 80, 95),
    made every fold throw away its CRPS, n_draws and coverage in
    silence. Columns are now named at full precision (coverage_97.5; the
    default 50/80/95 names are unchanged), a level outside (0, 1) or
    given twice is an error that suggests dividing by 100 where that
    fits, and the result carries coverage_levels, the nominal level of
    each column.

  - An error in a fold_info_fn threw away all of that fold's extras
    without a word. cv_spatial() caught it and dropped the whole list,
    so the columns were NA with fold_status "ok" and nothing logged; in
    cv_bayes() one failing quantile cost gp_k, n_draws, CRPS and every
    coverage column. The failure is now logged and named in
    fold_status$message, and cv_bayes() computes coverage and CRPS in a
    step of their own, so gp_k, n_draws and yhat_sd survive it. A
    fold_info_fn that returns the wrong shape (a vector where a value
    belongs, an unnamed or duplicated element, or a name fold_metrics
    already has, such as RMSE) is an error that says so; RMSE = -5 used
    to overwrite the fold's real RMSE.

  - cv_spatial() failed after fitting every fold when ..per_row came
    back from some folds only. The prediction rows were stacked with
    rbind(), which died on "numbers of columns of arguments do not
    match" once every fold had been fitted, so the work was lost and the
    message did not point at the cause. A fold without ..per_row
    (returned conditionally, of the wrong length, or from a fold_info_fn
    that threw) now gets NA in those columns, and one of the wrong
    length is logged.

  - Parallel cross-validation discarded folds that shared a core with a
    failure. mclapply() ran prescheduled, handing each core a chunk of
    folds: an error that escaped one fold was copied to every fold of
    its chunk, and a worker killed for lack of memory took all of its
    folds with it. With four folds on two cores, a failure on fold 2
    lost fold 4 as well, and overall pooled 40 of 80 rows. Each fold now
    runs in its own worker, so a failure costs that fold only, and an
    error that stops a sequential run (a metrics or fold_info_fn return
    value of the wrong shape) stops a parallel one too, naming the fold.
    A run in which nothing fails gives the same numbers as before, since
    the per-fold seeds are drawn before forking.

  - model_metrics(newdata = ) measured R-squared against a different
    baseline from every cv_*() function. It took the total sum of
    squares about the new rows' own mean, where cross-validation takes
    it about the training mean, so the same predictions on a split
    across a trend scored R-squared -0.89 here and 0.35 from
    cv_spatial(). model_metrics(), evaluate_insample() and
    compare_models() with newdata now use the training mean too, the
    out-of-sample convention, and the help pages say which baseline
    R-squared uses. In-sample numbers do not change.

  - R-squared and MAPE depended on the units of the response. The
    thresholds below which a total sum of squares or a percentage-error
    denominator counted as zero were absolute, so a response with
    standard deviation below about 1.5e-8 got R-squared NA while RMSE
    and MAE were fine (and select_features_forward(metric = "R2")
    selected nothing), and one on a 1e-15 scale lost MAPE and SMAPE as
    well. "Zero" is now 100 machine epsilons of the data's own magnitude
    for every metric: a rescaled response gets the same R-squared and
    MAPE, a constant one still gets NA, and on a response spanning many
    orders of magnitude a row whose denominator is no larger than 100
    epsilons times the largest (such as 1e-9 against 1e6) no longer
    enters MAPE.

  - compare_models() set out-of-bag random-forest metrics beside
    in-sample ones without saying so. Without newdata an rf_fit's fitted
    values are out-of-bag and every other backend's are in-sample, and
    the table could rank the models the wrong way round: GWR RMSE 0.77
    in-sample against RF 0.82 out-of-bag, where the forest's in-sample
    RMSE was 0.40. evaluate_insample() and compare_models() now carry a
    metric_basis column ("in-sample", "out-of-bag" or "newdata"),
    compare_models() logs a note when a table mixes them, and both help
    pages say what the metrics are computed on.

  - compare_models() put LOOIC and AICc side by side for fits on
    different rows. Both are sums over the rows a model was fitted to,
    so a model that lost 20 rows to a predictor's missing values showed
    LOOIC 32.2 against 54.8 for the model on all 70, and looked 22.6
    better while it was worse on the rows they share. A column whose
    models were fitted to different rows is now set to NA, with a
    warning naming each model's n. The same happens to fits of the same
    rows with different responses (a response and its log, say), since
    an information criterion compares models of one response only, and
    the warning now says so. It used to say they were "fitted to
    different rows (raw: n = 80, logged: n = 80)" and to "Refit them on
    the same rows".

  - compare_models() read significantly negative residual
    autocorrelation as missed spatial structure. The caution fired on a
    two-sided p-value whatever the sign, so the alternating in-sample
    residuals of a GP or a small-bandwidth GWR (Moran's I -0.13, p
    = 0.02) were logged as "may not fully capture the spatial
    structure", the opposite diagnosis. A negative z is now logged as
    what it usually means, a model tracking its data closely.

  - compare_models_cv() placed polygon rows in its shared blocks by a
    different point than every model is fitted at. The shared folds
    reduced polygons and lines to their point-on-surface whatever
    pointize said, while each backend fitted them at the pointize point:
    with pointize = "centroid", 119 of 150 L-shaped parcels were in a
    different fold from a standalone cv_gwr() run. The shared blocks now
    use pointize; with the default "auto" nothing changes.

  - Saved folds on lon/lat polygons were refused after sf_use_s2() was
    toggled. The provenance check located each probed row by its
    centroid, which on a geographic CRS is spherical with s2 on and
    planar with it off; the two differ by up to 5e-4 degrees on county
    polygons, 500 times the tolerance. Folds built before
    sf_use_s2(FALSE) (a common workaround for invalid polygons), or
    saved and read in a session set the other way, were rejected by
    every cv_*() as "built from different data". The probe now always
    takes the planar centroid; folds saved by an older version are
    checked the way they were made.

  - residual_morans_i() gave the wrong reason when it could not use a
    fit's residuals. An error from residuals() was thrown away, and it
    and a fit with no residuals() method (the ?new_spatial_fit example
    has none) were both reported as "could not extract enough residuals
    (n < 4)", on a 100-row fit; a residual vector of the wrong length
    was reported as "coordinate extraction failed". Each now has its own
    warning, quoting the error where there is one, except that a fit
    with no residuals() method is now scored on the response minus
    fitted() instead (below).

  - A GWR whose local regressions interpolate the data won on AICc.
    GWmodel's AICc is defined only while the effective number of
    parameters, tr(S), is below n - 2; past that its penalty turns
    negative. A small adaptive bandwidth reached it, and so did any
    bandwidth raised to the old floor, which for the bisquare and
    tricube kernels (they give the farthest neighbour in a window
    weight 0) fitted every window exactly. On 100 points
    compare_models() listed R2 = 1 and an AICc of -15033 against 292 for
    the automatic bandwidth, and at n = 60 gwr_model_selection()
    selected the real predictor plus three noise variables with -66689.
    fit_gwr_model() now reports such an AICc as NA and
    gwr_model_selection() ranks such models last, each with a warning
    giving tr(S). The adaptive floor is one neighbour higher for
    bisquare and tricube, and raising a supplied bandwidth to it is now
    a warning, not a log line. bandwidth = NULL was not affected
    above 20 points. The raised floor is enough unless several
    neighbours tie at the kernel's edge (a regular grid); the warning
    now says so. Where an adaptive bandwidth is already every
    observation, the undefined-AICc warning suggests fewer predictors,
    more observations or a gaussian or exponential kernel instead of a
    larger bandwidth.

  - fit_gwr_model() never checked a one-predictor model for local
    collinearity. The check ran only with two or more numeric
    predictors, but every local design includes the intercept, and a
    predictor nearly constant inside a window is collinear with it. A
    regional covariate nearly constant within each of four clusters gave
    local slopes from -97 to 221 around a true 3 with no warning, while
    adding a noise predictor to the same data warned at every location.
    One numeric predictor is now enough.

  - GWR said a singular window came back as NaN coefficients; it stops
    the fit. GWmodel's matrix inverse throws on an exactly singular
    window, so a 0/1 indicator constant within clusters failed the whole
    fit with a bare "inv(): matrix is singular", while the help page and
    the collinearity warning promised masked NaN coefficients. The
    non-finite coefficients GWmodel does return come from co-located
    points, where an adaptive bandwidth no larger than the number of
    observations at a site gives the kernel zero width (160 of 160 at 40
    sites of 4 observations, 4 neighbours), and the warning blamed
    singular windows for those. The fit error now says a window is
    singular and how many the collinearity check found, and the
    non-finite warning names co-located points when they are the cause.

  - An adaptive GWR bandwidth above the number of points was capped in
    silence. bandwidth = 1500, meant as metres with adaptive left at
    TRUE, became a 200-neighbour, near-global fit on 200 points without
    a word; its local slopes varied less than half as much as the
    intended fixed-distance fit's. fit_gwr_model() and
    gwr_model_selection() now warn, naming n and pointing to adaptive =
    FALSE. cv_gwr() repeats the warning in each fold whose training set
    is smaller than the bandwidth.

  - Below 20 points, bandwidth = NULL did not fit at the bandwidth
    bw.gwr() chose. GWmodel searches adaptive bandwidths from 20
    neighbours up to n, so with fewer points its choice exceeds n (18
    for 12 points) and was capped at n, a different kernel with a worse
    AICc (18.3 against 9.7), without a word. The cap stays and now
    raises a warning in fit_gwr_model() and gwr_model_selection();
    supply bandwidth for data this small.

  - A GWR predictor named twice raised false collinearity warnings.
    fit_gwr_model(predictor_vars = c("a", "b", "a")) fitted correctly,
    but its collinearity checks ran on the duplicated column and warned
    "exactly singular" and "100% of locations collinear", once in every
    cv_gwr() fold, and the doubled count raised the bandwidth floor.
    Names are now collapsed on entry, as gwr_model_selection() already
    did.

  - A gp_c you set did not change the gp_k derived for it. The basis
    count was always sized for the boundary factor the package would
    have chosen, so a wider boundary got the same number of basis
    functions and could no longer resolve the lower length-scale bound
    it was sized for. The advice under gp_c is to raise it for a
    long-range surface, which is exactly the case that coarsened the
    basis. On 200 uniform points, gp_c = 3 fitted with gp_k = 23 where
    the rule gives 43, and gp_c = 5 with 23 where the rule gives 70
    (capped at 50); the cap warning could never fire on this path. With
    gp_k = NULL the derived gp_k is now sized for the gp_c actually
    used, and a capped value is logged. An explicit gp_k still passes
    through untouched.

  - fit_bayesian_spatial_model() could not fit brms::categorical() or a
    mixture() family. Those families give each distributional parameter
    its own GP under the same coefficient names, and the automatic
    length-scale prior kept only the names, so every coefficient got two
    identical rows and brms stopped with "Duplicated prior
    specifications are not allowed" before sampling. A user's lscale
    prior restricted to one dpar failed the same way, because it was
    copied onto every category's coefficients. The prior now carries
    each coefficient's dpar, nlpar and resp, and a global or dpar-level
    lscale prior is expanded only onto the coefficients it addresses and
    never over a coefficient-level one the user already gave. With
    standardize_predictors = TRUE they still failed, on the automatic
    normal(0, 5) slope prior, which carried no dpar and so matched no
    slope of either family (brms: "The following priors do not
    correspond to any model parameter: b ~ normal(0, 5)", a prior the
    user never wrote). That prior is now set on each distributional
    parameter's slopes, as the length-scale prior is; a family with one
    mu gets the same single row as before.

  - A two-level factor response under brms::bernoulli() fitted, and then
    nothing could score it. The response check refused a non-numeric
    response only under gaussian, and the gaussian refusal itself
    pointed at bernoulli(). brms fits the factor, but residuals() came
    back all NA, summary(), model_metrics() and compare_models() stopped
    on "response is factor", and cv_bayes() ran a full fold of MCMC
    before aborting in the fold scoring (with parallel = 2, every fold
    ended as worker_error). A factor or character response is now
    refused, before anything is compiled, under every family except
    categorical() and the ordinal ones (cumulative, sratio, cratio,
    acat), with a message saying to convert it to 0/1. Numeric and
    logical 0/1 responses are unaffected.

  - predict() on an ordinal or categorical bayesian_fit returned all NA
    as a "posterior draw failed". brms returns posterior_epred() for
    those families as a draws x rows x categories array, which the
    method took for a failed draw: a real cumulative() fit returned NA
    for all five new rows, with only a log line, while fitted(),
    summary() and model_metrics() said merely that they got an array.
    predict() under its default type = "epred" and fitted() now stop,
    saying the family has a probability per category and pointing at
    type = "predict", draws = TRUE, whose share of draws in each
    category estimates its probability for any rows, and at
    brms::posterior_epred(<fit>$engine) for the training rows
    (posterior_epred(<fit>$engine, newdata = ), which the message used
    to suggest, refuses new rows without the scaled coordinates the
    method builds). The message now counts the caller's rows, where it
    counted the two GP-boundary rows as well ("150 x 7 x 3" for five
    rows). type = "predict" without draws = TRUE on a
    brms::categorical() fit, which returned the mean of unordered
    category indices (1.46, 1.97, ...), is now an error; for an ordinal
    family it is the expected category index, as documented. A genuinely
    failed draw still returns NA as documented, and the log line now
    carries the cause.

  - predict() on an rf_fit turned every ranger error into an all-NA
    vector. type = "quantiles" on a forest grown without quantreg =
    TRUE, and type = "se" without keep.inbag = TRUE, which the help page
    says are rejected, returned ten NAs for ten rows with no R
    condition, and model_metrics(newdata =, type = "se") then reported n
    = 0. A failure in ranger's predict method is now an error naming
    ranger's reason; the cv_*() fold loop records it as the fold's cause
    and predict_surface() stops naming the rows, as they already did for
    other backends. A newdata with no complete row still returns all NA
    with a log line, as for the other backends, rather than reaching
    ranger as a zero-row frame; that had made a predict_surface() chunk
    outside the covariates' coverage abort the whole surface.

  - check_convergence = FALSE returned convergence_ok = TRUE. The flag
    started out TRUE, so a fit whose checks never ran (its max R-hat
    was 1.28) claimed to have passed them over an empty diagnostics
    list, and print() had nothing to caveat. It is now NA when nothing
    was checked, and print() on the fit and on its summary() says
    "Convergence: NOT CHECKED"; summary()'s printout also repeats the
    "Convergence warnings present" flag, which it carried and never
    showed. A failed PSIS-LOO is now logged with its cause instead of
    "LOO computation failed." alone, which had left compare_models()
    showing LOOIC NA with nothing saying why.

  - The convergence check raised dozens of "The ESS has been capped"
    warnings. brms::neff_ratio() runs posterior's ESS over every GP
    basis weight, and posterior warns once per well-mixed one: an n = 80
    fit raised 44 R warnings, 41 of them this one. R keeps only the
    first 50 warnings, so a warning that mattered and came later, loo's
    Pareto-k among them, could be dropped. That one message is now
    muffled around the R-hat and ESS accessors; every other warning
    passes through, and the ratios are unchanged.

  - A saved rf_fit or bayesian_fit carried its engine twice. The model
    formula was built in the fitting function's frame and so captured
    it, and that frame holds the forest or the brmsfit itself; a formula
    serialises its environment, so saveRDS() wrote the engine a second
    time (1.62 MB for a 100-tree forest of 0.72 MB; about 80 MB for a 40
    MB brmsfit, which brms's own copy of the formula doubled even in
    saveRDS(fit$engine)). The formulas now carry the global environment,
    as a formula typed at the console does.

  - A forest with rows out of every tree's bag said nothing. ranger
    returns NaN as the out-of-bag prediction of a row every tree
    sampled, so with num_trees = 5 20 of 200 rows had NaN fitted values
    and summary() printed "n = 200" over an R-squared computed on 180.
    fit_rf_model() now warns with the count, and summary() prints
    "(computed on 180 of 200 rows ...)" when its metrics use fewer rows
    than the fit has. cv_rf() does not use its fold forests' out-of-bag
    predictions, so it warns once per run with the number of fold
    forests affected, instead of once per fold (each of which told the
    user to score the forest with cv_rf()).

  - area_of_applicability() counted rows that differ on a dropped
    zero-variance predictor as inside the AOA. A predictor constant in
    the training data --- a land-cover dummy absent from the training
    region --- is dropped from the distance, so new rows taking another
    value there were judged on the other predictors alone: 37 of 40
    urban rows came out inside, where a single urban training row would
    have kept the predictor and put 1 inside. Such rows now get DI = Inf
    (their scaled distance along that predictor is infinite), are
    counted outside, and a warning gives the count; print() says how
    many.

  - area_of_applicability() returned a threshold of 0 from duplicated
    training rows without saying why. Each training row's reference is
    its nearest other row, so exact duplicates in predictor space ---
    repeat visits to a site with static covariates, covariates from a
    raster coarser than the sampling --- have a training DI of 0; past
    about three quarters of the rows the threshold is 0 and only exact
    copies count as inside (30 sites visited four times: 0 of 200 new
    points inside, against 197 after deduplication). The rule is
    unchanged; the result is now logged with the remedy
    (leave-location-out folds, or deduplication) and print() shows the
    count of zero training DI.

  - area_of_applicability() refused all-zero weights, which the advice
    pmax(importance, 0) produces whenever the model found no useful
    predictor. A one-predictor forest gave all-zero weights in 9 fits
    of 20 when its predictor carried no signal, so the AOA was lost in
    the folds where it mattered most. With one predictor the weight
    cannot change the index (it is scale-invariant) and zero is accepted
    silently; with several, all are weighted equally, as weights = NULL
    would, with a warning.

  - A fractional chunk_size in area_of_applicability() marked
    extrapolation as inside the AOA. On the dense path (use_fnn = FALSE,
    or FNN not installed) block starts became fractional and the rows
    between blocks kept an initial DI of 0: with chunk_size = 2.5, 5
    of 25 far-out points were reported inside and 16 training DI of 0
    moved the threshold. chunk_size is now validated by name and
    truncated to whole rows.

  - area_of_applicability(model = fit, folds = folds) stopped with
    "fold 1 refers to rows outside 1:n" whenever prep_model_data() had
    dropped a row, for the same folds cv_*() accepted. One missing
    response in 200 rows was enough: make_folds() numbers the rows of
    the layer it is given, the fit keeps only the 199 rows
    prep_model_data() returned, and the fold IDs were read as positions
    in those. The documented workflow ("pass the same make_folds()
    result you passed to cv_spatial()") therefore failed on any layer
    with a missing or non-finite modelling value or an empty geometry.
    The rows the fit's "dropped" record names are now taken out of the
    folds, as cv_*() take them out, with a log line giving the count; a
    label vector with one label per row of the layer fitted from loses
    those labels. The threshold is the one you get by removing the rows
    from the folds by hand. Folds built on the model's own training data
    (the prep_model_data() output) are still read as positions in it,
    and a fold ID naming a row the data never had is still an error.

  - area_of_applicability() applied folds built on other rows without a
    word. Fold splits are row IDs, and cv_*() compare the sample of row
    locations make_folds() records against the data, refusing folds
    built on another layer. area_of_applicability() did not: a model
    fitted on the same rows in another order took the folds anyway and
    moved the threshold (0.2432 against 0.2404). It now makes the same
    check on the training data and refuses such folds with the cv_*()
    message. The check is skipped (logged) when the folds were built on
    polygons and the training data are the points a fit reduced them to,
    so a model fitted on polygons with folds built on those polygons
    keeps working.

  - predict_surface() filled a polygon grid with covariates from an
    arbitrary point inside each cell. st_nearest_feature() returns the
    first zero-distance match the spatial index yields, so a
    create_grid_polygons() grid got covariates that did not match the
    location predicted at (predictions off by up to 2.6 on a 0-30
    response) and that changed with the row order of covariates (by up
    to 4.7), and the result was a polygon layer where the manual
    promises points. The grid is now reduced to one representative point
    per cell first.

  - predict_surface(..., draws = TRUE) flattened the draw matrix into
    .pred. The argument was forwarded to predict(), and a backend that
    honours it returned an n_draws x n matrix that became .pred column
    by column (correlation with the right values: 0.03); with se = TRUE
    the duplicated argument was reported as "backend does not expose
    posterior draws". draws is now refused by name, and a predict()
    returning the wrong number of values is an error.

  - predict_surface()'s automatic grid could lose a whole column or row.
    When the extent is an exact multiple of the cell size, floor() of
    the ratio landed one short through rounding (0.3 / 0.1 gives 2
    cells), leaving a cell-wide strip uncovered; at the default n_cells
    this hit 1197 of 10000 random squares. The cell count now has a
    relative tolerance, so an exact multiple gets exactly that many
    cells.

  - predict_surface() kept a reused grid's old .pred_se. Passing an
    earlier surface as grid left that model's .pred_se beside the new
    .pred unless this call replaced it, even after logging "returning
    predictions only". It is now removed unless this call computes it.

  - A logger configuration made before loading the package could still
    abort its functions, and received its log lines. logger seeds a new
    namespace by copying every index of the user's global configuration,
    and 2.0.0 pinned the formatter on index 1 only. A user with two
    global indices set up before library(spatialkit) got
    formatter_sprintf or formatter_glue on the console echo, so a % or a
    { in a message ("fold 2 skipped: object 'cov_{x' not found") aborted
    the function that logged it and the R warning the manual promises
    never arrived; a third global index kept the user's own appender and
    received spatialkit's WARN and INFO lines in the user's log file.
    Both indices now have formatter, layout, appender and threshold
    pinned, every message is marked logger::skip_formatter(), and copied
    indices beyond the second are deleted (on logger 0.2.2, which cannot
    delete one, switched off). The global configuration is still never
    touched.

  - Deleting the session temp directory made every function that logs
    fail until the package was reloaded. The trace file's path was fixed
    in tempdir() at load time, so after an OS cleaner or
    unlink(tempdir()) every log call failed with "cannot open the
    connection", and a documented R warning (ensure_projected()'s CRS
    assumption, say) became that error; tempdir(check = TRUE), R's own
    recovery, did not help. The trace now resolves its path when a line
    is written, recreates the directory if it has gone, and drops a line
    it cannot write; no logging failure aborts the caller any more, so
    the warning always arrives.

  - Logged cautions were missing from knitted documents. The console
    echo wrote to stderr, which knitr does not capture, so an R
    Markdown, Quarto or pkgdown document showed the package's R warnings
    but none of its logged cautions (compare_models()'s significant
    residual autocorrelation, for one). While knitr is running the line
    is now also sent as an R message, so it appears in the output and
    message = FALSE hides it; a line that is raised as a warning too is
    not repeated. stderr gets exactly what it got before, and nothing
    changes outside knitr.

  - plot_tessellation_map(labels = TRUE) drew no labels on any layer
    this package builds. label_col defaulted to "grid_id", a column no
    function produces: Voronoi and Delaunay cells carry cell_id, grids
    poly_id and cell_id, and summarize_by_cell() output poly_id, so the
    map came back unlabelled with only a log line to say why. label_col
    now defaults to NULL, which takes the first of grid_id, cell_id,
    poly_id, polygon_id and id the layer has; grid_id stays first, so a
    layer that has one is labelled as before, and naming a column still
    works.

  - plot_tessellation_map() failed at print on a units, Date, POSIXct or
    difftime fill column. The fill scale was chosen with is.numeric(),
    so Date, POSIXct and difftime columns got a discrete scale
    ("Continuous value supplied to a discrete scale"), and an st_area()
    column (class units) passed the test and then broke the viridis
    scale's arithmetic. The function returned normally and the error
    came only when the plot was drawn. Date and POSIXct now get the
    continuous scale on a date or time axis, and units and difftime
    columns are drawn as numbers with the unit in the legend title
    ("area [m^2]") unless legend_title is given.

  - plot_folds() failed at print when its layers disagreed on having a
    CRS. A CRS-less boundary beside projected points, or the reverse,
    aborted inside coord_sf() with sf's "cannot transform sfc object
    with missing crs". Since plot_folds() began drawing the block
    outlines it also failed on the very layer the folds were built from:
    make_folds() projects CRS-less lon/lat points to a UTM zone, and
    stamps CRS-less points with a boundary's CRS, so the stored blocks
    carry a CRS the points do not. A CRS-less layer is now brought into
    the points' CRS, or failing that the folds' own, or the first
    layer's that has one: reprojected when its coordinates look like
    lon/lat, and stamped with a warning otherwise, as make_folds() does.

  - plot() on a custom spatial_fit without a residuals() method stopped
    with "could not extract residuals". ?new_spatial_fit calls that
    method optional, but residuals.default() returns NULL for a
    spatial_fit, so the residual map, the observed-against-predicted
    plot and the residual variogram all refused a backend that had the
    required fitted() method. The residuals are now the response minus
    fitted(), which is what the built-in backends return; without a
    fitted() method the error names the method to define.

  - citation("spatialkit") gave the year it was called in, not the year
    of the release. DESCRIPTION has no Date field, so inst/CITATION fell
    back to Sys.Date(). CRAN installs got the right year only by
    accident: meta$Date partially matched Date/Publication. The year is
    now looked up the way utils::citation() does it, by exact field
    name: the CRAN publication date, then Date, then the date R CMD
    build packaged the source (which a GitHub install via remotes or pak
    has). Only an install straight from a source directory records none
    of these, and only then does the current year appear.

  - The plots' size arguments did nothing on ggplot2 older than 3.4.0.
    The line layers pass linewidth =, which ggplot2 3.4.0 introduced;
    older versions warn "Ignoring unknown parameters" and draw at the
    default width, so outline_size, boundary_size and the other size
    arguments were silently ignored. Suggests now asks for ggplot2
    (>= 3.4.0).

  - The test suite calls local_mocked_bindings(), which testthat added
    in 3.1.7, but Suggests allowed 3.1.5. On 3.1.5 or 3.1.6 every test
    that mocks a function failed with "could not find function".
    Suggests now asks for testthat (>= 3.1.7).

  - determine_optimal_levels(criterion = "combined") could put first a
    cell count that Moran's I never scored, chosen by the rule the elbow
    had stopped using. On eight separated clusters (800 points,
    max_levels = 40, four seeds) it returned 10 7 6, 10 7 6, 6 5 10 and
    6 10 7. The geometric axis was still the chord on linear axes across
    the elbow's window, which ranked 6 or 7 above the elbow of 8. The
    candidates below the nine-cell floor, which Moran's I cannot score,
    shared an average rank that shrank as more of them went unscored (6
    of 9 for seven of them), although the help page said they ranked
    last. The geometric axis is now the log-log sag the elbow is read
    from, and it is flat when the curve has no elbow, so Moran's z alone
    orders the window. Every unscored candidate takes the last place on
    the Moran axis, and exact ties go to the k nearest the elbow. And an
    elbow below ten cells, a count Moran's I cannot score, is no longer
    ranked against the counts it can: ranking it put ten, the smallest
    count Moran's I scores, first whatever the response did (a response
    of noise and one varying by cluster gave the same answer). There
    "combined" returns the geometric ranking, logs why, and records it
    in the diagnostics (criterion = "geometric", fallback). The same
    layers now give 8 7 9, 8 7 9, 7 6 8 and 8 7 9, the geometric answer,
    and 800 uniform points, which have no elbow and are ordered by
    Moran's z, 10 11 7 where they gave 5 4 10. Supplying both
    response_var and predictor_vars selects this criterion by default.

  - build_tessellation(method = "hex") or "square" laid its lattice over
    a near-global lon/lat boundary in Web Mercator. The points' CRS is
    chosen for distances, and handed on as the grid's CRS it skipped the
    area check create_grid_polygons() makes: on a boundary from 170W
    to 170E and 60S to 70N, the full hexagons differed 5.75-fold in true
    area, where create_grid_polygons() on the same boundary used Equal
    Earth (0.7 percent). The lattice is now laid where
    create_grid_polygons() lays it: in the CRS picked for the points
    unless that CRS distorts areas across the boundary by more than 1
    percent, and otherwise in the equal-area CRS
    ensure_projected(purpose = "area") picks for the boundary, with a
    logged warning, the points indexed in the same CRS. This applies
    with no crs and with a geographic one; a local extent keeps its UTM
    zone.

  - A CRS-less study area given with lon/lat points could be read as a
    one-metre square. build_tessellation() resolved a boundary without a
    CRS against the UTM zone picked for the points, after projecting
    them, so a one-degree tile with integer corners (which the lon/lat
    heuristic declines) was stamped with that zone: one Voronoi cell,
    or 27 hexagons, and all 50 points indexed NA. A boundary in British
    National Grid metres was stamped with the UTM zone too. Such a
    boundary is now read in the points' own CRS when its coordinates fit
    the lon/lat envelope, with a warning, and refused with an error
    naming both layers when they do not; create_voronoi_polygons() and
    clip_target_for() read it the same way.

  - CRS-less lon/lat points with a CRS-less boundary in metres failed
    with "boundary must be polygonal". build_tessellation() stamped
    EPSG:4326 on the boundary without looking at its coordinates, so a
    UTM polygon was transformed to nothing and the error named its
    geometry type. It now stops with an error saying the two layers
    cannot be placed in one space.

  - create_voronoi_polygons() tessellated CRS-less lon/lat points in
    degrees, silently. It projected only when a CRS said lon/lat, so 60
    CRS-less points at 55N got cells in which 17.7 percent of sampled
    locations were not nearest to their cell's point, while
    build_tessellation() on the same points warned, took them as
    EPSG:4326 and projected them. It now applies the same lon/lat
    heuristic, with its warning, and returns what build_tessellation()
    returns (0.1 percent, at the cell edges). ?build_tessellation no
    longer says a CRS-less pair "stays in the same unnamed planar space"
    whatever its coordinates.

  - Delaunay triangles returned in a geographic crs did not contain
    their own points. build_tessellation(method = "triangles", crs
    = 4326) on 150 points left 9 of them touching no returned triangle,
    so a spatial join on the result did not reproduce index. The corners
    of the returned triangles are now put back on the input points after
    the round trip through the working projection, and every point lies
    in its indexed triangle.

  - The tessellation builders failed on an sf layer as crs.
    build_tessellation(), create_voronoi_polygons() and
    create_grid_polygons() stopped with an error from sf ("the condition
    has length > 1", or "cannot create a crs from an object of class
    sf"). They now take the layer's CRS, as ensure_projected(target_crs
    =) and harmonize_crs() do.

  - Random and k-means seeding on a lon/lat boundary warned "install
    package lwgeom" on every call. sf raises "coordinate ranges not
    computed along great circles" for each lon/lat draw when lwgeom,
    which this package does not depend on, is absent: one R warning per
    voronoi_seeds_random() or get_voronoi_seeds(method = "random") call
    and two per k-means call. That warning is muffled, and the draw is
    unchanged. The boundary's union is taken on the sphere as well, so
    with sf_use_s2(FALSE) sf no longer prints its planar st_union()
    message and the seeds are the ones an s2-on session gets.

  - A transect with sub-millimetre scatter got a sliver study area and a
    166,536-cell grid for approx_n_cells = 25. clip_target_for() called
    a bounding box degenerate only when its two ends were equal to
    rounding, so 30 points along 1000 m with a y scatter of 1e-6 got
    a 1000 x 9e-7 rectangle, over which the square grid had 166,536
    cells (16 seconds) and the hexagonal one 154,980 (28 seconds); a
    little thinner, and create_grid_polygons() stopped at max_cells
    telling the user to check the units of a cellsize they had not
    passed. A box whose short side is below a millionth of the long side
    is now degenerate too (a buffer around the points: 34 squares or 45
    hexagons for 25), and the max_cells error names the argument the
    size came from, target_cells (approx_n_cells), n or cellsize.

  - A whole build_tessellation() result passed as boundary failed with
    sf's no applicable method for 'st_geometry' applied to an object of
    class "list". build_tessellation(), create_voronoi_polygons() and
    clip_target_for() now say that the object looks like a
    build_tessellation() result and to pass its $boundary, as
    make_folds(blocks = ) already did for $cells;
    create_grid_polygons(), create_grid_polygons_cached() and
    ensure_stable_poly_id() add the same hint to their type errors.

  - build_tessellation(method = "triangles") recorded an approx_n_cells
    it had ignored. params$approx_n_cells held the ignored count where
    the help page says the count used is kept; it is now NULL, as for
    Voronoi, and the warning says so for both methods.

  - summarize_by_cell(deff = "variogram") said it was "Falling back to
    deff = 1" for a rejected sac, and then applied a design effect. A
    sac whose fit estimate_sac_range() had rejected was set aside with
    that warning, after which a variogram estimated from response_var
    was fitted and applied: in one check every row came back corrected,
    with a median design effect of 5.2, and in the development version
    the warning carried the fallback class, so
    tryCatch(spatialkit_deff_fallback = ) threw the corrected result
    away. When the estimate was rejected too, the one fallback raised
    two R warnings. A sac that an estimate replaces now gets a plain
    warning saying what replaced it, and the fallback warning is raised
    once per call, only when the standard errors really are the
    uncorrected ones, naming every reason, the rejected sac included.

  - summarize_by_cell(deff = "variogram") ignored a sac that carried no
    variogram model without saying so. A plain number, or
    units::set_units(1.5, "km"), was passed over and the design effect
    came from a variogram estimated from response_var, with no
    condition; the caller could not tell that the value given had not
    been used. Such a sac is now set aside with a warning, as a rejected
    sac is: a plain warning when a variogram is estimated instead, and
    the classed spatialkit_deff_fallback warning, naming it, when deff
    falls back to 1.

  - summarize_by_cell(deff = "kish") recorded no correction when only
    the predictor standard errors were corrected. The "deff_applied"
    attribute followed the primary variable's ICC alone, so with an
    unclustered response (ICC 0) and a clustered predictor (ICC 0.82) no
    attribute was attached, and in the development version every row
    said deff_applied = FALSE, while the predictor standard errors had
    been inflated elevenfold. The attribute is now attached whenever
    either ICC is positive; its deff stays the primary variable's (all 1
    in that case).

  - summarize_by_cell() corrected the response's standard errors with a
    residual variogram without a word. A sac from
    estimate_sac_range(..., predictor_vars = ) describes what the
    predictors leave unexplained, a weaker correlation than the
    response's own: the response standard errors came out at 0.60 of
    those from the response variogram in one check, while
    kriging_adequacy() warned about the same object. Such a sac is still
    used as given, as documented, but a warning now says the response
    standard errors are understated.

  - assign_features_to_polygons() dropped the features' own id column
    when the cells were keyed by id. The polygons' ID went through the
    spatial join under its own name, so a site id collided with the
    cells' id and was dropped, with a warning about a collision the
    result never had: it only gains polygon_id_col. Only a column named
    polygon_id_col is replaced now.

  - assign_features_to_polygons(largest = TRUE) let the polygon row
    order decide an exact tie in overlap. sf keeps the first of equal
    largest overlaps, so a 40 m square split evenly across the edge of
    cells 1 and 2 went to cell 1, or to cell 2 with the polygon rows
    reversed. An exact tie (areas equal to 9 significant digits) is now
    decided by tie_break among the equally largest polygons, and counted
    in attr(, "ties").

  - Every largest-overlap assignment of polygon features raised sf's
    "attribute variables are assumed to be spatially constant" warning.
    sf::st_join(largest = TRUE) adds grouping columns of its own before
    intersecting, so the warning came with every ordinary call (the
    reporting vignette hid it with warning = FALSE), said nothing about
    the data, and could not be avoided with st_agr(). That warning alone
    is now muffled.

  - summarize_by_cell(deff = "variogram") accepted a deff_max_n of less
    than 2. A value of 1 or 0 subsampled every cell to one point or
    none, so the mean correlation came out NaN, the design effect 1 and
    the standard errors uncorrected (4.1 times smaller than with the
    default in one check) with nothing to say so; NA stopped the call
    with "missing value where TRUE/FALSE needed". It must now be a
    single number of at least 2, and anything else is an error that
    names it.

  - The rows with no cell ID got a cell_weight of 0. A layer assigned
    with keep_unassigned = TRUE is summarised with its unassigned rows
    as a group whose ID is NA, and summarize_by_cell() lost that group's
    count of non-missing values: cell_weight was 0 beside n = 5 and a
    finite standard error. The group is now counted like any other.

  - attr(, "deff_applied")$deff turned into a vector of NAs for a fixed
    deff and one populated cell. With cells_sf, the realignment to the
    joined rows took a fixed deff = 2 for a per-cell vector whenever
    exactly one cell was summarised, and recorded c(2, NA, NA, ...). A
    fixed design effect is now recorded as the number.

  - make_folds(method = "block_kfold")'s refusal of a grid
    above 1,000,000 blocks told the caller to check a block_size they
    never passed, and a block_size hundreds of orders of magnitude too
    small slipped past it. With block_nx = 2000, block_ny = 1000 (or
    block_multiplier = 1e6) the error read "Check that block_size
    (unset) is expressed in the data's CRS units". With block_size
    = 1e-200 the cell count overflowed to Inf, which the guard let
    through, and st_make_grid() failed with "result would be too long a
    vector". The message now names what produced the grid: block_size,
    the range auto_range estimated, block_nx/block_ny, or
    block_multiplier x k. A count too large to represent is refused like
    any other.

  - make_folds() accepted any block_multiplier, and a units object for
    phi or min_train failed with an error that named no argument.
    block_multiplier = NA died inside st_make_grid() with "'length.out'
    must be a non-negative number". c(1, 3) silently used 3, and
    units::set_units(100, m) was read as 100, giving a 32 x 16 grid
    aimed at 500 blocks. phi = units::set_units(100, m) failed inside
    the units package with 'both operands of the expression should be
    "units" objects'. block_multiplier must now be a single positive
    number, and phi and min_train refuse a units object by name, as
    block_size, buffer and block_nx/block_ny already did.

  - residual_morans_i() and compare_models() had no residual Moran's I
    for a custom fit without a residuals() method. ?new_spatial_fit says
    residuals() is optional, and such a fit falls through to
    residuals.default() and gets NULL, so residual_morans_i() returned
    NULL with a warning and compare_models() reported all-NA
    resid_morans_* columns. On 80 points with an east-west trend the
    predictor could not explain, that hid a residual Moran's I of 0.754
    (z = 15.5, p = 4e-54). Such a fit is now scored on the observed
    response minus fitted(), which is what the built-in backends'
    residuals are and what plot() already used for it. NULL is returned
    only when that cannot be formed either, and the warning quotes the
    reason.

  - Re-using a cv_*() result's $folds renumbered every fold after a
    dropped one. The result's $folds holds the splits that survived,
    each carrying the fold_id it was reported under, so five folds with
    fold 3 dropped are labelled 1, 2, 4 and 5. Handed to a second cv_*()
    call, to score another learner on the same splits, they were
    numbered by position as 1, 2, 3 and 4, so the second run's fold 3
    was the first run's fold 4 (RMSE 2.90 in both), and a join on fold
    with the first run or with fold_separation() paired different folds.
    Splits that carry a distinct whole-number fold_id now keep it.
    Splits without one are still numbered by position, as are splits
    whose ids repeat.

  - A fold_info_fn whose ..per_row reused a predictions column name
    corrupted predictions. A ..per_row column called yhat (or y, fold,
    ..row_id or y_train_mean) was bound on beside the original, and
    predictions came back with two yhat columns. In the development
    version, dplyr::bind_rows() then renamed both (yhat...4, yhat...6)
    with only a message, so overall found no yhat and was all NA with
    n_pred = 0, beside per-fold RMSEs of 2.4 to 3.4 and no warning. Such
    a name is now an error naming the column, as a clashing scalar extra
    already was. Duplicated or empty ..per_row names are errors too, in
    a parallel run as in a sequential one.

  - A fold_info_fn that returned a named vector lost its values without
    a word. c(slope = 1.9) in place of list(slope = 1.9) (the shape
    metrics accepts) added no column, and fold_status read "ok" with an
    empty message on every fold. A named vector is now taken as the list
    it stands for. Any other return value that is not a list or NULL (a
    function, an environment) is now an error.

  - Folds built on a pointized copy of a polygon layer were refused as
    "built from different data". make_folds(coerce_to_points(parcels))
    applied to parcels itself, with the same rows and IDs, was refused
    on 80 L-shaped parcels because "64 of 64 checked row IDs sit at a
    different location here". The error blamed "folds from another
    dataset". The provenance check compares each polygon's centroid with
    the point coerce_to_points() gave it, and under "auto" that is
    st_point_on_surface(), which differs for every non-convex feature.
    The folds' params$row_probe now records whether they were built on
    POINT geometry. When that differs from the data, the error says so
    and tells you to build the folds with make_folds() on the layer
    passed, which reduces polygons to points itself. Folds from data
    that really differs get the old message.

  - evaluate_insample() returned NULL when no element of fits was a
    spatial_fit. Its help promises a data frame with one row per model,
    and the only sign that every element had been skipped was a log line
    per element, which spatialkit_quiet() hides and tryCatch(warning = )
    never sees. It is now an error that names fits.

  - cv_rf(parallel = ) printed its core-count message twice when asked
    for more workers than the machine has. cv_rf(parallel = 16) on
    a 4-core machine printed "cv parallel: 16 workers requested on a
    machine with 4 cores; using 4." twice. It resolved the count once
    for its thread policy and once more when the folds ran. It now
    prints the message once.

  - fit_gwr_model() called local regressions unstable because of where a
    predictor's units start. A temperature field in kelvin beside a
    second predictor drew "global design (intercept + predictors) has
    scaled condition index 230 ... (collinearity risk)" and "100% of 30
    sampled locations have a collinear local design ... Local
    regressions there are unstable", while the same field in degrees C
    drew no global warning and 6 of 30; the local slopes of the two fits
    agree to 1e-10. Both indices were computed on the uncentred design,
    where a predictor far from 0 against its spread (kelvin, a year,
    elevation in feet) is collinear with the intercept. That makes the
    local intercept an extrapolation to 0, but a slope's precision does
    not depend on where the predictor's origin is. The global index is
    now computed on the centred predictors, so it is 1 for a single
    predictor and a change of origin does not move it. The local survey
    keeps Belsley's uncentred index with the intercept as
    info$local_collinearity$cn and adds cn_slopes: the predictors
    centred at their weighted mean in the window, each divided by its
    standard deviation over the study area. A window's slopes count as
    collinear when cn_slopes is above 30 or singular, or when cn is
    above 1e6, where GWmodel's uncentred solve starts losing precision
    in the slopes. n_local_collinear and the warning count those
    windows. The kelvin and degrees C fits now get the same verdict (0
    of 200 locations), and a regional covariate nearly constant within
    clusters is still flagged everywhere (200 of 200). A window where
    only cn is above 30 is logged, not warned about.

  - gwr_model_selection() checked an adaptive bandwidth less strictly
    than fit_gwr_model(). bandwidth = 3e9 with adaptive = TRUE, a
    distance passed as a neighbour count, failed with "NAs introduced by
    coercion to integer range" and then a bare "missing value where
    TRUE/FALSE needed". bandwidth = 0.5 was rounded to 0 neighbours and
    quietly raised to the floor, where fit_gwr_model() refuses it.
    gwr_model_selection() now runs fit_gwr_model()'s check before it
    prepares anything: an adaptive count must be at least 1 and at most
    R's largest integer. cv_gwr() runs the same check once, up front,
    instead of failing it in every fold and returning "all folds
    failed".

  - gwr_model_selection() with a fixed bandwidth in the wrong units
    failed with a bare "inv(): matrix is singular". On lon/lat data
    projected to EPSG:32617, bandwidth = 0.2, adaptive = FALSE is 0.2
    metres against an extent of 44720 metres, so every window is empty.
    fit_gwr_model() warned about exactly this and explained the singular
    window, but the sweep said nothing. It now raises the same "less
    than a ten-thousandth of the data's extent" warning, and its error
    explains a singular window as fit_gwr_model()'s does.

  - print() on a GWR fit showed a fixed bandwidth in scientific notation
    and without its unit. A fixed bandwidth of 122372 m printed as
    "Bandwidth: 1.224e+05 (fixed, bisquare kernel)", so the value was
    rounded to four digits and its unit was missing, although the help
    page tells the reader to check it. It now prints "122,372 metre
    (fixed, bisquare kernel)", and an adaptive one as "42 neighbours
    (adaptive, ...)".

  - fit_gwr_model() did not refuse a character or factor response, as
    the README says every model function does. A response read from a
    CSV as text went into GWmodel::bw.gwr(), which failed twice with
    "Not compatible with requested type". That drew the
    arbitrary-fallback bandwidth warning, and the fit then stopped with
    "'x' must contain finite values only", naming neither the column nor
    the cause. It now stops first with "response 'zc' is not numeric (it
    is character)", as fit_rf_model() and fit_bayesian_spatial_model()
    do.

  - cv_bayes() under an ordinal or categorical family sampled every fold
    and then scored none of them. Such a family predicts a probability
    per response category, and cross-validation scores one number per
    row, so each fold compiled and sampled a full model and was then
    discarded: under brms::cumulative(), k = 2 on 50 rows took 4.35
    minutes to return "all 2 folds failed to produce predictions".
    cv_bayes() now refuses brms::categorical() and the ordinal families
    (cumulative, sratio, cratio, acat) before fitting anything, with a
    message naming the family, and the "Which metrics survive a
    non-Gaussian response" section (in ?model_metrics, ?cv_bayes and
    ?compare_models_cv) no longer calls CRPS and interval coverage
    meaningful for "any family the backend accepts".

  - predict_surface()'s automatic grid left up to a cell of the training
    extent uncovered on the east and north. It took floor(extent /
    cell_size) cells from the lower-left corner, so whenever the extent
    was not a whole number of cells the rest of it got no prediction.
    cell_size = 100 on a 980 x 956 extent covered 900 x 900: 13.5
    percent of the box was uncovered, and 14 of 120 training points lay
    in no cell. cell_size = 334 gave 2 x 2 cells and left 65 of the 120
    points out. The grid now has enough cells to cover the box and is
    centred on it. It overhangs the box by less than one cell, split
    evenly between the two sides, so every cell centre still lies inside
    the box. An extent that is an exact multiple of the cell size gets
    the same grid as before. At the default n_cells, a grid usually
    gains one column or row (102 x 99 cells instead of 101 x 98 on the
    extent above), and its centres move by less than half a cell.

  - plot_tessellation_map(labels = TRUE) warned at print on every
    lon/lat layer, and failed at print on a units label column. The
    label points were computed with st_point_on_surface(), and
    geom_sf_text() ran it again on those points when the plot was drawn.
    On longitude/latitude cells that raised "st_point_on_surface may not
    give correct results for longitude/latitude data", although nothing
    was wrong. A label_col holding st_area() values (class units) failed
    with "units package is not attached", as a units fill column did.
    The label points are now used as computed, and a units or difftime
    label is drawn formatted, with its unit ("15073.393 [m^2]"). A
    fill_col or label_col naming more than one column, or NA, is now
    refused by name. It used to fail with R's "'length = 2' in coercion
    to 'logical(1)'".

Documentation

  - Every figure in the vignettes carries alt text, which is what a
    screen reader announces and the only thing a reader gets when an
    image fails to load. Each one states what the picture shows and what
    it is there to demonstrate, rather than naming the axes.

  - Every help page now ends with a "See also" that leads somewhere;
    more than a third had none, every S3 method among them. Four
    families are new: spatial data preparation (ensure_projected(),
    harmonize_crs(), coerce_to_points(), clip_target_for(),
    prep_model_data()), package options and caches (spatialkit_quiet(),
    the two cache clearers, create_grid_polygons_cached()), methods on a
    fitted model (the predict(), fitted(), residuals(), coef(), print()
    and summary() methods of the three backends, which document one
    contract between them) and print methods for the other result
    objects. The seed generators, the cached grid constructor and
    ensure_stable_poly_id() join tessellation, model_metrics() joins
    model evaluation, sac_nugget() joins cross-validation,
    gp_lengthscale_bounds() joins model fitting, and
    area_of_applicability() is now in prediction as well as
    cross-validation. The website's reference index moves the two cache
    clearers into the same group, so it and the help pages agree.

  - Ten numbered scripts are installed with the package, in
    system.file("scripts", package = "spatialkit"). 00-run-all.R runs
    them in order; each of 01- to 10- is self-contained and covers one
    topic: tessellations, resolution, fold schemes, block sizing,
    fitting and diagnosing, model comparison, prediction surfaces and
    the area of applicability, feature selection, GWR and the Bayesian
    GP. They run on a simulated field with known structure, print what
    they are doing, and state what to look for in a figure before
    drawing it. No script asserts a number it did not compute, so a run
    that finishes is one whose claims held on the machine that ran it.
    SPATIALKIT_TOUR_OUTPUT writes the figures to a folder instead of the
    device; SPATIALKIT_TOUR_PAUSE = "no" skips the per-figure pause. A
    script skips the part that needs an absent package with a message
    naming it, and scripts 02, 09 and 10 skip themselves when gstat,
    GWmodel or brms is absent.

  - Four vignettes join spatialkit_nc_demo: getting-started (installing,
    what the coordinates are in, the pipeline from points to a scored
    model and a map, and a glossary of the terms that recur), resolution
    (the ladder, the four criteria and why they disagree),
    spatial-cross-validation (the five fold schemes, block sizing, and
    reading a CV result down to the last row) and diagnostics (residual
    autocorrelation, aggregation standard errors, kriging adequacy, area
    of applicability, and two ways to leak). Each is executed at build
    time and gates itself on the optional packages it needs.

  - The README is generated from README.Rmd, so its figures and every
    number in it are computed when it is built rather than pasted in. It
    is about half its previous length: the material that had accumulated
    in it moved to the vignettes above, and what remains is what the
    package is for, a quick start that shows a real cross-validation
    gap, the troubleshooting list, and pointers to the rest.

  - New vignette reporting: what leaves the session at the end of a run.
    The first half is the regions as a file someone else can use.
    Grouping a layer you already have with
    assign_features_to_polygons(largest = TRUE) (100 North Carolina
    counties into a hex grid of 18 regions, 12 of them populated, one
    row per county), answering whether a location falls in one with
    keep_unassigned = TRUE, IDs that survive a reprojection, and why the
    aggregates go into a GeoPackage instead of a shapefile: of the 10
    columns summarize_by_cell() produces, 2 come back from a shapefile
    with their names intact. The second half is a worked report of six
    numbers, each read out of an object the run already produced, with
    what each one is there to stop a reader believing (Roberts et
    al. 2017; Meyer and Pebesma 2021; Heaton et al. 2019). Its example
    reports a run that fails its own checks, which is the case the
    section exists for.

  - The documentation is published as a website at
    https://elkronos.github.io/gis_modeling_toolkit/: the README, every
    help page grouped by pipeline step, the six vignettes as articles,
    and this changelog. It is rebuilt from main on every push, so it
    describes the development version; the "development version" heading
    at the top of this file lists what the CRAN release does not have
    yet.

  - select_features_forward() now says what $score is: the
    cross-validated metric of the winning set at the final step, which
    is the selection criterion and is optimistically biased by the
    selection itself (Cawley and Talbot 2010) --- not a performance
    estimate of the selected model. The honest estimate comes from
    running the selection inside cv_spatial()'s fit_fn, which the page
    now spells out.

  - fit_bayesian_spatial_model() documents which brms families family =
    takes --- zero-inflated and hurdle counts, negative binomial,
    Bernoulli, beta, ordinal, categorical and mixture families all reach
    brms::brm() with the spatial GP term intact, and a factor response
    is taken only under categorical() and the ordinal families --- and
    gains a worked zero-inflated Poisson example. The same page now
    records a trap: a family object brms cannot name skips the
    response-type check entirely rather than falling back to the
    gaussian rule, so a malformed family buys less validation, not more.

  - fit_bayesian_spatial_model() gains a "Spatial confounding" section:
    a coefficient estimated beside a spatial random effect is a
    different estimand from the non-spatial one (Zimmerman and Ver
    Hoef 2022), can shrink toward zero when the response is smoother
    than the covariate (Bolin and Wallin 2025), and the honest
    diagnostic is to report both side by side. The section names both
    sides of the restricted-spatial-regression dispute (Hughes and
    Haran 2013; Hanks et al. 2015; Khan and Calder 2022) and the
    remedies expressible in a brms formula (Marques, Kneib and
    Klein 2022; Guan et al. 2023).

  - model_metrics() gains a "Which metrics survive a non-Gaussian
    response" section, inherited by cv_bayes() and compare_models_cv():
    RMSE and MAE are defined for any numeric response; MAPE, SMAPE and
    R-squared are Gaussian-shaped; CRPS and interval coverage from
    cv_bayes() are the proper scores for a count or bounded response. It
    also records that the all-folds-failed fold_metrics frame omits the
    coverage_* columns.

  - estimate_sac_range() documents what a count or other
    mean-variance-linked response does to the variogram, why detrending
    with predictor_vars helps but does not fix it, and what to do
    instead.

  - ?spatialkit now reads its nine-step pipeline as an argument ---
    steps 1 to 4 are the claim, steps 5, 7 and 9 the evidence that makes
    it checkable --- and gains a "Defaults and their sources" section
    listing which defaults rest on a citation and which were chosen, so
    the two are not mistaken for each other.

  - cv_spatial() documents its name collision with
    blockCV::cv_spatial(), which builds folds where this one runs them,
    and that blockCV's $folds_ids is accepted directly as folds
    everywhere.

  - ?determine_optimal_levels no longer calls the Cliff and Ord moments
    behind Moran's z exact: they assume cell means of equal variance,
    and with single-point cells beside cells of 70 or more points z ran
    slightly high (mean 0.2--0.34, 7--8% rejection at 5%);
    resolution_profile() says cell_diam_median is about 0.8 of the side
    of an equal-area square cell, not a width (and
    vignette("resolution") no longer calls it one); the Post-selection
    inference section and the vignette say that the estimation rows of a
    split leave cells in the selection half empty and cells across the
    border estimated from part of their points (coverage 0.88 against
    0.97 in simulations), and how to find the cells to trust;
    sac_nugget() says the nugget is extrapolated from the first lag bin
    (about max_dist / 30 at the default cutoff), not observed below the
    closest pair, and can be exactly 0 at the fit's bound.

  - The tessellation help pages are corrected: clip_target_for() returns
    the points' bounding box, not their convex hull, and is not the
    target build_tessellation(method = "voronoi") derives;
    create_voronoi_polygons() and build_tessellation() say that a
    multi-vertex MULTIPOINT feature gets one cell per vertex and that
    keep_duplicates has no effect; voronoi_seeds_random() tops a short
    draw up to exactly k; voronoi_seeds_kmeans() and get_voronoi_seeds()
    say that their stats::kmeans() partition is not the k-means++ one
    resolution_profile() scored, and no longer promise equal counts per
    cell.

  - Documentation corrected: summarize_by_cell()'s sac no longer
    identifies a rejected fit by a status attribute no sac_range carries
    (it is an NA value with a rejected_reason), and says a rejected sac
    is set aside and the variogram estimated as if none had been given;
    ensure_stable_poly_id() and the reporting vignette say IDs are the
    same across projections except for fine cells within the rounding
    step of one longitude (36 of 2,500 100 m cells via EPSG:3035), not
    always; summarize_by_cell()'s "Confidence intervals" says ..neff_*
    is NA for a single observation; the resolution vignette says the
    blocked split reduces the leak between halves rather than removing
    it.

  - ?estimate_sac_range now says that the exponential model is kept
    whenever it converges and so overestimates the range on smoother
    fields (about 1.8--2.1 times a Gaussian practical range, 1.3--1.4
    times a spherical one), and why the family is not chosen by fit
    (that biases exponential fields low, to 0.82); that 15 lag bins make
    a range spanning one or two of them come out long (60 m
    returned 89--102 m); that nothing tests for spatial structure (white
    noise gave a finite range in 8 of 30 draws); that the "decreases
    with distance" refusal also fires on small samples of ordinary
    fields (7--9 of 60 at n = 30, 0--1 at n = 100), which its warning
    now says below 100 points; that the anisotropy note goes to the log
    file at INFO, not the console, and that the longest directional
    range is read from directional_fitted after checking
    directional_status, not max(attr(, "directional")), which is NA
    exactly when the major axis ran past the fitted lags (the note
    itself, the nc_demo vignette and the help now say so); that an
    accepted range can still exceed half the width of the layer, leaving
    make_folds(auto_range = TRUE) room for one block (range_frac); and a
    flat variogram's refusal message no longer claims the data "never
    reached a sill" without mentioning that a structureless variogram
    ends there too. The README counts six refusals, not five.

  - ?make_folds: the NNDM details now say how the procedure differs from
    CAST::nndm() (a strict removal rule, one point more conservative per
    distance value, and ties broken by coordinates rather than row
    index) instead of calling it the same; the n > 5000 refusal gives
    the worst-case cost, O(n^3) time where min_train binds (about nine
    minutes at n = 3000), not O(n^2); drop_empty_blocks, boundary and
    block_multiplier say what they do on the cases above.

  - ?cv_rf no longer says a seed passed through ... overrides the
    per-fold forest seed (seed is the function's own argument and never
    reaches ...); ?cv_bayes says parallel = n compiles n Stan models at
    once, at several GB each, and what a fold killed for memory reports;
    ?residual_morans_i and ?compare_models say their methodological
    cautions are logged, not raised as warnings.

  - ?fit_gwr_model's "Collinearity diagnostics" section described the
    30-location unweighted spot check the survey replaced, a global
    index on the predictors alone, and a caveat that only a subset of
    locations is examined; it now describes the code (every location,
    kernel-weighted, a global index on the centred predictors, and at
    each location one index with the intercept and one for the slopes
    alone). ?fit_gwr_model and ?gwr_model_selection now state the
    adaptive bandwidth's floor and cap.

  - ?fit_bayesian_spatial_model: check_convergence says the checks write
    WARN log lines and set $info$convergence_ok rather than "issue
    warnings", and that under cmdstanr nothing is raised as an R
    warning; the basis adequacy check is described as logged, and as
    changing neither convergence_ok nor print() (the argument's text
    listed it among the checks that set convergence_ok to FALSE); the
    family argument and the non-Gaussian section no longer promise that
    any response type brms can fit works here; the spatial confounding
    section warns that coefficients under standardize_predictors = TRUE
    are per standard deviation before comparing with lm().
    ?coef.bayesian_fit gains a section on standardised predictors, and
    print() on such a fit names them. ?fitted.bayesian_fit no longer
    says a failed posterior draw returns NA (it has been an error since
    before this release). ?gp_lengthscale_bounds says its bounds are the
    prior's calibration range and do not shrink with n.

  - ?area_of_applicability states the outlier rule (type-7 quartiles)
    and how the threshold differs from CAST's (the fence itself, capped
    at the largest training DI, which is larger whenever n_outliers
    > 0), with the threshold = value that reproduces it; the internal
    note that CAST uses boxplot.stats() was out of date, and the package
    page no longer implies the threshold matches CAST. n_new is
    documented as every row of newdata (n_inside + n_outside + n_na),
    not the rows that passed the finite-value filter. ?predict_surface
    and the README say that se = TRUE on a bayesian_fit gives the SD of
    the mean surface and that type = "predict" gives the predictive SD.

  - The diagnostics vignette's area-of-applicability figure alt text no
    longer calls the training curve cross-validated (no folds are passed
    there); ?plot.aoa and ?plot.feature_selection describe the training
    DI and the rejected last step as drawn.

  - The examples that stamp EPSG:32632 (UTM zone 32N) on simulated
    points now put the points inside that zone. The README quick start,
    the fold_separation() example, the resolution, diagnostics and
    spatial cross-validation vignettes, and the fixture the scripts in
    inst/scripts/ share placed them at x and y between 0 and 1000, which
    is on the equator about 4.5 degrees east, outside the zone. They are
    now offset by 500000 m east and 5000000 m north, as the other
    examples already were. Every printed result is unchanged, since the
    package works in planar units; only the coordinates themselves, and
    the graticule on the maps, differ.

  - The examples of resolution_profile(), select_resolution(), summary()
    and plot() on a profile use set.seed(4), on which their comments
    hold (C_p interior at 28 cells, reliability on the range floor);
    with set.seed(2) C_p had moved to the support ceiling. The plot()
    example no longer promises four panels on points that have no elbow.

  - Tour script 02 runs to the end again: step 02.6 drew the elbow's
    pick, which its evenly spread points no longer have, and now draws
    the C_p pick. It prints the bound each criterion's optimum sits on,
    points to min_cell_n and range_floor rather than n_levels for moving
    one, and computes what the select_on = "split" comparison shows
    (both picks on the range floor) instead of calling the difference
    tuning. Script 08 no longer calls its held-out half untouched or its
    score gap the selection effect, and says its block size of 300 is
    under the range of about 357.

  - The README's resolution figure labels its middle count as the one
    resolution_profile() rates most reliable, with the width of that
    criterion's flat region; it was determine_optimal_levels()'s count,
    which on those evenly spread sites the function now says the ladder
    chose. The troubleshooting entry quotes the fallback messages
    determine_optimal_levels() now gives, with their causes.

  - The resolution vignette's criteria table describes the elbow as the
    log-log sag it is, NA on points with no cluster structure, and its
    introduction says determine_optimal_levels() still returns a count
    there, with a warning.

  - ?coerce_to_points (tmp_project) states the rule a CRS-less layer is
    read by: the lon/lat heuristic of ensure_projected(), not merely
    lying inside the lon/lat envelope.

  - ?summarize_by_cell now gives the derivation and the measured
    coverage of the small-sample rescaling that goes with every
    data-derived design effect, which the package page said were there;
    the package page lists it among the defaults that were chosen, not
    cited.

  - ?estimate_sac_range now says when the directional maximum is logged
    (only when it stands in for a singular or non-converged all-pairs
    fit, at WARN only above a ratio of 1.5, naming the directional
    ranges), and that a range shorter than the first lag bin can come
    back several times too long (a true 24 m returned 93--479 m at n
    = 1500, 19--32 m at cutoff = 0.1), so a variogram at its sill in the
    first one or two bins calls for a smaller cutoff whatever range was
    fitted.

  - ?make_folds no longer says auto_range "fits directional variograms
    to account for anisotropy". It sizes blocks from the omnidirectional
    range, and for a field known to be anisotropic the page now points
    to directional_fitted. An accepted range too wide for two blocks
    makes make_folds() stop with an error, which the page now says; it
    had claimed the grid "does not collapse to a single block". The log
    line announcing a lowered k just before that error is gone. The page
    no longer says NNDM never pushes a point's nearest-neighbour
    distance past phi: the last exclusion can take it past by one
    neighbour step, as in CAST::nndm(). The description lists all five
    methods. The page now says that k = 1 is raised to 2 by the three
    k-fold methods, and that only block_kfold returns
    params$blocks_supplied and params$boundary_supplied.

  - The README's entry for "response 'y' is not numeric" notes that
    fit_bayesian_spatial_model() takes a factor or character response
    under brms::categorical() and the ordinal families.

  - The README, ?summarize_by_cell and the getting-started, diagnostics,
    reporting and North Carolina vignettes now say which estimand the
    design-effect-corrected standard error is for. It is the standard
    error of a cell mean as an estimate of the population mean. For a
    cell's own mean, which is what a map reports, the default deff = 1
    standard error is the right one when the points are spread through
    the cell. Several of these pages had presented the correction as the
    right standard error for the cells themselves, and there it is too
    wide by sqrt(deff / (1 - rho)) (4.6 at 20 points a cell and rho
    = 0.5). The getting-started pipeline, which maps the cells, now
    aggregates at deff = 1.

  - Smaller corrections: the README's troubleshooting list adds
    "compare_models_cv(): no recognised model requested.", which is the
    error when no requested name is recognised, and says that "no viable
    models." means every recognised backend is uninstalled. The
    getting-started install table no longer says patchwork is needed for
    the plot_*() functions. The North Carolina vignette explains why
    some local designs of its GWR fit have a high condition index with
    the intercept (an uncentred elevation, nearly collinear with the
    intercept inside each window), and why the fit raises no
    collinearity warning: the slope index stays below 30, so the slopes
    are well determined and only the local intercepts are not. Its
    fold-map alt text now describes each blocked fold as whole blocks in
    separate parts of the state, not as one contiguous area.

  - vignette("diagnostics"): the "Two ways to leak" example of selection
    inside the folds leaked itself (blocks about 250 m across against an
    autocorrelation range of about 330 m) and its learner could not fit
    the intercept-only model, so tol did not apply to the first
    variable. It now passes block_size = 400, fits z ~ 1 for an empty
    predictor set, and shows sel$history.

  - ?fit_rf_model recommended area_of_applicability(weights =
    pmax(fit$info$importance, 0)) without condition. It now says this
    works only when that importance is finite, which it is not when no
    row is out of bag.

                 Changes in version 2.0.0 (2026-09-11)                  

Everything below is relative to 1.0.0 (published on CRAN 2026-08-07).
The major bump is warranted: three exported functions are removed, and
several defaults change the result of a fit or a comparison, so the same
script can get a different answer. Both are under "Breaking changes" —
read that section before upgrading a running analysis.

These notes describe what changed for a user of 1.0.0. A good deal of
this release was written after 1.0.0 and then revised before shipping;
defects that existed only between those points are not listed, since no
released version behaved that way. The commit history has that record in
full.

Throughout, raises a warning means a genuine R warning() — one
tryCatch(warning = ) catches, suppressWarnings() suppresses and
options(warn = 2) escalates. Logs a warning means a logger message in
the "spatialkit" namespace, which none of those touch.

Breaking changes

Corrections that change results

Each item below was measured and independently reproduced before it was
touched; the figures quoted are from those reproductions, so you can
judge whether an item affects an analysis you have already run.

  - residual_morans_i() no longer puts weight on a point's own residual.
    FNN::get.knn() reports a point's OWN index among its neighbours
    whenever exact duplicate coordinates are present, which put 1/k on
    the diagonal of a matrix Moran's I is only defined for with a zero
    diagonal. On 40 sites x 4 repeats with a response carrying no
    spatial structure, 120 of 160 rows gained a self-weight, mean I came
    out at +0.086 against E[I] = -0.0063, and 77% of samples were
    "significant" at p < 0.05 against a nominal 5%. Repeat observations
    at one site are exactly what make_folds(method =
    "leave_location_out") is for, so this was a mainstream input. The
    dense fallback never had the fault, so the statistic also depended
    silently on whether FNN happened to be installed; the two paths now
    share one neighbour lookup and agree exactly.
    
    Requesting k + 1 neighbours and dropping self is not sufficient on
    its own — the slot self occupied displaced a genuine co-located
    neighbour and left a farther point standing in for it (75 of 400
    retained pairs sat at distance 121 where a neighbour at distance 0
    existed). Duplicate coordinates are now grouped and answered
    exactly.

  - residual_morans_i() gains a null argument, defaulting to "auto".
    Model residuals are not exchangeable — they are orthogonal to the
    design matrix — so the classical randomisation moments are wrong for
    them. At n = 120 with six smooth covariates and independent errors,
    OLS residuals had mean I = -0.031 against the exchangeable E[I] =
    -0.008, and the z-score averaged -0.54 with sd 0.90 instead of 0
    and 1. The Cliff & Ord (1981, sec. 8.3) regression-residual moments
    restore mean z = -0.09, sd 1.03 and a 4.3% rejection rate against a
    nominal 5%, and agree with spdep::lm.morantest() to machine
    precision (verified at 1e-16 through the public function). "auto"
    applies them only when the fit's residuals really are the OLS
    residuals on the rebuilt design, which a forest's and a working
    GWR's are not; the null actually used is reported in the return
    value.

  - summarize_by_cell() standard errors under a design effect were too
    small. s / sqrt(n / deff) corrects the mean's variance for
    clustering but leaves s^2 biased low by the same clustering: for
    exchangeable correlation rho, E[s^2] = sigma^2 (n - deff)/(n - 1).
    The two errors compound. Measured 95% CI coverage at n = 20: 0.905
    at rho = 0.3, 0.796 at rho = 0.6, 0.632 at rho = 0.8; after
    rescaling by sqrt((n-1)/(n-deff)), 0.952 / 0.952 / 0.953. Applies to
    deff = "kish", deff = "variogram" and a fixed numeric deff. The
    default deff = 1 path is bit-identical to before.

  - determine_optimal_levels() ranks model-aware candidates on the
    standardised deviate, not on |Moran's I|. E[I] and Var(I) both
    depend on the cell count, so |I| shrinks as k grows whether or not
    the finer tessellation captures anything. Over 300 replicates of a
    response with no spatial structure, mean |I| fell monotonically
    from 0.114 at k = 10 to 0.050 at k = 60 — an |I| ranking prefers the
    largest candidate for arithmetic reasons alone. Candidates are now
    ordered by |z| using the Cliff & Ord residual moments (exact here,
    since the cell-level residuals are OLS residuals by construction);
    over the same runs z had mean ~0, sd ~1 and a 5% rejection rate
    of 0.040-0.057 at every k. The "diagnostics" attribute now carries
    moran_z alongside moran_i.

  - estimate_sac_range() sweeps four azimuths, not two. A +/-22.5 degree
    window around 0 and 90 covers exactly 90 of the 180 distinct
    azimuths — every direction between 23 and 67 degrees, and
    between 113 and 157, fell into neither. On simulated fields with 3:1
    anisotropy and a true major-axis range of 300, the estimate came
    back at 255 and 249 for major axes at 0 and 90 degrees but 151
    and 147 at 45 and 135. Since make_folds(auto_range = TRUE) sizes
    blocks from this number, a diagonally oriented field silently got
    blocks half as wide as the correlation they were meant to separate.
    c(0, 45, 90, 135) tiles all 180 azimuths; the same fields now
    return 255 / 245 / 249 / 228. A direction whose variogram never
    reaches a sill is excluded rather than taken as a long range, and
    the directional attribute now has four named entries.

  - fit_bayesian_spatial_model()'s calibrated length-scale prior never
    reached Stan. brms::set_prior(spec, class = "lscale") with no coef
    is a global prior, and brms applies a global prior only to
    coefficients with no individual prior of their own — every lscale
    coefficient always has one. brms dropped it with a note and Stan
    received brms's defaults, which made gp_lengthscale_bounds(), the
    tail calibration and $info$gp_lscale_prior dead weight. Confirmed
    with brms::make_stancode(): the requested prior is absent under the
    global form and present under the coefficient-level form, which is
    now used. $info$gp_lscale_prior is read back from
    brms::validate_prior(), so it records what brms will actually use.

  - The GP basis was sized against the wrong domain measure. brms builds
    the boundary as choose_L(x, c) = c * max(1, max(x) - min(x)) over
    the pooled, column-centred covariates — the full range, not the
    per-axis half-range in which Riutort-Mayol et al. state their
    inequalities. Recovering the boundary from make_standata()'s
    eigenvalues confirms L = c * full range exactly at every c, so the
    old convention built a boundary twice as wide as gp_k was sized for:
    the GP was under-resolved, and $info$gp_ell_min — the diagnostic
    meant to catch exactly that — was twice too lenient to fire. The c
    floor is now brms's own default 1.25 rather than 1.2.

  - fitted() on a gwr_fit could return a coefficient surface. The search
    for GWmodel's fitted-value column matched the whole candidate name
    vector with %in% and took the first hit in the SDF's column order —
    and the local coefficients come first. A predictor named fit, pred,
    prediction, fitted or yhat therefore returned its own coefficient
    column, silently: executed in-sample R^2 was -1.18 against a
    true 0.981, and residuals(), summary(), model_metrics(),
    compare_models() and every cv_gwr() fold consumed it without a
    warning. The search now runs in preference order and excludes any
    candidate that is also a model term; all five colliding names now
    give R^2 = 0.981, identical to the renamed control.

  - coef.gwr_fit() returned GWmodel's whole SDF data slot — 15 columns
    for a two-predictor fit, of which 3 are coefficients and the rest
    are standard errors, t-values, the response, the fitted values,
    residuals and Local_R2. It now returns the model terms only; reach
    for object$engine$SDF for the rest.

  - estimate_sac_range() is reproducible, and no longer disturbs the
    caller's RNG. seed now defaults to 123L rather than NULL. The n_max
    subsample is an internal approximation, not part of the answer, and
    leaving it unseeded made the returned range differ between runs on
    identical input (19531 / 19589 / 19605 on three calls) while
    silently advancing the caller's stream — and make_folds(auto_range =
    TRUE) sizes its blocks from that number. Pass seed = NULL for the
    old behaviour.

  - estimate_sac_range() rejects a non-numeric response. as.numeric() on
    a factor returns its level codes, so a factor response produced a
    variogram of an arbitrary integer relabelling of the categories and
    the estimated range changed when the levels were reordered (3700
    against 2497 on the same data). Factors and character columns are
    now an error naming the column; logicals are read as 0/1.

  - Fold sets built from a different dataset are refused. Fold splits
    are lists of ..row_id values, and row IDs are seq_len(nrow()) unless
    supplied, so passing cv_gwr() a folds object built from another
    dataset of the same size applied cleanly — every ID matched, every
    fold was populated, and the model was scored on splits describing
    other observations. make_folds() now records a small
    projection-invariant row fingerprint in params$row_probe, and
    cv_gwr(), cv_bayes(), cv_spatial() and cv_rf() error rather than
    proceed. Fold objects from earlier versions carry no fingerprint and
    are passed through unchecked.

  - evaluate_insample() rejects duplicated names in fits. model is the
    key compare_models() joins its metric and Moran's I tables on, so
    two fits called "GWR" produced a 2x2 cross-join: four rows, every
    one carrying the first fit's numbers, with the second fit never
    scored at all.

  - fit_gwr_model() rejects a non-numeric predictor. gwr.basic() expands
    contrasts via model.matrix() and fits, but gwr.predict() does not
    and fails, so the model appeared to fit and then silently predicted
    all NA.

  - fit_gwr_model() no longer rejects a two-valued continuous response.
    The "binary" error is now gated on the response being integer-like.
    A left-censored or saturated measurement (every observation at a
    detection limit or a ceiling) has two distinct values and is
    perfectly continuous; it now warns instead. The guard also runs once
    per fold inside cv_gwr(), where a small training fold can
    legitimately hold only two distinct values.

  - fitted() returning the wrong length, or nothing, is now an error in
    summary() and model_metrics() rather than a plausible row count over
    an all-NA comparison. new_spatial_fit() is the documented extension
    point, so a subclass with a missing or mis-sized fitted() method is
    user-reachable.

  - The cached fitted() on a bayesian_fit is stamped with the n and a
    digest of the data it was computed from. The cache environment has
    reference semantics — which is what makes it survive copy-on-modify
    — so fit2 <- fit gave both objects the same cache, and assigning
    different data to the copy returned the original's values at the
    original's length.

  - make_folds() drops rows with empty or non-finite coordinates, with a
    logged warning naming the count, rather than letting an EMPTY POINT
    reach block_kfold's nearest-block rescue and die with "replacement
    has length zero".

  - When every fold fails, the warning now names the first underlying
    error. Previously "all 5 folds failed" was the whole diagnosis even
    when the cause was simply that brms or GWmodel was not installed.

  - make_folds() records the CRS the folds were built in as params$crs.
    Geographic input is projected by ensure_projected() to a CRS the
    caller never chose, and block_size, sac_range and buffer are lengths
    in that CRS.

  - .morans_i_for_k() returns NA at or below nine cells, where every
    cell neighbours every other and Moran's I collapses to exactly
    -1/(k-1) for any residual vector — a function of the cell count
    alone.

  - residual_morans_i(fit, k = 1) works on machines without FNN. apply()
    simplified the length-1 result to a vector, making the neighbour
    index a 1 x n matrix and every row after the first out of bounds.

  - summarize_by_cell(deff = "kish") under-estimated the predictor ICC
    by about a factor of m. The pooled one-way ANOVA grouped the m
    z-scored predictor columns under the same cell label, so independent
    per-column cell effects averaged away in the shared cell mean and
    the between-cell sum of squares shrank by ~1/m. Measured at true rho
    = 0.5 with m = 4: pooled ICC 0.12 against 0.495 per column, so every
    predictor SE was ~44% too small. The pooled group is now (variable,
    cell), which recovers 0.49.

  - Design effects are built from what a column actually observes. A
    cell of 10 rows with 2 finite responses had its response SE formed
    at the 10-row design effect, then applied to a 2-observation mean:
    adding 8 NA-response rows moved the SE from 8.46 to 26.38. Each
    column's design effect now uses its own non-missing count, with the
    mean pairwise correlation recomputed over the observed locations
    when a cell has NAs; cell_weight is the effective count of the
    primary variable, not of rows (n still counts rows).

  - CRS-less coordinates get ONE interpretation, wherever they enter.
    prep_model_data() assumed EPSG:4326 for CRS-less data that looked
    like lon/lat and projected it, while every predict() method passed
    the fit's CRS as a target — a branch that stamped it onto the raw
    numbers. The same rows sat in two different places, and predict(fit,
    newdata = training rows) disagreed with fitted(fit) by up to one
    response SD (R² 0.98 in-sample, 0.64 via newdata). The heuristic is
    now a single function used by both branches; a fit records the
    assumption it was built under and replays it on CRS-less newdata.
    Two further symptoms of the same split — CRS-less LINESTRINGs
    aborting in coerce_to_points() ("crs not found"), and hex/square
    build_tessellation() refusing input voronoi accepted — are fixed
    with it. The assumption is now announced with a real R warning.

  - residual_morans_i() refuses a malformed weights matrix instead of
    silently substituting the default k-NN(8) matrix (I = 0.874 returned
    for four malformed shapes against 0.805 for the weights actually
    supplied).

  - The fold-provenance fingerprint no longer refuses the caller's own
    data. Three defects in the version introduced last pass: a character
    ..row_id was coerced to all-NA (and matched row 1 everywhere);
    coordinates were compared as "%.7g" strings and flipped on the ~1
    in 5000 that a reprojection moved by 5e-9°; and polygon input was
    probed after pointization, so a different pointize in the cv call
    read as different data. Both sides now probe the geometry as
    supplied, keep IDs in their own type, and compare numerically
    within 1e-6°.

  - n_folds_attempted counts the folds supplied. A fold whose test rows
    were all removed as incomplete vanished before fitting and was
    absent from both counts, so five supplied folds reported 4/4. It is
    announced with a real warning.

  - determine_optimal_levels()'s elbow uses the signed deviation below
    the chord. abs() let a concave bump above the chord — a k where
    k-means fell into a worse local optimum — win with the same
    magnitude.

  - estimate_sac_range() returns NA for a constant variable (an exactly
    explained response, or a constant one) instead of a fitted "range"
    of 168 or 673 from a variogram that is identically zero.
    make_folds(auto_range = TRUE) no longer re-opens the unseeded
    subsample by forwarding its own seed = NULL.

  - fit_bayesian_spatial_model() attaches a user-supplied global lscale
    prior at coefficient level, the same way it does its own, so
    set_prior(..., class = "lscale") reaches Stan instead of being
    discarded.

  - fitted() on a gwr_fit returns the prediction when a predictor is
    named prediction. The model-term exclusion added last pass was
    applied to gwr.predict()'s SDF too, where coefficients are suffixed
    _coef and the column literally named prediction is the prediction;
    predict() returned all NA.

  - Smaller: a logical response meets the same binary-response guard as
    0/1; model_metrics() errors on a non-numeric response instead of
    returning n = 0; fitted.bayesian_fit() errors when the posterior
    cannot be drawn instead of returning silent NA;
    create_grid_polygons() refuses a grid above max_cells (default 1e6)
    up front; make_folds(method = "buffered_loo") states its guard in
    bytes (splits are ~4n² bytes; the old n = 20000 cap admitted 1.6
    GB); -0 and 0 are the same coordinate in the duplicate-aware k-NN;
    .morans_i_for_k() returns the NA pair whenever the moments are
    unavailable.

  - Documented warnings are now R warnings. Eight paths the manual
    described as warning only wrote a logger line, invisible to
    tryCatch(warning = ), expect_warning() and options(warn = 2):
    residual_morans_i() returning NULL, the CRS assumption in
    ensure_projected() and stamping in harmonize_crs(), GWR
    collinearity, seed clamping in voronoi_seeds_kmeans() /
    get_voronoi_seeds(), dropped rows in ensure_stable_poly_id(), and
    the dropped column in assign_features_to_polygons(). Three
    deliberate methodological cautions (include_coords = TRUE,
    random_kfold feature selection, non-standardised Moran weights) stay
    logged and their documentation now says so.

  - predict() on a Bayesian GP fit depended on which other rows shared
    the call. brms 2.x stores Xgp, dmax and cmeans in a fit's GP basis
    but not the Hilbert-space boundary L, so brms:::.data_gp()
    recomputes it from whatever rows predict() is handed. Every
    eigenfunction of the approximation therefore moved with the newdata
    bounding box while the fitted basis coefficients stayed put.
    Measured: L was 5.57 at fit time, 4.02 for a five-row newdata
    and 3.63 for one row; predict_surface() on the same 10,000-cell grid
    differed by 1.78 between chunk_size = 5000 (the documented default)
    and a single call; and cv_bayes(), which predicts each test fold
    separately, scored every fold against a basis the model was never
    fitted with. predict.bayesian_fit() now pins the boundary by
    appending the training coordinate extrema and dropping them again,
    so the chunked, fold-wise and single-call answers are identical.

  - A predictor whose name is not a syntactic R name fitted a different
    model. Every backend builds its formula from these names, so a
    column called "B5-B4" was fitted as B5 - B4 — a different model,
    silently, with $predictor_vars still reporting "B5-B4" — and
    "band 4" died inside str2lang() with a parser error naming no
    column. prep_model_data() now refuses such names up front and says
    which column and what to do about it. (Backticking is not a fix: the
    sp coercion GWmodel needs runs the names through make.names()
    anyway, after which the formula and the data disagree.)

  - summarize_by_cell(deff = "kish") weighted the response with a
    predictor's ICC. The fallback to the predictor ICC is for the case
    where no response was supplied; applying it whenever the response's
    own ICC came out non-positive meant that merely adding a predictor
    to the call changed the response's regression weight by a factor of
    20, while the response SE was (correctly) left at deff = 1.

  - A column literally named n was summarised as the row count.
    dplyr::summarise() makes each new column visible to the expressions
    after it, so the row count shadowed a response or predictor named n:
    resp_mean_n came back equal to n, with sd and se NA, and no message.

  - Variogram models are fitted with a nugget. A nugget-free model
    forces the curve through the origin, and gstat's default N/h^2
    weights buy that constraint by collapsing the range: with a 50%
    nugget the fitted range came back at about 0.45 of the truth, so
    make_folds(auto_range = TRUE) built blocks less than half the
    correlation length it reported.

  - The directional maximum is used only when anisotropy is established.
    Splitting 180 degrees four ways leaves each directional variogram
    about a quarter of the point pairs, and the maximum of four noisy
    estimates is biased upward: on isotropic simulated fields the
    returned range ran about 40% above the truth and the "notable
    anisotropy" warning fired on the majority of them. An
    omnidirectional variogram is now fitted alongside and is the default
    answer; the directional maximum is used when all four directions
    fit, their ratio exceeds 1.5, and the widest stands more than 1.5x
    above the all-pairs estimate. On the package's own test field (true
    range 80) the old rule returned 248 and the new one returns 84;
    genuine 3:1 anisotropy is still recovered.

  - predict() replays the CRS decision the fit was made under, including
    a negative one. When CRS-less training data were passed through as
    planar, nothing recorded that, so every predict() re-ran the lon/lat
    heuristic on newdata alone — and a subset of those same training
    rows, whose own bounding box sits inside the lon/lat envelope, was
    judged differently from the whole: taken for degrees, reprojected,
    and predicted about 1e6 m from where it was fitted. predict(fit,
    training_subset) disagreed with fitted(fit)[subset] by more than the
    response's standard deviation while predict(fit, full_data) agreed
    exactly.

  - Antimeridian data are no longer flattened onto EPSG:3857. A bounding
    box cannot distinguish global coverage from a layer straddling
    +/-180 degrees, but the coordinates can (one very large gap in the
    sorted longitudes). Web Mercator splits such a layer: two
    stations 41 km apart came out 40,068 km apart, destroying every
    distance downstream. Only genuinely global coverage now falls back
    to EPSG:3857.

  - The projection for a wide extent is chosen by measurement, not by
    rule of thumb. Albers standard parallels came from the bounding box
    while the conic/azimuthal branch came from the centroid latitude, so
    a trans-equatorial extent put the parallels either side of the
    equator: at lat_1 = -lat_2 PROJ refused the string outright and
    ensure_projected(), make_folds() and build_tessellation() all
    aborted with an internal-looking "invalid crs", and just short of it
    the projection distorted distances by 15.6% against the single UTM
    zone's 1.65%. The candidates are now scored by projecting a sample
    of the data's own points and comparing planar with geodesic
    distances — WGS84 ellipsoidal distances (Vincenty), not
    sf::st_distance()'s s2 sphere, whose 0.24–0.56% gap from the
    ellipsoid is the size of the errors being ranked and mis-ordered the
    candidates on 16 of 40 random wide extents; the message reports both
    figures.

  - fitted.gwr_fit() could return a coefficient surface. Which GWmodel
    call produced an SDF was inferred from its column names, so a
    predictor ending in _coef made a gwr.basic() SDF look like a
    gwr.predict() one and switched off the model-term exclusion; with a
    second predictor named yhat, fitted() then returned that predictor's
    local coefficient surface (R2 -2.28 against 0.986) and residuals(),
    summary() and model_metrics() followed it. The mode is now passed
    explicitly.

  - A logical predictor supplied as text was silently mis-coded.
    is.numeric, is.factor and is.character are all FALSE for a logical,
    so the type guard skipped it entirely: a "TRUE"/"FALSE" character
    column — what a CSV round trip produces — was factor-coded 1/2
    against splits built on 0/1, sending every row to the TRUE side
    (correlation 0.21 with the correct predictions).

  - compare_models_cv() scores every model on identical folds. The
    documented guarantee only held when folds was supplied: with folds =
    NULL each backend built its own from k, seed, block_size and the
    rest, so a per-model block_size put 104 of 150 rows in different
    folds for GWR and RF, and seed = NULL did the same with no overrides
    at all. The folds are now built once, and the arguments that decide
    the split are protected.

  - User weights with a non-zero diagonal. Every moment of Moran's I
    assumes no observation is its own neighbour. A self-inclusive
    row-standardised kNN matrix — an easy thing to build by hand —
    rejected the null on 75% of white-noise residuals at a nominal 5%,
    with no condition raised. The diagonal is now zeroed with a warning.

  - build_tessellation()$index no longer snaps outside points to the
    nearest cell. Points outside every cell were assigned to whichever
    was closest, so a summary built from the index counted all 40 points
    of a layer whose study area held 10, while
    assign_features_to_polygons() on the same cells correctly
    reported 30 misses. Only a point within a thousandth of a cell width
    is snapped now; the rest are NA, and the count is logged.

  - ensure_stable_poly_id() is stable across projections again. The sort
    key was the raw double centroid, so cells sharing an exact x in one
    CRS differed by ~1e-11 degrees after a round trip through another:
    14 of 16 cells got a different ID depending on which CRS the layer
    arrived in — the exact failure the function exists to prevent. The
    key is now rounded to 7 decimals (about a centimetre).

  - A numeric deff is applied as the uniform inflation it is documented
    to be. The E[s^2] correction the estimated design effects apply is
    derived from within-cell correlation and is unjustified for a
    constant the caller chose: deff = 2 doubled a 3-point cell's SE
    instead of multiplying it by sqrt(2), and returned NA for every cell
    with n <= deff.

  - The Kish ICC guard matches its documentation. The docs promise "at
    least 2 cells with 2+ observations; falls back to deff = 1
    otherwise", and the check was only k >= 2, N >= 4: all-singleton
    cells give a within-group sum of squares of 0 and therefore an ICC
    of exactly 1 — a design effect of n from data carrying no
    within-cell information at all.

  - gp_lengthscale_bounds() and the GP basis use two dimensions and
    unique locations. stats::dist() uses every column, so POINT Z
    geometry gave 3-D length-scale bounds in mixed units; and brms::gp()
    defaults to gr = TRUE and reduces its covariates to unique rows
    before taking the boundary, so replicated locations made the
    package's S — and with it gp_c and the basis-adequacy threshold
    gp_ell_min — wrong by that factor (0.68x with one heavily-sampled
    station).

  - estimate_sac_range() and make_folds() drop the Z dimension.
    gstat::variogram() and sf::st_distance() use every coordinate
    dimension, so an XYZ layer had its elevation folded into each lag:
    the returned "range" was a length in 3-D (413.6 against 136.8 for
    the same stations) while the block grid, the buffered-LOO buffer,
    NNDM's neighbour distances and summarize_by_cell() all work in 2-D
    map distance.

  - determine_optimal_levels() returns the elbow first. Under the
    geometric criterion — the default, and the fallback every
    model-aware call takes below the nine-cell floor — the candidates
    came back sorted ascending, so the knee sat in the middle: k[1] and
    top_n = 1, both documented as "the top-ranked candidate", returned
    knee − 1 on every such call, and the help example answered 1 for two
    clearly separated clusters. The vector is now knee first, then its
    lower and upper neighbours, so position 1 means the same thing on
    both paths. The quick-start data's answer moves from 3 4 5 to 4 3 5.
    Same code in 1.0.0.

  - estimate_sac_range() is invariant to rotating the layer. Two things
    were not. The lag cutoff was a fraction of the bounding-box
    diagonal, a property of the axes rather than of the points, which
    grows by up to √2 when the same layer is rotated 45°; it is now a
    fraction of the farthest pair (found on the convex hull), as the
    documentation always said. And the directional maximum is never
    preferred over a usable all-pairs fit. The 1.0.0 code returned the
    widest of two axis-aligned directions; the three hurdles later put
    in front of it (all four directions fitted, ratio above 1.5, maximum
    above 1.5× the all-pairs fit) still let the noise through — on one
    isotropic field rotated in 10° steps they "established" anisotropy
    in 14 of 18 orientations and returned ranges from 225 to 529,
    a 2.35× spread produced by nothing but the direction the axes
    pointed. The four directional ranges are still reported, as a
    diagnostic, in the directional attribute; anisotropy_used is TRUE
    only when the all-pairs fit itself failed. A field known to be
    anisotropic should size its blocks from max(attr(range,
    "directional")) explicitly, and the log line says so.
    
    The variogram fit itself no longer depends on gstat's single
    starting value. fit.variogram() starts the optimiser at a third of
    the longest lag; for a field whose range is a small fraction of the
    extent that start is ten times too long, and whether the iteration
    landed or collapsed to a singular model depended on floating-point
    details — the same 250-point field fitted on Linux and came back
    singular on an arm64 Mac, where the estimate then fell through to
    the directional maximum (315 against a true range of 80). Shorter
    and longer starting ranges are tried as well, the converged,
    non-singular fit with the smallest weighted sum of squares (gstat's
    own criterion) is kept, and a converged spherical fit is preferred
    over an exponential one that did not converge. When no model fits at
    all — a flat, nugget-only variogram — the NA now carries the
    empirical variogram, so plot(fit, type = "variogram") draws it and
    says why there is no range instead of refusing.

  - residual_morans_i() splits tied neighbour distances. The k-nearest
    neighbour weights broke ties at the k-th distance by whichever point
    came first — in row order on the dense path, in kd-tree order under
    FNN — so the same data in a different row order gave a different I,
    z and p whenever observations shared a location (repeat visits), and
    on gridded data the answer depended on whether FNN was installed.
    Points tied at the k-th distance now share the remaining weight
    equally, on both paths, which is the only rule that is a function of
    the geometry alone; the matrix is row-standardised as before. Where
    no distances tie the weights equal spdep's k-NN weights exactly; on
    repeat-visit data every co-located twin is retained with its share,
    and the statistic no longer changes when the rows are shuffled or
    when FNN is installed. Same code in 1.0.0.

  - fit_gwr_model()'s collinearity diagnostic is the scaled condition
    index, thresholded at 30. It was kappa() on the raw, unscaled
    predictor matrix at 1e6 — a number that depends on the predictors'
    units, so rescaling a column changed it and the threshold was a
    threshold on nothing: a design with a scaled condition index
    of 1322, whose local coefficients ran from −86 to +150 around a true
    value of 2, raised nothing. Belsley's index (every column scaled to
    unit length, intercept included, the ratio of the largest to the
    smallest singular value) is now computed for the global design and
    for each sampled local window, and the conventional 30 is the
    threshold in both. Expect the warning on designs that were silent
    before.

API and default changes

  - estimate_sac_range() returns NA instead of a number when the fitted
    range runs past the longest observed lag, or the optimiser stopped
    at its iteration limit. Such a range is unidentified, not long — the
    empirical variogram never reached a sill — and sizing blocks or
    design effects with it is worse than declining to. The refusal
    carries rejected_range and rejected_reason attributes and keeps the
    fitted variogram, so plot(type = "variogram") can draw exactly the
    case worth looking at; the result is classed sac_range, so it prints
    as a bare NA rather than dumping the fit to the console. Callers
    that fed the old number straight into make_folds(block_size = ) now
    need to handle NA — that is the point.

  - Removed the legacy wrappers evaluate_models(), evaluate_models_cv()
    and phi_prior_bounds(). Use compare_models(), compare_models_cv()
    and gp_lengthscale_bounds().

  - compare_models_cv() gains an "RF" branch and an rf_args argument, so
    a ranger forest can be compared against GWR and the Bayesian GP on
    identical folds. Unrecognised model names now raise a warning and
    are dropped, and a request with nothing recognised left is an error.
    Previously a bare intersect() discarded anything outside c("GWR",
    "Bayesian") and fell back to GWR, so models = "RF" silently ran GWR
    and reported it as the answer.

  - coef() on a spatial_fit now either returns coefficients or errors;
    it never returns NULL. coef.gwr_fit() and coef.bayesian_fit() used
    to return NULL on failure, which is indistinguishable from "this
    model has no fixed effects", so lapply(fits, coef) quietly produced
    a short answer. coef.rf_fit() errors as before — a forest has no
    coefficients; use fit$info$importance. See ?new_spatial_fit.

  - Every predict() method errors when the number of predictions does
    not match the number of rows that survived cleaning. It used to
    recycle silently: two predictions for four clean rows produced a
    four-row answer.

  - create_grid_polygons_cached()'s default type is now "square",
    matching create_grid_polygons(). The two disagreed, so cached and
    uncached calls built different grids from the same arguments.

  - build_tessellation() keeps cell_id on "hex" and "square" grids
    instead of deleting it after indexing, so all four methods return
    the same ID column and plot_tessellation_map(fill_col = "cell_id")
    works on a grid. poly_id is retained alongside it. params$expand now
    echoes the value passed rather than a hard-coded 0; the grid and
    triangle methods still ignore expand, which is now documented.

  - clip_target_for() projects lon/lat input before applying expand. The
    fraction-of-extent form was computed from a bounding box in degrees
    and handed to sf::st_buffer(), which reads dist as metres. The
    returned clip target is therefore in the projected CRS, not the
    input CRS, and says so.

  - Voronoi, grid and triangle clipping union the boundary first.
    Against a multi-feature boundary, st_intersection() split every
    straddling cell into one row per boundary feature and grafted the
    boundary's attribute columns onto the result.

  - predict.bayesian_fit(newdata = NULL) honours summary, type and
    draws. It short-circuited to fitted(), which caches epred column
    means and nothing else, so summary = "median" returned means and
    type = "predict" returned expected values — silently.

  - prep_model_data() accepts predictor_vars = character(0), making
    intercept-only spatial GP models reachable from
    fit_bayesian_spatial_model(). fit_gwr_model() and fit_rf_model()
    reject an empty set explicitly.

  - fit_bayesian_spatial_model(control = ) is merged over the package
    defaults (adapt_delta = 0.9, max_treedepth = 12) rather than
    replacing them. Passing list(max_treedepth = 15) silently dropped
    adapt_delta — the setting the divergence warning tells you to raise.

  - cv_bayes()'s predictive_coverage is averaged across folds weighted
    by each fold's n_pred. The per-fold values are means over that
    fold's test rows, so an unweighted average is not the pooled
    quantity once fold sizes differ — and block_kfold tolerates a 3:1
    imbalance before it even logs a warning.

  - The three seeding functions (get_voronoi_seeds(),
    voronoi_seeds_kmeans(), voronoi_seeds_random()) all emit seed_id and
    method columns, so they are drop-in interchangeable.
    get_voronoi_seeds() returns seeds in the boundary's CRS for every
    method; per-branch alignment had made the final alignment block
    unreachable.

  - Geometry-type checks require every geometry to be an accepted type,
    not merely one of them. A mixed POINT/POLYGON layer passed a
    POINT-only check.

  - ensure_projected() errors on a target_crs that does not resolve to a
    usable CRS, rather than returning the input unchanged and letting
    unprojected coordinates flow into distance and area computations.

  - fit_bayesian_spatial_model() derives the GP basis count (gp_k) and
    boundary factor (gp_c) from the ratio of the estimated length-scale
    to the domain size, rather than from the number of observations.
    
    brms::gp() builds a full tensor grid, so gp(..x, ..y, k = gp_k)
    carries gp_k^2 basis functions — gp_k is the count per dimension.
    The previous rule reduced to max(15, floor(sqrt(n))) for any n
    above 45, making gp_k^2 identically n: at n = 10,000 the model
    carried 10,000 basis functions and an n × n design matrix, at which
    point the approximation was no longer approximating anything.
    
    Across the scenarios in dev/baseline-structural.rds the derived
    value is 22–24 per dimension and largely independent of n: at n
    = 2,000 the basis count falls from 1,936 to 576, at n = 10,000
    from 10,000 to 529, and at n = 200 it rises, 15 to 23 — a
    correction, not an optimisation, so results move in both directions.
    gp_c was hard-coded at 1.5, too small whenever the length-scale
    exceeds roughly half the domain half-range; the derived value
    ranges 2.85–3.59 over the same scenarios. Pass gp_k and gp_c
    explicitly to restore the old behaviour. For scale, a 4-fold
    cross-validated fit at n = 2,000 with 2 chains × 1,000 iterations
    takes 1,186 s on the reference machine at cross-validated R² 0.927
    (dev/baseline-accuracy.rds, 2026-08-20); no comparable timing was
    captured for 1.0.0, so none is quoted.

  - The GP term is built with scale = FALSE, changing the default result
    of every Bayesian fit. brms::gp() otherwise rescales its covariates
    so the maximum pairwise distance is 1 and reports lscale in that
    space, while this package standardises the coordinates itself and
    expresses the length-scale prior, gp_c and the basis adequacy
    threshold in those units. The two normalisations differed by roughly
    the maximum pairwise distance (~4.9 for standardised 2D
    coordinates), leaving the automatic prior about five times too
    diffuse — a likely contributor to divergent transitions and rejected
    initial values. There is now exactly one coordinate scaling.

  - The GP fits one length-scale per coordinate axis (gp_iso = FALSE), a
    second change to default Bayesian results. Coordinates are
    standardised per axis, so a single shared length-scale made the
    kernel anisotropic in the original CRS by the ratio sd(X)/sd(Y) — a
    property of the sampling layout, not of the process. Pass gp_iso =
    TRUE for the previous behaviour; cost is unchanged, since the tensor
    grid is gp_k^2 either way.

  - The automatic GP length-scale prior is a calibrated inverse-gamma
    rather than normal(0, sd), a third change to default Bayesian
    results. A half-normal on a positive parameter puts its mode at
    zero, so most of its mass sat at length-scales shorter than the
    basis can resolve — where the Hilbert-space approximation develops a
    funnel and the sampler diverges. The replacement pins 1% of its mass
    below the estimated lower bound and 1% above the upper. The two tail
    conditions have one exact solution — a one-dimensional root in the
    shape — and that is how it is found, so the calibration succeeds for
    any bounds with upper > lower; degenerate bounds fall back to the
    half-normal with a logged note. The prior applied is recorded in
    $info$gp_lscale_prior.

  - ensure_projected() no longer forces continental-extent data into a
    single UTM zone. Transverse Mercator scale error grows quadratically
    with distance from the central meridian, so data spanning the
    contiguous United States carried distance errors of roughly 7.5%
    near the extent edge, propagating silently into
    estimate_sac_range(), make_folds(block_kfold) block sizing, GWR
    bandwidth selection and the GP length-scale. Extents reaching more
    than 5° from the candidate zone's central meridian now receive an
    equal-area projection centred on the data, with a logged
    explanation. Only longitude offset triggers the switch, since
    cos(lat) shrinks the distance from the central meridian and a tall
    narrow north-south extent is UTM's design case. (Which equal-area
    projection is no longer decided by latitude band — see The
    projection for a wide extent is chosen by measurement below, which
    also replaced the EPSG:3857 fallback for wide bounding boxes with
    antimeridian detection.) Pass target_crs to override.

  - Core counts follow the session's mc.cores opt-in, and are capped.
    fit_bayesian_spatial_model()'s documented default was cores =
    max(1L, parallel::detectCores() - 1L), and cv_*(parallel = TRUE)
    auto-detected the same way with no cap — 63 workers on a 64-core
    host, and a hard error wherever _R_CHECK_LIMIT_CORES_ is set. The
    Bayesian default is now getOption("mc.cores", 1L); the auto-detect
    path is capped by that option when it is set; and every worker
    count, explicit or not, is capped at the machine's core count (with
    a message) and at two under R CMD check.

Guards, messages and stricter input handling

  - prep_model_data() now drops rows whose geometry is empty or whose
    coordinates are not finite, and counts them in its existing log
    line. st_geometry_type() calls an EMPTY POINT a "POINT", so nothing
    ever looked at the coordinates: such a row reached GWmodel as a raw
    sp coercion error naming no row, ranger(include_coords = TRUE) as
    "Missing data in columns", and brms at predict time as an infinite
    GP boundary that made every basis function on the surface NaN. It
    also refuses a response listed among its own predictors, which was
    leakage in the forest, a silently reduced model in GWR, and
    duplicated rows plus a phantom <none> entry in the GWR selection
    table.

  - make_folds() validates k. A non-integer k truncated inside
    rep(floor(n/k), k): only floor(k) folds were built, the last rows of
    the permutation landed in no test set, and k was echoed back
    unchanged, so length(folds) != k.

  - make_folds(block_kfold) refuses a grid it cannot build. There was no
    cap on nx * ny, so block_size in the wrong unit asked for 1e8-1e11
    cells and exhausted memory; the message now names the implied cell
    count, the extent and the CRS units, as create_grid_polygons()
    already did. predict_surface() gained the same guard on cell_size
    and n_cells.

  - Hand-built folds are checked. Train and test overlapping is not
    cross-validation — the model is fitted and scored on the same rows
    and the result is reported as a CV score (RMSE 0.50 against 0.97 for
    the same data properly split) — and is now refused; fold IDs that
    name no row were dropped silently by na.omit(match()) and are now
    counted and logged.

  - The fold provenance probe tolerates a row it could not measure. It
    was taken before the empty-geometry filter, so it carried NaN
    coordinates and every later cv_*() call on the same data died with
    R's internal "missing value where TRUE/FALSE needed".

  - estimate_sac_range() drops rows with unusable coordinates instead of
    returning NA for the whole layer under the message "variogram model
    fit failed", which blamed the fit rather than the row and disagreed
    with make_folds() on the same data.

  - GWR's local collinearity spot-check now runs after the bandwidth is
    chosen, includes the intercept column, and treats a non-finite
    condition number as extreme. On the default path (bandwidth = NULL)
    it used a stand-in window, so the documented warning never fired;
    without the intercept it could not see an indicator constant inside
    a window; and is.finite(cn) && cn > 1e6 discarded exactly singular
    designs. A new post-fit warning counts local regressions that
    returned non-finite coefficients — previously fitted(), summary()
    and model_metrics() silently reported metrics computed from the
    survivors (n = 18 of 200, R2 = 0.96).

  - fit_gwr_model() warns when a fixed bandwidth is implausibly small
    for the data's extent. The bandwidth is a distance in the CRS the
    fit runs in, which prep_model_data() may have chosen: 0.2 supplied
    for lon/lat data is 0.2 metres, and every local window came back
    empty with nothing raised. The argument's documentation now says so.

  - summary.spatial_fit() applies the same response-type guard as
    model_metrics(). A character response still produced n = 0 and
    all-NA metrics there, and a factor died inside abs().

  - print.spatial_fit() prints the CRS (which its documentation has
    always promised) and one Formula line. sprintf() vectorised over a
    multi-element deparse(), so every bayesian_fit printed its formula
    as two mangled fields — the gp() term this package builds is always
    long.

  - summarize_by_cell() warns when cells_sf carries duplicated IDs (each
    summary row is then repeated per matching cell, so sum(n) exceeds
    the number of points), and realigns every per-cell vector in
    deff_applied after the join, not only deff: $rbar was left in
    pre-join order, so deff[i] and rbar[i] described different cells.

  - assign_features_to_polygons() warns when no feature falls in any
    polygon instead of returning an empty layer silently.

  - create_grid_polygons(target_cells = ) builds square cells for type =
    "square". cellsize = c(w/nx, h/ny) forced an exact bbox tiling, so
    the cells were rectangles (aspect 10.5 on a 1000:1 strip). For type
    = "hex" a differing cellsize[2] is now collapsed with a warning
    before the max_cells estimate, which used both components and was
    therefore off by their ratio.

  - create_grid_polygons_cached()'s renumbering is documented: it
    applies ensure_stable_poly_id() and create_grid_polygons() does not,
    so the same cell carries a different poly_id from the two builders.

  - get_voronoi_seeds(method = "kmeans") clusters in two dimensions and
    drops rows with unusable coordinates, matching
    voronoi_seeds_kmeans(). It used every column of st_coordinates(), so
    POINT Z geometry was clustered in 3-D with elevation dominating, and
    an EMPTY POINT crashed inside kmeans().

  - determine_optimal_levels() refuses a factor or character response,
    which as.numeric() silently turned into level codes: re-ordering the
    levels of the same factor changed the chosen number of levels and
    every moran_z.

  - plot_tessellation_map() brings CRS-less layers into the plot's CRS
    instead of passing them through to fail inside ggplot_build() at
    print time with sf's message, naming no layer; and it borrows a CRS
    from an overlay when the tessellation itself has none.

  - Documentation corrected where it did not match behaviour:
    ensure_projected() states how it chooses a projection for a wide
    extent (by centroid latitude when that entry was written, by
    measured distortion in this release) and that the choice is
    announced rather than silent; .looks_like_lonlat()'s two tests are a
    disjunction and the extent test decides first, so a small planar
    survey inside the lon/lat envelope IS taken for degrees — the trade
    and its reasoning are now stated; harmonize_crs() no longer claims
    to match ensure_projected() while doing something else;
    summarize_by_cell() states which estimand its standard errors are
    for (the grand mean, where measured coverage is 0.95, not the cell's
    own mean, where the naive SE is the better estimate) and that the
    variogram path applies one correlation function to every column.

  - A misspelt newdata is an error, not an in-sample answer.
    model_metrics(), evaluate_insample() and compare_models() forward
    ... to predict(), which checks it only on the out-of-sample branch,
    so model_metrics(fit, newdta = hold) silently took the in-sample
    branch and returned an RMSE of 1.086 where the held-out answer was
    25.24, with the same return shape. Arguments in ... with no newdata
    are now refused by name. predict() on a gwr_fit or a bayesian_fit
    likewise refuses unknown arguments instead of swallowing them;
    predict.rf_fit() accepts only ranger's own predict arguments through
    ....

  - predict() enforces one newdata contract on all three fit classes. A
    bare sfc died inside two of them with R's "argument must be
    coercible to non-negative integer"; a numeric-at-fit predictor that
    arrived as character (a CSV round-trip) was refused by name by
    rf_fit and bayesian_fit and returned all-NA with a generic backend
    warning from gwr_fit; a missing column was reported by two different
    functions in two wordings. All three now run the same check first,
    so the message is the same whichever fit is behind it.

  - make_folds() validates block_size the way it validates k: a single
    finite positive number, else an error naming the argument. NA and a
    length-2 vector used to die as internal R errors, and a negative,
    zero or character value was silently ignored — yet echoed back in
    params$block_size as if it had been used. The grid-size guard also
    formats its own message: a block_size in the wrong unit could ask
    for a grid past 2³¹ cells on a side, which %d refused with "invalid
    format" instead of the documented refusal.

  - estimate_sac_range() says which variable it modelled. A predictor
    name absent from the data was dropped silently by intersect(), and
    with every name unknown the raw response was modelled — so
    make_folds(auto_range = TRUE) sized blocks from a range of 88.5
    instead of the residual range of 362.2 with nothing said. An unknown
    predictor is now an error, as it is everywhere else; a detrending
    fit that fails (an all-NA column) raises a warning and falls back to
    the raw response; and a new detrended attribute records which was
    used.

  - make_folds(block_kfold, boundary = ) refuses a boundary containing
    none of the points — almost always two layers in different places, a
    CRS that could only be stamped — and raises a warning counting the
    points that fall outside a boundary that contains some, since the
    region is silently extended to cover them. The single-block error
    names what produced the grid (the block size, block_nx/block_ny, or
    the automatic grid) rather than always blaming block_nx/block_ny,
    which the caller may never have passed.

  - Rows that no fold names are reported. The folds ↔ data guard was
    one-directional: fold IDs naming no row were counted and dropped,
    but rows in the data that appear in no fold's train or test set
    passed with no condition at any level — a folds object built on
    site[1:45, ] and applied to all 90 rows scored 45 of them and
    reported n_folds_attempted = n_folds_succeeded = 3. Every cv_*() now
    raises a warning with the count and an example row ID.

  - User-facing conditions name the function the user called. The cross-
    validation path leaked two internal names into ordinary console
    output — .remap_folds(): and .cv_run_folds(): — and cv_rf(), a
    wrapper around cv_spatial(), reported every message, warning and
    error in cv_spatial()'s name. The sf-input assertion shared by the
    tessellation and seeding functions said "Expected an sf object" with
    no function named; it now names the caller and, when handed a whole
    build_tessellation() result, says to pass its $cells.
    coerce_to_points(mode = "line_midpoint")'s MULTILINESTRING refusal
    is prefixed like every other.

  - clear_grid_cache(cache_env = ) removes only its own entries. It
    removed every binding in the environment it was handed and counted
    them all as "entries removed", so a user who passed a project
    environment lost unrelated objects. Cache keys now carry a
    spatialkit_grid:: prefix and nothing else is touched.

  - summarize_by_cell(deff = "variogram") now uses every structured
    component of a nested variogram model, each weighted by its partial
    sill, which is the correlation the model implies (1 - gamma(h) /
    sill). It read the single largest component, so a user-built Nug +
    Exp + Sph model gave a correlation of 0.108 at 200 m where
    gstat::variogramLine() implies 0.197, and the design effects and
    standard errors with it. Models from estimate_sac_range() are
    single-component and are unaffected. A model of a family the
    function does not implement (Matern, power, circular, ...) was
    silently read as exponential; it now falls back to deff = 1 with a
    warning that names the family. Exponential, spherical and Gaussian
    are supported.

  - cv_bayes()$predictions$yhat_sd is the posterior predictive standard
    deviation of each held-out row, from the same draws that give the
    coverage columns. It was an unconditional NA placeholder. It stays
    NA when compute_pred_intervals = FALSE or the draws failed for a
    fold, and it is documented.

  - fit_gwr_model() refuses n <= p + 1 observations with one error,
    before touching the backend. Such a fit cannot have a residual
    degree of freedom in any window; n = 2 used to warn three times
    ("only 2 observations", "fallback bandwidth", "2 of 2 local
    regressions singular") and return a fit whose fitted values were all
    NA.

  - create_grid_polygons_cached() gains max_entries (default 50): once
    the cache holds that many grids, adding one evicts the
    earliest-added. It never evicted, so a loop over a thousand
    boundaries held every grid (about 2 MB per 2,500 cells) for the life
    of the session. Nothing but the grids is written into a
    caller-supplied cache_env; the insertion order lives inside the
    package. The cache key also hashes the package version, so a cache
    that outlives an upgrade cannot serve a grid built by an older
    create_grid_polygons().

  - spatialkit_quiet() accepts a logger threshold as well as TRUE/FALSE,
    and the value it returns can be passed back: old <-
    spatialkit_quiet(); spatialkit_quiet(old) restores exactly the level
    that was in force. spatialkit_quiet(FALSE) put back the package
    default (WARN) whatever had been set, and the returned value was
    refused as "must be TRUE or FALSE".

  - .onUnload() disarms both logger appenders, so the spatialkit logger
    namespace no longer keeps pointing at the session's temp-file path
    after unloadNamespace("spatialkit"). .onLoad() re-registers them.

  - residual_morans_i() documents the second reason "residual" (and
    "auto") falls back to the randomisation null: fewer than four
    residual degrees of freedom, where the residual variance formula
    divides by (n - p)(n - p + 2). "residual" logs a warning when it
    does, "auto" does not, and df is then n - 1.

  - estimate_sac_range() states what "effective range" is for each model
    -- three times the range parameter for the exponential fit (95% of
    the sill), the range parameter itself for the spherical fallback
    (100%) -- and its return-value documentation matches the code: the
    no-model case (both fits singular) returns the classed NA with
    rejected_reason set, and only the cannot-even-start cases (no gstat,
    too few values, no variance, a degenerate extent) return a bare NA.
    Its example now demonstrates a fitted range on a simulated field
    with a known one (3 x 100 = 300), and a refusal on a field whose
    range the data cannot pin down.

  - make_folds() documents what block_multiplier does (the automatic
    grid aims for block_multiplier * k blocks, so each fold holds out
    about that many; 3 is a compromise, not a published constant), cites
    Roberts et al. (2017) and blockCV (Valavi et al. 2019) for sizing
    blocks from the autocorrelation range, and notes that blockCV takes
    the fitted variogram's range parameter where auto_range takes the
    effective range -- three times that parameter for an exponential
    fit, so larger blocks. phi for method = "nndm" is explained as Mila
    et al. (2022) define it: the autocorrelation range beyond which
    matching is unnecessary, which estimate_sac_range() supplies.

  - summarize_by_cell(deff = "kish") says which ICC estimator it is
    (ANOVA with Donner's n0, not REML) and how far the two can differ on
    an unbalanced draw, so the difference is not read as a defect;
    area_of_applicability()'s training-DI sentence now says what the
    code does (each fold's actual training rows, not "everything outside
    the fold"); ?spatialkit no longer lists coef() among the methods all
    three backends share (a forest has none); the internal elbow helper
    no longer claims to match Kneedle, which it does not on shouldered
    curves; the GP basis diagnostic's code comment attributes its 10%
    posterior-mass trigger to this package rather than to Riutort-Mayol
    et al. (2023).

  - summary() on a fit prints R^2 and Adj R^2 in ASCII with aligned
    labels (the superscript two rendered as R<U+00B2> on non-UTF-8
    consoles, and Adj R²= had no space); print() on a random forest
    likewise; the six console messages that carried an em dash use --.
    coef() on an rf_fit prefixes its error coef.rf_fit(): like its
    siblings.

  - get_voronoi_seeds(method = "kmeans") sizes its candidate cloud in
    double precision; 50L * as.integer(n) overflowed to NA
    above 42,949,672 seeds.

  - README: the opening leakage example says it needs ranger; the
    no-viable-models example is assigned so that it does not print every
    fold table on a machine that has the backends; the auto_range and
    cv_bayes() failure examples show every line the console actually
    prints; the installation table no longer suggests installing loo
    separately (brms installs and calls it); the "logged note" from
    determine_optimal_levels() is identified as an INFO-level line in
    the session log file, not console output; the roxygen2 sentence no
    longer names a version or a RoxygenNote field.

  - The memory-guard test for make_folds(block_kfold) stubs out
    sf::st_make_grid() for the two refused calls, so a regression of the
    guard fails fast by name instead of attempting an 8000 x 8000 grid.

  - Every cv_*() refuses a fold list whose splits carry no train/test
    element, with an error that names the problem. It read f$train /
    f$test straight, got NULL for both, and built empty folds -- so the
    row-coverage warning fired and blamed folds "built on a different or
    subsetted layer", which was not the cause, and the run returned an
    all-NA overall with n_folds_succeeded = 0. area_of_applicability()
    had always refused the same input by name; the two now agree. When
    the splits look positional (two unnamed vectors each) the error says
    so and points at the fold label vector instead.

  - The folds argument of every cv_*() documents all three shapes it has
    accepted since the label vector was added earlier in this pass -- a
    make_folds() result, a list of list(train =, test =) splits, or a
    vector of fold labels -- where the help listed only the first two.
    The label vector is what makes folds from another package usable
    directly: blockCV::cv_spatial() returns one as $folds_ids (its
    $folds_list holds two unnamed vectors per fold, which is the shape
    now refused by name).

  - MAPE and SMAPE are documented as what they are: averages over the
    rows whose denominator is non-zero. Both have a denominator that can
    vanish -- MAPE where the observation is zero, SMAPE where
    observation and prediction are both zero -- and each drops those
    rows rather than returning Inf, which is the right arithmetic but
    was reported nowhere. On a response taking exact zeros (counts,
    rainfall, abundance) the consequence is material: with 62 zeros out
    of 120, MAPE is an average over 58 rows presented as though it
    covered 120, and SMAPE drops precisely the rows a well-fitted model
    got right, so it reads worse than the fit deserves. The new
    "Percentage errors on responses with zeros" section on
    model_metrics() -- inherited by evaluate_insample(),
    compare_models(), compare_models_cv(), summary() and all four cv_*()
    -- says so, notes that the n column is the finite-pair count and not
    the row count either percentage error used, and points at
    RMSE/MAE/R-squared (and, for a Bayesian fit, CRPS and interval
    coverage) as the metrics unaffected by it. No computed value
    changes; returning the per-metric row count would alter the metric
    frame's column set and is deferred.

Bug fixes

  - Data carrying no CRS works again throughout. ensure_projected() now
    rejects a target_crs that does not resolve to a usable CRS
    (previously a typo silently made the call a no-op), but internal
    callers derive that target from another object —
    st_crs(training_data) — and that object is allowed to have no CRS.
    Passing NA_crs_ through turned every CRS-less workflow into a hard
    error: predict() on all three backends, make_folds(method = "nndm"),
    make_folds(boundary = ), prep_model_data(boundary = ) and
    predict_surface(). cv_rf() was worse than an error — the per-fold
    predict() threw, so every fold "failed" and $overall came back with
    n_pred = 0 and RMSE = NA behind a generic warning. Internal call
    sites now pass NULL ("choose one automatically") when the source has
    no CRS; the user-facing validation is unchanged.

  - build_tessellation(crs = ), create_voronoi_polygons(crs = ) and
    create_grid_polygons(crs = ) errored with sf's "cannot transform sfc
    object with missing crs" whenever the input had no CRS — exactly the
    users most likely to pass crs =. Reprojection is impossible there,
    but assumption is not: the target CRS is now stamped on with a loud
    warning, matching what ensure_projected() already documents. Input
    that does carry a CRS is still reprojected, not relabelled.

  - compare_models_cv() built its argument list with c(list(...),
    rf_args), so any gwr_args/rf_args entry whose name collided with one
    the function sets itself produced two entries of that name and
    do.call() died with "formal argument 'seed' matched by multiple
    actual arguments". Since cv_rf() has both k and seed as formals,
    rf_args = list(seed = 3) — straight from the documented usage — was
    enough to trigger it. Extras now replace base entries by name.
    data_sf, response_var, predictor_vars and folds are protected and
    dropped with a warning, because a per-model override of those would
    silently make the models incomparable.

  - compare_models() given a single bare spatial_fit reported all-NA
    Moran's I columns and logged "'fit' is not a spatial_fit object"
    once per component. A spatial_fit is itself a list, so it passed the
    is.list() check and the loops then iterated the fit's own components
    as though they were models. It is now wrapped into a one-element
    named list, exactly as evaluate_insample() already did.

  - residual_morans_i() failed on its own documented fast path. With FNN
    and Matrix installed the weights are a sparse Matrix, and
    base::crossprod() does not S4-dispatch on the dgeMatrix that W %*%
    resid produces, so the call died with "requires numeric/complex
    matrix/vector arguments" — taking compare_models(), which calls it
    automatically, down with it. Rewritten as sum(resid_c * (W %*%
    resid_c)), which is numerically identical and uses only dispatching
    primitives. determine_optimal_levels() carried the same bug.

  - residual_morans_i() no longer errors on constant non-zero residuals:
    the degeneracy guard tested the raw sum of squares where Moran's I
    is a function of the centred residuals, so VI came out NaN and if
    (VI > 0) raised "missing value where TRUE/FALSE needed". A
    non-finite VI is handled too.

  - .build_knn_weights()'s n > 5,000 guard tests for both FNN and
    Matrix. Keyed on FNN alone, an unbounded dense n × n allocation went
    through whenever FNN was present but Matrix was not.

  - assign_features_to_polygons() drops columns of features_sf that
    would collide with the polygon ID column, with a logged warning.
    sf::st_join() suffixed them (poly_id.x / poly_id.y), which defeated
    the rename afterwards and left the result with no rows — reachable
    simply by re-assigning already-assigned points. A join that still
    fails to produce the ID column now errors and names the columns it
    did produce.

  - summarize_by_cell() keeps the "deff_applied" attribute when cells_sf
    is supplied; dplyr::left_join() rebuilds attributes from its x
    template and dropped it. The per-cell vector is remapped onto the
    joined row order, NA for cells holding no observations.

  - summarize_by_cell() joins on the native ID type when both sides
    agree. Coercing unconditionally made the returned ID type depend on
    an unrelated argument and turned integer IDs into "1", "10", "2", ….
    A genuine class mismatch still coerces both to character and logs
    why.

  - summarize_by_cell() coerces non-POINT geometry before computing a
    variogram-based design effect. sf::st_coordinates() returns one row
    per vertex, so a POLYGON or multi-vertex MULTIPOINT feature
    misaligned the coordinate matrix with the data and fed the wrong
    points into every cell.

  - make_folds(method = "buffered_loo") errors when the buffer excludes
    so much of the data that no fold retains two training points. Those
    folds used to sail through and be dropped one at a time inside the
    CV loop, so the only symptom was a generic "all folds failed"
    warning at the very end.

  - make_folds(method = "block_kfold") refuses a block size yielding a
    single block covering the whole extent — one fold with an empty
    training set, reported as a run that merely happened to score NA. An
    accepted autocorrelation range could trigger it:
    estimate_sac_range() rejects ranges above half the bounding-box
    diagonal while block construction needs half the width.

  - make_folds() coerces MULTIPOINT geometry rather than merely
    accepting it, for the st_coordinates() reason above; every fold was
    misaligned silently.

  - Cross-validation no longer renumbers folds. .remap_folds() dropped
    unusable folds from a list, shifting every later fold's index, so
    fold_metrics$fold and predictions$fold stopped lining up with
    make_folds()$assignment$fold. The original index is carried through.
    Folds left with fewer than two training rows are detected there and
    logged, instead of failing one at a time deeper in.

  - cv_spatial() rejects a fit_fn whose predict() returns the wrong
    number of values. Both the metric computation and the prediction
    frame recycled silently, so two predictions against four test rows
    yielded a four-row frame with metrics computed against fabricated
    pairs.

  - Cross-validation under parallel = TRUE reports a fold that died in a
    worker. parallel::mclapply() returns a try-error rather than NULL,
    which the NULL filter kept, and the failure surfaced as "subscript
    out of bounds". conditionMessage() has no method for a try-error, so
    the diagnostic branch itself threw; the condition is now taken from
    the object's attribute.

  - Cross-validation under parallel = TRUE is reproducible from seed and
    gives results identical to parallel = FALSE. .cv_run_folds() called
    parallel::mclapply() without seeding the fork streams, so each
    worker seeded itself from the clock and process ID. One seed per
    fold is now drawn in the parent, making each fold's stream a
    function of (seed, fold index) alone.

  - estimate_sac_range() rejects a singular variogram fit.
    gstat::fit.variogram() signals failure by setting attr(.,
    "singular") and returning normally, so testing only for a try-error
    made the spherical fallback unreachable and let a singular fit's
    range flow out as the estimated autocorrelation range — which
    make_folds(auto_range = TRUE) then sizes spatial blocks from.

  - .extract_gwr_values() requires every model-matrix column to match a
    column of GWmodel's SDF before multiplying the local coefficients
    through. A partial match reconstructed a linear predictor missing
    one or more terms and returned it as the fitted value — plausible
    numbers that were simply wrong, feeding fitted(), residuals(),
    summary() and every metric with no warning. A non-numeric
    coefficient column is refused rather than coerced.

  - fit_gwr_model() separates the three degenerate response cases.
    Folded together, an all-dropped dataset was reported as "binary (0
    unique values)" and a constant response as "binary (1 unique
    value)", while a genuinely binary non-integer response (1.5 / 2.5)
    failed the integer-like gate and passed unremarked.

  - fit_gwr_model() and gwr_model_selection() validate bandwidth.
    Unvalidated, NA gave "missing value where TRUE/FALSE needed", a
    length-2 vector gave "the condition has length > 1", and with
    adaptive = FALSE a zero or negative distance reached GWmodel
    untouched.

  - fit_gwr_model()'s local-collinearity spot-check no longer touches
    the RNG at all. It sampled its 30 locations from the global stream
    and fires only when n > 30 with at least two numeric predictors, so
    the same script produced different fold assignments depending on how
    many predictors a model happened to carry; cv_gwr() calls it once
    per fold. The 30 locations are now evenly spaced ranks of the
    observations ordered by x, then y -- reproducible, independent of
    the row order, and drawing no random numbers.

  - predict.gwr_fit() returns an all-NA vector when every row of newdata
    is dropped as incomplete, matching the other two backends, rather
    than surfacing a raw sf-to-Spatial coercion error.

  - predict.bayesian_fit() transforms newdata to the training CRS before
    cleaning it, and derives the surviving rows from one sentinel column
    instead of a second, separately-maintained copy of the cleaning
    rules. It errors when a predictor standardised at fit time is absent
    from newdata or has arrived as character — silently skipping it
    handed brms an unscaled column against a model fitted on a scaled
    one. Its failure path returns a matrix when draws = TRUE, honouring
    the documented return shape.

  - plot() on a spatial_fit errors when there are no finite residuals,
    instead of producing a uniformly grey map from limits = c(Inf,
    -Inf). A perfect fit is handled too: all-zero residuals gave limits
    = c(0, 0), a degenerate diverging scale whose breaks collapse onto
    one value.

  - plot_tessellation_map() logs a warning for a fill_col that is not
    present, instead of drawing an unfilled outline map with nothing to
    say anything had gone wrong — a mistyped label_col already warned.
    xlim/ylim are validated, and the theme default moved out of the
    formals so a Suggests package never appears in an exported
    function's default arguments.

  - harmonize_crs() announces when it stamps a CRS rather than
    reprojecting. sf::st_set_crs() only relabels; the coordinates do not
    move. ensure_projected() already made that assumption loudly.

  - coerce_to_points() rejects an EMPTY LINESTRING rather than
    misaligning the result. st_line_sample() yields no midpoint for one
    (and segfaults in sf 1.0.x), so the sampled midpoints stopped
    corresponding 1:1 with the rows they are scattered back into. A
    count check backstops any other divergence.

  - evaluate_insample() errors on an unnamed list. The loop is over
    names(fits), so an unnamed list iterated zero times and returned
    NULL silently; compare_models() then died in seq_len(nrow(...))
    nowhere near the cause.

  - determine_optimal_levels() coerces MULTIPOINT geometry rather than
    admitting it, and errors on a factor or character predictor by name
    instead of dying inside colMeans() with "'x' must be numeric".

  - create_grid_polygons() passes both cellsize and n to
    sf::st_make_grid() when both are known; st_make_grid() does not
    ignore n in the presence of cellsize for square grids, and omitting
    it made sf recompute nx = ceiling(w / cellsize), which
    floating-point division pushes one past the intended count. n is
    parsed and validated once, up front, instead of being silently
    coerced to NULL in one branch and erroring in the other.

  - The grid cache key no longer truncates target_cells. as.integer()
    made 25.2 and 25.7 collide on one key, so the second call silently
    received the first one's grid, and a NULL target_cells collapsed
    paste0() to character(0), crashing the lookup.

  - build_tessellation() normalises a CRS-less points_sf to NULL rather
    than NA_crs_, which is a list and so was not treated as "no CRS
    supplied" downstream. Hex and square grids are built in the points'
    CRS, so the grid and the points no longer end up in different CRSs
    and break the point-to-cell index.

  - build_tessellation(method = "triangles") triangulates the point set
    when geometry is unavailable, via sf::st_triangulate() on the
    unioned points. The fallback previously triangulated the convex hull
    polygon, discarding every interior point. The result is still the
    Delaunay triangulation of the input; only the resolution of
    degenerate configurations can differ from qhull's, and the logged
    warning now says so.

  - ensure_stable_poly_id() logs a warning naming the geometry types
    when it drops non-polygonal rows, which it silently did before.

  - voronoi_seeds_kmeans() and voronoi_seeds_random() validate their
    inputs (voronoi_seeds_random() also accepts the sfc its
    documentation always promised), and clamping k to the number of
    distinct positions is logged. get_voronoi_seeds(method = "provided")
    logs a warning when n disagrees with nrow(seeds), which it ignores.

  - gp_lengthscale_bounds() validates coords_xy and q_small. A vector
    coords_xy failed inside .safe_dist() with "argument is of length
    zero" and an out-of-range q_small inside quantile(), neither naming
    the argument.

  - fit_bayesian_spatial_model() validates the response before handing
    it to Stan, where nothing points back at the column, and validates
    gp_k, gp_c and control. The inverse-gamma prior is written with
    %.10g rather than %.6f: a small scale rounded to the literal
    "0.000000" and Stan rejected inv_gamma(a, 0) from deep inside the
    model block. Tightly clustered coordinates get there. The
    half-normal fallback's scale is guarded the same way.

  - compare_models_cv() names, in a warning, any gwr_args entry it
    drops. cv_gwr() has no ..., so entries meant for fit_gwr_model()
    alone (e.g. longlat) were discarded silently and simply had no
    effect.

  - .compute_reg_metrics() errors on a y_train_mean that is neither a
    scalar baseline nor one value per observation, instead of recycling
    it against the filtered response and silently distorting R².

  - create_grid_polygons() no longer truncates the grid when cellsize
    and n are both supplied. This changes results. sf::st_make_grid()
    does not ignore n when cellsize is given: for square grids it takes
    the cell dimensions from cellsize and the counts from nx = n[1], ny
    = n[2], anchored at the bounding-box corner. cellsize = 25 with n
    = 2 on a 100 × 100 boundary therefore produced 4 cells
    covering 2,500 of 10,000 square units and silently left three
    quarters of the study area with no cells at all — and because clip =
    TRUE had nothing outside the boundary to discard, the result looked
    like an ordinary, complete grid. cellsize now wins, n is dropped
    with a logged warning naming what it would have done, and the same
    call returns 16 cells covering the whole boundary. n is still
    forwarded when the package derived cellsize from it or from
    target_cells, which is what the original code was written for:
    omitting it there lets sf recompute ceiling(w / cellsize) and
    floating-point division pushes the count one past the intended
    value.

  - fit_gwr_model() no longer refuses a continuous response that happens
    to take two values. This changes results: fits that used to error
    now run. The guard rejected any response with exactly two distinct
    finite values as "binary" and pointed at GWmodel::ggwr.basic(family
    = "binomial"). Two distinct values is not the same thing as binary:
    a measurement censored at a detection limit or saturated at a
    ceiling (0.0031 / 12.7401) is perfectly continuous, Gaussian GWR on
    it is a well-defined least-squares problem, and the advice to switch
    to a binomial family is nonsense for such values. The hard stop is
    now gated on the response also being integer-like, which is what the
    surrounding code already used to separate coded categories from
    measurements. A two-valued non-integer response raises a warning()
    naming the two values and asking you to confirm it is genuinely
    continuous, then fits. This also mattered inside cv_gwr(), where the
    guard runs once per fold and a small training fold can legitimately
    hold only two distinct values.

  - determine_optimal_levels() no longer reports a Moran's I that is
    arithmetically fixed. This changes which cell counts it returns.
    .morans_i_for_k() builds a min(8, n_cells - 1)-nearest-neighbour
    weight matrix, so at nine cells or fewer every cell neighbours every
    other one. The row-standardised matrix is then complete, W %*% e =
    -e/(n - 1) for any mean-zero residual vector, and Moran's I
    collapses to exactly -1/(n_cells - 1) whatever the data are. That is
    not merely uninformative: |I| = 1/(n_cells - 1) falls monotonically
    in the number of cells, so criterion = "morans_i" ranked the largest
    evaluated candidate first every time, and "combined" carried the
    same tilt at half weight. Candidates below the floor now return
    NA_real_ and are excluded from the model-aware ranking; when none
    clears it — the usual outcome at the default max_levels = 12, since
    the search evaluates a window around the elbow — the call falls back
    to the geometric ranking and logs a warning. Raise max_levels above
    roughly 10 for the model-aware criteria to contribute at all.
    predictor_vars also accepts logical columns now, read as 0/1,
    matching fit_rf_model()/cv_rf()/predict(); factor and character
    predictors are still refused by name.

  - residual_morans_i(fit, k = 1) no longer errors with "subscript out
    of bounds" on a machine without FNN. In the dense fallback the inner
    function returns a scalar at k = 1, so apply() simplified the
    neighbour table to a length-n vector and t() made it a 1 × n matrix;
    indexing nn_idx[i, ] then failed for every i > 1. The result is now
    forced to n × k.

  - make_folds() no longer dies on an empty or non-finite geometry.
    st_coordinates() yields one all-NA row per EMPTY POINT rather than
    zero rows, so a row-count check let them through: block_kfold's
    st_intersects() returned integer(0), ..block_id went NA, and the
    nearest-block rescue aborted with "replacement has length zero".
    Unusable rows are now dropped with a warning naming the count, after
    ..row_id is stamped so the survivors keep their original row
    identities, and for every method rather than just block_kfold —
    random_kfold would otherwise put an unplottable point in a fold, and
    nndm and buffered_loo both feed the coordinates to distance code.
    The rescue itself uses vapply() rather than apply(), so a point
    whose distances are all NA keeps its NA instead of collapsing the
    assignment. points_sf with no usable coordinates at all is an error
    naming that, not a downstream one.

  - Every cross-validation wrapper names the cause when folds fail.
    .cv_run_folds() returns each fold's error text rather than a bare
    NULL, and cv_gwr(), cv_bayes() and cv_spatial() append First error:
    ... to both the logged and the R-level "all N folds failed" message.
    Running cv_bayes() without brms installed previously produced five
    fold N fit failed warnings and an all-NA $overall with n_pred = 0 in
    which the word "brms" never appeared.

  - The package's log lines no longer depend on the user's global logger
    configuration. logger seeds a new namespace from the global one, so
    the "spatialkit" namespace inherited whatever formatter the user had
    set before loading — and every logging helper hands logger an
    already formatted string. Under a user's formatter_sprintf, every
    package message containing a literal % — the CRS distortion figures
    in ensure_projected(), the local collinearity percentage in
    fit_gwr_model() — hard-errored with "too few arguments", and because
    the helper logs before it raises the R warning, the warning the
    manual promises died with it. Under the default formatter_glue a
    {...} inside a fold error was re-evaluated. The namespace's
    formatter is now pinned to formatter_paste, so the message logged is
    the message written.

  - fit_bayesian_spatial_model(backend = "auto") chooses cmdstanr only
    when a CmdStan build is actually available. It chose it whenever the
    cmdstanr package could be loaded — a thin interface that is often
    installed without the toolchain it drives — so on such a machine
    every fit died inside the sampler with "CmdStan path has not been
    set yet. See ?set_cmdstan_path". The package's own weekly check-brms
    job was one such machine: it installs cmdstanr to satisfy Suggests
    and never builds CmdStan, and every scheduled run since the Stan
    smoke tests landed failed there. "auto" now falls back to rstan,
    which brms always brings, and logs the choice; an explicit backend =
    "cmdstanr" with no usable build is an error that says to run
    cmdstanr::install_cmdstan().

  - cv_bayes(seed = ) reaches the sampler. fit_bayesian_spatial_model()
    carries seed = 123 and the per-fold fit_args never set it, so every
    fold of every run sampled from Stan seed 123 and changing seed
    changed nothing on fixed folds. Each fold now draws its own sampler
    seed from the fold's seeded stream, as cv_rf() does for the forest;
    a seed in fit_args still overrides it for every fold.

  - cv_*(seed = NULL) is reproducible from set.seed() under parallel
    > 1, as the README promised without qualification. With seed = NULL
    no per-fold seeds were drawn, so each forked worker was seeded by
    mclapply() from the clock and the process ID — three runs after the
    same set.seed(777) gave 0.5687, 0.5649 and 0.5623 while the
    sequential call was reproducible. The per-fold seeds are now drawn
    from the caller's current stream (advancing it, as any RNG-consuming
    call would), so the sequential and parallel paths are the same
    function of the state set.seed() left.

  - Warnings raised inside a fold reach the caller from the parallel
    path. R conditions do not cross a fork, so under parallel > 1 every
    warning raised by the model — including fit_gwr_model()'s documented
    integer-response warning, raised in every fold — reached nobody,
    while the numbers came back identical and the run looked like a
    clean version of the same analysis. The worker now collects them and
    the parent re-raises each distinct message once.

  - cv_spatial() (and therefore cv_rf()) returns the same typed,
    zero-row fold_metrics frame as cv_gwr() and cv_bayes() when every
    fold failed, so subset(fold_metrics, RMSE < 5) works instead of
    erroring on a missing column.

  - residual_morans_i() refuses k large enough to make the neighbour
    matrix dense. The only size guard (n > 5000) applied to the dense
    fallback; with FNN and Matrix present, k >= n - 1 allocated n (n
    − 1) pairs unguarded. Requests above 2e7 pairs are now an error
    naming k and n.

New features

  - New fit_rf_model() and cv_rf(): a ranger random forest as a
    first-class backend, returning an rf_fit that works with
    cv_spatial(), predict_surface(), area_of_applicability() and plot()
    like any other model. Three defaults are opinionated: include_coords
    = FALSE (a forest given the coordinates memorises location and fails
    wherever it has not been — Meyer et al. 2019,
    https://doi.org/10.1016/j.ecolmodel.2019.108815 — and random CV does
    not catch it); fitted() returns out-of-bag predictions, so summary()
    on an rf_fit is not comparable with the other backends and says so
    ($info$fitted_are_oob); and importance defaults to permutation
    rather than impurity, which is biased toward continuous and
    high-cardinality predictors (Strobl et al. 2007,
    https://doi.org/10.1186/1471-2105-8-25). Compare backends with
    compare_models_cv(), which now has an RF branch.
    
    predict() on an rf_fit refuses the type confusions ranger would
    otherwise absorb silently: a numeric-at-fit predictor supplied as
    text (which ranger factor-codes, then applies numeric split
    thresholds to the codes), a logical-at-fit predictor supplied as
    text, and a categorical level the forest was never grown with — the
    level set is not enough, since a spatial fold holding out a whole
    class leaves a level with no training rows. Arguments that make
    ranger return a matrix (predict.all = TRUE, type = "quantiles") are
    rejected rather than flattened column-major. A constant seed is
    supplied to ranger's predict unless the caller passes one, so
    prediction does not consume the global RNG and
    predict_surface(chunk_size = ) — a performance knob — cannot shift
    later random draws. cv_rf(seed = ) reaches the forest in every fold,
    and gains pointize. Passing ranger's own spelling of an argument the
    wrapper already sets (num.trees, min.node.size, num.threads, mtry,
    importance, seed) through ... is an error naming the wrapper
    argument to use, rather than reaching ranger() twice. num_threads
    defaults to getOption("mc.cores", 1L) — one thread unless the
    session has opted in — for both the fit and predict(), rather than
    ranger's own default of every core on the machine, and
    cv_rf(parallel = ) runs each forked fold's forest on one thread
    unless told otherwise, so the worker count is never multiplied by a
    thread count. See ?fit_rf_model.

  - Every cv_*() and compare_models_cv() accept folds as a vector of
    fold labels, one per row — make_folds()$assignment$fold, the object
    most naturally to hand — in addition to a make_folds() result and a
    list of train/test splits, the three shapes area_of_applicability()
    already took. The label vector used to fail with R's "$ operator is
    invalid for atomic vectors".

  - New area_of_applicability(), implementing the dissimilarity index of
    Meyer & Pebesma (2021, https://doi.org/10.1111/2041-210X.13650).
    Predictors are centred and scaled on the training data's own
    statistics, optionally weighted by variable importance — by the
    importance itself, not its square root, matching CAST. A prediction
    point's DI is its distance to the nearest training point in that
    space over the mean pairwise training distance, and the threshold is
    the outlier-removed maximum of the training data's own DI. Pass the
    make_folds() result you actually validated with — the area is
    defined relative to a performance estimate, and a blocked estimate
    is a claim about predicting further away.
    
    A model fitted with include_coords = TRUE is measured in coordinate
    space, since an index that ignores location would report a point far
    outside the training extent as inside on ordinary covariate values
    alone; weights for the two coordinate columns default to the mean of
    those supplied, as the caller has never seen them. Non-POINT newdata
    is reduced to points, and a CRS present on one side is applied to
    the other. The zero-variance test is relative to each column's
    magnitude rather than an absolute tolerance, so a predictor is not
    dropped for the unit it was recorded in. Categorical predictors are
    refused rather than dummy-coded; logicals are read as 0/1. A
    make_folds() result is resolved by its ..row_id values, which
    coincide with row positions only when the input carried no prior
    IDs. See ?area_of_applicability.

  - New select_features_forward(): greedy forward feature selection with
    spatially blocked inner folds, which is the whole point of having
    it. Random inner folds inside blocked outer folds select variables
    that look predictive only because nearby points leak between train
    and test, and the outer loop then reports honest-looking numbers for
    a dishonestly chosen feature set. method defaults to "block_kfold"
    and logs a warning if set to "random_kfold". The empty set is scored
    first where the backend can fit it, so the first variable is judged
    against a null-model baseline rather than accepted unconditionally,
    and history carries that baseline as a step = 0 row. Every candidate
    set is scored on the same observations — the completeness filter
    matches prep_model_data() exactly, finiteness test included, so a
    candidate carrying a single Inf cannot be preferred for having an
    easier subset — and the inner folds are built once, before the
    sweep, rather than rebuilt per candidate. A max_fits budget guards
    against nesting a sweep inside leave-one-out outer folds. Where the
    backend cannot fit the empty set at all — fit_rf_model() and
    fit_gwr_model() both refuse a zero-length predictor_vars — the probe
    is silent on the console: its per-fold failures go to the file trace
    only, rather than printing the same lines a genuinely failed run
    prints.

  - New gwr_model_selection(): wraps GWmodel::gwr.model.selection() (Lu
    et al. 2014, https://doi.org/10.1080/10095020.2014.917453) and
    returns a ranked table instead of two loosely-coupled lists. It is
    the fast, in-sample counterpart to select_features_forward() — the
    same forward search scored by AICc, read from the documented
    c(bandwidth, AIC, AICc, RSS) layout of GWmodel's GWR.df, which
    carries no column names; the result records whether the table
    arrived in that shape. Candidates must be numeric, and dmat_max_n =
    Inf means always precompute the distance matrix. Both limitations
    are documented rather than papered over: one bandwidth is shared by
    every candidate (which is what makes the criteria comparable), and
    the null model is never evaluated, so the result always names at
    least one predictor. When it disagrees with the blocked estimate,
    believe the blocked estimate. See ?gwr_model_selection.

  - New predict_surface(): builds a regular grid over the training
    extent (or a grid you supply), joins covariates, predicts in chunks
    and returns sf. Supports boundary clipping, cell_size or approximate
    n_cells, and se = TRUE for a posterior-SD surface where the backend
    exposes draws.

  - New plot() method for spatial_fit, with type = "residuals",
    "observed_predicted" and "variogram" (the empirical residual
    variogram with the fitted model and effective range overlaid, so the
    fit can be judged rather than trusted). The variogram's distance
    axis is labelled in the units of the CRS it was actually fitted in —
    metres of an auto-chosen zone for a lon/lat fit, not the caller's
    degrees — it names the azimuth when a single direction is drawn, and
    a fit that did not converge says so in the caption. New plot_folds()
    maps a fold scheme, which is the fastest way to see whether spatial
    blocks separate the data or are smaller than the autocorrelation
    range and therefore leaking.

  - make_folds() gains method = "leave_location_out", which keeps every
    observation from a location (named by the new group_var) in the same
    fold. Repeated measurements at one site were previously
    unrepresentable: random k-fold splits them across folds, so the
    model is scored partly on sites it trained on.

  - make_folds() gains method = "nndm", implementing the
    distance-matching principle of Milà et al. (2022,
    https://doi.org/10.1111/2041-210X.13851), as in CAST::nndm(). Rather
    than choosing a buffer with nothing to justify it, the exclusion
    around each held-out point is sized so the training-to-test distance
    distribution reproduces the distances from your actual prediction
    locations (the new prediction_points) to the training data. The
    procedure follows the paper's iterative exclusion removal for
    removal and is deterministic: no random numbers are drawn, so the
    caller's RNG is untouched, and ties in the nearest-neighbour
    distance — every mutual-nearest-neighbour pair, all of a regular
    grid — are broken by the point's position rather than by its row
    index, so identical data give identical folds whatever order the
    rows arrive in (CAST breaks them by row). params$target_median,
    params$realised_median and params$max_ecdf_excess record how close
    the match came, and min_train (default 0.5) and phi control it.
    Matching is as close as the training configuration permits — the
    achievable distances are discrete order statistics. When prediction
    locations sit no further from the training data than training points
    sit from each other, plain leave-one-out already reproduces the
    target and nothing is excluded; that is the correct outcome. A
    non-POINT prediction_points layer is reduced to points first, since
    point-to-polygon distances are zero for any cell containing a
    training point and would collapse the scheme towards plain LOO.

  - summarize_by_cell() gains deff = "variogram", computing a per-cell
    design effect from a fitted variogram rather than one pooled
    intra-class correlation. For n points in a cell with correlation
    matrix R the effective sample size of the mean is n^2 / sum(R), so
    deff = sum(R) / n. This generalises the Kish option — a constant
    off-diagonal correlation recovers 1 + (n - 1) * rho exactly — but
    lets correlation decay with distance, which is what having fitted a
    variogram is for. Pass the fit via the new sac argument, or it is
    estimated when response_var is supplied — on the response, not on
    OLS residuals, even when predictor_vars are listed: the ..se_resp_*
    columns estimate the SE of the cell mean as an estimate of the
    response's grand mean, so the correlation to correct for is the
    response's own (measured grand-mean coverage 0.93 with the response
    variogram against 0.51 with the residual one on a field with a
    smooth predictor). Pass a residual variogram through sac if that is
    the field you want. Large cells are subsampled at deff_max_n
    (default 500), with the correlation scaled back to the cell's own
    size. A sac_range whose fit was rejected carries no usable
    correlation function, so both the supplied and the internally
    estimated path fall back to deff = 1 and say so rather than
    saturating the correlation at every within-cell distance. One
    correlation function is fitted and applied to every numeric column,
    response and predictors alike, because a variogram is a property of
    the field rather than of a variable type.

  - fit_bayesian_spatial_model() supports intercept-only models
    (predictor_vars = character(0)): the response is explained by the
    intercept and the spatial GP alone, the natural null for asking how
    much of a surface is spatial structure rather than covariate effect.

  - fit_bayesian_spatial_model() checks the posterior length-scale
    against the smallest scale the chosen basis can resolve and logs a
    warning when more than 10% of the posterior mass falls below it —
    the adequacy diagnostic recommended by Riutort-Mayol et al. (2023,
    https://doi.org/10.1007/s11222-022-10167-2), and what makes the
    smaller default gp_k safe rather than merely cheaper. $info gains
    gp_c, gp_n_basis, gp_ell_min and gp_lengthscale_bounds, and print()
    on a bayesian_fit and cv_bayes()'s fold_metrics report the total
    basis count alongside the per-dimension rank.

  - cv_spatial() raises a condition when folds fail, matching cv_gwr()
    and cv_bayes(); an all-failing fit_fn previously returned an all-NA
    overall and an empty fold_metrics with nothing at R condition level.
    The result records n_folds_attempted and n_folds_succeeded — compare
    them before trusting overall.

  - make_folds() records the CRS the folds were built in as params$crs
    ("EPSG:32632", an input string, or a WKT). block_size and sac_range
    are lengths in that CRS, which is not necessarily the one the caller
    passed: geographic input is projected by ensure_projected() to a CRS
    chosen for the extent. Without the label the units of a recorded
    block size were not recoverable from the result.

  - spatialkit_quiet() is a new exported helper. Both
    logger::log_appender() and logger::log_threshold() default to index
    = 1, which is the temp-file trace, so the two-line recipe in the
    README could not redirect or quieten the console echo (index 2) —
    there was no documented way to silence the package. The README now
    says so too.

Documentation

  - estimate_sac_range() documents its three return shapes (a range, a
    rejected range, and no fit at all) and which attributes each
    carries.

  - make_folds() documents that k is not always honoured: buffered_loo
    and nndm always return k = n, and block_kfold and leave_location_out
    lower it when the geometry or the grouping cannot support the
    request. Read folds$k.

  - new_spatial_fit() documents the coef() contract;
    summary.spatial_fit() and model_metrics() document that their
    metrics are in-sample for a gwr_fit and a bayesian_fit but
    out-of-bag for an rf_fit; prep_model_data() documents that the
    projected CRS is not an unconditional guarantee, since
    ensure_projected() passes a CRS-less dataset through unchanged when
    its coordinates do not look like lon/lat.

  - The vignette and inst/scripts/example_nc_demo.R read fit quality
    from fit$metrics$r_squared and CV results from cv$summary$rmse.
    Neither field has ever existed. Because sprintf() returns
    character(0) when any argument has length zero, the reporting lines
    printed nothing rather than erroring, so the shipped vignette
    silently omitted every number it claimed to show. Both now use
    model_metrics() and $overall.

  - The demo's Voronoi tessellation was built from all 300 observations
    rather than from the 40 k-means seeds it computed one line earlier —
    one cell per observation, a nearest-neighbour interpolation rather
    than an aggregation, compared side by side against two ~50-cell
    grids. The seeds are now used.

  - The vignette builds as rmarkdown::html_vignette rather than
    html_document, guards its ggplot2 and geometry use, demonstrates
    summarize_by_cell() instead of reimplementing it with
    group_by()/summarise(), and adds a spatial cross-validation section
    contrasting block_kfold against random_kfold on the same data.

  - The package-level help page (?spatialkit) gains "The pipeline, in
    order" and "Where to start" sections, so help(package =
    "spatialkit") leads somewhere rather than presenting 40 exports in
    alphabetical order.

  - Every exported function's description now says when to reach for it
    rather than only what it does, and @family / @seealso links connect
    each step of the pipeline to the one before and after it —
    assign_features_to_polygons() to summarize_by_cell(),
    determine_optimal_levels() to build_tessellation(),
    new_spatial_fit() to cv_spatial(), and the two seeding functions to
    each other. create_voronoi_polygons() versus create_grid_polygons(),
    and voronoi_seeds_kmeans() versus voronoi_seeds_random(), each say
    which to pick and why.

  - build_tessellation() documents that boundary is required for method
    = "hex" and method = "square" — the grid methods have no extent of
    their own — and optional for "voronoi" and "triangles", which derive
    one from the points. The error existed; the requirement was not
    written down anywhere.

  - create_grid_polygons() documents that target_cells, cellsize and n
    are three ways of sizing one grid and that exactly one should be
    supplied, that cellsize is in the units of the working CRS, and that
    cellsize takes precedence over n.

  - determine_optimal_levels() documents the nine-cell resolution floor
    on the model-aware criteria, why it exists, and that the whole call
    falls back to the geometric ranking when no candidate clears it.

  - compare_models_cv() documents that dropping every requested backend
    is an error ("no viable models.") rather than an empty comparison,
    and that the returned frame carries only the models that actually
    ran, so callers should check which names are present rather than
    assuming one row per request.

  - new_spatial_fit() documents the two obligations on a custom backend:
    return an object built by the constructor, and define a
    predict.<subclass>() method — cv_spatial() scores folds through the
    predict() generic, so without one every fold fails.

  - README. A new "Your own data" section shows both entry points —
    st_read() for a spatial file and read.csv() + st_as_sf() for a table
    of coordinates — using the nc.shp demo shapefile shipped with sf so
    it runs anywhere. The README previously manufactured every example
    inline with a hard-coded crs = 32632 and never showed data entering
    the package at all. A companion "CRS: what the numbers are in"
    subsection states that block sizes, buffers, bandwidths, variogram
    ranges and expand distances are in the units of the working CRS;
    that geographic input is projected automatically to a CRS chosen for
    the extent; and how to pin one.

  - README. New guidance where none existed: how to choose among the
    four tessellation methods, how k and block_size trade off against
    the autocorrelation range, what to do when estimate_sac_range()
    returns NA, how to read a design effect, which model backend to
    reach for (with the recorded cost of each), and a "Troubleshooting"
    section covering the errors a new user actually hits first. A worked
    hex-grid example replaces the previous picture-only coverage of the
    grid methods.

  - README. Three corrections. The estimate_sac_range() example showed a
    rejected range printing its attributes, which print.sac_range() has
    not done since the attribute dump was removed; it now shows the bare
    NA and reads the attributes explicitly. The
    determine_optimal_levels() passage claimed the
    residual-autocorrelation criterion was doing work at cell counts
    where it is arithmetically degenerate. The test-suite paragraph said
    "exactly one" test guards on brms; six do, five of them additionally
    gated behind SPATIALKIT_TEST_BRMS so they never run in the matrix.

  - inst/scripts/example_nc_demo.R said EPSG:2264 was projected "so
    distances are metric". Its unit is the US survey foot, which is what
    the script's own "Autocorrelation range: %.0f ft" line reports. The
    comment now says planar, and names the unit every distance,
    bandwidth and block size in the script is in.

  - Vignette. print(rf_fit) and summary(rf_fit) report the same OOB RMSE
    but different R² (0.4733 against 0.4715). The vignette now explains
    why: print.rf_fit() echoes ranger's r.squared (1 - MSE/var(y),
    unbiased n − 1 variance) while summary() recomputes 1 -
    SS_res/SS_tot from the same out-of-bag predictions with an n
    denominator, so the unexplained fractions differ by exactly n/(n
    − 1).

  - Every exported function has runnable examples: the eleven that
    shipped without any — clear_fitted_cache(), clear_grid_cache(),
    clip_target_for(), compare_models(), create_grid_polygons_cached(),
    ensure_stable_poly_id(), evaluate_insample(), harmonize_crs(),
    model_metrics(), voronoi_seeds_kmeans() and voronoi_seeds_random() —
    gained one, and the two \dontrun{} blocks say why they cannot be run
    (a Stan toolchain and minutes of MCMC).

  - residual_morans_i()'s default null is described correctly: null =
    "auto" uses the Cliff & Ord regression-residual moments whenever the
    residuals are OLS residuals on the rebuilt design, and the
    randomisation null otherwise; the README said "the randomisation
    variance" without qualification. The type-I error of a random
    forest's residual test is attributed to what the package actually
    feeds it — out-of-bag residuals, which are honest out-of-sample
    errors with their own spatial structure — not to "shrunk in-sample
    residuals". determine_optimal_levels() gives the real reason for its
    nine-cell floor (the standardised deviate is 0/0 there, so the
    criterion carries no information) rather than an argument from |I|
    that its own details section had just called wrong.
    summarize_by_cell() notes that the "use the naive SE for the cell's
    own mean" advice is calibrated under uniform within-cell sampling.
    cv_gwr(bandwidth = ) states its units and semantics like its
    siblings. quiet is documented as "suppress this function's progress
    messages" everywhere, with a pointer to spatialkit_quiet() for the
    console log echo it does not touch.

  - README: the square-grid call returns 36 cells, not 32; the
    installation section no longer promises a specific version from
    CRAN; the resolution figure and the quick-start output are
    regenerated for the elbow-first ordering (4 3 5, k = 4).

  - The DESCRIPTION now cites the methods it implements -- Lu et al.
    (2014) for GWmodel, Riutort-Mayol et al. (2023) for the Hilbert
    space Gaussian process, Strobl et al. (2007) for the permutation
    importance, Mila et al. (2022) for NNDM folds and Meyer and Pebesma
    (2021) for the area of applicability -- each with its DOI, and
    quotes only software names. The Riutort-Mayol reference on
    fit_bayesian_spatial_model()'s help page gives the article number
    (33, 17) rather than "33, 1".

  - No example is wrapped in \donttest{} any more: the fifteen that were
    -- every fit, cross-validation, comparison, plotting and surface
    example that needs a Suggests package -- run unconditionally behind
    their requireNamespace() guards, the slowest in under 2 s. The two
    Stan examples (fit_bayesian_spatial_model(), cv_bayes()) keep
    \dontrun{} because they need a C++ toolchain and minutes of MCMC,
    and say so in a leading comment.

                 Changes in version 1.0.0 (2026-08-07)                  

First CRAN release, published 2026-08-07.

  - CRS management: ensure_projected(), harmonize_crs(),
    coerce_to_points(), prep_model_data().
  - Voronoi, hexagonal, square and Delaunay tessellation
    (build_tessellation() and the create_*_polygons() functions), with
    boundary clipping, stable reproducible cell IDs
    (ensure_stable_poly_id()) and a memoised grid builder
    (create_grid_polygons_cached()).
  - Seeding (get_voronoi_seeds(), voronoi_seeds_kmeans(),
    voronoi_seeds_random()) and resolution selection
    (determine_optimal_levels()).
  - Feature-to-polygon assignment (assign_features_to_polygons()) and
    cell-level aggregation with design-effect-corrected standard errors
    (summarize_by_cell()).
  - GWR (fit_gwr_model()) and Bayesian spatial Gaussian process
    (fit_bayesian_spatial_model()) backends behind a common spatial_fit
    S3 class.
  - Spatial cross-validation: make_folds() with random_kfold,
    block_kfold and buffered_loo; estimate_sac_range(); cv_gwr(),
    cv_bayes() and cv_spatial().
  - Model comparison and diagnostics: compare_models(),
    compare_models_cv(), evaluate_insample(), residual_morans_i().
  - Tessellation mapping (plot_tessellation_map()) and scoped logging.
