Parsing ISOXML and Shapefile Field Boundaries: A Reference & Troubleshooting Guide

Problem statement

A field boundary looks like a solved problem until two exports of the same field disagree. A task controller writes the field as ISO 11783-10 ISOXML (TASKDATA.XML), a desktop GIS exports the same field as an ESRI shapefile, and the two geometries differ by tens of metres — or one of them lands in the Gulf of Guinea. Boundary decoding is the step inside Equipment Telemetry Parsing that resolves every raw GPS fix to a real field zone, so a boundary that is silently wrong misattributes an entire chemical pass to the wrong acreage and only surfaces during an audit.

Three concrete corruption scenarios recur, and this page exists to make each one impossible to commit unnoticed:

• Undefined coordinate reference system. A shapefile ships without its .prj sidecar, or an ISOXML PNT block carries locally projected metres while the schema promises WGS84 (EPSG:4326) degrees. Interpreting metres as degrees renders the field thousands of kilometres from its true centroid; the geometry is still is_valid, so nothing downstream complains.
• Broken ring winding on multipart fields. ISOXML serializes a field split by a waterway as disjoint polygon parts. Merge them without enforcing the OGC Simple Features winding rule — counter-clockwise exterior, clockwise holes — and the union produces self-intersecting rings that later raise TopologyException deep inside prescription generation.
• Precision compression drift. An OEM writes vertices at five decimal places to shrink the archive. At that resolution a WGS84 vertex can shift by roughly 0.11 m, which is enough to trip false overlap flags when the boundary feeds the Buffer Zone Calculations engine.

The remedy is a single deterministic normalization path: refuse any geometry whose CRS is undefined, reproject to WGS84 with a fixed axis order, union multipart parts while repairing validity, enforce OGC winding, and quarantine anything below a precision floor rather than snapping it into place. The canonical field record the result is written into is governed by Field Schema Design; the boundary itself is never allowed to define its own compliance fields.

Parameter reference table

Every value below changes whether a decoded boundary faithfully reproduces the field on the ground. Recommended values assume WGS84 storage and intermittently connected edge controllers.

Hardcode min_precision_dp and the winding rule in a version-stamped registry rather than scattering literals through controller code, so a single parser upgrade re-normalizes every historical boundary the same way.

Runnable implementation

The module below decodes either format into one canonical, WGS84, OGC-wound FieldBoundary. It targets Python 3.10+ (uses match/case and X | None unions), is fully typed, and depends on pyshp, lxml, pyproj, and shapely. The ISOXML attribute codes (PFD, PLN, LSG, PNT, and the A/C/D attributes) follow ISO 11783-10; note that TASKDATA.XML has no XML namespace, so the tags are matched bare.

from __future__ import annotations

import logging
from collections.abc import Iterator
from dataclasses import dataclass
from pathlib import Path

import shapefile  # pyshp
from lxml import etree
from pyproj import CRS, Transformer
from shapely.geometry import MultiPolygon, Polygon, shape
from shapely.geometry.polygon import orient
from shapely.ops import transform, unary_union
from shapely.validation import make_valid

logger = logging.getLogger("boundary.normalize")

WGS84 = CRS.from_epsg(4326)
MIN_PRECISION_DP = 6  # source decimals below this imply >0.11 m compression drift


@dataclass(frozen=True)
class FieldBoundary:
    field_id: str
    geometry: MultiPolygon  # always WGS84, OGC-wound (CCW exterior, CW holes)
    source_epsg: int
    repaired: bool


def _below_precision(poly: Polygon, min_dp: int) -> bool:
    """True if no coordinate carries more than (min_dp - 1) real decimals."""
    coords = list(poly.exterior.coords)
    for hole in poly.interiors:
        coords.extend(hole.coords)
    # If rounding to one place shy of the floor never changes a value, the
    # source was compressed below the precision we require -> quarantine it.
    return all(round(v, min_dp - 1) == v for xy in coords for v in xy)


def normalize_boundary(field_id: str, parts: list[Polygon], source: CRS) -> FieldBoundary:
    """Union multipart rings into one WGS84, OGC-wound geometry, or raise."""
    if source.to_epsg() is None:
        raise ValueError(f"{field_id}: source CRS is undefined; refusing to ingest")
    if not parts:
        raise ValueError(f"{field_id}: no polygon parts decoded")

    # 1. Precision guard runs on SOURCE coords, before any lossy reprojection.
    if all(_below_precision(p, MIN_PRECISION_DP) for p in parts):
        raise ValueError(f"{field_id}: coordinates below {MIN_PRECISION_DP}-dp floor")

    # 2. Reproject to WGS84 with a fixed lon/lat axis order.
    if source != WGS84:
        project = Transformer.from_crs(source, WGS84, always_xy=True).transform
        parts = [transform(project, p) for p in parts]

    # 3. Repair each part, then union disjoint parts into one geometry.
    repaired = any(not p.is_valid for p in parts)
    merged = unary_union([make_valid(p) for p in parts])

    # 4. Enforce OGC winding: sign=1.0 => CCW exterior, CW interior holes.
    match merged:
        case Polygon():
            polys = [orient(merged, sign=1.0)]
        case MultiPolygon():
            polys = [orient(p, sign=1.0) for p in merged.geoms]
        case _:
            raise ValueError(f"{field_id}: union produced non-areal geometry")

    logger.info(
        '{"event":"boundary_normalized","field_id":"%s","parts":%d,'
        '"source_epsg":%d,"repaired":%s}',
        field_id, len(polys), source.to_epsg(), str(repaired).lower(),
    )
    return FieldBoundary(field_id, MultiPolygon(polys), source.to_epsg(), repaired)


