How to clean yield monitor data in Python

Problem statement

Raw yield monitor data is a track log with a yield column, and the artefacts in it are structural rather than random: start and stop transients at the end of every pass, inflated or deflated values wherever the header was not full, a ten-to-twenty-second sensor lag that smears every boundary, and a scattering of physically impossible values.

Cleaning is therefore not outlier removal. It is a sequence of rules, each targeting a known mechanism, applied in an order that matters — and the count removed by each rule is part of the result, because a pipeline that silently discards a quarter of the points has made a decision nobody reviewed.

Quick answer

import numpy as np, pandas as pd

def clean_yield(y, lag_s=12, min_speed=1.5, max_speed=12, trim_points=8,
                full_swath=None, yield_range=(0.5, 25)):
    n0 = len(y)
    log = {}

    y = y.sort_values("timestamp").copy()
    y = add_pass_id(y)
    y = correct_lag(y, lag_s)

    keep = y["from_end"].ge(trim_points) & y["point_in_pass"].ge(trim_points)
    log["start_stop"] = int((~keep).sum()); y = y[keep]

    keep = y["speed_kmh"].between(min_speed, max_speed)
    log["speed"] = int((~keep).sum()); y = y[keep]

    keep = y["yield_lag_corrected"].notna()        # the lag shift has no value to give these
    log["lag_edge"] = int((~keep).sum()); y = y[keep]

    if full_swath is not None:
        keep = y["swath_m"] >= 0.95 * full_swath
        log["partial_width"] = int((~keep).sum()); y = y[keep]

    keep = y["yield_lag_corrected"].between(*yield_range)
    log["impossible"] = int((~keep).sum()); y = y[keep]

    print(f"{n0:,} → {len(y):,} points ({len(y)/n0:.1%} kept)")
    print(pd.Series(log).to_string())
    return y, log

Every rule reports, and each rule reports only what it removed. The first version of this function had no lag_edge line: shifting yield back by the lag leaves the last seconds of each pass with nothing to shift in, those points became NaN, and the range rule silently absorbed them — on the Ferrara harvest, 224 of the 225 points it attributed to "impossible" values were those gaps. If partial_width removes 30% of the points, that is a fact about the field's shape, not a bug — but it is a fact that belongs in the output, under its own name.

18,444 points; median speed 5.7 km/h; yield 0.02–18.0 t/ha, 1,367 points at the maximum

Step-by-step solution

1. Split by harvest date and by machine first

A file can contain two harvests or two combines. Each has its own calibration, its own lag and its own swath, and mixing them makes every later rule wrong.

2. Reconstruct the passes

A pass ends where there is a time gap or a large heading change. Everything afterwards — trimming, headland detection, lag correction — operates within a pass.

3. Correct the lag before any spatial rule

The lag displaces values along the track by tens of metres. Applying a spatial filter first smooths the displaced values and then displaces them again, which makes the error worse rather than better.

4. Trim the start and stop of every pass

Flow builds and decays over several seconds. Dropping the first and last six to ten points of each pass removes the transients; the exact number depends on the logging interval and the machine.

5. Filter speed

Very low speeds give division by a small number; very high speeds are transport between fields. Both produce spikes.

6. Handle partial width

Where the logged swath is below the full header width, the area is wrong and so is the yield. If the true swath is logged, correct it; if it is not, remove the points.

7. Remove physically impossible values, then statistical outliers

A negative yield or 60 t/ha of wheat is impossible. After the impossible values are gone, a local outlier rule — a median absolute deviation within a moving neighbourhood — catches the rest without being dominated by them.

8. Scale to the weighbridge

The sum of the cleaned map should match the delivered tonnage. Applying a single scale factor to match it is standard practice and is the only external check available.

Code examples

The data used on this page

The 2017 maize harvest of a 22 ha field at Porto Felloni near Ferrara, published on Zenodo (CC BY 4.0) — the field the yield monitor explainer takes apart. The monitor logged one point a second from a combine with a 6 m header, and the authors had already removed border effects. The file has a record number, easting, northing, yield in kg/ha and moisture, and nothing else: no timestamp, speed, heading or swath. So the timestamp here is the record number in seconds, heading and speed come from consecutive positions, and the coordinates are read as the Italian Gauss-Boaga grid they really are, whatever the file declares.

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

raw = gpd.read_file("FIELD_YIELD/yield2017/yield2017.shp")
y = gpd.GeoDataFrame(raw.drop(columns="geometry"), geometry=gpd.points_from_xy(raw["EST"], raw["NORD"]),
                     crs="EPSG:3003")
y["timestamp"] = pd.to_datetime(y["NUMBER"], unit="s")          # 1 Hz logging
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
print(f"{len(y):,} points; median speed {y['speed_kmh'].median():.1f} km/h; "
      f"yield {y['yield_t_ha'].min():.2f}–{y['yield_t_ha'].max():.1f} t/ha, "
      f"{(y['yield_t_ha'] == y['yield_t_ha'].max()).sum():,} points at the maximum")

