API Satellitari e Meteorologiche per AgriTech: Dati Predittivi

EVI = G * (NIR - Red) / (NIR + C1*Red - C2*Blue + L)

Su Sentinel-2:
EVI = 2.5 * (B8 - B4) / (B8 + 6*B4 - 7.5*B2 + 1)

Costanti standard:
  G  = 2.5  (guadagno)
  C1 = 6.0  (coefficiente resistenza aerosol)
  C2 = 7.5  (coefficiente resistenza aerosol)
  L  = 1.0  (fattore aggiustamento canopy)

Range: -1.0 a +1.0 (tipicamente 0 a 0.9 per vegetazione)

SAVI = (NIR - Red) * (1 + L) / (NIR + Red + L)

Su Sentinel-2:
SAVI = (B8 - B4) * (1 + 0.5) / (B8 + B4 + 0.5)

L = 0.5 (valore standard per vegetazione media)
L = 1.0 per vegetazione molto rada
L = 0.25 per vegetazione densa

Vantaggioso rispetto a NDVI quando:
- Copertura vegetale < 40%
- Suoli chiari e brillanti (calcari, sabbie)
- Areali semi-aridi del Sud Italia

# Stima LAI empirica da NDVI (relazione di Boegh et al., corretta per Sentinel-2)
# Valida per cereali e colture erbacee europee

LAI_stima = -0.3 + 10.2 * NDVI  # valida per NDVI in [0.2, 0.8]

# Per la vite (relazione specifica da letteratura italiana):
LAI_vite = 0.57 * EXP(3.28 * NDVI)  # Dalla et al., 2019

# Soglie LAI indicative per colture italiane:
# Grano: LAI ottimale 3-5 m2/m2 alla spigatura
# Mais:  LAI ottimale 3-5 m2/m2 alla fioritura
# Vite:  LAI ottimale 1.5-2.5 m2/m2 durante invaiatura
# Pomodoro: LAI ottimale 3-4 m2/m2 durante copertura massima

# Installa le librerie necessarie
pip install sentinelhub rasterio numpy pandas geopandas shapely \
            requests python-dotenv matplotlib scipy scikit-learn \
            openmeteo-requests retry-requests xarray

# requirements.txt con versioni validate per Python 3.11+
# sentinelhub==3.11.4      - interfaccia CDSE e SentinelHub
# rasterio==1.4.0          - lettura/scrittura raster geospaziali
# numpy==2.1.0             - calcolo matriciale per bande
# geopandas==1.0.1         - dati vettoriali e spatial operations
# shapely==2.0.6           - geometrie per AOI (Area of Interest)
# openmeteo-requests==0.3.3 - client Open-Meteo API
# scikit-learn==1.5.2      - modello predittivo stress idrico

import os
import requests
from datetime import datetime, timedelta
from dotenv import load_dotenv

load_dotenv()

# Configurazione credenziali CDSE (da variabili d'ambiente)
CDSE_USERNAME = os.getenv("CDSE_USERNAME")
CDSE_PASSWORD = os.getenv("CDSE_PASSWORD")
CDSE_CLIENT_ID = "cdse-public"
CDSE_TOKEN_URL = "https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token"

def get_cdse_access_token(username: str, password: str) -> str:
    """
    Ottieni access token OAuth2 dal CDSE Identity Server.
    Il token ha durata di 600 secondi (10 minuti).
    """
    response = requests.post(
        CDSE_TOKEN_URL,
        data={
            "grant_type": "password",
            "client_id": CDSE_CLIENT_ID,
            "username": username,
            "password": password,
        },
        headers={"Content-Type": "application/x-www-form-urlencoded"},
        timeout=30,
    )
    response.raise_for_status()
    token_data = response.json()
    return token_data["access_token"]


class CDSESession:
    """Sessione autenticata con gestione automatica del refresh token."""

    def __init__(self, username: str, password: str):
        self._username = username
        self._password = password
        self._token: str | None = None
        self._token_expiry: datetime | None = None

    def _refresh_token(self) -> None:
        self._token = get_cdse_access_token(self._username, self._password)
        # Margine di 60 secondi prima della scadenza reale (600s)
        self._token_expiry = datetime.utcnow() + timedelta(seconds=540)

    def get_headers(self) -> dict:
        if self._token is None or datetime.utcnow() >= self._token_expiry:
            self._refresh_token()
        return {"Authorization": f"Bearer {self._token}"}


# Inizializza sessione
session = CDSESession(CDSE_USERNAME, CDSE_PASSWORD)

from shapely.geometry import box
import json

ODATA_BASE_URL = "https://catalogue.dataspace.copernicus.eu/odata/v1"

