衛星画像NDVIを無料・自動取得する――Planetary Computer × Python実践ガイド

「特定エリアの植生変化を衛星データで追いたいが、NASAのダウンロードサイトを手作業でたどるのは辛い」――GIS実務でよく聞く悩みです。本記事では Microsoft Planetary Computer の STAC API を Python から叩いて、Landsat／Sentinel-2 の NDVI（正規化差植生指数）GeoTIFF を完全自動生成するスクリプト landsat_ndvi.py を解説します。認証不要・無償、config.ini を1行書き換えるだけで複数年度の差分比較まで行えます。

なぜ Planetary Computer なのか
STAC API の検索モデル：3つの軸で絞り込む
スクリプト全体構成
衛星別バンドと反射率変換の自動判定
NDVI の算出とノイズ処理
2年度差分（ΔNDVI）と伐採検知マスク
config.ini の設定と実行方法
タイムアウト対策：StacApiIO へのリトライ注入
出力ファイルと後続パイプラインへの接続
インストールと動作確認
まとめ

なぜ Planetary Computer なのか

地表反射率データを手に入れる手段は複数ありますが、実務で使い続けるには「STAC API で機械的に検索できる」「COG（Cloud Optimized GeoTIFF）形式で配信される」「認証トークンが自動付与される」の3点が揃っていることが重要です。

Planetary Computer はこれをすべて満たしています。planetary_computer.sign_inplace を渡すだけで SAS Token が自動注入され、アカウント登録すら不要でスクリプトが動きます。Landsat は 1972 年以降、Sentinel-2 は 2015 年以降のアーカイブが対象で、森林・農地の長期変動解析にも対応できます。

STAC API の検索モデル：3つの軸で絞り込む

STAC（SpatioTemporal Asset Catalog）の検索は本質的にシンプルです。「どこ（bbox / ポリゴン）」「いつ（datetime）」「何のデータ（collection）」の3軸を組み合わせて Item のリストを返してもらい、Item が持つ Asset URL から実データを取得するだけです。

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=[137.5, 34.5, 138.5, 35.2],   # 浜松周辺
    datetime="2024-05-01/2024-06-30",
    query={"eo:cloud_cover": {"lt": 50}},
    sortby="+eo:cloud_cover",
    max_items=5,
)

sortby="+eo:cloud_cover" で雲量の少ないシーンを先頭に並べ替えているのがポイントです。先頭1件を採用することで、その期間の「最良シーン」を自動選択できます。

スクリプト全体構成

処理は大きく4ステップに分かれています。

ステップ	処理内容
[1/4] 入力読み込み	GeoJSON または KML → Shapely ジオメトリ → bbox 生成
[2/4] STAC 接続	Planetary Computer に接続、StacApiIO でリトライ設定
[3/4] シーン確認	全年度の候補シーンを事前チェック。1件でも欠ければ中断
[4/4] DL・NDVI算出	Red・NIR バンドをダウンロード → 反射率変換 → NDVI 保存

事前チェックを分離しているのは、「2024年は取得できたが 2025年が0件だった」という状況を途中で気づかず時間を無駄にしないためです。全年度の揃いを確認してからダウンロードに進みます。

衛星別バンドと反射率変換の自動判定

Landsat のバンド番号はシリーズごとに異なります。シーン ID の接頭辞（LC08 LC09 / LT04 LT05）から自動判定するため、ユーザーが意識する必要はありません。

衛星	接頭辞	赤バンド	NIR バンド	解像度
Landsat 9 / 8	LC09 / LC08	SR_B4	SR_B5	30m
Landsat 5 / 4	LT05 / LT04	SR_B3	SR_B4	30m
Landsat 7	LE07	—	—	除外
Sentinel-2A/2B	S2A / S2B	B04	B08	10m

Landsat 7 を除外しているのは SLC-off 問題（2003 年以降、縞状のデータ欠損）のためです。伐採検知などの差分解析では誤検知の原因になるため、シーン ID が LE07 で始まる場合は自動的にスキップします。

反射率への変換式は衛星ごとに固定です。

# Landsat Collection 2 L2（USGS 公式スケールファクター）
reflectance = DN * 0.0000275 - 0.2   # 0〜1 にクリップ

# Sentinel-2 L2A
reflectance = DN * 0.0001             # 0〜1 にクリップ

NDVI の算出とノイズ処理

反射率に変換した Red・NIR バンドから NDVI を算出します。分母がゼロになる画素（水面の暗い画素など）は NaN を割り当てて nodata 扱いにします。

with np.errstate(invalid="ignore", divide="ignore"):
    ndvi = np.where(
        (nir + red) == 0,
        np.nan,
        (nir - red) / (nir + red)
    )

