Variable-Rate Prescription Maps

A variable-rate application (VRA) prescription is the artifact that turns an agronomic intention into machine instructions: a set of field polygons, each carrying the rate a sprayer or spreader should apply while its section controller is inside that polygon. Producing one is deceptively easy to get subtly wrong. The zones can be agronomically perfect and still export a rate that violates a product label, or reference acreage that a waterway setback forbids treating, or disagree with the field outline the task controller actually holds. This subsystem, part of Season Planning & Crop Rotation Optimization, owns the path from management zones to a defensible, regulation-clamped prescription — and its non-negotiable rule is that no zone rate ever leaves the pipeline above the regulated ceiling for the product being applied.

Zone delineation and rate agronomy pull in one direction — put more nitrogen where the yield potential justifies it, less where it does not. Regulation pulls in another. A product label carries a maximum application rate derived from EPA registration and 40 CFR tolerances, water features impose setback exclusions, and the field the prescription is written against must be the same polygon the machine believes it is working. A prescription generator that optimizes rate without enforcing all three produces a file that looks correct in a desktop GIS and fails on the ground: the rate controller either applies an illegal rate or the operator overrides the map and the whole exercise is discarded.

The failure is quiet and expensive. An agronomic model asked for 260 kg/ha of nitrogen on a high-yield zone will happily emit it; if the label ceiling is 246 kg/ha, that 14 kg/ha excess is a reportable over-application on every hectare of that zone, discovered only when a season’s records are reconciled against the label. The remedy is structural, not advisory: every per-zone rate is passed through a single clamp against the regulated maximum before it can be serialized, the sprayable envelope is computed by subtracting buffer exclusions from the canonical field boundary, and any zone left with no treatable area is dropped rather than exported as a phantom instruction. This topic covers that pipeline end to end; two focused guides go deeper on the hardest legs — generating the nitrogen rate itself and exporting the result to ISOXML for task controllers.

Prerequisites and dependencies

The prescription generator sits downstream of boundary synchronization and the compliance registry, and upstream of the export format. Confirm four contracts before wiring it in.

  • A canonical field boundary. Zones are reconciled against the authoritative field polygon, not against whatever outline shipped with the soil samples. Drift reconciliation and boundary versioning are owned by Field Boundary Synchronization; this pipeline consumes the boundary it publishes and never redraws it.
  • Buffer-zone exclusions. The un-sprayable setback geometry around water features, wells, and sensitive areas is computed by the Buffer Zone Calculations engine. The prescription subtracts that geometry from the field before any rate is attributed to acreage.
  • Regulated maximum rates. The per-product ceiling used by the clamp originates in EPA/USDA Rule Mapping and is resolved at runtime from a version-stamped registry — never hard-coded per client, so a label revision re-clamps every future export.
  • A soil or yield sample layer. Management zones are derived from a per-point productivity index (soil test, historical yield, or a fused stability index) carried in a WGS84 point layer.

The stack is pinned for reproducibility across planning workers: geopandas>=0.14 and shapely>=2.0 for the spatial algebra, numpy>=1.26 for the zoning quantiles, and Python 3.10+ for match/case and X | None unions. Regulated rates are handled as Decimal end to end so that rounding for export never nudges a value across the ceiling. The governing external references are the EPA pesticide label and use-rate requirements, EPA 40 CFR pesticide tolerances, and USDA NRCS nutrient management (Code 590) guidance for agronomic rate ceilings.

Architecture of this subsystem

The generator is a five-stage transform. A sample layer is binned into management zones; each zone is assigned an agronomic rate; every rate is clamped to the regulated maximum; the zones are reconciled against the canonical boundary minus buffer exclusions; and the surviving zones are serialized as a prescription. The clamp and the reconciliation both draw on the compliance registry, and every stage writes an entry to an append-only audit ledger recording the requested rate, the clamped rate, and the boundary version it was reconciled against.

