from __future__ import annotations from collections import deque from dataclasses import dataclass from typing import Any import numpy as np INDEX_NAMES = ("ndvi", "ndre", "gndvi") def _normalised_difference(high: np.ndarray, low: np.ndarray) -> np.ndarray: numerator = high - low denominator = high + low return np.divide(numerator, denominator, out=np.full_like(numerator, np.nan), where=np.abs(denominator) > 1e-8) def compute_indices(green: np.ndarray, red: np.ndarray, red_edge: np.ndarray, nir: np.ndarray) -> dict[str, np.ndarray]: arrays = [np.asarray(value, dtype=np.float32) for value in (green, red, red_edge, nir)] if len({array.shape for array in arrays}) != 1 or arrays[0].ndim != 2: raise ValueError("All four reflectance bands must be aligned two-dimensional arrays") green_value, red_value, red_edge_value, nir_value = arrays return { "ndvi": _normalised_difference(nir_value, red_value), "ndre": _normalised_difference(nir_value, red_edge_value), "gndvi": _normalised_difference(nir_value, green_value), } def quality_check(bands: dict[str, np.ndarray], metadata: dict[str, Any]) -> dict[str, Any]: required = {"green", "red", "red_edge", "nir"} failures: list[str] = [] if set(bands) != required: failures.append("required_band_set") shapes = {np.asarray(value).shape for value in bands.values()} if len(shapes) != 1 or not shapes or len(next(iter(shapes))) != 2: failures.append("band_alignment") finite_fraction = min((float(np.isfinite(value).mean()) for value in bands.values()), default=0.0) out_of_range_fraction = max((float(((value < 0) | (value > 1)).mean()) for value in bands.values()), default=1.0) saturation_fraction = max((float((value >= 0.98).mean()) for value in bands.values()), default=1.0) if finite_fraction < 0.995: failures.append("finite_fraction") if out_of_range_fraction > 0.001: failures.append("reflectance_range") if saturation_fraction > 0.01: failures.append("saturation_fraction") for key in ( "dataProvenance", "captureId", "scaleFactor", "pixelSizeM", "orthomosaicCoRegistered", "bandMapValidated", "processingLevel", "measurementScale", "calibrationMethod", "crs", "gridTransform", ): if key not in metadata: failures.append(f"metadata:{key}") if metadata.get("orthomosaicCoRegistered") is not True: failures.append("orthomosaic_co_registration") if metadata.get("bandMapValidated") is not True: failures.append("band_map_validation") if metadata.get("processingLevel") != "orthomosaic_reflectance": failures.append("processing_level_not_reflectance") if metadata.get("measurementScale") not in {"relative_reflectance", "surface_reflectance"}: failures.append("measurement_scale") if metadata.get("dataProvenance") not in {"synthetic", "field"}: failures.append("data_provenance") return { "passed": not failures, "failures": sorted(set(failures)), "metrics": { "finiteFraction": round(finite_fraction, 6), "outOfRangeFraction": round(out_of_range_fraction, 6), "saturationFraction": round(saturation_fraction, 6), "shape": list(next(iter(shapes))) if len(shapes) == 1 else None, }, "interpretation": "QC pass permits index calculation, not agronomic diagnosis.", } def _robust_low_score(array: np.ndarray, valid: np.ndarray) -> np.ndarray: values = array[valid] median = float(np.nanmedian(values)) mad = float(np.nanmedian(np.abs(values - median))) scale = max(1.4826 * mad, 1e-4) return np.clip((median - array) / scale, 0.0, 6.0) def anomaly_score(indices: dict[str, np.ndarray], field_mask: np.ndarray | None = None) -> np.ndarray: if set(indices) != set(INDEX_NAMES): raise ValueError(f"indices must be exactly {INDEX_NAMES}") shape = indices["ndvi"].shape if any(value.shape != shape for value in indices.values()): raise ValueError("Index arrays are not aligned") valid = np.ones(shape, dtype=bool) if field_mask is None else np.asarray(field_mask, dtype=bool) valid &= np.logical_and.reduce([np.isfinite(value) for value in indices.values()]) if valid.sum() < 100: raise ValueError("At least 100 valid field pixels are required") stacked = np.stack([_robust_low_score(indices[name], valid) for name in INDEX_NAMES]) score = np.nanmean(stacked, axis=0) / 6.0 score[~valid] = np.nan return score.astype(np.float32) @dataclass(frozen=True) class Zone: zone_id: str pixel_count: int area_m2: float mean_score: float peak_score: float centroid_row: float centroid_col: float def find_priority_zones( score: np.ndarray, pixel_size_m: float, threshold: float = 0.38, minimum_area_m2: float = 2.0, ) -> list[Zone]: if pixel_size_m <= 0 or not 0 < threshold < 1 or minimum_area_m2 <= 0: raise ValueError("Invalid zone extraction parameters") active = np.isfinite(score) & (score >= threshold) visited = np.zeros(active.shape, dtype=bool) zones: list[Zone] = [] for row, col in zip(*np.where(active)): if visited[row, col]: continue queue: deque[tuple[int, int]] = deque([(int(row), int(col))]) visited[row, col] = True pixels: list[tuple[int, int]] = [] while queue: current_row, current_col = queue.popleft() pixels.append((current_row, current_col)) for next_row, next_col in ((current_row - 1, current_col), (current_row + 1, current_col), (current_row, current_col - 1), (current_row, current_col + 1)): if 0 <= next_row < active.shape[0] and 0 <= next_col < active.shape[1] and active[next_row, next_col] and not visited[next_row, next_col]: visited[next_row, next_col] = True queue.append((next_row, next_col)) area_m2 = len(pixels) * pixel_size_m * pixel_size_m if area_m2 < minimum_area_m2: continue values = np.array([score[r, c] for r, c in pixels]) zones.append(Zone( zone_id="", pixel_count=len(pixels), area_m2=round(area_m2, 3), mean_score=round(float(values.mean()), 4), peak_score=round(float(values.max()), 4), centroid_row=round(float(np.mean([p[0] for p in pixels])), 2), centroid_col=round(float(np.mean([p[1] for p in pixels])), 2), )) zones.sort(key=lambda zone: (-zone.mean_score, -zone.area_m2)) return [Zone(zone_id=f"P{index:02d}", **{key: value for key, value in zone.__dict__.items() if key != "zone_id"}) for index, zone in enumerate(zones, 1)]