How to write a variable-rate prescription map

Problem statement

A prescription map is a polygon layer with a rate column, and the machine that reads it has opinions about all three words. The polygons must be simple, non-overlapping, and no finer than the implement can follow; the rate must be in the units the controller expects; and the file must be in a format and a coordinate system the terminal accepts.

Getting the agronomy right and the file wrong is the common outcome. The map loads, the controller applies a default rate everywhere, and nobody notices until the invoice for the fertiliser arrives.

This guide builds a prescription from zones, and gets the mechanical details right.

Quick answer

import geopandas as gpd, shapely

RATES = {"low": 120.0, "mid": 160.0, "high": 200.0}   # kg N/ha by zone — placeholders, see step 1
METRES = "EPSG:3003"                                   # a grid in metres that suits the field

def absorb_fragments(rx, min_area):
    """Hand each part under min_area to the neighbour it shares the longest edge with."""
    parts = rx.explode(index_parts=False).reset_index(drop=True)
    while (parts.area < min_area).any():
        i = parts.area.idxmin()                        # smallest first
        g = parts.geometry[i]
        near = parts.drop(i)
        near = near[near.intersects(g)]
        edge = near.boundary.intersection(g.boundary).length
        parts = parts.drop(i)
        if len(edge) and edge.max() > 0:
            j = edge.idxmax()
            parts.loc[j, "geometry"] = parts.geometry[j].union(g)
    return parts.dissolve("zone").explode(index_parts=False).reset_index()

rx = zones.to_crs(METRES)
rx["geometry"] = rx.intersection(field)                # nothing outside the field
# half a 24 m working width; coverage_simplify moves each shared edge once, for both sides
rx["geometry"] = shapely.coverage_simplify(rx.geometry.values, 12)
rx = absorb_fragments(rx, min_area=1000)               # 0.1 ha
rx["rate"] = rx["zone"].map(RATES)
rx["rate_units"] = "kg/ha"
rx["product"] = "urea 46N"

overlap = rx.area.sum() - rx.union_all().area
gap = field.area - rx.union_all().area
assert rx["rate"].notna().all(), "some polygons have no rate"
assert overlap < 1, f"{overlap:.0f} m² of overlap"
assert gap < 0.01 * field.area, f"{gap / 1e4:.2f} ha of the field has no rate"
rx.to_crs(4326).to_file("prescription.shp")           # most terminals want WGS84
print(f"{len(zones.explode())} zone parts → {len(rx)} polygons, smallest {rx.area.min() / 1e4:.2f} ha; "
      f"overlap {overlap:.2f} m², uncovered {gap:.0f} m²")

561 zone parts → 17 polygons, smallest 0.10 ha; overlap 0.00 m², uncovered 601 m²

The zones here come straight from a kriged yield raster, so each is a scatter of 2.2 m pixel clumps: 561 parts, half of them under 5 m². The first version of this snippet simplified each polygon on its own with simplify(0.0001).buffer(0), which moves the two copies of every shared edge independently — on this field the overlap assertion stopped it, and Example 1 measures the damage at 12 m: 1.54 ha of overlap and 1.52 ha uncovered, out of 19.5. shapely.coverage_simplify treats the zones as one coverage and moves each shared edge once. Fragments are then handed to a neighbour rather than dropped, because dropping them is how a prescription gets holes. The measurements are in the field's own grid; the snippet once used the Dutch national grid, EPSG:28992, which reads this field as 19.48 ha against 19.51 — only 0.2% out, but a UTM zone one off (33N) makes it 19.97.

Three assertions before writing. Every one of them corresponds to a real failure in the field: a missing rate becomes zero or a default, an overlap makes the controller choose unpredictably, and a gap falls through to whatever default the terminal holds.

Step-by-step solution

1. Decide the rates from something other than the zones

Zones say where the field differs; they do not say what rate to apply. The rate comes from an agronomic model, a nutrient balance, a soil test, or — best — an on-farm trial that measured the response. How to analyse an on-farm strip trial in Python covers the last.

2. Check the rates against the machine's limits

Every applicator has a minimum and a maximum rate and a rate resolution. A prescription outside those is clipped silently. Ask, and clip deliberately.

3. Simplify to the implement's resolution

A spreader with a 24 m boom cannot follow a 5 m boundary, and a spinner's distribution pattern is several metres wide in any case. Simplifying at roughly half the working width removes detail that cannot be applied.

4. Remove fragments and slivers

A polygon smaller than the machine's response distance — speed times the actuator's lag, typically 5–15 m of travel — causes the rate to change and change back before it takes effect. Merge them into their neighbours; deleting them leaves holes.

