Variable Rate Export to ISOXML: Python Workflow for Precision Agriculture

import geopandas as gpd
from shapely.validation import make_valid
from pyproj import CRS
import logging

logger = logging.getLogger(__name__)

def normalize_prescription(input_path: str, metric_epsg: int = 32632) -> gpd.GeoDataFrame:
    """
    Load, repair, and validate a prescription GeoDataFrame.
    Repairs topology in the supplied metric CRS, then reprojects to EPSG:4326.

    Args:
        input_path: Path to a shapefile or GeoPackage.
        metric_epsg: EPSG code for the local metric CRS used during repair.
                     Default 32632 = UTM zone 32N (central Europe).
                     Adjust to the field's UTM zone before use.
    Returns:
        GeoDataFrame in EPSG:4326, topology-clean, coordinate-rounded.
    Raises:
        ValueError: If the dataset is empty or coordinates fall outside WGS84 bounds.
    """
    gdf = gpd.read_file(input_path)
    if gdf.empty:
        raise ValueError("Input dataset contains no features.")

    # Repair in metric CRS so buffer distances are meaningful
    if gdf.crs is None:
        raise ValueError("Input has no CRS. Assign EPSG before calling this function.")
    metric_crs = CRS.from_epsg(metric_epsg)
    if gdf.crs != metric_crs:
        logger.info("Reprojecting to metric CRS %s for topology repair", metric_epsg)
        gdf = gdf.to_crs(metric_crs)

    gdf["geometry"] = gdf["geometry"].apply(make_valid)
    gdf = gdf[~gdf["geometry"].is_empty].reset_index(drop=True)

    # Reproject to WGS84 for ISOXML output
    wgs84 = CRS.from_epsg(4326)
    gdf = gdf.to_crs(wgs84)

    # Validate geographic bounds
    min_lon, min_lat, max_lon, max_lat = gdf.total_bounds
    assert -180 <= min_lon <= max_lon <= 180, "Longitude out of WGS84 range"
    assert -90 <= min_lat <= max_lat <= 90, "Latitude out of WGS84 range"

    # Round to 7 d.p. (~1 cm) to avoid firmware truncation artifacts
    from shapely import set_precision
    gdf["geometry"] = gdf["geometry"].apply(lambda g: set_precision(g, grid_size=1e-7))

    logger.info("Normalized %d features, CRS: EPSG:4326", len(gdf))
    return gdf

gdf = normalize_prescription("prescriptions/corn_n_zones.gpkg", metric_epsg=32615)
assert gdf.crs.to_epsg() == 4326, "CRS must be EPSG:4326 after normalization"
assert gdf.geometry.is_valid.all(), "Topology repair failed on at least one feature"
print(f"Features ready for export: {len(gdf)}, CRS: {gdf.crs}")

from lxml import etree
import uuid

NS = "urn:iso:std:iso:11783:-10:taskdata:1.0"

def build_pdt_element(product_id: str, name: str, unit: str = "mg1ha-1") -> etree._Element:
    """
    Build a single PDT (Product) XML element.

    ISOXML unit strings follow the DDI registry.  Common values:
      mg1ha-1  →  kg/ha (mass per area, scaled ×1e6 internally)
      ml1ha-1  →  L/ha  (volume per area)
    Always confirm the DDI string against the terminal's implementation guide.
    """
    pdt = etree.Element("{%s}PDT" % NS, attrib={
        "A": product_id,           # Unique product reference ID
        "B": name,                 # Human-readable name (max 32 chars)
        "C": unit,                 # DDI unit string
        "D": str(uuid.uuid4())[:8] # Optional manufacturer-assigned UID
    })
    return pdt


