Model Validation Demo (v0.5.0+)¶

TerraFlow ships spatial-block cross-validation (Roberts et al. 2017) as the primary validation surface. The Cohen's κ and Moran's I wrappers that lived in terraflow.validation through v0.4.x were removed in v0.5.0 — they didn't earn Methods-section citations (downstream papers cited sklearn / esda, not TerraFlow). This notebook shows:

How to invoke spatial-block CV via terraflow validate.
How to compute Cohen's κ via sklearn.metrics.cohen_kappa_score directly on the run's features.parquet.
How to compute Moran's I via esda.Moran on score residuals.

All three patterns are documented in Migration v0.4 → v0.5.

1. Build synthetic inputs + run pipeline¶

In [1]:

Copied!





import json, tempfile, textwrap
from pathlib import Path
import numpy as np
import pandas as pd
import rasterio
from rasterio.transform import from_origin
from terraflow.pipeline import run_pipeline
from terraflow.validation import run_validation

tmp = Path(tempfile.mkdtemp())
raster_path = tmp / 'raster.tif'
arr = np.arange(25, dtype='float32').reshape(5, 5)
with rasterio.open(
    raster_path, 'w', driver='GTiff', height=5, width=5, count=1,
    dtype='float32', crs='EPSG:4326',
    transform=from_origin(-100.0, 40.0, 0.01, 0.01),
) as ds:
    ds.write(arr, 1)

climate_path = tmp / 'climate.csv'
pd.DataFrame({
    'lat':        [40.0,  40.01,  40.02],
    'lon':        [-100.0, -99.99, -99.98],
    'mean_temp':  [18.0,   19.0,   20.0],
    'total_rain': [100.0, 120.0,  140.0],
}).to_csv(climate_path, index=False)

config_path = tmp / 'config.yml'
config_path.write_text(textwrap.dedent(f'''
    raster_path: {raster_path}
    climate_csv: {climate_path}
    output_dir: {tmp}/outputs
    roi: {{ type: bbox, xmin: -101.0, ymin: 39.0, xmax: -99.0, ymax: 41.0 }}
    model_params: {{ v_min: 0.0, v_max: 25.0, t_min: 0.0, t_max: 40.0,
                    r_min: 0.0, r_max: 300.0, w_v: 0.4, w_t: 0.3, w_r: 0.3 }}
    max_cells: 50
    validation:
      n_blocks_side: 2
      buffer_deg: 0.01
'''))

df = run_pipeline(config_path)
run_dir = Path(df.attrs['run_dir'])
print('run_dir:', run_dir)
import json, tempfile, textwrap
from pathlib import Path
import numpy as np
import pandas as pd
import rasterio
from rasterio.transform import from_origin
from terraflow.pipeline import run_pipeline
from terraflow.validation import run_validation

tmp = Path(tempfile.mkdtemp())
raster_path = tmp / 'raster.tif'
arr = np.arange(25, dtype='float32').reshape(5, 5)
with rasterio.open(
    raster_path, 'w', driver='GTiff', height=5, width=5, count=1,
    dtype='float32', crs='EPSG:4326',
    transform=from_origin(-100.0, 40.0, 0.01, 0.01),
) as ds:
    ds.write(arr, 1)

climate_path = tmp / 'climate.csv'
pd.DataFrame({
    'lat':        [40.0,  40.01,  40.02],
    'lon':        [-100.0, -99.99, -99.98],
    'mean_temp':  [18.0,   19.0,   20.0],
    'total_rain': [100.0, 120.0,  140.0],
}).to_csv(climate_path, index=False)

config_path = tmp / 'config.yml'
config_path.write_text(textwrap.dedent(f'''
    raster_path: {raster_path}
    climate_csv: {climate_path}
    output_dir: {tmp}/outputs
    roi: {{ type: bbox, xmin: -101.0, ymin: 39.0, xmax: -99.0, ymax: 41.0 }}
    model_params: {{ v_min: 0.0, v_max: 25.0, t_min: 0.0, t_max: 40.0,
                    r_min: 0.0, r_max: 300.0, w_v: 0.4, w_t: 0.3, w_r: 0.3 }}
    max_cells: 50
    validation:
      n_blocks_side: 2
      buffer_deg: 0.01
'''))

df = run_pipeline(config_path)
run_dir = Path(df.attrs['run_dir'])
print('run_dir:', run_dir)

INFO:terraflow:Loaded config from /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/config.yml

INFO:terraflow:Computed run fingerprint 2YOlX6x-yV8Ffp-5xpEk6_QhtwG-4-ydksUnEiA0KMc (config=b942a8b7ec0cc565d4bb30d577fa7af626929b13b7da2ccb1d7473ce3d1647b4, inputs=2)

