A yield map has impossible values and stripes

Problem statement

The map shows a flat plateau of identical maximum values, a ring of low yields around the headland, passes that do not match their neighbours and a drift from one end of the day to the other. None of it is agronomy, and each pattern has a specific mechanical cause.

The distinction that matters is between values that are impossible and values that are merely wrong. An impossible value — negative, or several times the crop's physical maximum — is a broken record. A value clipped at the monitor's ceiling is censored: real, high, and flattened. A stripe is a systematic error in a term of the yield calculation, and removing the points does not fix the bias in the rest of the pass.

This page diagnoses a real yield monitor file: 18,444 one-second records from a combine harvesting 6 m of maize at a time in a field near Ferrara, Italy, in 2017, published with the study by Kayad and colleagues on Zenodo.

Quick answer

Diagnose before cleaning, because the pattern names the cause:

print(y["yield_t_ha"].describe(percentiles=[.001, .01, .5, .99, .999]).round(2).to_dict())
print(f"negative or zero: {(y.yield_t_ha <= 0).sum():,}")
print(f"above 3x the median: {(y.yield_t_ha > 3 * y.yield_t_ha.median()).sum():,}")
print(f"exactly at the maximum, {y.yield_t_ha.max()} t/ha: {(y.yield_t_ha == y.yield_t_ha.max()).sum():,}")
print(f"speed below 1 km/h: {(y.speed_kmh < 1).sum():,}; above 12 km/h: {(y.speed_kmh > 12).sum():,}")
print("swath recorded per point:", "swath_m" in y)

# is the variation between passes or within them? (full passes only)
full = y[y["pass_len"] >= 20]
between = full.groupby("pass_id")["yield_t_ha"].mean().var()
within = full.groupby("pass_id")["yield_t_ha"].var().mean()
print(f"between-pass variance {between:.2f}, within-pass {within:.2f}, ratio {between / within:.2f}")

{'count': 18444.0, 'mean': 14.49, 'std': 2.84, 'min': 0.02, '0.1%': 0.17, '1%': 3.03, '50%': 14.97, '99%': 18.0, '99.9%': 18.0, 'max': 18.0}
negative or zero: 0
above 3x the median: 0
exactly at the maximum, 18.0 t/ha: 1,367
speed below 1 km/h: 12; above 12 km/h: 120
swath recorded per point: False
between-pass variance 2.42, within-pass 5.43, ratio 0.45

No record here is impossible in the negative-or-absurd sense. What the file has instead is a ceiling: 1,367 points — 7.4% of the harvest — sit at exactly 18.0 t/ha, because the monitor clips there. The between-to-within ratio says the full passes differ from each other by about half as much as the yield varies along them; whether that is field or machine is what the stripe test below tries to separate. The file records no swath width, so the classic partial-width stripe cannot be checked from this data at all.

Impossible values, speed extremes, pass stripes

1. Impossible values first

Negative yield, zero yield with a non-zero flow, and anything above about three times the crop's realistic maximum. These are broken records and they distort every statistic used to find the subtler problems. Just as important is the value that is not impossible but is not real either: a ceiling.

ceiling = y["yield_t_ha"].max()
share = full.groupby("pass_id")["yield_t_ha"].agg(lambda s: (s == ceiling).mean())
print(f"points at {ceiling} t/ha by year:", {yr: int((d.yield_t_ha == d.yield_t_ha.max()).sum()) for yr, d in years.items()})
print(f"share of a full pass at the ceiling: median {share.median():.1%}, worst pass {share.max():.1%}")

points at 18.0 t/ha by year: {2016: 356, 2017: 1367, 2018: 1424}
share of a full pass at the ceiling: median 3.6%, worst pass 41.8%

All three years were clipped at the same 18.0 t/ha, and in the worst pass of 2017 two fifths of the points sit there. They are the high-yielding part of the field with its values flattened, so deleting them biases the map downwards; flag them as censored instead.