NDVI の値域は −1.0〜+1.0 で、健全な植生は +0.6 以上、裸地・市街地は 0 付近、水面は負値になります。出力は Float32 GeoTIFF で、入力シーンの CRS（UTM 帯）をそのまま継承します。後続の PostGIS インポートや QGIS での重ね合わせをそのまま行えます。

2年度差分（ΔNDVI）と伐採検知マスク

2年度分のシーンが揃った場合、スクリプトは自動的に ΔNDVI を算出します。

# ΔNDVI = NDVI(後年) − NDVI(前年)
# 負値が大きいほど植生が減少（伐採・枯死の可能性）
ndvi_diff = ndvi_2025 - ndvi_2024

# 伐採マスク：ΔNDVI が閾値を下回る画素を 1 に
DEFOR_THRESHOLD = -0.2
defor_mask = np.where(ndvi_diff < DEFOR_THRESHOLD, 1.0, 0.0)

閾値 −0.2 は経験的な値です。森林の伐採では ΔNDVI が −0.3〜−0.6 に達することが多く、−0.2 でほぼ確実に拾えます。草地の季節変動（±0.1 程度）との分離も可能です。浜松・天竜の人工林を対象に試行したところ、視覚的に把握できる伐採跡地との対応が良好でした。

グリッドサイズが年度間でわずかに異なる場合（Sentinel-2 の MGRS タイル境界など）は rasterio.warp.reproject でバイリニアリサンプリングして自動整合します。

config.ini の設定と実行方法

動作パラメータはすべて config.ini に集約されています。スクリプト本体を触らずに対象・期間・衛星を切り替えられます。

[DEFAULT]
input_file   = ../Input/Target_Range.geojson
collection   = landsat-c2-l2
years        = 2024,2025
month_from   = 05
month_to     = 06
max_cloud    = 50
out_dir      = ../Output

years をカンマ区切りで3件以上指定することも可能で、その場合は最若年と2番目の年度ペアで pipeline.ini が生成されます。collection を sentinel-2-l2a に変えるだけで 10m 解像度の Sentinel-2 に切り替わります。

Windows 環境では auto.bat をダブルクリックするだけで実行できます。

@echo off
cd /d "%~dp0Python_src"
python landsat_ndvi.py
pause

タイムアウト対策：StacApiIO へのリトライ注入

Planetary Computer の STAC API は時間帯によってタイムアウト（504）を返すことがあります。デフォルトの pystac_client.Client.open() はリトライ機構を持たないため、1回の失敗でスクリプトが止まります。以下の設定でリトライを注入してください。

from pystac_client.stac_api_io import StacApiIO
from urllib3.util.retry import Retry

retry = Retry(
    total=5,
    backoff_factor=2,                                  # 2→4→8→16→32秒
    status_forcelist=[408, 429, 500, 502, 503, 504],
    allowed_methods=["GET", "POST"],
)
stac_io = StacApiIO(max_retries=retry, timeout=180)

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
    stac_io=stac_io,
)

広域 bbox や複数年度を連続して投げる場合は、検索ループ自体にも自前リトライを加えておくと安心です。max_items=5 に絞ることで1リクエストあたりの負荷も下がります。

出力ファイルと後続パイプラインへの接続

スクリプトが生成するファイルは次のとおりです。

ファイル	内容
`{シーンID}_NDVI.TIF`	年度ごとの NDVI ラスター（Float32、nodata=NaN）
`ndvi_diff.tif`	ΔNDVI（後年 − 前年）
`deforestation_mask.tif`	伐採推定マスク（1=伐採疑い、0=変化なし）
`pipeline.ini`	後続処理（伐採抽出・林分クリップ）用設定ファイル

pipeline.ini は後続の伐採抽出スクリプトが読み込む設定ファイルで、取得したバンド TIF のパスや対象年度が自動で書き込まれます。NDVI 生成から PostGIS インポート・MapLibre 可視化までをひとつのパイプラインとして連結するための接着剤です。

インストールと動作確認

Python 3.9 以上が前提です。3.13 での動作確認済みです。

pip install pystac-client planetary-computer requests numpy rasterio shapely pyproj tqdm

パッケージを入れたら、まず max_cloud=80・max_items=1 で試し取りして接続と出力を確認してください。大きな TIF（赤バンド1ファイルで 200〜400 MB）のダウンロードが走るため、最初の実行は時間がかかります。再実行時は既存ファイルをスキップするので、中断後のリカバリーも安全です。

まとめ

Planetary Computer の STAC API は、衛星データを「URL を知っていれば誰でも取れる構造化されたカタログ」として公開しています。pystac-client との組み合わせで、GeoJSON を渡すだけで複数年度の NDVI と伐採マスクが出力される自動化パイプラインを構築できます。

#NDVI #衛星画像 #PlanetaryComputer #Python #STACAPI #Landsat #Sentinel2 #GeoTIFF #QGIS #PostGIS #GIS #植生解析 #リモートセンシング #森林解析 #伐採検知