def map_zones_to_rates(gdf: gpd.GeoDataFrame, rate_col: str, product_id: str) -> list[dict]:
    """
    Extract zone geometry and rate metadata for TZN generation.

    Args:
        gdf: GeoDataFrame in EPSG:4326 with a numeric rate column.
        rate_col: Column name containing application rates in the PDT's unit.
        product_id: Must match the 'A' attribute of the corresponding PDT element.
    Returns:
        List of dicts with keys: zone_id, geometry, rate, product_id.
    Raises:
        KeyError: If rate_col is absent.
        ValueError: If any rate is negative or non-finite.
    """
    if rate_col not in gdf.columns:
        raise KeyError(f"Column '{rate_col}' not found. Available: {list(gdf.columns)}")

    import numpy as np
    rates = gdf[rate_col].to_numpy(dtype=float)
    if not np.all(np.isfinite(rates)):
        raise ValueError("Non-finite values detected in rate column. Fill or drop before export.")
    if np.any(rates < 0):
        raise ValueError("Negative application rates are invalid. Clamp to 0 before export.")

    return [
        {
            "zone_id": f"Z{i+1:04d}",
            "geometry": row.geometry,
            "rate": float(row[rate_col]),
            "product_id": product_id,
        }
        for i, row in gdf.iterrows()
    ]

def build_taskdata_xml(zones: list[dict], pdt_element: etree._Element,
                       task_name: str = "VRA_Prescription") -> bytes:
    """
    Assemble a complete, namespace-correct TASKDATA.XML document.

    Returns raw UTF-8 bytes (no BOM) ready for file write.
    Element order: TASKDATA → PDT → TSK → TZN* → CTR* → LSG → LSG/PNT
    """
    nsmap = {None: NS}
    root = etree.Element("{%s}TASKDATA" % NS, nsmap=nsmap, attrib={
        "VersionMajor": "4",
        "VersionMinor": "3",
        "ManagementSoftwareManufacturer": "ag-gis-exporter",
        "ManagementSoftwareVersion": "1.0",
        "DataTransferOrigin": "1",  # 1 = FMIS (farm management information system)
    })

    # Attach product definition first — TSK references it by A-attribute
    root.append(pdt_element)

    # TSK (Task) element
    tsk = etree.SubElement(root, "{%s}TSK" % NS, attrib={
        "A": "TSK1",
        "B": task_name,
        "G": "1",  # TaskStatus: 1 = planned
    })

    for zone in zones:
        # TZN = Treatment Zone — one per management zone polygon
        tzn = etree.SubElement(tsk, "{%s}TZN" % NS, attrib={
            "A": zone["zone_id"],
        })

        # ProcessDataVariable: links zone to product and rate
        etree.SubElement(tzn, "{%s}PDV" % NS, attrib={
            "A": "0007",               # DDI 0x0007 = Setpoint Application Rate
            "B": str(int(round(zone["rate"] * 1e6))),  # Value in mg1ha-1 units ×1e6
            "C": zone["product_id"],   # References PDT.A
        })

        # CTR (Polygon boundary) — serialize exterior ring
        ctr = etree.SubElement(tzn, "{%s}CTR" % NS)
        coords = list(zone["geometry"].exterior.coords)
        for lon, lat in coords:
            etree.SubElement(ctr, "{%s}PNT" % NS, attrib={
                "A": "2",              # Point type: 2 = field access/boundary
                "C": f"{lat:.7f}",     # Note: ISOXML uses lat (C) before lon (D)
                "D": f"{lon:.7f}",
            })

    return etree.tostring(root, pretty_print=True,
                          xml_declaration=True, encoding="UTF-8")

from lxml import etree
from pathlib import Path

def validate_isoxml(xml_bytes: bytes, xsd_path: str) -> None:
    """
    Validate serialised ISOXML bytes against the ISO 11783-10 XSD.

    Args:
        xml_bytes: Output of build_taskdata_xml().
        xsd_path: Path to the XSD file obtained from the terminal manufacturer
                  or the AEF ISOBUS implementation guide.
    Raises:
        etree.XMLSchemaError: With a descriptive message listing all violations.
    """
    xsd_doc = etree.parse(xsd_path)
    schema = etree.XMLSchema(xsd_doc)
    doc = etree.fromstring(xml_bytes)
    if not schema.validate(doc):
        errors = "\n".join(str(e) for e in schema.error_log)
        raise etree.XMLSchemaError(
            f"ISOXML validation failed with {len(schema.error_log)} error(s):\n{errors}"
        )
    print("XSD validation passed.")