5. Make the coverage clean

No overlaps, no gaps, valid geometry, single-part polygons. A gap gets the controller's default rate, which is usually not what you intended.

6. Write in the format the terminal reads

Shapefile is still the most widely accepted interchange. ISOXML (ISO 11783-10) is the standard for task data and is what most modern terminals prefer. Several manufacturers have their own formats. Ask before exporting, and test with a small file.

7. Report the totals before the file leaves

Area by rate, total product, and the average rate. A prescription that applies 40% more nitrogen than the uniform plan is a decision, and it should be a visible one.

Code examples

The data used on this page

A real field: the 22 ha corn field at Porto Felloni near Ferrara whose three yield-monitor harvests are published on Zenodo (CC BY 4.0), used on yield monitor data explained and laying out a soil sampling grid. The field outline is the concave hull of the 2017 yield points — 19.5 ha of it, because the authors had trimmed the borders — and the zones are the mean of the three kriged yield maps in the archive, split into thirds. The files' declared CRS is wrong; they are in Gauss-Boaga (EPSG:3003).

import os
import zipfile

import geopandas as gpd
import numpy as np
import rasterio
import requests
import shapely
from rasterio.features import shapes
from shapely.geometry import shape

if not os.path.exists("FIELD_YIELD"):
    reply = requests.get("https://zenodo.org/records/14216165/files/FIELD_YIELD.zip?download=1", timeout=300)
    reply.raise_for_status()
    open("FIELD_YIELD.zip", "wb").write(reply.content)
    zipfile.ZipFile("FIELD_YIELD.zip").extractall(".")

CRS = "EPSG:3003"
points = gpd.read_file("FIELD_YIELD/yield2017/yield2017.shp").set_crs(CRS, allow_override=True)
field = shapely.concave_hull(points.union_all(), ratio=0.05)

maps = []
for year in (2016, 2017, 2018):
    with rasterio.open(f"FIELD_YIELD/yield{year}krig/{year}yieldkrig.tif") as src:
        maps.append(src.read(1, masked=True))
        transform = src.transform
mean = np.ma.mean(maps, axis=0)
breaks = np.quantile(mean.compressed(), [1 / 3, 2 / 3])
classes = np.digitize(mean.filled(np.nan), breaks).astype("int16")
classes[mean.mask] = -1
zones = gpd.GeoDataFrame(
    [{"zone": ("low", "mid", "high")[int(v)], "geometry": shape(g)}
     for g, v in shapes(classes, mask=classes >= 0, transform=transform)], crs=CRS).dissolve("zone").reset_index()
print(f"field {field.area / 1e4:.1f} ha; zone areas, ha: "
      f"{dict(zip(zones['zone'], (zones.area / 1e4).round(1)))}")

field 19.5 ha; zone areas, ha: {'high': 6.7, 'low': 6.7, 'mid': 6.7}

Example 1 — build the prescription with the checks

import geopandas as gpd, shapely

def build_prescription(zones, field, rates, product, units="kg/ha",
                       working_width_m=24, min_area_ha=0.1,
                       rate_limits=(0, 400), rate_step=None, crs_out=4326,
                       area_crs=3003):
    """Zones in any CRS, field outline in area_crs. Uses absorb_fragments from the quick answer."""
    rx = zones.to_crs(area_crs)
    rx["geometry"] = rx.intersection(field)
    rx["rate"] = rx["zone"].map(rates)
    missing = rx["rate"].isna().sum()
    if missing:
        raise ValueError(f"{missing} zones have no rate")

    lo, hi = rate_limits
    clipped = ((rx["rate"] < lo) | (rx["rate"] > hi)).sum()
    rx["rate"] = rx["rate"].clip(lo, hi)
    if clipped:
        print(f"clipped {clipped} zones to the machine limits {lo}–{hi}")
    if rate_step:
        rx["rate"] = (rx["rate"] / rate_step).round() * rate_step

    # one rate per zone, so dissolving on zone keeps the rates with their zones
    rx["geometry"] = shapely.coverage_simplify(rx.geometry.values, working_width_m / 2)
    before = rx.explode().shape[0]
    rx = absorb_fragments(rx, min_area_ha * 1e4)
    rx["rate"] = rx["zone"].map(rx.groupby("zone")["rate"].first())
    print(f"{before} parts → {len(rx)} polygons after absorbing fragments under {min_area_ha} ha")

    rx["product"] = product
    rx["rate_units"] = units
    rx["area_ha"] = rx.area / 1e4

    covered = rx.union_all()
    overlap = rx.area.sum() - covered.area
    gap = field.area - covered.area
    if overlap > 1:
        raise ValueError(f"{overlap:.0f} m² of overlapping polygons")
    print(f"uncovered {gap:.0f} m² ({gap / field.area:.2%} of the field)")

    total = float((rx["rate"] * rx["area_ha"]).sum())
    print(rx.groupby("rate")["area_ha"].sum().round(2).to_string())
    print(f"total {total:,.0f} {units.split('/')[0]}; "
          f"area-weighted mean rate {total / rx['area_ha'].sum():.1f} {units}")
    return rx.to_crs(crs_out)

