AgriTech 向けの衛星および天気 API: Sentinel-2 と Python による予測データ
2024 年に、欧州宇宙機関は 1.5ペタバイトのデータ 毎日 コペルニクス計画のセンチネル衛星から。このうち農業モニタリングに特化した部分 ヨーロッパでは、これは最大かつ最も戦略的に関連性の高いセグメントを表します。それでも大多数は 最も技術的に進んだ企業を含むイタリアの農業企業のうち、まだ農業生産を行っていない企業の割合 この公共、無料、科学的に検証されたインフラストラクチャ。
その理由はデータの不足ではなく、データへのアクセス、処理、変換の複雑さです。 具体的な農業上の決定。 Sentinel-2 衛星は 5 日ごとにブドウ畑の上空を通過します。 解像度はピクセルあたり 10 メートルです。各パスは健康をコード化する 13 のスペクトル バンドを生成します 植生、水ストレス、クロロフィル含有量、寄生虫の存在。データと組み合わせる 天気、IoT 地上センサー、予測モデルなど、このデータは根本的に変革をもたらす可能性があります 農業に関する意思決定の質を高め、投入コストを 15 ~ 25% 削減し、収量を増加させます。 最大10〜20%。
世界の地球観測市場には価値がある 2025年には100.7億ドル Straits Research によると、2033 年までに 172 億人 (CAGR 6.92%) に成長すると予想されています。 農業部門はアプリケーションの 21% を占め、最も高い割合で成長しています。 精密農業、収量監視、灌漑の最適化に対する需要が原動力となっています。 この記事では、Sentinel-2 にアクセスするための完全な Python パイプラインを段階的に構築します。 CDSE (Copernicus Data Space Ecosystem) を介して、NDVI およびその他の植生指数を計算し、統合します Open-Meteo からの気象データを利用し、水ストレスの予測モデルをフィードします。
この記事で学べること
- コペルニクス プログラムのアーキテクチャと CDSE 経由で Sentinel-2 に無料でアクセスする方法
- 植生指数: NDVI、EVI、SAVI、LAI - イタリアの作物の数式、解釈、およびしきい値
- 完全な Python 実装: CDSE 認証、バンドのダウンロード、rasterio と numpy による NDVI 計算
- 天気 API の比較: Open-Meteo、OpenWeatherMap、Tomorrow.io、衛星データ統合
- GeoPandas、rasterio、PostGIS を使用したベクトルおよびラスター データ用の地理空間パイプライン
- 水ストレス予測モデル: Sentinel-2 + 気象 + IoT と scikit-learn の組み合わせ
- 実践的なケーススタディ: プーリア州のブドウ園、季節ごとの NDVI モニタリング、灌漑の最適化
- イタリアの法律: AGEA、SIAN、オープンデータ CAP、デジタル CAP
フードテックシリーズ - すべての記事
| # | アイテム | レベル | Stato |
|---|---|---|---|
| 1 | Python と MQTT を使用した精密農業用の IoT パイプライン | 高度な | 利用可能 |
| 2 | 作物監視のための ML Edge: 圃場でのコンピュータ ビジョン | 高度な | 利用可能 |
| 3 | アグリテック向けの衛星および天気 API: 予測データ (ここにいます) | 高度な | 現在 |
| 4 | 食品におけるブロックチェーンのトレーサビリティ: 現場からスーパーマーケットまで | 中級 | 近日公開 |
| 5 | 食品産業における品質管理のためのコンピュータービジョン | 高度な | 近日公開 |
| 6 | FSMA とデジタル コンプライアンス: 規制プロセスの自動化 | 中級 | 近日公開 |
| 7 | 垂直農法: IoT と ML による環境制御 | 高度な | 近日公開 |
| 8 | Prophet と LightGBM を使用した食品小売の需要予測 | 中級 | 近日公開 |
| 9 | ファーム インテリジェンス ダッシュボード: Grafana を使用したリアルタイム分析 | 中級 | 近日公開 |
| 10 | サプライチェーンの食品の最適化: 廃棄物削減のための ML | 中級 | 近日公開 |
コペルニクス計画: 地球観測のためのヨーロッパのインフラストラクチャ
欧州連合のコペルニクス計画は、世界最大の観測インフラです。 土地は決して建設されませんでした。 ESA (欧州宇宙機関) と欧州委員会によって管理されています。 EUMETSAT の運用サポート、Copernicus は名前付き衛星群を運用します 「センチネル」はそれぞれ特定の監視任務のために設計されています。農業の場合、 絶対基準衛星e センチネル-2.
Sentinel-2 コンステレーションは、2 つの双衛星 (Sentinel-2A と Sentinel-2B) で構成されています。 高度786kmの太陽同期軌道。 2 つの衛星構成により、 時間 5日間のレビュー 赤道上、ヨーロッパ緯度では 2 ~ 3 日 (イタリア) 含まれています)。各衛星には MSI (MultiSpectral Instrument) センサー、手押しほうきが搭載されています。 290kmのスワイプ幅と13のスペクトルバンドでデータを取得するスキャナー 10 メートル、20 メートル、60 メートルの 3 つの空間解像度に分散されます。
Sentinel-2 の最も戦略的な特徴は、 完全無料:2014年から コペルニクスのすべてのデータは、商用または非商用を問わず、あらゆる用途にオープンアクセスで利用できます。 コマーシャル。 2023 年からの新しいアクセス ポータルと コペルニクス データ空間エコシステム (CDSE)、で廃止された以前の Copernicus Open Access Hub (SciHub) に代わるものです。 2023. CDSE は、取得レイテンシーを伴いながら、完全な Sentinel-2 アーカイブへの即時アクセスを提供します 最新のイメージを利用するには約 3 時間かかります。
Satellite Sentinel-2: MSI センサーの技術仕様
| バンド | 波長(nm) | 解決 | 主な用途 |
|---|---|---|---|
| B1 - 海岸エアロゾル | 443 | 60メートル | 空気の質、大気補正 |
| B2 - 青 | 490 | 10メートル | 植生/土壌の識別、マッピング |
| B3 - 緑 | 560 | 10メートル | 草勢、葉色 |
| B4 -赤 | 665 | 10メートル | クロロフィル、NDVI (赤いバンド) |
| B5 - 植生レッドエッジ | 705 | 20メートル | 植生ストレス、レッドエッジ NDVI |
| B6 - 植生レッドエッジ | 740 | 20メートル | クロロフィル含有量、種分類 |
| B7 - 植生レッドエッジ | 783 | 20メートル | バイオマス、LAI |
| B8 - 近赤外線 | 842 | 10メートル | NDVI(近赤外帯)、バイオマス、LAI |
| B8a - 植生レッドエッジ | 865 | 20メートル | キャノピー構造、NDVI ナロー |
| B9 - 水蒸気 | 940 | 60メートル | 大気中の水蒸気補正 |
| B10 - SWIR - 巻雲 | 1375年 | 60メートル | 巻雲の検出 |
| B11 - SWIR | 1610 | 20メートル | 水分ストレス、葉水分量 |
| B12 - SWIR | 2190 | 20メートル | 高度な水分ストレス、老化 |
農業に関連する Sentinel-2 データ処理には主に 2 つのレベルがあります。 レベル 1C (L1C) - 幾何学補正によるトップ大気の輝き、e レベル2A(L2A) - 大気底での反射率 (大気底、BOA) Sen2Cor アルゴリズムを介して大気補正が適用されます。直接的な農業用途の場合、 レベル2Aは大気の影響が補正されているのでおすすめです。 異なる日付と異なる地理的位置の間で比較可能な反射率値。
アグリテック向け衛星プラットフォーム: 2025 年の完全な比較
Sentinel-2 が唯一の利用可能なオプションではありません。農業向け衛星データ市場 Planet Labs などの民間プロバイダー、Landsat (NASA/USGS) などの既存プロバイダーが含まれます。 Google Earth Engineなどの統合クラウドプラットフォーム。選択は複雑なトレードオフに依存します 空間解像度、取得頻度、管理されたクラウド カバー、利用可能な予算の間で調整します。
精密農業向け衛星プラットフォームの比較 - 2025 年
| プラットフォーム | 空間解像度 | 再訪問の頻度 | スペクトルバンド | 料金 | データ遅延 | API |
|---|---|---|---|---|---|---|
| センチネル 2 (ESA/コペルニクス) | 10m (可視/近赤外)、20m (SWIR) | 5日間(EUでは2~3日) | 13のMSIバンド | 無料(オープンデータ) | 取得から約 3 時間 | CDSE OData、STAC、SentinelHub、OpenEO |
| プラネット プラネットスコープ | 3.7メートル | 毎日(準毎日) | 8バンド(PS2.SD) | 商用サブスクリプション;研究のために無料で | 24時間 | Planet API v1、サブスクリプション API |
| Landsat 8/9 (NASA/USGS) | 30 m (マルチスペクトル)、15 m (パン) | 16日 | 11バンドOLI/TIRS | 無料(オープンデータ) | 24~48時間 | USGS EarthExplorer、Google Earth エンジン |
| MODIS (NASA アース/アクア) | 250m / 500m / 1km | 1~2日 | 36バンド | 無料(オープンデータ) | 24~48時間 | NASA EarthData STAC、Google Earth Engine |
| Google Earth エンジン | データセットに依存 (10m ~ 1km) | データセットによって異なります | 統合された複数のデータセット | 無料、非営利。有料コマーシャル | インスタントクラウド処理 | Python ee、JavaScript API |
| マクサー ワールドビュー-3 | 0.3m(パン)、1.24m(マルチ) | 1~4日 | 29バンド(CAVIS) | ~25-40 USD/km2 | 4~8時間 | Maxar ストリーミング API |
| エアバス プレアデス ネオ | 0.3m | 1~2日 | 6つのマルチスペクトルバンド | ~15-30 ユーロ/km2 | 24~48時間 | OneAtlas API |
イタリアの農業用途の大部分では、 センチネル 2 とその選択 最適な: 0.5 ヘクタールを超える区画には 10 メートルの解像度で十分です (ピクセル統計が信頼できなくなるしきい値)、再訪問頻度 イタリアでは 2 ~ 3 日で、季節のモニタリングには十分です。また、完全に無料です。 導入に対する経済的障壁を排除します。 Planet Labs は次の場合にのみ関連性を持ちます。 ほぼ毎日の監視を必要とするアプリケーション (例: 早期発症の検出など) 集中的な果物や野菜などの高価値作物における水ストレス)または 10 メートル未満の解像度。 MODIS と Landsat は、大規模なエリアや数十年にわたる時系列の分析に引き続き役立ちます。
植生指数: イタリアの作物に対する NDVI、EVI、SAVI、LAI
植生指数 (VI) は、以下から数学的に導出される指標です。 特定の植生特性を定量化するスペクトル バンドの組み合わせ。彼らは搾取します クロロフィルは赤色帯域(665nm)を強く吸収し、強く反射するという事実 近赤外線 (842 nm) では、これらのバンド間の比率が、その密度と活力をコード化します。 数学的な優雅さを備えた植物。
NDVI - 正規化された差異植生指数
NDVI は、Rouse らによって導入された、世界で最も使用されている植生指数です。 1973年に。 数式は次のとおりです。
NDVIの式
NDVI = (NIR - Red) / (NIR + Red)
Su Sentinel-2:
NDVI = (B8 - B4) / (B8 + B4)
Range: da -1.0 a +1.0
NDVI 値は標準化されたスケールに従って解釈されますが、最適なしきい値は異なります 作物および季節学的段階ごとに。次の表は、次の基準しきい値を示しています。 科学文献と経験的検証に基づいたイタリアの主要作物 典型的な地域 (ピアヌーラ パダナ、プーリア、シチリア島、トスカーナ):
イタリア作物に対する NDVI の解釈
| NDVI 範囲 | 一般的な解釈 | コムギ(コムギ) | つる (ブドウ) | トマト(ナス) | トウモロコシ(ゼアメイズ) | オリーブ(オレア) |
|---|---|---|---|---|---|---|
| < 0.1 | 裸地、岩、雪 | 播種されていない土地 | 冬、芽吹き前 | 収穫後 | 播種前 | 列間の裸地 |
| 0.1~0.2 | 非常にまばらな植物または乾燥した植物 | 初期の緊急事態 | 休眠/重度のストレス | 最近の移植 | 緊急時(VE) | 深刻な水/熱ストレス |
| 0.2~0.4 | 植生がまばら、適度なストレス | 初期分げつ | 開花前、軽いストレス | 初期の栄養成長 | ステージ V1 ~ V3 | 貧弱な植生、ストレス |
| 0.4~0.6 | 適度な植生、健康状態良好 | フルライジング | 開花前・結実 | 活発な栄養発育 | ステージ V4 ~ V8 | 通常の葉、よく水やり |
| 0.6~0.8 | 植物が密生し、樹勢が優れている | ヘディング(季節のピーク) | フルフォレーション、ベレゾン | 最大被覆率(開花) | ステージ V10~VT(開花) | 葉が濃く、状態は良好 |
| > 0.8 | 非常に密集した植物、森林 | レア(畑の端の森) | 生垣と境界木 | 全体をカバーする温室 | まれな最適条件 | 高密度のオリーブ畑 |
EVI - 強化された植生インデックス
Hueteらによって開発されたEVI。 MODIS センサーの場合、領域内の NDVI 制限を修正します。 植生密度が高く(NDVI飽和度が0.7以上)、土壌が露出している地域 (地面の反射率によって生じる NDVI 誤差)。式は次のとおりです。
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 - 土壌調整植生指数
1988 年に Huete によって提案された SAVI は、次のような半乾燥環境で特に役立ちます。 露出した土壌がスペクトル応答に大きな影響を与えるプーリアまたはシチリア。 補正係数 L は地面の影響を軽減します。
SAVI = (NIR - Red) * (1 + L) / (NIR + Red + L)
Su Sentinel-2:
SAVI = (B8 - B4) * (1 + 0.5) / (B8 + B4 + 0.5)
L = 0.5 (valore standard per vegetazione media)
L = 1.0 per vegetazione molto rada
L = 0.25 per vegetazione densa
Vantaggioso rispetto a NDVI quando:
- Copertura vegetale < 40%
- Suoli chiari e brillanti (calcari, sabbie)
- Areali semi-aridi del Sud Italia
LAI - 葉面積指数
LAI (葉面積指数) は、単位表面積あたりの葉の総表面積を定量化します。 土。これは潜在的な生産性と関連する基本的な農業パラメータであり、 汗。 Sentinel-2 を使用すると、経験的アルゴリズムまたは物理的アルゴリズムを使用して LAI を推定できます。 (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 massima
衛星植生指数の限界
- 曇り: Sentinel-2 と光学雲により画像が使用できなくなります。イタリアの平均雲量は冬で 48%、夏で 15% です。 QA60 バンドでは自動クラウド マスキングが可能です。
- NDVI飽和: NDVI ~0.7 ~ 0.8 を超えると、感度は大幅に低下します。高密度の作物にはEVIを使用してください。
- ピクセルの混合: 解像度 10 m では、ピクセルに植生と土壌の両方を含めることができます。小さな区画 (< 0.5 ha) の場合、平均 NDVI はあまり代表的ではありません。
- 季節学的季節性: 異なる日付間の NDVI を比較するには、絶対値だけでなく季節相を考慮する必要があります。
- 年ごとの変動: NDVI は、農業慣行に関係なく、季節の気象条件 (気温、降水量) によっても変化します。
Python の実装: CDSE にアクセスして Sentinel-2 をダウンロードする
Copernicus Data Space Ecosystem (CDSE) は、Sentinel データにアクセスするための単一のポータルになりました 2023 年以降。次の 4 つの主要な API を提供します。 OData API (製品の検索とダウンロード)、 STAC API (時空間資産カタログ)、 OpenEO API (処理中 標準化された)e Sentinel ハブ API (サービスとして処理)。始めるには、 dataspace.copernicus.eu で無料アカウントを作成し、OAuth2 資格情報を生成します。
Python 環境のセットアップと依存関係
# Installa le librerie necessarie
pip install sentinelhub rasterio numpy pandas geopandas shapely \
requests python-dotenv matplotlib scipy scikit-learn \
openmeteo-requests retry-requests xarray
# requirements.txt con versioni validate per Python 3.11+
# sentinelhub==3.11.4 - interfaccia CDSE e SentinelHub
# rasterio==1.4.0 - lettura/scrittura raster geospaziali
# numpy==2.1.0 - calcolo matriciale per bande
# geopandas==1.0.1 - dati vettoriali e spatial operations
# shapely==2.0.6 - geometrie per AOI (Area of Interest)
# openmeteo-requests==0.3.3 - client Open-Meteo API
# scikit-learn==1.5.2 - modello predittivo stress idrico
OAuth2 による CDSE 認証
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)
OData API 経由で Sentinel-2 製品を検索
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]}")
Sentinel Hub Python を使用してバンドをダウンロードしてアクセスする
シーン全体をダウンロードせずにバンドを処理する場合 (シーンごとに 800 MB ~ 1.2 GB の重さ) L2A 製品)、 対象領域で必要なバンド:
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}%)")
NDVI 計算と 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")
アグリテック向け気象 API: 比較と実装
気象データは、精密農業の 2 番目に重要なデータ ソースです。 地上の NDVI および IoT データと統合されているため、予測モデルを構築できます。 水ストレス、病気のリスク、蒸発散量、収量予測。パノラマ 2025 年の気象 API は、無料のオープンソース プロバイダーと商用プロバイダーに分けられます イタリア ARPA の無料枠と地域サービスを利用できます。
アグリテック向け気象 API - 2025 年の比較
| プロバイダー | 価格 | 空間解像度 | 予報 | 農業変数 | 歴史的アーカイブ | API |
|---|---|---|---|---|---|---|
| オープンウェザー | 無料(非営利)。コマーシャル月額 $15 | 1~11km(マルチモデル) | 16日 | ET0、土壌温度、土壌水分、蒸気圧 | 1940年~現在(ERA5) | REST JSON、Python クライアント |
| 天気マップを開く | 1 分あたり最大 60 件の通話が無料。プロ月額 $40 | 2.5km (アイコン) | 5 日間 (無料)、16 日間 (プロ) | 制限あり、ET0なし | 有料履歴のみ | REST JSON、Python SDK |
| 明日.io | 1 日あたり 500 回の通話が無料。月額 199 ドルのコア | 1km | 21日 | 土壌水分、ET0、噴霧条件、害虫リスク | 6年 | REST、WebSocket、gRPC |
| ARPA ロンバルディア / ピエモンテ / ヴェネト | 無料(地域オープンデータ) | 定時駅 (イタリア国内に最大 1500 駅のネットワーク) | いいえ (観察のみ) | すべての農業気象変数 | 完全なアーカイブ (数十年) | WMS、REST、CSV/JSONダウンロード |
| メテオアム / CFS/ECMWF | ECMWF の無料オープンデータ。無料の NOAA CFS | 9~25km | 10~45日 | 気象ベースの変数 | 1940年からERA5経由で現在 | ECMWF API (Python)、CDS API |
| ビジュアルクロッシング | 1 日あたり 1000 件のレコードが無料。月額 35 ドルの標準 | ステーションから補間 | 15日 | ET0、成長度日数、霜の危険性 | 無制限の有料 | REST JSON、CSV、SQL ライク |
実稼働環境での AgriTech アプリケーションの場合は、次の使用をお勧めします。 オープンウェザー どうやって ベース (無料、オープンソース、歴史家向けの ERA5 データ、16 日間の予測) と統合します。 ローカル校正用の ARPA ステーションのネットワーク。 Open-Meteo には農業気象変数が含まれています ET0 (基準蒸発散量)、複数の深さの土壌水分、蒸気などの重要な問題 水性ですが、商用プロバイダーから常に無料で入手できるわけではありません。
キャッシュと再試行を備えた Open-Meteo クライアントの実装
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")
完全な地理空間パイプライン: GeoPandas、rasterio、PostGIS
専門的な衛星データ管理には包括的な地理空間インフラストラクチャが必要です ラスター データ (NDVI 衛星画像、天気図) とベクター データ (境界線) を統合します。 区画、管理区域、サンプリングポイント)。ここで説明するパイプラインでは、 地理空間の事実上の標準 Python スタック: ラステリオ ラスターの場合、 ジオパンダ ベクトルの場合、 ポストGIS 空間データベースとして GDAL 基礎となるフォーマット変換エンジンとして。
地理空間パイプラインのアーキテクチャ
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 │
└─────────────────────────────────────────────────────────────────────────┘
ゾーン統計: rasterstats を使用したプロットごとの平均 NDVI
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 クエリ
-- 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;
ML による水ストレス予測モデル
衛星 NDVI と気象データおよび IoT データを組み合わせることで、予測モデルを構築できます 3 ~ 7 日前に水ストレスを把握できるため、事前の灌漑計画が可能になります。 ここで説明するモデルは、Sentinel-2 の NDVI 時系列、 ET0 と降水量から計算された水分バランス、および IoT センサーからの土壌水分 (または 物理センサーが利用できない場合は、Open-Meteo ERA5-Land を使用します)。
モデルの特徴エンジニアリング
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"]]
ケーススタディ: マンドゥリアのプリミティーヴォブドウ園での NDVI モニタリング
説明した概念を具体的にするために、アプリケーションからインスピレーションを得たケーススタディを紹介します。 DOC Primitivo di Manduria (プーリア州ターラント) の領土内にあります。エリアには条件があります 地中海のブドウ栽培の典型的な気候: 暑くて乾燥した夏 (夏の降水量 < 50 mm)、 粘土質石灰岩土壌、カウンターエスパリエ、苗木の育成システム。
ケーススタディのパラメータ - Vigna Primitivo Manduria
| パラメータ | 価値 |
|---|---|
| 位置 | マンドゥーリア (TA)、プーリア州 - 緯度 40.38、経度 17.47 |
| 総表面積 | 35ヘクタールを12の区画に分割 |
| 主な品種 | プリミティーボ (80%)、ネグロアマーロ (20%) |
| 繁殖システム | プーリアの苗木(8本)、カウンターエスパリエ(4本) |
| 植物年 | 2003-2015 (若いブドウと成熟したブドウの混合畑) |
| 灌漑 | 一滴ずつ (すべてのプロットではありません) |
| 追跡衛星 | Sentinel-2A/B、L2A、解像度 10m |
| 使用された取得の頻度 | 2024 年の田舎の 28 シーン (雲 < 20%) |
| 統合型IoTセンサー | 気象観測所 6 台、土壌湿度センサー 18 台 (30/60 cm) |
| 分析期間 | 2024年4月~10月 |
季節ごとの NDVI 結果と灌漑に関する決定
季節ごとのNDVI分析により、著しく異なる水ストレスパターンが明らかになった 灌漑管理と製品の品質に直接影響する区画間:
プロット別の NDVI 結果 - 2024 年キャンペーン
| プロット | バラエティ | NDVI 4 月 | NDVI ダウン | NDVI 7 月 (ピーク) | NDVI 8月 | NDVIセット | 灌漑に関する警報 |
|---|---|---|---|---|---|---|---|
| パークA1 | 原始的/苗木 | 0.28 | 0.42 | 0.58 | 0.41 | 0.31 | 適度な針ストレス |
| パークA2 | 原始的/苗木 | 0.22 | 0.36 | 0.47 | 0.29 | 0.24 | 針をセットする際の重度のストレス |
| パークB1 | プリミティブ/カウンタースパリア | 0.31 | 0.53 | 0.66 | 0.57 | 0.44 | アラートはありません |
| パークB2 | プリミティブ/カウンタースパリア | 0.29 | 0.49 | 0.63 | 0.52 | 0.40 | アラートはありません |
| パークC1 | ネグロアマーロ/木 | 0.25 | 0.44 | 0.55 | 0.38 | 0.28 | 適度な針ストレス |
| パークC2 | ネグロアマーロ/カウンタースパリエラ | 0.30 | 0.51 | 0.61 | 0.50 | 0.38 | アラートはありません |
現れたパターンは明確であり、農学的に重要です。 緊急灌漑は、8 月から 9 月にかけて、例年に比べてより顕著な水ストレスを示します。 ドリップシステムを備えたカウンターエスパリエプロット。特にプロット A2 には、 は、2024 年 7 月から 8 月にかけて NDVI が 0.47 から 0.29 に低下したことを示し、次の点と時間的に相関しています。 降水がなく、最高気温が36度を超える日が23日間続きました。
自動アラートスクリプト
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
衛星監視の投資収益率を定量化することは、 ワイン生産会社への技術投資を正当化する。マンドゥリアのブドウ畑については、 2024 年のキャンペーンの費用対効果分析では次のことがわかります。
ROI 分析 - 衛星監視システム Vigna Manduria (35 ヘクタール、2024 年キャンペーン)
| Voce | 価値 | 注意事項 |
|---|---|---|
| 費用 | ||
| Sentinel-2 データ アクセス (CDSE) | 0ユーロ | 無料のオープンデータ |
| オープンウェザー天気 API | 0ユーロ | 非営利無料 |
| Pythonパイプライン開発 | 2,500ユーロ | 30 開発時間 x 83 ユーロ/時間 (1 回限り、3 年間で償却) |
| クラウドホスティング (VM + ストレージ) | 600ユーロ/年 | VM 4 vCPU + 50 GB S3 ストレージ |
| 統合農業コンサルティング | 1,000ユーロ | 一回限りの作物固有のしきい値キャリブレーション |
| 1年目の総コスト | 4,100ユーロ | 以降の年: ~1,400 ユーロ/年 |
| 推定される便益 (35 ヘクタール、年間 100,000 本) | ||
| 灌漑の削減 (タイミングの最適化) | 2,800ユーロ | 15% 節水、8,500 m3 x 0.33 EUR/m3 |
| 植物検疫処理の削減 | 1,750ユーロ | リスク領域のみに的を絞った介入 (治療量の 20% 削減) |
| ブドウの品質向上 | 8,500ユーロ | 12 ヘクタールの平均糖度 +0.5 ポイント、品質プレミアム +0.20 ユーロ/kg |
| 水ストレスによる損失の削減 | 3,200ユーロ | 予想される厳しいストレス下で 2 つのプロットで 8% の収量回復 |
| 1 年目の合計給付金 | 16,250ユーロ | |
| 1年目のROI | 296% | 約3ヶ月で返済可能 |
Google Earth Engine: 大規模分析のためのクラウド処理
地域または国家規模での分析 (例: レベルでの品種マッピング) DOC 地区、州規模での春の霜による被害の評価、モニタリング AGEA などの支払い機関向けの CAP)、Google Earth Engine (GEE) がコンピューティング能力を提供 クラウドネイティブでは専用の 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}")
イタリアのアグリテックに関する規制とオープンデータ
イタリアの農業オープンデータ エコシステムは、次の 3 つの主要な主体を中心に構造化されています。 アゲア (農業支払庁)、 シアン (全国農業情報システム) と接続された地域ポータル。 2025 年 2 月には、 AGEAは国家戦略極へ向けたSIANの移行を完了し、最新情報を更新した クラウドインフラストラクチャと地理空間データ処理機能の拡大。
イタリアの主な農業オープンデータソース
| ソース | URL | 利用可能なデータ | 形式 | アップデート |
|---|---|---|---|---|
| SIANオープンデータ | データ.sian.it | 高解像度オルソ写真、土地利用地図 (AI ベース)、企業グラフィック プラン | GeoTIFF、シェープファイル、WFS | 年間 |
| AGEA EO サービス | シアン.it | マルチスペクトルオルソフォト、Xバンドレーダー画像、DEM、DSM | ECW、GeoTIFF | 定期(キャンペーン) |
| ISTAT - 農業統計 | istat.it/アグリコルトゥーラ | UAA、生産、企業構造、農業センサス | CSV、JSON、SDMX | 年間 |
| MATTM ナショナル ジオポータル | geoportale.gov.it | カルタ ナチュラ、コリーヌ土地被覆、IDT、DBT | WMS、WFS、シェープファイル | 変数 |
| ARPA 地域 (15 地域) | ハープ.[地域].it | 気象気候観測所データ、地域予測モデル | CSV、JSON、NetCDF | ほぼリアルタイム |
| コペルニクス ランド サービス | land.copernicus.eu | CORINE 土地被覆、草原、水と湿り気、都市アトラス | GeoTIFF、シェープファイル | トリエンナーレ/年次 |
アグリテック向け PNRR 税額控除移行 5.0
衛星監視システムと精密農業への投資により、 移行計画 5.0 (法律 207/2024 - 2025 予算法) の恩恵を受ける資格があります。 2025年から認められる農業ビジネスは投資に対する税額控除を利用できる 有形および無形の資本財 4.0 には、 人工知能コンポーネントと地理空間データ分析。
農業事業に対する移行 5.0 税額控除率 (2025 年)
| 投資範囲 | エネルギー節約 3 ~ 6% | 6 ~ 10% の節約 | 10% 以上の節約 |
|---|---|---|---|
| 最大250万ユーロ | 35% | 40% | 45% |
| 250万ユーロから1,000万ユーロまで | 15% | 20% | 25% |
| 1,000万ユーロから5,000万ユーロまで | 5% | 10% | 15% |
NDVI と天気予報に基づく高精度灌漑システムは適格となる可能性があります 「省エネ >6%」の帯域については、水の消費量の削減が証明されている場合 (ポンピングエネルギー)ベースラインと比較。宣誓された技術的専門知識の要求。
実稼働アーキテクチャ: ベスト プラクティスとアンチパターン
コンポーネントを説明した後は、コンポーネントのベスト プラクティスを要約することが重要です。 本番環境での堅牢な実装を行い、一般的なアンチパターンを特定します。 システムが脆弱または不正確になる可能性があります。
ベストプラクティス - AgriTech 衛星監視システム
- プログレッシブクラウドフィルタリング: 製品の雲量率 (レベル 1 メタデータ) だけを使用しないでください。 QA60 (Sentinel-2) または BQA (Landsat) を使用して、常にピクセル レベルでマスキングを適用します。 5% クラウド製品では、特定のプロットが完全に隠れている場合があります。
- クラウドに最適化された GeoTIFF (COG): S3/MinIO から効率的に範囲リクエストにアクセスできるように、NDVI ラスターを COG 形式で保存します。クラウド ストレージ用に最適化されていない従来の GeoTIFF は避けてください。
- 一貫した投影: 最終データの場合はすべてを EPSG:4326 (WGS84) に、イタリアの場合は UTM ゾーン 32N (EPSG:32632) に標準化して、正しいメトリック測定を維持します。明示的な再投影なしに CRS を混合しないでください。
- IoT グラウンド トゥルースによるデータ検証: 現場測定(土壌湿度センサー、ポータブル LAI メーター、収集された農業データ)を使用して NDVI しきい値を校正します。ローカルキャリブレーションを行わない NDVI だけでは、15 ~ 30% の偽陽性/陰性アラートが発生する可能性があります。
- タイムアーカイブ管理: プロットごとに少なくとも 3 ~ 5 年の NDVI 時系列を維持します。過去の平均 (同じ月、同じ作物) と比較した NDVI の異常は、絶対値よりもはるかに有益です。
- レート制限とキャッシュ API: CDSE と Open-Meteo にはレート制限があります。不必要な再ダウンロードを避けるために、常にローカル キャッシュ (ファイルベースまたは Redis) を実装してください。 Sentinel-2 のダウンロードには 30 ~ 120 秒かかります。同じデータを 2 回リクエストしないでください。
- ロギングと監査証跡: 各 NDVI 計算は、ソース衛星シーン、取得日、雲量率、アルゴリズム バージョンまで追跡可能である必要があります。 PAC/AGEA 監査および異常のデバッグの基本。
一般的なアンチパターン - AgriTech サテライト システム
- 大気補正を無視した場合: レベル 2A (大気の底部) の代わりにレベル 1C (大気の上部) を使用すると、大気湿度が高い条件 (ポー渓谷の霧、アドリア海) では NDVI に 10 ~ 20% 程度の系統誤差が生じます。定量的なアプリケーションには常に L2A を使用してください。
- 唯一の指標としての NDVI: EVI (密集した植生)、NDWI (水ストレス)、レッドエッジ NDVI (NDVI に表示される前の初期ストレス) を無視し、NDVI のみに依存します。プロフェッショナル システムでは、少なくとも 3 ~ 4 つのインデックスを組み合わせて使用します。
- ピクセルサイズとパーセル: 10m のピクセルで 0.3 ヘクタール未満のプロットに NDVI を適用すると、統計的に信頼できない平均値 (30 ピクセル未満) が生成されます。小規模な敷地の場合は、Planet Labs (3.7m) またはドローンセンサーを検討してください。
- 時間的合成の欠如: ある日付で利用可能な単一の取得を取得すると、雲量が 19% であっても、部分的な雲領域で NDVI アーティファクトが発生します。 medoid 合成は常に 10 ~ 15 日の時間枠で使用してください。
- データのバージョン管理を行わないパイプライン: 異なるバージョンの Python/rasterio/sentinelhub を使用して NDVI を再計算すると、同じソース データに対してわずかに異なる値が生成されます。 DVC または同等のものを使用して、コンピューティング パイプラインのバージョンを設定します。
完全なパイプライン: Prefect とのオーケストレーション
説明したすべてのコンポーネントを統合することで、完全な衛星追跡パイプラインが完成します。 Prefect を使用して調整できます (シリーズを参照) パイプライン オーケストレーション) 5 日ごとに自動的に実行されます (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()
結論と次のステップ
CDSE 経由で無料でアクセスできる Sentinel-2 衛星データは、 おそらくイタリアのアグリテックで最も活用されていないデータソースです。パイソンでは、 rasterio、sentinelhub、CDSE アカウントを使用して、 品質と粒度において多くのソリューションを上回る NDVI 監視システム 有料でコマーシャル。
成功の鍵は技術の洗練ではなく、 校正 地元の農学: NDVI しきい値、アラート時間枠、重み 予測モデルの特徴の一部は特定の作物に合わせて調整する必要があります。 地元の気候と同社の農業慣行について。校正されたシステム プーリア州のブドウの木では、それなしではポー渓谷の小麦でも同じようにうまく機能しません。 人間による校正介入。
同様のシステムを実装したい人が覚えておくべき重要なポイントは次のとおりです。
- 常に使用する センチネル-2 レベル-2A (大気補正を適用) 異なる日付間の定量的な比較に使用します。
- NDVI を少なくとも 1 つの補完指標と組み合わせます (密集したブドウ畑には EVI、初期の水ストレスには NDWI、クロロフィル モニタリングにはレッド エッジ NDVI)。
- 統合する オープンウェザー ERA5-ランド 水バランスについては、モデルからの ET0 + 降水量 + 土壌水分、および現場に IoT センサーがない場合の優れた代替手段となります。
- 1 つ保管してください 少なくとも 3 年以上の歴史のあるシリーズ 各区画について: 過去の平均と比較した異常の方が、絶対的な NDVI 値よりもはるかに役立ちます。
- スケール調整 Google Earth エンジン 地域規模または全国規模で処理する必要がある場合にのみ、個々の企業にとっては、CDSE 上のローカル Python パイプラインの方が制御しやすく、ロックインもありません。
フードテックシリーズは続く
衛星パイプラインを構築しました。このシリーズの他の記事をご覧ください。 アグリテック アーキテクチャを完成させます。
- 第1条: Python と MQTT を使用した精密農業用の IoT パイプライン - 地上センサーを衛星パイプラインと統合する方法。
- 第2条: 作物監視のための ML Edge: 圃場でのコンピュータ ビジョン - 植物の病気を早期に検出するためのコンピューター ビジョン モデル。
- 第4条: 食品におけるブロックチェーンのトレーサビリティ: 現場からスーパーマーケットまで - NDVI データがオンチェーンの品質証明書をどのように強化するか。
この記事で説明されている MLOps と予測モデルの展開について詳しくは、 シリーズを見る MLflow を使用したビジネス向け MLOps そしてシリーズ ビジネスにおける LLM: RAG エンタープライズ.







