GeoVersa logo

GeoVersa

Deep Learning + Geostatistics for Spatial Prediction

Overview · Install · DeepSCORPAN · Architecture · Mathematics · Auto-Configuration · Research-Family · Benchmark · Usage

Overview

GeoVersa is a research line on spatial prediction that combines machine learning with geostatistical residual correction, evaluated under the design-based paradigm of Wadoux et al. (2021).

It has two distinct faces, and this README documents both:

What it is Status
DB-TARF Design-Based Targeted Adaptive Residual Forest — the shipped, CRAN-published estimator on CRAN (v0.3.0)
ConvKrigingNet2D the deep-kriging research backbone (CNN + tabular + coordinate encoders with a differentiable anisotropic kriging layer) research code, in-repo

DB-TARF is the only model in this repository that beats the Wadoux reference RF on the design-based protocol; the deep backbone does not (see Benchmark). The deep-learning research family is preserved here in full, with its benchmark history.

The R package (DB-TARF)

The DB-TARF estimator plus its design-based conformal prediction intervals ship as the GeoVersa R package, available on CRAN:

install.packages("GeoVersa")

library(GeoVersa)

fit <- dbtarf_train_predict(
  train_df       = train,
  test_df        = test,
  response_name  = "y",
  predictor_names = c("elev", "slope", "ndvi"),
  coord_names    = c("x", "y")
)

fit$pred_test   # point predictions
fit$pi_q90      # design-based conformal 90% half-widths

The package deliberately contains only the DB-TARF core (R/dbtarf.R). The deep-learning research code (ConvKrigingNet2D and the torch-based learners) stays in code/ and is excluded from the package build — it is not CRAN-installable.

The deep-kriging research backbone

The ConvKrigingNet2D / GeoVersa configuration is the deep path documented throughout the architecture and mathematics sections below:

From SCORPAN to DeepSCORPAN

Digital Soil Mapping is grounded in the SCORPAN framework:

S(\mathbf{s}) = f(S, C, O, R, P, A, N) + \varepsilon(\mathbf{s})

where $S, C, O, R, P, A, N$ represent the classical soil-forming factors and $\varepsilon(\mathbf{s})$ is the spatially structured residual.

Classical regression-kriging keeps these two parts separate:

\hat{S}(\mathbf{s}_0) =
\hat{f}(\mathbf{x}_0) +
\sum_{i=1}^{n} \lambda_i(\mathbf{s}_0)\,\hat{\varepsilon}(\mathbf{s}_i)

GeoVersa reframes that idea as a DeepSCORPAN program: the trend model and the residual geostatistical correction are learned inside one trainable system.

For the active GeoVersa benchmark path, the final standardized prediction is:

\hat{y}_i^{(s)} =
\hat{y}_{i}^{\mathrm{base}} + \beta\,\delta_i

with:

This is the current concrete realization of the broader DeepSCORPAN idea in the repository.

SCORPAN Mapping in GeoVersa

SCORPAN factor Current GeoVersa realization
$S, C, O, R, P, A$ Tabular encoder over point-level covariates
$O, R$ as local landscape texture 2D CNN patch encoder
$N$ Coordinate encoder
$\varepsilon(\mathbf{s})$ Differentiable anisotropic residual-kriging layer

Active Architecture

The benchmarked architecture has four coupled components:

Component Input Output Role
Tabular encoder point covariates $\mathbf{x}_i$ $\mathbf{e}_i^{\mathrm{tab}}$ nonlinear trend from SCORPAN-style attributes
Patch encoder raster patch $\mathbf{P}_i$ $\mathbf{e}_i^{\mathrm{patch}}$ local terrain and remote-sensing texture
Coordinate encoder location $\mathbf{s}_i$ $\mathbf{e}_i^{\mathrm{coord}}$ smooth spatial trend component
Residual kriging layer neighbor residual bank $\delta_i$ anisotropic spatial correction 
flowchart LR
    X["Tabular covariates x_i"] --> TAB["Tabular MLP"]
    P["Raster patch P_i"] --> PATCH["Patch CNN + projection"]
    C["Coordinates s_i"] --> COORD["Coordinate MLP + projection"]

    TAB --> FUSE["Fusion block"]
    PATCH --> FUSE
    COORD --> FUSE

    FUSE --> BASE["Scalar head: yhat_base"]
    BASE --> BANK["Residual bank refresh"]
    BANK --> KRIG["Anisotropic residual kriging"]
    C --> KRIG

    BASE --> OUT["Final predictor"]
    KRIG --> OUT
    OUT["yhat = yhat_base + beta delta"]