Example 2 — compare against the uniform plan

import numpy as np, pandas as pd

def compare_with_uniform(rx, uniform_rate, units="kg/ha"):
    area = rx["area_ha"].sum()
    variable = float((rx["rate"] * rx["area_ha"]).sum())
    flat = uniform_rate * area
    rows = rx.groupby("rate")["area_ha"].sum().reset_index()
    rows["share"] = rows["area_ha"] / area
    rows["vs_uniform"] = rows["rate"] - uniform_rate
    print(rows.round(2).to_string(index=False))
    print(f"\n{area:.1f} ha; variable {variable:,.0f} vs uniform {flat:,.0f} "
          f"({variable / flat - 1:+.1%})")
    return {"area_ha": area, "variable_total": variable, "uniform_total": flat,
            "difference_pct": variable / flat - 1}

rx = build_prescription(zones, field, {"low": 120.0, "mid": 160.0, "high": 200.0}, "urea 46N",
                        rate_limits=(0, 180), rate_step=10)
summary = compare_with_uniform(rx, uniform_rate=160)

# what the first version did: simplify each polygon on its own
naive = zones.intersection(field).simplify(12)
print(f"\nper-polygon simplify: {(naive.area.sum() - naive.union_all().area) / 1e4:.2f} ha overlap, "
      f"{(field.area - naive.union_all().area) / 1e4:.2f} ha uncovered")
# and an overlap test by predicate, after reprojecting to WGS84
pairs = gpd.sjoin(rx, rx, predicate="overlaps")
print(f"{len(pairs) // 2} pairs 'overlap' in WGS84 by predicate; "
      f"their shared area is {sum(rx.geometry[a].intersection(rx.geometry[b]).area for a, b in zip(pairs.index, pairs.index_right)):.1e} square degrees")

field_area = {epsg: gpd.GeoSeries([field], crs=CRS).to_crs(epsg).area.iloc[0] / 1e4
              for epsg in (3003, 32633, 28992)}
print("the field in three grids:", ", ".join(f"EPSG:{k} {v:.3f} ha" for k, v in field_area.items()))

clipped 1 zones to the machine limits 0–180
559 parts → 17 polygons after absorbing fragments under 0.1 ha
uncovered 601 m² (0.31% of the field)
rate
120.0    5.82
160.0    6.95
180.0    6.68
total 3,013 kg; area-weighted mean rate 154.9 kg/ha
 rate  area_ha  share  vs_uniform
120.0     5.82   0.30       -40.0
160.0     6.95   0.36         0.0
180.0     6.68   0.34        20.0

19.5 ha; variable 3,013 vs uniform 3,113 (-3.2%)

per-polygon simplify: 1.54 ha overlap, 1.52 ha uncovered
6 pairs 'overlap' in WGS84 by predicate; their shared area is 1.6e-15 square degrees
the field in three grids: EPSG:3003 19.513 ha, EPSG:32633 19.965 ha, EPSG:28992 19.477 ha

Here the variable plan uses 3.2% less than a flat 160 kg/ha, and not for agronomic reasons: the high zone asked for 200 and the spreader's 180 limit clipped it. The 561 pixel clumps became 17 polygons, 601 m² of the field edge is left without a polygon, and the overlap check had to be made by area — the overlaps predicate after reprojecting to WGS84 flags six pairs whose shared area is 1.6 × 10⁻¹⁵ square degrees, a few square millimetres.

Publishing the comparison is what turns a prescription from a black box into a decision. A variable plan that costs 15% more than the flat one may well be right; it should not be a surprise.

Example 3 — fill the gaps so the controller never uses a default

import geopandas as gpd

