Skip to content

0.3 roadmap

Jump to bottom
Ivan Kolomiets edited this page May 15, 2026 · 1 revision

0.3 Roadmap — uncertainty-aware soil analysis

Status: planning. This is a discussion document, not a commitment. Iterate freely.

Headline

0.3 makes uncertainty a first-class citizen across pedotri's surface.

The motivation: SoilGrids (and most modern digital soil maps) publish Q0.05 / Q0.50 / Q0.95 quantile layers alongside the mean. Treating those rasters as deterministic — classifying the mean, averaging the mean over a region — throws away the most scientifically valuable thing the upstream model produced. Two concrete workflows are crippled by this:

  1. Classification near class boundaries. A point at (clay=27, sand=45) sits comfortably in USDA clay loam if you classify the mean — but if clay's 90% CI is [18, 36] and sand's is [35, 55], the uncertainty box overlaps four classes. The honest answer is a distribution over classes.
  2. Regional aggregation for stock / inventory work. When you only have a coarse AOI (a village, a commune, a farm address) and need a number for SOC stock, the right output is a posterior distribution over the regional mean, masked to the relevant land use — not a single point pulled from village-centroid data.

0.3 addresses both with a small, coherent set of additions.

Out of scope (and why)

  • Sentinel imagery, crop classification, productivity zoning. Real product, wrong library. Belongs in a companion package built on top of pedotri (+ sentinelhub-py, pystac-client, scikit-learn). Forcing this into core would explode the dep tree and break pedotri's mental model.
  • Climate data. Different shape (time series), different sources (ERA5, WorldClim), mature libraries already exist. Companion package territory.
  • DEM-derived topographic indices. Pedology-adjacent and a good fit eventually — pushed to 0.4 to keep 0.3 focused.

Module additions

src/pedotri/
├── uncertainty.py        # distribution helpers           (NEW)
├── sources/                                               (NEW)
│   ├── soilgrids.py      # mean + Q05/Q50/Q95 stack
│   └── worldcover.py     # optional dep, mask source
├── zonal.py              # regional aggregation           (NEW)
├── classifier.py         # extend with uncertainty kwargs
└── raster.py             # classify_geotiff with uncertainty + mask

Functional style throughout, consistent with existing raster.py / interp.py — (array, profile) tuples flow between primitives. No RasterLayer base class yet; add only if ≥3 places in 0.3 itself earn it.

Three data flows

1. Point classification with uncertainty

(sand_mean, clay_mean) + (sand_q05, sand_q95, clay_q05, clay_q95)
        |
        v
fit marginal distributions      (uncertainty.fit_from_quantiles)
        |
        v
project onto compositional simplex   (sand + silt + clay = 100)
        |
   +----+----+
   v         v
"distance"   "monte_carlo"
classify     sample N compositions,
mean, cmp    classify each, tally
|dist| vs σ  class frequencies
        |
        v
ClassificationResult { key, probabilities, entropy, confidence }

2. Raster classification with uncertainty + mask

soil_stack (mean + Q05 + Q95 per axis as GeoTIFFs)
        |
        v
align grids (existing _check_rasters_aligned)
        |
        v
[optional] mask: vector polygon OR raster (e.g. WorldCover)
                 + class allowlist (mask_include=[...])
        |
        v
chunk pixels (windowed read — no full materialization)
        |
        v
vectorized uncertainty-aware classify per chunk
        |
        v
write: class_raster.tif + confidence_raster.tif
       (+ optional probability stack)

3. Regional aggregation — the SOC stock case

inputs:
  region:     AOI polygon (village, commune, field)
  properties: { "soc": (mean, q05, q95),
                "bd":  (mean, q05, q95), ... }
  mask:       worldcover raster OR user raster/vector
  include:    [worldcover.CROPLAND]
        |
        v
rasterize region polygon              -> region_mask
intersect with land-cover mask        -> effective_mask
        |
        v
for n in 1..N_samples:
   draw realization of each pixel from its (Q05, Q50, Q95) dist
   compute weighted regional mean over effective_mask
        |
        v
AggregateDistribution {
   mean, std, q05, q50, q95, samples,
   n_pixels_used, mask_coverage
}
        |
        v