Important implementation note

The active kriging layer is purely spatial in its weighting rule. It stores latent vectors in the bank for compatibility with the broader GeoVersa research program, but the current benchmarked path does not use latent similarity inside the kriging weights. That distinction matters and is reflected below in the mathematics.

Mathematical Formulation

Notation

Symbol Meaning
$\mathbf{x}_i \in \mathbb{R}^p$ tabular covariates
$\mathbf{P}_i \in \mathbb{R}^{C \times H \times W}$ raster patch
$\mathbf{s}_i \in \mathbb{R}^2$ coordinates
y_i observed target
T(y_i) optional target transform
y_i^{(s)} standardized training target

The standardized target is:

y_i^{(s)} = \frac{T(y_i) - \mu_y}{\sigma_y}

Encoder Fusion

Each modality is encoded independently and projected into a shared latent space:

\begin{aligned}
\mathbf{e}_i^{\mathrm{tab}}   &= f_{\mathrm{tab}}(\mathbf{x}_i), \\
\mathbf{e}_i^{\mathrm{patch}} &= W_{\mathrm{patch}}\,f_{\mathrm{cnn}}(\mathbf{P}_i), \\
\mathbf{e}_i^{\mathrm{coord}} &= W_{\mathrm{coord}}\,f_{\mathrm{coord}}(\mathbf{s}_i).
\end{aligned}

These embeddings are fused into a joint representation:

\mathbf{z}_i =
f_{\mathrm{fuse}}
\left[
\mathbf{e}_i^{\mathrm{tab}},
\mathbf{e}_i^{\mathrm{patch}},
\mathbf{e}_i^{\mathrm{coord}}
\right]

and mapped to the base trend:

\hat{y}_i^{\mathrm{base}} = h(\mathbf{z}_i)

Residual Memory Bank

At each bank refresh, the current backbone defines residuals over the training set:

r_j = y_j^{(s)} - \hat{y}_j^{\mathrm{base}}

so the memory bank can be written as:

\mathcal{B} = \left\{ \left(\mathbf{z}_j, \mathbf{s}_j, r_j\right) \right\}_{j=1}^{n_{\mathrm{train}}}

Anisotropic Residual Kriging

For query point $i$ and neighbor $j$, GeoVersa rotates coordinate offsets by a learned anisotropy angle $\theta$:

\begin{aligned}
\Delta x_{ij} &= x_i - x_j, \\
\Delta y_{ij} &= y_i - y_j, \\
u_{ij} &= \cos(\theta)\,\Delta x_{ij} + \sin(\theta)\,\Delta y_{ij}, \\
v_{ij} &= -\sin(\theta)\,\Delta x_{ij} + \cos(\theta)\,\Delta y_{ij}.
\end{aligned}

The anisotropic distance is:

d_{ij}^{\mathrm{aniso}} =
\sqrt{
\left(\frac{u_{ij}}{\ell_{\mathrm{maj}}}\right)^2 +
\left(\frac{v_{ij}}{\ell_{\mathrm{min}}}\right)^2 + \varepsilon
}

The active weighting rule is the exponential covariance softmax used in AnisotropicExpCovKrigingLayer_Auto:

w_{ij} =
\frac{\exp\!\left(-3\,d_{ij}^{\mathrm{aniso}}\right)}
{\sum_{k \in \mathcal{N}(i)} \exp\!\left(-3\,d_{ik}^{\mathrm{aniso}}\right)}

and the residual correction is:

\delta_i = \sum_{j \in \mathcal{N}(i)} w_{ij}\,r_j

This is intentionally different from earlier latent-attention variants explored in the GeoVersa research line. The benchmarked path described here uses pure anisotropic spatial weighting.

Final Predictor

The global gate is learned end to end:

\beta = \sigma(\mathrm{logit}_{\beta})

