diff --git a/DESCRIPTION b/DESCRIPTION
index 1bd0629..f6b6a93 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,6 +22,7 @@ Depends:
 Imports:
     datamatch,
     dplyr,
+    fancyfx,
     ggplot2,
     parallel,
     parsnip,
@@ -38,7 +39,6 @@ Imports:
     yardstick
 Suggests:
     derivoce,
-    fancyfx,
     knitr,
     leaflet,
     ranger,
diff --git a/NAMESPACE b/NAMESPACE
index c903863..7637f4c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -21,6 +21,7 @@ export(default_copernicus_datasets)
 export(default_species_catalog)
 export(derivoce_choices)
 export(derivoce_covariates)
+export(derivoce_dependency_columns)
 export(derivoce_required_inputs)
 export(derivoce_steps_for)
 export(ensemble_rules)
@@ -69,11 +70,13 @@ export(power_curve)
 export(prejoin_steps)
 export(project_patch_model)
 export(projection_map)
+export(projection_overlap)
 export(raw_abundance_suffix)
 export(run_taupatch)
 export(run_taupatch_app)
 export(save_config)
 export(single_stages)
+export(spatial_bias)
 export(species_catalog_from)
 export(split_dates)
 export(stage_suffix_pattern)
diff --git a/R/derivoce.R b/R/derivoce.R
index 8fea5d7..d65e8b5 100644
--- a/R/derivoce.R
+++ b/R/derivoce.R
@@ -542,9 +542,10 @@ derivoce_choices <- function(selected, bathymetry = character(),
   available <- union(selected, as.character(fetchable))
 
   candidate <- function(id, label, group, step, expensive = FALSE,
-                        requires = character()) {
+                        requires = character(), depends = character()) {
     list(id = id, label = label, group = group, expensive = expensive,
-         step = step, requires = setdiff(requires, selected))
+         step = step, requires = setdiff(requires, selected),
+         depends = depends)
   }
 
   # One per fetched covariate. These are the cheap ones, and the ones a habitat
@@ -575,6 +576,16 @@ derivoce_choices <- function(selected, bathymetry = character(),
          step = function(v) list(type = "distance_to_front", var = v))
   )
 
+  # Selected covariates first, then the rest of the catalogue. Deriving from a
+  # covariate you are not modelling is an ordinary thing to want - integrated
+  # chlorophyll without chlorophyll itself, because the accumulated bloom is
+  # what feeds the animals and the instantaneous value is not - and there was
+  # no way to ask for it: the derived covariate was only offered once its
+  # source had been made a predictor, and then there was no way to take the
+  # source back out.
+  #
+  # The extra ones say so in their group rather than their label, so the cost
+  # is visible in the picker without every label carrying a parenthesis.
   out <- list()
   for (entry in per_covariate) {
     for (v in selected) {
@@ -583,6 +594,13 @@ derivoce_choices <- function(selected, bathymetry = character(),
         entry$step(v), entry$expensive
       )
     }
+    for (v in setdiff(available, selected)) {
+      out[[length(out) + 1]] <- candidate(
+        paste0(v, entry$suffix), sprintf(entry$label, v),
+        paste0(entry$group, " (downloads ", v, ")"),
+        entry$step(v), entry$expensive, requires = v
+      )
+    }
   }
 
   # Steps that read particular covariates, offered only when those were fetched.
@@ -602,6 +620,28 @@ derivoce_choices <- function(selected, bathymetry = character(),
       "EKE", "Eddy kinetic energy", "Flow",
       list(type = "eke"), requires = c("UO", "VO")
     )
+    # Steps see what earlier steps produced, so a gradient can be taken of the
+    # speed rather than of the two components it came from. That is the
+    # quantity a front in the flow actually is - the components can each be
+    # changing steeply while the speed is constant, which is a turn and not a
+    # shear - and it is the original pipeline's uv_grad.
+    #
+    # `depends` rather than `requires` because what is needed is another step,
+    # not another download: picking this pulls current_speed in whether or not
+    # the speed itself was asked for as a predictor.
+    for (entry in per_covariate) {
+      if (entry$type == "distance_to_front") next
+      out[[length(out) + 1]] <- candidate(
+        paste0("speed", entry$suffix), sprintf(entry$label, "current speed"),
+        "Flow", entry$step("speed"), entry$expensive,
+        requires = c("UO", "VO"), depends = "speed"
+      )
+    }
+    out[[length(out) + 1]] <- candidate(
+      "EKE_grad", "Spatial gradient of eddy kinetic energy (per km)", "Flow",
+      list(type = "horizontal_gradient", vars = "EKE"),
+      requires = c("UO", "VO"), depends = "EKE"
+    )
     # Backward rather than forward, because backward finds the attracting
     # structures where water converges and plankton accumulate, which is the
     # question a habitat model is asking. Forward finds transport barriers, and
@@ -647,10 +687,40 @@ derivoce_steps_for <- function(ids, selected, bathymetry = character()) {
   if (length(ids) == 0) return(list())
 
   choices <- derivoce_choices(selected, bathymetry)
-  chosen <- Filter(function(x) x$id %in% ids, choices)
+  chosen <- Filter(function(x) x$id %in% with_dependencies(ids, choices),
+                   choices)
   lapply(chosen, function(x) x$step)
 }
 
+#' The chosen derived covariates, plus the ones they are computed from
+#'
+#' A gradient of current speed needs the speed, and the speed is itself a
+#' derived covariate rather than a download. Choosing the gradient therefore
+#' has to pull in the step that produces what it reads — otherwise the config
+#' asks derivoce for a gradient of a column that was never computed.
+#'
+#' Order is not this function's problem. `derivoce_choices()` lists a
+#' dependency before anything that depends on it, and steps are emitted in that
+#' order, so `current_speed` runs before the gradient of `speed` without
+#' anything having to sort them.
+#'
+#' @param ids chosen derived covariate ids
+#' @param choices the [derivoce_choices()] list
+#' @return `ids` with any dependencies added
+#' @keywords internal
+with_dependencies <- function(ids, choices) {
+  needed <- ids
+  repeat {
+    depends <- unlist(lapply(
+      Filter(function(x) x$id %in% needed, choices),
+      function(x) x$depends %||% character()
+    ))
+    grown <- union(needed, depends %||% character())
+    if (length(grown) == length(needed)) return(grown)
+    needed <- grown
+  }
+}
+
 #' Covariates a set of derived choices needs fetching
 #'
 #' A derived covariate is computed from others, and those have to be downloaded
@@ -670,6 +740,47 @@ derivoce_required_inputs <- function(ids, selected, bathymetry = character()) {
   if (length(ids) == 0) return(character())
 
   choices <- derivoce_choices(selected, bathymetry)
-  chosen <- Filter(function(x) x$id %in% ids, choices)
+  chosen <- Filter(function(x) x$id %in% with_dependencies(ids, choices),
+                   choices)
   unique(unlist(lapply(chosen, function(x) x$requires))) %||% character()
 }
+
+#' Derived columns pulled in only as ingredients
+#'
+#' Choosing the gradient of current speed computes the speed on the way, and
+#' that column then sits in the modelling data looking exactly like one that was
+#' asked for. It was not: the ingredient of a derived covariate is no more a
+#' predictor than the velocity components behind an FSLE are.
+#'
+#' This names the derived columns a selection produced without anyone choosing
+#' them, so they can go into `covariates.exclude` alongside the downloads that
+#' [derivoce_required_inputs()] finds. A column named here is still computed and
+#' still available to anything later that reads it; it just does not become a
+#' predictor.
+#'
+#' @param ids chosen derived covariate ids
+#' @param selected time-varying covariate names
+#' @param bathymetry static seafloor covariate names
+#' @return character vector of derived column names, possibly empty
+#' @examples
+#' # Asking for the gradient of current speed computes the speed too, and that
+#' # is an ingredient rather than a request.
+#' derivoce_dependency_columns("speed_grad", c("SST", "SSS"))
+#'
+#' # Asking for both makes the speed a request, so it is not excluded.
+#' derivoce_dependency_columns(c("speed", "speed_grad"), c("SST", "SSS"))
+#' @seealso [derivoce_required_inputs()], which does the same for downloads
+#' @export
+derivoce_dependency_columns <- function(ids, selected,
+                                        bathymetry = character()) {
+  if (length(ids) == 0) return(character())
+
+  choices <- derivoce_choices(selected, bathymetry)
+  pulled <- setdiff(with_dependencies(ids, choices), ids)
+  if (length(pulled) == 0) return(character())
+
+  # The id of a per-covariate candidate is the column it produces, which is
+  # what `exclude` has to name.
+  chosen <- Filter(function(x) x$id %in% pulled, choices)
+  unique(vapply(chosen, function(x) x$id, character(1)))
+}
diff --git a/R/model.R b/R/model.R
index 49f1860..bc98227 100644
--- a/R/model.R
+++ b/R/model.R
@@ -35,6 +35,10 @@ fit_patch_model <- function(dat, config) {
   }
 
   model_data <- dat[c(predictors, "patch")] |> as.data.frame()
+  # Kept beside the modelling data rather than in it: the coordinates are not
+  # predictors and must not become any, but spatial_bias() needs to know where
+  # each station was, and `.row` in the held-out predictions indexes this.
+  coordinates <- station_coordinates(dat)
 
   rec <- build_recipe(model_data, config)
   spec <- build_model_spec(config)
@@ -92,6 +96,12 @@ fit_patch_model <- function(dat, config) {
                                  times = bootstrap_times(config),
                                  seed = config$model$seed)
 