[optional] combine(): propagate through user formula
   stock = soc * bd * depth * area  ->  distribution of stock

Proposed API surface

# 1. Point — backward-compatible extension
result = pedotri.classify(
    sand=27, clay=45,
    sand_uncertainty=(20, 34),     # (q05, q95) or stddev
    clay_uncertainty=(38, 52),
    method="distance",             # or "monte_carlo"
    detailed=True,
)
result.probabilities   # {'clay_loam': 0.62, 'sandy_clay_loam': 0.28, ...}
result.entropy         # in [0, log(n_classes)]

# 2. Raster — extends classify_geotiff
pedotri.classify_geotiff(
    sand=("sand_mean.tif", "sand_q05.tif", "sand_q95.tif"),
    clay=("clay_mean.tif", "clay_q05.tif", "clay_q95.tif"),
    classification="USDA",
    mask="worldcover_2021.tif",
    mask_include=[40],             # WorldCover cropland
    out_class="class.tif",
    out_confidence="confidence.tif",
)

# 3. Regional aggregation — new
agg = pedotri.zonal_aggregate(
    region=village_polygon,
    properties={
        "soc": ("soc_mean.tif", "soc_q05.tif", "soc_q95.tif"),
        "bd":  ("bd_mean.tif",  "bd_q05.tif",  "bd_q95.tif"),
    },
    mask="worldcover_2021.tif",
    mask_include=[40],
    n_samples=1000,
)
agg["soc"].q05, agg["soc"].mean, agg["soc"].q95
agg["soc"].mask_coverage          # fraction of region pixels kept

# Propagate through a user formula
stock = agg.combine(
    lambda soc, bd, depth=0.3, area=village_area:
        soc * bd * depth * area
)
stock.mean, stock.q05, stock.q95

Design decisions to nail down

# Decision Default Notes
1 Compositional constraint Sample marginals, renormalize on simplex Independent normal sampling produces invalid compositions (sand+clay>100). Logit-normal / Dirichlet as future option.
2 Distribution shape from Q05/Q50/Q95 Truncated normal; log-normal for SOC-like skewed properties Doesn't uniquely determine a distribution — make the choice an explicit parameter, document the assumption.
3 Spatial autocorrelation in aggregation Independent pixels in 0.3 with prominent caveat Naive MC underestimates regional uncertainty since neighbors are correlated. Proper modeling needs SoilGrids' covariance or geostatistical fitting → 0.4.
4 region vs mask semantics Two separate inputs region defines where you want an answer; mask defines which pixels are valid. Cleaner than overloading one parameter.
5 WorldCover packaging Behind pedotri[landcover] extra mask= accepts any raster/vector. WorldCover fetcher is convenience, not core.
6 SoilGrids fetcher caching On-disk cache from day one, keyed by (lon, lat, property, depth) Agents calling the MCP classify_point tool will hammer the endpoint; un-cached calls risk rate limits and slow agent loops. Default cache dir under $XDG_CACHE_HOME/pedotri/soilgrids/, configurable, with a TTL (SoilGrids releases are infrequent).

Phased release plan

  • 0.3.0 — uncertainty.py; extend classify() and classify_geotiff() with uncertainty kwargs (point + raster); distance-method default, Monte Carlo opt-in.
  • 0.3.x — sources/soilgrids.py fetcher; zonal.zonal_aggregate + AggregateDistribution.combine(); optional sources/worldcover.py behind [landcover] extra.
  • 0.4 — see below.

Each release ships independently. 0.3.0 is enough to be useful on its own (uncertainty-aware classification on existing data); 0.3.x adds the data-acquisition and aggregation primitives that make end-to-end workflows easy.

0.4 — geospatial plumbing

Two themes for 0.4, both deferred from 0.3 to keep that release tight:

Reprojection & multi-resolution alignment