2. Very low or very high speed

Yield is mass flow divided by speed times width. At 0.2 km/h the divisor is almost zero and the yield is enormous; above about 12 km/h the machine is usually in transport rather than harvesting. The file records no speed, so it is derived from successive positions — 12 points below 1 km/h and 120 above 12 km/h in 2017.

3. Stripes aligned with the passes

Three causes, distinguishable by their pattern:

The moisture reading is in this file, so the second can be looked at:

per_pass = full.groupby("pass_id").agg(yield_t_ha=("yield_t_ha", "mean"), moisture=("moisture_pct", "mean"),
                                       start=("timestamp", "min"))
hourly = y.set_index("timestamp")[["yield_t_ha", "moisture_pct"]].resample("1h").mean().round(1)
hourly.index = [f"hour {i + 1}" for i in range(len(hourly))]
print(hourly.T.to_string())
print(f"pass mean yield vs moisture: r = {per_pass.yield_t_ha.corr(per_pass.moisture):.2f}; "
      f"moisture vs time: r = {per_pass.moisture.corr(per_pass.start.astype('int64')):.2f}")

              hour 1  hour 2  hour 3  hour 4  hour 5  hour 6
yield_t_ha      13.7    13.7    14.1    14.3    15.4    15.6
moisture_pct    18.9    18.8    21.0    20.3    22.6    22.7
pass mean yield vs moisture: r = 0.47; moisture vs time: r = 0.70

Over the six hours of logging — record numbers at one a second — the grain moisture reading rose from about 19% to 23%, and the mapped yield rose with it. That could be the field — the combine worked across it through the day — or it could be a moisture correction applied to a reading that drifted. The data alone cannot say; the pattern is what tells you to go and ask whether the moisture sensor was checked. It is also a warning about rules of thumb: grain moisture is often said to fall through a harvest day, and here it went the other way.

4. A halo around the headland

Start and stop transients, turning, and partial width all concentrate at the ends of passes:

by_position = full.groupby(pd.cut(full["point_in_pass"], [-1, 2, 5, 10, 20, 50, 10_000],
                                  labels=["1-3", "4-6", "7-11", "12-21", "22-51", "52+"]))["yield_t_ha"].median()
print("median t/ha by position in the pass:", by_position.round(1).to_dict())

median t/ha by position in the pass: {'1-3': 9.7, '4-6': 12.0, '7-11': 13.3, '12-21': 14.6, '22-51': 15.3, '52+': 15.0}

The first three seconds of every full pass have a median of 9.7 t/ha against about 15 in the middle — a ring of low values around the field that is the machine filling up, not the crop thinning out. Reconstruct the passes and trim the first and last few points of each.

5. Smearing along the direction of travel

The grain takes time to travel from the header to the flow sensor, so each value is recorded some seconds after the crop was cut, which displaces it along the direction of travel — at this field's median 5.7 km/h, each second of lag is about 1.6 m. It shows as a boundary that is blurred in the direction of travel and sharp across it. How to clean yield monitor data estimates and corrects the lag for this same file.

6. Two populations in the histogram

A bimodal distribution usually means two combines, two harvest dates, or two crops in one file. Split before cleaning. The three Ferrara harvests are one crop in one field with similar medians — 14.7, 15.0 and 14.3 t/ha — so stacking them adds points without adding a population.

7. Check the total against the weighbridge

The one external check. A map total that differs from the delivered tonnage by more than about 10% means the calibration, the swath assumption or the moisture basis is wrong — a bias that cleaning will not remove. The published data does not include the delivered tonnage, so the check cannot be run for this field; Example 3 shows it.

A diagnostic report on a real yield dataset

The data used on this page