Variable-rate prescription generation pipeline Left-to-right pipeline of five stages. Field data sources feed Derive Zones, which quantile-bins a productivity index into management zones. Assign Rates maps each zone to a Decimal rate in kilograms per hectare. Clamp caps each rate at the EPA maximum and flags clamped zones. A compliance registry of EPA maximum rates and buffer exclusions feeds the clamp and the reconcile stage. Reconcile clips zones to the canonical field boundary minus buffer exclusions and drops empty zones. The final stage emits a clamped, clipped VRA prescription for ISOXML export. Each stage writes a dashed branch to an append-only audit ledger recording the requested rate, the clamped rate, and the boundary version. Field data sources soil-test grid yield history Compliance registry EPA max application rates buffer-zone exclusions 1 · Derive Zones soil + yield index numpy quantile bins → GeoDataFrame 2 · Assign Rates agronomic target per-zone Decimal kg/ha 3 · Clamp min(rate, EPA max) flag clamped zones never exceed cap 4 · Reconcile clip to boundary subtract buffers drop empty zones VRA prescription clamped + clipped → ISOXML export Append-only audit ledger requested rate · clamped rate · boundary version · excluded acreage

The inputs are a WGS84 sample layer, the canonical field boundary, the buffer exclusion geometry, and an agronomic target; the outputs are a list of clamped, reconciled zone prescriptions ready for serialization. The integration is one-directional: this service consumes boundaries and registry values and produces a prescription — it never edits the boundary registry or the compliance thresholds it reads.

Step-by-step implementation

The reference generator below is a single module with four composable functions plus an export step. It targets Python 3.10+ and is fully typed.

Step 1 — Derive management zones

Management zones are equal-count bins of a per-point productivity index. Quantile edges track the yield or soil gradient without assuming a particular distribution, and the CRS is asserted up front because every later area and reconciliation calculation depends on it.

Parameter Type Default Purpose
samples GeoDataFrame — WGS84 point layer carrying the productivity index.
n_zones int 4 Number of management zones; equal-count quantile bins.
index_col str — Column holding the soil/yield productivity index.

Step 2 — Assign per-zone rates and clamp to the regulated maximum

Each zone’s mean index is mapped onto an agronomic response between a target floor and ceiling, in Decimal. The proposed rate is then passed through clamp_to_regulated_max, which replaces any value above the label ceiling with the ceiling and records that the zone was clamped. This is the guarantee the whole subsystem exists to provide: nothing downstream can raise a rate back above the cap.

Parameter Type Default Purpose
target.floor Decimal — Minimum agronomic rate (kg/ha) for the lowest zone.
target.ceiling Decimal — Maximum agronomic rate before the regulatory clamp.
product str — Registry key selecting the regulated ceiling to clamp against.

Step 3 — Reconcile against boundaries and export

The sprayable envelope is the canonical field boundary with buffer exclusions subtracted. Each zone is intersected with that envelope; a zone that keeps no area is dropped, not exported as a zero-area instruction. The survivors are written to a GIS layer that the ISOXML exporter consumes.

python 
from __future__ import annotations

import logging
from dataclasses import dataclass
from decimal import Decimal, ROUND_HALF_UP

import geopandas as gpd
import numpy as np
from shapely.geometry import base

logger = logging.getLogger("vra.prescription")

# Regulated per-product ceilings are resolved from the compliance registry at
# runtime; the literals here stand in for that lookup. Units are kg/ha of a.i.
EPA_MAX_KG_HA: dict[str, Decimal] = {
    "nitrogen": Decimal("246.0"),
    "phosphorus": Decimal("90.0"),
}
QUANTIZE = Decimal("0.1")  # export resolution the task controller accepts


@dataclass(frozen=True)
class AgronomicTarget:
    product: str
    floor: Decimal    # lowest agronomic rate (kg/ha)
    ceiling: Decimal  # highest agronomic rate BEFORE the regulatory clamp


@dataclass(frozen=True)
class ZonePrescription:
    zone_id: int
    rate_kg_ha: Decimal
    clamped: bool
    area_ha: float
    geometry: base.BaseGeometry  # WGS84 (EPSG:4326)


def derive_zones(samples: gpd.GeoDataFrame, n_zones: int, index_col: str) -> gpd.GeoDataFrame:
    """Bin a per-point productivity index into ordered zones via numpy quantiles."""
    if samples.crs is None or samples.crs.to_epsg() != 4326:
        raise ValueError("sample layer must be WGS84 (EPSG:4326)")
    values = samples[index_col].to_numpy(dtype=float)
    # Quantile edges give equal-count zones that follow the yield/soil gradient.
    edges = np.quantile(values, np.linspace(0.0, 1.0, n_zones + 1))
    edges[0], edges[-1] = -np.inf, np.inf  # keep the open ends inclusive
    zone_id = np.digitize(values, edges[1:-1], right=False)
    zoned = samples.assign(zone_id=zone_id)
    # Dissolve the points of each zone into one polygon (Voronoi/buffer upstream).
    return zoned.dissolve(by="zone_id", aggfunc={index_col: "mean"}).reset_index()