def search_sentinel2_products(
    bbox: tuple[float, float, float, float],  # (lon_min, lat_min, lon_max, lat_max)
    start_date: str,   # formato: "2024-05-01T00:00:00.000Z"
    end_date: str,     # formato: "2024-05-31T23:59:59.000Z"
    cloud_cover_max: float = 20.0,
    product_type: str = "S2MSI2A",  # L2A = Bottom of Atmosphere, preferire per agricoltura
    max_results: int = 20,
) -> list[dict]:
    """
    Cerca prodotti Sentinel-2 nel catalogo CDSE.

    Args:
        bbox: Bounding box nell'ordine (lon_min, lat_min, lon_max, lat_max)
        start_date: Data inizio (formato ISO 8601)
        end_date: Data fine (formato ISO 8601)
        cloud_cover_max: Percentuale massima di copertura nuvolosa (0-100)
        product_type: S2MSI2A (L2A, consigliato) o S2MSI1C (L1C)
        max_results: Numero massimo di risultati

    Returns:
        Lista di dizionari con metadata prodotto
    """
    lon_min, lat_min, lon_max, lat_max = bbox
    aoi_wkt = f"POLYGON(( {lon_min} {lat_min}, {lon_max} {lat_min}, {lon_max} {lat_max}, {lon_min} {lat_max}, {lon_min} {lat_min} ))"

    filter_query = (
        f"Collection/Name eq 'SENTINEL-2' "
        f"and Attributes/OData.CSC.DoubleAttribute/any(att:att/Name eq 'cloudCover' "
        f"    and att/OData.CSC.DoubleAttribute/Value le {cloud_cover_max}) "
        f"and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' "
        f"    and att/OData.CSC.StringAttribute/Value eq '{product_type}') "
        f"and ContentDate/Start ge {start_date} "
        f"and ContentDate/Start lt {end_date} "
        f"and OData.CSC.Intersects(area=geography'SRID=4326;{aoi_wkt}')"
    )

    params = {
        "$filter": filter_query,
        "$orderby": "ContentDate/Start desc",
        "$top": max_results,
        "$expand": "Attributes",
    }

    response = requests.get(
        f"{ODATA_BASE_URL}/Products",
        params=params,
        timeout=60,
    )
    response.raise_for_status()
    data = response.json()
    return data.get("value", [])


# Esempio: cerca immagini per la vigna di Manduria (Puglia)
# Campagna estiva: maggio-settembre 2024
bbox_manduria = (17.4, 40.3, 17.7, 40.5)  # Zona DOC Primitivo di Manduria

products = search_sentinel2_products(
    bbox=bbox_manduria,
    start_date="2024-05-01T00:00:00.000Z",
    end_date="2024-09-30T23:59:59.000Z",
    cloud_cover_max=10.0,
    product_type="S2MSI2A",
)

print(f"Trovati {len(products)} prodotti Sentinel-2 L2A con copertura nuvolosa <= 10%")
for p in products[:5]:
    cloud = next(
        (a["Value"] for a in p.get("Attributes", []) if a["Name"] == "cloudCover"),
        "N/A",
    )
    print(f"  {p['Name']} | Cloud: {cloud:.1f}% | Date: {p['ContentDate']['Start'][:10]}")

from sentinelhub import (
    SHConfig,
    SentinelHubRequest,
    DataCollection,
    MimeType,
    BBox,
    CRS,
    bbox_to_dimensions,
)
import numpy as np

# Configurazione SentinelHub via CDSE
config = SHConfig()
config.sh_client_id = os.getenv("SH_CLIENT_ID")       # da apps.sentinel-hub.com
config.sh_client_secret = os.getenv("SH_CLIENT_SECRET")
config.sh_base_url = "https://sh.dataspace.copernicus.eu"  # endpoint CDSE

# Definizione Area of Interest (AOI) - Vigna Manduria 100 ha circa
aoi_bbox = BBox(
    bbox=[17.42, 40.34, 17.52, 40.42],
    crs=CRS.WGS84,
)
# Calcola dimensioni raster a 10m di risoluzione
size = bbox_to_dimensions(aoi_bbox, resolution=10)
print(f"Dimensioni raster: {size[0]} x {size[1]} pixel a 10m")

# Evalscript per scaricare B04 (Red) e B08 (NIR) per calcolo NDVI
evalscript_ndvi_bands = """
//VERSION=3
function setup() {
    return {
        input: [{
            bands: ["B04", "B08", "CLM"],  // Red, NIR, Cloud Mask
            units: "REFLECTANCE"
        }],
        output: {
            bands: 3,
            sampleType: "FLOAT32"
        }
    };
}

function evaluatePixel(sample) {
    // Restituisce [Red, NIR, CloudMask]
    return [sample.B04, sample.B08, sample.CLM];
}
"""

# Richiesta dati
request = SentinelHubRequest(
    evalscript=evalscript_ndvi_bands,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A.define_from(
                "s2l2a",
                service_url=config.sh_base_url,
            ),
            time_interval=("2024-07-15", "2024-07-20"),
            other_args={"dataFilter": {"maxCloudCoverage": 10}},
        )
    ],
    responses=[SentinelHubRequest.output_response("default", MimeType.TIFF)],
    bbox=aoi_bbox,
    size=size,
    config=config,
)

# Download dati (restituisce array numpy [H, W, Bande])
images = request.get_data()
band_data = images[0]  # shape: (height, width, 3)

red_band  = band_data[:, :, 0].astype(np.float32)  # B04 - Red
nir_band  = band_data[:, :, 1].astype(np.float32)  # B08 - NIR
cloud_mask = band_data[:, :, 2]                     # CLM - 0=clear, 1=cloud

print(f"Bande scaricate - Shape: {red_band.shape}")
print(f"Red B04 - Min: {red_band.min():.4f}, Max: {red_band.max():.4f}, Mean: {red_band.mean():.4f}")
print(f"NIR B08 - Min: {nir_band.min():.4f}, Max: {nir_band.max():.4f}, Mean: {nir_band.mean():.4f}")
print(f"Pixel nuvolosi: {cloud_mask.sum()} su {cloud_mask.size} totali ({100*cloud_mask.mean():.1f}%)")