The FIELD_YIELD dataset from Kayad et al. (2019), Zenodo record 14216165 (CC BY 4.0): maize yield monitor points for a 22 ha field at Porto Felloni near Ferrara, 2016–2018, with a record number, position, yield in kg/ha and grain moisture. The monitor logged one point a second, so the record number stands in for time; per the paper, the authors had already removed some border effects. The coordinates are Italian Gauss–Boaga (EPSG:3003) whatever the files' own labels say; the setup derives time, heading, speed and pass numbers the same way as the yield monitor guide.

import os
import zipfile

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

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


def load(path):
    raw = gpd.read_file(path)
    y = gpd.GeoDataFrame(raw.drop(columns="geometry"), geometry=gpd.points_from_xy(raw["EST"], raw["NORD"]),
                         crs="EPSG:3003").sort_values("NUMBER").reset_index(drop=True)
    y["timestamp"] = pd.to_datetime(y["NUMBER"], unit="s")             # record number: one a second
    dx, dy = y["EST"].diff(), y["NORD"].diff()
    y["heading_deg"] = np.degrees(np.arctan2(dx, dy)) % 360
    y["speed_kmh"] = np.hypot(dx, dy) * 3.6 / y["NUMBER"].diff().clip(lower=1)
    y["yield_t_ha"] = y["YIELD"] / 1000
    y["moisture_pct"] = y["HUMID"]
    turn = y["heading_deg"].diff().abs().mod(360)
    y["pass_id"] = ((y["NUMBER"].diff() > 10) | (np.minimum(turn, 360 - turn) > 60)).cumsum()
    y["point_in_pass"] = y.groupby("pass_id").cumcount()
    y["pass_len"] = y.groupby("pass_id")["point_in_pass"].transform("size")
    return y


years = {2016: load("FIELD_YIELD/yield2016/yield2016.gpkg"),
         2017: load("FIELD_YIELD/yield2017/yield2017.shp"),
         2018: load("FIELD_YIELD/yield2018/yield2018.shp")}
y = years[2017]
print({yr: f"{len(d):,} points, median {d['yield_t_ha'].median():.1f} t/ha" for yr, d in years.items()})
print(f"2017: {y['pass_id'].nunique()} passes, {y.loc[y['pass_len'] >= 20, 'pass_id'].nunique()} of them 20 s or longer; "
      f"median speed {y['speed_kmh'].median():.1f} km/h; {(y['timestamp'].max() - y['timestamp'].min())} of logging")

{2016: '21,169 points, median 14.7 t/ha', 2017: '18,444 points, median 15.0 t/ha', 2018: '22,616 points, median 14.3 t/ha'}
2017: 144 passes, 59 of them 20 s or longer; median speed 5.7 km/h; 0 days 05:44:00 of logging

Example 1 — the diagnostic report

def diagnose_yield(y, crop_max=20.0, min_pass=20):
    v = y["yield_t_ha"]
    r = {"points": len(y),
         "negative_or_zero": int((v <= 0).sum()),
         "above_crop_max": int((v > crop_max).sum()),
         "at_the_maximum": int((v == v.max()).sum()),
         "maximum": float(v.max()),
         "median": round(float(v.median()), 2)}
    if "speed_kmh" in y:
        r["speed_under_1"] = int((y.speed_kmh < 1).sum())
        r["speed_over_12"] = int((y.speed_kmh > 12).sum())
    if "swath_m" in y:
        r["partial_swath_share"] = round(float((y.swath_m < 0.95 * y.swath_m.max()).mean()), 3)
    if "pass_id" in y:
        full = y[y.groupby("pass_id")["pass_id"].transform("size") >= min_pass]
        between = float(full.groupby("pass_id")[v.name].mean().var())
        within = float(full.groupby("pass_id")[v.name].var().mean())
        r["between_within_ratio"] = round(between / within, 3) if within else None
        if "moisture_pct" in y:
            m = full.groupby("pass_id")["moisture_pct"].mean()
            r["moisture_range_between_passes"] = round(float(m.max() - m.min()), 2)
    if "timestamp" in y:
        hourly = y.set_index("timestamp")[v.name].resample("1h").mean()
        r["largest_hourly_step"] = round(float(hourly.diff().abs().max()), 2)
    return r