and the final standardized prediction is:

\hat{y}_i^{(s)} = \hat{y}_i^{\mathrm{base}} + \beta\,\delta_i

Training Objective

During warmup, only the backbone is optimized:

\mathcal{L}_{\mathrm{warmup}} =
\mathrm{Huber}\!\left(y^{(s)}, \hat{y}^{\mathrm{base}}\right)

After warmup, the full model is trained with:

\begin{aligned}
\mathcal{L}_{\mathrm{full}}
&=
\mathrm{Huber}\!\left(y^{(s)}, \hat{y}\right)
 + \lambda_{\mathrm{base}}\,
\mathrm{Huber}\!\left(y^{(s)}, \hat{y}^{\mathrm{base}}\right) \\
&\quad
 + \alpha_{\mathrm{ME}}
\left(
\mathrm{mean}\!\left(\hat{y}^{\mathrm{base}}\right) -
\mathrm{mean}\!\left(y^{(s)}\right)
\right)^2 \\
&\quad
 + \lambda_{\mathrm{cov}}
\left(
\frac{\mathrm{sd}\!\left(\hat{y}^{\mathrm{base}}\right)}
{\mathrm{sd}\!\left(y^{(s)}\right) + \varepsilon}
 - 1
\right)^2
\end{aligned}

For the current standalone GeoVersa benchmark:

Automatic Configuration

The active auto-config logic lives in code/ConvKrigingNet2D_Auto_v5.R and derives the model from the fold geometry, the fitted variogram, and device constraints.

Variogram-Derived Spatial Initialization

From the empirical variogram, the code initializes:

The active rules are:

\begin{aligned}
K &= \mathrm{clamp}
\left(
\mathrm{round}
\left(
\frac{n_{\mathrm{train}} \pi \,\mathrm{range}_{\mathrm{major}}^2}
{\mathrm{area}}
\right),
6, 30
\right), \\
\mathrm{logit}_{\beta,0} &= 2 - 6r.
\end{aligned}

Capacity Rules

The latent widths and patch geometry follow:

\begin{aligned}
d &=
\mathrm{clamp}
\left(
64 \left\lceil \frac{\sqrt{n_{\mathrm{train}}}}{8} \right\rceil,
128, 512
\right), \\
\mathrm{patch\_size} &=
\mathrm{clamp}
\left(
\left\lfloor \sqrt{n_{\mathrm{train}}} \right\rfloor,
8, 31
\right), \\
\mathrm{patch\_dim} &=
\mathrm{clamp}
\left(
\left\lceil \sqrt{C H W} \right\rceil,
\frac{d}{4}, d
\right), \\
\mathrm{coord\_dim} &=
\mathrm{clamp}
\left(
32 + 24(1 - \rho_{\mathrm{aniso}}) + 8(1-r),
32, 64
\right),
\end{aligned}

where:

\rho_{\mathrm{aniso}} = \frac{\ell_{\mathrm{min}}}{\ell_{\mathrm{maj}}}

If patches are already cached, the actual cached tensor size overrides the heuristic patch_size.

Loss Weights And Warmup

The current v5 rules are:

\begin{aligned}
\lambda_{\mathrm{base}} &= \max(0.05,\; 0.10r), \\
\alpha_{\mathrm{ME}} &= 0.375r, \\
\lambda_{\mathrm{cov}} &= 0.025(1-r), \\
\mathrm{max\_warmup} &= \mathrm{clamp}\!\left(\mathrm{round}(4 + 16r), 4, 20\right).
\end{aligned}

The $0.05$ floor on $\lambda_{\mathrm{base}}$ is active because the standalone benchmark needs a minimum direct backbone signal in strongly spatial folds. Without that floor, the backbone under-learns.

Optimization Rules

The initial learning rate is estimated with a Polyak-style probe:

\alpha_{\mathrm{init}} =
\mathrm{clamp}
\left(
0.01 \frac{\mathcal{L}}{\lVert \nabla \mathcal{L} \rVert_2^2},
10^{-5},
10^{-3}
\right)

Weight decay scales with model size:

\mathrm{wd} =
\mathrm{clamp}
\left(
\frac{10^{-3}}{\sqrt{(n_{\mathrm{params}}/10^6)/5}},
10^{-4},
10^{-2}
\right)