18,444 points; median speed 5.7 km/h; yield 0.02–18.0 t/ha, 1,367 points at the maximum

Example 1 — passes, lag and trimming

import numpy as np, pandas as pd

def add_pass_id(y, gap_s=10, turn_deg=60):
    y = y.sort_values("timestamp").copy()
    dt = y["timestamp"].diff().dt.total_seconds()
    turn = y["heading_deg"].diff().abs().mod(360)
    turn = np.minimum(turn, 360 - turn)
    y["pass_id"] = ((dt > gap_s) | (turn > turn_deg)).cumsum()
    y["point_in_pass"] = y.groupby("pass_id").cumcount()
    y["from_end"] = (y.groupby("pass_id")["point_in_pass"].transform("max")
                     - y["point_in_pass"])
    return y

def correct_lag(y, lag_s=12.0, col="yield_t_ha"):
    out = []
    for _, g in y.groupby("pass_id", sort=False):
        g = g.sort_values("timestamp").copy()
        t = (g["timestamp"] - g["timestamp"].iloc[0]).dt.total_seconds().values
        g["yield_lag_corrected"] = np.interp(t, t - lag_s, g[col].values,
                                             left=np.nan, right=np.nan)
        out.append(g)
    return pd.concat(out)

passes = add_pass_id(y)
sizes = passes.groupby("pass_id").size()
print(f"{sizes.size} passes; {(sizes >= 200).sum()} of 200+ points hold "
      f"{sizes[sizes >= 200].sum():,} of {len(y):,} points; {(sizes == 1).sum()} are a single point")

144 passes; 43 of 200+ points hold 16,847 of 18,444 points; 52 are a single point

The single-point "passes" are the turns at the ends of the field: the heading swings more than 60° between two fixes, and the rule counts each swing as a pass of its own. They are harmless here, because the trimming rule removes them, but they are why pass counts from this function should not be quoted as the number of times the combine crossed the field.

The quick answer's pipeline, run on the harvest — there is no swath column, so the partial-width rule is skipped, and the lag is the 12 s default because this field offers no boundary to estimate it from (Example 2):

clean, log = clean_yield(y, lag_s=12, full_swath=None, yield_range=(0.5, 25))

18,444 → 16,952 points (91.9% kept)
start_stop    1201
speed           66
lag_edge       224
impossible       1

Example 2 — estimate the lag rather than assuming it

import numpy as np

def estimate_lag(y, boundary_geom, candidates=np.arange(0, 25, 1.0)):
    """Choose the lag that makes the yield step align with a known boundary."""
    best, best_score = None, -np.inf
    inside = y.geometry.within(boundary_geom)
    for lag in candidates:
        shifted = correct_lag(y, lag)["yield_lag_corrected"]
        a, b = shifted[inside.values], shifted[~inside.values]
        a, b = a[np.isfinite(a)], b[np.isfinite(b)]
        if len(a) < 50 or len(b) < 50:
            continue
        # t-like separation: a sharp boundary maximises it
        score = abs(a.mean() - b.mean()) / np.sqrt(a.var() / len(a) + b.var() / len(b))
        if score > best_score:
            best, best_score = lag, score
    print(f"best lag {best:.1f} s (separation {best_score:.1f})")
    return best

A trial strip, a variety change or a field boundary all give a known step. The Ferrara field has none of them in its published data — the authors removed the borders, and there is no strip or variety change — so this function has nothing to run on here, and the explainer's attempt with opposite-direction passes gave a flat curve with no clear best lag. That is the normal situation, and it is the reason to run a deliberate test on a known edge once per combine: without one, the lag is a default nobody has checked.

Example 3 — local outliers and the weighbridge scale

import numpy as np, geopandas as gpd
from scipy.spatial import cKDTree

def local_outliers(y, col="yield_lag_corrected", k=20, mad_threshold=3.5):
    xy = np.c_[y.geometry.x, y.geometry.y]
    tree = cKDTree(xy)
    _, idx = tree.query(xy, k=k + 1, workers=-1)
    v = y[col].values
    neigh = v[idx[:, 1:]]
    med = np.nanmedian(neigh, axis=1)
    mad = np.nanmedian(np.abs(neigh - med[:, None]), axis=1)
    score = 0.6745 * (v - med) / np.where(mad == 0, np.nan, mad)
    flag = np.abs(score) > mad_threshold
    print(f"local outliers: {np.nansum(flag):,} of {len(y):,} "
          f"({np.nansum(flag)/len(y):.1%})")
    return y[~np.nan_to_num(flag, nan=False)]

