Shapefile Validation for Farm Equipment

import geopandas as gpd
import pandas as pd
from pathlib import Path


def inspect_shapefile(input_path: str) -> gpd.GeoDataFrame:
    gdf = gpd.read_file(input_path, engine="pyogrio")

    # Enforce single geometry type — mixed types break most firmware parsers
    geom_types = gdf.geom_type.unique()
    assert len(geom_types) == 1, f"Mixed geometry types: {geom_types}. Separate before validation."

    # Mandatory agronomic columns
    required_cols = {"zone_id", "rate_kg_ha", "crop_code"}
    missing = required_cols - set(gdf.columns)
    assert not missing, f"Missing required columns: {missing}"

    # Detect field name collisions after 10-char truncation
    truncated = [c[:10] for c in gdf.columns if c != "geometry"]
    collisions = [c for c in truncated if truncated.count(c) > 1]
    assert not collisions, (
        f"Column names collide after 10-char truncation: {set(collisions)}. "
        "Rename before export."
    )

    print(f"Schema OK — {len(gdf)} features, CRS: {gdf.crs}")
    return gdf

from pyproj import CRS


def harmonize_crs(gdf: gpd.GeoDataFrame, target_epsg: int) -> gpd.GeoDataFrame:
    target = CRS.from_epsg(target_epsg)

    if gdf.crs is None:
        raise ValueError(
            "Shapefile has no CRS definition. Assign the source CRS manually "
            "with gdf.set_crs(epsg=...) before calling this function."
        )

    if not gdf.crs.equals(target):
        source_epsg = gdf.crs.to_epsg()
        print(f"Reprojecting EPSG:{source_epsg} → EPSG:{target_epsg}")
        gdf = gdf.to_crs(epsg=target_epsg)

    assert gdf.crs.to_epsg() == target_epsg, "CRS after reprojection does not match target."
    return gdf