def assign_rate(zone_frac: float, target: AgronomicTarget) -> Decimal:
    """Map a zone's normalized index (0..1) to an agronomic rate in kg/ha."""
    span = target.ceiling - target.floor
    frac = Decimal(str(max(0.0, min(1.0, zone_frac))))  # guard the interpolation
    return (target.floor + span * frac).quantize(QUANTIZE, rounding=ROUND_HALF_UP)


def clamp_to_regulated_max(rate: Decimal, product: str) -> tuple[Decimal, bool]:
    """CLAMP a proposed rate to the regulated ceiling. Returns (rate, was_clamped)."""
    ceiling = EPA_MAX_KG_HA[product]
    if rate > ceiling:
        logger.warning(
            '{"event":"rate_clamped","product":"%s","requested":"%s","ceiling":"%s"}',
            product, rate, ceiling,
        )
        return ceiling, True
    return rate, False


def _area_hectares(geom: base.BaseGeometry) -> float:
    """Area in hectares via an equal-area reprojection (EPSG:6933)."""
    projected = gpd.GeoSeries([geom], crs=4326).to_crs(6933)
    return float(projected.area.iloc[0]) / 10_000.0


def build_prescription(
    zones: gpd.GeoDataFrame,
    field_boundary: base.BaseGeometry,
    buffer_exclusions: base.BaseGeometry,
    target: AgronomicTarget,
    index_col: str,
) -> list[ZonePrescription]:
    """Zone -> rate -> clamp -> reconcile against the canonical boundary and buffers."""
    plan_area = field_boundary.difference(buffer_exclusions)  # sprayable acreage only
    idx = zones[index_col].to_numpy(dtype=float)
    lo, hi = float(idx.min()), float(idx.max())
    out: list[ZonePrescription] = []
    for row in zones.itertuples(index=False):
        # Uniform index (hi == lo) degrades to a flat plan instead of div-by-zero.
        frac = 0.0 if hi == lo else (getattr(row, index_col) - lo) / (hi - lo)
        rate = assign_rate(frac, target)
        rate, clamped = clamp_to_regulated_max(rate, target.product)
        geom = row.geometry.intersection(plan_area)  # keep only sprayable part
        if geom.is_empty:
            logger.info(
                '{"event":"zone_dropped","zone_id":%d,"reason":"outside_plan_area"}',
                int(row.zone_id),
            )
            continue
        out.append(
            ZonePrescription(int(row.zone_id), rate, clamped, _area_hectares(geom), geom)
        )
    logger.info(
        '{"event":"prescription_built","zones":%d,"clamped":%d}',
        len(out), sum(p.clamped for p in out),
    )
    return out


def export_prescription(rx: list[ZonePrescription], path: str) -> None:
    """Serialize the clamped, reconciled zones to a GIS layer for the controller."""
    frame = gpd.GeoDataFrame(
        {
            "zone_id": [p.zone_id for p in rx],
            "rate_kg_ha": [float(p.rate_kg_ha) for p in rx],
            "clamped": [p.clamped for p in rx],
            "area_ha": [round(p.area_ha, 4) for p in rx],
        },
        geometry=[p.geometry for p in rx],
        crs="EPSG:4326",
    )
    frame.to_file(path, driver="GeoJSON")
    logger.info('{"event":"prescription_exported","path":"%s","zones":%d}', path, len(rx))

Expected log output when one high-yield zone is clamped to the label ceiling:

text 
{"event": "rate_clamped", "product": "nitrogen", "requested": "261.4", "ceiling": "246.0"}
{"event": "prescription_built", "zones": 4, "clamped": 1}
{"event": "prescription_exported", "path": "field_12_rx.geojson", "zones": 4}

Expected log output when a zone falls entirely inside a buffer exclusion:

text 
{"event": "zone_dropped", "zone_id": 0, "reason": "outside_plan_area"}
{"event": "prescription_built", "zones": 3, "clamped": 0}

Edge cases and known failure modes