import rasterio
from rasterio.transform import from_bounds
from rasterio.crs import CRS as RasterioCRS
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

def calculate_ndvi(
    red: np.ndarray,
    nir: np.ndarray,
    cloud_mask: np.ndarray | None = None,
) -> np.ndarray:
    """
    Calcola NDVI dalla banda rossa e NIR.
    Gestisce divisione per zero con np.errstate.
    Maschera i pixel nuvolosi con np.nan.

    Returns:
        Array NDVI con float32, range [-1, 1], NaN per pixel invalidi
    """
    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = np.where(
            (nir + red) == 0,
            np.nan,
            (nir - red) / (nir + red),
        ).astype(np.float32)

    # Maschera nuvole
    if cloud_mask is not None:
        ndvi = np.where(cloud_mask == 1, np.nan, ndvi)

    # Clip valori fisicamente impossibili (artefatti strumentali)
    ndvi = np.clip(ndvi, -1.0, 1.0)
    return ndvi


def save_ndvi_geotiff(
    ndvi: np.ndarray,
    bbox_wgs84: tuple[float, float, float, float],
    output_path: str,
) -> None:
    """
    Salva l'array NDVI come GeoTIFF georeferenziato in WGS84.

    Args:
        ndvi: Array NDVI 2D (H, W)
        bbox_wgs84: (lon_min, lat_min, lon_max, lat_max)
        output_path: Path del file output
    """
    height, width = ndvi.shape
    lon_min, lat_min, lon_max, lat_max = bbox_wgs84

    transform = from_bounds(lon_min, lat_min, lon_max, lat_max, width, height)

    with rasterio.open(
        output_path,
        "w",
        driver="GTiff",
        height=height,
        width=width,
        count=1,
        dtype=rasterio.float32,
        crs=RasterioCRS.from_epsg(4326),
        transform=transform,
        compress="lzw",
        nodata=np.nan,
    ) as dst:
        dst.write(ndvi, 1)
        dst.update_tags(
            NDVI_COMPUTED_AT=datetime.utcnow().isoformat(),
            SENTINEL2_PRODUCT_TYPE="S2MSI2A",
            FORMULA="(B08-B04)/(B08+B04)",
        )


# Calcola e salva NDVI
ndvi_map = calculate_ndvi(red_band, nir_band, cloud_mask)

# Statistiche NDVI per l'AOI (escluso NaN)
valid_pixels = ndvi_map[~np.isnan(ndvi_map)]
print(f"\nStatistiche NDVI - Vigna Manduria (luglio 2024)")
print(f"  Pixel validi: {len(valid_pixels)} su {ndvi_map.size}")
print(f"  NDVI medio:   {valid_pixels.mean():.3f}")
print(f"  NDVI mediano: {np.median(valid_pixels):.3f}")
print(f"  NDVI p10:     {np.percentile(valid_pixels, 10):.3f} (zone a stress)")
print(f"  NDVI p90:     {np.percentile(valid_pixels, 90):.3f} (zone ottimali)")

bbox_wgs84 = (17.42, 40.34, 17.52, 40.42)
save_ndvi_geotiff(ndvi_map, bbox_wgs84, "ndvi_manduria_20240717.tif")
print("GeoTIFF NDVI salvato: ndvi_manduria_20240717.tif")

import openmeteo_requests
import requests_cache
import pandas as pd
from retry_requests import retry
from dataclasses import dataclass
from typing import Optional

@dataclass
class WeatherForecast:
    """Dati meteo giornalieri per decisioni agronomiche."""
    date: pd.DatetimeIndex
    temperature_max: np.ndarray    # Celsius
    temperature_min: np.ndarray    # Celsius
    precipitation: np.ndarray      # mm
    et0: np.ndarray                # mm/giorno - evapotraspirazione di riferimento (FAO-56)
    wind_speed_max: np.ndarray     # km/h
    solar_radiation: np.ndarray    # MJ/m2/giorno
    soil_moisture_0_7cm: np.ndarray   # m3/m3 (0-7 cm profondità)
    soil_moisture_7_28cm: np.ndarray  # m3/m3 (7-28 cm profondità)
    soil_temp_0_7cm: np.ndarray    # Celsius


