Interfejsy API satelitów i pogody dla AgriTech: dane predykcyjne za pomocą Sentinel-2 i Python
W 2024 roku Europejska Agencja Kosmiczna przejęła ponad 1,5 petabajta danych każdego dnia z satelitów Sentinel programu Copernicus. Spośród nich część poświęcona monitorowaniu rolnictwa w Europie stanowi największy i najbardziej strategicznie istotny segment. Jednak zdecydowana większość włoskich firm rolniczych, w tym tych najbardziej zaawansowanych technologicznie, jeszcze nie eksploatuje tę publiczną, bezpłatną i potwierdzoną naukowo infrastrukturę.
Powodem nie jest brak danych, lecz złożoność ich dostępu, przetwarzania i przekształcania konkretne decyzje agronomiczne. Satelita Sentinel-2 przelatuje nad Twoją winnicą co 5 dni rozdzielczość 10 metrów na piksel. Każde przejście wytwarza 13 pasm widmowych, które kodują zdrowie roślinności, stres wodny, zawartość chlorofilu, obecność pasożytów. W połączeniu z danymi pogodę, czujniki naziemne IoT i modele predykcyjne – dane te mogą radykalnie przekształcić jakość decyzji agronomicznych, zmniejszając koszty środków produkcji o 15-25% i zwiększając plony do 10-20%.
Globalny rynek obserwacji Ziemi jest tego wart 10,07 miliardów dolarów w 2025 roku według Straits Research, przy oczekiwanym wzroście do 17,20 miliardów do 2033 r. (CAGR 6,92%). Segment rolnictwa reprezentuje 21% wniosków i rozwija się w najszybszym tempie, napędzany popytem na rolnictwo precyzyjne, monitorowanie plonów i optymalizację nawadniania. W tym artykule budujemy krok po kroku kompletny potok Pythona umożliwiający dostęp do Sentinel-2 poprzez CDSE (Ekosystem przestrzeni danych Copernicus), obliczyć NDVI i inne wskaźniki roślinności, zintegrować dane meteorologiczne z Open-Meteo i posłużą do stworzenia modelu predykcyjnego stresu wodnego.
Czego dowiesz się w tym artykule
- Architektura programu Copernicus i darmowy dostęp do Sentinel-2 poprzez CDSE
- Wskaźniki wegetacyjne: NDVI, EVI, SAVI, LAI - wzory matematyczne, interpretacja i progi dla upraw włoskich
- Kompletna implementacja Pythona: uwierzytelnianie CDSE, pobieranie pasma, obliczanie NDVI za pomocą rasterio i numpy
- Porównanie interfejsów API pogody: Open-Meteo, OpenWeatherMap, Tomorrow.io i integracja danych satelitarnych
- Potok geoprzestrzenny z GeoPandas, rasterio i PostGIS dla danych wektorowych i rastrowych
- Model predykcyjny stresu wodnego: połączenie Sentinel-2 + pogody + IoT z scikit-learn
- Praktyczne studium przypadku: winnica w Apulii, sezonowy monitoring NDVI, optymalizacja nawadniania
- Ustawodawstwo włoskie: AGEA, SIAN, WPR otwartych danych i WPR cyfrowa
Seria FoodTech - Wszystkie artykuły
| # | Przedmiot | Poziom | Państwo |
|---|---|---|---|
| 1 | Rurociąg IoT dla rolnictwa precyzyjnego z Pythonem i MQTT | Zaawansowany | Dostępny |
| 2 | ML Edge do monitorowania upraw: wizja komputerowa na polach | Zaawansowany | Dostępny |
| 3 | Interfejsy API satelitów i pogody dla AgriTech: dane predykcyjne (tutaj jesteś) | Zaawansowany | Aktualny |
| 4 | Identyfikowalność Blockchain w żywności: od pola do supermarketu | Mediator | Już wkrótce |
| 5 | Wizja komputerowa w kontroli jakości w przemyśle spożywczym | Zaawansowany | Już wkrótce |
| 6 | FSMA i zgodność cyfrowa: automatyzacja procesów regulacyjnych | Mediator | Już wkrótce |
| 7 | Rolnictwo pionowe: kontrola środowiska za pomocą IoT i ML | Zaawansowany | Już wkrótce |
| 8 | Prognozowanie popytu na sprzedaż detaliczną żywności za pomocą Prophet i LightGBM | Mediator | Już wkrótce |
| 9 | Pulpit nawigacyjny Farm Intelligence: analityka w czasie rzeczywistym za pomocą Grafany | Mediator | Już wkrótce |
| 10 | Optymalizacja żywności w łańcuchu dostaw: ML na rzecz redukcji odpadów | Mediator | Już wkrótce |
Program Copernicus: Europejska infrastruktura obserwacji Ziemi
Unijny program Copernicus to największa infrastruktura obserwacyjna na świecie Ziemia nigdy nie zbudowana. Zarządzane przez ESA (Europejską Agencję Kosmiczną) i Komisję Europejską Wsparcie operacyjne EUMETSAT, Copernicus obsługuje konstelację nazwanych satelitów „Sentinel”, każdy przeznaczony do konkretnych misji monitorujących. Dla rolnictwa, absolutny satelita referencyjny, np Strażnik-2.
Konstelacja Sentinel-2 składa się z dwóch bliźniaczych satelitów (Sentinel-2A i Sentinel-2B) orbita synchroniczna ze Słońcem na wysokości 786 km. Konfiguracja dwóch satelitów gwarantuje: czas Przegląd 5-dniowy na równiku i 2-3 dni na europejskich szerokościach geograficznych (Włochy zawarte). Każdy satelita ma na pokładzie czujnik MSI (MultiSpectral Instrument), czyli miotłę skaner zbierający dane z szerokością przesuwania 290 km i 13 pasmami widmowymi dystrybuowany w trzech rozdzielczościach przestrzennych: 10 metrów, 20 metrów i 60 metrów.
Najbardziej strategiczną cechą Sentinel-2 jest jego całkowicie za darmo: od 2014 r wszystkie dane programu Copernicus są dostępne w otwartym dostępie do dowolnego wykorzystania, komercyjnego lub niekomercyjnego komercyjny. Nowy portal dostępowy od 2023 roku i Ekosystem przestrzeni danych programu Copernicus (CDSE), który zastąpił wycofany z użytku poprzedni Copernicus Open Access Hub (SciHub). 2023. CDSE oferuje natychmiastowy dostęp do pełnego archiwum Sentinel-2, z opóźnieniem w pozyskiwaniu danych najnowsze obrazy dostępne przez około 3 godziny.
Satelita Sentinel-2: Dane techniczne czujnika MSI
| Zespół | Długość fali (nm) | Rezolucja | Główna aplikacja |
|---|---|---|---|
| B1 - Aerozole przybrzeżne | 443 | 60 m | jakość powietrza, korekta atmosferyczna |
| B2 - Niebieski | 490 | 10 m | Rozróżnianie roślinności/gleby, kartowanie |
| B3 - Zielony | 560 | 10 m | Wigor roślinności, kolor liści |
| B4 -Czerwony | 665 | 10 m | Chlorofil, NDVI (czerwony pasek) |
| B5 - Roślinność czerwona krawędź | 705 | 20 m | Stres wegetacyjny, Red Edge NDVI |
| B6 - Roślinność czerwona krawędź | 740 | 20 m | Zawartość chlorofilu, klasyfikacja gatunkowa |
| B7 - Roślinność czerwona krawędź | 783 | 20 m | Biomasa, LAI |
| B8 - NIR | 842 | 10 m | NDVI (pasmo NIR), biomasa, LAI |
| B8a - Roślinność czerwona krawędź | 865 | 20 m | Konstrukcja czaszy, NDVI wąska |
| B9 - Para wodna | 940 | 60 m | Korekta atmosferycznej pary wodnej |
| B10 - SWIR - Cirrus | 1375 | 60 m | Wykrywanie chmur Cirrus |
| B11 - SWIR | 1610 | 20 m | Stres wodny, zawartość wody w liściach |
| B12 - SWIR | 2190 | 20 m | Zaawansowany stres wodny, starzenie się |
Istnieją głównie dwa poziomy przetwarzania danych Sentinel-2 istotne dla rolnictwa: Poziom 1C (L1C) - jasność górnej atmosfery z korekcją geometryczną, np Poziom 2A (L2A) - współczynnik odbicia na dnie atmosfery (Dół atmosfery, BOA) z korekcją atmosferyczną zastosowaną poprzez algorytm Sen2Cor. Do bezpośrednich zastosowań w rolnictwie, Poziom 2A jest zalecany, ponieważ jest już skorygowany pod kątem wpływu atmosfery porównywalne wartości współczynnika odbicia dla różnych dat i różnych lokalizacji geograficznych.
Platformy satelitarne dla AgriTech: pełne porównanie 2025
Sentinel-2 nie jest jedyną dostępną opcją. Rynek danych satelitarnych dla rolnictwa obejmuje dostawców komercyjnych, takich jak Planet Labs, operatorów zasiedziałych, takich jak Landsat (NASA/USGS) oraz zintegrowane platformy chmurowe, takie jak Google Earth Engine. Wybór zależy od złożonego kompromisu pomiędzy rozdzielczością przestrzenną, częstotliwością akwizycji, zarządzanym zachmurzeniem i dostępnym budżetem.
Porównanie platform satelitarnych dla rolnictwa precyzyjnego - 2025
| Platforma | Rozdzielczość przestrzenna | Ponownie odwiedź częstotliwość | Pasma widmowe | Koszt | Opóźnienie danych | Pszczoła |
|---|---|---|---|---|---|---|
| Sentinel-2 (ESA/Copernicus) | 10m (vis/NIR), 20m (SWIR) | 5 dni (2-3 w UE) | 13 pasm MSI | Bezpłatne (otwarte dane) | ~3 godziny od nabycia | CDSE OData, STAC, SentinelHub, OpenEO |
| Planeta PlanetSkop | 3,7 m | Codziennie (quasi-codziennie) | 8 pasm (PS2.SD) | Subskrypcja komercyjna; bezpłatnie na badania | 24 godziny | Planet API v1, API subskrypcji |
| Landsat 8/9 (NASA/USGS) | 30 m (wielospektralny), 15 m (panoramiczny) | 16 dni | 11 pasm OLI/TIRS | Bezpłatne (otwarte dane) | 24-48 godzin | USGS EarthExplorer, Google Earth Engine |
| MODIS (NASA Ziemia/Woda) | 250m / 500m / 1km | 1-2 dni | 36 zespołów | Bezpłatne (otwarte dane) | 24-48 godzin | NASA EarthData STAC, Google Earth Engine |
| Silnik Google Earth | Zależy od zbioru danych (10m-1km) | To zależy od zbioru danych | Zintegrowane wiele zbiorów danych | Bezpłatny, niekomercyjny; płatna reklama | Natychmiastowe przetwarzanie w chmurze | Python, API JavaScript |
| Maxar WorldView-3 | 0,3 m (panorama), 1,24 m (multi) | 1-4 dni | 29 pasm (CAVIS) | ~25-40 USD/km2 | 4-8 godzin | Maxar Streaming API |
| Airbusa Plejady Neo | 0,3 m | 1-2 dni | 6 pasm wielospektralnych | ~15-30 EUR/km2 | 24-48 godzin | API OneAtlas |
W przypadku zdecydowanej większości zastosowań rolniczych we Włoszech Sentinel-2 i wybór optymalny: rozdzielczość 10 metrów jest wystarczająca dla działek większych niż 0,5 ha (próg, poniżej którego statystyki pikseli stają się niewiarygodne), częstotliwość ponownych odwiedzin 2-3 dni we Włoszech i odpowiednie do sezonowego monitorowania, a także całkowicie bezpłatne eliminuje wszelkie bariery ekonomiczne w adopcji. Planet Labs staje się istotne tylko dla aplikacje wymagające niemal codziennego monitorowania (na przykład wykrywanie wczesnego pojawienia się infekcji). stres wodny w uprawach o wysokiej wartości, takich jak intensywne owoce i warzywa) lub rozdzielczość poniżej 10 m. MODIS i Landsat pozostają przydatne do analiz dużych obszarów i wielodekadowych szeregów czasowych.
Wskaźniki roślinności: NDVI, EVI, SAVI i LAI dla upraw włoskich
Wskaźniki roślinności (VI) to wskaźniki wyprowadzone matematycznie kombinacja pasm widmowych, które określają ilościowo określone właściwości roślinności. Wykorzystują fakt, że chlorofil silnie absorbuje w paśmie czerwonym (665 nm) i silnie odbija w bliskiej podczerwieni (842 nm): stosunek między tymi pasmami koduje gęstość i żywotność roślinność z matematyczną elegancją.
NDVI – znormalizowany różnicowy wskaźnik roślinności
NDVI to najczęściej stosowany na świecie wskaźnik roślinności, wprowadzony przez Rouse i in. w 1973 r. Wzór matematyczny to:
Formuła NDVI
NDVI = (NIR - Red) / (NIR + Red)
Su Sentinel-2:
NDVI = (B8 - B4) / (B8 + B4)
Range: da -1.0 a +1.0
Wartości NDVI są interpretowane zgodnie ze znormalizowaną skalą, ale optymalne progi są różne według fazy uprawnej i fenologicznej. Poniższa tabela przedstawia progi referencyjne dla głównych włoskich upraw, zaczerpnięte z literatury naukowej i walidacji empirycznej na typowych obszarach (Pianura Padana, Apulia, Sycylia, Toskania):
Interpretacja NDVI dla włoskich upraw
| Zakres NDVI | Interpretacja ogólna | Pszenica (Triticum) | Winorośl (Vitis) | Pomidor (Solanum) | Kukurydza (Zea Mays) | Oliwka (Olea) |
|---|---|---|---|---|---|---|
| < 0,1 | Goła ziemia, skały, śnieg | Niezasiana ziemia | Zima, przed pączkowaniem | Po zbiorach | Przedsiewnie | Goła ziemia między rzędami |
| 0,1 - 0,2 | Bardzo rzadka lub sucha roślinność | Początkowa sytuacja awaryjna | Uśpienie/silny stres | Niedawny przeszczep | Awaryjne (VE) | Silny stres wodny/termiczny |
| 0,2 - 0,4 | Rzadka roślinność, umiarkowany stres | Początkowe krzewienie | Przed kwitnieniem, lekki stres | Początkowy wzrost wegetatywny | Etapy V1-V3 | Słaba roślinność, stres |
| 0,4 - 0,6 | Umiarkowana roślinność, dobre zdrowie | Pełne powstanie | Przed kwitnieniem/zawiązaniem owoców | Aktywny rozwój wegetatywny | Etapy V4-V8 | Liście normalne, dobrze podlewane |
| 0,6 - 0,8 | Gęsta roślinność, doskonały wigor | Kurs (szczyt sezonowy) | Pełna foliacja, veraison | Maksymalne krycie (kwitnienie) | Etapy V10-VT (kwitnienie) | Gęste liście, stan doskonały |
| > 0,8 | Bardzo gęsta roślinność, las | Rzadkie (lasy na skraju pola) | Żywopłoty i drzewa graniczne | Szklarnie z całkowitym pokryciem | Rzadkie optymalne warunki | Gaje oliwne o dużym zagęszczeniu |
EVI – ulepszony wskaźnik roślinności
EVI, opracowany przez Huete i in. dla czujnika MODIS koryguje limity NDVI w obszarach o dużym zagęszczeniu roślinności (nasycenie NDVI powyżej 0,7) oraz na obszarach o odsłoniętych glebach (Błędy NDVI spowodowane odbiciem od podłoża). Formuła jest następująca:
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 – wskaźnik roślinności dostosowany do gleby
SAVI, zaproponowany przez Huete w 1988 roku, jest szczególnie przydatny w środowiskach półsuchych, takich jak Apulia lub Sycylia, gdzie odsłonięta gleba znacząco wpływa na odpowiedź widmową. Współczynnik korygujący L zmniejsza wpływ gruntu:
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 – wskaźnik powierzchni liści
LAI (wskaźnik powierzchni liści) określa całkowitą powierzchnię liścia na jednostkę powierzchni gleba. Jest to podstawowy parametr agronomiczny związany z potencjalną produktywnością i wydajnością pot. Sentinel-2 pozwala oszacować LAI za pomocą algorytmów empirycznych lub fizycznych (Procesor biofizyczny SNAP):
# 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 massimaOgraniczenia satelitarnych wskaźników roślinności
- Zachmurzenie: Sentinel-2 i chmury optyczne sprawiają, że obraz jest bezużyteczny. We Włoszech średnie zachmurzenie wynosi 48% zimą i 15% latem. Opaska QA60 umożliwia automatyczne maskowanie chmur.
- Nasycenie NDVI: Powyżej NDVI ~0,7-0,8 czułość drastycznie spada. W przypadku upraw o dużej gęstości użyj EVI.
- Mieszanka pikseli: Przy rozdzielczości 10 m piksel może zawierać zarówno roślinność, jak i glebę. Dla małych działek (< 0,5 ha) średni NDVI jest mniej reprezentatywny.
- Sezonowość fenologiczna: Porównywanie NDVI pomiędzy różnymi datami wymaga uwzględnienia fazy fenologicznej, a nie tylko wartości bezwzględnej.
- Zmienność międzyroczna: NDVI zmienia się również w zależności od sezonowych warunków pogodowych (temperatura, opady), niezależnie od praktyk agronomicznych.
Implementacja Pythona: Dostęp do CDSE i pobieranie Sentinel-2
Ekosystem przestrzeni danych programu Copernicus (CDSE) stał się pojedynczym portalem umożliwiającym dostęp do danych Sentinel od 2023. Oferuje cztery główne API: API OData (wyszukiwanie i pobieranie produktów), API STAC-a (Katalog zasobów przestrzenno-czasowych), API OpenEO (przetwarzanie standaryzowane) tj API Sentinel Hub (przetwarzanie jako usługa). Aby rozpocząć, utwórz darmowe konto na dataspace.copernicus.eu i wygeneruj dane uwierzytelniające OAuth2.
Konfiguracja środowiska Python i zależności
# 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 idricoUwierzytelnianie CDSE za pomocą 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)Wyszukaj produkty Sentinel-2 za pośrednictwem interfejsu API OData
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]}")Pobieraj i uzyskuj dostęp do pasm za pomocą Sentinel Hub Python
Do przetwarzania zespołu bez pobierania całej sceny (która waży 800 MB - 1,2 GB na plik). Produkt L2A), bardziej efektywne jest użycie interfejsu API procesu Sentinel Hub, który zwraca tylko pasma wymagane w interesującym Cię regionie:
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}%)")Obliczanie NDVI i zapisywanie 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")Interfejsy API pogodowe dla AgriTech: porównanie i wdrożenie
Dane pogodowe to drugie krytyczne źródło danych w rolnictwie precyzyjnym. Zintegrowane z danymi NDVI i IoT w terenie, umożliwiają budowanie modeli predykcyjnych deficyt wody, ryzyko chorób, ewapotranspiracja i prognoza plonów. Panorama interfejsów API pogody w 2025 r. dzieli się na bezpłatnych dostawców open source i dostawców komercyjnych z bezpłatnymi poziomami i usługami regionalnymi włoskiej ARPA.
Interfejsy API pogodowe dla AgriTech – porównanie 2025
| Dostawcy | Cena | Rozdzielczość przestrzenna | Prognoza | Zmienne agro | Archiwum Historyczne | Pszczoła |
|---|---|---|---|---|---|---|
| Otwarta pogoda | Bezpłatny (niekomercyjny); Reklama ~15 USD miesięcznie | 1-11 km (wiele modeli) | 16 dni | ET0, temperatura gleby, wilgotność gleby, ciśnienie pary | 1940-obecnie (ERA5) | REST JSON, klient Pythona |
| Otwórz mapę pogody | Bezpłatne do 60 połączeń/min; 40 USD/miesiąc pro | 2,5 km (IKONA) | 5 dni (bezpłatnie), 16 dni (pro) | Ograniczona, bez ET0 | Tylko płatna historia | REST JSON, SDK Pythona |
| Jutro.io | Bezpłatne 500 połączeń dziennie; Rdzeń 199 USD miesięcznie | 1 km | 21 dni | Wilgotność gleby, ET0, warunki oprysku, ryzyko szkodników | 6 lat | REST, WebSocket, gRPC |
| ARPA Lombardia / Piemont / Veneto | Bezpłatnie (regionalne otwarte dane) | Stacje punktualne (sieć ~1500 stacji we Włoszech) | Nie (tylko zaobserwowane) | Wszystkie zmienne rolno-pogodowe | Pełne archiwum (dekady) | Pobieranie WMS, REST, CSV/JSON |
| Meteoam/CFS/ECMWF | bezpłatne otwarte dane ECMWF; Bezpłatne NOAA CFS | 9-25 km | 10-45 dni | Zmienne pogodowe | 1940-obecnie za pośrednictwem ERA5 | API ECMWF (Python), API CDS |
| Przejście wizualne | Bezpłatne 1000 rekordów dziennie; Standardowo 35 USD miesięcznie | Interpolowane ze stacji | 15 dni | ET0, liczba dni uprawy, ryzyko przymrozków | Nielimitowane płatne | REST JSON, CSV, podobny do SQL |
W przypadku zastosowań AgriTech w produkcji zaleca się stosowanie Otwarta pogoda jak baza danych (bezpłatna, open source, dane ERA5 dla historyków, prognoza na 16 dni) integrująca ją z sieci stacji ARPA do lokalnej kalibracji. Open-Meteo uwzględnia zmienne agrometeorologiczne krytyczne kwestie, takie jak ET0 (ewapotranspiracja referencyjna), wilgotność gleby na różnych głębokościach i para wodna wodnych, które nie zawsze są dostępne bezpłatnie u dostawców komercyjnych.
Wdrożenie klienta Open-Meteo z funkcją Cache i Retry
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")Kompletny potok geoprzestrzenny: GeoPandas, rasterio i PostGIS
Profesjonalne zarządzanie danymi satelitarnymi wymaga kompleksowej infrastruktury geoprzestrzennej integrujący dane rastrowe (zdjęcia satelitarne NDVI, mapy pogodowe) z danymi wektorowymi (granice działki, obszary zagospodarowane, punkty poboru próbek). Potok, który tu opisujemy, przyjmuje de facto standardowy stos Pythona dla zastosowań geoprzestrzennych: raster dla rastrów, GeoPanda dla wektorów, PocztaGIS jako przestrzenna baza danych, np GDAL jako podstawowy silnik translacji formatów.
Architektura rurociągów geoprzestrzennych
Stos technologii rurociągów geoprzestrzennych AgriTech
┌─────────────────────────────────────────────────────────────────────────┐
│ 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 │
└─────────────────────────────────────────────────────────────────────────┘
Statystyki strefowe: Średnie NDVI na działkę ze statystykami rastrowymi
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)")Zapytania PostGIS do zaawansowanej analizy przestrzennej
-- 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;Model predykcyjny stresu wodnego z ML
Połączenie satelitarnego NDVI z danymi meteorologicznymi i IoT pozwala na budowanie modeli predykcyjnych niedoboru wody z 3–7-dniowym wyprzedzeniem, co pozwala na proaktywne planowanie nawadniania. Opisany tutaj model integruje trzy źródła danych: szeregi czasowe NDVI z Sentinel-2, bilans wodny wyliczony z ET0 i opadów oraz wilgotność gleby z czujników IoT (lub z Open-Meteo ERA5-Land, jeśli czujniki fizyczne nie są dostępne).
Inżynieria cech modelu
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"]]Studium przypadku: Monitoring NDVI w winnicy Primitivo w Mandurii
Aby ukonkretnić opisane koncepcje, przedstawiamy studium przypadku inspirowane aplikacją prawdziwe na terytorium DOC Primitivo di Manduria (Taranto, Apulia). Teren ma warunki typowe dla śródziemnomorskiej uprawy winorośli: gorące i suche lata (opady w lecie < 50 mm), gleby gliniasto-wapienne, przeciwszpalery i systemy do sadzenia drzewek.
Parametry studium przypadku - Vigna Primitivo Manduria
| Parametr | Wartość |
|---|---|
| Lokalizacja | Manduria (TA), Apulia – szer. 40,38, dł. 17,47 |
| Całkowita powierzchnia | 35 ha podzielone na 12 działek |
| Główne odmiany | Primitivo (80%), Negroamaro (20%) |
| System hodowlany | Sadzonka Apulii (8 szt.), Przeciwespalier (4 szt.) |
| Rok rośliny | 2003-2015 (winnica mieszana młoda/dorosła) |
| Nawadnianie | Kropla po kropli (nie wszystkie działki) |
| Śledzenie satelity | Sentinel-2A/B, L2A, rozdzielczość 10m |
| Stosowana częstotliwość przejęć | 28 scen na wsi 2024 (chmura < 20%) |
| Zintegrowane czujniki IoT | 6 stacji pogodowych, 18 czujników wilgotności gleby na wysokości 30/60 cm |
| Okres analizy | kwiecień – październik 2024 r |
Sezonowe wyniki NDVI i decyzje dotyczące nawadniania
Sezonowa analiza NDVI ujawniła znacząco zróżnicowane wzorce stresu wodnego między działkami, co ma bezpośredni wpływ na zarządzanie nawadnianiem i jakość produktu:
Wyniki NDVI według fabuły – kampania 2024
| Działka | Różnorodność | NDVI kwiecień | NDVI w dół | NDVI lipiec (szczyt) | NDVI sierpień | Zestaw NDVI | Alarm dotyczący nawadniania |
|---|---|---|---|---|---|---|---|
| Parkuj A1 | Prymitywne / drzewko | 0,28 | 0,42 | 0,58 | 0,41 | 0,31 | Umiarkowane naprężenie igły |
| Zaparkuj A2 | Prymitywne / drzewko | 0,22 | 0,36 | 0,47 | 0,29 | 0,24 | Silny stres związany z igłą |
| Park B1 | Prymitywny/kontrapalier | 0,31 | 0,53 | 0,66 | 0,57 | 0,44 | Brak alertów |
| Park B2 | Prymitywny/kontrapalier | 0,29 | 0,49 | 0,63 | 0,52 | 0,40 | Brak alertów |
| Zaparkuj C1 | Negroamaro/drzewo | 0,25 | 0,44 | 0,55 | 0,38 | 0,28 | Umiarkowane naprężenie igły |
| Zaparkuj C2 | Negroamaro/kontrpalliera | 0,30 | 0,51 | 0,61 | 0,50 | 0,38 | Brak alertów |
Powstały wzór jest wyraźny i ma znaczenie agronomiczne: działki sadzonek są pozbawione nawadnianie awaryjne wykazuje bardziej wyraźny niedobór wody w okresie od sierpnia do września w porównaniu z nawadnianiem awaryjnym Działki przeciwespalierowe z systemem kroplowym. W szczególności działka A2 ma wykazał spadek NDVI z 0,47 do 0,29 w okresie od lipca do sierpnia 2024 r., czasowo skorelowany z sekwencja 23 dni bez opadów i maksymalnych temperatur powyżej 36 stopni.
Automatyczny skrypt ostrzegający
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 Satelitarnego Systemu Monitoringu
Kwantyfikacja zwrotu z inwestycji w monitorowanie satelitarne ma fundamentalne znaczenie uzasadniają inwestycje technologiczne w przedsiębiorstwa zajmujące się uprawą winorośli. Za winnicę Mandurii, analiza kosztów i korzyści kampanii na rok 2024 pokazuje:
Analiza ROI - System Monitoringu Satelitarnego Vigna Manduria (35 ha, kampania 2024)
| Głos | Wartość | Notatki |
|---|---|---|
| KOSZTY | ||
| Dostęp do danych Sentinel-2 (CDSE) | 0 EUR | Bezpłatne otwarte dane |
| Interfejs API pogody typu Open Weather | 0 EUR | Niekomercyjne, bezpłatne |
| Rozwój potoku w Pythonie | 2500 EUR | 30 godzin programistycznych x 83 EUR/h (jednorazowe, amortyzowane przez 3 lata) |
| Hosting w chmurze (VM + pamięć masowa) | 600 EUR/rok | VM 4 vCPU + 50 GB pamięci S3 |
| Integracyjne doradztwo agronomiczne | 1000 EUR | Jednorazowa kalibracja progowa dla konkretnej uprawy |
| Całkowity koszt za rok 1 | 4100 EUR | Kolejne lata: ~1400 EUR/rok |
| SZACOWANE KORZYŚCI (35 ha, 100 000 butelek/rok) | ||
| Redukcja nawadniania (optymalizacja czasu) | 2800 EUR | 15% oszczędności wody, 8500 m3 x 0,33 EUR/m3 |
| Ograniczenie zabiegów fitosanitarnych | 1750 EUR | Ukierunkowana interwencja tylko w obszarach ryzyka (-20% zabiegów) |
| Poprawa jakości winogron | 8500 EUR | +0,5 pkt. średniego stopnia Brixa na 12 ha, +0,20 EUR/kg premii jakościowej |
| Redukcja strat spowodowanych stresem wodnym | 3200 EUR | 8% odzysk plonów na 2 działkach przy oczekiwanym silnym stresie |
| Suma świadczeń rok 1 | 16 250 EUR | |
| Rok 1 ROI | 296% | Zwrot w ciągu około 3 miesięcy |
Google Earth Engine: przetwarzanie w chmurze do analiz na dużą skalę
Do analiz w skali regionalnej lub krajowej (np. mapowanie odmian na poziomie DOC powiat, ocena szkód spowodowanych przymrozkami wiosennymi w skali województwa, monitoring WPR dla agencji płatniczych, takich jak AGEA), Google Earth Engine (GEE) oferuje moc obliczeniową natywne w chmurze, które w przeciwnym razie wymagałyby dedykowanej 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}")Regulacje i otwarte dane dla AgriTech we Włoszech
Rolniczy ekosystem otwartych danych we Włoszech opiera się na trzech głównych aktorach: WIEK (Agencja Płatności Rolnych), ul SIAN (Krajowy System Informacji Rolniczej) i powiązane z nim portale regionalne. W lutym 2025 r. AGEA zakończyła migrację SIAN w kierunku Narodowego Bieguna Strategicznego, aktualizacja infrastrukturę chmurową i poszerzanie możliwości przetwarzania danych geoprzestrzennych.
Główne włoskie otwarte źródła danych dotyczące rolnictwa
| Źródło | Adres URL | Dostępne dane | Format | Aktualizacja |
|---|---|---|---|---|
| Otwarte dane SIAN | data.sian.it | Ortofotomapy w wysokiej rozdzielczości, mapa zagospodarowania przestrzennego (oparta na sztucznej inteligencji), korporacyjne plany graficzne | GeoTIFF, plik kształtu, WFS | Coroczny |
| Usługi AGEA EO | sian.it | Ortofotomapy wielospektralne, obrazy radarowe w paśmie X, DEM, DSM | GeoTIFF, ECW | Czasopismo (kampania) |
| ISTAT – Statystyka Rolnicza | ist.it/agricoltura | Powierzchnie użytków rolnych, produkcja, struktura firmy, spis rolny | CSV, JSON, SDMX | Coroczny |
| Krajowy Geoportal MATTM | geoportale.gov.it | Carta Natura, Corine Land Cover, IDT, DBT | WMS, WFS, Shapefile | Zmienny |
| Regionalne ARPA (15 regionów) | harfa.[region].it | Dane stacji metoklimatycznych, regionalne modele prognostyczne | CSV, JSON, NetCDF | Prawie w czasie rzeczywistym |
| Służba Ziemi Kopernika | land.copernicus.eu | CORINE Pokrycie terenu, użytki zielone, woda i wilgoć, Atlas miejski | GeoTIFF, plik kształtu | Trzyletnie/roczne |
Ulgi podatkowe PNRR – przejście 5.0 dla AgriTech
Inwestycje w satelitarne systemy monitorowania i rolnictwo precyzyjne mogą zakwalifikować się do świadczeń Planu Przejściowego 5.0 (Ustawa 207/2024 – Ustawa budżetowa 2025). Przedsiębiorcy rolni, dopuszczeni od 2025 r., będą mogli skorzystać z ulg podatkowych na inwestycje w materialne i niematerialne dobra kapitałowe 4.0, w tym oprogramowanie do zarządzania przedsiębiorstwem komponenty sztucznej inteligencji i analiza danych geoprzestrzennych.
Przejście 5.0 Stawki ulg podatkowych dla przedsiębiorstw rolniczych (2025)
| Zakres inwestycji | Oszczędność energii 3-6% | Oszczędności 6-10% | Oszczędności >10% |
|---|---|---|---|
| Do 2,5 mln EUR | 35% | 40% | 45% |
| Od 2,5 mln do 10 mln EUR | 15% | 20% | 25% |
| Od 10M do 50M EUR | 5% | 10% | 15% |
Kwalifikują się precyzyjne systemy nawadniające sterowane przez NDVI i prognozy pogody dla przedziału „oszczędność energii >6%”, jeżeli dokumentują zmniejszenie zużycia wody (pompowanie energii) w porównaniu z wartością bazową. Zapytanie o ekspertyzę techniczną pod przysięgą.
Architektura produkcji: najlepsze praktyki i anty-wzorce
Po opisaniu komponentów ważne jest podsumowanie najlepszych praktyk dotyczących solidną implementację w produkcji i identyfikację wspólnych anty-wzorców prowadzić do kruchości lub niedokładności systemów.
Najlepsze praktyki - System monitoringu satelitarnego AgriTech
- Progresywne filtrowanie chmur: Nie używaj tylko procentu zachmurzenia produktu (metadane poziomu 1). Zawsze stosuj maskowanie na poziomie pikseli za pomocą QA60 (Sentinel-2) lub BQA (Landsat). Produkt chmurowy o stężeniu 5% może mieć określoną działkę całkowicie przesłoniętą.
- Zoptymalizowany w chmurze GeoTIFF (COG): Przechowuj rastry NDVI w formacie COG, aby zapewnić efektywny dostęp do żądań zasięgu z S3/MinIO. Unikaj klasycznych plików GeoTIFF, które nie są zoptymalizowane pod kątem przechowywania w chmurze.
- Spójna projekcja: Standaryzuj WSZYSTKO zgodnie z EPSG:4326 (WGS84) w celu uzyskania danych końcowych lub ze strefą UTM 32N (EPSG:32632) dla Włoch, aby zachować prawidłowe pomiary metryczne. Nigdy nie mieszaj CRS bez wyraźnej nagany.
- Walidacja danych w oparciu o prawdę IoT: Kalibracja progów NDVI za pomocą pomiarów terenowych (czujniki wilgotności gleby, przenośny miernik LAI, zebrane dane agronomiczne). Samo NDVI bez lokalnej kalibracji może dawać fałszywe alarmy dodatnie/ujemne w zakresie 15–30%.
- Zarządzanie archiwum czasu: Utrzymuj szeregi czasowe NDVI przez co najmniej 3–5 lat dla każdego wykresu. Anomalia NDVI w porównaniu ze średnią historyczną (ten sam miesiąc, ten sam zbiór) dostarcza znacznie więcej informacji niż wartość bezwzględna.
- API ograniczania szybkości i buforowania: CDSE i Open-Meteo mają limity stawek. Zawsze wdrażaj lokalną pamięć podręczną (opartą na plikach lub Redis), aby uniknąć niepotrzebnego ponownego pobierania. Pobieranie Sentinel-2 zajmuje 30–120 sekund: nie żądaj dwukrotnie tych samych danych.
- Rejestrowanie i ścieżka audytu: Każde obliczenie NDVI musi być powiązane ze źródłową sceną satelitarną, datą pozyskania, procentem zachmurzenia i wersją algorytmu. Podstawa audytów PAC/AGEA i debugowania anomalii.
Powszechne anty-wzorce - systemy satelitarne AgriTech
- Pomijając korektę atmosferyczną: Stosowanie poziomu 1C (góra atmosfery) zamiast poziomu 2A (dno atmosfery) wprowadza systematyczne błędy w NDVI rzędu 10-20% w warunkach dużej wilgotności powietrza (mgła w dolinie Padu, Morze Adriatyckie). ZAWSZE używaj L2A do zastosowań ilościowych.
- NDVI jako jedyny wskaźnik: Polegaj wyłącznie na NDVI, ignorując EVI (dla gęstej roślinności), NDWI (stres wodny), Red Edge NDVI (wczesny stres, zanim będzie widoczny na NDVI). Profesjonalny system wykorzystuje kombinację co najmniej 3-4 indeksów.
- Rozmiar piksela i paczka: Zastosowanie NDVI do działek mniejszych niż 0,3 ha z pikselami oddalonymi o 10 m daje statystycznie niewiarygodne średnie (poniżej 30 pikseli). W przypadku małych działek rozważ Planet Labs (3,7 m) lub czujniki dronów.
- Brak kompozycji czasowej: Przyjmowanie pojedynczego akwizycji dostępnego dla danej daty, nawet przy zachmurzeniu 19%, prowadzi do artefaktów NDVI w obszarach częściowego zachmurzenia. Zawsze używaj kompozycji medoidowej w oknie czasowym 10-15 dni.
- Potok bez wersjonowania danych: Ponowne obliczenie NDVI z różnymi wersjami Pythona/rasterio/sentinelhub daje nieco inne wartości dla tych samych danych źródłowych. Wersjonuj swoje potoki obliczeniowe za pomocą DVC lub równoważnego.
Kompletny potok: orkiestracja z prefektem
Dzięki integracji wszystkich opisanych komponentów powstaje kompletny potok śledzenia satelitarnego można zaaranżować za pomocą Prefect (zobacz serię Orkiestracja rurociągów) uruchamiać się automatycznie co 5 dni (z częstotliwością ponownych odwiedzin 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()Wnioski i dalsze kroki
Dane satelitarne Sentinel-2 są dostępne bezpłatnie za pośrednictwem CDSE prawdopodobnie najbardziej niedostatecznie wykorzystywane źródło danych we włoskim AgriTech. Z Pythonem, rasterio, Sentinelhub i konto CDSE, możliwe jest zbudowanie konta System monitorowania NDVI, który przewyższa wiele rozwiązań pod względem jakości i szczegółowości komercyjne za opłatą.
Kluczem do sukcesu nie jest zaawansowanie technologiczne, ale kalibracja lokalna agronomia: progi NDVI, okna czasowe alertów, wagi cechy modelu predykcyjnego muszą być skalibrowane pod kątem konkretnej uprawy, na lokalny klimat i praktyki agronomiczne firmy. Kalibrowany system w przypadku winorośli w Apulii bez niego nie sprawdzi się równie dobrze w przypadku pszenicy w Dolinie Padu interwencja człowieka w zakresie kalibracji.
Kluczowe punkty do zapamiętania dla tych, którzy chcą wdrożyć podobny system:
- Zawsze używaj Strażnik-2 Poziom-2A (zastosowano korektę atmosferyczną) do porównań ilościowych pomiędzy różnymi datami.
- Połącz NDVI z co najmniej jednym indeksem uzupełniającym (EVI dla gęstych winnic, NDWI dla wczesnego stresu wodnego, Red Edge NDVI dla monitorowania chlorofilu).
- Zintegrować Otwarta pogoda ERA5-Land dla bilansu wodnego: ET0 + opady + wilgotność gleby z modelu i doskonałe zastępstwo, gdy nie masz czujników IoT w terenie.
- Zachowaj jeden seria historyczna trwająca co najmniej 3 lata dla każdej działki: anomalia porównana ze średnią historyczną jest znacznie bardziej użyteczna niż bezwzględna wartość NDVI.
- Skaluj za pomocą Silnik Google Earth tylko wtedy, gdy trzeba przetwarzać na skalę regionalną lub krajową: w przypadku indywidualnych firm lokalny potok Pythona na CDSE jest bardziej kontrolowany i nie wymaga blokowania.
Seria FoodTech trwa
Zbudowałeś rurociąg satelitarny. Teraz zapoznaj się z innymi artykułami z tej serii uzupełnij architekturę AgriTech:
- Artykuł 1: Rurociąg IoT dla rolnictwa precyzyjnego z Pythonem i MQTT - Jak zintegrować czujniki naziemne z rurociągiem satelitarnym.
- Artykuł 2: ML Edge do monitorowania upraw: wizja komputerowa na polach - Komputerowe modele wizyjne do wczesnego wykrywania chorób roślin.
- Artykuł 4: Identyfikowalność Blockchain w żywności: od pola do supermarketu - W jaki sposób dane NDVI zasilają certyfikaty jakości w łańcuchu.
Aby dowiedzieć się więcej o MLOps i wdrażaniu modeli predykcyjnych opisanych w tym artykule, zobacz serial MLOps dla firm z MLflow i serial LLM w biznesie: RAG Enterprise.