Threshold Mapping for Crop Health

pip install rasterio numpy scikit-image geopandas scipy

import rasterio
import numpy as np
from pathlib import Path

def open_and_validate_vi(input_path: str) -> rasterio.DatasetReader:
    """Open a vegetation index raster and enforce minimum quality criteria."""
    p = Path(input_path)
    if not p.exists():
        raise FileNotFoundError(f"Raster not found: {input_path}")

    src = rasterio.open(input_path)

    # Must be a projected CRS — geographic CRS produces inaccurate zone areas
    if src.crs is None:
        raise ValueError("Raster has no CRS — reproject before classifying.")
    if src.crs.is_geographic:
        raise ValueError(
            f"CRS {src.crs.to_epsg()} is geographic. "
            "Reproject to a UTM zone (e.g., EPSG:32632) before classifying."
        )

    # Must be single-band float
    assert src.count == 1, f"Expected single-band raster, got {src.count} bands."
    assert src.dtypes[0] in ("float32", "float64"), (
        f"Expected float32/64 index raster, got {src.dtypes[0]}. "
        "Cast to float32 after radiometric conversion."
    )

    # Nodata must be declared so windowed reads handle border tiles correctly
    if src.nodata is None:
        src.close()
        raise ValueError("Raster has no nodata value set. Set nodata in rasterio profile.")

    print(f"Validated: {p.name} | CRS: {src.crs.to_epsg()} | "
          f"Size: {src.width}×{src.height} | Nodata: {src.nodata}")
    return src

from skimage.filters import threshold_multiotsu

def compute_thresholds(
    src: rasterio.DatasetReader,
    strategy: str = "fixed",
    fixed_breaks: list = None,
    n_classes: int = 4,
    sample_fraction: float = 0.05
) -> list[float]:
    """
    Derive threshold breakpoints for an n-class classification.

    For adaptive strategies, a random spatial sample avoids loading the full raster.
    sample_fraction controls what fraction of pixels to sample (default 5%).
    """
    if strategy == "fixed":
        breaks = fixed_breaks if fixed_breaks else [0.35, 0.55, 0.75]
        assert len(breaks) == n_classes - 1, (
            f"Need {n_classes-1} breakpoints for {n_classes} classes, "
            f"got {len(breaks)}."
        )
        return breaks

    if strategy in ("otsu", "jenks"):
        # Stratified random sample to avoid loading full raster into RAM
        total_pixels = src.width * src.height
        sample_size = max(10_000, int(total_pixels * sample_fraction))
        nodata = src.nodata

        rows = np.random.randint(0, src.height, sample_size)
        cols = np.random.randint(0, src.width, sample_size)

        # Read individual pixel windows (vectorised via rasterio sample)
        xy = list(zip(
            src.transform.c + cols * src.transform.a,
            src.transform.f + rows * src.transform.e
        ))
        values = np.array([v[0] for v in src.sample(xy)], dtype=np.float32)
        valid = values[~np.isnan(values) & (values != nodata) & np.isfinite(values)]

        if len(valid) < n_classes * 50:
            raise ValueError(
                f"Insufficient valid pixels for adaptive thresholding "
                f"(need ≥{n_classes*50}, got {len(valid)})."
            )

        if strategy == "otsu":
            thresholds = threshold_multiotsu(valid, classes=n_classes)
            return sorted(thresholds.tolist())

        if strategy == "jenks":
            # Jenks via Fisher-Jenks algorithm from mapclassify (optional dep)
            try:
                import mapclassify
                jnk = mapclassify.JenksCaspall(valid, k=n_classes)
                return sorted(jnk.bins[:-1].tolist())
            except ImportError:
                raise ImportError(
                    "mapclassify is required for Jenks thresholds: "
                    "pip install mapclassify"
                )

    raise ValueError(f"Unknown strategy '{strategy}'. Choose 'fixed', 'otsu', or 'jenks'.")

from rasterio.windows import Window

def classify_windowed(
    src: rasterio.DatasetReader,
    output_path: str,
    thresholds: list[float],
    window_size: int = 1024
) -> None:
    """
    Classify a vegetation index raster into discrete health zones using
    windowed I/O. Writes a uint8 GeoTIFF where 0 = nodata, 1–N = class.
    """
    nodata = src.nodata
    n_classes = len(thresholds) + 1
    breaks = sorted(thresholds)

    profile = src.profile.copy()
    profile.update(dtype="uint8", count=1, nodata=0, compress="lzw", predictor=2)

    classified_arr = np.zeros((src.height, src.width), dtype=np.uint8)

    for row_off in range(0, src.height, window_size):
        for col_off in range(0, src.width, window_size):
            h = min(window_size, src.height - row_off)
            w = min(window_size, src.width - col_off)
            win = Window(col_off, row_off, w, h)

            data = src.read(1, window=win).astype(np.float32)
            valid = (data != nodata) & np.isfinite(data)

            if not np.any(valid):
                continue

            classes = np.zeros_like(data, dtype=np.uint8)

            # Assign classes 1..N based on sorted breakpoints
            prev_upper = -np.inf
            for cls_id, upper in enumerate(breaks, start=1):
                mask = valid & (data > prev_upper) & (data <= upper)
                classes[mask] = cls_id
                prev_upper = upper
            # Last class captures everything above the highest breakpoint
            classes[valid & (data > breaks[-1])] = n_classes

            classified_arr[row_off:row_off+h, col_off:col_off+w] = classes

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(classified_arr, 1)

    print(f"Classified raster written to {output_path}")
    # Verify class distribution
    unique, counts = np.unique(classified_arr[classified_arr > 0], return_counts=True)
    for u, c in zip(unique, counts):
        print(f"  Class {u}: {c:,} px ({100*c/classified_arr.size:.1f}%)")