for year, data in years.items():
    print(year, diagnose_yield(data))

2016 {'points': 21169, 'negative_or_zero': 0, 'above_crop_max': 0, 'at_the_maximum': 356, 'maximum': 18.0, 'median': 14.66, 'speed_under_1': 22, 'speed_over_12': 95, 'between_within_ratio': 0.455, 'moisture_range_between_passes': 4.69, 'largest_hourly_step': 1.47}
2017 {'points': 18444, 'negative_or_zero': 0, 'above_crop_max': 0, 'at_the_maximum': 1367, 'maximum': 18.0, 'median': 14.97, 'speed_under_1': 12, 'speed_over_12': 120, 'between_within_ratio': 0.447, 'moisture_range_between_passes': 6.05, 'largest_hourly_step': 1.09}
2018 {'points': 22616, 'negative_or_zero': 0, 'above_crop_max': 0, 'at_the_maximum': 1424, 'maximum': 18.0, 'median': 14.32, 'speed_under_1': 24, 'speed_over_12': 96, 'between_within_ratio': 0.483, 'moisture_range_between_passes': 3.1, 'largest_hourly_step': 1.63}

The same report for three harvests of one field shows what is the machine and what is the season: the 18.0 t/ha ceiling is there every year, and the speed and moisture figures move from year to year. With a 20 t/ha crop maximum the ceiling never registers as "above the maximum" — only the count of points at the exact maximum finds it.

Example 2 — separate a stripe from a real pattern

def stripe_test(y, value="yield_t_ha", pass_col="pass_id", min_pass=20):
    """How much do neighbouring passes' means agree, in harvest order and across the field?"""
    full = y[y.groupby(pass_col)[pass_col].transform("size") >= min_pass]
    g = full.groupby(pass_col)
    means = g[value].mean()
    heading = np.radians(full["heading_deg"].median())
    across = g["EST"].mean() * np.cos(heading) - g["NORD"].mean() * np.sin(heading)
    result = {"passes": len(means), "between_pass_sd": round(float(means.std()), 2),
              "lag1_harvest_order": round(float(means.loc[g["timestamp"].min().sort_values().index].autocorr(1)), 3),
              "lag1_across_field": round(float(means.loc[across.sort_values().index].autocorr(1)), 3)}
    print(result)
    return result


for year, data in years.items():
    print(year, end=" ")
    stripe_test(data)

2016 {'passes': 61, 'between_pass_sd': 1.31, 'lag1_harvest_order': 0.37, 'lag1_across_field': -0.025}
2017 {'passes': 59, 'between_pass_sd': 1.56, 'lag1_harvest_order': 0.278, 'lag1_across_field': 0.072}
2018 {'passes': 69, 'between_pass_sd': 1.85, 'lag1_harvest_order': 0.502, 'lag1_across_field': 0.135}

A real field pattern varies smoothly across adjacent passes, so the pass means are autocorrelated when the passes are put in order across the field; a machine artefact is independent from pass to pass. Here the correlation of neighbouring passes across the field is weak in every year, and stronger in harvest order — which is what a drift through the day, like the moisture one above, would produce. Neither number is a verdict on its own: a combine that harvests in lands rather than back and forth breaks the link between harvest order and position, and a pass that crosses two soil types is an average of both. Treat the test as a prompt to look at the map, not as a classifier.

Example 3 — the weighbridge check

def weighbridge_check(y, delivered_t, area_ha, value="yield_t_ha"):
    mapped = float(y[value].mean() * area_ha)
    factor = delivered_t / mapped
    verdict = ("good" if 0.95 <= factor <= 1.05
               else "acceptable" if 0.90 <= factor <= 1.10
               else "the calibration, swath or moisture basis is wrong")
    print(f"map {mapped:.1f} t, delivered {delivered_t:.1f} t, factor {factor:.3f} — {verdict}")
    return factor

