From 54239b3e3602a0a7f09aa522692a376454416e57 Mon Sep 17 00:00:00 2001
From: Brian Maitner <bmaitner@usf.edu>
Date: Tue, 22 Sep 2026 10:43:57 -0400
Subject: [PATCH] working on local version

---
 DESCRIPTION                     |   8 +-
 R/GVS_local.R                   | 659 ++++++++++++++++++++++++++++++++
 R/local_build.R                 | 349 +++++++++++++++++
 R/local_cache.R                 | 287 ++++++++++++++
 R/local_index.R                 | 239 ++++++++++++
 tests/testthat/test-local-gvs.R | 143 +++++++
 6 files changed, 1684 insertions(+), 1 deletion(-)
 create mode 100644 R/GVS_local.R
 create mode 100644 R/local_build.R
 create mode 100644 R/local_cache.R
 create mode 100644 R/local_index.R
 create mode 100644 tests/testthat/test-local-gvs.R

diff --git a/DESCRIPTION b/DESCRIPTION
index 4d410bb..56c3d67 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -26,5 +26,11 @@ Suggests:
         testthat,
         devtools,
         BIEN,
-        vcr (>= 0.6.0)
+        vcr (>= 0.6.0),
+        sf,
+        terra,
+        nanoparquet,
+        geosphere,
+        countrycode,
+        GNRS
 VignetteBuilder: knitr
