Multispectral Sensors: Beyond Visible Light

class MultispectralAnalyzer:
    """Compute vegetation indices from multispectral bands."""

    def __init__(self, red, nir, red_edge=None, green=None):
        self.red = red.astype(float)
        self.nir = nir.astype(float)
        self.red_edge = red_edge.astype(float) if red_edge is not None else None
        self.green = green.astype(float) if green is not None else None

    def ndvi(self):
        """Normalized Difference Vegetation Index."""
        denominator = self.nir + self.red
        denominator[denominator == 0] = 1e-10  # Avoid division by zero
        return (self.nir - self.red) / denominator

    def ndre(self):
        """Normalized Difference Red Edge."""
        if self.red_edge is None:
            raise ValueError("Red edge band required for NDRE")
        denominator = self.nir + self.red_edge
        denominator[denominator == 0] = 1e-10
        return (self.nir - self.red_edge) / denominator

    def savi(self, L=0.5):
        """Soil-Adjusted Vegetation Index."""
        denominator = self.nir + self.red + L
        denominator[denominator == 0] = 1e-10
        return ((self.nir - self.red) * (1 + L)) / denominator

    def evi(self):
        """Enhanced Vegetation Index (2.5)."""
        denominator = self.nir + 6 * self.red - 7.5 * (self.green or self.red) + 1
        denominator[denominator == 0] = 1e-10
        return 2.5 * (self.nir - self.red) / denominator

    def classify_health(self, ndvi_values):
        """Classify vegetation health from NDVI."""
        classification = np.zeros_like(ndvi_values, dtype=int)
        classification[ndvi_values < 0.1] = 0   # Bare soil
        classification[(ndvi_values >= 0.1) & (ndvi_values < 0.3)] = 1  # Sparse
        classification[(ndvi_values >= 0.3) & (ndvi_values < 0.5)] = 2  # Moderate
        classification[(ndvi_values >= 0.5) & (ndvi_values < 0.7)] = 3  # Healthy
        classification[ndvi_values >= 0.7] = 4  # Dense/vigorous
        return classification

# Demo: Simulated multispectral data
np.random.seed(42)
height, width = 100, 100

# Create synthetic vegetation bands
red = np.random.uniform(0.05, 0.15, (height, width))
nir = np.random.uniform(0.3, 0.6, (height, width))

# Simulate a stressed area (lower NIR, higher red)
red[40:60, 40:60] = 0.25
nir[40:60, 40:60] = 0.15

analyzer = MultispectralAnalyzer(red=red, nir=nir)
ndvi = analyzer.ndvi()
health = analyzer.classify_health(ndvi)

print(f"NDVI range: {ndvi.min():.3f} to {ndvi.max():.3f}")
print(f"Mean NDVI: {ndvi.mean():.3f}")
print(f"Health classes: {np.bincount(health.flatten().astype(int))}")
print(f"Stressed pixels: {np.sum(health < 2)}")

from collections import defaultdict

class CropHealthMapper:
    """Generate crop health maps from drone multispectral data."""

    def __init__(self, field_boundary=None):
        self.field_boundary = field_boundary
        self.zone_stats = defaultdict(list)

    def create_health_zones(self, ndvi_map, grid_size=10):
        """Divide field into zones and compute statistics."""
        h, w = ndvi_map.shape
        zones = {}

        for i in range(0, h, grid_size):
            for j in range(0, w, grid_size):
                zone = ndvi_map[i:i+grid_size, j:j+grid_size]
                if zone.size > 0:
                    zone_id = f"zone_{i//grid_size}_{j//grid_size}"
                    zones[zone_id] = {
                        "mean_ndvi": float(np.mean(zone)),
                        "std_ndvi": float(np.std(zone)),
                        "min_ndvi": float(np.min(zone)),
                        "max_ndvi": float(np.max(zone)),
                        "pixels": int(zone.size),
                        "health_status": self._classify_zone(np.mean(zone))
                    }

        return zones

    def _classify_zone(self, mean_ndvi):
        """Classify zone health status."""
        if mean_ndvi >= 0.6:
            return "Excellent"
        elif mean_ndvi >= 0.4:
            return "Good"
        elif mean_ndvi >= 0.2:
            return "Fair"
        elif mean_ndvi >= 0.1:
            return "Poor"
        else:
            return "Critical"

    def generate_recommendations(self, zones):
        """Generate management recommendations."""
        recommendations = []

        poor_zones = [z for z, data in zones.items() if data["mean_ndvi"] < 0.3]
        if poor_zones:
            recommendations.append({
                "priority": "HIGH",
                "action": "Irrigate or fertilize",
                "zones": poor_zones,
                "reason": "NDVI below 0.3 indicates stressed vegetation"
            })

        variable_zones = [z for z, data in zones.items() if data["std_ndvi"] > 0.15]
        if variable_zones:
            recommendations.append({
                "priority": "MEDIUM",
                "action": "Investigate variability",
                "zones": variable_zones,
                "reason": "High NDVI variance within zone"
            })

        return recommendations

# Demo
mapper = CropHealthMapper()

np.random.seed(42)
ndvi_map = np.random.normal(0.5, 0.15, (100, 100))
ndvi_map[20:40, 60:80] = 0.15  # Stressed area

zones = mapper.create_health_zones(ndvi_map, grid_size=20)
recommendations = mapper.generate_recommendations(zones)

print(f"Total zones: {len(zones)}")
for zone_id, data in list(zones.items())[:5]:
    print(f"  {zone_id}: {data['health_status']} (NDVI: {data['mean_ndvi']:.3f})")

print(f"\nRecommendations: {len(recommendations)}")
for rec in recommendations:
    print(f"  [{rec['priority']}] {rec['action']}")