weighbridge_check(y, delivered_t=..., area_ha=...)      # the tonnage is not in the published data

A factor well away from 1 is not a scaling problem to be divided away; it means the machine was measuring something other than what was delivered, and the same error will recur next year unless the calibration is fixed.

Why a ceiling is worse than an outlier

Why a ceiling is worse than an outlier

An outlier is one bad point among good ones, and any robust filter finds it. A ceiling changes every high value in the same way: the points at 18.0 t/ha are real harvest, but their true values are anywhere above that. A mean over a field with 7% of its points clipped is biased low, a map of the best areas is flat-topped, and no filter can recover the lost values. The honest record says how many points are censored.

Why partial width produces a one-sided stripe

The monitor divides measured mass by an assumed area. On a half-width pass the true area is half the assumed one, so the recorded yield is half the true value. That is a 50% error on one pass — far larger than any sensor error — and it appears as a single consistently low stripe rather than as scattered noise, because it affects the whole pass. A file that does not record the working width, like this one, cannot show it except as a pass that is low end to end.

Why moisture produces stripes and trends

Yield is reported at a standard moisture. If the correction differs between passes — a drifting sensor, or a manual entry changed part way through the day — each pass carries a different multiplicative offset. Since passes are harvested in order, a drift appears as a trend with the harvest sequence rather than with the field.

Why the headland halo is unavoidable without trimming

At the end of every pass the machine slows, stops, turns and restarts. Flow decays and rebuilds over several seconds, the swath is partial on the turn, and the position is changing quickly. All of the affected points are at the ends of passes, so they form a ring. No statistical filter can distinguish them from real headland variation; only the pass geometry can.

Why the weighbridge check cannot be skipped

Calibration error, a wrong swath setting and a wrong moisture basis all scale the entire map. Nothing internal to the map reveals them, because every internal statistic is affected equally. The delivered tonnage is the only independent measurement in the whole system.

Crop maxima and zero-flow points

• Split by date and machine first. A bimodal histogram is usually two of something.
• Crop maximum is crop-specific. 20 t/ha is absurd for wheat and low for sugar beet.
• Zero yield with zero flow is correct — the machine was not harvesting.
• Real headland effects exist. Trimming removes the artefact and the signal together.
• Moisture basis must be stated. Dry-matter and standard-moisture yields differ by 15%.
• A stripe that repeats every year is the machine, not the field.
• Keep the removed points as evidence.
• Fix the calibration rather than scaling the map, or it recurs.

Internal links

• How to clean yield monitor data in Python — the cleaning pipeline
• Yield monitor data explained: what the numbers really are — where each artefact comes from
• How to remove spatial outliers in Python — the local rule
• Management zones explained — why uncleaned yield makes bad zones
• How to clean a GPS track in Python — reconstructing the passes
• Track has impossible speeds — the speed filter
• Repair, reject or flag — deciding what to do with bad rows
• How to build a data cleaning report in GeoPandas — reporting what was removed

FAQ

Why does my yield map have negative values?

Broken records — a sensor fault, a zero-speed division, or a calibration with a negative intercept. Remove them before computing anything.

What causes stripes down a yield map?

Partial header width recorded as full (one low stripe), an inconsistent moisture correction (alternating stripes), or a calibration change mid-harvest (a step).

How do I tell a stripe from a real field pattern?

Check whether the pass means are autocorrelated in spatial order. A real pattern varies smoothly between adjacent passes; an artefact does not.

Why is the headland always brighter or darker?

Start and stop transients, turning and partial width all concentrate there. Trim the first and last few points of each pass.

My map total is 25% above the weighbridge. Should I scale it?

No. That size of error means the calibration, the swath setting or the moisture basis is wrong, and scaling hides a problem that will recur.

Why is my histogram bimodal?

Usually two combines, two harvest dates or two crops in one file. Split before cleaning.