def get_weather_forecast(
    latitude: float,
    longitude: float,
    start_date: Optional[str] = None,   # "YYYY-MM-DD" per archivio
    end_date: Optional[str] = None,
    forecast_days: int = 16,
) -> WeatherForecast:
    """
    Recupera dati meteo da Open-Meteo con cache locale (1 ora) e retry automatico.
    Combina forecast (16 giorni) con archivio storico ERA5 se start_date specificato.

    Args:
        latitude: Latitudine WGS84
        longitude: Longitudine WGS84
        start_date: Se fornito, richiede dati storici (non forecast)
        end_date: Fine periodo storico
        forecast_days: Giorni di forecast (1-16)

    Returns:
        WeatherForecast con tutti i parametri agrometeorologici
    """
    # Cache HTTP locale (1 ora di validita) + retry su errori transitori
    cache_session = requests_cache.CachedSession(".cache/open_meteo", expire_after=3600)
    retry_session = retry(cache_session, retries=5, backoff_factor=0.2)
    client = openmeteo_requests.Client(session=retry_session)

    base_url = (
        "https://archive-api.open-meteo.com/v1/archive"
        if start_date
        else "https://api.open-meteo.com/v1/forecast"
    )

    params = {
        "latitude": latitude,
        "longitude": longitude,
        "daily": [
            "temperature_2m_max",
            "temperature_2m_min",
            "precipitation_sum",
            "et0_fao_evapotranspiration",
            "wind_speed_10m_max",
            "shortwave_radiation_sum",
            "soil_moisture_0_to_7cm",
            "soil_moisture_7_to_28cm",
            "soil_temperature_0_to_7cm",
        ],
        "timezone": "Europe/Rome",  # Fuso orario italiano
        "wind_speed_unit": "kmh",
    }

    if start_date:
        params["start_date"] = start_date
        params["end_date"] = end_date
    else:
        params["forecast_days"] = forecast_days

    # Open-Meteo usa ERA5 per archivio storico (1940-oggi), eccellente per agricoltura
    if start_date:
        # Archivio ERA5-Land: alta risoluzione 9 km per Italia
        params["models"] = "era5_land"

    responses = client.weather_api(base_url, params=params)
    r = responses[0]  # primo (e unico) punto

    daily = r.Daily()

    return WeatherForecast(
        date=pd.date_range(
            start=pd.to_datetime(daily.Time(), unit="s", utc=True),
            end=pd.to_datetime(daily.TimeEnd(), unit="s", utc=True),
            freq=pd.Timedelta(seconds=daily.Interval()),
            inclusive="left",
        ),
        temperature_max      = daily.Variables(0).ValuesAsNumpy(),
        temperature_min      = daily.Variables(1).ValuesAsNumpy(),
        precipitation        = daily.Variables(2).ValuesAsNumpy(),
        et0                  = daily.Variables(3).ValuesAsNumpy(),
        wind_speed_max       = daily.Variables(4).ValuesAsNumpy(),
        solar_radiation      = daily.Variables(5).ValuesAsNumpy(),
        soil_moisture_0_7cm  = daily.Variables(6).ValuesAsNumpy(),
        soil_moisture_7_28cm = daily.Variables(7).ValuesAsNumpy(),
        soil_temp_0_7cm      = daily.Variables(8).ValuesAsNumpy(),
    )


# Esempio: meteo vigna Manduria (campagna 2024)
lat_manduria, lon_manduria = 40.38, 17.47

weather_storico = get_weather_forecast(
    latitude=lat_manduria,
    longitude=lon_manduria,
    start_date="2024-04-01",
    end_date="2024-09-30",
)
weather_forecast = get_weather_forecast(
    latitude=lat_manduria,
    longitude=lon_manduria,
    forecast_days=16,
)

df_meteo = pd.DataFrame({
    "data":              weather_storico.date,
    "t_max_C":           weather_storico.temperature_max,
    "t_min_C":           weather_storico.temperature_min,
    "pioggia_mm":        weather_storico.precipitation,
    "et0_mm":            weather_storico.et0,
    "rad_MJ":            weather_storico.solar_radiation,
    "sw_0_7cm":          weather_storico.soil_moisture_0_7cm,
    "sw_7_28cm":         weather_storico.soil_moisture_7_28cm,
})

print(df_meteo.describe())
print(f"\nET0 totale campagna apr-set 2024: {df_meteo['et0_mm'].sum():.1f} mm")
print(f"Pioggia totale campagna apr-set 2024: {df_meteo['pioggia_mm'].sum():.1f} mm")
deficit_idrico = df_meteo["et0_mm"].sum() - df_meteo["pioggia_mm"].sum()
print(f"Deficit idrico stagionale (ET0 - pioggia): {deficit_idrico:.1f} mm")

import geopandas as gpd
from rasterstats import zonal_stats
import rasterio
from sqlalchemy import create_engine, text
from geoalchemy2 import Geometry
import pandas as pd

# Carica shapefile appezzamenti vigna
gdf_parcelle = gpd.read_file("dati/parcelle_vigna_manduria.shp")
gdf_parcelle = gdf_parcelle.to_crs("EPSG:4326")  # Riproietta in WGS84

print(f"Appezzamenti caricati: {len(gdf_parcelle)}")
print(gdf_parcelle[["id_parcella", "varieta", "ettari", "anno_impianto"]].head(10))

def compute_zonal_ndvi(
    ndvi_raster_path: str,
    parcelle_gdf: gpd.GeoDataFrame,
    acquisition_date: str,
) -> gpd.GeoDataFrame:
    """
    Calcola statistiche NDVI zonali per ogni appezzamento.

    Per ogni parcella calcola:
    - NDVI medio (indicatore di vigor generale)
    - NDVI mediano (robusto agli outlier)
    - NDVI std (indica variabilità intra-parcella)
    - Percentile 10 (zone a stress all'interno della parcella)
    - Percentile 90 (zone ottimali)
    - Percentuale pixel validi (esclude NaN da nuvole)

    Returns:
        GeoDataFrame con colonne NDVI aggiunte
    """
    stats = zonal_stats(
        vectors=parcelle_gdf.geometry,
        raster=ndvi_raster_path,
        stats=["mean", "median", "std", "percentile_10", "percentile_90", "count", "nodata"],
        nodata=np.nan,
        geojson_out=False,
    )

    df_stats = pd.DataFrame(stats)
    df_stats.columns = [
        "ndvi_mean", "ndvi_median", "ndvi_std",
        "ndvi_p10", "ndvi_p90", "pixel_count", "pixel_nodata",
    ]
    df_stats["acquisition_date"] = pd.to_datetime(acquisition_date)
    df_stats["pixel_valid_pct"] = (
        df_stats["pixel_count"] /
        (df_stats["pixel_count"] + df_stats["pixel_nodata"].fillna(0)) * 100
    )

    # Classificazione stress basata su NDVI medio (soglie per vite in estate)
    df_stats["stress_category"] = pd.cut(
        df_stats["ndvi_mean"],
        bins=[-1.0, 0.25, 0.40, 0.55, 0.70, 1.0],
        labels=["Stress Severo", "Stress Moderato", "Normale", "Buono", "Ottimale"],
    )

    return gpd.GeoDataFrame(
        pd.concat([parcelle_gdf.reset_index(drop=True), df_stats], axis=1),
        geometry="geometry",
        crs=parcelle_gdf.crs,
    )