from scipy.ndimage import label, generic_filter

def sieve_and_smooth(
    classified_path: str,
    min_patch_px: int = 25,
    smooth_radius: int = 1
) -> np.ndarray:
    """
    Remove small isolated patches (sieve) and apply majority smoothing.

    min_patch_px: minimum connected patch size in pixels (e.g. 25 px ≈ 2.5 m²
                  at 10 cm GSD — roughly the size of one soybean row).
    smooth_radius: neighbourhood radius for majority filter (pixels).
    """
    with rasterio.open(classified_path) as src:
        arr = src.read(1).copy()
        profile = src.profile.copy()

    nodata_mask = arr == 0

    # Sieve: label connected components per class, remove undersized patches
    sieved = arr.copy()
    for cls in np.unique(arr[arr > 0]):
        binary = (arr == cls)
        labelled, n_feat = label(binary)
        for feat_id in range(1, n_feat + 1):
            patch = labelled == feat_id
            if patch.sum() < min_patch_px:
                sieved[patch] = 0  # mark for reassignment

    # Fill sieved-out pixels with nearest valid neighbour (majority in 3×3)
    def majority(values):
        vals = values[values > 0]
        if len(vals) == 0:
            return 0
        counts = np.bincount(vals.astype(int))
        return int(np.argmax(counts))

    kernel_size = 2 * smooth_radius + 1
    smoothed = generic_filter(
        sieved.astype(np.float32),
        majority,
        size=kernel_size,
        mode="nearest"
    ).astype(np.uint8)

    smoothed[nodata_mask] = 0  # restore nodata pixels

    with rasterio.open(classified_path, "r+") as dst:
        dst.write(smoothed, 1)

    return smoothed

import geopandas as gpd
from rasterio.features import shapes
from shapely.geometry import shape

def vectorise_zones(
    classified_path: str,
    output_gpkg: str,
    class_labels: dict = None
) -> gpd.GeoDataFrame:
    """
    Convert a classified zone raster to a GeoPackage of clean zone polygons.
    Requires that the input CRS is projected (metres) for accurate area fields.
    """
    labels = class_labels or {1: "Stress", 2: "Moderate", 3: "Healthy", 4: "High Vigour"}

    with rasterio.open(classified_path) as src:
        assert not src.crs.is_geographic, (
            "Raster CRS must be projected for area calculations. "
            f"Got: {src.crs}"
        )
        arr = src.read(1)
        transform = src.transform
        crs = src.crs

    geoms, class_vals, areas = [], [], []
    for geom_dict, val in shapes(arr, mask=(arr > 0), transform=transform):
        geom = shape(geom_dict)
        geoms.append(geom)
        class_vals.append(int(val))
        areas.append(geom.area)  # m² when CRS is projected

    gdf = gpd.GeoDataFrame(
        {
            "class_id": class_vals,
            "label": [labels.get(v, f"Class {v}") for v in class_vals],
            "area_ha": [a / 10_000 for a in areas],
        },
        geometry=geoms,
        crs=crs,
    )

    # Remove slivers (< 1 m²)
    gdf = gdf[gdf["area_ha"] > 0.0001].reset_index(drop=True)

    gdf.to_file(output_gpkg, layer="crop_health_zones", driver="GPKG")
    print(f"Vectorised zones saved: {output_gpkg}")
    print(gdf.groupby("label")["area_ha"].sum().round(2))
    return gdf

def run_threshold_pipeline(
    vi_raster: str,
    classified_out: str,
    vector_out: str,
    strategy: str = "fixed",
    fixed_breaks: list = None,
    n_classes: int = 4,
    min_patch_px: int = 25,
) -> gpd.GeoDataFrame:
    """
    End-to-end crop health threshold mapping pipeline.
    Returns a GeoDataFrame of management zone polygons.
    """
    src = open_and_validate_vi(vi_raster)
    try:
        breaks = compute_thresholds(src, strategy=strategy,
                                    fixed_breaks=fixed_breaks,
                                    n_classes=n_classes)
        print(f"Thresholds ({strategy}): {breaks}")
        classify_windowed(src, classified_out, breaks)
    finally:
        src.close()

    sieve_and_smooth(classified_out, min_patch_px=min_patch_px)
    return vectorise_zones(classified_out, vector_out)


# Example call — NDVI maize mid-season with fixed agronomic thresholds
if __name__ == "__main__":
    zones = run_threshold_pipeline(
        vi_raster="ndvi_field_utm32.tif",
        classified_out="crop_health_classified.tif",
        vector_out="management_zones.gpkg",
        strategy="fixed",
        fixed_breaks=[0.35, 0.55, 0.75],
    )

# After classify_windowed() — sanity check class distribution
with rasterio.open(classified_out) as src:
    arr = src.read(1)
    valid = arr[arr > 0]
    if len(valid) == 0:
        raise ValueError("No valid classified pixels — check nodata handling.")
    for cls in range(1, n_classes + 1):
        pct = 100 * (arr == cls).sum() / len(valid)
        print(f"Class {cls}: {pct:.1f}%")
        # Warn if any single class captures >90% of the field
        assert pct < 90, f"Class {cls} dominates ({pct:.1f}%) — verify thresholds."