INFO:terraflow:DataCatalog built: raster=raster.tif (CRS EPSG:4326, shape (5, 5)), climate=climate.csv (3 rows, vars=['mean_temp', 'total_rain'])

INFO:terraflow:Loaded raster from /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/raster.tif

INFO:terraflow:Loaded climate CSV from /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/climate.csv with 3 rows

INFO:terraflow:Climate variables: ['mean_temp', 'total_rain']

INFO:terraflow:Climate CSV validated successfully: 3 valid records

INFO:terraflow:Loaded raster: /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/raster.tif (CRS: EPSG:4326)

INFO:terraflow:Loaded climate data: /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/climate.csv

INFO:terraflow:Clipped raster to ROI

INFO:terraflow.climate:ClimateInterpolator initialised: strategy='spatial', interpolation_method='linear', variogram_mode='standard', records=3, variables=['mean_temp', 'total_rain']

INFO:terraflow:Initialized climate interpolator with strategy='spatial', method='linear'

INFO:terraflow:Sampled 25 cells from 25 valid cells in ROI

INFO:terraflow.climate:Linear interpolation failed (QH6154 Qhull precision error: Initial simplex is flat (facet 2 is coplanar with the interior point)

While executing:  | qhull d Qz Qc Qbb Q12 Qt
Options selected for Qhull 2020.2.r 2020/08/31:
  run-id 97641640  delaunay  Qz-infinity-point  Qcoplanar-keep  Qbbound-last
  Q12-allow-wide  Qtriangulate  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _max-width 0.02  Error-roundoff 1.4e-13  _one-merge 9.7e-13
  Visible-distance 2.8e-13  U-max-coplanar 2.8e-13  Width-outside 5.5e-13
  _wide-facet 1.7e-12  _maxoutside 1.1e-12

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p3(v4):    40 -1e+02 1e+02
- p1(v3):    40 -1e+02   0.1
- p2(v2):    40 -1e+02 -1.9e-14
- p0(v1):    40 -1e+02  0.21

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 1.4e-13.  The center point, facets and distances
to the center point are as follows:

center point    40.01   -99.99    25.08

facet p1 p2 p0 distance= -1.1e-08
facet p3 p2 p0 distance= -3.2e-15
facet p3 p1 p0 distance= -6.7e-15
facet p3 p1 p2 distance= -1.5e-14

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:        40     40.02  difference= 0.02
  1:      -100  -2.225e-308  difference=  100
  2:  -1.876e-14       100  difference=  100

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 1.4e-13.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.
), falling back to nearest-neighbor

INFO:terraflow.climate:Linear interpolation failed (QH6154 Qhull precision error: Initial simplex is flat (facet 2 is coplanar with the interior point)

While executing:  | qhull d Qz Qc Qbb Q12 Qt
Options selected for Qhull 2020.2.r 2020/08/31:
  run-id 97641640  delaunay  Qz-infinity-point  Qcoplanar-keep  Qbbound-last
  Q12-allow-wide  Qtriangulate  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _max-width 0.02  Error-roundoff 1.4e-13  _one-merge 9.7e-13
  Visible-distance 2.8e-13  U-max-coplanar 2.8e-13  Width-outside 5.5e-13
  _wide-facet 1.7e-12  _maxoutside 1.1e-12

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p3(v4):    40 -1e+02 1e+02
- p1(v3):    40 -1e+02   0.1
- p2(v2):    40 -1e+02 -1.9e-14
- p0(v1):    40 -1e+02  0.21

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 1.4e-13.  The center point, facets and distances
to the center point are as follows:

center point    40.01   -99.99    25.08

facet p1 p2 p0 distance= -1.1e-08
facet p3 p2 p0 distance= -3.2e-15
facet p3 p1 p0 distance= -6.7e-15
facet p3 p1 p2 distance= -1.5e-14

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:        40     40.02  difference= 0.02
  1:      -100  -2.225e-308  difference=  100
  2:  -1.876e-14       100  difference=  100

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 1.4e-13.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.
), falling back to nearest-neighbor

INFO:terraflow:Interpolated climate for 25 cells using strategy='spatial'

INFO:terraflow:Closed raster dataset

INFO:terraflow:Artifacts written to /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/outputs/runs/2YOlX6x-yV8Ffp-5xpEk6_QhtwG-4-ydksUnEiA0KMc (fingerprint=2YOlX6x-yV8Ffp-5xpEk6_QhtwG-4-ydksUnEiA0KMc, cells=25, total=0.05s)

run_dir: /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/outputs/runs/2YOlX6x-yV8Ffp-5xpEk6_QhtwG-4-ydksUnEiA0KMc

2. Spatial block cross-validation¶

Run via the public API. Results land in report.json under the validation key — block-level scores + buffer-aware fold construction (Roberts et al. 2017).

In [2]:

Copied!