# Calcola NDVI zonale per tutte le parcelle
gdf_ndvi = compute_zonal_ndvi(
    ndvi_raster_path="ndvi_manduria_20240717.tif",
    parcelle_gdf=gdf_parcelle,
    acquisition_date="2024-07-17",
)

# Parcelle in stress - priorità irrigazione
parcelle_in_stress = gdf_ndvi[
    gdf_ndvi["stress_category"].isin(["Stress Severo", "Stress Moderato"])
]
print(f"\nParcelle in stress ({len(parcelle_in_stress)}/{len(gdf_ndvi)}):")
print(parcelle_in_stress[["id_parcella", "varieta", "ndvi_mean", "stress_category", "ettari"]])

# Salva in PostGIS per persistenza e query spaziali
engine = create_engine("postgresql://agritech:pass@localhost:5432/vigna_db")

gdf_ndvi.to_postgis(
    name="ndvi_storico",
    con=engine,
    if_exists="append",
    index=False,
    dtype={"geometry": Geometry("MULTIPOLYGON", srid=4326)},
)
print("\nDati NDVI salvati in PostGIS (tabella: ndvi_storico)")

-- Schema PostGIS per sistema AgriTech
CREATE TABLE IF NOT EXISTS parcelle (
    id_parcella    SERIAL PRIMARY KEY,
    codice_catasto VARCHAR(20) UNIQUE NOT NULL,  -- codice catastale italiano
    varieta        VARCHAR(50),
    ettari         NUMERIC(8, 4),
    anno_impianto  INTEGER,
    sistema_allevamento VARCHAR(30),
    geom           GEOMETRY(MULTIPOLYGON, 4326)
);

CREATE INDEX idx_parcelle_geom ON parcelle USING GIST(geom);

CREATE TABLE IF NOT EXISTS ndvi_storico (
    id            SERIAL PRIMARY KEY,
    id_parcella   INTEGER REFERENCES parcelle(id_parcella),
    data_satellite DATE NOT NULL,
    ndvi_mean     NUMERIC(6, 4),
    ndvi_median   NUMERIC(6, 4),
    ndvi_std      NUMERIC(6, 4),
    ndvi_p10      NUMERIC(6, 4),
    ndvi_p90      NUMERIC(6, 4),
    pixel_valid_pct NUMERIC(5, 2),
    stress_category VARCHAR(20),
    satellite     VARCHAR(20) DEFAULT 'Sentinel-2',
    created_at    TIMESTAMPTZ DEFAULT NOW()
);

CREATE INDEX idx_ndvi_parcella_data ON ndvi_storico(id_parcella, data_satellite);

-- Query: tendenza NDVI ultime 4 acquisizioni per parcella
SELECT
    p.codice_catasto,
    p.varieta,
    n.data_satellite,
    n.ndvi_mean,
    LAG(n.ndvi_mean, 1) OVER (PARTITION BY n.id_parcella ORDER BY n.data_satellite) AS ndvi_prev,
    n.ndvi_mean - LAG(n.ndvi_mean, 1) OVER (PARTITION BY n.id_parcella ORDER BY n.data_satellite) AS ndvi_delta
FROM parcelle p
JOIN ndvi_storico n ON p.id_parcella = n.id_parcella
WHERE n.data_satellite >= NOW() - INTERVAL '60 days'
ORDER BY p.id_parcella, n.data_satellite;

-- Query: parcelle con NDVI in calo > 0.1 nell'ultimo mese (alert stress)
SELECT
    p.codice_catasto,
    p.varieta,
    p.ettari,
    latest.ndvi_mean AS ndvi_corrente,
    prev.ndvi_mean   AS ndvi_30gg_fa,
    (latest.ndvi_mean - prev.ndvi_mean) AS variazione
FROM parcelle p
JOIN ndvi_storico latest ON p.id_parcella = latest.id_parcella
    AND latest.data_satellite = (SELECT MAX(data_satellite) FROM ndvi_storico WHERE id_parcella = p.id_parcella)
JOIN ndvi_storico prev ON p.id_parcella = prev.id_parcella
    AND prev.data_satellite = (SELECT MAX(data_satellite) FROM ndvi_storico
        WHERE id_parcella = p.id_parcella AND data_satellite < NOW() - INTERVAL '25 days')
WHERE (latest.ndvi_mean - prev.ndvi_mean) < -0.10
ORDER BY variazione ASC;

import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report
import joblib
from typing import Tuple