def _polygon_from_pln(pln: etree._Element) -> Polygon | None:
    """Build one Polygon from an ISOXML PLN: LSG type 1 = exterior, 2 = hole."""
    exterior: list[tuple[float, float]] = []
    holes: list[list[tuple[float, float]]] = []
    for lsg in pln.iterfind("LSG"):
        # PNT: C = PointNorth (latitude), D = PointEast (longitude); shapely = (x, y).
        pts = [(float(p.get("D")), float(p.get("C"))) for p in lsg.iterfind("PNT")]
        if len(pts) < 4:  # a closed ring needs at least four coordinates
            continue
        match lsg.get("A"):
            case "1":
                exterior = pts
            case "2":
                holes.append(pts)
    return Polygon(exterior, holes) if exterior else None


def load_isoxml(path: Path) -> Iterator[tuple[str, list[Polygon], CRS]]:
    root = etree.parse(str(path)).getroot()
    for pfd in root.iterfind(".//PFD"):  # PFD = Partfield
        field_id = pfd.get("A") or "unknown"  # A = PartfieldId
        parts = [poly for pln in pfd.iterfind("PLN")
                 if (poly := _polygon_from_pln(pln)) is not None]
        yield field_id, parts, WGS84  # ISO 11783-10 mandates WGS84


def load_shapefile(path: Path) -> Iterator[tuple[str, list[Polygon], CRS]]:
    prj = path.with_suffix(".prj")
    if not prj.exists():
        raise ValueError(f"{path.name}: missing .prj sidecar; CRS undefined")
    source = CRS.from_wkt(prj.read_text())
    reader = shapefile.Reader(str(path), encoding="cp437")  # dBASE III+ default codec
    for i, sr in enumerate(reader.iterShapeRecords()):
        geom = shape(sr.shape.__geo_interface__)
        parts = list(geom.geoms) if geom.geom_type == "MultiPolygon" else [geom]
        field_id = str(sr.record.as_dict().get("FARM_ID", f"row{i}"))
        yield field_id, parts, source


def parse_field_boundaries(path: Path) -> list[FieldBoundary]:
    match path.suffix.lower():
        case ".xml":
            rows = load_isoxml(path)
        case ".shp":
            rows = load_shapefile(path)
        case other:
            raise ValueError(f"unsupported boundary format: {other!r}")
    return [normalize_boundary(fid, parts, crs) for fid, parts, crs in rows]

The precision guard deliberately runs on the source coordinates before reprojection, because reprojecting compressed metres into WGS84 manufactures long decimal tails that hide the original loss. The source.to_epsg() is None check and the .prj existence check are the two guards that stop an undefined-CRS field from being committed as though it were WGS84.

Log patterns and observable signals

Every decode emits structured JSON so any stored boundary can be traced back to the exact archive and CRS it came from.

Success path:

{"event": "boundary_normalized", "field_id": "PFD1", "parts": 2, "source_epsg": 4326, "repaired": false}

Warning (multipart input needed topology repair before union — usually benign self-healing):

{"event": "boundary_normalized", "field_id": "PFD1", "parts": 1, "source_epsg": 32615, "repaired": true}

Warning (a .dbf attribute failed the UTF-8 probe and fell back to CP437):

{"event": "dbf_encoding_fallback", "file": "field_04.dbf", "field": "FARM_ID", "codec": "cp437", "action": "transcode_to_utf8"}

Error (source below the precision floor — the compression drift this guide prevents):

{"event": "boundary_quarantined", "field_id": "PFD1", "reason": "precision_below_floor", "min_dp": 6, "action": "manual_review"}

Error (shapefile without projection metadata):

{"event": "boundary_rejected", "file": "field_04.shp", "reason": "prj_missing", "action": "halt_ingest"}

When triaging, filter on boundary_rejected and boundary_quarantined first: they mean a field never entered the registry and any downstream job keyed on it will find nothing. A high repaired: true rate on ISOXML sources points upstream at a firmware revision emitting malformed multipart parts — correlate it against the machine and export version captured by Connecting John Deere API to Python backend.

Safe override protocol

Occasionally a legitimate boundary must be admitted despite a missing .prj or an implicit CRS — a historical export from a decommissioned controller, or a regional dataset whose UTM zone is known out of band. The override must never bypass normalize_boundary; it only supplies the missing CRS so the same validation still runs.