Condition Symptom Fix
Agronomic model requests a rate above the label ceiling Export would carry an illegal over-application on every hectare of the zone clamp_to_regulated_max caps the value at the ceiling, sets clamped=True, and logs rate_clamped.
A zone lies entirely within a buffer exclusion Zone would export a rate for un-sprayable acreage The intersection empties; the zone is dropped with zone_dropped and never serialized.
Sample layer is not WGS84 Areas and reconciliation are silently wrong derive_zones asserts EPSG:4326 and raises before any zoning happens.
Zones and boundary come from different vintages Sliver gaps/overlaps at field edges Reconcile against the current boundary from Field Boundary Synchronization; snap zones to the active version.
All samples share one index value Quantile edges collapse to a single zone The hi == lo guard yields a uniform rate — a flat-rate plan rather than a divide-by-zero.
Buffer exclusion geometry is self-intersecting difference raises TopologyException Repair with make_valid upstream; the buffer engine must publish valid geometry.
Rounding to export resolution nudges a rate above the ceiling Controller rejects the zone Quantize first, then clamp — the clamp is always the final gate before export.

Compliance and audit integration

Every prescription is a compliance artifact, and the clamp is where the compliance boundary is enforced. The regulated maximum for each product is resolved from the versioned threshold set in EPA/USDA Rule Mapping, which draws the ceiling from the product’s EPA pesticide label and the applicable 40 CFR tolerance. Because the ceiling is looked up rather than hard-coded, a label revision that lowers a maximum re-clamps every prescription generated afterward, and a prescription generated last season replays against the ceiling that was in force at the time.

Buffer exclusions carry the second half of the compliance guarantee. The setback geometry the reconciliation subtracts is produced by the Buffer Zone Calculations engine from EPA Worker Protection Standard and label setback rules; treating acreage inside a setback is a violation regardless of the rate applied. Each zone’s audit entry records the requested rate, the clamped rate, whether the clamp fired, the boundary version reconciled against, and the excluded acreage removed — a continuous trail from an exported prescription back to the agronomic target and the regulatory ceiling that constrained it. Following USDA NRCS Code 590 nutrient-management expectations, the clamped-zone count and total excluded area should be surfaced as planning-health metrics.

Verification

Confirm correct operation in staging with a reproducible fault-injection scenario before promoting a change:

  1. Golden prescription. Feed a four-zone sample layer whose highest agronomic rate sits just under the ceiling and assert build_prescription returns four zones, clamped=False on all of them, and areas that sum to the sprayable envelope.
  2. Inject an over-ceiling zone. Set the target ceiling above the regulated maximum so the top zone would exceed it, and assert that zone’s exported rate equals the registry ceiling exactly, clamped=True, and a rate_clamped entry is written — the over-application never reaches the file.
  3. Inject a fully-excluded zone. Place a zone entirely inside the buffer geometry and assert it is dropped with zone_dropped, the exported zone count drops by one, and no zero-area geometry is serialized.
  4. Inject a wrong CRS. Pass a sample layer in a projected CRS and assert derive_zones raises before any rate is computed, so a mis-projected plan can never be exported.

A run passes only when each injected fault produces the expected safe result and the audit ledger holds one traceable entry per zone. The rate agronomy and the export format each have their own verification, covered in the guides below.

Frequently Asked Questions

Why clamp every zone instead of validating the prescription after it is built? Because a post-hoc validator can reject a bad file but cannot make a good one — it leaves the pipeline able to emit an illegal rate right up to the final check. Clamping inside the rate assignment means an over-ceiling value is impossible to serialize: the clamp replaces it with the ceiling, flags the zone, and logs the event, so the export is correct by construction rather than by a gate that a refactor could bypass.

What happens to a management zone that falls entirely inside a buffer exclusion? It is dropped, not exported with a zero rate. The zone is intersected with the sprayable envelope (boundary minus buffers); if nothing survives, a zone_dropped entry is written and the zone never appears in the prescription. Exporting a zero-area or zero-rate instruction would clutter the task controller and obscure the fact that the acreage is legally un-treatable.

Why reconcile against the canonical boundary instead of the outline that came with the samples? Sample layers are frequently drawn against a stale or approximate field outline. Attributing rates to that outline produces slivers and phantom acreage where it disagrees with the boundary the machine actually holds. Reconciling against the versioned canonical boundary keeps the prescription geometry identical to the field the controller works, which is what makes the applied-area records defensible.

Why use Decimal for the rates rather than floats? Regulated rates are compared against a hard ceiling and rounded to an export resolution. Binary floats make both operations imprecise — a value that should equal the ceiling can land a hair above it and trip a false violation, or a rounding step can drift. Decimal with an explicit quantization keeps the clamp exact and the exported value reproducible.

Up: Season Planning & Crop Rotation Optimization

Subtopics