+  # Computed before the evaluation table so it can go in it. It is a property
+  # of the folds rather than of the fit, so every model type on these folds
+  # gets the same number - which is what makes it a caveat on a comparison
+  # rather than a score in one.
+  bias <- spatial_bias(list(coordinates = coordinates, predictions = predictions))
+
   # The point estimate stays the model fitted on everything. The ensemble only
   # ever adds columns beside it, so turning uncertainty on never moves a map.
   ensemble <- if (is.null(uncertainty)) {
@@ -103,7 +113,9 @@ fit_patch_model <- function(dat, config) {
   list(
     workflow = fitted,
     metrics = cv_metrics,
-    evaluation = evaluation_table(predictions, cv_metrics, cutoff, bounds),
+    evaluation = evaluation_table(predictions, cv_metrics, cutoff, bounds,
+                                  ssb = overall_ssb(bias)),
+    spatial_bias = bias,
     predictions = predictions,
     classification_threshold = cutoff,
     # The cutoff is estimated, and how much it moves decides whether a binarised
@@ -117,11 +129,25 @@ fit_patch_model <- function(dat, config) {
     # can be produced without refitting. It is the station table, so hundreds to
     # a few thousand rows.
     model_data = model_data,
+    coordinates = coordinates,
     predictors = predictors,
     threshold = attr(dat, "threshold")
   )
 }
 
+#' Station coordinates, in the order the modelling data is in
+#'
+#' `NULL` when the station table has no coordinates, which is the case for a
+#' hand-built frame in a test rather than for anything a run produces.
+#'
+#' @param dat labeled modeling data
+#' @return a two-column data frame of `lon` and `lat`, or `NULL`
+#' @keywords internal
+station_coordinates <- function(dat) {
+  if (!all(c("lon", "lat") %in% names(dat))) return(NULL)
+  as.data.frame(dat[c("lon", "lat")])
+}
+
 #' Assemble a self-explanatory evaluation table
 #'
 #' The cross-validated metrics table reports sensitivity, specificity and kappa
@@ -171,7 +197,8 @@ fit_patch_model <- function(dat, config) {
 #' \doi{10.1111/2041-210X.13140} — the same argument for rare events in species
 #' distribution models, which is what a patch is
 #' @keywords internal
-evaluation_table <- function(predictions, cv_metrics, cutoff, bounds = NULL) {
+evaluation_table <- function(predictions, cv_metrics, cutoff, bounds = NULL,
+                             ssb = NULL) {
   threshold_free <- c("roc_auc", "pr_auc")
 
   ranking <- data.frame(
@@ -195,7 +222,17 @@ evaluation_table <- function(predictions, cv_metrics, cutoff, bounds = NULL) {
   at_best$std_err <- NA_real_
   at_best$note <- "TSS-optimal cutoff; use this one to binarise a projection"
 
-  out <- rbind(ranking, at_default, at_best)
+  # Reported beside the ranking metrics rather than in a file of its own,
+  # because it is the number that says how much to believe them. A reader who
+  # sees roc_auc 0.86 and has to go looking for the caveat will not.
+  bias <- if (is.null(ssb) || is.na(ssb)) NULL else data.frame(
+    metric = "ssb", threshold = NA_real_, value = ssb, std_err = NA_real_,
+    note = paste("spatial sorting bias of the folds, not a model score;",
+                 "1 is fair, near 0 means the metrics above are optimistic"),
+    stringsAsFactors = FALSE
+  )
+
+  out <- rbind(ranking, at_default, at_best, bias)
 
   # Every row gets an interval, including the ones that never had a standard
   # error - which is the point. Column order puts the two uncertainty measures
diff --git a/R/pipeline.R b/R/pipeline.R
index 0241f8e..c9b384c 100644
--- a/R/pipeline.R
+++ b/R/pipeline.R
@@ -117,6 +117,10 @@ run_taupatch <- function(config_path, project = TRUE, keep_covariates = 50000) {
   # Read off the evaluation table rather than the metrics one, since that is
   # the table both a single model and an ensemble fill in the same way.
   message("  ROC AUC: ", signif(evaluation_value(model, "roc_auc"), 4))
+  # Said out loud rather than left in a column, because it is the caveat on the
+  # line above it and a reader who has to go looking for it will not.
+  ssb <- overall_ssb(model$spatial_bias)
+  if (!is.na(ssb)) message("  ", spatial_bias_note(ssb))
 
   projections <- NULL
   if (project) {
@@ -160,6 +164,9 @@ write_model_outputs <- function(model, config) {
   readr::write_csv(model$evaluation, file.path(out, "evals.csv"))
   readr::write_csv(model$metrics, file.path(out, "cv_metrics.csv"))
   readr::write_csv(model$importance, file.path(out, "var_importance.csv"))
+  if (!is.null(model$spatial_bias)) {
+    readr::write_csv(model$spatial_bias, file.path(out, "spatial_bias.csv"))
+  }
   plot_importance(model$importance, file.path(out, "var_importance.png"))
   write_diagnostic_plots(model, out)
 
diff --git a/R/plot_effects.R b/R/plot_effects.R
index 0997a89..e27ba3d 100644
--- a/R/plot_effects.R
+++ b/R/plot_effects.R
@@ -239,30 +239,6 @@ gam_smooth_terms <- function(fitted) {
   out[order(-out$edf), ]
 }
 
-#' Whether fancyfx is available to draw smooths
-#'
-#' Its own function so the optional path can be exercised in tests without
-#' mocking `requireNamespace()` itself, which every package that loads a
-#' graphics device also goes through.
-#'
-#' @section It used to be called fancygam:
-#' The package was renamed when it grew past GAMs. The rename is why this
-#' matters more than a find-and-replace: `chross22/fancygam` still resolves on
-#' GitHub, so `Remotes: chross22/fancygam` kept installing — but what it
-#' installs now declares `Package: fancyfx`, so `requireNamespace("fancygam")`
-#' returned `FALSE` on every fresh install and the smooth plots were skipped in
-#' silence. Anyone with the old package still sitting in their library saw
-#' nothing wrong.
-#'
-#' That is the failure mode to watch for here: this function gates a diagnostic
-#' rather than the run, so a wrong answer costs a plot and no error.
-#'
-#' @return `TRUE` when fancyfx is installed
-#' @keywords internal
-has_fancyfx <- function() {
-  requireNamespace("fancyfx", quietly = TRUE)
-}
-
 #' Variables a fitted GAM gave a smooth to
 #'
 #' Not every predictor gets one — [model_formula()] gives a linear term to any
@@ -287,8 +263,7 @@ gam_smoothed_variables <- function(fitted) {
 #' rather than reconstructed by prediction, so it carries uncertainty, which a
 #' partial dependence curve cannot.
 #'
-#' Drawn by [fancyfx](https://github.com/chross22/fancyfx), which is a Suggests
-#' — a run without it still gets the generic partial effect curves.
+#' Drawn by [fancyfx](https://github.com/chross22/fancyfx).
 #'
 #' @section Why the axes read in standard deviations:
 #' The smooths belong to the model, and the model was fitted on the recipe's
@@ -309,11 +284,6 @@ gam_smoothed_variables <- function(fitted) {
 #' @seealso [gam_smooth_terms()] for the numbers behind these
 #' @export
 plot_gam_smooths <- function(model, vars = NULL, path = NULL) {
-  if (!has_fancyfx()) {
-    stop("The 'fancyfx' package is required to plot GAM smooths. ",
-         "Install it with remotes::install_github('chross22/fancyfx').",
-         call. = FALSE)
-  }
   workflow <- if (inherits(model, "workflow")) model else model$workflow
 
   vars <- vars %||% gam_smoothed_variables(workflow)
@@ -382,16 +352,12 @@ write_effect_plots <- function(model, out) {
       written <- c(written, "smooth_terms.csv")
     }
 
-    # The fitted smooths themselves, with their uncertainty. Skipped without a
-    # word when fancyfx is absent: it is a Suggests, and the generic partial
-    # effect curves above already cover the question.
-    if (has_fancyfx()) {
-      smooths <- try_diagnostic(
-        plot_gam_smooths(model, path = file.path(out, "gam_smooths.png")),
-        "GAM smooth plots"
-      )
-      if (!is.null(smooths)) written <- c(written, "gam_smooths.png")
-    }
+    # The fitted smooths themselves, with their uncertainty.
+    smooths <- try_diagnostic(
+      plot_gam_smooths(model, path = file.path(out, "gam_smooths.png")),
+      "GAM smooth plots"
+    )
+    if (!is.null(smooths)) written <- c(written, "gam_smooths.png")
   }
 
   invisible(written)
diff --git a/R/plotting.R b/R/plotting.R
index 0e3cdf9..a239eac 100644
--- a/R/plotting.R
+++ b/R/plotting.R
@@ -23,23 +23,27 @@ plot_projection <- function(predicted, year, month, species, path) {
   invisible(path)
 }
 
-#' Plot how far a monthly projection can be trusted
+#' Draw the layers that say how far to trust a projection
 #'
-#' The companion to [plot_projection()]: the same cells, panelled by what is
-#' uncertain about them rather than by what is predicted. Whichever of the
-#' uncertainty surfaces the run produced are drawn, and nothing else.
+#' Two or three panels above the map, one per question the projection can be
+#' doubted on. Only the ones a run actually produced are drawn.
 #'
-#' The two panels are deliberately not merged into a single "confidence" layer.
+#' The panels are deliberately not merged into a single "confidence" layer.
 #' They measure different things and are free to disagree in either direction,
 #' and the disagreement is the informative part: a cell can be stable across
 #' every member and still be extrapolated, because agreement between members
-#' trained on the same data is not evidence about ground the data never covered.
-#' Blending the two would average that case away instead of showing it.
-#'
-#' Diverging colour on the novelty panel, centred at zero, because zero is the
-#' meaningful break: above it the model is interpolating, below it the cell is
-#' outside the training range on some predictor and the model has no evidence
-#' for what it says there.
+#' trained on the same data is not evidence about ground the data never
+#' covered. Blending them would average that case away instead of showing it.
+#'
+#' @section Drawn by fancyfx:
+#' The panels themselves come from [fancyfx::plotUncertainty()] and
+#' [fancyfx::plotExtrapolation()] rather than being built here. That is worth a
+#' sentence because it used to be ninety lines of `ggplot2` in this file, and
+#' the replacement is better in two ways this package would otherwise have had
+#' to write for itself: cells are downsampled above `max.cells`, which a real
+#' Copernicus grid needs and the hand-built version did not do, and the
+#' extrapolation ramp diverges about zero in colours that survive the common
+#' colour vision deficiencies.
 #'
 #' @param predicted a projection from `predict_grid()` with uncertainty columns
 #' @param year year being projected
@@ -50,83 +54,43 @@ plot_projection <- function(predicted, year, month, species, path) {
 #' @seealso [novelty_surface()] for how the novelty panel is computed
 #' @export
 plot_projection_uncertainty <- function(predicted, year, month, species, path) {
+  drawable <- intersect(c("suitability_sd", "algorithm_sd", "novelty"),
+                        names(predicted))
+  if (length(drawable) == 0) return(invisible(NULL))
+
   panels <- list()
 
-  if ("suitability_sd" %in% names(predicted)) {
-    panels$spread <- ggplot2::ggplot(
-      predicted, ggplot2::aes(x = .data$lon, y = .data$lat,
-                              fill = .data$suitability_sd)) +
-      ggplot2::geom_raster() +
-      ggplot2::scale_fill_viridis_c(option = "magma", na.value = "white",
-                                     name = "SD") +
-      ggplot2::labs(subtitle = "Spread across the ensemble")
+  if ("suitability_sd" %in% drawable) {
+    panels$spread <- fancyfx::plotUncertainty(
+      projection_raster(predicted, "suitability_sd"), legend.lab = "SD"
+    ) + ggplot2::labs(subtitle = "Spread across the ensemble")
   }
 
-  if ("algorithm_sd" %in% names(predicted)) {
-    # A third panel only when a multi-algorithm ensemble was fitted, and kept
-    # separate from the spread panel above it on purpose: that one is one
-    # algorithm refitted on resampled stations, this one is different
-    # algorithms on the same stations. A cell where the forest and the GLM
-    # disagree is not the same worry as a cell where the forest is unstable.
-    panels$algorithms <- ggplot2::ggplot(
-      predicted, ggplot2::aes(x = .data$lon, y = .data$lat,
-                              fill = .data$algorithm_sd)) +
-      ggplot2::geom_raster() +
-      ggplot2::scale_fill_viridis_c(option = "cividis", na.value = "white",
-                                     name = "SD") +
-      ggplot2::labs(subtitle = "Disagreement between algorithms")
+  # A multi-algorithm ensemble only. Kept separate from the panel above on
+  # purpose: that one is a single algorithm refitted on resampled stations,
+  # this one is different algorithms on the same stations, and a cell where the
+  # forest and the GLM disagree is not the same worry as a cell where the
+  # forest is unstable.
+  if ("algorithm_sd" %in% drawable) {
+    panels$algorithms <- fancyfx::plotUncertainty(
+      projection_raster(predicted, "algorithm_sd"), legend.lab = "SD",
+      option = "cividis"
+    ) + ggplot2::labs(subtitle = "Disagreement between algorithms")
   }
 
-  if ("novelty" %in% names(predicted)) {
-    # Extrapolated cells are usually a small minority, and a scale stretched
-    # over the whole range renders them invisible - which defeats the panel.
-    # So the visual weight goes to the negatives: saturated red below zero,
-    # against a muted ramp above it. The alarming colour marks the cells to
-    # distrust rather than, as a plain diverging scale does, the safest ones.
-    outside <- sum(!is.na(predicted$novelty) & predicted$novelty < 0)
-    scored <- sum(!is.na(predicted$novelty))
-
-    subtitle <- if (outside == 0) {
-      "Every cell is inside the training range"
-    } else {
-      culprit <- names(sort(table(
-        predicted$novel_variable[!is.na(predicted$novelty) & predicted$novelty < 0]
-      ), decreasing = TRUE))[1]
-      sprintf("%d of %d cells (%d%%) are outside the training range, mostly %s",
-              outside, scored, max(1, round(100 * outside / scored)), culprit)
-    }
-
-    panels$novelty <- ggplot2::ggplot(
-      predicted, ggplot2::aes(x = .data$lon, y = .data$lat,
-                              fill = .data$novelty)) +
-      ggplot2::geom_raster() +
-      ggplot2::scale_fill_gradientn(
-        colours = c("#67001f", "#d6604d", "#fddbc7", "#e8eef2", "#7fa8c4",
-                    "#3a6b8a"),
-        # Zero sits where the warm colours end. rescale puts the break at the
-        # true zero of this month's range rather than at the midpoint of it,
-        # so the boundary means the same thing on every map in a run.
-        values = novelty_scale_positions(predicted$novelty),
-        na.value = "white", name = "Similarity"
-      ) +
-      ggplot2::labs(subtitle = subtitle)
+  if ("novelty" %in% drawable) {
+    # Already a MESS surface, so it is passed as one rather than recomputed
+    # from the covariates - which is what `training = NULL` means here.
+    panels$novelty <- fancyfx::plotExtrapolation(
+      projection_raster(predicted, "novelty"), legend.lab = "Similarity"
+    ) + ggplot2::labs(subtitle = novelty_subtitle(predicted))
   }
 
-  if (length(panels) == 0) return(invisible(NULL))
-
-  panels <- lapply(panels, function(p) {
-    p + ggplot2::coord_quickmap() +
-      ggplot2::labs(x = NULL, y = NULL) +
-      ggplot2::theme_bw() +
-      ggplot2::theme(panel.grid = ggplot2::element_blank(),
-                     legend.position = "bottom")
-  })
-
   title <- paste0(species, " - ", month.name[month], " ", year,
                   ": how far to trust it")
 
   # Stacked with a shared title rather than side by side, so each panel keeps
-  # the aspect ratio the coastline has and neither is squashed. patchwork is a
+  # the aspect ratio the coastline has and none is squashed. patchwork is a
   # Suggests, and its `/` operator is the whole reason it is wanted here - so
   # without it, fall back to the panel that carries more of the answer rather
   # than failing. Novelty is that panel: the spread cannot tell you a cell is
@@ -142,37 +106,47 @@ plot_projection_uncertainty <- function(predicted, year, month, species, path) {
   }
 
   ggplot2::ggsave(path, plot = combined, width = 7,
-                  height = if (stacked) 11 else 7, dpi = 150)
+                  height = if (stacked) 4 + 3.5 * length(panels) else 7,
+                  dpi = 150)
   invisible(path)
 }
 
-#' Where zero falls on a novelty colour ramp
+#' One column of a projection, as a raster fancyfx can draw
 #'
-#' `scale_fill_gradientn()` places its colours at positions in `[0, 1]` across
-#' the data range, so zero — the only value on a novelty surface that means
-#' anything fixed — lands somewhere different on every map unless it is put
-#' there deliberately. This returns positions that pin the warm-to-cool break to
-#' the true zero, so red means extrapolated on every month of a run rather than
-#' meaning "low for this month".
+#' The projection is a table of cells at this point in the pipeline, and every
+#' fancyfx map takes a `SpatRaster`. The cells are already on a regular grid —
+#' `project_patch_model()` builds the GeoTIFF from the same frame — so this is
+#' a change of container rather than any resampling.
 #'
-#' @param novelty the novelty values being plotted
-#' @return a vector of six positions in `[0, 1]`, for the six ramp colours
+#' @param predicted a projection from `predict_grid()`
+#' @param column which column to rasterize
+#' @return a single-layer `SpatRaster`
 #' @keywords internal
-novelty_scale_positions <- function(novelty) {
-  finite <- novelty[is.finite(novelty)]
-  # Nothing to anchor: an all-positive or empty surface gets an even ramp.
-  if (length(finite) == 0) return(seq(0, 1, length.out = 6))
-
-  low <- min(finite)
-  high <- max(finite)
-  if (high <= low) return(seq(0, 1, length.out = 6))
-
-  zero <- (0 - low) / (high - low)
-  # Clamped off the ends so the break stays visible when everything is on one
-  # side of zero, which is the common case.
-  zero <- min(max(zero, 0.02), 0.98)
+projection_raster <- function(predicted, column) {
+  out <- terra::rast(as.data.frame(predicted[c("lon", "lat", column)]),
+                     type = "xyz", crs = "EPSG:4326")
+  names(out) <- column
+  out
+}
 
-  c(0, zero * 0.5, zero * 0.95, zero, zero + (1 - zero) * 0.5, 1)
+#' The one line that says whether a novelty panel needed looking at
+#'
+#' @param predicted a projection carrying `novelty` and `novel_variable`
+#' @return a subtitle string
+#' @keywords internal
+novelty_subtitle <- function(predicted) {
+  outside <- sum(!is.na(predicted$novelty) & predicted$novelty < 0)
+  scored <- sum(!is.na(predicted$novelty))
+  if (outside == 0) return("Every cell is inside the training range")
+
+  culprit <- if ("novel_variable" %in% names(predicted)) {
+    names(sort(table(
+      predicted$novel_variable[!is.na(predicted$novelty) & predicted$novelty < 0]
+    ), decreasing = TRUE))[1]
+  }
+  sprintf("%d of %d cells (%d%%) are outside the training range%s",
+          outside, scored, max(1, round(100 * outside / scored)),
+          if (is.null(culprit)) "" else paste(", mostly", culprit))
 }
 
 #' Build a leaflet map of a projection GeoTIFF
diff --git a/R/power.R b/R/power.R
index 711a250..fae105d 100644
--- a/R/power.R
+++ b/R/power.R
@@ -485,3 +485,117 @@ power_point <- function(model_data, configs, predictors, types, fraction,
     stringsAsFactors = FALSE
   )
 }
+
+#' How much two projected surfaces agree about the habitat
+#'
+#' [compare_runs()] asks which run scores better. This asks something the
+#' scores cannot: whether the two maps are the *same map*. Two runs can rank
+#' stations equally well and disagree completely about where the habitat is,
+#' and a comparison of AUCs would report them as equivalent.
+#'
+#' Schoener's *D* and Warren's *I*, both on 0 to 1. **1** is identical
+#' surfaces; **0** is no overlap at all. There is no threshold at which two
+#' niches become "the same" — the number is a description, and what counts as
+#' close depends on what the maps are for.
+#'
+#' @section Read them as high numbers:
+#' Both statistics compare the suitability *distributions* cell by cell after
+#' normalising them, so they run high. Two surfaces of independent random noise
+#' over the same grid score around 0.7 on *D* and 0.9 on *I* — not because they
+#' agree about anything, but because both spread their probability over the
+#' same cells in similar proportions. Treat those as the floor for surfaces of
+#' this shape rather than reading 0.9 as near-identity, and compare overlaps
+#' against each other rather than against 1.
+#'
+#' @section What it is good for:
+#' Two questions, mostly. *Do two species occupy the same habitat?* — run
+#' `cfin` and `ctyp` over the same months and compare their surfaces. And *does
+#' the choice of algorithm change the map?* — the more useful one when a
+#' comparison of scores came back inconclusive, because a small difference in
+#' AUC with a low overlap means the two models are doing genuinely different
+#' things and the tie in performance is hiding it.
+#'
+#' @section Matching the cells:
+#' The two surfaces must be the same cells in the same order, or the overlap is
+#' measured between places that have nothing to do with each other. Cells are
+#' matched on their coordinates rather than on row position, and cells present
+#' in one projection and not the other are dropped with a warning — a
+#' projection that covers less ground usually has a reason, and averaging over
+#' the difference would hide it.
+#'
+#' @param x,y projections to compare: data frames with `lon`, `lat` and
+#'   `suitability`, as [project_patch_model()] writes to `suitability.csv`
+#' @param statistic `"both"`, `"D"`, or `"I"`
+#' @param digits how many decimal places the coordinates are matched on.
+#'   Projections from the same grid agree exactly; this is for two that came
+#'   off slightly different pipelines
+#' @return a named numeric vector, with an `n_cells` attribute
+#' @examples
+#' \dontrun{
+#' rf <- readr::read_csv("output/rf/projections/suitability.csv")
+#' gam <- readr::read_csv("output/gam/projections/suitability.csv")
+#'
+#' projection_overlap(rf[rf$month == 7, ], gam[gam$month == 7, ])
+#' }
+#' @references
+#' Schoener TW (1968). The *Anolis* lizards of Bimini: resource partitioning in
+#' a complex fauna. *Ecology* **49**(4), 704-726. \doi{10.2307/1935534} — *D*
+#'
+#' Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus
+#' conservatism: quantitative approaches to niche evolution. *Evolution*
+#' **62**(11), 2868-2883. \doi{10.1111/j.1558-5646.2008.00482.x} — *I*
+#' @seealso [compare_runs()], which asks which surface is *better* rather than
+#'   whether they are the same
+#' @export
+projection_overlap <- function(x, y, statistic = c("both", "D", "I"),
+                               digits = 6) {
+  statistic <- match.arg(statistic)
+  x <- projection_cells(x, "x")
+  y <- projection_cells(y, "y")
+
+  key <- function(cells) {
+    paste(round(cells$lon, digits), round(cells$lat, digits), sep = "\r")
+  }
+  shared <- intersect(key(x), key(y))
+  dropped <- length(union(key(x), key(y))) - length(shared)
+
+  if (length(shared) < 2) {
+    stop("The two projections share ", length(shared), " cells, so there is ",
+         "no common ground to measure overlap over.\nThey were probably ",
+         "projected onto different grids; `covariates.grid` decides that.",
+         call. = FALSE)
+  }
+  if (dropped > 0) {
+    warning("The two projections do not cover the same cells: ", dropped,
+            " of ", length(shared) + dropped, " are in one and not the other, ",
+            "and the overlap uses the ", length(shared), " they share.",
+            call. = FALSE)
+  }
+
+  out <- fancyfx::niche_overlap(
+    x$suitability[match(shared, key(x))],
+    y$suitability[match(shared, key(y))],
+    statistic = statistic
+  )
+  attr(out, "n_cells") <- length(shared)
+  out
+}
+
+#' The cells of a projection, checked
+#'
+#' @param projection a projection table
+#' @param label which argument it came from, for the error
+#' @return a data frame of `lon`, `lat` and `suitability`
+#' @keywords internal
+projection_cells <- function(projection, label) {
+  projection <- as.data.frame(projection)
+  needed <- c("lon", "lat", "suitability")
+  missing <- setdiff(needed, names(projection))
+  if (length(missing) > 0) {
+    stop("`", label, "` is missing: ", paste(missing, collapse = ", "),
+         ".\nprojection_overlap() reads the table project_patch_model() ",
+         "writes to suitability.csv; one month of it at a time, since two ",
+         "months stacked are two surfaces rather than one.", call. = FALSE)
+  }
+  projection[needed]
+}
diff --git a/R/spatial_bias.R b/R/spatial_bias.R
new file mode 100644
index 0000000..b76611c
--- /dev/null
+++ b/R/spatial_bias.R
@@ -0,0 +1,133 @@
+#' How much of the AUC is an artefact of where the stations are
+#'
+#' A cross-validated AUC is optimistic when the held-out patches sit closer to
+#' the training patches than the held-out non-patches do. The model can then
+#' score well by recognising *where* it has already been rather than *what*
+#' makes a patch, and nothing in the ROC curve distinguishes the two.
+#'
+#' This is the spatial sorting bias of Hijmans (2012): the mean distance from
+#' each held-out patch to the nearest training patch, over the same distance
+#' for the held-out non-patches.
+#'
+#' \deqn{SSB = \frac{\overline{d}(\mathrm{held\ out\ patch},\ \mathrm{training\ patch})}{\overline{d}(\mathrm{held\ out\ non\ patch},\ \mathrm{training\ patch})}}
+#'
+#' * **Near 1** — held-out patches and non-patches are equally far from the
+#'   training patches. The split is spatially fair and the AUC means what it
+#'   appears to.
+#' * **Near 0** — held-out patches are much closer to training patches than the
+#'   non-patches are. **The AUC is inflated**, and by an amount this number
+#'   does not tell you.
+#'
+#' It is a property of the *split*, not of the model. Every model type fitted on
+#' the same folds gets the same value, which is the point: it says how much to
+#' trust the comparison between them, not which of them won.
+#'
+#' @section Why an ECOMON run should look at it:
+#' Survey stations are not scattered at random. They sit on transects, revisited
+#' season after season, so a randomly held-out station usually has a near
+#' neighbour in the training folds — often the same station in another year.
+#' That is the situation this measures, and random `model.cv_folds` cannot avoid
+#' it. A low value is an argument for spatially blocked folds, not for a
+#' different model.
+#'
+#' @param model a fitted model from [fit_patch_model()] or
+#'   [fit_patch_ensemble()], carrying `coordinates` and held-out `predictions`
+#' @param geo whether to measure great-circle distances. `TRUE` by default,
+#'   because the coordinates are longitude and latitude, and a degree of
+#'   longitude is not a degree of latitude anywhere but the equator
+#' @return a data frame with one row per fold and a `fold` of `"overall"` for
+#'   the mean, carrying `patch_distance`, `non_patch_distance` and `ssb`
+#' @examples
+#' \dontrun{
+#' model <- fit_patch_model(dat, config)
+#' spatial_bias(model)
+#' }
+#' @references
+#' Hijmans RJ (2012). Cross-validation of species distribution models: removing
+#' spatial sorting bias and calibration with a null model. *Ecology* **93**(3),
+#' 679-688. \doi{10.1890/11-0826.1}
+#' @seealso [evaluation_table()], which reports the summary beside the AUC it
+#'   qualifies
+#' @export
+spatial_bias <- function(model, geo = TRUE) {
+  coordinates <- model$coordinates
+  predictions <- model$predictions
+
+  if (is.null(coordinates) || is.null(predictions) ||
+      !all(c(".row", "id", "patch") %in% names(predictions))) {
+    return(NULL)
+  }
+
+  folds <- split(seq_len(nrow(predictions)), predictions$id)
+  rows <- lapply(names(folds), function(fold) {
+    at <- folds[[fold]]
+    held_out <- predictions$.row[at]
+    is_patch <- predictions$patch[at] == "patch"
+
+    # The reference is the training patches: every patch station except the
+    # ones held out in this fold. A held-out patch that is close to one of them
+    # is a patch the model has, in effect, already seen.
+    all_patches <- predictions$.row[predictions$patch == "patch"]
+    training <- setdiff(all_patches, held_out[is_patch])
+
+    if (!any(is_patch) || !any(!is_patch) || length(training) == 0) {
+      return(NULL)
+    }
+
+    out <- fancyfx::spatial_sorting_bias(
+      presence = coordinates[held_out[is_patch], , drop = FALSE],
+      absence = coordinates[held_out[!is_patch], , drop = FALSE],
+      reference = coordinates[training, , drop = FALSE],
+      geo = geo
+    )
+    data.frame(fold = fold, patch_distance = out[["presence"]],
+               non_patch_distance = out[["absence"]], ssb = out[["ssb"]],
+               stringsAsFactors = FALSE)
+  })
+
+  out <- do.call(rbind, Filter(Negate(is.null), rows))
+  if (is.null(out) || nrow(out) == 0) return(NULL)
+
+  # One number to report beside the AUC, and the folds it came from kept so a
+  # single odd fold is visible rather than averaged into the answer.
+  overall <- data.frame(
+    fold = "overall",
+    patch_distance = mean(out$patch_distance, na.rm = TRUE),
+    non_patch_distance = mean(out$non_patch_distance, na.rm = TRUE),
+    ssb = mean(out$ssb, na.rm = TRUE),
+    stringsAsFactors = FALSE
+  )
+  rbind(out, overall)
+}
+
+#' The one-line reading of a spatial sorting bias
+#'
+#' Written for the run log, where a bare ratio would be ignored. The thresholds
+#' are conventional rather than derived — Hijmans (2012) offers no cutoff, and
+#' presenting one as though he did would be worse than a rule of thumb labelled
+#' as one.
+#'
+#' @param ssb the overall spatial sorting bias
+#' @return a single string
+#' @keywords internal
+spatial_bias_note <- function(ssb) {
+  if (is.na(ssb)) return("spatial sorting bias could not be computed")
+  if (ssb >= 0.8) {
+    return(paste0("spatial sorting bias ", signif(ssb, 2),
+                  ": the folds are spatially fair"))
+  }
+  paste0("spatial sorting bias ", signif(ssb, 2),
+         ": held-out patches sit closer to training patches than non-patches ",
+         "do, so the AUC above is optimistic. Consider spatially blocked folds")
+}
+
+#' The single spatial sorting bias from a per-fold table
+#'
+#' @param bias the result of [spatial_bias()], or `NULL`
+#' @return the overall ratio, or `NA_real_`
+#' @keywords internal
+overall_ssb <- function(bias) {
+  if (is.null(bias) || !("fold" %in% names(bias))) return(NA_real_)
+  value <- bias$ssb[bias$fold == "overall"]
+  if (length(value) == 1) value else NA_real_
+}
diff --git a/R/uncertainty.R b/R/uncertainty.R
index b618789..e75e35a 100644
--- a/R/uncertainty.R
+++ b/R/uncertainty.R
@@ -231,6 +231,14 @@ ensemble_spread <- function(ensemble, newdata, level = 0.9) {
 #' because its chlorophyll is higher than anything a station saw" is a decision
 #' about whether to widen the training window or clip the map.
 #'
+#' @section Where the numbers come from:
+#' [fancyfx::mess()] computes both columns. taupatch used to implement MESS
+#' itself, and the two implementations agreed to the last decimal on every case
+#' tested -- which is the argument for there being one of them rather than two.
+#' What survives here is the interface: the `novelty` and `novel_variable`
+#' names, which the projection CSV and the GeoTIFF layers are written under and
+#' which a run\'s output should not lose to an internal tidy-up.
+#'
 #' @param grid the cells to score, with one column per predictor
 #' @param model_data the data the model was fitted on
 #' @param predictors predictor column names
@@ -256,67 +264,16 @@ novelty_surface <- function(grid, model_data, predictors) {
                       stringsAsFactors = FALSE))
   }
 
-  similarity <- vapply(predictors, function(v) {
-    variable_similarity(grid[[v]], model_data[[v]])
-  }, numeric(nrow(grid)))
-  # vapply drops to a vector when the grid has one row, which then indexes
-  # wrongly below.
-  similarity <- matrix(similarity, nrow = nrow(grid),
-                       dimnames = list(NULL, predictors))
-
-  worst <- apply(similarity, 1, function(row) {
-    if (all(is.na(row))) return(NA_integer_)
-    which.min(row)
-  })
-
-  data.frame(
-    novelty = apply(similarity, 1, min, na.rm = FALSE),
-    novel_variable = ifelse(is.na(worst), NA_character_, predictors[worst]),
-    stringsAsFactors = FALSE
-  )
-}
-
-#' Similarity of values to a training distribution, for one predictor
-#'
-#' The per-variable half of [novelty_surface()]. Negative below the training
-#' minimum and above its maximum, scaled by the training range; inside, twice
-#' the distance to the nearer tail in percentile terms, so the median scores 100.
-#'
-#' @param values the values to score
-#' @param train the training values for the same predictor
-#' @return a numeric vector the length of `values`
-#' @keywords internal
-variable_similarity <- function(values, train) {
-  train <- sort(train[is.finite(train)])
-  n <- length(train)
-  if (n == 0) return(rep(NA_real_, length(values)))
-
-  low <- train[1]
-  high <- train[n]
-  span <- high - low
-
-  # A predictor that never varied in training carries no information about
-  # what is ordinary, so it can only say same-or-not.
-  if (span == 0) {
-    return(ifelse(is.na(values), NA_real_, ifelse(values == low, 100, -Inf)))
-  }
-
-  # Percentile of each value within the training distribution. findInterval
-  # with left.open counts the training values strictly below, which is the
-  # count this needs and is a binary search rather than a full comparison
-  # against every training row - the difference between a grid that scores in
-  # a moment and one that does not.
-  below <- findInterval(values, train, left.open = TRUE)
-  percentile <- 100 * below / n
+  # The arithmetic is fancyfx's. This function is the name taupatch's output
+  # columns are written under, and the place the scale is explained; keeping it
+  # as a thin call means the projection CSV and the GeoTIFF layer keep the
+  # names they have always had.
+  out <- fancyfx::mess(as.data.frame(grid[predictors]),
+                       as.data.frame(model_data[predictors]),
+                       vars = predictors, limiting = TRUE)
 
-  out <- ifelse(
-    percentile == 0, 100 * (values - low) / span,
-    ifelse(percentile <= 50, 2 * percentile,
-           ifelse(percentile < 100, 2 * (100 - percentile),
-                  100 * (high - values) / span))
-  )
-  out[is.na(values)] <- NA_real_
-  out
+  data.frame(novelty = out$mess, novel_variable = out$mess_variable,
+             stringsAsFactors = FALSE)
 }
 
 #' Summarise a novelty surface for the run log
diff --git a/README.md b/README.md
index 7f89191..f7ec6c2 100644
--- a/README.md
+++ b/README.md
@@ -145,6 +145,9 @@ remotes::install_github("chross22/datamatch")
 `datamatch` no longer depends on `BigelowLab/copernicus` — it calls the
 Copernicus Marine Toolbox directly — so that package no longer needs installing.
 
+`fancyfx` draws the effect and uncertainty figures and is installed with the
+package, from the `Remotes` field, so it needs no separate step.
+
 Derived covariates (`covariates.derivoce`) additionally need:
 
 ```r
@@ -878,14 +881,10 @@ On top of that, each model contributes what only it can:
 | `gam` | Effective degrees of freedom per smooth. An `edf` of 1 means the smooth collapsed to a line |
 | `rf` / `brt` | None. The partial effect curve *is* their answer |
 
-With [`fancyfx`](https://github.com/chross22/fancyfx) installed, a GAM also
-gets its **fitted smooths** drawn — each term with its standard error band and a
-rug showing where the data actually is. Those carry uncertainty, which a partial
-dependence curve cannot:
-
-```r
-remotes::install_github("chross22/fancyfx")
-```
+A GAM also gets its **fitted smooths** drawn by
+[`fancyfx`](https://github.com/chross22/fancyfx) — each term with its standard
+error band and a rug showing where the data actually is. Those carry
+uncertainty, which a partial dependence curve cannot.
 
 Their x axes read in standard deviations, because the smooths belong to the model
 and the model was fitted on the recipe's output. Set `covariates.normalize: false`
@@ -1256,7 +1255,7 @@ diagnostics/cv_predictions.csv     held-out predictions, for any metric not tabu
 diagnostics/partial_effects.png    what each predictor does to patch probability
 diagnostics/coefficients.png       glm only: signed effects with intervals
 diagnostics/smooth_terms.csv       gam only: effective degrees of freedom per smooth
-diagnostics/gam_smooths.png        gam only, with fancyfx: fitted smooths with error bands
+diagnostics/gam_smooths.png        gam only: fitted smooths with error bands, drawn by fancyfx
 diagnostics/members/<type>/        with model.ensemble: the above, one directory per algorithm
 projections/suitability.csv       every cell of every month: species, year, month, lon, lat, probability
                                   plus the interval and novelty columns, with projection.uncertainty
@@ -1514,6 +1513,19 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`.
   Learning* **52**(3), 239–281.
   [doi:10.1023/A:1024068626366](https://doi.org/10.1023/A:1024068626366) — the
   variance correction behind the jackknife's `p_value` and `compare_runs()`
+- Hijmans RJ (2012). Cross-validation of species distribution models: removing
+  spatial sorting bias and calibration with a null model. *Ecology* **93**(3),
+  679–688. [doi:10.1890/11-0826.1](https://doi.org/10.1890/11-0826.1) — the
+  `ssb` row in `evals.csv`
+- Schoener TW (1968). The *Anolis* lizards of Bimini: resource partitioning in
+  a complex fauna. *Ecology* **49**(4), 704–726.
+  [doi:10.2307/1935534](https://doi.org/10.2307/1935534) — Schoener's *D* in
+  `projection_overlap()`
+- Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus
+  conservatism: quantitative approaches to niche evolution. *Evolution*
+  **62**(11), 2868–2883.
+  [doi:10.1111/j.1558-5646.2008.00482.x](https://doi.org/10.1111/j.1558-5646.2008.00482.x)
+  — Warren's *I*
 - Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of
   power calculations for data analysis. *The American Statistician* **55**(1),
   19–24.
diff --git a/docs/fancyfx_plan.md b/docs/fancyfx_plan.md
new file mode 100644
index 0000000..61c971f
--- /dev/null
+++ b/docs/fancyfx_plan.md
@@ -0,0 +1,220 @@
+# Taking what fancyfx now offers
+
+`fancyfx` has gone from 0.2.0 to 0.8.0 and now exports 26 functions, a good
+number of which overlap what taupatch built for itself. This is the plan for
+adopting it: what to take, what to leave, what becomes redundant, and the one
+thing that cannot be taken as it stands.
+
+Written against `fancyfx 0.8.0` (`d429991`) and taupatch 0.2.0 (`983cf55`).
+
+## The finding that shapes everything else
+
+**fancyfx's four evaluation plots cannot consume taupatch's cross-validated
+predictions, and adopting them as they stand would replace honest diagnostics
+with optimistic ones.**
+
+This was worth checking rather than assuming, and the answer is not the one the
+function signatures suggest. `plotROC()`, `plotThreshold()`, `plotImportance()`
+and `plotCalibration()` all route through `fancyfx:::evaluation_pairs()`, which
+does this and only this:
+
+```r
+observed  <- binary_response(model, newdata)
+predicted <- predict_probability(model, newdata, ...)
+```
+
+There is no path to hand it predictions that already exist. It re-predicts,
+always.
+
+taupatch's ROC, threshold, calibration and PR curves are drawn from **pooled
+out-of-fold predictions**: every station is scored by the one fold model that
+never saw it, collected in `model$predictions` with `.row`, `id`, `patch` and
+`.pred_patch`. That is cross-validated performance. Handing the final fitted
+model and its own training stations to fancyfx instead produces something
+different — and fancyfx knows it, which is to its credit:
+
+* `evaluation_pairs()` warns: *"Evaluating on the data the model was fitted to.
+  These metrics are optimistic and are not validation."*
+* `in_sample_caption()` prints it onto the figure: *"In-sample: scored on the
+  data the model was fitted to. Optimistic, and not validation."*
+
+The `folds` argument does not rescue this. It splits the *scoring* of one
+already-fitted model by fold; it does not make each fold's predictions come from
+a model that held that fold out. On taupatch's data `is_training_data()` returns
+`TRUE`, so every plot would carry the caption above.
+
+A second, blunter obstacle sits on top of the first. taupatch fits through a
+`recipes` workflow, so the engine's response is the baked `..y`, not `patch`.
+Handing `model_data` to `threshold_metrics()` does not produce an optimistic
+number — it errors outright:
+
+```
+newdata has no column '..y', the model's response.
+```
+
+So these four cannot be adopted **as fancyfx currently stands**. The fix is not
+in taupatch, and it is not a workaround: it is a generalisation fancyfx wants
+anyway.
+
+### Upstream change 1: let evaluation take predictions
+
+An entry point that accepts predictions rather than a model — the shape any
+cross-validated workflow already has, taupatch included:
+
+```r
+threshold_metrics(observed = , predicted = , folds = )
+```
+
+`evaluation_pairs()` returns early when handed those, with `in.sample = FALSE`,
+and every plot built on it works unchanged. This is the whole change. It costs
+fancyfx nothing, it removes the `..y` problem along with the in-sample one, and
+it makes fancyfx usable by *any* package that cross-validates rather than only
+by ones that keep a single fitted model around.
+
+Only with this in place do the four evaluation plots become a streamlining:
+taupatch deletes four hand-built ggplots and keeps its out-of-fold numbers.
+Without it they are a downgrade, and would not be worth taking.
+
+## What to take now
+
+### New capability — nothing in taupatch does these
+
+| function | what it adds | why it belongs here |
+|---|---|---|
+| `spatial_sorting_bias()` | Hijmans (2012). How much of the AUC is an artefact of presences and absences being differently distributed in space | ECOMON stations are clustered along transects. A spatially sorted sample inflates AUC, and nothing in the current evaluation says by how much. This is a genuine gap |
+| `niche_overlap()` | Schoener's *D* and Warren's *I* between two projected surfaces | taupatch models three species and can now compare runs. "Do *cfin* and *ctyp* occupy the same habitat?" and "do rf and gam draw the same map?" are the natural next questions, and `compare_runs()` deliberately answers neither — it compares performance, not surfaces |
+| `niche_equivalency()` | A permutation test for whether two niches are the same | The inferential half of the above |
+
+### Better versions of what exists
+
+| function | replaces | why it is a clean swap |
+|---|---|---|
+| `plotUncertainty()` | the spread panel of `plot_projection_uncertainty()` | Takes a raster, not a model, so the in-sample problem does not arise. Adds `cv`/`range`/`iqr` beside `sd`, and `max.cells` downsampling, which the current panel needs on a real grid |
+| `plotExtrapolation()` | the novelty panel of the same | Same. Adds `novel.only` |
+| `comparePlots()` | nothing yet | Overlays several models' effect curves for one covariate — exactly the per-member view the new ensemble wants |
+| `mess()` | `novelty_surface()` | The same method implemented twice. Requires upstream change 2 below |
+
+### Upstream change 2: let mess() take a data frame, and name the culprit
+
+`fancyfx::mess()` requires a `SpatRaster` and returns the surface alone:
+
+```
+Error: x must be a SpatRaster of covariates, not a <data.frame>.
+```
+
+taupatch's projections are data frames at that point in the pipeline, and its
+`novelty_surface()` returns a second column — `novel_variable`, the predictor
+that put the cell outside the training range. That column is the actionable
+half: "this shelf is extrapolated" is a shrug, "extrapolated because its
+chlorophyll is higher than any station saw" is a decision. It is used in the run
+log, the projection CSV, and the plot subtitle.
+
+Two small generalisations, then, and both make `mess()` better on its own terms:
+accept a data frame as well as a raster, and return which variable was
+responsible. With those, `novelty_surface()` deletes and taupatch calls
+`mess()`.
+
+### Deliberately not taken
+
+Excluded by the "streamlining **and** improvement" test rather than by any fault
+of theirs:
+
+* **`thin_points()`** — spatial thinning changes which stations are modelled.
+  That is a methodological decision for a study to make, not a default for a
+  package to acquire, and it streamlines nothing here.
+* **`calc_deviance()`** — defined for `glm` and `gam` and not for `rf` or
+  `brt`. taupatch reports metrics that mean the same thing for all four types
+  on purpose; a column that is `NA` for half of them works against that.
+* **`plotHexbin()`** — the station map is not currently a problem worth a new
+  code path.
+
+## Redundancy assessment
+
+To be done **after** the integration lands, against the code as it then stands,
+not against this plan. The candidates, in the order they are worth examining:
+
+1. **`novelty_surface()` against `fancyfx::mess()`.** These are the same method
+   — the multivariate environmental similarity surface of Elith et al. (2010) —
+   implemented twice. One should go. taupatch's returns `novel_variable`, the
+   predictor responsible, which is the actionable half and must survive the
+   merge whichever way it goes.
+2. **`plot_projection_uncertainty()`** becomes a thin arrangement of two
+   fancyfx panels rather than two hand-built ggplots.
+3. **`ensemble_spread()` against `fancyfx::ensemble_summary()`.** Same
+   statistics over a stack of member predictions. taupatch's returns a data
+   frame aligned to the projection rows; fancyfx's takes a raster. Whether they
+   can be one function depends on which side the alignment lives.
+4. **`permutation_importance()`** — taupatch's is deliberately model-agnostic
+   and computed by prediction so all four model types are comparable. fancyfx's
+   takes `(model, newdata)` and has the same in-sample problem as the plots.
+   Blocked with them, and possibly not worth merging even after.
+5. **`plot_importance()`, `plot_calibration()`, `plot_roc_curve()`,
+   `plot_threshold_performance()`** — not candidates. Their fancyfx
+   counterparts answer a different question on different evidence.
+
+The rule for the assessment: **a duplicate is only redundant if the replacement
+answers the same question on the same evidence.** `mess()` and
+`novelty_surface()` do. `plotROC()` and `plot_roc_curve()` currently do not,
+because one is in-sample and the other is not, and that is a difference in what
+the number means rather than in how it is drawn.
+
+## The app
+
+Into the existing **Diagnostics** tab rather than a new one, since that is where
+a user already goes for this material:
+
+* spatial sorting bias beside the evaluation table, with the one-line
+  explanation of what a high value means for the AUC above it
+* explained deviance in the same block
+* the ensemble's per-member effect curves, via `comparePlots()`, when the run
+  fitted an ensemble
+* niche overlap between the active species and one other, when more than one has
+  been run
+
+The app must not gain a dependency the package does not have: `fancyfx` is a
+Suggests, so every panel here is behind `has_fancyfx()` and degrades to the
+current view without it.
+
+## Order of work
+
+Two of these are changes to fancyfx. They come first, because the taupatch side
+of each is small and pointless without them.
+
+**In fancyfx**
+
+1. `threshold_metrics(observed=, predicted=, folds=)` — upstream change 1.
+2. `mess()` on a data frame, returning the responsible variable — upstream
+   change 2.
+
+**In taupatch**
+
+3. `plotUncertainty()` / `plotExtrapolation()` into
+   `plot_projection_uncertainty()`. Needs nothing upstream; two hand-built
+   ggplots delete.
+4. `comparePlots()` for ensemble members.
+5. `spatial_sorting_bias()` into the evaluation table.
+6. `niche_overlap()` / `niche_equivalency()` as a projection comparison,
+   alongside `compare_runs()`.
+7. `novelty_surface()` replaced by `mess()`, once (2) lands.
+8. The four evaluation plots replaced, once (1) lands.
+9. App: Diagnostics tab.
+10. Redundancy pass, against the code as it then stands.
+The four evaluation plots are not on this list, and are not expected to join
+it.
+
+## References
+
+Elith J, Kearney M, Phillips S (2010). The art of modelling range-shifting
+species. *Methods in Ecology and Evolution* **1**(4), 330–342.
+doi:10.1111/j.2041-210X.2010.00036.x — MESS
+
+Hijmans RJ (2012). Cross-validation of species distribution models: removing
+spatial sorting bias and calibration with a null model. *Ecology* **93**(3),
+679–688. doi:10.1890/11-0826.1 — spatial sorting bias
+
+Schoener TW (1968). The *Anolis* lizards of Bimini: resource partitioning in a
+complex fauna. *Ecology* **49**(4), 704–726. doi:10.2307/1935534 — *D*
+
+Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus
+conservatism. *Evolution* **62**(11), 2868–2883.
+doi:10.1111/j.1558-5646.2008.00482.x — *I*, and the equivalency test
diff --git a/inst/shiny/app.R b/inst/shiny/app.R
index ac3e084..325184e 100644
--- a/inst/shiny/app.R
+++ b/inst/shiny/app.R
@@ -920,7 +920,15 @@ server <- function(input, output, session) {
                                             input$bathymetry %||% character())
     config$covariates$selected <- union(input$covariates %||% character(),
                                         ingredients)
-    config$covariates$exclude <- ingredients
+    # A derived covariate can also be an ingredient: the gradient of current
+    # speed computes the speed on the way, and that column is no more a
+    # predictor than the velocity components behind it are.
+    config$covariates$exclude <- union(
+      ingredients,
+      derivoce_dependency_columns(input$derived %||% character(),
+                                  input$covariates %||% character(),
+                                  input$bathymetry %||% character())
+    )
     config$covariates$bathymetry <- input$bathymetry %||% character()
     config$covariates$climate <- input$climate %||% character()
 
@@ -1210,9 +1218,9 @@ server <- function(input, output, session) {
   # A GAM has its own partial effects, and they are better than the generic
   # ones: read out of the fitted model rather than reconstructed by prediction,
   # so they carry the uncertainty a partial dependence curve cannot.
+  # fancyfx is an Imports now, so this is only ever about the model type.
   use_fancyfx <- reactive({
-    identical(run_result()$model$type, "gam") &&
-      requireNamespace("fancyfx", quietly = TRUE)
+    identical(run_result()$model$type, "gam")
   })
 
   output$partial_effects_note <- renderUI({
@@ -1256,8 +1264,8 @@ server <- function(input, output, session) {
       ))
     }
     if (identical(model$type, "gam")) {
-      # The smooths themselves are the partial effects panel above when
-      # fancyfx is present, so they are not repeated here.
+      # The smooths themselves are the partial effects panel above, so they
+      # are not repeated here.
       return(tagList(
         h4("Smooth terms"),
         helpText("Effective degrees of freedom per smooth. An edf of 1 means",
diff --git a/man/derivoce_dependency_columns.Rd b/man/derivoce_dependency_columns.Rd
new file mode 100644
index 0000000..d9a6dd3
--- /dev/null
+++ b/man/derivoce_dependency_columns.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/derivoce.R
+\name{derivoce_dependency_columns}
+\alias{derivoce_dependency_columns}
+\title{Derived columns pulled in only as ingredients}
+\usage{
+derivoce_dependency_columns(ids, selected, bathymetry = character())
+}
+\arguments{
+\item{ids}{chosen derived covariate ids}
+
+\item{selected}{time-varying covariate names}
+
+\item{bathymetry}{static seafloor covariate names}
+}
+\value{
+character vector of derived column names, possibly empty
+}
+\description{
+Choosing the gradient of current speed computes the speed on the way, and
+that column then sits in the modelling data looking exactly like one that was
+asked for. It was not: the ingredient of a derived covariate is no more a
+predictor than the velocity components behind an FSLE are.
+}
+\details{
+This names the derived columns a selection produced without anyone choosing
+them, so they can go into \code{covariates.exclude} alongside the downloads that
+\code{\link[=derivoce_required_inputs]{derivoce_required_inputs()}} finds. A column named here is still computed and
+still available to anything later that reads it; it just does not become a
+predictor.
+}
+\examples{
+# Asking for the gradient of current speed computes the speed too, and that
+# is an ingredient rather than a request.
+derivoce_dependency_columns("speed_grad", c("SST", "SSS"))
+
+# Asking for both makes the speed a request, so it is not excluded.
+derivoce_dependency_columns(c("speed", "speed_grad"), c("SST", "SSS"))
+}
+\seealso{
+\code{\link[=derivoce_required_inputs]{derivoce_required_inputs()}}, which does the same for downloads
+}
diff --git a/man/evaluation_table.Rd b/man/evaluation_table.Rd
index 3aaa506..fd32c5e 100644
--- a/man/evaluation_table.Rd
+++ b/man/evaluation_table.Rd
@@ -4,7 +4,7 @@
 \alias{evaluation_table}
 \title{Assemble a self-explanatory evaluation table}
 \usage{
-evaluation_table(predictions, cv_metrics, cutoff, bounds = NULL)
+evaluation_table(predictions, cv_metrics, cutoff, bounds = NULL, ssb = NULL)
 }
 \arguments{
 \item{predictions}{held-out predictions from resampling}
diff --git a/man/has_fancyfx.Rd b/man/has_fancyfx.Rd
deleted file mode 100644
index a35d08d..0000000
--- a/man/has_fancyfx.Rd
+++ /dev/null
@@ -1,31 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plot_effects.R
-\name{has_fancyfx}
-\alias{has_fancyfx}
-\title{Whether fancyfx is available to draw smooths}
-\usage{
-has_fancyfx()
-}
-\value{
-\code{TRUE} when fancyfx is installed
-}
-\description{
-Its own function so the optional path can be exercised in tests without
-mocking \code{requireNamespace()} itself, which every package that loads a
-graphics device also goes through.
-}
-\section{It used to be called fancygam}{
-
-The package was renamed when it grew past GAMs. The rename is why this
-matters more than a find-and-replace: \code{chross22/fancygam} still resolves on
-GitHub, so \code{Remotes: chross22/fancygam} kept installing — but what it
-installs now declares \code{Package: fancyfx}, so \code{requireNamespace("fancygam")}
-returned \code{FALSE} on every fresh install and the smooth plots were skipped in
-silence. Anyone with the old package still sitting in their library saw
-nothing wrong.
-
-That is the failure mode to watch for here: this function gates a diagnostic
-rather than the run, so a wrong answer costs a plot and no error.
-}
-
-\keyword{internal}
diff --git a/man/novelty_scale_positions.Rd b/man/novelty_scale_positions.Rd
deleted file mode 100644
index 966855b..0000000
--- a/man/novelty_scale_positions.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plotting.R
-\name{novelty_scale_positions}
-\alias{novelty_scale_positions}
-\title{Where zero falls on a novelty colour ramp}
-\usage{
-novelty_scale_positions(novelty)
-}
-\arguments{
-\item{novelty}{the novelty values being plotted}
-}
-\value{
-a vector of six positions in \verb{[0, 1]}, for the six ramp colours
-}
-\description{
-\code{scale_fill_gradientn()} places its colours at positions in \verb{[0, 1]} across
-the data range, so zero — the only value on a novelty surface that means
-anything fixed — lands somewhere different on every map unless it is put
-there deliberately. This returns positions that pin the warm-to-cool break to
-the true zero, so red means extrapolated on every month of a run rather than
-meaning "low for this month".
-}
-\keyword{internal}
diff --git a/man/novelty_subtitle.Rd b/man/novelty_subtitle.Rd
new file mode 100644
index 0000000..5832b16
--- /dev/null
+++ b/man/novelty_subtitle.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plotting.R
+\name{novelty_subtitle}
+\alias{novelty_subtitle}
+\title{The one line that says whether a novelty panel needed looking at}
+\usage{
+novelty_subtitle(predicted)
+}
+\arguments{
+\item{predicted}{a projection carrying \code{novelty} and \code{novel_variable}}
+}
+\value{
+a subtitle string
+}
+\description{
+The one line that says whether a novelty panel needed looking at
+}
+\keyword{internal}
diff --git a/man/novelty_surface.Rd b/man/novelty_surface.Rd
index 9295ed5..4c076a1 100644
--- a/man/novelty_surface.Rd
+++ b/man/novelty_surface.Rd
@@ -40,6 +40,16 @@ half: "this coast is extrapolated" is a shrug, and "this coast is extrapolated
 because its chlorophyll is higher than anything a station saw" is a decision
 about whether to widen the training window or clip the map.
 }
+\section{Where the numbers come from}{
+
+\code{\link[fancyfx:mess]{fancyfx::mess()}} computes both columns. taupatch used to implement MESS
+itself, and the two implementations agreed to the last decimal on every case
+tested -- which is the argument for there being one of them rather than two.
+What survives here is the interface: the \code{novelty} and \code{novel_variable}
+names, which the projection CSV and the GeoTIFF layers are written under and
+which a run\'s output should not lose to an internal tidy-up.
+}
+
 \examples{
 train <- data.frame(SST = c(4, 8, 12, 16), CHL = c(0.2, 0.5, 1.0, 2.0))
 grid <- data.frame(SST = c(10, 25), CHL = c(0.6, 0.6))
diff --git a/man/overall_ssb.Rd b/man/overall_ssb.Rd
new file mode 100644
index 0000000..177ed53
--- /dev/null
+++ b/man/overall_ssb.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spatial_bias.R
+\name{overall_ssb}
+\alias{overall_ssb}
+\title{The single spatial sorting bias from a per-fold table}
+\usage{
+overall_ssb(bias)
+}
+\arguments{
+\item{bias}{the result of \code{\link[=spatial_bias]{spatial_bias()}}, or \code{NULL}}
+}
+\value{
+the overall ratio, or \code{NA_real_}
+}
+\description{
+The single spatial sorting bias from a per-fold table
+}
+\keyword{internal}
diff --git a/man/plot_gam_smooths.Rd b/man/plot_gam_smooths.Rd
index e0831bb..30cd0ee 100644
--- a/man/plot_gam_smooths.Rd
+++ b/man/plot_gam_smooths.Rd
@@ -24,8 +24,7 @@ rather than reconstructed by prediction, so it carries uncertainty, which a
 partial dependence curve cannot.
 }
 \details{
-Drawn by \href{https://github.com/chross22/fancyfx}{fancyfx}, which is a Suggests
-— a run without it still gets the generic partial effect curves.
+Drawn by \href{https://github.com/chross22/fancyfx}{fancyfx}.
 }
 \section{Why the axes read in standard deviations}{
 
diff --git a/man/plot_projection_uncertainty.Rd b/man/plot_projection_uncertainty.Rd
index 86f7717..163c8f4 100644
--- a/man/plot_projection_uncertainty.Rd
+++ b/man/plot_projection_uncertainty.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/plotting.R
 \name{plot_projection_uncertainty}
 \alias{plot_projection_uncertainty}
-\title{Plot how far a monthly projection can be trusted}
+\title{Draw the layers that say how far to trust a projection}
 \usage{
 plot_projection_uncertainty(predicted, year, month, species, path)
 }
@@ -21,23 +21,29 @@ plot_projection_uncertainty(predicted, year, month, species, path)
 \code{path} invisibly, or \code{NULL} when there is nothing to draw
 }
 \description{
-The companion to \code{\link[=plot_projection]{plot_projection()}}: the same cells, panelled by what is
-uncertain about them rather than by what is predicted. Whichever of the
-uncertainty surfaces the run produced are drawn, and nothing else.
+Two or three panels above the map, one per question the projection can be
+doubted on. Only the ones a run actually produced are drawn.
 }
 \details{
-The two panels are deliberately not merged into a single "confidence" layer.
+The panels are deliberately not merged into a single "confidence" layer.
 They measure different things and are free to disagree in either direction,
 and the disagreement is the informative part: a cell can be stable across
 every member and still be extrapolated, because agreement between members
-trained on the same data is not evidence about ground the data never covered.
-Blending the two would average that case away instead of showing it.
+trained on the same data is not evidence about ground the data never
+covered. Blending them would average that case away instead of showing it.
+}
+\section{Drawn by fancyfx}{
 
-Diverging colour on the novelty panel, centred at zero, because zero is the
-meaningful break: above it the model is interpolating, below it the cell is
-outside the training range on some predictor and the model has no evidence
-for what it says there.
+The panels themselves come from \code{\link[fancyfx:plotUncertainty]{fancyfx::plotUncertainty()}} and
+\code{\link[fancyfx:plotExtrapolation]{fancyfx::plotExtrapolation()}} rather than being built here. That is worth a
+sentence because it used to be ninety lines of \code{ggplot2} in this file, and
+the replacement is better in two ways this package would otherwise have had
+to write for itself: cells are downsampled above \code{max.cells}, which a real
+Copernicus grid needs and the hand-built version did not do, and the
+extrapolation ramp diverges about zero in colours that survive the common
+colour vision deficiencies.
 }
+
 \seealso{
 \code{\link[=novelty_surface]{novelty_surface()}} for how the novelty panel is computed
 }
diff --git a/man/projection_cells.Rd b/man/projection_cells.Rd
new file mode 100644
index 0000000..710269a
--- /dev/null
+++ b/man/projection_cells.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{projection_cells}
+\alias{projection_cells}
+\title{The cells of a projection, checked}
+\usage{
+projection_cells(projection, label)
+}
+\arguments{
+\item{projection}{a projection table}
+
+\item{label}{which argument it came from, for the error}
+}
+\value{
+a data frame of \code{lon}, \code{lat} and \code{suitability}
+}
+\description{
+The cells of a projection, checked
+}
+\keyword{internal}
diff --git a/man/projection_overlap.Rd b/man/projection_overlap.Rd
new file mode 100644
index 0000000..5e587bf
--- /dev/null
+++ b/man/projection_overlap.Rd
@@ -0,0 +1,84 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{projection_overlap}
+\alias{projection_overlap}
+\title{How much two projected surfaces agree about the habitat}
+\usage{
+projection_overlap(x, y, statistic = c("both", "D", "I"), digits = 6)
+}
+\arguments{
+\item{x, y}{projections to compare: data frames with \code{lon}, \code{lat} and
+\code{suitability}, as \code{\link[=project_patch_model]{project_patch_model()}} writes to \code{suitability.csv}}
+
+\item{statistic}{\code{"both"}, \code{"D"}, or \code{"I"}}
+
+\item{digits}{how many decimal places the coordinates are matched on.
+Projections from the same grid agree exactly; this is for two that came
+off slightly different pipelines}
+}
+\value{
+a named numeric vector, with an \code{n_cells} attribute
+}
+\description{
+\code{\link[=compare_runs]{compare_runs()}} asks which run scores better. This asks something the
+scores cannot: whether the two maps are the \emph{same map}. Two runs can rank
+stations equally well and disagree completely about where the habitat is,
+and a comparison of AUCs would report them as equivalent.
+}
+\details{
+Schoener's \emph{D} and Warren's \emph{I}, both on 0 to 1. \strong{1} is identical
+surfaces; \strong{0} is no overlap at all. There is no threshold at which two
+niches become "the same" — the number is a description, and what counts as
+close depends on what the maps are for.
+}
+\section{Read them as high numbers}{
+
+Both statistics compare the suitability \emph{distributions} cell by cell after
+normalising them, so they run high. Two surfaces of independent random noise
+over the same grid score around 0.7 on \emph{D} and 0.9 on \emph{I} — not because they
+agree about anything, but because both spread their probability over the
+same cells in similar proportions. Treat those as the floor for surfaces of
+this shape rather than reading 0.9 as near-identity, and compare overlaps
+against each other rather than against 1.
+}
+
+\section{What it is good for}{
+
+Two questions, mostly. \emph{Do two species occupy the same habitat?} — run
+\code{cfin} and \code{ctyp} over the same months and compare their surfaces. And \emph{does
+the choice of algorithm change the map?} — the more useful one when a
+comparison of scores came back inconclusive, because a small difference in
+AUC with a low overlap means the two models are doing genuinely different
+things and the tie in performance is hiding it.
+}
+
+\section{Matching the cells}{
+
+The two surfaces must be the same cells in the same order, or the overlap is
+measured between places that have nothing to do with each other. Cells are
+matched on their coordinates rather than on row position, and cells present
+in one projection and not the other are dropped with a warning — a
+projection that covers less ground usually has a reason, and averaging over
+the difference would hide it.
+}
+
+\examples{
+\dontrun{
+rf <- readr::read_csv("output/rf/projections/suitability.csv")
+gam <- readr::read_csv("output/gam/projections/suitability.csv")
+
+projection_overlap(rf[rf$month == 7, ], gam[gam$month == 7, ])
+}
+}
+\references{
+Schoener TW (1968). The \emph{Anolis} lizards of Bimini: resource partitioning in
+a complex fauna. \emph{Ecology} \strong{49}(4), 704-726. \doi{10.2307/1935534} — \emph{D}
+
+Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus
+conservatism: quantitative approaches to niche evolution. \emph{Evolution}
+\strong{62}(11), 2868-2883. \doi{10.1111/j.1558-5646.2008.00482.x} — \emph{I}
+}
+\seealso{
+\code{\link[=compare_runs]{compare_runs()}}, which asks which surface is \emph{better} rather than
+whether they are the same
+}
diff --git a/man/projection_raster.Rd b/man/projection_raster.Rd
new file mode 100644
index 0000000..24d05e2
--- /dev/null
+++ b/man/projection_raster.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plotting.R
+\name{projection_raster}
+\alias{projection_raster}
+\title{One column of a projection, as a raster fancyfx can draw}
+\usage{
+projection_raster(predicted, column)
+}
+\arguments{
+\item{predicted}{a projection from \code{predict_grid()}}
+
+\item{column}{which column to rasterize}
+}
+\value{
+a single-layer \code{SpatRaster}
+}
+\description{
+The projection is a table of cells at this point in the pipeline, and every
+fancyfx map takes a \code{SpatRaster}. The cells are already on a regular grid —
+\code{project_patch_model()} builds the GeoTIFF from the same frame — so this is
+a change of container rather than any resampling.
+}
+\keyword{internal}
diff --git a/man/spatial_bias.Rd b/man/spatial_bias.Rd
new file mode 100644
index 0000000..0c6224c
--- /dev/null
+++ b/man/spatial_bias.Rd
@@ -0,0 +1,70 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spatial_bias.R
+\name{spatial_bias}
+\alias{spatial_bias}
+\title{How much of the AUC is an artefact of where the stations are}
+\usage{
+spatial_bias(model, geo = TRUE)
+}
+\arguments{
+\item{model}{a fitted model from \code{\link[=fit_patch_model]{fit_patch_model()}} or
+\code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}}, carrying \code{coordinates} and held-out \code{predictions}}
+
+\item{geo}{whether to measure great-circle distances. \code{TRUE} by default,
+because the coordinates are longitude and latitude, and a degree of
+longitude is not a degree of latitude anywhere but the equator}
+}
+\value{
+a data frame with one row per fold and a \code{fold} of \code{"overall"} for
+the mean, carrying \code{patch_distance}, \code{non_patch_distance} and \code{ssb}
+}
+\description{
+A cross-validated AUC is optimistic when the held-out patches sit closer to
+the training patches than the held-out non-patches do. The model can then
+score well by recognising \emph{where} it has already been rather than \emph{what}
+makes a patch, and nothing in the ROC curve distinguishes the two.
+}
+\details{
+This is the spatial sorting bias of Hijmans (2012): the mean distance from
+each held-out patch to the nearest training patch, over the same distance
+for the held-out non-patches.
+
+\deqn{SSB = \frac{\overline{d}(\mathrm{held\ out\ patch},\ \mathrm{training\ patch})}{\overline{d}(\mathrm{held\ out\ non\ patch},\ \mathrm{training\ patch})}}
+\itemize{
+\item \strong{Near 1} — held-out patches and non-patches are equally far from the
+training patches. The split is spatially fair and the AUC means what it
+appears to.
+\item \strong{Near 0} — held-out patches are much closer to training patches than the
+non-patches are. \strong{The AUC is inflated}, and by an amount this number
+does not tell you.
+}
+
+It is a property of the \emph{split}, not of the model. Every model type fitted on
+the same folds gets the same value, which is the point: it says how much to
+trust the comparison between them, not which of them won.
+}
+\section{Why an ECOMON run should look at it}{
+
+Survey stations are not scattered at random. They sit on transects, revisited
+season after season, so a randomly held-out station usually has a near
+neighbour in the training folds — often the same station in another year.
+That is the situation this measures, and random \code{model.cv_folds} cannot avoid
+it. A low value is an argument for spatially blocked folds, not for a
+different model.
+}
+
+\examples{
+\dontrun{
+model <- fit_patch_model(dat, config)
+spatial_bias(model)
+}
+}
+\references{
+Hijmans RJ (2012). Cross-validation of species distribution models: removing
+spatial sorting bias and calibration with a null model. \emph{Ecology} \strong{93}(3),
+679-688. \doi{10.1890/11-0826.1}
+}
+\seealso{
+\code{\link[=evaluation_table]{evaluation_table()}}, which reports the summary beside the AUC it
+qualifies
+}
diff --git a/man/spatial_bias_note.Rd b/man/spatial_bias_note.Rd
new file mode 100644
index 0000000..bd90bd4
--- /dev/null
+++ b/man/spatial_bias_note.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spatial_bias.R
+\name{spatial_bias_note}
+\alias{spatial_bias_note}
+\title{The one-line reading of a spatial sorting bias}
+\usage{
+spatial_bias_note(ssb)
+}
+\arguments{
+\item{ssb}{the overall spatial sorting bias}
+}
+\value{
+a single string
+}
+\description{
+Written for the run log, where a bare ratio would be ignored. The thresholds
+are conventional rather than derived — Hijmans (2012) offers no cutoff, and
+presenting one as though he did would be worse than a rule of thumb labelled
+as one.
+}
+\keyword{internal}
diff --git a/man/station_coordinates.Rd b/man/station_coordinates.Rd
new file mode 100644
index 0000000..7008018
--- /dev/null
+++ b/man/station_coordinates.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/model.R
+\name{station_coordinates}
+\alias{station_coordinates}
+\title{Station coordinates, in the order the modelling data is in}
+\usage{
+station_coordinates(dat)
+}
+\arguments{
+\item{dat}{labeled modeling data}
+}
+\value{
+a two-column data frame of \code{lon} and \code{lat}, or \code{NULL}
+}
+\description{
+\code{NULL} when the station table has no coordinates, which is the case for a
+hand-built frame in a test rather than for anything a run produces.
+}
+\keyword{internal}
diff --git a/man/variable_similarity.Rd b/man/variable_similarity.Rd
deleted file mode 100644
index 6728270..0000000
--- a/man/variable_similarity.Rd
+++ /dev/null
@@ -1,22 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/uncertainty.R
-\name{variable_similarity}
-\alias{variable_similarity}
-\title{Similarity of values to a training distribution, for one predictor}
-\usage{
-variable_similarity(values, train)
-}
-\arguments{
-\item{values}{the values to score}
-
-\item{train}{the training values for the same predictor}
-}
-\value{
-a numeric vector the length of \code{values}
-}
-\description{
-The per-variable half of \code{\link[=novelty_surface]{novelty_surface()}}. Negative below the training
-minimum and above its maximum, scaled by the training range; inside, twice
-the distance to the nearer tail in percentile terms, so the median scores 100.
-}
-\keyword{internal}
diff --git a/man/with_dependencies.Rd b/man/with_dependencies.Rd
new file mode 100644
index 0000000..873a98d
--- /dev/null
+++ b/man/with_dependencies.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/derivoce.R
+\name{with_dependencies}
+\alias{with_dependencies}
+\title{The chosen derived covariates, plus the ones they are computed from}
+\usage{
+with_dependencies(ids, choices)
+}
+\arguments{
+\item{ids}{chosen derived covariate ids}
+
+\item{choices}{the \code{\link[=derivoce_choices]{derivoce_choices()}} list}
+}
+\value{
+\code{ids} with any dependencies added
+}
+\description{
+A gradient of current speed needs the speed, and the speed is itself a
+derived covariate rather than a download. Choosing the gradient therefore
+has to pull in the step that produces what it reads — otherwise the config
+asks derivoce for a gradient of a column that was never computed.
+}
+\details{
+Order is not this function's problem. \code{derivoce_choices()} lists a
+dependency before anything that depends on it, and steps are emitted in that
+order, so \code{current_speed} runs before the gradient of \code{speed} without
+anything having to sort them.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-derivoce.R b/tests/testthat/test-derivoce.R
index 183bf2a..6d665d6 100644
--- a/tests/testthat/test-derivoce.R
+++ b/tests/testthat/test-derivoce.R
@@ -274,3 +274,66 @@ test_that("a run with no derivoce block is untouched", {
   env <- fetch_covariates(config, years = 2018, months = 6)
   expect_identical(add_derivoce_covariates(env, config), env)
 })
+
+test_that("a derived covariate is offered whether or not its source is modelled", {
+  # Integrated chlorophyll without chlorophyll itself is an ordinary thing to
+  # want - the accumulated bloom is what feeds the animals, the instantaneous
+  # value is not - and it used to be unaskable: the derived form appeared only
+  # once its source was a predictor, and then the source could not be removed.
+  ids <- vapply(derivoce_choices(c("SST", "SSS")), function(x) x$id,
+                character(1))
+
+  expect_true("CHL_int" %in% ids)
+  expect_true("CHL_lag1" %in% ids)
+  expect_true("SST_int" %in% ids)
+  # And the source is fetched, so it can be excluded from the predictors.
+  expect_equal(derivoce_required_inputs("CHL_int", c("SST", "SSS")), "CHL")
+  # A source that *is* selected needs no extra download.
+  expect_length(derivoce_required_inputs("SST_int", c("SST", "SSS")), 0)
+})
+
+test_that("the extra choices say in their group that they cost a download", {
+  choices <- derivoce_choices(c("SST"))
+  group_of <- function(id) {
+    Filter(function(x) x$id == id, choices)[[1]]$group
+  }
+
+  expect_equal(group_of("SST_int"), "Temporal")
+  expect_match(group_of("CHL_int"), "downloads CHL")
+})
+
+test_that("current speed can be differentiated, not just its components", {
+  # The gradient of the speed is the quantity a front in the flow is. The
+  # components can each be changing steeply while the speed is constant, which
+  # is a turn rather than a shear.
+  ids <- vapply(derivoce_choices(c("SST")), function(x) x$id, character(1))
+  expect_true("speed_grad" %in% ids)
+  expect_true("EKE_grad" %in% ids)
+
+  steps <- derivoce_steps_for("speed_grad", c("SST"))
+
+  # The speed has to be computed before it can be differentiated, and the
+  # config is read in order.
+  expect_equal(vapply(steps, function(s) s$type, character(1)),
+               c("current_speed", "horizontal_gradient"))
+  expect_equal(steps[[2]]$vars, "speed")
+  expect_setequal(derivoce_required_inputs("speed_grad", c("SST")),
+                  c("UO", "VO"))
+})
+
+test_that("a dependency pulled in on the way is not made a predictor", {
+  expect_equal(derivoce_dependency_columns("speed_grad", c("SST")), "speed")
+  # Unless it was asked for in its own right.
+  expect_length(derivoce_dependency_columns(c("speed", "speed_grad"), c("SST")),
+                0)
+  # A download is not a derived column; derivoce_required_inputs covers those.
+  expect_length(derivoce_dependency_columns("CHL_int", c("SST")), 0)
+  expect_length(derivoce_dependency_columns(character(), c("SST")), 0)
+})
+
+test_that("asking for a derived covariate twice over does not duplicate a step", {
+  steps <- derivoce_steps_for(c("speed", "speed_grad"), c("SST"))
+
+  expect_equal(sum(vapply(steps, function(s) s$type, character(1)) ==
+                     "current_speed"), 1)
+})
diff --git a/tests/testthat/test-effects.R b/tests/testthat/test-effects.R
index 8df1153..ee3342e 100644
--- a/tests/testthat/test-effects.R
+++ b/tests/testthat/test-effects.R
@@ -176,7 +176,6 @@ test_that("the smoothed variables are read back off the fitted model", {
 test_that("fancyfx draws the fitted smooths for a GAM", {
   skip_on_cran()
   skip_if_not_installed("mgcv")
-  skip_if_not_installed("fancyfx")
   model <- fitted_of_type("gam")
 
   plot <- suppressMessages(plot_gam_smooths(model))
@@ -194,7 +193,6 @@ test_that("fancyfx draws the fitted smooths for a GAM", {
 test_that("smooth plots are drawn against the data the model actually saw", {
   skip_on_cran()
   skip_if_not_installed("mgcv")
-  skip_if_not_installed("fancyfx")
   model <- fitted_of_type("gam")
 
   # gratia reports the smooths in the recipe's output units, so the rug has to
@@ -207,21 +205,6 @@ test_that("smooth plots are drawn against the data the model actually saw", {
   expect_gt(abs(mean(model$model_data$SST)), 1)
 })
 
-test_that("without fancyfx the run still gets its generic curves", {
-  skip_on_cran()
-  skip_if_not_installed("mgcv")
-  # fancyfx is a Suggests: its absence removes the extra plot, not the
-  # diagnostics.
-  local_mocked_bindings(has_fancyfx = function() FALSE)
-  out <- tempfile("diag"); dir.create(out)
-  suppressMessages(suppressWarnings(write_effect_plots(fitted_of_type("gam"), out)))
-  written <- list.files(out)
-
-  expect_false("gam_smooths.png" %in% written)
-  expect_true("partial_effects.png" %in% written)
-  expect_true("smooth_terms.csv" %in% written)
-})
-
 test_that("maps draw, not just build", {
   skip_on_cran()
   # geom_sf() refuses to work under any coord but coord_sf(), and refuses when
@@ -266,3 +249,41 @@ test_that("abundance over the record is one continuous series", {
   expect_false("GeomLine" %in% geoms)
   expect_true(all(geoms %in% c("GeomPoint", "GeomJitter")))
 })
+
+test_that("the uncertainty panels are drawn for whichever surfaces exist", {
+  skip_on_cran()
+  cells <- expand.grid(lon = seq(-70, -66, by = 0.5),
+                       lat = seq(41, 44, by = 0.5))
+  set.seed(1)
+  cells$suitability_sd <- runif(nrow(cells), 0, 0.1)
+  cells$novelty <- runif(nrow(cells), -20, 90)
+  cells$novel_variable <- "SST"
+
+  path <- tempfile(fileext = ".png")
+  expect_equal(plot_projection_uncertainty(cells, 2018, 6, "cfin", path), path)
+  expect_true(file.exists(path))
+
+  # An algorithm ensemble adds a third panel; a run without one does not get it.
+  cells$algorithm_sd <- runif(nrow(cells), 0, 0.2)
+  three <- tempfile(fileext = ".png")
+  plot_projection_uncertainty(cells, 2018, 6, "cfin", three)
+  expect_gt(file.info(three)$size, file.info(path)$size)
+})
+
+test_that("a projection with no uncertainty surfaces draws nothing", {
+  cells <- data.frame(lon = c(-70, -69), lat = c(41, 42),
+                      suitability = c(0.2, 0.8))
+
+  expect_null(plot_projection_uncertainty(cells, 2018, 6, "cfin",
+                                          tempfile(fileext = ".png")))
+})
+
+test_that("the novelty subtitle counts the cells and names the culprit", {
+  clean <- data.frame(novelty = c(10, 40, 90))
+  expect_match(novelty_subtitle(clean), "Every cell is inside")
+
+  mixed <- data.frame(novelty = c(-5, 40, -2, 90),
+                      novel_variable = c("CHL", NA, "CHL", NA))
+  expect_match(novelty_subtitle(mixed), "2 of 4 cells")
+  expect_match(novelty_subtitle(mixed), "mostly CHL")
+})
diff --git a/tests/testthat/test-power.R b/tests/testthat/test-power.R
index 3e32779..c6abd3c 100644
--- a/tests/testthat/test-power.R
+++ b/tests/testthat/test-power.R
@@ -211,3 +211,77 @@ test_that("a power curve is reproducible", {
   expect_equal(once$difference, twice$difference)
   expect_equal(once$std_err, twice$std_err)
 })
+
+# ---- projection overlap ------------------------------------------------------
+
+overlap_cells <- function(seed = 1) {
+  cells <- expand.grid(lon = seq(-70, -66, by = 0.5),
+                       lat = seq(41, 44, by = 0.5))
+  set.seed(seed)
+  cells$suitability <- stats::runif(nrow(cells))
+  cells
+}
+
+test_that("a surface overlaps itself completely", {
+  cells <- overlap_cells()
+
+  out <- projection_overlap(cells, cells)
+
+  expect_equal(unname(out[["D"]]), 1)
+  expect_equal(unname(out[["I"]]), 1)
+  expect_equal(attr(out, "n_cells"), nrow(cells))
+})
+
+test_that("a surface overlaps its own inverse less than itself", {
+  cells <- overlap_cells()
+  flipped <- cells
+  flipped$suitability <- 1 - cells$suitability
+
+  expect_lt(projection_overlap(cells, flipped)[["D"]],
+            projection_overlap(cells, cells)[["D"]])
+})
+
+test_that("one statistic can be asked for on its own", {
+  cells <- overlap_cells()
+
+  expect_named(projection_overlap(cells, cells, statistic = "D"), "D")
+  expect_named(projection_overlap(cells, cells, statistic = "I"), "I")
+  expect_setequal(names(projection_overlap(cells, cells)), c("D", "I"))
+})
+
+test_that("cells are matched on coordinates rather than on row order", {
+  # Two projections written by different runs need not be in the same order,
+  # and comparing row i of one with row i of the other would be comparing two
+  # different places.
+  cells <- overlap_cells()
+  shuffled <- cells[order(stats::runif(nrow(cells))), ]
+
+  expect_equal(unname(projection_overlap(cells, shuffled)[["D"]]), 1)
+})
+
+test_that("projections covering different ground are intersected, loudly", {
+  cells <- overlap_cells()
+  partial <- cells[cells$lat <= 43, ]
+
+  expect_warning(out <- projection_overlap(cells, partial),
+                 "do not cover the same cells")
+  expect_equal(attr(out, "n_cells"), nrow(partial))
+})
+
+test_that("projections with no common ground are refused", {
+  cells <- overlap_cells()
+  elsewhere <- cells
+  elsewhere$lon <- elsewhere$lon + 100
+
+  expect_error(suppressWarnings(projection_overlap(cells, elsewhere)),
+               "share 0 cells")
+})
+
+test_that("a table that is not a projection says what it is missing", {
+  cells <- overlap_cells()
+
+  expect_error(projection_overlap(cells[, c("lon", "lat")], cells),
+               "missing: suitability")
+  expect_error(projection_overlap(cells, cells[, c("lon", "suitability")]),
+               "missing: lat")
+})
diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R
index 19cb0f8..5314c68 100644
--- a/tests/testthat/test-uncertainty.R
+++ b/tests/testthat/test-uncertainty.R
@@ -42,26 +42,22 @@ test_that("malformed uncertainty settings are refused at config load", {
 
 # ---- novelty ----------------------------------------------------------------
 
-test_that("similarity is negative outside the training range and scaled by it", {
-  train <- c(0, 10)   # range of 10
+test_that("the MESS scale still reads the way the documentation says", {
+  # The arithmetic belongs to fancyfx::mess() now, and is tested there. What is
+  # checked here is the property this package documents and its readers rely
+  # on: 100 at the median, falling to 0 at the edge of the training range, and
+  # negative outside it in proportion to how far.
+  train <- data.frame(x = 0:100)
 
-  # Half a range below the minimum, and one range above the maximum.
-  expect_equal(unname(variable_similarity(-5, train)), -50)
-  expect_equal(unname(variable_similarity(20, train)), -100)
+  middle <- novelty_surface(data.frame(x = 50), train, "x")$novelty
+  edge <- novelty_surface(data.frame(x = 99), train, "x")$novelty
+  outside <- novelty_surface(data.frame(x = c(-50, 150)), train, "x")$novelty
 
-  # Inside, it is positive. The endpoints are the edge of the range, not
-  # outside it, so they are not negative.
-  expect_gte(variable_similarity(5, train), 0)
-})
-
-test_that("similarity peaks at the middle of the training data", {
-  train <- 1:101
-
-  middle <- variable_similarity(51, train)
-  edge <- variable_similarity(95, train)
-
-  expect_gt(middle, edge)
   expect_lte(middle, 100)
+  expect_gt(middle, edge)
+  expect_gte(edge, 0)
+  # Half a training range below the minimum, and half a range above the max.
+  expect_equal(outside, c(-50, -50))
 })
 
 test_that("a cell is as novel as its worst predictor, and says which", {
@@ -173,3 +169,88 @@ test_that("the extra surfaces reach the GeoTIFF as named layers", {
   expect_true(any(grepl("_uncertainty\\.png$", plots)))
   expect_true(any(!grepl("_uncertainty\\.png$", plots)))
 })
+
+# ---- spatial sorting bias ---------------------------------------------------
+
+test_that("spatially fair folds score near 1 and sorted ones score near 0", {
+  # Built rather than fitted, so the answer is known in advance. Patches and
+  # non-patches drawn from the same places is the fair case; patches clustered
+  # away from the non-patches is the sorted one.
+  set.seed(1)
+  n <- 200
+  fair <- data.frame(lon = runif(n, -70, -66), lat = runif(n, 41, 44))
+  # Shuffled rather than alternating: `rep` of a 2-cycle across a 4-fold cycle
+  # puts every patch in the odd folds, leaving each fold with one class and
+  # nothing to compare.
+  labels <- sample(rep(c("patch", "non_patch"), length.out = n))
+  predictions <- data.frame(
+    .row = seq_len(n),
+    id = rep(paste0("Fold", 1:4), length.out = n),
+    patch = factor(labels, levels = c("patch", "non_patch"))
+  )
+
+  unbiased <- spatial_bias(list(coordinates = fair, predictions = predictions))
+  expect_equal(overall_ssb(unbiased), 1, tolerance = 0.35)
+
+  # Now put every patch in one corner and every non-patch in another.
+  sorted <- fair
+  is_patch <- predictions$patch == "patch"
+  sorted$lon[is_patch] <- runif(sum(is_patch), -70, -69.5)
+  sorted$lon[!is_patch] <- runif(sum(!is_patch), -66.5, -66)
+
+  biased <- spatial_bias(list(coordinates = sorted, predictions = predictions))
+  expect_lt(overall_ssb(biased), overall_ssb(unbiased))
+  expect_lt(overall_ssb(biased), 0.3)
+})
+
+test_that("spatial bias reports every fold as well as the overall", {
+  set.seed(2)
+  n <- 120
+  coordinates <- data.frame(lon = runif(n, -70, -66), lat = runif(n, 41, 44))
+  predictions <- data.frame(
+    .row = seq_len(n),
+    id = rep(paste0("Fold", 1:3), length.out = n),
+    patch = factor(sample(rep(c("patch", "non_patch"), length.out = n)),
+                   levels = c("patch", "non_patch"))
+  )
+
+  out <- spatial_bias(list(coordinates = coordinates, predictions = predictions))
+
+  expect_equal(nrow(out), 4)
+  expect_equal(out$fold, c("Fold1", "Fold2", "Fold3", "overall"))
+  expect_true(all(out$ssb > 0))
+})
+
+test_that("spatial bias declines to answer without coordinates", {
+  # A model fitted on a hand-built frame has none, and that is not an error.
+  predictions <- data.frame(.row = 1:4, id = "Fold1",
+                            patch = factor(c("patch", "non_patch"),
+                                           levels = c("patch", "non_patch")))
+
+  expect_null(spatial_bias(list(coordinates = NULL, predictions = predictions)))
+  expect_null(spatial_bias(list(coordinates = data.frame(lon = 1, lat = 1),
+                                predictions = NULL)))
+  expect_true(is.na(overall_ssb(NULL)))
+})
+
+test_that("the spatial bias note says what a low value means for the AUC", {
+  expect_match(spatial_bias_note(0.95), "spatially fair")
+  expect_match(spatial_bias_note(0.2), "optimistic")
+  expect_match(spatial_bias_note(NA_real_), "could not be computed")
+})
+
+test_that("the fitted model carries its coordinates and its spatial bias", {
+  skip_on_cran()
+  config <- mock_config()
+  config$model$trees <- 50
+  config$model$cv_folds <- 5
+  dat <- labeled_mock_data(config)
+
+  model <- fit_patch_model(dat, config)
+
+  expect_equal(nrow(model$coordinates), nrow(model$model_data))
+  expect_setequal(names(model$coordinates), c("lon", "lat"))
+  expect_false(is.null(model$spatial_bias))
+  # And it reaches the evaluation table, beside the metrics it qualifies.
+  expect_true("ssb" %in% model$evaluation$metric)
+})
diff --git a/tools/citations.csv b/tools/citations.csv
index 3ece160..c9c1d7d 100644
--- a/tools/citations.csv
+++ b/tools/citations.csv
@@ -41,3 +41,6 @@ elith2011,10.1111/j.1472-4642.2010.00725.x,,Elith,2011,A statistical explanation
 nadeau2003,10.1023/A:1024068626366,,Nadeau,2003,Inference for the Generalization Error,Machine Learning,52,239-281,crossref,
 dietterich1998,10.1162/089976698300017197,,Dietterich,1998,Approximate Statistical Tests for Comparing Supervised Classification Learning Algorithms,Neural Computation,10,1895-1923,crossref,
 hoenig2001,10.1198/000313001300339897,,Hoenig,2001,The Abuse of Power,The American Statistician,55,19-24,crossref,
+hijmans2012,10.1890/11-0826.1,,Hijmans,2012,Cross‐validation of species distribution models: removing spatial sorting bias and calibration with a null model,Ecology,93,679-688,crossref,
+schoener1968,10.2307/1935534,,Schoener,1968,The Anolis Lizards of Bimini: Resource Partitioning in a Complex Fauna,Ecology,49,704-726,crossref,
+warren2008,10.1111/j.1558-5646.2008.00482.x,,Warren,2008,ENVIRONMENTAL NICHE EQUIVALENCY VERSUS CONSERVATISM: QUANTITATIVE APPROACHES TO NICHE EVOLUTION,Evolution,62,2868-2883,crossref,