def build_stress_features(
    df_ndvi: pd.DataFrame,      # colonne: data, ndvi_mean, ndvi_std, ndvi_p10
    df_meteo: pd.DataFrame,     # colonne: data, t_max_C, t_min_C, pioggia_mm, et0_mm, sw_0_7cm
    lookback_days: int = 14,
) -> pd.DataFrame:
    """
    Costruisce feature set per modello predittivo stress idrico.
    Combina NDVI satellite con bilancio idrico e meteo.

    Features generate:
    - NDVI corrente e variazione a 5/10 giorni
    - Deficit idrico cumulato (ET0 - pioggia) a 7/14/21 giorni
    - Growing Degree Days (GDD) da inizio stagione
    - Soil moisture media e tendenza
    - Variabilità termica (t_max - t_min)
    """
    df = df_meteo.merge(df_ndvi, on="data", how="left").sort_values("data")
    df["ndvi_mean"] = df["ndvi_mean"].interpolate(method="linear", limit=5)

    # Bilancio idrico cumulato (mm) - positivo = surplus, negativo = deficit
    df["bilancio_idrico"] = df["pioggia_mm"] - df["et0_mm"]
    df["deficit_7gg"]  = df["bilancio_idrico"].rolling(7).sum()
    df["deficit_14gg"] = df["bilancio_idrico"].rolling(14).sum()
    df["deficit_21gg"] = df["bilancio_idrico"].rolling(21).sum()

    # NDVI trend (variazione puntuale e su finestra)
    df["ndvi_delta_5gg"]  = df["ndvi_mean"] - df["ndvi_mean"].shift(5)
    df["ndvi_delta_10gg"] = df["ndvi_mean"] - df["ndvi_mean"].shift(10)

    # Growing Degree Days (base 10 gradi - standard per vite)
    df["t_media"] = (df["t_max_C"] + df["t_min_C"]) / 2
    df["gdd_giornaliero"] = np.maximum(df["t_media"] - 10, 0)
    df["gdd_cumulato"] = df["gdd_giornaliero"].cumsum()

    # Variabilità termica (indicatore stress termico)
    df["escursione_termica"] = df["t_max_C"] - df["t_min_C"]

    # Soil moisture media e trend
    df["sm_media_7gg"] = df["sw_0_7cm"].rolling(7).mean()
    df["sm_trend"]     = df["sw_0_7cm"] - df["sw_0_7cm"].shift(3)

    feature_cols = [
        "ndvi_mean", "ndvi_std", "ndvi_p10",
        "ndvi_delta_5gg", "ndvi_delta_10gg",
        "deficit_7gg", "deficit_14gg", "deficit_21gg",
        "et0_mm", "pioggia_mm",
        "t_max_C", "t_min_C", "escursione_termica",
        "gdd_cumulato", "gdd_giornaliero",
        "sm_media_7gg", "sm_trend",
    ]

    return df[["data"] + feature_cols].dropna()


def train_stress_model(
    features_df: pd.DataFrame,
    labels: pd.Series,  # 0=no stress, 1=stress moderato, 2=stress severo
) -> Tuple[Pipeline, dict]:
    """
    Addestra modello GradientBoosting per classificazione stress idrico.
    Usa pipeline sklearn con StandardScaler integrato.

    Returns:
        (pipeline_addestrata, metriche_validazione)
    """
    feature_cols = [c for c in features_df.columns if c != "data"]
    X = features_df[feature_cols].values
    y = labels.values

    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42, stratify=y
    )

    pipeline = Pipeline([
        ("scaler", StandardScaler()),
        ("clf", GradientBoostingClassifier(
            n_estimators=200,
            learning_rate=0.05,
            max_depth=4,
            subsample=0.8,
            random_state=42,
        )),
    ])

    pipeline.fit(X_train, y_train)
    y_pred = pipeline.predict(X_test)

    cv_scores = cross_val_score(pipeline, X, y, cv=5, scoring="f1_weighted")

    metriche = {
        "f1_cv_mean": cv_scores.mean(),
        "f1_cv_std":  cv_scores.std(),
        "classification_report": classification_report(
            y_test, y_pred,
            target_names=["No Stress", "Stress Moderato", "Stress Severo"],
        ),
    }

    return pipeline, metriche


# Salva modello per deployment
def save_model(pipeline: Pipeline, path: str) -> None:
    joblib.dump(pipeline, path)
    print(f"Modello salvato: {path}")


def predict_stress_7days(
    pipeline: Pipeline,
    features_forecast: pd.DataFrame,  # features calcolate su dati forecast 7 giorni
) -> pd.DataFrame:
    """
    Predici stress idrico per i prossimi 7 giorni.
    Usa le features calcolate su dati forecast meteo (Open-Meteo, 16 gg).

    Returns:
        DataFrame con data, probabilità per classe, classe predetta
    """
    feature_cols = [c for c in features_forecast.columns if c != "data"]
    X_forecast = features_forecast[feature_cols].values

    proba = pipeline.predict_proba(X_forecast)
    pred_class = pipeline.predict(X_forecast)

    labels = ["no_stress", "stress_moderato", "stress_severo"]
    result_df = pd.DataFrame(proba, columns=[f"prob_{l}" for l in labels])
    result_df["classe_predetta"] = pred_class
    result_df["data"] = features_forecast["data"].values
    result_df["label_predetto"] = [labels[c] for c in pred_class]

    return result_df[["data", "label_predetto", "prob_no_stress",
                       "prob_stress_moderato", "prob_stress_severo"]]

