Satelitní a meteorologická rozhraní API pro AgriTech: Prediktivní data s Sentinel-2 a Python
V roce 2024 Evropská kosmická agentura získala přes 1,5 petabajtu dat každý den ze satelitů Sentinel programu Copernicus. Z toho část věnovaná monitorování zemědělství v Evropě představuje největší a strategicky nejvýznamnější segment. Přesto velká většina italských zemědělských společností, včetně těch technologicky nejvyspělejších, zatím nevyužívá tuto veřejnou, bezplatnou a vědecky ověřenou infrastrukturu.
Důvodem není nedostatek dat: jde o složitost přístupu k nim, jejich zpracování a transformace konkrétní agronomická rozhodnutí. Satelit Sentinel-2 prochází nad vaší vinicí každých 5 dní rozlišení 10 metrů na pixel. Každý průchod vytváří 13 spektrálních pásem, která kódují zdraví vegetace, vodní stres, obsah chlorofylu, přítomnost parazitů. V kombinaci s daty počasí, pozemní senzory IoT a prediktivní modely mohou být tato data radikálně transformační kvalita agronomických rozhodnutí, snížení vstupních nákladů o 15-25 % a zvýšení výnosů až 10-20%.
Globální trh pozorování Země za to stojí 10,07 miliardy dolarů v roce 2025 podle Straits Research, s růstem očekávaným na 17,20 miliard do roku 2033 (CAGR 6,92 %). Segment zemědělství představuje 21 % aplikací a roste nejvyšším tempem, řízena poptávkou po precizním zemědělství, sledování výnosů a optimalizaci zavlažování. V tomto článku vytvoříme krok za krokem kompletní kanál Pythonu pro přístup k Sentinel-2 přes CDSE (Copernicus Data Space Ecosystem), vypočítat NDVI a další vegetační indexy, integrovat meteorologická data z Open-Meteo a poskytují předpovědní model pro vodní stres.
Co se dozvíte v tomto článku
- Architektura programu Copernicus a jak získat bezplatný přístup k Sentinel-2 přes CDSE
- Vegetační indexy: NDVI, EVI, SAVI, LAI - matematické vzorce, interpretace a prahové hodnoty pro italské plodiny
- Kompletní implementace Pythonu: CDSE autentizace, stahování pásma, výpočet NDVI s rasteriem a numpy
- Porovnaná rozhraní API pro počasí: Open-Meteo, OpenWeatherMap, Tomorrow.io a integrace satelitních dat
- Geoprostorové potrubí s GeoPandas, rasteriem a PostGIS pro vektorová a rastrová data
- Prediktivní model vodního stresu: kombinace Sentinel-2 + počasí + IoT se scikit-learn
- Praktická případová studie: vinice v Puglii, sezónní monitorování NDVI, optimalizace zavlažování
- Italská legislativa: AGEA, SIAN, open data CAP a digital CAP
FoodTech Series – všechny články
| # | Položka | Úroveň | Stát |
|---|---|---|---|
| 1 | IoT Pipeline pro přesné zemědělství s Pythonem a MQTT | Moderní | K dispozici |
| 2 | ML Edge pro monitorování plodin: Počítačové vidění v polích | Moderní | K dispozici |
| 3 | Satelitní a meteorologická rozhraní API pro AgriTech: Prediktivní data (jste zde) | Moderní | Proud |
| 4 | Sledovatelnost blockchainu v potravinách: od pole až po supermarket | Střední | Již brzy |
| 5 | Počítačová vize pro kontrolu kvality v potravinářském průmyslu | Moderní | Již brzy |
| 6 | FSMA a Digital Compliance: Automatizace regulačních procesů | Střední | Již brzy |
| 7 | Vertikální zemědělství: Kontrola životního prostředí s IoT a ML | Moderní | Již brzy |
| 8 | Prognóza poptávky pro maloobchod s potravinami s Prophet a LightGBM | Střední | Již brzy |
| 9 | Farm Intelligence Dashboard: Analýza v reálném čase s Grafanou | Střední | Již brzy |
| 10 | Optimalizace potravin v dodavatelském řetězci: ML pro snížení odpadu | Střední | Již brzy |
Program Copernicus: Evropská infrastruktura pro pozorování Země
Program Evropské unie Copernicus je největší pozorovací infrastrukturou v zemi Pozemek nikdy nezastavěný. Spravuje ESA (Evropská kosmická agentura) a Evropská komise s Provozní podpora EUMETSAT, Copernicus provozuje konstelaci pojmenovaných satelitů „Sentinel“, každý je určen pro konkrétní monitorovací mise. Pro zemědělství, absolutní referenční satelit e Sentinel-2.
Souhvězdí Sentinel-2 se skládá ze dvou dvojčat (Sentinel-2A a Sentinel-2B). sluneční synchronní oběžná dráha ve výšce 786 km. Konfigurace dvou satelitů zaručuje a čas 5denní recenze na rovníku a 2-3 dny v evropských zeměpisných šířkách (Itálie v ceně). Každý satelit nese na palubě MSI (MultiSpectral Instrument) senzor, násadu skener, který získává data o šířce swipe 290 km a 13 spektrálních pásmech distribuovány ve třech prostorových rozlišeních: 10 metrů, 20 metrů a 60 metrů.
Nejstrategičtější vlastností Sentinelu-2 je jeho zcela zdarma: od roku 2014 všechna data Copernicus jsou k dispozici v otevřeném přístupu pro jakékoli použití, komerční i nekomerční komerční. Nový přístupový portál z roku 2023 a Ekosystém datového prostoru Copernicus (CDSE), který nahradil předchozí Copernicus Open Access Hub (SciHub) vyřazený z provozu v r 2023. CDSE nabízí okamžitý přístup ke kompletnímu archivu Sentinel-2 s latencí akvizice pro nejnovější snímky jsou k dispozici přibližně 3 hodiny.
Satellite Sentinel-2: Technické specifikace snímače MSI
| Kapela | Vlnová délka (nm) | Rezoluce | Hlavní aplikace |
|---|---|---|---|
| B1 - Pobřežní aerosoly | 443 | 60 m | kvalita vzduchu, korekce atmosféry |
| B2 - Modrá | 490 | 10 m | Diskriminace vegetace/půdy, mapování |
| B3 - Zelená | 560 | 10 m | Vegetační vitalita, barva listů |
| B4 -Červený | 665 | 10 m | Chlorofyl, NDVI (červený pruh) |
| B5 - Vegetace Red Edge | 705 | 20 m | Stres vegetace, Red Edge NDVI |
| B6 - Vegetace Red Edge | 740 | 20 m | Obsah chlorofylu, druhová klasifikace |
| B7 - Vegetace Red Edge | 783 | 20 m | Biomasa, LAI |
| B8 - NIR | 842 | 10 m | NDVI (pásmo NIR), biomasa, LAI |
| B8a - Vegetace Red Edge | 865 | 20 m | Struktura vrchlíku, NDVI úzký |
| B9 - Vodní pára | 940 | 60 m | Korekce atmosférické vodní páry |
| B10 - SWIR - Cirrus | 1375 | 60 m | Detekce cirrusových mraků |
| B11 - ZAVÍREJTE | 1610 | 20 m | Vodní stres, obsah vody v listech |
| B12 - ZAVÍREJTE | 2190 | 20 m | Pokročilý vodní stres, stárnutí |
Pro zemědělství jsou relevantní především dvě úrovně zpracování dat Sentinel-2: Úroveň 1C (L1C) - nejvyšší vyzařování atmosféry s geometrickou korekcí, např Úroveň 2A (L2A) - odrazivost na dně atmosféry (Bottom of Atmosphere, BOA) s atmosférickou korekcí aplikovanou pomocí algoritmu Sen2Cor. Pro přímé zemědělské aplikace, Úroveň 2A je doporučená, protože je již korigována na vliv atmosféry srovnatelné hodnoty odrazivosti mezi různými daty a různými geografickými polohami.
Satelitní platformy pro AgriTech: Úplné srovnání 2025
Sentinel-2 není jedinou dostupnou možností. Trh satelitních dat pro zemědělství zahrnuje komerční poskytovatele, jako je Planet Labs, zavedené společnosti jako Landsat (NASA/USGS) a integrované cloudové platformy, jako je Google Earth Engine. Výběr závisí na složitém kompromisu mezi prostorovým rozlišením, frekvencí akvizic, spravovaným cloudovým pokrytím a dostupným rozpočtem.
Srovnání satelitních platforem pro přesné zemědělství – 2025
| Platforma | Prostorové rozlišení | Znovu navštívit frekvenci | Spektrální pásma | Náklady | Latence dat | API |
|---|---|---|---|---|---|---|
| Sentinel-2 (ESA/Copernicus) | 10 m (vis/NIR), 20 m (SWIR) | 5 dní (2-3 v EU) | 13 pásem MSI | Zdarma (otevřená data) | ~3 hodiny od pořízení | CDSE OData, STAC, SentinelHub, OpenEO |
| Planeta PlanetScope | 3,7 m | Denně (kvazidenní) | 8 pásem (PS2.SD) | Komerční předplatné; zdarma pro výzkum | 24 hodin | Planet API v1, Subscriptions API |
| Landsat 8/9 (NASA/USGS) | 30 m (multispektrální), 15 m (panev) | 16 dní | 11 pásem OLI/TIRS | Zdarma (otevřená data) | 24-48 hodin | USGS EarthExplorer, Google Earth Engine |
| MODIS (NASA Earth/Aqua) | 250 m / 500 m / 1 km | 1-2 dny | 36 kapel | Zdarma (otevřená data) | 24-48 hodin | NASA EarthData STAC, Google Earth Engine |
| Google Earth Engine | Závisí na datovém souboru (10 m-1 km) | Záleží na datové sadě | Integrované multi-datové sady | Zdarma nekomerční; placená reklama | Okamžité cloudové zpracování | Python ee, JavaScript API |
| Maxar WorldView-3 | 0,3 m (panev), 1,24 m (multi) | 1-4 dny | 29 kapel (CAVIS) | ~25-40 USD/km2 | 4-8 hodin | Maxar Streaming API |
| Airbus Plejády Neo | 0,3 m | 1-2 dny | 6 multispektrálních pásem | ~15-30 EUR/km2 | 24-48 hodin | OneAtlas API |
Pro velkou většinu zemědělských aplikací v Itálii, Sentinel-2 a volba optimální: rozlišení 10 metrů je dostatečné pro pozemky větší než 0,5 hektaru (prahová hodnota, pod kterou se statistika pixelů stává nespolehlivou), frekvenci opakovaných návštěv 2-3 dny v Itálii a dostačující pro sezónní sledování a zcela zdarma eliminuje jakoukoli ekonomickou překážku přijetí. Planet Labs se stává relevantní pouze pro aplikace, které vyžadují téměř každodenní monitorování (například detekce raného nástupu vodní stres u plodin s vysokou hodnotou, jako je intenzivní ovoce a zelenina) nebo rozlišení pod 10 m. MODIS a Landsat zůstávají užitečné pro analýzy velkých oblastí a vícedekádových časových řad.
Vegetační indexy: NDVI, EVI, SAVI a LAI pro italské plodiny
Vegetační indexy (VI) jsou metriky odvozené matematicky kombinace spektrálních pásů, které kvantifikují specifické vlastnosti vegetace. Vykořisťují skutečnost, že chlorofyl silně absorbuje v červeném pásu (665 nm) a silně odráží v blízké infračervené oblasti (842 nm): poměr mezi těmito pásy kóduje hustotu a vitalitu vegetace s matematickou elegancí.
NDVI - Normalized Difference Vegetation Index
NDVI je nejpoužívanější vegetační index na světě, který zavedli Rouse et al. v roce 1973. Matematický vzorec je:
Vzorec NDVI
NDVI = (NIR - Red) / (NIR + Red)
Su Sentinel-2:
NDVI = (B8 - B4) / (B8 + B4)
Range: da -1.0 a +1.0
Hodnoty NDVI jsou interpretovány podle standardizované stupnice, ale optimální prahové hodnoty se liší podle plodinové a fenologické fáze. V následující tabulce jsou uvedeny referenční prahové hodnoty pro hlavní italské plodiny, odvozené z vědecké literatury a empirického ověření v typických oblastech (Pianura Padana, Puglia, Sicílie, Toskánsko):
Interpretace NDVI pro italské plodiny
| Rozsah NDVI | Obecný výklad | Pšenice (Triticum) | Réva (Vitis) | Rajče (Solanum) | kukuřice (Zea mays) | Olivový (Olea) |
|---|---|---|---|---|---|---|
| < 0,1 | Holá půda, kámen, sníh | Neosetá půda | Zima, předpučení | Posklizňová | Předsetí | Holá půda mezi řadami |
| 0,1 - 0,2 | Velmi řídký nebo suchý porost | Počáteční pohotovost | Dormance/silný stres | Nedávná transplantace | Nouzové (VE) | Silný vodní/tepelný stres |
| 0,2 - 0,4 | Řídký porost, střední stres | Počáteční odkopávání | Předkvět, lehký stres | Počáteční vegetativní růst | Fáze V1-V3 | Špatná vegetace, stres |
| 0,4 - 0,6 | Mírná vegetace, dobrý zdravotní stav | Plné stoupání | Předkvětová / ovocná sada | Aktivní vegetativní vývoj | Etapy V4-V8 | Normální listy, dobře zalévané |
| 0,6 - 0,8 | Hustá vegetace, vynikající vitalita | Okruh (sezónní špička) | Plná foliace, veraison | Maximální krytí (kvetoucí) | Fáze V10-VT (kvetoucí) | Husté olistění, výborný stav |
| > 0,8 | Velmi hustý porost, les | Vzácné (lesy na okraji pole) | Živé ploty a hraniční stromy | Skleníky s celkovým pokrytím | Vzácné optimální podmínky | Olivové háje s vysokou hustotou |
EVI - Enhanced Vegetation Index
EVI, vyvinutý Huete et al. pro snímač MODIS koriguje limity NDVI v oblastech s vysokou hustotou vegetace (saturace NDVI nad 0,7) a v oblastech s obnaženými půdami (chyby NDVI způsobené odrazem země). Vzorec je:
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 - Soil Adjusted Vegetation Index
SAVI, navržený Huete v roce 1988, je zvláště užitečný v semiaridních prostředích, jako je např. Puglia nebo Sicílie, kde exponovaná půda významně ovlivňuje spektrální odezvu. Korekční faktor L snižuje vliv země:
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 ItaliaLAI - Index listové plochy
LAI (Leaf Area Index) kvantifikuje celkový povrch listu na jednotku plochy povrchu půda. Je to základní agronomický parametr související s potenciální produktivitou a pocení. Sentinel-2 vám umožňuje odhadnout LAI pomocí empirických nebo fyzikálních algoritmů (SNAP BioPhysical Processor):
# 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 massimaOmezení indexů satelitní vegetace
- oblačnost: Sentinel-2 a optické - mraky dělají obraz nepoužitelným. V Itálii je průměrná oblačnost 48 % v zimě a 15 % v létě. Pásmo QA60 umožňuje automatické maskování oblačnosti.
- Saturace NDVI: Nad NDVI ~0,7-0,8 je citlivost drasticky snížena. Pro plodiny s vysokou hustotou použijte EVI.
- Mix pixelů: Při rozlišení 10 m může pixel obsahovat vegetaci i půdu. U malých pozemků (< 0,5 ha) je průměrná NDVI méně reprezentativní.
- Fenologická sezónnost: Porovnání NDVI mezi různými daty vyžaduje vzít v úvahu fenologickou fázi, nikoli pouze absolutní hodnotu.
- Meziroční variabilita: NDVI se také liší v závislosti na sezónních povětrnostních podmínkách (teplota, srážky) bez ohledu na agronomické postupy.
Implementace Pythonu: CDSE Access and Download Sentinel-2
Copernicus Data Space Ecosystem (CDSE) se stal jediným portálem pro přístup k datům Sentinel od roku 2023. Nabízí čtyři hlavní rozhraní API: OData API (vyhledávání a stahování produktů), STAC API (SpatioTemporal Asset Catalog), OpenEO API (zpracování standardizované) e Sentinel Hub API (zpracování jako služba). Chcete-li začít, vytvořte si bezplatný účet na dataspace.copernicus.eu a vygenerujte přihlašovací údaje OAuth2.
Nastavení prostředí Python a závislosti
# 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 idricoAutentizace CDSE s OAuth2
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)Vyhledávejte produkty Sentinel-2 přes OData API
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]}")Stahujte pásma a přistupujte k nim pomocí Sentinel Hub Python
Pro zpracování pásma bez stahování celé scény (která váží 800 MB - 1,2 GB na L2A produkt), je efektivnější používat Sentinel Hub Process API, které vrací pouze pásma požadovaná v oblasti zájmu:
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}%)")Výpočet NDVI a ukládání GeoTIFF
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")Počasí API pro AgriTech: Porovnání a implementace
Údaje o počasí jsou druhým kritickým zdrojem dat precizního zemědělství. Díky integraci s daty NDVI a IoT v terénu vám umožňují vytvářet prediktivní modely vodní stres, riziko onemocnění, evapotranspirace a předpověď výnosů. Panorama meteorologických API v roce 2025 se dělí na bezplatné open-source poskytovatele, komerční poskytovatele s bezplatnými úrovněmi a regionálními službami italské ARPA.
Počasí API pro AgriTech – srovnání 2025
| Poskytovatelé | Cena | Prostorové rozlišení | Předpověď | Agro proměnné | Historický archiv | API |
|---|---|---|---|---|---|---|
| Otevřené počasí | Zdarma (nekomerční); Reklama ~ 15 $ měsíčně | 1-11 km (více modelů) | 16 dní | ET0, teplota půdy, vlhkost půdy, tlak par | 1940–současnost (ERA5) | REST JSON, klient Pythonu |
| OpenWeatherMap | Zdarma až 60 hovorů/min; 40 $/měsíc pro | 2,5 km (ICON) | 5 dní (zdarma), 16 dní (pro) | Omezené, bez ET0 | Pouze placená historie | REST JSON, Python SDK |
| Zítra.io | Zdarma 500 hovorů/den; 199 $ měsíčně jádro | 1 km | 21 dní | SoilMoisture, ET0, podmínky postřiku, riziko škůdců | 6 let | REST, WebSocket, gRPC |
| ARPA Lombardie / Piemont / Veneto | Zdarma (regionální otevřená data) | Přesné stanice (síť ~1500 stanic v Itálii) | Ne (pouze pozorováno) | Všechny agro-počasí proměnné | Kompletní archiv (desítky let) | Stažení WMS, REST, CSV/JSON |
| Meteoam / CFS / ECMWF | ECMWF zdarma otevřená data; Zdarma NOAA CFS | 9-25 km | 10-45 dní | Proměnné založené na počasí | 1940-dosud prostřednictvím ERA5 | ECMWF API (Python), CDS API |
| Vizuální přejezd | Zdarma 1000 záznamů/den; Standardně 35 $ měsíčně | Interpolováno ze stanic | 15 dní | ET0, dny růstu, riziko mrazu | Neomezené placené | REST JSON, CSV, jako SQL |
Pro aplikace AgriTech ve výrobě se doporučuje použít Otevřené počasí jak základna (volná, open-source, data ERA5 pro historiky, 16 dní předpovědi), která ji integruje s sítě stanic ARPA pro místní kalibraci. Open-Meteo zahrnuje agrometeorologické proměnné kritické problémy, jako je ET0 (referenční evapotranspirace), vlhkost půdy ve více hloubkách a pára vodné, které nejsou vždy dostupné zdarma od komerčních poskytovatelů.
Implementace klienta Open-Meteo s mezipamětí a opakováním
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")Kompletní Geospatial Pipeline: GeoPandas, rasterio a PostGIS
Profesionální správa satelitních dat vyžaduje komplexní geoprostorovou infrastrukturu která integruje rastrová data (satelitní snímky NDVI, mapy počasí) s vektorovými daty (hranice pozemky, hospodářské oblasti, místa odběru vzorků). Potrubí, které zde popisujeme, přijímá de facto standardní zásobník Pythonu pro geoprostor: rasterio pro rastry, GeoPandas pro vektory, PostGIS jako prostorová databáze e GDAL jako základní stroj pro překlad formátů.
Architektura geoprostorového potrubí
AgriTech Geospatial Pipeline Technology Stack
┌─────────────────────────────────────────────────────────────────────────┐
│ DATI DI INPUT │
│ [Sentinel-2 GeoTIFF] [Shapefile Appezzamenti] [CSV Meteo/IoT] │
│ │ │ │ │
│ rasterio geopandas pandas │
└─────────┼──────────────────────┼───────────────────────┼────────────────┘
│ │ │
┌─────────▼──────────────────────▼───────────────────────▼────────────────┐
│ PROCESSING LAYER (Python) │
│ [NDVI Calculation] [Zonal Statistics] [Spatial Join] │
│ numpy/rasterio rasterstats geopandas │
│ │ │ │ │
└─────────┼──────────────────────┼───────────────────────┼────────────────┘
│ │ │
└──────────────────────▼───────────────────────┘
│
┌────────────────────────────────▼────────────────────────────────────────┐
│ STORAGE LAYER (PostGIS + S3) │
│ [PostGIS - geometrie + attributi] [S3/MinIO - GeoTIFF raster] │
│ - Tabella parcelle con NDVI medio - File raster originali │
│ - Serie storica per parcella - Archivio scene satellite │
│ - Spatial index GIST - Formato Cloud Optimized GeoTIFF │
└────────────────────────────────┬────────────────────────────────────────┘
│
┌────────────────────────────────▼────────────────────────────────────────┐
│ SERVING LAYER │
│ [REST API - FastAPI] [Dashboard - Grafana] [ML Models] │
│ NDVI per appezzamento Mappe NDVI tempo reale Predizione stress │
└─────────────────────────────────────────────────────────────────────────┘
Zonální statistika: Průměrné NDVI na pozemek s rastrovými statistikami
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)")PostGIS dotazy pro pokročilou prostorovou analýzu
-- 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;Prediktivní model vodního stresu s ML
Kombinace satelitního NDVI s meteorologickými a IoT daty vám umožňuje vytvářet prediktivní modely vodního stresu 3-7 dní předem, což umožňuje proaktivní plánování zavlažování. Zde popsaný model integruje tři zdroje dat: časovou řadu NDVI ze Sentinel-2, vodní bilance vypočítaná z ET0 a srážek a půdní vlhkost ze senzorů IoT (nebo z Open-Meteo ERA5-Land, pokud nejsou k dispozici fyzické senzory).
Funkce inženýrství pro model
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"]]Případová studie: Monitorování NDVI ve vinici Primitivo v Mandurii
Aby byly popsané pojmy konkrétní, uvádíme případovou studii inspirovanou aplikací skutečné na území DOC Primitivo di Manduria (Taranto, Puglia). Oblast má podmínky typické pro středomořské vinařství: horká a suchá léta (letní srážky < 50 mm), jílovito-vápencové půdy, tréninkové systémy protiespalier a stromků.
Parametry případové studie - Vigna Primitivo Manduria
| Parametr | Hodnota |
|---|---|
| Umístění | Manduria (TA), Puglia - šířka 40,38, délka 17,47 |
| Celková plocha povrchu | 35 hektarů rozdělených na 12 parcel |
| Hlavní odrůdy | Primitivo (80 %), Negroamaro (20 %) |
| Systém chovu | Apulský stromek (8 kusů), Counter-espalier (4 kusy) |
| Rok rostliny | 2003-2015 (smíšené mladé/dospělé vinice) |
| Zavlažování | Po kapkách (ne všechny pozemky) |
| Sledovací satelit | Sentinel-2A/B, L2A, rozlišení 10m |
| Frekvence používaných akvizic | 28 scén na venkově 2024 (oblačno < 20 %) |
| Integrované IoT senzory | 6 meteostanic, 18 senzorů vlhkosti půdy na 30/60 cm |
| Analytické období | duben–říjen 2024 |
Sezónní výsledky NDVI a rozhodnutí o zavlažování
Sezónní analýza NDVI odhalila významně diferencované vzorce vodního stresu mezi pozemky, s přímými důsledky pro řízení zavlažování a kvalitu produktu:
Výsledky NDVI podle grafu – kampaň 2024
| Spiknutí | Odrůda | NDVI duben | NDVI dolů | NDVI Jul (vrchol) | NDVI srpen | sada NDVI | Upozornění na zavlažování |
|---|---|---|---|---|---|---|---|
| Park A1 | Primitivní/stromek | 0,28 | 0,42 | 0,58 | 0,41 | 0,31 | Mírné namáhání jehly |
| Park A2 | Primitivní/stromek | 0,22 | 0,36 | 0,47 | 0,29 | 0,24 | Silné namáhání jehlou |
| Park B1 | Primitivní/kontrapalér | 0,31 | 0,53 | 0,66 | 0,57 | 0,44 | Žádná upozornění |
| Park B2 | Primitivní/kontrapalér | 0,29 | 0,49 | 0,63 | 0,52 | 0,40 | Žádná upozornění |
| Park C1 | Negroamaro/strom | 0,25 | 0,44 | 0,55 | 0,38 | 0,28 | Mírné namáhání jehly |
| Park C2 | Negroamaro/counterspalliera | 0,30 | 0,51 | 0,61 | 0,50 | 0,38 | Žádná upozornění |
Vzor, který se objevil, je jasný a agronomicky významný: stromkové pozemky bez nouzové zavlažování vykazují výraznější vodní stres v srpnu až září ve srovnání s pult-espalier parcely s odkapávacím systémem. Pozemek A2 má zejména vykázal pokles NDVI z 0,47 na 0,29 mezi červencem a srpnem 2024, dočasně koreloval s sled 23 dnů bez srážek a maximální teploty nad 36 stupňů.
Skript automatického upozornění
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['varieta']}, {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")ROI satelitního monitorovacího systému
Kvantifikace návratnosti investic satelitního sledování je zásadní pro odůvodnit technologické investice do vinařských společností. Pro vinici Manduria, analýza nákladů a přínosů pro kampaň 2024 ukazuje:
Analýza návratnosti investic – satelitní monitorovací systém Vigna Manduria (35 ha, kampaň 2024)
| Hlas | Hodnota | Poznámky |
|---|---|---|
| NÁKLADY | ||
| Sentinel-2 Data Access (CDSE) | 0 EUR | Otevřená data zdarma |
| Open-Weather Weather API | 0 EUR | Nekomerční zdarma |
| Vývoj Python pipeline | 2 500 EUR | 30 vývojářských hodin x 83 EUR/h (jednorázově, amortizováno po dobu 3 let) |
| Cloud hosting (VM + úložiště) | 600 EUR/rok | VM 4 vCPU + 50 GB úložiště S3 |
| Integrační agronomické poradenství | 1 000 EUR | Jednorázová kalibrace prahu pro konkrétní plodinu |
| Celkové náklady za rok 1 | 4 100 EUR | Následující roky: ~1 400 EUR/rok |
| ODHADOVANÉ PŘÍNOSY (35 ha, 100 000 lahví/rok) | ||
| Snížení zavlažování (optimalizace načasování) | 2 800 EUR | 15% úspora vody, 8 500 m3 x 0,33 EUR/m3 |
| Omezení fytosanitárních ošetření | 1 750 EUR | Cílená intervence pouze v rizikových oblastech (-20 % ošetření) |
| Zlepšení kvality hroznů | 8 500 EUR | +0,5 bodu průměrný stupeň Brix na 12 ha, +0,20 EUR/kg prémie za jakost |
| Snížení ztrát z vodního stresu | 3 200 EUR | 8% výtěžnost na 2 plochách při očekávaném silném stresu |
| Celkové výhody rok 1 | 16 250 EUR | |
| Rok 1 ROI | 296 % | Návratnost cca 3 měsíce |
Google Earth Engine: Cloudové zpracování pro analýzu ve velkém měřítku
Pro analýzu v regionálním nebo celostátním měřítku (např. mapování odrůd na úrovni Okres DOC, hodnocení škod jarními mrazy v zemském měřítku, monitoring CAP pro platební agentury, jako je AGEA), Google Earth Engine (GEE) nabízí výpočetní kapacitu cloud-native, které by jinak vyžadovaly vyhrazené infrastruktury HPC.
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}")Předpisy a otevřená data pro AgriTech v Itálii
Zemědělský ekosystém otevřených dat v Itálii je strukturován kolem tří hlavních aktérů: AGEA (Zemědělská platební agentura), SIAN (Národní zemědělský informační systém) a navazující regionální portály. V únoru 2025 AGEA dokončila migraci SIAN směrem k Národnímu strategickému pólu, aktualizaci cloudovou infrastrukturu a rozšiřující možnosti zpracování geoprostorových dat.
Hlavní italské zemědělské otevřené zdroje dat
| Zdroj | URL | Data jsou k dispozici | Formát | Aktualizovat |
|---|---|---|---|---|
| Otevřená data SIAN | data.sian.it | Ortofota s vysokým rozlišením, mapa využití půdy (na základě AI), firemní grafické plány | GeoTIFF, Shapefile, WFS | Výroční |
| Služby AGEA EO | sian.it | Multispektrální ortofota, X-pásmové radarové snímky, DEM, DSM | GeoTIFF, ECW | Periodikum (kampaň) |
| ISTAT - Zemědělská statistika | istat.it/agricoltura | UAA, výroba, struktura společnosti, zemědělský census | CSV, JSON, SDMX | Výroční |
| Národní geoportál MATTM | geoportale.gov.it | Carta Natura, Corine Land Cover, IDT, DBT | WMS, WFS, Shapefile | Variabilní |
| Regionální ARPA (15 regionů) | harfa.[region].it | Meteoklimatické stanice, regionální předpovědní modely | CSV, JSON, NetCDF | Téměř v reálném čase |
| Pozemní služba Copernicus | land.copernicus.eu | CORINE Land Cover, travní porosty, voda a vlhko, atlas měst | GeoTIFF, Shapefile | Trienále/ročník |
Přechod daňových kreditů PNRR 5.0 pro AgriTech
Investice do satelitních monitorovacích systémů a precizního zemědělství mohou mít nárok na výhody plánu přechodu 5.0 (zákon 207/2024 – zákon o rozpočtu 2025). Zemědělské podniky, které budou přijaty od roku 2025, mají přístup k daňovým úlevám na investice v hmotném a nehmotném investičním majetku 4.0, včetně softwaru pro řízení podniku s komponenty umělé inteligence a analýza geoprostorových dat.
Přechodové sazby 5.0 daňových kreditů pro zemědělské podniky (2025)
| Investiční rozsah | Úspora energie 3–6 % | Úspora 6–10 % | Úspora > 10 % |
|---|---|---|---|
| Až 2,5 milionu EUR | 35 % | 40 % | 45 % |
| Od 2,5 mil. do 10 mil. EUR | 15 % | 20 % | 25 % |
| Od 10 milionů do 50 milionů EUR | 5% | 10 % | 15 % |
Mohou se kvalifikovat přesné zavlažovací systémy řízené NDVI a předpověďmi počasí pro pásmo „úspora energie >6 %“, pokud doloží snížení spotřeby vody (čerpací energie) ve srovnání se základní linií. Žádost o přísežnou technickou expertizu.
Produkční architektura: Best Practices a Anti-Patterns
Po popisu komponent je důležité shrnout osvědčené postupy pro robustní implementace ve výrobě a identifikovat společné anti-vzory, které vést ke křehkým nebo nepřesným systémům.
Osvědčené postupy – Satelitní monitorovací systém AgriTech
- Progresivní cloudové filtrování: Nepoužívejte pouze procento cloudového pokrytí produktu (metadata úrovně 1). Vždy aplikujte maskování na úrovni pixelů pomocí QA60 (Sentinel-2) nebo BQA (Landsat). 5% cloudový produkt může mít konkrétní graf zcela zakrytý.
- GeoTIFF optimalizovaný pro cloud (COG): Ukládejte rastry NDVI ve formátu COG pro efektivní přístup k požadavku na rozsah z S3/MinIO. Vyhněte se klasickým GeoTIFFům, které nejsou optimalizovány pro cloudové úložiště.
- Konzistentní projekce: Standardizujte VŠECHNO na EPSG:4326 (WGS84) pro konečná data nebo na UTM zónu 32N (EPSG:32632) pro Itálii, aby byla zachována správná metrická měření. Nikdy nemíchejte CRS bez explicitní reprojekce.
- Ověření dat pomocí základní pravdy IoT: Kalibrujte prahové hodnoty NDVI pomocí měření v terénu (senzory vlhkosti půdy, přenosný měřič LAI, shromážděná agronomická data). Samotné NDVI bez místní kalibrace může poskytnout falešně pozitivní/negativní výstrahy o 15–30 %.
- Správa časového archivu: Udržujte alespoň 3-5 let časové řady NDVI pro každý pozemek. Anomálie NDVI ve srovnání s historickým průměrem (stejný měsíc, stejná plodina) je mnohem vypovídající než absolutní hodnota.
- Rozhraní API pro omezení rychlosti a ukládání do mezipaměti: CDSE a Open-Meteo mají rychlostní limity. Vždy implementujte místní mezipaměť (založenou na souborech nebo Redis), abyste se vyhnuli zbytečnému opětovnému stahování. Stažení Sentinel-2 trvá 30–120 sekund: nepožadujte stejná data dvakrát.
- Protokolování a auditní záznam: Každý výpočet NDVI musí být sledovatelný ke zdrojové satelitní scéně, datu akvizice, procentu oblačnosti, verzi algoritmu. Zásadní pro audity PAC/AGEA a pro ladění anomálií.
Společné anti-vzory – satelitní systémy AgriTech
- Ignorování atmosférické korekce: Použití úrovně 1C (horní část atmosféry) namísto úrovně 2A (spodní část atmosféry) zavádí systematické chyby v NDVI řádově 10–20 % v podmínkách vysoké vzdušné vlhkosti (mlha v údolí Pádu, Jaderské moře). VŽDY používejte L2A pro kvantitativní aplikace.
- NDVI jako jediný indikátor: Spolehněte se pouze na NDVI a ignorujte EVI (pro hustou vegetaci), NDWI (vodní stres), Red Edge NDVI (brzký stres, než je viditelný na NDVI). Profesionální systém používá minimálně 3-4 indexy v kombinaci.
- Velikost pixelu a balíček: Použití NDVI na pozemky menší než 0,3 ha s pixely ve vzdálenosti 10 m vytváří statisticky nespolehlivé průměry (méně než 30 pixelů). Pro malé pozemky zvažte Planet Labs (3,7 m) nebo senzory dronů.
- Nedostatek časového skládání: Vezmeme-li jedinou dostupnou akvizici na určité datum, dokonce i při 19% oblačnosti, vede k artefaktům NDVI v dílčích oblastech oblačnosti. Vždy používejte medoid compositing v časovém okně 10-15 dnů.
- Potrubí bez verzování dat: Přepočet NDVI s různými verzemi Pythonu/rasterio/sentinelhub vytváří mírně odlišné hodnoty pro stejná zdrojová data. Verze vašich výpočetních kanálů pomocí DVC nebo ekvivalentu.
Complete Pipeline: Orchestrace s prefektem
Díky integraci všech popsaných komponent je kompletní satelitní sledovací potrubí lze zorganizovat pomocí Prefect (viz série Pipeline Orchestrace) spouštět automaticky každých 5 dní (s frekvencí opakovaných návštěv Sentinel-2).
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()Závěry a další kroky
Představují družicová data Sentinel-2, dostupná zdarma přes CDSE pravděpodobně nejvíce málo využívaný zdroj dat v italském AgriTech. S Pythonem, rasterio, sentinelhub a účet CDSE, je možné sestavit a Monitorovací systém NDVI, který kvalitou a granularitou předčí mnohá řešení komerční za poplatek.
Klíčem k úspěchu není technologická vyspělost, ale kalibrace místní agronomické: prahové hodnoty NDVI, časová okna výstrah, váhy prvky v prediktivním modelu musí být kalibrovány pro konkrétní plodinu, na místní klima a agronomické postupy společnosti. Kalibrovaný systém na vinné révě v Puglii to bez něj stejně dobře nefunguje ani na pšenici v údolí Pádu lidský kalibrační zásah.
Klíčové body, které je třeba pamatovat pro ty, kteří chtějí implementovat podobný systém:
- Vždy používejte Sentinel-2 Úroveň-2A (aplikovaná atmosférická korekce) pro kvantitativní srovnání mezi různými daty.
- Kombinujte NDVI s alespoň jedním doplňkovým indexem (EVI pro husté vinice, NDWI pro raný vodní stres, Red Edge NDVI pro monitorování chlorofylu).
- Integrovat Open-Weather ERA5-Land pro vodní bilanci: ET0 + srážky + vlhkost půdy z modelu a vynikající proxy, když nemáte IoT senzory v terénu.
- Nechte si jeden historická řada minimálně 3 roky pro každý pozemek: anomálie ve srovnání s historickým průměrem je mnohem užitečnější než absolutní hodnota NDVI.
- Měřítko s Google Earth Engine pouze tehdy, když potřebujete zpracovat regionální nebo národní měřítka: pro jednotlivé společnosti je místní Python pipeline na CDSE lépe ovladatelný a bez lock-in.
FoodTech Series pokračuje
Postavili jste satelitní potrubí. Nyní prozkoumejte další články ze série pro dokončete architekturu AgriTech:
- článek 1: IoT Pipeline pro přesné zemědělství s Pythonem a MQTT - Jak integrovat pozemní senzory se satelitním potrubím.
- Článek 2: ML Edge pro monitorování plodin: Počítačové vidění v polích - Modely počítačového vidění pro včasnou detekci chorob rostlin.
- Článek 4: Sledovatelnost blockchainu v potravinách: od pole až po supermarket - Jak NDVI data pohání certifikáty kvality v řetězci.
Chcete-li se dozvědět více o MLOps a nasazení prediktivních modelů popsaných v tomto článku, viz seriál MLOps for Business s MLflow a seriál LLM v podnikání: RAG Enterprise.