def fill_to_field(rx, field, field_crs, default_rate, area_crs=3003):
    """Any part of the field with no polygon gets an explicit rate."""
    rx_m = rx.to_crs(area_crs)
    # the field comes in its own CRS: labelling it with rx.crs (WGS84 after export)
    # would read metres as degrees
    f = gpd.GeoDataFrame(geometry=[field], crs=field_crs).to_crs(area_crs)
    gaps = gpd.overlay(f, rx_m[["geometry"]], how="difference").explode(index_parts=False)
    gaps = gaps[gaps.area > 100]
    if len(gaps):
        print(f"{len(gaps)} gaps totalling {gaps.area.sum()/1e4:.2f} ha "
              f"→ filling at {default_rate}")
        gaps["rate"] = default_rate
        gaps["product"] = rx_m["product"].iloc[0]
        gaps["rate_units"] = rx_m["rate_units"].iloc[0]
        rx_m = gpd.pd.concat([rx_m, gaps], ignore_index=True)
    return gpd.GeoDataFrame(rx_m, crs=area_crs).to_crs(rx.crs)

filled = fill_to_field(rx, field, CRS, default_rate=160)
print(f"{len(rx)} polygons → {len(filled)}; covered {filled.to_crs(3003).area.sum() / 1e4:.2f} "
      f"of {field.area / 1e4:.2f} ha")

1 gaps totalling 0.02 ha → filling at 160
17 polygons → 18; covered 19.47 of 19.51 ha

On this field that is one 0.02 ha gap along the edge; the rest of the 601 m² is slivers under the 100 m² threshold. An explicit fill is always better than relying on the controller's default, because the default is a machine setting that somebody may have changed.

Explanation

Why detail finer than the implement is worse than useless

A boom sprayer with a 24 m width applies an average over that width; a spinner's pattern is wider still. A prescription boundary at 5 m resolution cannot be applied — and every boundary crossing costs the controller a rate change, which has a response lag. The result is a machine oscillating between rates while applying an average that has nothing to do with the map.

Why gaps are more dangerous than overlaps

An overlap is usually resolved by the controller taking the first or the last polygon, which is wrong in a small area. A gap falls through to the terminal's default rate, which may be zero, may be the last rate used, or may be a value somebody set for a different field. It is the one failure that can apply a completely unintended rate across a large area.

Why the shapefile persists

It is universally readable, it is simple enough that a terminal's embedded software can parse it, and every manufacturer has supported it for decades. Its limitations — ten-character field names, no real date type — barely matter for a file with three columns. ISOXML is the proper standard and is what modern terminals prefer, but shapefile remains the safest thing to hand somebody whose terminal you have not seen.

Why rate resolution matters

Applicators change rate in steps, and the step can be coarse — 10 kg/ha is common on older spreaders. A prescription with rates of 137 and 142 kg/ha is applied as the same rate, so the zones it distinguishes are not distinguished on the ground. Rounding to the machine's step before writing makes the map honest about what will happen.

Edge cases or notes

• Ask for the terminal's format first. Exporting blind wastes a field operation.
• Units differ by product. kg/ha, l/ha, seeds/m², thousand seeds/ha.
• Product form matters. 120 kg N/ha is not 120 kg/ha of urea.
• Shapefile truncates field names to ten characters.
• Most terminals want WGS84. Some want the local grid; ask.
• Headlands are often a separate rate, or excluded entirely.
• Test with one small file before a whole farm.
• Keep the as-applied data. It is the only evidence of what happened.

Internal links

• A prescription map is rejected by the machine controller — when the file will not load
• Management zones explained — where the polygons come from
• How to build management zones with clustering in Python — producing them
• How to analyse an on-farm strip trial in Python — where the rates should come from
• How to simplify geometry in GeoPandas — simplifying to the implement
• How to fix gaps and overlaps in a polygon coverage — making the coverage clean
• Shapefile column names are truncated — the ten-character limit
• How to package GIS deliverables in Python — what to hand over

FAQ

What format should a prescription map be in?

Ask the terminal. Shapefile is the most widely accepted interchange; ISOXML is the proper standard and what most modern terminals prefer.

How much should I simplify the polygons?

To about half the implement's working width. Finer detail cannot be applied and only makes the controller cross more boundaries.

What happens if my polygons have gaps?

The controller applies its default rate there, which may be zero or a value left from another field. Fill gaps explicitly.

Why does the machine apply the same rate to two different zones?

Because the applicator's rate resolution is coarser than the difference. Round the rates to the machine's step before writing.

Should I include the headland?

Usually as its own rate or excluded entirely — it is compacted, differently fertilised and turned on.

What should I check before exporting?

Every polygon has a rate, no overlaps, no uncovered area, no fragments below the machine's response distance, valid single-part geometry, and the totals compared with the uniform plan.