Guard conditions, all mandatory:

1. Explicit zone, never a guess. The operator supplies a concrete EPSG (for example UTM zone 15N, EPSG:32615); the parser never heuristically infers a zone and commits silently.
2. Centroid sanity gate. After reprojection the field centroid must fall within a known regional bounding box; if it does not, the override is rejected and the archive returns to quarantine.
3. Dry-run first. The override runs with persistence disabled and must pass every topology check before it is allowed to write.
4. Immutable audit trail. The original archive hash, the assumed EPSG, and the operator identity are written to an append-only ledger and flagged for compliance review, consistent with EPA pesticide application recordkeeping expectations.

from pathlib import Path

from pyproj import CRS


def override_missing_crs(
    path: Path,
    assumed_epsg: int,
    region_bbox: tuple[float, float, float, float],  # minx, miny, maxx, maxy in WGS84
    approver: str,
    dry_run: bool = True,
) -> list[FieldBoundary]:
    source = CRS.from_epsg(assumed_epsg)  # explicit, never inferred
    reader = shapefile.Reader(str(path), encoding="cp437")
    boundaries: list[FieldBoundary] = []
    for i, sr in enumerate(reader.iterShapeRecords()):
        geom = shape(sr.shape.__geo_interface__)
        parts = list(geom.geoms) if geom.geom_type == "MultiPolygon" else [geom]
        field_id = str(sr.record.as_dict().get("FARM_ID", f"row{i}"))
        fb = normalize_boundary(field_id, parts, source)
        cx, cy = fb.geometry.centroid.x, fb.geometry.centroid.y
        minx, miny, maxx, maxy = region_bbox
        if not (minx <= cx <= maxx and miny <= cy <= maxy):
            raise ValueError(f"{field_id}: centroid {cx:.5f},{cy:.5f} outside region")
        boundaries.append(fb)
    logger.warning(
        '{"event":"crs_override","file":"%s","assumed_epsg":%d,'
        '"approver":"%s","dry_run":%s}',
        path.name, assumed_epsg, approver, str(dry_run).lower(),
    )
    return boundaries  # caller persists only when dry_run is False and review passes

The assumed_epsg is a required argument with no default: forcing the caller to name the CRS is what keeps a convenient override from degenerating into the very undefined-CRS bug it exists to work around.

Troubleshooting

• boundary_rejected with reason: prj_missing. Root cause: a shapefile arrived without its .prj sidecar, so the CRS is undefined. Remediation: request a re-export that includes the projection file; never assume EPSG:4326. If the zone is genuinely known, route the archive through override_missing_crs with an explicit EPSG and the centroid gate.
• Field renders in the ocean or thousands of km off. Root cause: locally projected metres were interpreted as WGS84 degrees, or a transform ran with always_xy=False and swapped latitude and longitude. Remediation: confirm the source EPSG, and pin always_xy=True on every Transformer; the centroid sanity gate catches this before persistence.
• TopologyException: found non-noded intersection downstream. Root cause: multipart parts were unioned without OGC winding, so an interior ring was read as filled area and self-intersects. Remediation: keep orient(sign=1.0) before storage and rely on make_valid; verify ISOXML LSG type codes (1 exterior, 2 interior) are being honored in _polygon_from_pln.
• boundary_quarantined with reason: precision_below_floor. Root cause: an OEM compressed vertices to five or fewer decimals, below the 0.11 m drift threshold. Remediation: re-pull the archive at full precision from the source system; do not lower min_precision_dp to force ingestion, as that admits the drift the floor exists to block.
• Garbled FARM_ID or truncated parcel IDs. Root cause: the .dbf was read with the wrong codec, or a UTF-8 identifier exceeded the dBASE 10-character field limit. Remediation: probe for a BOM and only fall back to cp437 on failure, and map long identifiers to an internal UUID registry rather than truncating them — downstream GIS and EPA/USDA rule mapping reject non-conforming manifests.

Frequently asked questions

Should a shapefile without a .prj file ever be ingested? Not on the normal path — a missing .prj means the CRS is undefined, and assuming WGS84 is how a field lands in the ocean. Reject it, and only admit it through the explicit CRS-override protocol above.

Why enforce OGC ring winding before unioning ISOXML parts? Because LSG type codes distinguish exterior rings from interior holes; merging without correct winding reads a hole as filled area and self-intersects, surfacing as a TopologyException far downstream instead of at ingest.

What is the precision floor protecting against? OEM compression to five decimal places can move a vertex ~0.11 m, enough to cause false overlaps in buffer and prescription logic. The parser quarantines sub-floor sources instead of snapping and hiding the loss.

Related

• Connecting John Deere API to Python backend — the OAuth and pagination layer that delivers the archives this parser decodes.
• Enforcing EPA buffer zones with geospatial Python — the setback engine that consumes these normalized boundaries.
• Field Schema Design — the canonical field record every decoded boundary is written into.

Up: Equipment Telemetry Parsing · Farm Data Ingestion & Field Boundary Synchronization