A prescription map is rejected by the machine controller

On this page

Problem statement

The file is on the stick, the terminal sees it, and one of four things happens: it refuses to import, it imports and shows an empty field, it imports and applies a constant rate, or it applies rates that bear no relation to the map.

All four are mechanical rather than agronomic, and they are diagnosed by inspecting the file rather than by reloading it. The usual causes are a coordinate system the terminal does not accept, a rate column whose name or type the terminal does not recognise, geometry that is invalid, overlapping or far too fine, and a format that the terminal supports in a different version from the one you wrote.

The examples use a real field: the 20 ha maize field near Ferrara whose kriged yield maps for 2016–2018 are published on Zenodo. Three management zones are clustered from the average of the three years and given rates of 180, 150 and 120 kg N/ha — illustrative rates, not a recommendation — and the prescription is then written the ways that go wrong in practice.

Quick answer

Check the file the way the terminal will read it:

import geopandas as gpd

def inspect(path, rate_col="rate"):
    rx = gpd.read_file(path)
    print(path)
    print("  CRS:", rx.crs.to_string() if rx.crs else None, "| bounds:", rx.total_bounds.round(3).tolist())
    print("  columns:", rx.columns.tolist())
    print("  rate dtype:", rx[rate_col].dtype if rate_col in rx else "no rate column")
    print("  geometry:", rx.geom_type.value_counts().to_dict(),
          "| invalid:", int((~rx.geometry.is_valid).sum()), "| empty:", int(rx.geometry.is_empty.sum()))
    print("  overlapping pairs:", len(gpd.sjoin(rx, rx, predicate="overlaps")) // 2)

for path in ["out/rx_good.shp", "out/rx_no_prj.shp", "out/rx_longname.shp", "out/rx_text.shp", "out/rx_simplified.shp"]:
    inspect(path)

out/rx_good.shp
  CRS: EPSG:4326 | bounds: [12.209, 44.691, 12.216, 44.698]
  columns: ['zone', 'rate', 'geometry']
  rate dtype: int64
  geometry: {'Polygon': 335} | invalid: 0 | empty: 0
  overlapping pairs: 14
out/rx_no_prj.shp
  CRS: None | bounds: [754285.083, 4953633.137, 754836.556, 4954391.963]
  columns: ['zone', 'rate', 'geometry']
  rate dtype: int64
  geometry: {'Polygon': 335} | invalid: 0 | empty: 0
  overlapping pairs: 0
out/rx_longname.shp
  CRS: EPSG:4326 | bounds: [12.209, 44.691, 12.216, 44.698]
  columns: ['zone', 'prescripti', 'geometry']
  rate dtype: no rate column
  geometry: {'Polygon': 335} | invalid: 0 | empty: 0
  overlapping pairs: 14
out/rx_text.shp
  CRS: EPSG:4326 | bounds: [12.209, 44.691, 12.216, 44.698]
  columns: ['zone', 'rate', 'geometry']
  rate dtype: str
  geometry: {'Polygon': 335} | invalid: 0 | empty: 0
  overlapping pairs: 14
out/rx_simplified.shp
  CRS: EPSG:4326 | bounds: [12.209, 44.691, 12.216, 44.698]
  columns: ['zone', 'rate', 'geometry']
  rate dtype: int64
  geometry: {'Polygon': 335} | invalid: 0 | empty: 0
  overlapping pairs: 80

Each broken file shows its fault in one line: no CRS and bounds in the hundreds of thousands, a rate column renamed prescripti, a rate column that is text, and 80 overlapping pairs after simplification. None of those files raised an error when it was written. And the "good" file has 14 overlapping pairs of its own, which step 4 explains.

Six failures; every one is visible in the file before it reaches the terminal.

Match symptom to cause, check CRS and the rate column

1. Match the symptom to the cause

symptom	usual cause
will not import at all	format or version the terminal does not read
imports, field is empty	wrong CRS, or geometry far from the machine's position
imports, constant rate applied	rate column missing, misnamed, or stored as text
rates look scrambled	overlapping polygons, or multipart geometry
rate applied only in part of the field	gaps falling through to the default
machine oscillates	polygons finer than the controller's response

2. Check the coordinate system

Terminals generally expect WGS 84 longitude and latitude (ISOXML, the ISO 11783-10 task format, uses it); some accept a projected grid. A shapefile with no .prj carries no CRS at all:

no_prj = gpd.read_file("out/rx_no_prj.shp")
print(no_prj.crs, no_prj.total_bounds.round())

None [ 754285. 4953633.  754837. 4954392.]

Eastings of 754,000 and northings of 4,953,000 are UTM zone 32N metres. Read as degrees they are nowhere on Earth, so a terminal that assumes WGS 84 has nothing to draw.

3. Check the rate column's name and type

Shapefile truncates field names to ten characters, so prescription_rate becomes prescripti — GeoPandas warns, the terminal just does not find the name it looks for.

The type problem is subtler than it is usually described. A missing number on its own does not make a column text: pandas stores it as NaN in a float column, which the shapefile writes as a numeric field. Text comes from the values themselves — rates copied from a spreadsheet as 180 kg/ha, or a placeholder like n/a:

import pandas as pd
import pyogrio

print(pd.Series([180, None, 150]).dtype, "->", pyogrio.read_info("out/rx_nan.shp")["dtypes"].tolist())
print(pd.Series(["180 kg/ha", "150 kg/ha"]).dtype, "->", pyogrio.read_info("out/rx_text.shp")["dtypes"].tolist())

text = gpd.read_file("out/rx_text.shp")
text["rate"] = pd.to_numeric(text["rate"].str.extract(r"([\d.]+)")[0], errors="coerce")
assert text["rate"].notna().all(), "non-numeric rates"
print(text["rate"].dtype, sorted(text["rate"].unique().tolist()))

float64 -> ['int64', 'float64']
str -> ['int64', 'object']
int64 [120, 150, 180]

4. Check the geometry

Valid, single-part, non-overlapping polygons with no holes where you did not intend them. Zones polygonised from a raster are valid and never overlap, but they are ragged and full of holes:

print(len(zones), "polygons;", int(zones.geometry.interiors.map(len).gt(0).sum()), "with holes;",
      int((zones.geom_type == "MultiPolygon").sum()), "multipart")

335 polygons; 19 with holes; 0 multipart

Polygonised zones do not overlap in the grid they were drawn in, yet the WGS 84 file above has 14 overlapping pairs. The cause is T-junctions: where three zones meet, one polygon's edge carries a vertex that its neighbour's straight edge runs through without. In UTM the two edges are the same line. Reprojection bends the straight edge into a slightly different line than the one through the extra vertex, and the neighbours now overlap by a sliver. Snapping every polygon to the layer's own vertices inserts the missing ones:

import shapely

def ov(frame):
    return len(gpd.sjoin(frame, frame, predicate="overlaps")) // 2

vertices = shapely.multipoints(shapely.get_coordinates(zones.geometry.values))
snapped = zones.set_geometry(shapely.snap(zones.geometry.values, vertices, 0.001), crs=zones.crs)
print("coverage valid in UTM:", shapely.coverage_is_valid(zones.geometry.values), "->",
      shapely.coverage_is_valid(snapped.geometry.values))
print("overlapping pairs after reprojecting:", ov(zones.to_crs(4326)), "->", ov(snapped.to_crs(4326)))
print("vertices added:", shapely.get_num_coordinates(snapped.geometry.values).sum()
      - shapely.get_num_coordinates(zones.geometry.values).sum())

coverage valid in UTM: False -> True
overlapping pairs after reprojecting: 14 -> 0
vertices added: 28

5. Check for gaps

Any part of the field with no polygon gets the controller's default rate — often zero, sometimes the last value used. Gaps are usually created by cleaning: dropping small polygons, or simplifying each polygon on its own, opens holes between neighbours.

field = zones.union_all()
kept = zones[zones.area >= 1000]
print(f"dropping polygons under 0.1 ha leaves {field.difference(kept.union_all()).area / 1e4:.2f} ha uncovered")

dropping polygons under 0.1 ha leaves 2.38 ha uncovered

6. Check the resolution against the machine

A polygon narrower than the distance the machine travels during the actuator's response — speed times lag, commonly quoted as 5–15 m — changes the rate and changes it back before it takes effect. A zone map polygonised straight from a 2.2 m raster is far finer than that:

print(f"{(zones.area < 1000).sum()} of {len(zones)} polygons are under 0.1 ha; "
      f"{(zones.area < 5).sum()} are a single 2.2 m cell")

326 of 335 polygons are under 0.1 ha; 135 are a single 2.2 m cell

7. Test with a small file first

One field, three zones, on the actual terminal, before a whole farm. The failure modes are terminal-specific and no amount of specification reading substitutes for a test.

Eight checks, all runnable in Python, none of which the terminal will report.

A pre-flight check and common repairs

The data used on this page

The kriged yield maps of the Ferrara field for 2016, 2017 and 2018 from the Zenodo archive 14216165, on a 2.2 m grid treated as UTM zone 32N (EPSG:32632; the archive's Gauss-Boaga eastings with the leading million dropped land within about 37 m of it). Each year is standardised, the three are averaged, and k-means splits the field into three zones, polygonised with rasterio.features.shapes into zones. Then the prescription is written six ways into out/: correctly (rx_good), without its .prj, with a long rate column name, with a missing rate, with rates as text, and with each polygon simplified to 12 m on its own.

import os
import warnings
import zipfile

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
import shapely
from rasterio import features
from sklearn.cluster import KMeans

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(".")

years = {}
for year in (2016, 2017, 2018):
    with rasterio.open(f"FIELD_YIELD/yield{year}krig/{year}yieldkrig.tif") as src:
        years[year] = src.read(1, masked=True).filled(np.nan)
        transform = src.transform
inside = np.isfinite(years[2016])
mean_z = np.mean([(a - np.nanmean(a)) / np.nanstd(a) for a in years.values()], axis=0)

kmeans = KMeans(3, n_init=10, random_state=0).fit(mean_z[inside].reshape(-1, 1))
rank = np.argsort(np.argsort(kmeans.cluster_centers_.ravel()))       # 1 = lowest-yielding zone
zone_grid = np.zeros(inside.shape, "int16")
zone_grid[inside] = rank[kmeans.labels_] + 1

shapes = list(features.shapes(zone_grid, mask=inside, transform=transform))
zones = gpd.GeoDataFrame({"zone": [int(v) for _, v in shapes]},
                         geometry=[shapely.geometry.shape(g) for g, _ in shapes], crs=32632)
RATES = {1: 180, 2: 150, 3: 120}                                       # kg N/ha, illustrative
zones["rate"] = zones["zone"].map(RATES)

os.makedirs("out", exist_ok=True)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    zones.to_crs(4326).to_file("out/rx_good.shp")
    zones.to_file("out/rx_no_prj.shp")
    os.remove("out/rx_no_prj.prj")                                     # the sidecar that went missing
    zones.rename(columns={"rate": "prescription_rate"}).to_crs(4326).to_file("out/rx_longname.shp")
    with_gap = zones.to_crs(4326)
    with_gap["rate"] = [None] + with_gap["rate"].tolist()[1:]
    with_gap.to_file("out/rx_nan.shp")
    as_text = zones.to_crs(4326)
    as_text["rate"] = as_text["rate"].map(lambda r: f"{r} kg/ha")
    as_text.to_file("out/rx_text.shp")
    simplified = zones.copy()
    simplified["geometry"] = simplified.geometry.simplify(12).buffer(0)
    simplified.to_crs(4326).to_file("out/rx_simplified.shp")

print(f"{inside.sum():,} cells = {zones.area.sum() / 1e4:.2f} ha; {len(zones)} zone polygons; "
      f"zone areas (ha): {zones.dissolve('zone').area.div(1e4).round(2).to_dict()}")

41,126 cells = 20.01 ha; 335 zone polygons; zone areas (ha): {1: 1.83, 2: 8.38, 3: 9.8}

Example 1 — the pre-flight check

import numpy as np, pandas as pd, geopandas as gpd

def preflight(path, rate_col="rate", expect_crs=4326, working_width_m=24,
              min_area_ha=0.1, area_crs=32632):
    rx = gpd.read_file(path)
    problems = []

    if rx.crs is None:
        problems.append("no CRS — a shapefile with no .prj is read as an assumption")
    elif expect_crs and rx.crs.to_epsg() != expect_crs:
        problems.append(f"CRS is EPSG:{rx.crs.to_epsg()}, expected EPSG:{expect_crs}")

    if rate_col not in rx.columns:
        near = [c for c in rx.columns if c.lower().startswith(rate_col[:4].lower())]
        problems.append(f"no '{rate_col}' column; closest are {near or rx.columns.tolist()}")
    else:
        if not pd.api.types.is_numeric_dtype(rx[rate_col]):
            problems.append(f"'{rate_col}' is {rx[rate_col].dtype}, not numeric")
        elif rx[rate_col].isna().any():
            problems.append(f"{int(rx[rate_col].isna().sum())} null rates")

    if (~rx.geometry.is_valid).any():
        problems.append(f"{int((~rx.geometry.is_valid).sum())} invalid geometries")
    if rx.geometry.is_empty.any():
        problems.append(f"{int(rx.geometry.is_empty.sum())} empty geometries")
    if (rx.geom_type == "MultiPolygon").any():
        problems.append(f"{int((rx.geom_type == 'MultiPolygon').sum())} multipart polygons")

    overlaps = len(gpd.sjoin(rx, rx, predicate="overlaps")) // 2
    if overlaps:
        problems.append(f"{overlaps} overlapping polygon pairs")

    if rx.crs is not None:
        m = rx.to_crs(area_crs)
        tiny = int((m.area < min_area_ha * 1e4).sum())
        if tiny:
            problems.append(f"{tiny} of {len(m)} polygons under {min_area_ha} ha")
        narrow = int((m.area / m.length < working_width_m / 4).sum())
        if narrow:
            problems.append(f"{narrow} polygons with area/perimeter under {working_width_m / 4:g} m")

    print(path)
    for p in problems:
        print("  !", p)
    if not problems:
        print("  preflight passed")
    return problems

for path in ["out/rx_good.shp", "out/rx_no_prj.shp", "out/rx_longname.shp", "out/rx_nan.shp",
             "out/rx_text.shp", "out/rx_simplified.shp"]:
    preflight(path)

out/rx_good.shp
  ! 14 overlapping polygon pairs
  ! 326 of 335 polygons under 0.1 ha
  ! 332 polygons with area/perimeter under 6 m
out/rx_no_prj.shp
  ! no CRS — a shapefile with no .prj is read as an assumption
out/rx_longname.shp
  ! no 'rate' column; closest are ['zone', 'prescripti', 'geometry']
  ! 14 overlapping polygon pairs
  ! 326 of 335 polygons under 0.1 ha
  ! 332 polygons with area/perimeter under 6 m
out/rx_nan.shp
  ! 1 null rates
  ! 14 overlapping polygon pairs
  ! 326 of 335 polygons under 0.1 ha
  ! 332 polygons with area/perimeter under 6 m
out/rx_text.shp
  ! 'rate' is str, not numeric
  ! 14 overlapping polygon pairs
  ! 326 of 335 polygons under 0.1 ha
  ! 332 polygons with area/perimeter under 6 m
out/rx_simplified.shp
  ! 80 overlapping polygon pairs
  ! 326 of 335 polygons under 0.1 ha
  ! 329 polygons with area/perimeter under 6 m

Even the "good" file fails, and it should: polygonised straight from the zone raster, it is hundreds of polygons too fine for any implement, and its T-junctions overlap once it is in WGS 84. area / length is a crude width proxy — on raster-edged polygons the stair-stepped perimeter drags it down — but it points at the ribbons and slivers that make a controller oscillate.

Example 2 — repair the common failures

The obvious repair works on the polygons: coerce the rate, simplify each polygon to the implement, drop the small ones, and dissolve by rate.

def repair_prescription(rx, rate_col="rate", crs_out=4326, working_width_m=24,
                        min_area_ha=0.1, area_crs=32632):
    rx = rx.copy()
    rx[rate_col] = pd.to_numeric(rx[rate_col], errors="coerce")
    rx["geometry"] = rx.geometry.buffer(0)
    rx = rx[rx.geometry.notna() & ~rx.geometry.is_empty].explode(index_parts=False)

    m = rx.to_crs(area_crs)
    m["geometry"] = m.geometry.simplify(working_width_m / 2).buffer(0)
    m = m[m.area >= min_area_ha * 1e4]
    m = m.dissolve(by=rate_col, as_index=False).explode(index_parts=False)
    m = m[m.area >= min_area_ha * 1e4]
    return m.to_crs(crs_out)[[rate_col, "geometry"]]


def coverage_summary(rx, field, crs=32632):
    m = rx.to_crs(crs)
    union = m.union_all()
    pairs = gpd.sjoin(m, m, predicate="overlaps")
    return (f"{len(m)} polygons; {len(pairs) // 2} overlapping pairs, "
            f"{max(m.area.sum() - union.area, 0) / 1e4:.2f} ha claimed twice, "
            f"{field.difference(union).area / 1e4:.2f} ha of the field uncovered")

field = zones.union_all()
print("polygon repair:", coverage_summary(repair_prescription(zones.to_crs(4326)), field))

polygon repair: 9 polygons; 6 overlapping pairs, 0.41 ha claimed twice, 1.99 ha of the field uncovered

It trades one set of faults for another. Simplifying each polygon separately moves shared edges apart, so neighbouring zones overlap in some places and leave gaps in others, and the dropped polygons leave holes that the dissolve does not fill. Those holes are exactly where the controller falls back to its default.

Clean the zones where they are still a raster instead. rasterio.features.sieve merges every patch smaller than a threshold into its largest neighbour, so the zones stay a complete, non-overlapping tiling when they are polygonised:

min_cells = int(np.ceil(1000 / abs(transform.a * transform.e)))      # 0.1 ha in 2.2 m cells
sieved = features.sieve(zone_grid, size=min_cells, mask=inside)
sieved_shapes = features.shapes(sieved, mask=inside, transform=transform)
clean = gpd.GeoDataFrame(
    [{"zone": int(v), "geometry": shapely.geometry.shape(g)} for g, v in sieved_shapes], crs=32632)
clean["rate"] = clean["zone"].map(RATES)
clean = clean.set_geometry(                                           # T-junctions, as in step 4
    shapely.snap(clean.geometry.values, shapely.multipoints(shapely.get_coordinates(clean.geometry.values)), 0.001),
    crs=clean.crs)

print("raster sieve:  ", coverage_summary(clean, field))
print(f"{((sieved != zone_grid) & inside).mean() / inside.mean():.1%} of cells moved to a neighbouring zone")

raster sieve:   9 polygons; 0 overlapping pairs, 0.00 ha claimed twice, 0.00 ha of the field uncovered
11.7% of cells moved to a neighbouring zone

The sieve changes the prescription too — about one cell in nine now gets its neighbour's rate — but it says so, and every square metre of the field still has exactly one rate.

Measured in Example 2. The sieve moved 11.7% of cells into a neighbouring zone; the polygon repair left a tenth of the field on the controller default.

Example 3 — write in more than one format

import pathlib

def export_prescription(rx, stem="deliver/prescription", rate_col="rate"):
    out = pathlib.Path(stem).parent
    out.mkdir(parents=True, exist_ok=True)

    # shapefile: the widest-supported interchange, ten-character field names
    rx.rename(columns={rate_col: rate_col[:10]}).to_file(f"{stem}.shp")

    # GeoJSON for inspection and for web tools
    rx.to_file(f"{stem}.geojson", driver="GeoJSON", COORDINATE_PRECISION=7)

    for p in sorted(out.glob(f"{pathlib.Path(stem).name}.*")):
        print(f"  {p.name:24} {p.stat().st_size / 1024:8.1f} KB")
    print("  check the terminal's required format before relying on any of these")

export_prescription(clean.to_crs(4326)[["zone", "rate", "geometry"]])
preflight("deliver/prescription.shp")

  prescription.cpg              0.0 KB
  prescription.dbf              0.4 KB
  prescription.geojson        118.1 KB
  prescription.prj              0.1 KB
  prescription.shp             67.9 KB
  prescription.shx              0.2 KB
  check the terminal's required format before relying on any of these
deliver/prescription.shp
  ! 4 polygons with area/perimeter under 6 m

The one flag left on the delivered file is the width proxy, on four of the nine zones whose stair-stepped raster edges make their perimeters long. Whether those zones are really too narrow for the implement is a question for a map and the machine's specification, not for the proxy.

Shipping a shapefile and a GeoJSON costs nothing and means the operator has an alternative when the first one fails at the gate. COORDINATE_PRECISION=7 keeps about a centimetre, which is finer than any implement can steer.

Text rates and missing CRS information

Why a text rate column produces a constant rate

A terminal reads the attribute table and looks for a numeric field. A column that pandas typed as object — because its values carried units, or a placeholder string — is written to the shapefile as a character field. A terminal that finds no numeric rate typically falls back to the default configured on the machine and applies it everywhere. Nothing errors in Python, and the map looks fine on the screen.

Why missing CRS information puts the field in the wrong place

A shapefile's coordinate system lives in the .prj sidecar. Without it, a terminal assumes something — usually WGS 84 — and a prescription written in UTM or a national grid, with coordinates in the hundreds of thousands or millions, is then interpreted as degrees, far outside any valid range. The field appears empty because the machine is nowhere near it.

Why overlaps scramble the rates

When two polygons cover the same ground, the controller has to choose. Implementations differ — the first, the last, the largest, the highest rate — and nothing in the file says which. The applied pattern then depends on the order features happen to be stored in, which is why an overlap can behave differently on two terminals reading the same file.

Why polygons finer than the response distance cause oscillation

A rate change is a mechanical action with a lag — a valve, a gate, a motor. During that lag the machine travels several metres. A polygon narrower than that distance is entered and left before the rate has changed, so the controller is permanently chasing a target it never reaches, and the applied rate is an average with no relation to the map.

A missing number alone stays numeric (NaN in a float column); units or placeholders in the values make the column text.

Terminal formats, truncation and small test files

Ask for the terminal's format before exporting anything.
Shapefile truncates names to ten characters, with only a Python warning.
Test with one small file on the actual machine.
Units are not in the file. kg/ha and l/ha look identical — and putting them in the values makes the column text.
Some terminals cap the number of polygons. Sieve the zones before polygonising, then dissolve by rate.
A .prj is not optional. Ship the whole shapefile set.
ISOXML is the standard (ISO 11783-10) and is what most modern terminals prefer.
Clean zones as a raster, not as polygons. Per-polygon simplification opens gaps and overlaps between neighbours.
Keep the as-applied data. It is the only evidence of what actually happened.

Internal links

How to write a variable-rate prescription map — building the file correctly
Management zones explained — where the polygons come from
Shapefile column names are truncated — the ten-character limit
How to fix invalid geometries in GeoPandas — repairing what the terminal would skip
How to fix gaps and overlaps in a polygon coverage — making the coverage clean
How to simplify geometry in GeoPandas — simplifying to the implement
Mixed geometry types cause a save error — multipart problems
How to package GIS deliverables in Python — shipping the whole set

FAQ

Why does my prescription apply a constant rate?

Usually the rate column is stored as text — values like 180 kg/ha — or its name was truncated to ten characters by the shapefile driver. A terminal that finds no numeric rate typically uses its default. A missing value on its own does not cause this: pandas keeps it as NaN in a numeric column.

Why does the field appear empty on the terminal?

A missing or wrong CRS. Without a .prj, a prescription in UTM or a national grid is read as degrees and lands outside any valid range — the Ferrara field's eastings of 754,000 are nowhere on Earth as longitudes.

Why are the applied rates scrambled?

Overlapping polygons. Different controllers resolve an overlap differently, and nothing in the file says which one wins.

Why is the machine oscillating between rates?

Polygons narrower than the distance travelled during the actuator's response — commonly 5–15 m. Remove the small patches before polygonising (a raster sieve), rather than simplifying each polygon, which opens gaps and overlaps between zones.

What format should I use?

Whatever the terminal wants. Shapefile is the most widely accepted interchange; ISOXML is the standard modern terminals prefer.

How do I avoid all of this?

Run a file-level pre-flight — CRS, column name and type, validity, single-part, no overlaps, no gaps, minimum size — clean the zones as a raster before polygonising, and test one small file on the actual machine first.