import smtplib
from email.mime.text import MIMEText
from email.mime.multipart import MIMEMultipart
from typing import List

def check_ndvi_alerts(
    gdf_ndvi: gpd.GeoDataFrame,
    stress_threshold: float = 0.35,
    min_valid_pixels_pct: float = 70.0,
) -> List[dict]:
    """
    Controlla le parcelle con NDVI sotto soglia di stress.
    Genera lista di alert per notifica agronomo.

    Args:
        gdf_ndvi: GeoDataFrame con colonna ndvi_mean e pixel_valid_pct
        stress_threshold: Soglia NDVI sotto cui scattare alert (default 0.35 per vite)
        min_valid_pixels_pct: Percentuale minima pixel validi per affidabilità

    Returns:
        Lista di dict con dettagli alert per parcella
    """
    alerts = []
    for _, row in gdf_ndvi.iterrows():
        if row["pixel_valid_pct"] < min_valid_pixels_pct:
            continue  # Troppi pixel nuvolosi, dato inaffidabile

        if row["ndvi_mean"] < stress_threshold:
            severity = "SEVERO" if row["ndvi_mean"] < 0.25 else "MODERATO"
            alerts.append({
                "id_parcella":   row["id_parcella"],
                "varieta":       row.get("varieta", "N/A"),
                "ettari":        row.get("ettari", 0),
                "ndvi_corrente": round(row["ndvi_mean"], 3),
                "ndvi_soglia":   stress_threshold,
                "severity":      severity,
                "data":          row.get("acquisition_date", "N/A"),
                "pixel_valid":   round(row["pixel_valid_pct"], 1),
            })

    return sorted(alerts, key=lambda x: x["ndvi_corrente"])


def send_alert_email(
    alerts: List[dict],
    recipient: str,
    sender: str,
    smtp_host: str = "smtp.gmail.com",
    smtp_port: int = 587,
) -> None:
    """Invia email di alert NDVI all'agronomo."""
    if not alerts:
        return

    body_lines = [
        f"Alert Monitoraggio NDVI Satellite - {alerts[0]['data']}\n",
        f"Parcelle in stress idrico: {len(alerts)}\n",
        "=" * 60,
        "",
    ]
    for a in alerts:
        body_lines.append(
            f"  [{a['severity']}] Parcella {a['id_parcella']} ({a['varietà]}, {a['ettari']} ha)\n"
            f"    NDVI: {a['ndvi_corrente']} (soglia: {a['ndvi_soglia']}) | Pixel validi: {a['pixel_valid']}%\n"
        )

    msg = MIMEMultipart()
    msg["Subject"] = f"[ALERT NDVI] {len(alerts)} parcelle in stress - {alerts[0]['data']}"
    msg["From"]    = sender
    msg["To"]      = recipient
    msg.attach(MIMEText("\n".join(body_lines), "plain"))

    smtp_password = os.getenv("SMTP_PASSWORD")
    with smtplib.SMTP(smtp_host, smtp_port) as server:
        server.starttls()
        server.login(sender, smtp_password)
        server.send_message(msg)
        print(f"Alert inviato a {recipient} per {len(alerts)} parcelle in stress")

import ee

# Autenticazione Google Earth Engine (richiede account GEE)
ee.Authenticate()
ee.Initialize(project="your-gee-project-id")

def calculate_ndvi_gee(
    aoi_geojson: dict,
    start_date: str,
    end_date: str,
    cloud_threshold: int = 20,
) -> ee.ImageCollection:
    """
    Calcola NDVI su larga scala usando Google Earth Engine.
    Adatto per analisi su comprensori DOC o scala regionale.

    Args:
        aoi_geojson: GeoJSON dell'area di interesse
        start_date: "YYYY-MM-DD"
        end_date:   "YYYY-MM-DD"
        cloud_threshold: Percentuale max copertura nuvolosa (0-100)

    Returns:
        ImageCollection con NDVI calcolato per ogni scena
    """
    aoi = ee.Geometry(aoi_geojson)

    def mask_s2_clouds(image: ee.Image) -> ee.Image:
        """Maschera nuvole usando QA60 band di Sentinel-2."""
        qa = image.select("QA60")
        cloud_bit_mask = 1 << 10   # bit 10: nuvole opache
        cirrus_bit_mask = 1 << 11  # bit 11: nuvole cirro
        mask = qa.bitwiseAnd(cloud_bit_mask).eq(0).And(
               qa.bitwiseAnd(cirrus_bit_mask).eq(0))
        return image.updateMask(mask).divide(10000)  # scala a [0,1]

    def add_ndvi(image: ee.Image) -> ee.Image:
        """Aggiunge banda NDVI all'immagine Sentinel-2."""
        ndvi = image.normalizedDifference(["B8", "B4"]).rename("NDVI")
        return image.addBands(ndvi)

    collection = (
        ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
        .filterBounds(aoi)
        .filterDate(start_date, end_date)
        .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", cloud_threshold))
        .map(mask_s2_clouds)
        .map(add_ndvi)
    )

    return collection


# Calcola NDVI mediano per stagione (mosaico cloud-free)
def get_seasonal_ndvi_mosaic(
    collection: ee.ImageCollection,
    aoi: ee.Geometry,
) -> ee.Image:
    """Crea mosaico NDVI stagionale (mediano = robusto a outlier residui)."""
    ndvi_mosaic = collection.select("NDVI").median().clip(aoi)
    return ndvi_mosaic