def validate_and_repair_geometry(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    invalid_mask = ~gdf.is_valid
    invalid_count = int(invalid_mask.sum())

    if invalid_count > 0:
        print(f"Repairing {invalid_count} invalid geometries with make_valid()...")
        gdf.loc[invalid_mask, "geometry"] = (
            gdf.loc[invalid_mask, "geometry"].make_valid()
        )

    # Remove sliver polygons under 10 m² that cause erratic nozzle switching
    pre_sliver_count = len(gdf)
    gdf = gdf[gdf.geometry.area >= 10.0].copy()
    removed = pre_sliver_count - len(gdf)
    if removed:
        print(f"Removed {removed} sliver features (area < 10 m²).")

    # Explode MultiPolygons to simple Polygons for maximum controller compatibility
    gdf = gdf.explode(index_parts=False).reset_index(drop=True)

    # Final validity assertion — fail loudly if repair did not succeed
    still_invalid = (~gdf.is_valid).sum()
    assert still_invalid == 0, (
        f"{still_invalid} geometries remain invalid after repair. "
        "See /yield-mapping-variable-rate-prescription-generation/"
        "shapefile-validation-for-farm-equipment/"
        "debugging-shapefile-geometry-errors-in-qgis-and-python/ "
        "for manual intervention steps."
    )

    print(f"Geometry validation passed — {len(gdf)} features.")
    return gdf

def enforce_attribute_constraints(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    # Coerce rate column to float; replace non-numeric strings with NaN then fill
    gdf["rate_kg_ha"] = pd.to_numeric(gdf["rate_kg_ha"], errors="coerce")
    nan_rate_count = int(gdf["rate_kg_ha"].isna().sum())
    if nan_rate_count:
        print(f"Warning: {nan_rate_count} NaN rate values replaced with 0.0")
    gdf["rate_kg_ha"] = gdf["rate_kg_ha"].fillna(0.0)

    # Enforce agronomic application bounds (0–500 kg/ha for granular products)
    gdf.loc[gdf["rate_kg_ha"] < 0, "rate_kg_ha"] = 0.0
    gdf.loc[gdf["rate_kg_ha"] > 500, "rate_kg_ha"] = 500.0

    # Normalise categorical identifiers
    gdf["zone_id"] = gdf["zone_id"].astype(str).str.strip()
    gdf["crop_code"] = gdf["crop_code"].astype(str).str.upper().str.strip()

    # Drop features missing critical identifiers — do not silently pass them through
    before = len(gdf)
    gdf = gdf.dropna(subset=["zone_id", "crop_code"])
    dropped = before - len(gdf)
    if dropped:
        print(f"Dropped {dropped} features with missing zone_id or crop_code.")

    # Post-condition: no NaN in rate or identifier columns
    assert gdf["rate_kg_ha"].notna().all()
    assert gdf["zone_id"].notna().all()

    return gdf

def export_validated_shapefile(gdf: gpd.GeoDataFrame, output_path: str) -> None:
    # Build explicit 10-character column mapping — fail if any collision detected
    col_mapping: dict[str, str] = {}
    seen_truncated: set[str] = set()
    for col in gdf.columns:
        if col == "geometry":
            continue
        short = col[:10]
        assert short not in seen_truncated, (
            f"Column name collision after truncation: '{col}' → '{short}'. "
            "Rename the column before exporting."
        )
        col_mapping[col] = short
        seen_truncated.add(short)

    gdf = gdf.rename(columns=col_mapping)

    # Write with pyogrio for structured error output on schema violations
    gdf.to_file(output_path, driver="ESRI Shapefile", encoding="UTF-8", engine="pyogrio")

    # Read back and verify feature count integrity
    written = gpd.read_file(output_path, engine="pyogrio")
    assert len(written) == len(gdf), (
        f"Export integrity failure: wrote {len(gdf)} features, "
        f"read back {len(written)}."
    )
    print(f"Validated shapefile written: {output_path} ({len(written)} features)")

def validate_prescription_shapefile(
    input_path: str,
    output_path: str,
    target_epsg: int,
) -> None:
    gdf = inspect_shapefile(input_path)
    gdf = harmonize_crs(gdf, target_epsg=target_epsg)
    gdf = validate_and_repair_geometry(gdf)
    gdf = enforce_attribute_constraints(gdf)
    export_validated_shapefile(gdf, output_path)


# Example: validate a wheat VRA prescription for UTM Zone 14N
validate_prescription_shapefile(
    input_path="raw/wheat_vra_2024.shp",
    output_path="validated/wheat_vra_2024_validated.shp",
    target_epsg=32614,
)

from shapely.ops import unary_union

def check_overlaps(gdf: gpd.GeoDataFrame, tolerance_m2: float = 1.0) -> None:
    merged = unary_union(gdf.geometry)
    total_individual = gdf.geometry.area.sum()
    overlap_area = total_individual - merged.area
    if overlap_area > tolerance_m2:
        raise ValueError(
            f"Zones overlap by {overlap_area:.1f} m². "
            "Resolve topology before exporting to equipment."
        )

gdf = gpd.read_file(input_path, engine="pyogrio", encoding="latin-1")

assert (gdf["zone_id"].str.len() > 0).all(), \
    "zone_id contains empty strings after normalisation."

def verify_output(output_path: str, expected_epsg: int) -> None:
    gdf = gpd.read_file(output_path, engine="pyogrio")

    # CRS integrity
    assert gdf.crs.to_epsg() == expected_epsg, \
        f"Output CRS mismatch: expected EPSG:{expected_epsg}, got {gdf.crs.to_epsg()}"

    # All geometries valid
    assert gdf.is_valid.all(), "Output contains invalid geometries."

    # No NaN rates
    assert gdf["rate_kg_ha"].notna().all(), "NaN values remain in rate_kg_ha."

    # Field name length
    long_names = [c for c in gdf.columns if c != "geometry" and len(c) > 10]
    assert not long_names, f"Field names exceed 10 chars: {long_names}"

    # Bounding box sanity — reject empty or near-zero extent files
    bounds = gdf.total_bounds  # (minx, miny, maxx, maxy)
    extent_m = max(bounds[2] - bounds[0], bounds[3] - bounds[1])
    assert extent_m > 10, f"Suspiciously small spatial extent: {extent_m:.1f} m"

    # Rate statistics — flag if max rate is zero (indicates attribute coercion failure)
    max_rate = gdf["rate_kg_ha"].max()
    assert max_rate > 0, "Maximum application rate is 0 — check attribute pipeline."

    print(f"Verification passed: {len(gdf)} features, "
          f"EPSG:{expected_epsg}, rate range "
          f"{gdf['rate_kg_ha'].min():.1f}–{gdf['rate_kg_ha'].max():.1f} kg/ha")

import logging

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("shapefile_validation")


def run_validation_pipeline(
    input_path: str,
    output_dir: str,
    target_epsg: int,
) -> str:
    from pathlib import Path
    import datetime

    ts = datetime.datetime.utcnow().strftime("%Y%m%dT%H%M%SZ")
    stem = Path(input_path).stem
    output_path = str(Path(output_dir) / f"{stem}_validated_{ts}.shp")

    logger.info("Starting shapefile validation: %s", input_path)
    validate_prescription_shapefile(input_path, output_path, target_epsg)
    verify_output(output_path, expected_epsg=target_epsg)
    logger.info("Validation complete: %s", output_path)
    return output_path