Today _check_rasters_aligned requires inputs to share grid + CRS exactly. Real workflows don't: SoilGrids ships in EPSG:4326 at ~250 m, WorldCover is 10 m in EPSG:4326, a user mask might be a vector in a local UTM zone, DEMs come at 30 m / 90 m. 0.4 introduces:

  • An explicit target grid concept (CRS + transform + shape) that every multi-raster op aligns to.
  • Reprojection helpers with sensible resampling defaults per data type — bilinear for continuous (clay %, SOC), nearest for categorical (class codes, land cover), average for downsampling continuous to coarser grids.
  • Resolution policy for multi-resolution inputs: upsample fine → coarse vs downsample coarse → fine, with the default driven by the use case (zonal aggregation typically downsamples fine inputs to the AOI's working resolution; pixel-level classification typically picks the finest common grid).
  • Caching of reprojection results — these are expensive, especially against SoilGrids at native resolution.

Spatially-correlated uncertainty sampling

The independent-pixel assumption in 0.3 aggregation is a known understatement. 0.4 should offer at least one correlated-sampling mode — likely a simple isotropic Gaussian field with a range parameter, with the option to use SoilGrids' published per-property correlation if/when accessible. Documented as a refinement, not a default switch, since it materially changes the output and the modeling assumption.

Also in 0.4 (lighter)

  • DEM-derived topographic indices: slope, TWI, curvature. Pedology-adjacent; useful primitives for downstream productivity-zone work in a companion package.
  • Companion-package interop guidance: documented protocol for the (array, profile) convention, example of building a productivity-zone pipeline on top of pedotri + sentinelhub-py externally.

Resolved interface decisions

Performance — document + pre-flight check, no dask

No automatic chunking via dask in 0.3. Instead:

  • Document MC cost characteristics in the wiki (rough rule: cost ≈ n_pixels × n_samples × cost(classify_one)).
  • Implement a pre-flight estimate before any MC raster run: compute the expected work and emit a warning above a soft threshold, abort above a hard threshold unless confirm=True is passed. This protects against accidental continental-scale MC runs at N=1000.
  • For huge AOIs, users chunk manually with rasterio windowed reads. Existing classify_geotiff already supports this pattern; extend the docstring with an MC-aware example.

Probability output schema — top-5, multi-band GeoTIFF

Per-class probability output is a 10-band GeoTIFF (top-5 by probability per pixel):

  • Bands 1–5: class codes for the top-5 classes at that pixel, ranked.
  • Bands 6–10: corresponding probabilities (as uint8 0–255 or float32 0–1, configurable).
  • Each band gets a description tag ("rank_1_class", "rank_1_prob", …).
  • A sidecar JSON (or GeoTIFF tag) records the code→name mapping for the active classification + locale.

Why 10 bands and not separate files: keeps the product as a single artifact, plays well with rasterio/GDAL workflows, easy to read top-1 quickly (just bands 1 + 6). KA5 (31 classes) would be unwieldy as 31 separate per-class files; top-5 covers the realistic decision space.

Adjacent idea, separate path: the "probability as alpha" suggestion is great for visualization — extend render_classified_png to optionally encode modal-class probability as the alpha channel, so confidence shows up as transparency on the rendered map. That's a rendering feature, not a data-product schema, and can ship independently.

CLI — no changes in 0.3

Python API only. Existing CLI stays as-is. Revisit if zonal aggregation becomes a common one-off script call.

MCP server — one new tool, smart input modes

Add classify_point to mcp_server.py. Single tool, multiple input modes:

@mcp.tool()
def classify_point(
    sand: float | None = None,
    clay: float | None = None,
    sand_uncertainty: tuple[float, float] | None = None,
    clay_uncertainty: tuple[float, float] | None = None,
    lon: float | None = None,
    lat: float | None = None,
    classification: str = "USDA",
    locale: str = "en",
) -> dict: ...

Behavior:

Inputs What happens
sand, clay only Deterministic classify. Return modal class.
sand, clay, *_uncertainty Uncertainty-aware classify. Return modal + probabilities + entropy.
lon, lat (no sand/clay) Fetch mean + Q05/Q95 from SoilGrids → uncertainty-aware classify automatically.
Both location and explicit values Error — ambiguous; prefer one.

Uncertainties stay optional everywhere — they're only used when provided explicitly or auto-populated from SoilGrids. Plain deterministic flow remains the default for users who hand in lab values.

No separate classify_with_uncertainty tool. No zonal_aggregate MCP tool in 0.3 — it's a heavier operation that fits the Python API better; revisit if there's user demand.

Home

Concepts

Compliance

Use cases

Development

External

Clone this wiki locally