doc = etree.fromstring(xml_bytes)
# Assert required elements are present
assert doc.find(".//{%s}PDT" % NS) is not None, "Missing PDT element"
assert doc.find(".//{%s}TSK" % NS) is not None, "Missing TSK element"
tzn_count = len(doc.findall(".//{%s}TZN" % NS))
assert tzn_count == len(zones), f"TZN count {tzn_count} ≠ zone count {len(zones)}"
print(f"Structural check passed: {tzn_count} treatment zones serialised.")

import zipfile
import os
import tempfile
from pathlib import Path

def package_isoxml(xml_bytes: bytes, output_zip: str,
                   extra_files: dict[str, bytes] | None = None) -> Path:
    """
    Package TASKDATA.XML (and any auxiliary files) into a terminal-ready ZIP.

    Args:
        xml_bytes: Validated XML as bytes.
        output_zip: Output path for the ZIP archive.
        extra_files: Optional dict mapping archive filename → bytes content
                     (e.g., {"TASKDATA.XSD": xsd_bytes} for self-validating packages).
    Returns:
        Path to the created ZIP file.
    """
    out_path = Path(output_zip)
    out_path.parent.mkdir(parents=True, exist_ok=True)

    with zipfile.ZipFile(out_path, "w", zipfile.ZIP_DEFLATED) as zf:
        # Write TASKDATA.XML — no BOM, LF line endings enforced by lxml
        zf.writestr("TASKDATA.XML", xml_bytes)
        if extra_files:
            for name, content in extra_files.items():
                zf.writestr(name, content)

    zip_size_kb = out_path.stat().st_size / 1024
    print(f"Package created: {out_path} ({zip_size_kb:.1f} KB)")
    return out_path

from pathlib import Path

# ---- Inputs ----
INPUT_SHAPEFILE = "prescriptions/corn_n_zones.gpkg"
RATE_COLUMN     = "n_kg_ha"
PRODUCT_NAME    = "Nitrogen"
OUTPUT_ZIP      = "exports/corn_n_vra.zip"
XSD_PATH        = "schemas/TASKDATA.xsd"   # Obtain from AEF or terminal OEM

# ---- Pipeline ----
# 1. Normalize
gdf = normalize_prescription(INPUT_SHAPEFILE, metric_epsg=32615)

# 2. Build product definition
pdt = build_pdt_element("PDT1", PRODUCT_NAME, unit="mg1ha-1")

# 3. Map zones
zones = map_zones_to_rates(gdf, RATE_COLUMN, product_id="PDT1")

# 4. Build XML
xml_bytes = build_taskdata_xml(zones, pdt, task_name="Corn_N_VRA_2024")

# 5. Validate
validate_isoxml(xml_bytes, XSD_PATH)

# 6. Package
package_isoxml(xml_bytes, OUTPUT_ZIP)
print("Export complete.")

import zipfile
from lxml import etree

def verify_package(zip_path: str, expected_zone_count: int) -> None:
    """Quick structural verification of a packaged ISOXML archive."""
    with zipfile.ZipFile(zip_path) as zf:
        names = zf.namelist()
        assert "TASKDATA.XML" in names, "TASKDATA.XML missing from archive"
        xml_bytes = zf.read("TASKDATA.XML")

    # Encoding: must be UTF-8, no BOM
    assert not xml_bytes.startswith(b'\xef\xbb\xbf'), \
        "UTF-8 BOM detected — will cause parser failure on most terminals"

    doc = etree.fromstring(xml_bytes)
    tzn_nodes = doc.findall(".//{%s}TZN" % NS)
    assert len(tzn_nodes) == expected_zone_count, \
        f"Expected {expected_zone_count} TZN elements, found {len(tzn_nodes)}"

    # Spot-check coordinate range on first zone
    first_lat = float(tzn_nodes[0].find(".//{%s}PNT" % NS).get("C"))
    first_lon = float(tzn_nodes[0].find(".//{%s}PNT" % NS).get("D"))
    assert -90 <= first_lat <= 90, f"First latitude out of range: {first_lat}"
    assert -180 <= first_lon <= 180, f"First longitude out of range: {first_lon}"

    print(f"Package verified: {len(tzn_nodes)} zones, first point lat={first_lat:.5f} lon={first_lon:.5f}")

verify_package(OUTPUT_ZIP, expected_zone_count=len(gdf))