diff --git a/R/GVS_local.R b/R/GVS_local.R
new file mode 100644
index 0000000..36a2315
--- /dev/null
+++ b/R/GVS_local.R
@@ -0,0 +1,659 @@
+#' Validate coordinates without an internet connection, with historical divisions
+#'
+#' Offline Geocoordinate Validation Service.  For each coordinate: the GADM
+#' country, state/province and county/parish containing it; its distance to the
+#' centroids of those divisions, and whether it is likely a centroid, as the
+#' service computes them; the same centroid test against the historical countries
+#' containing it (CShapes 2.0, 1886-2019), so that a record placed at the USSR's
+#' centroid is detected; and, when a GNRS resolution is supplied, whether the point
+#' lies in the political division the record names, now or at any time.
+#'
+#' @param occurrence_dataframe A data.frame with numeric \code{latitude} and
+#'   \code{longitude} (WGS84), and optionally \code{date} (years or ISO dates),
+#'   used when \code{history = "at_date"}.
+#' @param history How historical countries count towards geovalidity.
+#'   \code{"all"} (default): a point is geovalid if it lies in the named division
+#'   as it was at any time (1886 onwards).  \code{"at_date"}: only versions valid
+#'   at the record's date (± \code{tolerance_years}) count; records without a date
+#'   are treated as under \code{"all"}.  Centroid detection always uses every
+#'   historical version, since a centroid copied from an old gazetteer can appear
+#'   in a record of any date.
+#' @param gnrs Optional output of \code{GNRS::GNRS_local()} for the same rows, in
+#'   the same order (with \code{history = "all"} or \code{"at_date"}, so it carries
+#'   \code{entity_key}).  Needed for the geovalidity columns.
+#' @param dir Cache directory, shared with GNRS.
+#' @param gpkg Optional path to the GADM GeoPackage the index was built from.
+#'   With it, points in boundary cells are located exactly; without it they keep
+#'   the raster answer and are flagged in \code{locate_method}.
+#' @param thr_abs,thr_rel Centroid thresholds: named numeric vectors with any of
+#'   \code{country}, \code{state}, \code{county}, absolute in km and relative
+#'   (distance over the division's maximum distance from that centroid).  Either may
+#'   be \code{NULL}, to test only the other kind of distance; a level missing from a
+#'   vector is not tested with that kind of distance.  At least one must be given.
+#'   The defaults are \code{GVS_thresholds("spec995")}: a point is a centroid if it is
+#'   within both the absolute and the relative threshold of its country's or county's
+#'   centroid.  On BIEN's own records this catches 66\% of the occurrences BIEN declares
+#'   centroids while flagging 1.6\% of all occurrences, against 39\% and 0.53\% for the
+#'   service's historical rule (relative 0.002).  \code{\link{GVS_thresholds}} holds the
+#'   other sets, including the TSS-optimised ones the GVS manuscript reports.
+#' @param combine When a level has both an absolute and a relative threshold:
+#'   \code{"or"} flags a point meeting either, \code{"and"} only a point meeting
+#'   both.  Levels are always combined with OR.
+#' @param max_uncertainty_m Coordinates whose inherent uncertainty (from the decimal
+#'   places submitted) is at least this are flagged \code{is_low_precision}; \code{NULL}
+#'   switches the flag off.  The default, 11 km, is one decimal place or fewer.  This is
+#'   a separate flag from \code{is_centroid}: a rounded coordinate is imprecise, which is
+#'   not the same claim as sitting on a division's centroid, though in practice gazetteer
+#'   points are both.  On the GVS benchmark it catches 74\% of the state-level records
+#'   (rounded gazetteer points a median 23 km from any state centroid, which centroid
+#'   distance cannot catch) at specificity 0.998.
+#' @param tolerance_years Tolerance either side of a historical version's dates
+#'   when \code{date} is given.
+#' @return A data.frame, one row per input row.
+#' @export
+GVS_local <- function(occurrence_dataframe, gnrs = NULL, dir = gvs_cache_dir(), gpkg = NULL,
+                      thr_abs = c(country = 38.466, county = 0.24727),
+                      thr_rel = c(country = 0.026656, county = 0.0025027),
+                      combine = c("and", "or"), max_uncertainty_m = 11000,
+                      history = c("all", "at_date"), tolerance_years = 1) {
+  history <- match.arg(history)
+  combine <- match.arg(combine)
+  gvs_check_thresholds(thr_abs, thr_rel)
+  if (!inherits(occurrence_dataframe, "data.frame") ||
+      !all(c("latitude", "longitude") %in% names(occurrence_dataframe))) {
+    stop("occurrence_dataframe should be a data.frame with latitude and longitude columns", call. = FALSE)
+  }
+  if (!is.null(gnrs) && nrow(gnrs) != nrow(occurrence_dataframe)) {
+    stop("gnrs should have one row per row of occurrence_dataframe", call. = FALSE)
+  }
+  dp_lat <- gvs_decimal_places(occurrence_dataframe$latitude)
+  dp_lon <- gvs_decimal_places(occurrence_dataframe$longitude)
+  lon <- suppressWarnings(as.numeric(occurrence_dataframe$longitude))
+  lat <- suppressWarnings(as.numeric(occurrence_dataframe$latitude))
+  n <- length(lon)
+  dates <- if (history == "at_date" && "date" %in% names(occurrence_dataframe)) {
+    gvs_parse_record_date(occurrence_dataframe$date)
+  } else NULL
+
+  # ---- 1. current GADM divisions -------------------------------------------------
+  loc <- gvs_locate_current(lon, lat, dir = dir, gpkg = gpkg)
+
+  # ---- 2. centroids of the current divisions (the service's computation) ---------
+  ct <- lapply(stats::setNames(gvs_levels(), gvs_levels()), function(l)
+    as.data.frame(nanoparquet::read_parquet(gvs_centroid_path(l, dir))))
+  thr_of <- function(thr, level) if (is.null(thr) || !level %in% names(thr)) NULL else unname(thr[[level]])
+  dc <- gvs_level_distances(lon, lat, loc$gid_0, ct$country, thr = thr_of(thr_abs, "country"))
+  ds <- gvs_level_distances(lon, lat, loc$gid_1, ct$state, thr = thr_of(thr_abs, "state"))
+  dp <- gvs_level_distances(lon, lat, loc$gid_2, ct$county, thr = thr_of(thr_abs, "county"))
+  out <- data.frame(
+    latitude = lat, longitude = lon,
+    gid_0 = loc$gid_0, country = loc$country, gid_1 = loc$gid_1, state = loc$name_1,
+    gid_2 = loc$gid_2, county = loc$name_2, locate_method = loc$locate_method,
+    country_cent_dist = dc$dist, country_cent_dist_relative = dc$rel, country_cent_type = dc$type,
+    state_cent_dist = ds$dist, state_cent_dist_relative = ds$rel, state_cent_type = ds$type,
+    county_cent_dist = dp$dist, county_cent_dist_relative = dp$rel, county_cent_type = dp$type,
+    stringsAsFactors = FALSE
+  )
+  rel_m <- cbind(dc$rel, ds$rel, dp$rel)
+  dist_m <- cbind(dc$dist, ds$dist, dp$dist)
+  type_m <- cbind(dc$type, ds$type, dp$type)
+  allna <- rowSums(!is.na(rel_m)) == 0
+  k <- max.col(-replace(rel_m, is.na(rel_m), Inf), ties.method = "first")
+  ix <- cbind(seq_len(n), k)
+  out$centroid_dist <- ifelse(allna, NA, dist_m[ix])
+  out$centroid_dist_relative <- ifelse(allna, NA, rel_m[ix])
+  out$centroid_type <- ifelse(allna, NA, type_m[ix])
+  out$centroid_poldiv <- ifelse(allna, NA, c("country", "state", "county")[k])
+  out$is_centroid <- as.integer(
+    gvs_centroid_flag(dc$dist, dc$rel, thr_of(thr_abs, "country"), thr_of(thr_rel, "country"), combine) |
+    gvs_centroid_flag(ds$dist, ds$rel, thr_of(thr_abs, "state"), thr_of(thr_rel, "state"), combine) |
+    gvs_centroid_flag(dp$dist, dp$rel, thr_of(thr_abs, "county"), thr_of(thr_rel, "county"), combine))
+
+  # ---- 2b. precision of the submitted coordinates ---------------------------------
+  out$coordinate_decimal_places <- pmin(dp_lat, dp_lon)
+  out$coordinate_inherent_uncertainty_m <- gvs_coordinate_uncertainty_m(lat, dp_lat, dp_lon)
+  out$is_low_precision <- if (is.null(max_uncertainty_m)) NA_integer_ else {
+    z <- out$coordinate_inherent_uncertainty_m >= max_uncertainty_m
+    as.integer(ifelse(is.na(z), FALSE, z))
+  }
+
+  # ---- 3. centroids of historical countries (CShapes) -----------------------------
+  hc <- gvs_historical_centroids(lon, lat, NULL, dir, thr_of(thr_abs, "country"), thr_of(thr_rel, "country"),
+                                 tolerance_years, combine)
+  out$hist_centroid_version <- hc$version_id
+  out$hist_centroid_name <- hc$country_name
+  out$hist_centroid_valid_from <- hc$valid_from
+  out$hist_centroid_valid_to <- hc$valid_to
+  out$hist_centroid_dist <- hc$dist
+  out$hist_centroid_dist_relative <- hc$rel
+  out$hist_centroid_type <- hc$type
+  out$is_centroid_historical <- as.integer(!is.na(hc$version_id))
+  out$is_centroid_any <- as.integer(out$is_centroid == 1L | out$is_centroid_historical == 1L)
+
+  # ---- 4. geovalidity against the named division ----------------------------------
+  if (!is.null(gnrs)) {
+    gv <- gvs_geovalidity(lon, lat, dates, loc, gnrs, dir, tolerance_years)
+    out <- cbind(out, gv)
+  }
+  rownames(out) <- NULL
+  out
+}
+
+#' Group identical rows of several vectors without pasting them
+#'
+#' Internal.  \code{first} indexes one representative row per group (in sort
+#' order, not order of appearance); \code{group} maps every row to its
+#' representative's position in \code{first}.  NAs compare equal to each other.
+#' @keywords internal
+#' @noRd
+gvs_dedupe <- function(...) {
+  cols <- list(...)
+  n <- length(cols[[1]])
+  if (!n) return(list(first = integer(0), group = integer(0)))
+  o <- do.call(order, c(unname(cols), list(na.last = TRUE, method = "radix")))
+  same <- rep(TRUE, n - 1L)
+  for (x in cols) {
+    a <- x[o][-n]
+    b <- x[o][-1L]
+    eq <- a == b
+    eq[is.na(eq)] <- is.na(a[is.na(eq)]) & is.na(b[is.na(eq)])
+    same <- same & eq
+  }
+  start <- c(TRUE, !same)
+  gid <- cumsum(start)
+  group <- integer(n)
+  group[o] <- gid
+  list(first = o[start], group = group)
+}
+
+#' Centroid-detection threshold sets
+#'
+#' The threshold sets fitted to the GVS benchmark of 300,000 known centroids (BIEN
+#' records georeferenced as country, state or county centroids) and 300,000
+#' iNaturalist records, described in \code{gvs_ms/R_scripts/04b_centroid_threshold_approaches.R}.
+#' Use with \code{\link{GVS_local}} as \code{thr_abs}, \code{thr_rel} and
+#' \code{combine}.
+#'
+#' \describe{
+#'   \item{spec99}{Maximum sensitivity at specificity 0.99, requiring a point to be
+#'     within both the absolute and the relative threshold of a division's centroid
+#'     (\code{combine = "and"}).  The best of the rules tried at this specificity, by
+#'     both weightings.  Held-out sensitivity per record 0.55-0.65; per distinct
+#'     location 0.79-0.81, made up of 0.94-0.95 of county centroids, 0.34-0.40 of
+#'     country centroids and almost no state centroids, which at this specificity cost
+#'     more good records than they save, so the state test is off.  (Per record, country
+#'     detection reads 0.83-0.92, but the benchmark's 100,000 country records sit on 132
+#'     distinct coordinates, so that figure reflects a few well-replicated places; the
+#'     country thresholds are fitted on those 132 locations and should be treated as
+#'     provisional.)  Flags 2.9\% of BIEN occurrences.  State-level gazetteer points are
+#'     better caught by \code{max_uncertainty_m}, not by centroid distance.}
+#'   \item{spec995}{As \code{spec99} at specificity 0.995 (held-out sensitivity
+#'     0.51-0.54).}
+#'   \item{spec999}{Specificity 0.999; county centroids only (sensitivity 0.29-0.33).}
+#'   \item{tss_abs, tss_rel}{The TSS-optimised sets of the GVS manuscript, each used
+#'     on its own with \code{combine = "or"}: specificity 0.84 and 0.80, flagging 26\%
+#'     and 22\% of BIEN occurrences.  \code{tss_abs} with \code{tss_rel} and
+#'     \code{combine = "or"} (all six at once) was never validated: specificity 0.73.}
+#'   \item{bien_default}{The service's historical rule, relative distance 0.002 at any
+#'     level (specificity 0.997, sensitivity 0.29).}
+#' }
+#' @return A list with \code{thr_abs}, \code{thr_rel} and \code{combine}, or, with no
+#'   arguments, the names of the available sets.
+#' @param set One of "spec99", "spec995", "spec999", "tss_abs", "tss_rel", "bien_default".
+#' @export
+#' @examples
+#' GVS_thresholds()
+#' GVS_thresholds("spec99")
+GVS_thresholds <- function(set = NULL) {
+  sets <- list(
+    spec99 = list(thr_abs = c(country = 32.911, county = 1.0357),
+                  thr_rel = c(country = 0.070671, county = 0.0048643), combine = "and"),
+    spec995 = list(thr_abs = c(country = 38.466, county = 0.24727),
+                   thr_rel = c(country = 0.026656, county = 0.0025027), combine = "and"),
+    spec999 = list(thr_abs = NULL, thr_rel = c(county = 0.0066129), combine = "and"),
+    tss_abs = list(thr_abs = c(country = 69.99529, state = 31.75486, county = 1.01550),
+                   thr_rel = NULL, combine = "or"),
+    tss_rel = list(thr_abs = NULL,
+                   thr_rel = c(country = 0.01606376, state = 0.15369636, county = 0.03163571),
+                   combine = "or"),
+    bien_default = list(thr_abs = NULL, thr_rel = c(country = 0.002, state = 0.002, county = 0.002),
+                        combine = "or")
+  )
+  if (is.null(set)) return(names(sets))
+  set <- match.arg(set, names(sets))
+  sets[[set]]
+}
+
+#' Decimal places of submitted coordinates
+#'
+#' Internal.  As the service counts them: trailing zeros do not count, so "12.50" is one
+#' decimal place.  A coordinate passed as a number has already lost its trailing zeros,
+#' which is why \code{\link{GVS_local}} keeps the column as supplied; pass character
+#' coordinates to reproduce the service exactly.
+#' @keywords internal
+#' @noRd
+gvs_decimal_places <- function(x) {
+  if (is.numeric(x)) x <- formatC(x, format = "f", digits = 12, drop0trailing = FALSE)
+  x <- trimws(as.character(x))
+  sci <- grepl("[eE]", x)
+  if (any(sci, na.rm = TRUE)) {
+    x[which(sci)] <- formatC(suppressWarnings(as.numeric(x[which(sci)])), format = "f",
+                             digits = 12, drop0trailing = FALSE)
+  }
+  frac <- sub("^[^.]*[.]?", "", x)
+  frac[!grepl("[.]", x)] <- ""
+  frac <- sub("0+$", "", frac)
+  out <- nchar(frac)
+  out[is.na(x)] <- NA_integer_
+  as.integer(out)
+}
+
+#' Inherent uncertainty (m) of a coordinate, from its decimal places
+#'
+#' Internal.  The service's \code{coordinate_inherent_uncertainty_m}: the diagonal of the
+#' cell the rounding leaves the point in, each axis taking its own decimal places, at
+#' 111,320 m per degree with longitude converging as cos(latitude).  Reproduces the
+#' service on its benchmark to within 1\% for every record but one (a coordinate written
+#' in scientific notation).
+#' @keywords internal
+#' @noRd
+gvs_coordinate_uncertainty_m <- function(lat, dp_lat, dp_lon) {
+  m_per_deg <- 111320
+  sqrt((10^(-dp_lat) * m_per_deg)^2 + (10^(-dp_lon) * m_per_deg * cos(lat * pi / 180))^2)
+}
+
+#' Check user-supplied centroid thresholds
+#' @keywords internal
+#' @noRd
+gvs_check_thresholds <- function(thr_abs, thr_rel) {
+  if (is.null(thr_abs) && is.null(thr_rel)) {
+    stop("Give thr_abs, thr_rel or both.", call. = FALSE)
+  }
+  for (nm in c("thr_abs", "thr_rel")) {
+    thr <- get(nm)
+    if (is.null(thr)) next
+    if (!is.numeric(thr) || is.null(names(thr)) || !all(names(thr) %in% gvs_levels()) ||
+        anyDuplicated(names(thr)) || any(!is.finite(thr)) || any(thr < 0)) {
+      stop(nm, " should be NULL or a named vector of non-negative numbers with names among ",
+           "'country', 'state', 'county'.", call. = FALSE)
+    }
+  }
+  invisible(TRUE)
+}
+
+#' Centroid test at one level: absolute, relative, or both combined
+#'
+#' Internal.  A threshold of NULL means that kind of distance is not tested; with
+#' neither, nothing is flagged.  Missing distances never flag.
+#' @keywords internal
+#' @noRd
+gvs_centroid_flag <- function(dist, rel, a = NULL, r = NULL, combine = c("or", "and")) {
+  combine <- match.arg(combine)
+  le <- function(x, t) { z <- x <= t; z[is.na(z)] <- FALSE; z }
+  fa <- if (is.null(a)) NULL else le(dist, a)
+  fr <- if (is.null(r)) NULL else le(rel, r)
+  if (is.null(fa) && is.null(fr)) return(rep(FALSE, length(dist)))
+  if (is.null(fa)) return(fr)
+  if (is.null(fr)) return(fa)
+  if (combine == "or") fa | fr else fa & fr
+}
+
+#' Parse an optional date column (years, ISO dates, Dates)
+#' @keywords internal
+#' @noRd
+gvs_parse_record_date <- function(x) {
+  if (inherits(x, "Date")) return(x)
+  x <- trimws(as.character(x))
+  out <- rep(as.Date(NA), length(x))
+  year <- !is.na(x) & grepl("^[0-9]{3,4}$", x)
+  if (any(year)) out[year] <- as.Date(sprintf("%04d-07-01", as.integer(x[year])))
+  iso <- !is.na(x) & !year & grepl("^[0-9]{4}-[0-9]{2}", x)
+  if (any(iso)) out[iso] <- suppressWarnings(as.Date(substr(paste0(x[iso], "-01"), 1, 10)))
+  out
+}
+
+#' Distances from points to the six centroids of their division at one level
+#'
+#' Internal.  As the service: the nearest of the six centroid types by absolute
+#' (ellipsoidal) distance is reported, with its relative distance (planar degrees
+#' over that centroid's maximum distance within the division).  Where a base type
+#' and its _main counterpart coincide the _main label is reported.
+#'
+#' \code{m} indexes the row of \code{tab} for each point (default: matched on
+#' \code{gids}).  Distances within 0.1\% of an absolute threshold in \code{thr}
+#' are recomputed with \code{geosphere::distGeo}, so comparisons with the
+#' thresholds are those of the service; elsewhere the Andoyer-Lambert
+#' approximation is used (see \code{gvs_distance_km}).
+#' @keywords internal
+#' @noRd
+gvs_level_distances <- function(lon, lat, gids, tab, m = match(gids, tab$gid), thr = NULL) {
+  n <- length(lon)
+  eval_order <- c(4, 5, 6, 1, 2, 3)
+  eval_label <- c("std_main", "pos_main", "bb_main", "std", "pos", "bb")
+  D <- R <- matrix(NA_real_, n, 6)
+  has <- which(!is.na(m))
+  mh <- m[has]
+  for (col in seq_len(6)) {
+    i <- eval_order[col]
+    clon <- tab[[paste0("c", i, "_lon")]][mh]
+    clat <- tab[[paste0("c", i, "_lat")]][mh]
+    dm <- tab[[paste0("c", i, "_dmax")]][mh]
+    d <- gvs_distance_km(lon[has], lat[has], clon, clat)
+    thr <- thr[is.finite(thr)]
+    if (length(thr) && requireNamespace("geosphere", quietly = TRUE)) {
+      near <- which(Reduce(`|`, lapply(thr, function(t) abs(d - t) <= 1e-5 * t)))
+      if (length(near)) {
+        d[near] <- geosphere::distGeo(cbind(lon[has][near], lat[has][near]), cbind(clon[near], clat[near])) / 1000
+      }
+    }
+    D[has, col] <- d
+    R[has, col] <- sqrt((lon[has] - clon)^2 + (lat[has] - clat)^2) / dm
+  }
+  k <- max.col(-replace(D, is.na(D), Inf), ties.method = "first")
+  idx <- cbind(seq_len(n), k)
+  ok <- !is.na(m)
+  data.frame(dist = ifelse(ok, D[idx], NA), rel = ifelse(ok, R[idx], NA),
+             type = ifelse(ok, eval_label[k], NA_character_), stringsAsFactors = FALSE)
+}
+
+#' Ellipsoidal (WGS84) distance in km
+#'
+#' Internal.  Andoyer-Lambert: the central angle on the auxiliary sphere of
+#' reduced latitudes, corrected for flattening to first order.  Vectorised, and
+#' within about 1e-5 of the geodesic (Karney) distance at the distances that
+#' matter for centroids; exact comparisons with thresholds are made in
+#' \code{gvs_level_distances}.
+#' @keywords internal
+#' @noRd
+gvs_distance_km <- function(lon1, lat1, lon2, lat2) {
+  if (!length(lon1)) return(numeric(0))
+  a <- 6378.137
+  f <- 1 / 298.257223563
+  r <- pi / 180
+  b1 <- atan((1 - f) * tan(lat1 * r))
+  b2 <- atan((1 - f) * tan(lat2 * r))
+  h <- sin((b2 - b1) / 2)^2 + cos(b1) * cos(b2) * sin((lon2 - lon1) * r / 2)^2
+  sig <- 2 * asin(pmin(1, sqrt(h)))
+  P <- (b1 + b2) / 2
+  Q <- (b2 - b1) / 2
+  X <- (sig - sin(sig)) * sin(P)^2 * cos(Q)^2 / cos(sig / 2)^2
+  Y <- (sig + sin(sig)) * cos(P)^2 * sin(Q)^2 / sin(sig / 2)^2
+  d <- a * (sig - f / 2 * (X + Y))
+  d[sig == 0] <- 0
+  d
+}
+
+#' Pairs of points and centroids within a latitude/longitude window
+#'
+#' Internal.  First a coarse mask (0.1 degree cells) of the union of all windows
+#' discards the points that are near no centroid, in one vectorised lookup; the
+#' rest are sorted by latitude and each centroid takes the slice within
+#' \code{rlat} of it, keeping those within \code{rlon} in longitude.
+#' @return data.frame(p = point index, j = centroid index)
+#' @keywords internal
+#' @noRd
+gvs_pairs_within <- function(lon, lat, clon, clat, rlat, rlon, cell = 0.1) {
+  none <- data.frame(p = integer(0), j = integer(0))
+  good <- is.finite(clon) & is.finite(clat)
+  ok <- which(is.finite(lon) & is.finite(lat) & abs(lon) <= 180 & abs(lat) <= 90)
+  if (!length(ok) || !any(good)) return(none)
+
+  nx <- round(360 / cell)
+  ny <- round(180 / cell)
+  mask <- matrix(FALSE, nrow = nx, ncol = ny)
+  for (j in which(good)) {
+    x0 <- max(1L, floor((clon[j] - rlon[j] + 180) / cell))
+    x1 <- min(nx, floor((clon[j] + rlon[j] + 180) / cell) + 1L)
+    y0 <- max(1L, floor((clat[j] - rlat[j] + 90) / cell))
+    y1 <- min(ny, floor((clat[j] + rlat[j] + 90) / cell) + 1L)
+    if (x1 >= x0 && y1 >= y0) mask[x0:x1, y0:y1] <- TRUE
+  }
+  cx <- pmin(nx, floor((lon[ok] + 180) / cell) + 1L)
+  cy <- pmin(ny, floor((lat[ok] + 90) / cell) + 1L)
+  ok <- ok[mask[cbind(cx, cy)]]
+  if (!length(ok)) return(none)
+
+  o <- ok[order(lat[ok])]
+  sl <- lat[o]
+  ps <- js <- vector("list", length(clon))
+  for (j in which(good)) {
+    lo <- findInterval(clat[j] - rlat[j], sl, left.open = TRUE) + 1L
+    hi <- findInterval(clat[j] + rlat[j], sl)
+    if (hi < lo) next
+    idx <- o[lo:hi]
+    idx <- idx[abs(lon[idx] - clon[j]) <= rlon[j]]
+    if (length(idx)) {
+      ps[[j]] <- idx
+      js[[j]] <- rep.int(j, length(idx))
+    }
+  }
+  data.frame(p = as.integer(unlist(ps)), j = as.integer(unlist(js)))
+}
+
+gvs_session <- new.env(parent = emptyenv())
+
+#' Read the CShapes versions and geometry from the cache, or NULL
+#'
+#' Internal.  Cached for the session, keyed on the files' paths and times.
+#' @keywords internal
+#' @noRd
+gvs_read_cshapes <- function(dir, geometry = TRUE) {
+  vf <- file.path(dir, "cshapes-versions.gz.parquet")
+  gf <- file.path(dir, "cshapes-geometry.gpkg")
+  if (!file.exists(vf) || (geometry && !file.exists(gf))) return(NULL)
+  stamp <- paste(normalizePath(vf), file.mtime(vf), if (geometry) file.mtime(gf))
+  hit <- gvs_session$cshapes
+  if (!is.null(hit) && identical(hit$stamp, stamp)) return(hit$value)
+  v <- as.data.frame(nanoparquet::read_parquet(vf))
+  v$valid_from <- as.Date(v$valid_from)
+  v$valid_to <- as.Date(v$valid_to)
+  if (geometry) {
+    g <- sf::st_read(gf, quiet = TRUE)
+    g <- g[match(v$version_id, g$version_id), ]
+    attr(v, "geometry") <- sf::st_geometry(g)
+  }
+  gvs_session$cshapes <- list(stamp = stamp, value = v)
+  v
+}
+
+#' Which (point, version) pairs have the point inside the version's polygon
+#'
+#' Internal.  Planar (GEOS) point-in-polygon, as for GADM.
+#' @return logical, one per pair
+#' @keywords internal
+#' @noRd
+gvs_pairs_inside <- function(lon, lat, p, j, geom) {
+  if (!length(p)) return(logical(0))
+  old_s2 <- sf::sf_use_s2()
+  suppressMessages(sf::sf_use_s2(FALSE))
+  on.exit(suppressMessages(sf::sf_use_s2(old_s2)), add = TRUE)
+  up <- unique(p)
+  uj <- unique(j)
+  pts <- sf::st_as_sf(data.frame(lon = lon[up], lat = lat[up]), coords = c("lon", "lat"), crs = sf::st_crs(geom))
+  hits <- suppressMessages(sf::st_intersects(pts, geom[uj]))
+  hp <- rep(up, lengths(hits))
+  hj <- uj[unlist(hits, use.names = FALSE)]
+  # pair keys as doubles: p < 2^31 and j < 2^20
+  (as.numeric(p) * 1048576 + j) %in% (as.numeric(hp) * 1048576 + hj)
+}
+
+#' Centroid test against historical countries containing each point
+#'
+#' Internal.  A point is tested against a CShapes version only if it is within
+#' reach of one of that version's six centroids under the country thresholds
+#' (reach in latitude = the larger of the absolute threshold at 110.5 km per
+#' degree and the relative threshold times that centroid's maximum distance; in
+#' longitude the absolute part is widened by 1/cos(latitude)).  Qualifying pairs
+#' are then required to have the point inside the version's polygon and, when the
+#' record has a date, the version to be valid at that date (± tolerance).  Of
+#' several qualifying versions the earliest is reported.  As for current GADM
+#' divisions, windows do not wrap the antimeridian.
+#' @keywords internal
+#' @noRd
+gvs_historical_centroids <- function(lon, lat, dates, dir, thr_abs, thr_rel, tolerance_years,
+                                     combine = "or") {
+  n <- length(lon)
+  out <- data.frame(version_id = rep(NA_character_, n), country_name = NA_character_,
+                    valid_from = as.Date(NA), valid_to = as.Date(NA),
+                    dist = NA_real_, rel = NA_real_, type = NA_character_, stringsAsFactors = FALSE)
+  if (is.null(thr_abs) && is.null(thr_rel)) return(out)
+  v <- gvs_read_cshapes(dir)
+  if (is.null(v)) return(out)
+
+  # unique coordinates (and dates) only
+  dd0 <- if (is.null(dates)) gvs_dedupe(lon, lat) else gvs_dedupe(lon, lat, as.numeric(dates))
+  u <- dd0$first
+  ulon <- lon[u]
+  ulat <- lat[u]
+  udate <- if (is.null(dates)) rep(as.Date(NA), length(u)) else dates[u]
+
+  pairs <- list()
+  for (i in 1:6) {
+    clon <- v[[paste0("c", i, "_lon")]]
+    clat <- v[[paste0("c", i, "_lat")]]
+    rel_deg <- if (is.null(thr_rel)) 0 else thr_rel * v[[paste0("c", i, "_dmax")]]
+    abs_deg <- if (is.null(thr_abs)) 0 else thr_abs / 110.5
+    rlat <- pmax(abs_deg, rel_deg)
+    rlon <- pmax(abs_deg / pmax(cos(pmin(abs(clat) + rlat, 89.9) * pi / 180), 0.01), rel_deg)
+    pairs[[i]] <- gvs_pairs_within(ulon, ulat, clon, clat, rlat, rlon)
+  }
+  pr <- do.call(rbind, pairs)
+  pr <- pr[!duplicated(as.numeric(pr$p) * 1048576 + pr$j), , drop = FALSE]
+  if (!nrow(pr)) return(out)
+
+  # dates
+  tol <- 365.25 * tolerance_years
+  d <- udate[pr$p]
+  pr <- pr[is.na(d) | (v$valid_from[pr$j] - tol <= d & v$valid_to[pr$j] + tol >= d), , drop = FALSE]
+  if (!nrow(pr)) return(out)
+
+  # distances to that version's centroids, and the thresholds
+  dd <- gvs_level_distances(ulon[pr$p], ulat[pr$p], NULL, v, m = pr$j, thr = thr_abs)
+  q <- gvs_centroid_flag(dd$dist, dd$rel, thr_abs, thr_rel, combine)
+  pr <- cbind(pr, dd)[q, , drop = FALSE]
+  if (!nrow(pr)) return(out)
+
+  pr <- pr[gvs_pairs_inside(ulon, ulat, pr$p, pr$j, attr(v, "geometry")), , drop = FALSE]
+  if (!nrow(pr)) return(out)
+  pr <- pr[order(pr$p, v$valid_from[pr$j]), , drop = FALSE]
+  pr <- pr[!duplicated(pr$p), , drop = FALSE]
+
+  res <- out[seq_along(u), , drop = FALSE]
+  res$version_id[pr$p] <- v$version_id[pr$j]
+  res$country_name[pr$p] <- v$country_name[pr$j]
+  res$valid_from[pr$p] <- v$valid_from[pr$j]
+  res$valid_to[pr$p] <- v$valid_to[pr$j]
+  res$dist[pr$p] <- pr$dist
+  res$rel[pr$p] <- pr$rel
+  res$type[pr$p] <- pr$type
+  res <- res[dd0$group, , drop = FALSE]
+  rownames(res) <- NULL
+  res
+}
+
+#' Is each point inside the political division its record names, now or ever?
+#'
+#' Internal.  Uses the GNRS resolution of the same rows.
+#' \describe{
+#'   \item{geovalid_country_current}{the point's current GADM country is the
+#'     resolved country (GNRS gid_0).  NA when GNRS resolved a former country or
+#'     no country.}
+#'   \item{geovalid_country_any}{the point is in its current country as above, or
+#'     in a CShapes version of the named entity (only versions valid at the
+#'     record's date when it has one), or, for a former country without CShapes
+#'     geometry (e.g. the Netherlands Antilles), in one of its current successors.}
+#'   \item{geovalid_basis}{"current GADM", the CShapes version id, "successor:
+#'     <keys>", or NA.}
+#'   \item{geovalid_state_current, geovalid_county_current}{the point's current
+#'     GADM state / county are the resolved ones.}
+#'   \item{subnational_status}{"valid", "invalid", or "unverifiable": GNRS could
+#'     not place the named state/county in a successor, or they disagree with
+#'     today's divisions while the country is valid only historically.}
+#'   \item{geovalid}{the default for BIEN uses: the country was ever accurate and
+#'     the sub-national divisions are not contradicted (valid or unverifiable).}
+#' }
+#' @keywords internal
+#' @noRd
+gvs_geovalidity <- function(lon, lat, dates, loc, gnrs, dir, tolerance_years) {
+  n <- length(lon)
+  s <- function(x) { x <- as.character(x); x[is.na(x)] <- ""; x }
+  g0 <- s(gnrs$gid_0); g1 <- s(gnrs$gid_1); g2 <- s(gnrs$gid_2)
+  ek <- if ("entity_key" %in% names(gnrs)) s(gnrs$entity_key) else rep("", n)
+  hist <- if ("is_historical" %in% names(gnrs)) gnrs$is_historical %in% TRUE else rep(FALSE, n)
+  succ <- if ("successors" %in% names(gnrs)) s(gnrs$successors) else rep("", n)
+  gsub_status <- if ("subnational_status" %in% names(gnrs)) s(gnrs$subnational_status) else rep("", n)
+  named <- nzchar(ek) | nzchar(g0)
+
+  cur_country <- ifelse(hist | !nzchar(g0), NA, s(loc$gid_0) == g0)
+  any_country <- cur_country %in% TRUE
+  basis <- ifelse(any_country, "current GADM", NA_character_)
+
+  # CShapes versions of the named entity containing the point
+  need <- which(!any_country & nzchar(ek) & is.finite(lon) & is.finite(lat))
+  pf <- file.path(dir, "history-periods.gz.parquet")
+  v <- if (length(need) && file.exists(pf)) gvs_read_cshapes(dir) else NULL
+  if (!is.null(v)) {
+    per <- as.data.frame(nanoparquet::read_parquet(pf))
+    per$from <- as.Date(per$from)
+    per$to <- as.Date(per$to)
+    dd0 <- if (is.null(dates)) gvs_dedupe(lon[need], lat[need], ek[need]) else
+      gvs_dedupe(lon[need], lat[need], ek[need], as.numeric(dates[need]))
+    u <- need[dd0$first]
+    # candidate versions: those whose gwcode ever carried the entity
+    cand <- merge(data.frame(p = u, entity_key = ek[u]), per[, c("entity_key", "gwcode", "from", "to")],
+                  by = "entity_key")
+    if (nrow(cand)) {
+      cand <- merge(cand, data.frame(j = seq_len(nrow(v)), gwcode = v$gwcode), by = "gwcode")
+      # the version overlaps a period in which the gwcode was that entity
+      cand <- cand[v$valid_from[cand$j] <= cand$to & v$valid_to[cand$j] >= cand$from, , drop = FALSE]
+      if (!is.null(dates) && nrow(cand)) {
+        tol <- 365.25 * tolerance_years
+        d <- dates[cand$p]
+        cand <- cand[is.na(d) | (v$valid_from[cand$j] - tol <= d & v$valid_to[cand$j] + tol >= d), , drop = FALSE]
+      }
+      cand <- cand[!duplicated(as.numeric(cand$p) * 1048576 + cand$j), c("p", "j"), drop = FALSE]
+      if (nrow(cand)) {
+        cand <- cand[gvs_pairs_inside(lon, lat, cand$p, cand$j, attr(v, "geometry")), , drop = FALSE]
+        cand <- cand[order(cand$p, v$valid_from[cand$j]), , drop = FALSE]
+        cand <- cand[!duplicated(cand$p), , drop = FALSE]
+        hit <- match(u[dd0$group], cand$p)
+        ok <- !is.na(hit)
+        any_country[need[ok]] <- TRUE
+        basis[need[ok]] <- v$version_id[cand$j[hit[ok]]]
+      }
+    }
+  }
+
+  # a former country without CShapes geometry: its current successors
+  with_geometry <- if (file.exists(pf)) unique(nanoparquet::read_parquet(pf)$entity_key) else character(0)
+  via <- which(!any_country & hist & nzchar(succ) & !(ek %in% with_geometry) & !is.na(loc$gid_0))
+  if (length(via) && requireNamespace("countrycode", quietly = TRUE)) {
+    in_succ <- vapply(via, function(i) {
+      keys <- strsplit(succ[i], ";", fixed = TRUE)[[1]]
+      iso3 <- suppressWarnings(countrycode::countrycode(keys, "iso2c", "iso3c", warn = FALSE))
+      iso3[keys == "XK"] <- "XKO"
+      loc$gid_0[i] %in% iso3
+    }, logical(1))
+    any_country[via[in_succ]] <- TRUE
+    basis[via[in_succ]] <- paste("successor:", succ[via[in_succ]])
+  }
+
+  st_cur <- ifelse(!nzchar(g1), NA, s(loc$gid_1) == g1)
+  co_cur <- ifelse(!nzchar(g2), NA, s(loc$gid_2) == g2)
+  sub_given <- nzchar(g1) | nzchar(g2) | nzchar(gsub_status)
+  contradicted <- (st_cur %in% FALSE) | (co_cur %in% FALSE)
+  sub_status <- ifelse(!sub_given, NA_character_,
+                ifelse(gsub_status == "unverifiable", "unverifiable",
+                ifelse(!contradicted, "valid",
+                ifelse(any_country & !(cur_country %in% TRUE), "unverifiable", "invalid"))))
+  data.frame(
+    geovalid_country_current = cur_country,
+    geovalid_country_any = ifelse(named, any_country, NA),
+    geovalid_basis = basis,
+    geovalid_state_current = st_cur,
+    geovalid_county_current = co_cur,
+    subnational_status = sub_status,
+    geovalid = ifelse(named, any_country & !(sub_status %in% "invalid"), NA),
+    stringsAsFactors = FALSE
+  )
+}
diff --git a/R/local_build.R b/R/local_build.R
new file mode 100644
index 0000000..0471cc4
--- /dev/null
+++ b/R/local_build.R
@@ -0,0 +1,349 @@
+#' Install the reference data GVS_local() resolves against
+#'
+#' Builds the local copy of the reference data, so that coordinates can be
+#' validated with no calls to the GVS API.  Two components:
+#'
+#' \code{"gadm"} is the GADM release the divisions come from: the world
+#' GeoPackage (a large download), from which the attributes of levels 0 to 2 are
+#' kept.  This component is shared with the GNRS package and is skipped if that
+#' package has already built it.  A GeoPackage already on disk can be used
+#' instead of downloading one, with \code{gpkg}.
+#'
+#' \code{"centroids"} is what the resolver reads: for every country,
+#' state/province and county/parish, the six centroid types the service uses
+#' (centre of mass, point on surface and bounding-box centre, each for the whole
+#' division and for its largest part) and, for each, the greatest distance from
+#' that centroid to the division's boundary.  Relative distances are measured
+#' against those maxima, which is how a point 2 km from the centroid of a small
+#' county and a point 2 km from the centroid of a large country are told apart.
+#' Deriving them needs the GADM geometry, so the GeoPackage must be present (or
+#' downloadable) even though only these tables are kept.
+#'
+#' @param sources Components to build.  Defaults to both.
+#' @param dir Cache directory.  Defaults to the standard user cache location,
+#'   shared with GNRS.
+#' @param gpkg Optional path to a GADM world GeoPackage already on disk.  Used
+#'   instead of downloading, and not deleted afterwards.
+#' @param gadm_layer Layer within the GeoPackage.  GADM 4.1 ships one layer of
+#'   all levels.
+#' @param overwrite Rebuild components that are already built?
+#' @param keep_archive Keep the downloaded GADM archive after building?
+#' @param quiet Suppress progress messages?
+#' @return The status table, invisibly.
+#' @seealso \code{\link{GVS_local}}, \code{\link{GVS_local_status}}
+#' @export
+#' @examples \dontrun{
+#' # Using a GeoPackage already on disk
+#' GVS_local_build(gpkg = "gadm_410.gpkg")
+#'
+#' # Downloading GADM (about 1.4 GB)
+#' GVS_local_build()
+#' }
+GVS_local_build <- function(sources = c("gadm", "centroids"),
+                            dir = gvs_cache_dir(create = TRUE),
+                            gpkg = NULL,
+                            gadm_layer = "gadm_410",
+                            overwrite = FALSE,
+                            keep_archive = FALSE,
+                            quiet = FALSE) {
+  registry <- gvs_builtin_registry()
+  unknown <- setdiff(sources, names(registry))
+  if (length(unknown) > 0) {
+    stop("Unknown source(s): ", paste(unknown, collapse = ", "), call. = FALSE)
+  }
+  for (pkg in c("sf", "terra", "nanoparquet")) {
+    if (!requireNamespace(pkg, quietly = TRUE)) {
+      stop(
+        "Building the local reference data needs the '", pkg, "' package. ",
+        "Install it with install.packages(\"", pkg, "\") and build again.",
+        call. = FALSE
+      )
+    }
+  }
+  if (!dir.exists(dir)) dir.create(dir, recursive = TRUE, showWarnings = FALSE)
+
+  # The GeoPackage is needed by both components; it is downloaded once, used by
+  # whichever of them runs, and removed afterwards unless the caller supplied it
+  supplied <- !is.null(gpkg)
+  if (supplied && !file.exists(gpkg)) {
+    stop("No GeoPackage at: ", gpkg, call. = FALSE)
+  }
+  needs_geometry <- "centroids" %in% sources &&
+    (overwrite || !gvs_is_built("centroids", dir))
+  needs_attributes <- "gadm" %in% sources &&
+    (overwrite || !gvs_is_built("gadm", dir))
+
+  if (!supplied && (needs_geometry || needs_attributes)) {
+    gpkg <- gvs_fetch_gadm(dir = dir, quiet = quiet)
+    on.exit({
+      unlink(gpkg)
+      if (!keep_archive) unlink(gvs_gadm_archive_path(dir))
+    }, add = TRUE)
+  } else if (!needs_geometry && !needs_attributes && !quiet) {
+    message("Already built; nothing to do (use overwrite = TRUE to rebuild).")
+  }
+
+  if (needs_attributes) {
+    gvs_build_gadm_attributes(gpkg, gadm_layer, dir = dir, quiet = quiet)
+  }
+  if (needs_geometry) {
+    gvs_build_centroids(gpkg, gadm_layer, dir = dir, quiet = quiet)
+  }
+
+  invisible(GVS_local_status(dir))
+}
+
+#' Download and extract the GADM world GeoPackage
+#'
+#' Internal.  Returns the path of the extracted GeoPackage.
+#' @keywords internal
+#' @noRd
+gvs_fetch_gadm <- function(dir = gvs_cache_dir(), quiet = FALSE) {
+  spec <- gvs_builtin_registry()$gadm
+  archive <- gvs_gadm_archive_path(dir)
+
+  if (file.exists(archive)) {
+    if (!quiet) message("Using cached download of GADM ", spec$version)
+  } else {
+    if (!quiet) {
+      message(
+        "Downloading ", spec$full_name, " ", spec$version,
+        " (about ", spec$download_mb, " MB) ..."
+      )
+    }
+    partial <- paste0(archive, ".part")
+    old <- options(timeout = max(7200, getOption("timeout")))
+    on.exit(options(old), add = TRUE)
+    status <- utils::download.file(spec$url, partial, mode = "wb", quiet = quiet, cacheOK = FALSE)
+    if (status != 0 || !file.exists(partial)) {
+      unlink(partial)
+      stop("Download failed for GADM.", call. = FALSE)
+    }
+    file.rename(partial, archive)
+  }
+
+  if (!quiet) message("Extracting the GeoPackage ...")
+  members <- utils::unzip(archive, list = TRUE)$Name
+  member <- members[grepl("\\.gpkg$", members)][1]
+  if (is.na(member)) {
+    stop("No GeoPackage found in the GADM archive.", call. = FALSE)
+  }
+  utils::unzip(archive, files = member, exdir = dir, overwrite = TRUE)
+  file.path(dir, member)
+}
+
+#' Keep the attributes of GADM levels 0 to 2
+#'
+#' Internal.  The same table GNRS builds, written to the same path, so that
+#' either package can install it.  Read without touching the geometry.
+#' @keywords internal
+#' @noRd
+gvs_build_gadm_attributes <- function(gpkg, gadm_layer, dir = gvs_cache_dir(), quiet = FALSE) {
+  if (!quiet) message("Reading the GADM attributes ...")
+  fields <- c("GID_0", "COUNTRY", "GID_1", "NAME_1", "GID_2", "NAME_2")
+  d <- sf::st_read(
+    gpkg,
+    query = sprintf('SELECT %s FROM "%s"', paste(fields, collapse = ", "), gadm_layer),
+    quiet = TRUE
+  )
+  d <- sf::st_drop_geometry(d)
+  blank <- function(x) is.na(x) | !nzchar(x) | x == "NA"
+  divisions <- rbind(
+    unique(data.frame(
+      level = 0L, gid = d$GID_0, name = d$COUNTRY,
+      gid_0 = d$GID_0, gid_1 = NA_character_, gid_2 = NA_character_
+    )),
+    unique(data.frame(
+      level = 1L, gid = d$GID_1, name = d$NAME_1,
+      gid_0 = d$GID_0, gid_1 = d$GID_1, gid_2 = NA_character_
+    )[!blank(d$GID_1), ]),
+    unique(data.frame(
+      level = 2L, gid = d$GID_2, name = d$NAME_2,
+      gid_0 = d$GID_0, gid_1 = d$GID_1, gid_2 = d$GID_2
+    )[!blank(d$GID_2), ])
+  )
+  divisions <- divisions[!duplicated(divisions$gid), ]
+  nanoparquet::write_parquet(divisions, gvs_gadm_path(dir), compression = "gzip")
+
+  spec <- gvs_builtin_registry()$gadm
+  provenance <- list(
+    source = "gadm", full_name = spec$full_name, version = spec$version,
+    url = spec$url, license = spec$license, publisher = spec$publisher,
+    downloaded = as.character(Sys.Date()),
+    n_admin0 = sum(divisions$level == 0L),
+    n_admin1 = sum(divisions$level == 1L),
+    n_admin2 = sum(divisions$level == 2L)
+  )
+  saveRDS(provenance, gvs_provenance_path("gadm", dir))
+  if (!quiet) {
+    message(
+      "  ", provenance$n_admin0, " countries, ", provenance$n_admin1,
+      " level-1 and ", provenance$n_admin2, " level-2 divisions"
+    )
+  }
+  invisible(provenance)
+}
+
+#' The six centroids of one set of dissolved divisions
+#'
+#' Internal.  \code{geom} is one dissolved division per element.  Column order
+#' follows the service's own: centre of mass, point on surface, bounding-box
+#' centre, then the same three computed on the division's largest part.  For a
+#' single-part division the two halves coincide.
+#'
+#' Each centroid carries the greatest distance from it to the division's
+#' boundary, which relative distances are measured against.  Distances are
+#' Euclidean in degrees, as the service's relative distances are, and the
+#' farthest boundary point of a polygon is always a vertex.  They are measured
+#' to the WHOLE division even for the largest-part centroids, as the service
+#' measures them (checked against its own tables: for the United States 279.2
+#' degrees to the farthest vertex of the country, not 32.6 within the largest
+#' part).  Distances and centroids are deliberately planar in WGS84 longitude and
+#' latitude, including for divisions that span the antimeridian (Fiji's maximum is
+#' 357.97 degrees).  The service is built to detect centroids that users assigned
+#' naively, and a naive planar calculation for a dateline-straddling division puts
+#' its centroid far from the division (near longitude 0); correcting for the
+#' antimeridian here would stop those centroids from being recognised.  The
+#' point-on-surface and largest-part centroid types cover users who avoided that
+#' mistake.
+#' @keywords internal
+#' @noRd
+gvs_six_centroids <- function(geom) {
+  n <- length(geom)
+  out <- matrix(NA_real_, n, 18)
+
+  main <- gvs_largest_part(geom)
+
+  for (half in 1:2) {
+    g <- if (half == 1) geom else main
+    cent <- suppressWarnings(sf::st_coordinates(sf::st_centroid(g)))
+    pos <- suppressWarnings(sf::st_coordinates(sf::st_point_on_surface(g)))
+    bb <- vapply(g, function(x) {
+      b <- sf::st_bbox(x)
+      c((b[["xmin"]] + b[["xmax"]]) / 2, (b[["ymin"]] + b[["ymax"]]) / 2)
+    }, numeric(2))
+    types <- list(cent[, 1:2, drop = FALSE], pos[, 1:2, drop = FALSE], t(bb))
+    for (k in seq_along(types)) {
+      col <- (half - 1) * 9 + (k - 1) * 3 + 1
+      out[, col] <- types[[k]][, 1]
+      out[, col + 1] <- types[[k]][, 2]
+      out[, col + 2] <- gvs_max_vertex_distance(geom, types[[k]])
+    }
+  }
+
+  # c1..c6 as the service orders them: 1 std, 2 pos, 3 bb, 4 std_main, 5 pos_main, 6 bb_main
+  colnames(out) <- paste0("c", rep(1:6, each = 3), c("_lon", "_lat", "_dmax"))
+  as.data.frame(out)
+}
+
+#' The largest part of each (multi)polygon, by planar area
+#' @keywords internal
+#' @noRd
+gvs_largest_part <- function(geom) {
+  parts <- lapply(geom, function(g) {
+    pieces <- suppressWarnings(sf::st_cast(sf::st_sfc(g), "POLYGON"))
+    if (length(pieces) <= 1) {
+      return(g)
+    }
+    areas <- suppressWarnings(as.numeric(sf::st_area(sf::st_set_crs(pieces, NA))))
+    pieces[[which.max(areas)]]
+  })
+  sf::st_sfc(parts, crs = sf::st_crs(geom))
+}
+
+#' Greatest Euclidean-degree distance from each centroid to its division's vertices
+#' @keywords internal
+#' @noRd
+gvs_max_vertex_distance <- function(geom, centres) {
+  vapply(seq_along(geom), function(i) {
+    xy <- sf::st_coordinates(sf::st_sfc(geom[[i]]))
+    if (!nrow(xy) || any(is.na(centres[i, ]))) {
+      return(NA_real_)
+    }
+    max(sqrt((xy[, 1] - centres[i, 1])^2 + (xy[, 2] - centres[i, 2])^2))
+  }, numeric(1))
+}
+
+#' Derive the centroid tables from the GADM geometry
+#'
+#' Internal.  Countries are read one at a time, so peak memory is one country's
+#' geometry rather than the world's.  GADM ships one row per finest division, so
+#' the divisions of each level are dissolved from those rows: a level-2 division
+#' is the union of its level-3 and finer children, and a country the union of
+#' everything in it.  Dissolving matters for the largest-part centroids, where
+#' the largest part of a country is an island rather than whichever administrative
+#' piece of it happens to be biggest.
+#' @keywords internal
+#' @noRd
+gvs_build_centroids <- function(gpkg, gadm_layer, dir = gvs_cache_dir(), quiet = FALSE) {
+  old_s2 <- sf::sf_use_s2()
+  suppressMessages(sf::sf_use_s2(FALSE))
+  on.exit(suppressMessages(sf::sf_use_s2(old_s2)), add = TRUE)
+
+  countries <- sf::st_read(
+    gpkg,
+    query = sprintf('SELECT DISTINCT GID_0 FROM "%s" ORDER BY GID_0', gadm_layer),
+    quiet = TRUE
+  )
+  countries <- sf::st_drop_geometry(countries)$GID_0
+  if (!quiet) message("Deriving centroids for ", length(countries), " countries ...")
+
+  levels <- gvs_levels()
+  field <- c(country = "GID_0", state = "GID_1", county = "GID_2")
+  acc <- stats::setNames(vector("list", length(levels)), levels)
+  t0 <- Sys.time()
+
+  for (i in seq_along(countries)) {
+    g0 <- countries[i]
+    v <- terra::vect(
+      gpkg,
+      query = sprintf(
+        'SELECT GID_0, GID_1, GID_2, geom FROM "%s" WHERE GID_0 = \'%s\'',
+        gadm_layer, g0
+      )
+    )
+    if (!nrow(v)) next
+
+    for (lev in levels) {
+      key <- field[[lev]]
+      ids <- as.character(terra::values(v)[[key]])
+      keep <- !is.na(ids) & nzchar(ids) & ids != "NA"
+      if (!any(keep)) next
+      vk <- v[keep, ]
+      # dissolve to one geometry per division of this level
+      agg <- terra::aggregate(vk, by = key)
+      s <- sf::st_as_sf(agg)
+      gid <- as.character(s[[key]])
+      geom <- sf::st_geometry(s)
+      cents <- gvs_six_centroids(geom)
+      acc[[lev]][[length(acc[[lev]]) + 1L]] <- cbind(data.frame(gid = gid), cents)
+    }
+
+    if (!quiet && (i %% 25 == 0 || i == length(countries))) {
+      el <- as.numeric(difftime(Sys.time(), t0, units = "mins"))
+      message(sprintf(
+        "  %d/%d countries | %.1f min | %s counties so far",
+        i, length(countries), el,
+        format(sum(vapply(acc$county, nrow, integer(1))), big.mark = ",")
+      ))
+    }
+  }
+
+  provenance <- list(
+    source = "centroids",
+    full_name = gvs_builtin_registry()$centroids$full_name,
+    version = gvs_builtin_registry()$gadm$version,
+    gadm_version = gvs_builtin_registry()$gadm$version,
+    downloaded = as.character(Sys.Date()),
+    built_minutes = round(as.numeric(difftime(Sys.time(), t0, units = "mins")), 1)
+  )
+  for (lev in levels) {
+    tab <- do.call(rbind, acc[[lev]])
+    tab <- tab[!duplicated(tab$gid), ]
+    nanoparquet::write_parquet(tab, gvs_centroid_path(lev, dir), compression = "gzip")
+    provenance[[paste0("n_", lev)]] <- nrow(tab)
+    if (!quiet) message("  ", lev, ": ", format(nrow(tab), big.mark = ","), " divisions")
+  }
+  saveRDS(provenance, gvs_provenance_path("centroids", dir))
+  invisible(provenance)
+}
diff --git a/R/local_cache.R b/R/local_cache.R
new file mode 100644
index 0000000..bf3586c
--- /dev/null
+++ b/R/local_cache.R
@@ -0,0 +1,287 @@
+#' Directory holding the local GVS reference data
+#'
+#' The reference data is downloaded and derived on demand rather than shipped
+#' with the package, and is kept in the standard user cache directory so that it
+#' survives package updates and can be removed with \code{GVS_local_remove()}.
+#'
+#' GVS and GNRS derive their reference data from the same GADM release, so they
+#' share one cache directory: whichever package builds GADM first, the other
+#' finds it.  The location follows GNRS's, and either package can install it.
+#' Override with \code{options(GVS.cache_dir = )}, or for both packages at once
+#' with \code{options(GNRS.cache_dir = )}.
+#'
+#' @param create Should the directory be created if it does not exist?
+#' @return Path to the cache directory.
+#' @keywords internal
+#' @noRd
+gvs_cache_dir <- function(create = FALSE) {
+  dir <- getOption(
+    "GVS.cache_dir",
+    getOption("GNRS.cache_dir", tools::R_user_dir("GNRS", which = "cache"))
+  )
+  if (create && !dir.exists(dir)) {
+    dir.create(dir, recursive = TRUE, showWarnings = FALSE)
+  }
+  dir
+}
+
+#' Data sources a local build can fetch or derive
+#'
+#' Internal.  \code{"gadm"} is the release the divisions and centroids come
+#' from, shared with GNRS: the world GeoPackage, from which GNRS keeps the
+#' attributes of levels 0 to 2 and GVS derives its centroid tables.  It is a
+#' large download, so a GeoPackage already on disk can be passed to
+#' \code{GVS_local_build()} instead.
+#'
+#' \code{"centroids"} is what GVS itself resolves against: for every country,
+#' state/province and county/parish, the six centroid types the service uses and
+#' the maximum distance within the division that each is measured against.  It is
+#' derived from the GADM geometry and is a few megabytes.
+#'
+#' @return A named list of source definitions.
+#' @keywords internal
+#' @noRd
+gvs_builtin_registry <- function() {
+  list(
+    gadm = list(
+      source = "gadm",
+      full_name = "GADM administrative areas",
+      publisher = "GADM",
+      version = "4.1",
+      url = "https://geodata.ucdavis.edu/gadm/gadm4.1/gadm_410-gpkg.zip",
+      license = "Free for academic and other non-commercial use; redistribution and commercial use need permission (gadm.org/license.html)",
+      citation = "GADM (2022). Database of Global Administrative Areas, version 4.1. https://gadm.org/",
+      # The world GeoPackage. Only the attributes of levels 0 to 2 are kept, a
+      # few megabytes; extracting it needs about twice its size free.
+      download_mb = 1400,
+      disk_mb = 5
+    ),
+    centroids = list(
+      source = "centroids",
+      full_name = "GVS political division centroids",
+      publisher = "Botanical Information and Ecology Network",
+      url = "https://gvsapi.xyz/gvs_api.php",
+      license = "GPL-3 (software); GADM terms apply to the data",
+      citation = paste(
+        "Maitner B. S., Boyle B. L. & Enquist B. J. Geocoordinate Validation",
+        "Service (GVS). https://github.com/ojalaquellueva/gvs"
+      ),
+      # Derived from the GADM geometry, not downloaded
+      download_mb = 0,
+      disk_mb = 9
+    )
+  )
+}
+
+#' Paths of the files a built source consists of
+#'
+#' Internal.  Everything a source writes is prefixed with its name, so that a
+#' component can be rebuilt or removed on its own, and the \code{.gz.parquet}
+#' suffix records the codec.  The GADM files are the ones GNRS writes, and are
+#' read here if that package built them first.
+#' @keywords internal
+#' @noRd
+gvs_centroid_path <- function(level, dir = gvs_cache_dir()) {
+  file.path(dir, paste0("centroids-", level, ".gz.parquet"))
+}
+
+#' @keywords internal
+#' @noRd
+gvs_gadm_archive_path <- function(dir = gvs_cache_dir()) {
+  file.path(dir, paste0("gadm-", gvs_builtin_registry()$gadm$version, ".zip"))
+}
+
+#' @keywords internal
+#' @noRd
+gvs_gadm_path <- function(dir = gvs_cache_dir()) {
+  file.path(dir, "gadm-divisions.gz.parquet")
+}
+
+#' @keywords internal
+#' @noRd
+gvs_provenance_path <- function(source, dir = gvs_cache_dir()) {
+  file.path(dir, paste0(source, "-provenance.rds"))
+}
+
+#' The three levels the service resolves
+#' @keywords internal
+#' @noRd
+gvs_levels <- function() c("country", "state", "county")
+
+#' Has a source been built?
+#'
+#' Internal.  Built means the tables the resolver reads exist, not that a
+#' download completed: an interrupted build can leave one without the other.
+#' @keywords internal
+#' @noRd
+gvs_is_built <- function(source, dir = gvs_cache_dir()) {
+  switch(source,
+    gadm = file.exists(gvs_gadm_path(dir)),
+    centroids = all(file.exists(gvs_centroid_path(gvs_levels(), dir))),
+    FALSE
+  )
+}
+
+#' Files a source occupies in the cache
+#' @keywords internal
+#' @noRd
+gvs_source_files <- function(source, dir = gvs_cache_dir()) {
+  if (!dir.exists(dir)) {
+    return(character(0))
+  }
+  files <- list.files(dir, pattern = paste0("^", source, "-"), full.names = TRUE)
+  files[!dir.exists(files)]
+}
+
+#' How a set of source names would be written in a call
+#' @keywords internal
+#' @noRd
+gvs_source_arg <- function(sources) {
+  quoted <- paste0('"', sources, '"')
+  if (length(quoted) == 1) quoted else paste0("c(", paste(quoted, collapse = ", "), ")")
+}
+
+#' Report on the locally cached GVS reference data
+#'
+#' Shows each component of the local reference data, whether it has been built
+#' for offline use, and for those that have, which version it is and how much
+#' space it occupies.  A local result can be cited using the versions reported
+#' here.
+#'
+#' The cache is shared with the GNRS package, which derives its political
+#' division names from the same GADM release, so a GADM component built by
+#' either package is reported as built by both.
+#'
+#' @param dir Cache directory.  Defaults to the standard user cache location.
+#' @return A data.frame with one row per component.  \code{version} and
+#'   \code{downloaded} describe what is installed and are NA for a component
+#'   that has not been built.  \code{size_mb} is what the component occupies on
+#'   disk now; \code{download_mb} is what fetching it costs.
+#' @seealso \code{\link{GVS_local_build}}, \code{\link{GVS_local}}
+#' @export
+#' @examples {
+#'   status <- GVS_local_status()
+#' }
+GVS_local_status <- function(dir = gvs_cache_dir()) {
+  registry <- gvs_builtin_registry()
+  sources <- names(registry)
+
+  built <- vapply(sources, function(s) gvs_is_built(s, dir), logical(1))
+
+  provenance <- lapply(sources, function(s) {
+    path <- gvs_provenance_path(s, dir)
+    if (file.exists(path)) readRDS(path) else NULL
+  })
+  names(provenance) <- sources
+
+  from_record <- function(field, empty) {
+    vapply(sources, function(s) {
+      record <- provenance[[s]]
+      if (!built[[s]] || is.null(record)) {
+        return(empty)
+      }
+      value <- record[[field]]
+      if (is.null(value) || length(value) == 0) empty else as.character(value)[1]
+    }, empty)
+  }
+
+  size_mb <- vapply(
+    sources,
+    function(s) round(sum(file.size(gvs_source_files(s, dir))) / 1024^2, 1),
+    numeric(1)
+  )
+
+  out <- data.frame(
+    source = sources,
+    full_name = vapply(registry, function(x) x$full_name, character(1)),
+    built = unname(built),
+    version = unname(from_record("version", NA_character_)),
+    downloaded = unname(from_record("downloaded", NA_character_)),
+    size_mb = unname(size_mb),
+    download_mb = vapply(registry, function(x) as.numeric(x$download_mb), numeric(1)),
+    stringsAsFactors = FALSE,
+    row.names = NULL
+  )
+
+  absent <- out$source[!out$built]
+  if (length(absent) == length(sources)) {
+    message(
+      "No local reference data built yet in:\n  ", dir,
+      "\nRun GVS_local_build() to set it up."
+    )
+  } else if (length(absent) > 0) {
+    message(
+      "Not built: ", paste(absent, collapse = ", "),
+      ". Add with GVS_local_build(", gvs_source_arg(absent), ")."
+    )
+  }
+
+  out
+}
+
+#' Delete the locally cached GVS reference data
+#'
+#' Removes the reference data GVS derived, or one component of it.  It can be
+#' rebuilt at any time with \code{GVS_local_build()}.
+#'
+#' The GADM component is shared with the GNRS package: removing it leaves that
+#' package without its GADM layer too, and is refused unless
+#' \code{shared = TRUE}.
+#'
+#' @param dir Cache directory.  Defaults to the standard user cache location.
+#' @param sources NULL, the default, removes the centroid tables.  Otherwise the
+#'   components to remove.
+#' @param shared Allow removal of components shared with GNRS?
+#' @param ask Ask for confirmation before deleting? Defaults to TRUE in an
+#'   interactive session.
+#' @return TRUE if anything was removed, FALSE otherwise, invisibly.
+#' @export
+#' @examples \dontrun{
+#' GVS_local_remove()
+#' }
+GVS_local_remove <- function(dir = gvs_cache_dir(), sources = "centroids",
+                             shared = FALSE, ask = interactive()) {
+  if (!dir.exists(dir)) {
+    message("Nothing to remove; no cache directory at:\n  ", dir)
+    return(invisible(FALSE))
+  }
+  unknown <- setdiff(sources, names(gvs_builtin_registry()))
+  if (length(unknown) > 0) {
+    message("Unknown source(s): ", paste(unknown, collapse = ", "))
+    return(invisible(FALSE))
+  }
+  if ("gadm" %in% sources && !shared) {
+    message(
+      "The GADM component is shared with the GNRS package. ",
+      "Remove it with shared = TRUE, or with GNRS_local_remove(sources = \"gadm\")."
+    )
+    sources <- setdiff(sources, "gadm")
+    if (length(sources) == 0) {
+      return(invisible(FALSE))
+    }
+  }
+
+  files <- unlist(lapply(sources, gvs_source_files, dir = dir), use.names = FALSE)
+  size_mb <- round(sum(file.size(files), na.rm = TRUE) / 1024^2, 1)
+
+  if (ask) {
+    answer <- readline(paste0(
+      "Delete ", paste(sources, collapse = ", "), " (", size_mb, " MB) in\n  ",
+      dir, "\n? [y/N] "
+    ))
+    if (!tolower(trimws(answer)) %in% c("y", "yes")) {
+      message("Nothing removed.")
+      return(invisible(FALSE))
+    }
+  }
+
+  unlink(files)
+  for (s in sources) unlink(gvs_provenance_path(s, dir))
+  message("Removed ", size_mb, " MB from ", dir)
+  invisible(TRUE)
+}
+
+#' Default value for NULL
+#' @keywords internal
+#' @noRd
+`%||%` <- function(x, y) if (is.null(x)) y else x
diff --git a/R/local_index.R b/R/local_index.R
new file mode 100644
index 0000000..4894614
--- /dev/null
+++ b/R/local_index.R
@@ -0,0 +1,239 @@
+#' Paths of the location index
+#' @keywords internal
+#' @noRd
+gvs_index_path <- function(what, dir = gvs_cache_dir()) {
+  switch(what,
+    units = file.path(dir, "gadmindex-units.gz.parquet"),
+    raster = file.path(dir, "gadmindex-units-30s.tif"),
+    edges = file.path(dir, "gadmindex-edges-30s.tif"),
+    stop("unknown index file: ", what)
+  )
+}
+
+#' Build the raster location index for the current GADM release
+#'
+#' Internal.  Locating a point in its country, state/province and county/parish
+#' by exact point-in-polygon against GADM is too slow for large inputs.  This
+#' builds, once, a 30 arc-second raster of GADM level-2 divisions (level 1 or 0
+#' where a country has no level 2) and a mask of cells near a division boundary.
+#' A point in a cell well inside one division is located by raster lookup; a point
+#' in a boundary cell, or in a cell no division covers (coasts), is resolved
+#' exactly against the GeoPackage when one is available and flagged otherwise.
+#'
+#' Method, chosen to bound memory: GDAL burns each feature's unit number straight
+#' from the GeoPackage, the unit being numbered in SQL (\code{DENSE_RANK} over the
+#' unit key), so no features are loaded into R and no reclassification is
+#' needed.  Boundary cells are those whose (2 \code{edge_radius} + 1)-cell square
+#' neighbourhood holds more than one value (a unit, or no unit), computed in
+#' blocks of rows.  Because a cell takes the unit at its centre, a sliver of
+#' another division narrower than a cell can be missed; \code{edge_radius} widens
+#' the margin.
+#'
+#' Outputs: \code{gadmindex-units.gz.parquet} (unit -> gid_0, gid_1, gid_2 and
+#' names), \code{gadmindex-units-30s.tif} (INT32, 0 = no division) and
+#' \code{gadmindex-edges-30s.tif} (INT8, 1 = boundary cell).
+#' @keywords internal
+#' @noRd
+gvs_build_index <- function(gpkg, gadm_layer = "gadm_410", dir = gvs_cache_dir(create = TRUE),
+                            resolution = 1 / 120, edge_radius = 1, block_rows = 500, quiet = FALSE) {
+  for (pkg in c("sf", "terra", "nanoparquet")) {
+    if (!requireNamespace(pkg, quietly = TRUE)) stop("Building the index needs the '", pkg, "' package.", call. = FALSE)
+  }
+  t0 <- Sys.time()
+  el <- function() round(as.numeric(difftime(Sys.time(), t0, units = "mins")), 1)
+  blank <- function(col) sprintf("(%s IS NULL OR %s IN ('', 'NA'))", col, col)
+  key_sql <- sprintf("CASE WHEN NOT %s THEN GID_2 WHEN NOT %s THEN GID_1 ELSE GID_0 END",
+                     blank("GID_2"), blank("GID_1"))
+  unit_sql <- sprintf("DENSE_RANK() OVER (ORDER BY %s)", key_sql)
+
+  # ---- units: one per level-2 division (or the finest level a country has, to 2)
+  if (!quiet) message("Reading GADM attributes ...")
+  attrs <- sf::st_drop_geometry(sf::st_read(
+    gpkg, quiet = TRUE,
+    query = sprintf('SELECT %s AS unit, GID_0, COUNTRY, GID_1, NAME_1, GID_2, NAME_2 FROM "%s"', unit_sql, gadm_layer)
+  ))
+  isblank <- function(x) is.na(x) | !nzchar(x) | x == "NA"
+  first <- attrs[!duplicated(attrs$unit), , drop = FALSE]
+  first <- first[order(first$unit), , drop = FALSE]
+  stopifnot(identical(as.integer(first$unit), seq_len(nrow(first))))
+  units <- data.frame(
+    unit = as.integer(first$unit), gid_0 = first$GID_0, country = first$COUNTRY,
+    gid_1 = ifelse(isblank(first$GID_1), NA_character_, first$GID_1),
+    name_1 = ifelse(isblank(first$GID_1), NA_character_, first$NAME_1),
+    gid_2 = ifelse(isblank(first$GID_2), NA_character_, first$GID_2),
+    name_2 = ifelse(isblank(first$GID_2), NA_character_, first$NAME_2),
+    stringsAsFactors = FALSE
+  )
+  nanoparquet::write_parquet(units, gvs_index_path("units", dir), compression = "gzip")
+  if (!quiet) message("  ", nrow(attrs), " GADM features in ", nrow(units), " units (", el(), " min)")
+
+  # ---- burn unit numbers with GDAL, streaming from the GeoPackage
+  units_tif <- gvs_index_path("raster", dir)
+  unlink(units_tif)
+  if (!quiet) message("Rasterizing GADM units at ", round(resolution * 3600), " arc-seconds ...")
+  sf::gdal_utils("rasterize", gpkg, units_tif, options = c(
+    "-sql", sprintf('SELECT %s AS unit, geom FROM "%s"', unit_sql, gadm_layer),
+    "-a", "unit", "-a_nodata", "0", "-init", "0",
+    "-tr", format(resolution, digits = 15), format(resolution, digits = 15),
+    "-te", "-180", "-90", "180", "90", "-ot", "Int32",
+    "-co", "COMPRESS=DEFLATE", "-co", "TILED=YES", "-co", "BIGTIFF=YES"
+  ), quiet = TRUE)
+  if (!quiet) message("  done (", el(), " min)")
+
+  # ---- boundary cells
+  if (!quiet) message("Marking boundary cells ...")
+  gvs_mark_edges(units_tif, gvs_index_path("edges", dir), radius = edge_radius, block_rows = block_rows)
+  if (!quiet) message("  done (", el(), " min)")
+
+  provenance <- list(
+    source = "index", full_name = "GVS raster location index (GADM level 2 at 30 arc-seconds)",
+    version = gvs_builtin_registry()$gadm$version, downloaded = as.character(Sys.Date()),
+    resolution_arcsec = round(resolution * 3600), edge_radius = edge_radius,
+    n_features = nrow(attrs), n_units = nrow(units), built_minutes = el()
+  )
+  saveRDS(provenance, gvs_provenance_path("gadmindex", dir))
+  invisible(provenance)
+}
+
+#' Mark cells whose square neighbourhood holds more than one unit
+#'
+#' Internal.  Reads the unit raster in blocks of rows (with \code{radius} rows of
+#' overlap), takes running maxima and minima over the neighbourhood, and writes 1
+#' where they differ.  No-division cells count as a value (0), so coastal cells
+#' are boundary cells.  Columns do not wrap at the antimeridian.
+#' @keywords internal
+#' @noRd
+gvs_mark_edges <- function(units_tif, edges_tif, radius = 1, block_rows = 500) {
+  r <- terra::rast(units_tif)
+  nr <- terra::nrow(r)
+  nc <- terra::ncol(r)
+  shift_ext <- function(m, fun) {
+    # running fun over columns then rows, window 2 * radius + 1, edges replicated
+    out <- m
+    for (k in seq_len(radius)) {
+      left <- cbind(m[, -seq_len(k), drop = FALSE], m[, rep(nc, k), drop = FALSE])
+      right <- cbind(m[, rep(1L, k), drop = FALSE], m[, seq_len(nc - k), drop = FALSE])
+      out <- fun(out, left, right)
+    }
+    m2 <- out
+    nrm <- nrow(m2)
+    for (k in seq_len(radius)) {
+      up <- rbind(m2[-seq_len(k), , drop = FALSE], m2[rep(nrm, k), , drop = FALSE])
+      down <- rbind(m2[rep(1L, k), , drop = FALSE], m2[seq_len(nrm - k), , drop = FALSE])
+      out <- fun(out, up, down)
+    }
+    out
+  }
+  terra::readStart(r)
+  on.exit(terra::readStop(r), add = TRUE)
+  tmpl <- terra::rast(r)
+  unlink(edges_tif)
+  b <- terra::writeStart(tmpl, edges_tif, overwrite = TRUE,
+                         wopt = list(datatype = "INT1U", names = "edge",
+                                     gdal = c("COMPRESS=DEFLATE", "TILED=YES", "BIGTIFF=YES")))
+  for (s in seq(1L, nr, by = block_rows)) {
+    e <- min(nr, s + block_rows - 1L)
+    rs <- max(1L, s - radius)
+    re <- min(nr, e + radius)
+    v <- terra::readValues(r, row = rs, nrows = re - rs + 1L, col = 1L, ncols = nc, mat = FALSE)
+    v[is.na(v)] <- 0L
+    m <- matrix(as.integer(v), nrow = re - rs + 1L, ncol = nc, byrow = TRUE)
+    rm(v)
+    edge <- shift_ext(m, pmax) != shift_ext(m, pmin)
+    keep <- (s - rs + 1L):(e - rs + 1L)
+    terra::writeValues(tmpl, as.integer(t(edge[keep, , drop = FALSE])), s, length(keep))
+    rm(m, edge)
+    gc(verbose = FALSE)
+  }
+  terra::writeStop(tmpl)
+  invisible(edges_tif)
+}
+
+#' Locate points in current GADM divisions
+#'
+#' Internal.  Raster lookup, with exact point-in-polygon for points in boundary
+#' cells or in cells no division covers, when \code{gpkg} is given.  Points are
+#' processed as unique coordinates.
+#'
+#' @return data.frame with gid_0, gid_1, gid_2, country, name_1, name_2, and
+#'   \code{locate_method}: "raster", "exact", "exact: outside all divisions", or
+#'   "raster (boundary cell, not verified)" / "outside raster (not verified)" when
+#'   no GeoPackage was given.
+#' @keywords internal
+#' @noRd
+gvs_locate_current <- function(lon, lat, dir = gvs_cache_dir(), gpkg = NULL, gadm_layer = "gadm_410",
+                               tile_deg = 1) {
+  units <- as.data.frame(nanoparquet::read_parquet(gvs_index_path("units", dir)))
+  ru <- terra::rast(gvs_index_path("raster", dir))
+  re <- terra::rast(gvs_index_path("edges", dir))
+
+  dd0 <- gvs_dedupe(lon, lat)
+  ulon <- lon[dd0$first]
+  ulat <- lat[dd0$first]
+  ok <- is.finite(ulon) & is.finite(ulat) & abs(ulat) <= 90 & abs(ulon) <= 180
+  xy <- cbind(ulon[ok], ulat[ok])
+
+  unit <- rep(NA_integer_, length(ulon))
+  edge <- rep(NA, length(ulon))
+  unit[ok] <- as.integer(terra::extract(ru, xy)[, 1])
+  unit[ok & unit %in% 0L] <- NA_integer_
+  edge[ok] <- as.integer(terra::extract(re, xy)[, 1]) %in% 1L
+
+  method <- rep(NA_character_, length(ulon))
+  method[ok] <- "raster"
+  needs_exact <- ok & (is.na(unit) | edge)
+
+  if (any(needs_exact)) {
+    if (is.null(gpkg)) {
+      method[needs_exact & !is.na(unit)] <- "raster (boundary cell, not verified)"
+      method[needs_exact & is.na(unit)] <- "outside raster (not verified)"
+    } else {
+      idx <- which(needs_exact)
+      exact <- gvs_exact_units(ulon[idx], ulat[idx], units, gpkg, gadm_layer, tile_deg)
+      unit[idx] <- exact
+      method[idx] <- ifelse(is.na(exact), "exact: outside all divisions", "exact")
+    }
+  }
+
+  m <- dd0$group
+  un <- unit[m]
+  data.frame(
+    gid_0 = units$gid_0[un], country = units$country[un],
+    gid_1 = units$gid_1[un], name_1 = units$name_1[un],
+    gid_2 = units$gid_2[un], name_2 = units$name_2[un],
+    locate_method = method[m], stringsAsFactors = FALSE
+  )
+}
+
+#' Exact point-in-polygon against the GeoPackage, one tile at a time
+#' @keywords internal
+#' @noRd
+gvs_exact_units <- function(lon, lat, units, gpkg, gadm_layer, tile_deg = 1) {
+  old_s2 <- sf::sf_use_s2()
+  suppressMessages(sf::sf_use_s2(FALSE))
+  on.exit(suppressMessages(sf::sf_use_s2(old_s2)), add = TRUE)
+  ukey <- ifelse(!is.na(units$gid_2), units$gid_2, ifelse(!is.na(units$gid_1), units$gid_1, units$gid_0))
+  out <- rep(NA_integer_, length(lon))
+  tx <- floor(lon / tile_deg)
+  ty <- floor(lat / tile_deg)
+  tiles <- unique(data.frame(tx, ty))
+  for (t in seq_len(nrow(tiles))) {
+    i <- which(tx == tiles$tx[t] & ty == tiles$ty[t])
+    x0 <- tiles$tx[t] * tile_deg
+    y0 <- tiles$ty[t] * tile_deg
+    bb <- sf::st_as_sfc(sf::st_bbox(c(xmin = x0, ymin = y0, xmax = x0 + tile_deg, ymax = y0 + tile_deg),
+                                     crs = sf::st_crs(4326)))
+    g <- suppressMessages(suppressWarnings(sf::st_read(
+      gpkg, quiet = TRUE, wkt_filter = sf::st_as_text(bb),
+      query = sprintf('SELECT GID_0, GID_1, GID_2, geom FROM "%s"', gadm_layer)
+    )))
+    if (!nrow(g)) next
+    blank <- function(x) is.na(x) | !nzchar(x) | x == "NA"
+    gkey <- ifelse(!blank(g$GID_2), g$GID_2, ifelse(!blank(g$GID_1), g$GID_1, g$GID_0))
+    p <- sf::st_as_sf(data.frame(lon = lon[i], lat = lat[i]), coords = c("lon", "lat"), crs = 4326)
+    hit <- suppressMessages(sf::st_intersects(p, g))
+    first <- vapply(hit, function(h) if (length(h)) h[1] else NA_integer_, integer(1))
+    out[i] <- match(gkey[first], ukey)
+  }
+  out
+}
diff --git a/tests/testthat/test-local-gvs.R b/tests/testthat/test-local-gvs.R
new file mode 100644
index 0000000..2999d8b
--- /dev/null
+++ b/tests/testthat/test-local-gvs.R
@@ -0,0 +1,143 @@
+# Offline GVS: helpers that need no cache, then stages that need a built cache
+# (skipped unless the option GVS.test_cache names a directory holding one).
+
+test_that("gvs_dedupe groups identical rows, NAs equal", {
+  set.seed(2)
+  n <- 5000
+  a <- sample(c(1.5, 2.25, NA, 3), n, TRUE)
+  b <- sample(c("x", "y", NA), n, TRUE)
+  g <- GVS:::gvs_dedupe(a, b)
+  k <- paste(a, b)
+  expect_equal(length(g$first), length(unique(k)))
+  expect_true(all(k == k[g$first][g$group]))
+})
+
+test_that("Andoyer-Lambert distances agree with the geodesic", {
+  skip_if_not_installed("geosphere")
+  set.seed(3)
+  n <- 2000
+  lon1 <- runif(n, -180, 180)
+  lat1 <- runif(n, -85, 85)
+  p2 <- geosphere::destPoint(cbind(lon1, lat1), runif(n, 0, 360), runif(n, 0, 3e6))
+  g <- geosphere::distGeo(cbind(lon1, lat1), p2) / 1000
+  l <- GVS:::gvs_distance_km(lon1, lat1, p2[, 1], p2[, 2])
+  expect_lt(max(abs(l - g) / pmax(g, 1e-9)), 1e-5)
+  expect_equal(GVS:::gvs_distance_km(10, 10, 10, 10), 0)
+})
+
+test_that("record dates parse from years and ISO strings", {
+  d <- GVS:::gvs_parse_record_date(c("1985", "1991-12-01", "1991-12", "", NA, "c. 1900"))
+  expect_equal(d, as.Date(c("1985-07-01", "1991-12-01", "1991-12-01", NA, NA, NA)))
+})
+
+test_that("window pair search matches brute force", {
+  set.seed(4)
+  lon <- runif(3000, -180, 180)
+  lat <- runif(3000, -90, 90)
+  clon <- runif(40, -170, 170)
+  clat <- runif(40, -80, 80)
+  rlat <- runif(40, 0.5, 5)
+  rlon <- runif(40, 0.5, 8)
+  got <- GVS:::gvs_pairs_within(lon, lat, clon, clat, rlat, rlon)
+  bf <- which(outer(abs(lat), rep(1, 40)) >= 0 &
+                abs(outer(lat, clat, "-")) <= rep(rlat, each = 3000) &
+                abs(outer(lon, clon, "-")) <= rep(rlon, each = 3000), arr.ind = TRUE)
+  expect_setequal(paste(got$p, got$j), paste(bf[, 1], bf[, 2]))
+})
+
+test_that("level distances pick the nearest centroid and threshold comparisons are exact", {
+  skip_if_not_installed("geosphere")
+  tab <- data.frame(gid = "A")
+  for (i in 1:6) {
+    tab[[paste0("c", i, "_lon")]] <- 10 + i / 10
+    tab[[paste0("c", i, "_lat")]] <- 50
+    tab[[paste0("c", i, "_dmax")]] <- 2
+  }
+  r <- GVS:::gvs_level_distances(c(10.1, 10.62, 0), c(50, 50, 0), c("A", "A", "B"), tab, thr = 5)
+  expect_equal(r$type, c("std", "bb_main", NA))
+  expect_equal(r$dist[1], 0)
+  expect_equal(r$rel[2], 0.02 / 2, tolerance = 1e-9)
+  expect_true(is.na(r$dist[3]))
+  # a point exactly at the threshold distance is compared with the geodesic value
+  p <- geosphere::destPoint(c(10.1, 50), 180, 5000)
+  r2 <- GVS:::gvs_level_distances(p[1], p[2], "A", tab, thr = 5)
+  expect_equal(r2$dist, geosphere::distGeo(c(10.1, 50), p) / 1000)
+})
+
+cache <- getOption("GVS.test_cache", "")
+
+test_that("historical centroids and geovalidity (needs a built cache)", {
+  skip_if(!nzchar(cache) || !dir.exists(cache), "no test cache (option GVS.test_cache)")
+  skip_if_not_installed("sf")
+  v <- nanoparquet::read_parquet(file.path(cache, "cshapes-versions.gz.parquet"))
+  su <- v[v$gwcode == 365 & v$valid_from > as.Date("1950-01-01") & v$valid_to < as.Date("1992-01-01"), ][1, ]
+  hc <- GVS:::gvs_historical_centroids(c(su$c1_lon, 2.35), c(su$c1_lat, 48.86), NULL, cache,
+                                       69.99529, 0.01606376, 1)
+  expect_match(hc$version_id[1], "^cshapes:365:")
+  expect_true(is.na(hc$version_id[2]))
+
+  # GNRS-shaped input: USSR record in Vilnius; Czechoslovakia record in Dresden;
+  # Netherlands Antilles record on Curacao (no CShapes geometry: successors)
+  gn <- data.frame(gid_0 = c(NA, NA, NA), gid_1 = NA, gid_2 = NA,
+                   entity_key = c("SUHH", "CSHH", "ANHH"), is_historical = TRUE,
+                   successors = c("RU;LT", "CZ;SK", "CW;SX;BQ;AW"), subnational_status = NA)
+  loc <- data.frame(gid_0 = c("LTU", "DEU", "CUW"), gid_1 = NA, gid_2 = NA)
+  lon <- c(25.28, 13.74, -68.93)
+  lat <- c(54.69, 51.05, 12.12)
+  gv <- GVS:::gvs_geovalidity(lon, lat, NULL, loc, gn, cache, 1)
+  expect_equal(gv$geovalid, c(TRUE, FALSE, TRUE))
+  expect_match(gv$geovalid_basis[1], "^cshapes:365:")
+  expect_match(gv$geovalid_basis[3], "^successor:")
+  # at a date after the USSR, only under history = "all"
+  gd <- GVS:::gvs_geovalidity(lon, lat, as.Date(c("2005-07-01", NA, NA)), loc, gn, cache, 1)
+  expect_false(gd$geovalid[1])
+})
+
+test_that("centroid thresholds: absolute, relative, or both", {
+  d <- c(10, 50, 10, NA)
+  r <- c(0.5, 0.01, 0.01, 0.01)
+  expect_equal(GVS:::gvs_centroid_flag(d, r, a = 20), c(TRUE, FALSE, TRUE, FALSE))
+  expect_equal(GVS:::gvs_centroid_flag(d, r, r = 0.02), c(FALSE, TRUE, TRUE, TRUE))
+  expect_equal(GVS:::gvs_centroid_flag(d, r, a = 20, r = 0.02, combine = "or"), c(TRUE, TRUE, TRUE, TRUE))
+  expect_equal(GVS:::gvs_centroid_flag(d, r, a = 20, r = 0.02, combine = "and"), c(FALSE, FALSE, TRUE, FALSE))
+  expect_equal(GVS:::gvs_centroid_flag(d, r), rep(FALSE, 4))
+  expect_error(GVS:::gvs_check_thresholds(NULL, NULL), "thr_abs, thr_rel or both")
+  expect_error(GVS:::gvs_check_thresholds(c(province = 1), NULL), "named vector")
+  expect_error(GVS:::gvs_check_thresholds(c(1, 2), NULL), "named vector")
+  expect_true(GVS:::gvs_check_thresholds(NULL, c(county = 0.03)))
+})
+
+test_that("GVS_thresholds returns usable sets", {
+  expect_true(all(c("spec99", "tss_abs", "bien_default") %in% GVS_thresholds()))
+  s <- GVS_thresholds("spec99")
+  expect_equal(s$combine, "and")
+  expect_true(GVS:::gvs_check_thresholds(s$thr_abs, s$thr_rel))
+  expect_false("state" %in% names(s$thr_abs))
+  for (nm in GVS_thresholds()) {
+    x <- GVS_thresholds(nm)
+    expect_true(GVS:::gvs_check_thresholds(x$thr_abs, x$thr_rel))
+    expect_true(x$combine %in% c("or", "and"))
+  }
+})
+
+test_that("coordinate precision and inherent uncertainty follow the service", {
+  expect_equal(GVS:::gvs_decimal_places(c("12.50", "12.5", "12", "-9.1839", "0.0001", NA)),
+               c(1L, 1L, 0L, 4L, 4L, NA_integer_))
+  expect_equal(GVS:::gvs_decimal_places(c(12.5, 12, -9.1839)), c(1L, 0L, 4L))
+  # per-axis: latitude to 4 places, longitude to 2, at 9.18 S -> longitude cell dominates
+  u <- GVS:::gvs_coordinate_uncertainty_m(-9.1839, 4L, 2L)
+  expect_equal(round(u), 1099)
+  # whole degrees at 45 N
+  expect_equal(round(GVS:::gvs_coordinate_uncertainty_m(45, 0L, 0L)), 136339)
+  # more decimals means less uncertainty
+  expect_lt(GVS:::gvs_coordinate_uncertainty_m(45, 4L, 4L), GVS:::gvs_coordinate_uncertainty_m(45, 2L, 2L))
+})
+
+test_that("the default thresholds are the fitted spec995 set", {
+  f <- formals(GVS_local)
+  s <- GVS_thresholds("spec995")
+  expect_equal(eval(f$thr_abs), s$thr_abs)
+  expect_equal(eval(f$thr_rel), s$thr_rel)
+  expect_equal(eval(f$combine)[1], s$combine)
+  expect_equal(eval(f$max_uncertainty_m), 11000)
+})