def scale_to_weighbridge(y, delivered_t, area_ha, col="yield_lag_corrected"):
    mapped = y[col].mean() * area_ha
    factor = delivered_t / mapped
    print(f"map total {mapped:.1f} t, delivered {delivered_t:.1f} t, "
          f"scale factor {factor:.3f}")
    y = y.copy()
    y["yield_scaled"] = y[col] * factor
    return y, factor

kept = local_outliers(clean)
ceiling = kept["yield_t_ha"].max()
print(f"after cleaning, {(kept['yield_t_ha'] == ceiling).sum():,} points still sit at the {ceiling:.1f} t/ha maximum")

local outliers: 414 of 16,952 (2.4%)
after cleaning, 1,262 points still sit at the 18.0 t/ha maximum

The last line is a limit of every rule on this page. The monitor clipped yield at 18 t/ha, so the points at the ceiling are not outliers — they are the high-yielding part of the field with its values flattened — and neither a range rule nor a local one will remove them. They need flagging as censored, not deleting.

The weighbridge scale needs the delivered tonnage, which the published data does not include, so it cannot be run here:

y_scaled, factor = scale_to_weighbridge(kept, delivered_t=..., area_ha=22)

A scale factor far from 1 — outside roughly 0.9 to 1.1 — means the calibration or the swath assumption is wrong, and scaling it away hides a problem that will recur next year.

Explanation

Why the order of the rules matters

The lag correction moves values along the track, so it must come before anything that depends on position — local outlier detection, headland masking, rasterising. Pass reconstruction must come before trimming, because "the first eight points" is meaningless without a pass. Removing impossible values before the statistical rule stops a handful of 400 t/ha points from setting the threshold for everything else.

Why partial-width points are removed rather than corrected

The monitor divides measured mass by assumed area. If the true swath was logged, the correction is arithmetic. If it was not — which is common on older monitors — there is no information from which to recover the area, and the value is simply wrong by an unknown factor. Removing them loses coverage along one edge of the field and keeps the rest honest.

Why a local rule beats a global one

Yield varies across a field by a factor of two or more, so a global threshold either keeps outliers in the low-yielding part or removes real data from the high-yielding part. Comparing each point with its twenty nearest neighbours adapts to the local level, and the median absolute deviation makes it robust to the outliers it is looking for.

Why the weighbridge is the only external check

Everything else in a yield map is internally consistent by construction: the monitor's calibration, the swath assumption and the moisture correction all affect the whole map together, so no internal statistic reveals them. The delivered tonnage is measured independently, and the ratio between it and the map total is the one number that tests the whole chain.

Edge cases or notes

• Split by date and machine first. Everything downstream assumes one of each.
• The lag is a machine property. Estimate it once, reuse it.
• Headlands may deserve their own treatment rather than deletion.
• Overlap between passes double-counts. Detect it geometrically if the swath is logged.
• Moisture basis must be consistent before any comparison between passes.
• Keep the removed points. They are the evidence for the cleaning report.
• A scale factor outside 0.9–1.1 is a warning, not a correction.
• Report the kept fraction. Whatever it is, it should be visible, rule by rule.
• A clipped maximum is not an outlier. Points piled at the monitor's ceiling are censored values; flag them.

Internal links

• Yield monitor data explained: what the numbers really are — where each artefact comes from
• A yield map has impossible values and stripes — diagnosing before cleaning
• How to remove spatial outliers in Python — the local rule in general
• Management zones explained — what the cleaned map feeds
• How to clean a GPS track in Python — pass reconstruction
• How to calculate speed, distance and direction along a track — the speed filter
• How to analyse an on-farm strip trial in Python — where the lag limits resolution
• Repair, reject or flag — deciding what to do with bad rows

FAQ

How much yield data should I expect to remove?

It depends on the file, which is why each rule has to report. On the Ferrara harvest, which the authors had already cleaned of border effects, the pipeline still removed the share shown in Example 1 — mostly start-and-stop trimming of points that were no longer transients. Report the count removed by each rule, and check that each rule is removing what it is meant to.

What order should the cleaning rules go in?

Split by date and machine, reconstruct passes, correct the lag, trim starts and stops, filter speed, handle partial width, remove impossible values, then local outliers.

How do I find the sensor lag?

Sweep candidate lags and pick the one that maximises the separation across a known boundary — a trial strip, a variety change or a field edge.

Should I correct or remove partial-width points?

Correct them if the true swath was logged. If it was not, the area is unrecoverable and the points have to go.

Why use a local outlier rule?

Because yield varies across a field by a factor of two, so a global threshold removes real data in the good areas and keeps outliers in the poor ones.

What does the weighbridge check tell me?

Whether the whole chain — calibration, swath, moisture — is right. A scale factor outside about 0.9 to 1.1 is a problem to fix, not to scale away.