report_path = run_validation(config_path)
report = json.loads(Path(report_path).read_text(encoding='utf-8'))
print(json.dumps(report.get('validation', {}), indent=2, default=str))
report_path = run_validation(config_path)
report = json.loads(Path(report_path).read_text(encoding='utf-8'))
print(json.dumps(report.get('validation', {}), indent=2, default=str))

INFO:terraflow:Running validation on /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/outputs/runs/2YOlX6x-yV8Ffp-5xpEk6_QhtwG-4-ydksUnEiA0KMc

INFO:terraflow:Validation block written to /var/folders/lw/h_v37nds21q4xf2_q24hhz8h0000gn/T/tmpm5xw9jkb/outputs/runs/2YOlX6x-yV8Ffp-5xpEk6_QhtwG-4-ydksUnEiA0KMc/report.json

{
  "method": "spatial_block_cv",
  "citation": "Roberts et al. 2017, Ecography",
  "n_blocks_side": 2,
  "buffer_deg": 0.01,
  "n_folds": 4,
  "mean_fold_accuracy": 0.8194444444444444,
  "kriging_loocv_rmse": null,
  "note": "model has no free parameters; fold accuracy reflects spatial label consistency, not fit generalization"
}

3. Cohen's κ — external (v0.5.0 pattern)¶

Removed from terraflow.validation in v0.5.0. Compute directly via sklearn.metrics.cohen_kappa_score against any reference label set. Below: a synthetic reference DataFrame matched to the run's features.parquet.

In [3]:

Copied!





from sklearn.metrics import cohen_kappa_score
features = pd.read_parquet(run_dir / 'features.parquet')
rng = np.random.default_rng(7)
ref_labels = rng.choice(['low','medium','high'], size=len(features))
kappa = cohen_kappa_score(ref_labels, features['label'].to_numpy())
print(f'cohen_kappa (synthetic reference): {kappa:.3f}')
from sklearn.metrics import cohen_kappa_score
features = pd.read_parquet(run_dir / 'features.parquet')
rng = np.random.default_rng(7)
ref_labels = rng.choice(['low','medium','high'], size=len(features))
kappa = cohen_kappa_score(ref_labels, features['label'].to_numpy())
print(f'cohen_kappa (synthetic reference): {kappa:.3f}')

cohen_kappa (synthetic reference): 0.062

4. Moran's I on residuals — external (v0.5.0 pattern)¶

Removed from terraflow.validation in v0.5.0. Spatial autocorrelation of residuals is the meaningful diagnostic — autocorrelation of the raw score reflects the suitability surface itself, not the model error.

Below: integer-code both the predicted labels (from features.parquet) and the synthetic reference set, then compute residuals = predicted − reference and pass those to esda.Moran with a k-NN spatial weights matrix.

In [4]:

Copied!





try:
    import esda, libpysal
    label_to_int = {'low': 0, 'medium': 1, 'high': 2}
    pred_int = features['label'].map(label_to_int).to_numpy()
    ref_int = pd.Series(ref_labels).map(label_to_int).to_numpy()
    residuals = pred_int - ref_int  # integer-coded prediction error per cell
    coords = features[['lon', 'lat']].to_numpy()
    w = libpysal.weights.KNN.from_array(coords, k=min(8, len(features) - 1))
    w.transform = 'r'
    mi = esda.Moran(residuals, w)
    print(f"Moran's I on residuals: {mi.I:.4f} (p_norm={mi.p_norm:.4f})")
except ImportError:
    print('esda + libpysal not installed. Install with: pip install esda libpysal')
try:
    import esda, libpysal
    label_to_int = {'low': 0, 'medium': 1, 'high': 2}
    pred_int = features['label'].map(label_to_int).to_numpy()
    ref_int = pd.Series(ref_labels).map(label_to_int).to_numpy()
    residuals = pred_int - ref_int  # integer-coded prediction error per cell
    coords = features[['lon', 'lat']].to_numpy()
    w = libpysal.weights.KNN.from_array(coords, k=min(8, len(features) - 1))
    w.transform = 'r'
    mi = esda.Moran(residuals, w)
    print(f"Moran's I on residuals: {mi.I:.4f} (p_norm={mi.p_norm:.4f})")
except ImportError:
    print('esda + libpysal not installed. Install with: pip install esda libpysal')

esda + libpysal not installed. Install with: pip install esda libpysal

Why these wrappers were cut¶

Both Cohen's κ and Moran's I are one-call patterns against features.parquet. Maintaining thin wrappers in terraflow.validation earned no Methods-section citations (downstream users cite sklearn / esda directly), so they were removed in v0.5.0 to keep the validation surface focused on spatial-block CV — the one diagnostic for which TerraFlow's per-cell schema + run fingerprinting genuinely add value. See issue #136 + PR #142.