# Esempio: Comprensorio DOC Primitivo di Manduria (~21.000 ha)
aoi_doc_manduria = {
    "type": "Polygon",
    "coordinates": [[
        [17.2, 40.2], [17.8, 40.2], [17.8, 40.6],
        [17.2, 40.6], [17.2, 40.2],
    ]],
}

collection_estate = calculate_ndvi_gee(
    aoi_geojson=aoi_doc_manduria,
    start_date="2024-07-01",
    end_date="2024-08-31",
    cloud_threshold=15,
)
print(f"Scene disponibili: {collection_estate.size().getInfo()}")

# Export NDVI mosaic su Google Drive (per analisi successive in Python locale)
ndvi_mosaic = get_seasonal_ndvi_mosaic(collection_estate, ee.Geometry(aoi_doc_manduria))
task = ee.batch.Export.image.toDrive(
    image=ndvi_mosaic,
    description="NDVI_DOC_Manduria_Estate2024",
    folder="AgriTech_Export",
    fileNamePrefix="ndvi_manduria_estate2024",
    region=ee.Geometry(aoi_doc_manduria),
    scale=10,
    crs="EPSG:4326",
    maxPixels=1e10,
)
task.start()
print(f"Export avviato - ID task: {task.id}")

from prefect import flow, task
from prefect.schedules import CronSchedule
from datetime import datetime, timedelta
import logging

logger = logging.getLogger(__name__)

@task(retries=3, retry_delay_seconds=60, name="search-sentinel2-products")
def search_products_task(bbox: tuple, lookback_days: int = 7) -> list[dict]:
    """Task Prefect: ricerca prodotti Sentinel-2 recenti."""
    end_date = datetime.utcnow()
    start_date = end_date - timedelta(days=lookback_days)
    products = search_sentinel2_products(
        bbox=bbox,
        start_date=start_date.strftime("%Y-%m-%dT00:00:00.000Z"),
        end_date=end_date.strftime("%Y-%m-%dT23:59:59.000Z"),
        cloud_cover_max=20.0,
    )
    logger.info(f"Trovati {len(products)} prodotti Sentinel-2")
    return products


@task(retries=2, retry_delay_seconds=120, name="download-compute-ndvi")
def compute_ndvi_task(
    product: dict,
    aoi_bbox: tuple,
    output_dir: str = "/data/ndvi",
) -> str:
    """Task Prefect: scarica bande e calcola NDVI per un prodotto."""
    # Implementazione basata su funzioni definite nelle sezioni precedenti
    acquisition_date = product["ContentDate"]["Start"][:10]
    output_path = f"{output_dir}/ndvi_{acquisition_date}.tif"

    # Download + calcolo (codice omesso per brevita, usa funzioni precedenti)
    logger.info(f"NDVI calcolato per {acquisition_date}: {output_path}")
    return output_path


@task(name="zonal-stats-postgis")
def zonal_stats_task(ndvi_path: str, parcelle_shapefile: str) -> dict:
    """Task Prefect: calcola statistiche zonali e salva in PostGIS."""
    gdf_parcelle = gpd.read_file(parcelle_shapefile).to_crs("EPSG:4326")
    acquisition_date = ndvi_path.split("ndvi_")[1].replace(".tif", "")

    gdf_ndvi = compute_zonal_ndvi(ndvi_path, gdf_parcelle, acquisition_date)

    engine = create_engine(os.getenv("POSTGIS_URL"))
    gdf_ndvi.to_postgis("ndvi_storico", engine, if_exists="append", index=False)

    alerts = check_ndvi_alerts(gdf_ndvi)
    return {"parcelle_processate": len(gdf_ndvi), "alerts": len(alerts), "alert_details": alerts}


@task(name="send-alerts")
def alerts_task(stats: dict) -> None:
    """Task Prefect: invia alert email se presenti."""
    if stats["alerts"] > 0:
        send_alert_email(
            alerts=stats["alert_details"],
            recipient=os.getenv("AGRONOMO_EMAIL"),
            sender=os.getenv("ALERT_EMAIL"),
        )
        logger.info(f"Alert inviati per {stats['alerts']} parcelle")


@flow(
    name="satellite-monitoring-pipeline",
    schedule=CronSchedule(cron="0 8 */5 * *"),  # ogni 5 giorni alle 8:00
)
def satellite_monitoring_flow(
    bbox: tuple = (17.35, 40.28, 17.60, 40.50),
    parcelle_shapefile: str = "/data/parcelle_vigna.shp",
) -> None:
    """
    Pipeline completa di monitoraggio satellite per vigna.
    Si esegue ogni 5 giorni, allineata alla frequenza Sentinel-2.
    """
    products = search_products_task(bbox=bbox)

    if not products:
        logger.warning("Nessun prodotto trovato nell'ultimo periodo")
        return

    # Processo il prodotto più recente con minima copertura nuvolosa
    best_product = sorted(
        products,
        key=lambda p: next(
            (a["Value"] for a in p.get("Attributes", []) if a["Name"] == "cloudCover"), 100
        ),
    )[0]

    ndvi_path = compute_ndvi_task(best_product, bbox)
    stats = zonal_stats_task(ndvi_path, parcelle_shapefile)
    alerts_task(stats)

    logger.info(
        f"Pipeline completata: {stats['parcelle_processate']} parcelle, {stats['alerts']} alert"
    )


if __name__ == "__main__":
    satellite_monitoring_flow()