Batch size is chosen from the smaller of:

After warmup, the code adapts:

GeoVersa Research Family

The repository is not only the current ConvKrigingNet2D benchmark. It also records a broader family of GeoVersa research directions through archived searches, confirmations, and benchmark folders under results/.

Benchmarked Core

Functional And Basis-Oriented Families

Attention, Memory, And Set-Based Families

Structured Expert And Hybrid Families

Explicit Spatial-Structure Families

These families should be understood as part of the GeoVersa research program. The current benchmarked and documented production path, however, is the pure ConvKrigingNet2D / GeoVersa route described in this README.

Wadoux Reference And Current Results

The repository now keeps two separate comparison tracks:

1. GeoVersa Benchmark

code/run_wadoux_style_rf_conv_comparison.R

This runs GeoVersa itself on Wadoux-style validation splits.

2. Wadoux RF Reference Reproduction

code/run_wadoux_rf_reference.R

This reproduces the Random Forest reference under the same validation framework.

Tracked upstream reference:

Role Of Each Validation Protocol

The validation suite is not a collection of interchangeable scores. Each protocol estimates a different empirical risk functional, so a scientifically correct reading must respect the estimand behind the split design, in the sense emphasized by Wadoux et al. (2021).

Let $\ell(y,\hat{y})$ denote a prediction loss and let $\mathcal{D}$ be the full prediction domain.

For the optional full-domain benchmark, the population-style target is:

\widehat{\mathcal{R}}_{\mathrm{Population}}
=
\frac{1}{|\mathcal{D}|}
\sum_{\mathbf{s} \in \mathcal{D}}
\ell\!\left(y(\mathbf{s}), \hat{y}(\mathbf{s})\right)

The DesignBased protocol instead estimates map error from an independent simple-random test sample $\mathcal{S}_{\mathrm{SRS}} \subset \mathcal{D}$:

\widehat{\mathcal{R}}_{\mathrm{Design}}
=
\frac{1}{|\mathcal{S}_{\mathrm{SRS}}|}
\sum_{\mathbf{s} \in \mathcal{S}_{\mathrm{SRS}}}
\ell\!\left(y(\mathbf{s}), \hat{y}(\mathbf{s})\right)

This is why DesignBased is the main benchmark in this repository whenever the scientific claim is map accuracy rather than only local interpolation skill.

For RandomKFold, the repository evaluates the out-of-fold risk under random partitioning:

\widehat{\mathcal{R}}_{\mathrm{RandomKFold}}
=
\frac{1}{K}
\sum_{k=1}^{K}
\frac{1}{|F_k|}
\sum_{i \in F_k}
\ell\!\left(y_i, \hat{y}_i^{(-k)}\right)

For SpatialKFold, the same out-of-fold estimator is applied after replacing random folds by spatially separated blocks $B_k$:

\widehat{\mathcal{R}}_{\mathrm{SpatialKFold}}
=
\frac{1}{K}
\sum_{k=1}^{K}
\frac{1}{|B_k|}
\sum_{i \in B_k}
\ell\!\left(y_i, \hat{y}_i^{(-k)}\right)

For BLOOCV, each focal test location is evaluated after excluding a buffer of radius $h$ around that location from calibration:

\widehat{\mathcal{R}}_{\mathrm{BLOO}}(h)
=
\frac{1}{M}
\sum_{m=1}^{M}
\ell\!\left(
y(\mathbf{s}_m),
\hat{y}^{(-\mathcal{N}_h(\mathbf{s}_m))}(\mathbf{s}_m)
\right)

flowchart LR
    A["Population"] --> A1["Full-domain empirical risk"]
    B["DesignBased"] --> B1["Independent SRS estimate of map accuracy"]
    C["RandomKFold"] --> C1["Random out-of-fold predictive risk"]
    D["SpatialKFold"] --> D1["Blocked spatial generalization risk"]
    E["BLOOCV"] --> E1["Buffered local extrapolation/interpolation stress test"]

The correct interpretation is therefore:

GeoVersa's thesis is not weakened by this distinction. On the contrary: the model is designed to combine a deep covariate-driven backbone with spatial residual correction, so it should be read as a system with two complementary ambitions:

DB-TARF: the only model here that beats the RF reference

DB-TARF (Design-Based Targeted Adaptive Residual Forest) predicts RF_selected + lambda * correction, choosing an RF profile and a cross-fitted residual correction that is accepted only when it improves the out-of-fold RMSE. It is the estimator published on CRAN.

The canonical benchmark below is the v0.3.0 reference campaign: the exact code published on CRAN (R/dbtarf.R), runner defaults (ensemble_top_k = 3, n_folds = 5), n = 500, DesignBased, seed 123, and 30 repetitions in every one of the three Wadoux sampling designs. RMSE and ME are in Mg ha⁻¹.

Design RF (Wadoux default) RF_tuned Cubist DB-TARF ΔRMSE Paired wins Wilcoxon $p$
random 32.49 32.48 33.26 31.60 −0.895 28 / 30 6.5×10⁻⁹
regular_grid 32.28 32.27 33.17 31.38 −0.901 30 / 30 9.3×10⁻¹⁰
clustered_random 35.78 35.60 37.16 34.57 −1.205 28 / 30 1.8×10⁻⁸

95 % CIs for ΔRMSE: random [−1.08, −0.71], regular_grid [−1.05, −0.75], clustered_random [−1.47, −0.94].

Two readings matter here:

The gain is largest under the hardest clustered design (−1.21), the regime where a calibration sample least resembles the population. The reference is the untuned default RF, as in Wadoux et al. (2021); DB-TARF performs all of its selection on the calibration sample only.

The advantage is not an artifact of the estimand

DB-TARF selects its correction by an out-of-fold criterion and is scored on the design-based protocol it is named for, so a fair objection is that the gain might be specific to that yardstick. Repeating the comparison under three further protocols (same configuration, n = 500, seed 123, 10 repetitions each) shows it is not:

Estimand $n$ RF DB-TARF ΔRMSE Wins Wilcoxon $p$
Population (SRS 20 000) 10 32.18 31.51 −0.669 10/10 0.00098
DesignBased (SRS 500) 30 32.49 31.60 −0.895 28/30 6.5×10⁻⁹
Random 10-fold CV 10 32.24 31.12 −1.125 10/10 0.00098
Spatial block CV (350 km) 10 34.72 33.45 −1.275 10/10 0.00098

DB-TARF wins in all 10 draws under every one of the three additional protocols, and mean MEC rises under each.

Four honest qualifications:

The point estimates do order −0.67 → −0.90 → −1.13 → −1.28, but adjacent protocols are not statistically distinguishable (random vs spatial CV, p = 0.36) and this is not an ordering in spatial stress — random 10-fold CV is the least spatially demanding of the four. BLOOCV is excluded by cost.

Full method, honest caveats, reproduction commands and paired statistics: docs/dbtarf_wadoux_validation.md.

The deep backbone does not beat RF

For contrast, the standalone ConvKrigingNet2D / GeoVersa deep path on the same protocol (random, DesignBased, n = 500, but only 3 repetitions, auto profile, mps):

Model ME RMSE Pearson² MEC
RF benchmark 0.187 32.313 0.883 0.877
GeoVersa, standalone auto 0.403 33.033 0.873 0.870

The deep backbone does not beat the RF benchmark on DesignBased (ΔRMSE = +0.720). This is an honest negative result and is why the CRAN package ships DB-TARF rather than the deep model. Note this run uses 3 repetitions against DB-TARF's 30, so its RF baseline (32.313) is a noisier draw than the campaign's (32.493); the two RF columns are not directly comparable across tables.

The original saved Wadoux .Rdata outputs are still not present in the mirrored upstream checkout. Their availability is tracked in docs/wadoux2021-reference/official_rdata_manifest.csv.

Experimental DAD-RK Path

DAD-RK means Design-Aware Density-Reweighted Kriging. It keeps the active GeoVersa Auto v5 network but changes the training objective so the loss approximates a design-based map-accuracy target more closely:

\mathcal{L}_{\mathrm{DAD}}
=
\frac{\sum_i w_i\,\ell(y_i,\hat y_i)}{\sum_i w_i},
\qquad
w_i \propto
\left(d_i^{\mathrm{spatial}}\right)^{1-\rho}
\left(d_i^{\mathrm{environmental}}\right)^\rho

where the two distances are local k-nearest-neighbour spacing proxies in coordinate space and covariate space. Dense clusters get lower weight; isolated spatial/environmental support gets higher weight. The implementation is opt-in through:

Sys.setenv(
  GEOVERSA_AUTO_DESIGN_WEIGHTS = "1",
  GEOVERSA_AUTO_DESIGN_WEIGHT_MODE = "hybrid"
)

The ready-to-run entrypoint is code/Benchmark_DADRK.R.

Auto v5 also includes an internal module meta-controller (GEOVERSA_AUTO_META_CONTROLLER=1 by default). It decides whether to activate optional branches such as quantile basis, cross tab-patch fusion, DAM, VRM/CVRM, PedoMoE, NKO, local gates, local linear kriging, multiscale kriging, and DAD-RK from fold-level signals: sample size, nugget ratio, anisotropy, density imbalance, covariate redundancy, patch/tabular capacity ratio, and baseline residual variance. Explicit environment flags still act as manual overrides. The pure meta-controller benchmark entrypoint is code/Benchmark_GeoVersaMeta.R.

How To Run

GeoVersa benchmark

Sys.setenv(
  WADOUX_AUTO_V5_SCRIPT = normalizePath(file.path(getwd(), "code", "ConvKrigingNet2D_Auto_v5.R"), mustWork = FALSE),
  WADOUX_AUTO_SCRIPT    = normalizePath(file.path(getwd(), "code", "ConvKrigingNet2D_Auto.R"), mustWork = FALSE),
  WADOUX_MODELS         = "ConvKrigingNet2D",
  WADOUX_PROTOCOLS      = "DesignBased",
  WADOUX_N_ITER         = "10",
  WADOUX_SAMPLE_SIZE    = "500",
  WADOUX_MODEL_PROFILE  = "auto",
  WADOUX_DEVICE         = "mps",
  WADOUX_RESULTS_DIR    = "results/geoversa_run"
)

source("code/run_wadoux_style_rf_conv_comparison.R")

Useful override knobs exposed by the runner:

Wadoux RF reference reproduction

Sys.setenv(
  WADOUX_RF_SCENARIO    = "random",
  WADOUX_RF_PROTOCOLS   = "Population,DesignBased",
  WADOUX_RF_N_ITER      = "10",
  WADOUX_RF_SAMPLE_SIZE = "500",
  WADOUX_R2_METHOD      = "pearson",
  WADOUX_RF_RESULTS_DIR = "results/wadoux_rf_reference"
)

source("code/run_wadoux_rf_reference.R")

Import official Wadoux .Rdata outputs

source("code/import_wadoux_official_rdata.R")

DB-TARF DesignBased benchmark

Sys.setenv(
  DBT_MODELS      = "RF,Cubist,DBTARF",
  DBT_PROTOCOLS   = "DesignBased",
  DBT_N_ITER      = "30",
  DBT_SAMPLE_SIZE = "500",
  DBT_TRAIN_SEED  = "123",
  DBT_RESULTS_DIR = "results/canon_v030_full_random_confirm30"
)

source("code/run_dbtarf_wadoux_comparison.R")

Repository Layout

code/
  ConvKrigingNet2D_Auto.R
  ConvKrigingNet2D_Auto_v5.R
  run_wadoux_style_rf_conv_comparison.R
  run_wadoux_rf_reference.R
  import_wadoux_official_rdata.R
  wadoux2021_rf_reproduction_helpers.R

docs/
  wadoux2021-reference/

logo/
  GeoVersa Logo.png

results/
  geoversa_blw_confirm10_20260406/
  wadoux2021_rf_reference_random_designbased_10iter_pearson_20260405/
  many archived GeoVersa family searches and confirmations

Citation

@software{GeoVersa,
  author = {Rodrigues, Hugo},
  title  = {{GeoVersa}: Deep Learning + Geostatistics for Spatial Prediction},
  year   = {2026},
  url    = {https://github.com/HugoMachadoRodrigues/GeoVersa}
}

References