衛星データ取得の全工程を解説：openEO・STAC・xarray・rioxarrayを使いこなす【Sentinel-2データアクセス編】

キリン

対象プログラムについて
本記事は、書籍付属のサンプルプログラム（Apache License 2.0）の技術的な内容を解説することを目的としています。コードの直接転載は行わず、処理の流れと技術概念の理解を促すことを目的としています。
原著リポジトリ: https://github.com/tamanome/satelliteBook

はじめに

衛星データ解析の入門書では、「衛星データをどうやって取得するか」という部分が最初の壁になりがちです。本の解説を読む前に、実際のプログラムが何をやっているのかの大枠を掴んでおくと、本を読んだときの理解が格段に早くなります。

この記事では、第3章1節「衛星画像データを取得する」のノートブックが行っている処理を、なぜそうするのかという理由とともに丁寧に解説します。

本記事で扱う技術

$$
\begin{array}{|l|l|}
\hline
\textbf{ライブラリ / サービス} & \textbf{役割} \\
\hline
\text{CDSE (Copernicus Data Space Ecosystem)} & \text{ESAが提供するSentinel衛星データのプラットフォーム} \\
\hline
\text{openEO} & \text{衛星データ処理のオープンAPI標準} \\
\hline
\text{xarray} & \text{多次元配列の操作（NetCDF対応）} \\
\hline
\text{rioxarray} & \text{xarrayにGIS機能を追加する拡張} \\
\hline
\text{STAC} & \text{衛星データカタログの標準仕様} \\
\hline
\text{pystac-client} & \text{STACカタログへのPythonクライアント} \\
\hline
\text{rasterio} & \text{ラスターデータ（GeoTIFF等）の読み書き} \\
\hline
\text{GDAL} & \text{地理空間データ変換の定番ライブラリ} \\
\hline
\end{array}
$$

プログラム全体の構成（フロー）

このノートブックは大きく2つのルートでSentinel-2データを取得します。

衛星データ取得
├── ルート1: CDSE + openEO（本書メインルート）
│   ├── 1. CDSEバックエンドに接続
│   ├── 2. OIDC認証
│   ├── 3. DataCube構築（検索条件設定）
│   ├── 4. NetCDF形式でダウンロード
│   ├── 5. xarrayで読み込み・操作
│   ├── 6. SCL（雲マスク用）の可視化
│   ├── 7. RGBトゥルーカラー画像の表示
│   └── 8. rioxarrayでGeoTIFF出力
│
└── ルート2: STAC + pystac-client（参考ルート）
    ├── 1. STACカタログへの接続・検索
    ├── 2. GeoDataFrameで結果整理・雲量ソート
    ├── 3. サムネイル画像で目視確認
    ├── 4. COGファイルのダウンロード
    ├── 5. rasterioでGeoTIFF作成・マスク処理
    └── 6. GDALで8bit変換・表示用画像作成

ルート1：CDSE + openEOによるデータ取得

なぜCDSEを使うのか

2023年10月、Sentinel衛星データの旧プラットフォーム「Copernicus Open Access Hub」が閉鎖されました。それに伴い、従来広く使われていた sentinelsat ライブラリによるデータ取得ができなくなりました。

新しいプラットフォームである CDSE（Copernicus Data Space Ecosystem） では、以下のデータが無料で利用可能です（2024年4月現在）。

  • Sentinel-1 / 2 / 3 / 5P / 6（SAR・光学・大気など）

  • Landsat-5 / 7 / 8

  • Copernicus DEM（数値標高モデル）

  • MODIS（Terra / Aqua）

  • その他多数

ユーザー登録は公式ドキュメントに従って行います。メールアドレスの確認が必要な点に注意してください。

openEOとは

openEO（open Earth Observation） は、衛星データの検索・処理・ダウンロードを統一されたAPIで行うための国際標準です。バックエンド（サーバー側の処理環境）を意識せずに、同じコードで複数のプラットフォームにアクセスできるのが特長です。

プログラムでは openeo.connect("openeo.dataspace.copernicus.eu") でCDSEのバックエンドに接続します。

接続・認証の仕組み

バックエンドへの接続

接続後、connection.describe_collection("SENTINEL2_L2A") を呼ぶと、そのコレクション（データセット）のメタデータが返ってきます。これで接続が成功しているかを確認できます。返ってくる情報には「どの地域をカバーしているか」「どんなバンドがあるか」「時間範囲はいつからか」などが含まれます。

OIDC認証

データを実際に取得するためには認証が必要です。プログラムでは connection.authenticate_oidc() を実行します。

OIDC（OpenID Connect） はOAuth2をベースにした認証プロトコルです。実行するとURLが出力されるので、ブラウザでそこにアクセスしてCDSEアカウントでログイン・許可するとトークンが発行されます。一度認証するとトークンがローカルにキャッシュされるため、次回以降は自動的に利用されます。

DataCubeの概念とload_collection

openEOの核心的な概念が DataCube です。

connection.load_collection() を呼ぶ時点ではデータはまだダウンロードされません。代わりに、「どのデータが欲しいか」という指示書（DataCube）が構築されるだけです。実際のデータ転送はダウンロード命令を出した時点で行われます（遅延評価）。

プログラムで設定している検索条件は以下のとおりです。


なぜLevel-2Aを使うのか

Sentinel-2のデータにはLevel-1C（大気補正なし・トップオブアトモスフィア反射率）とLevel-2A（大気補正済み・地表反射率）があります。植生解析・土地被覆分類など多くの応用では、大気の影響を除いたLevel-2Aが標準的に使用されます。

バンドの意味


NetCDFでのダウンロードと日本語パス問題

DataCubeを download() するとNetCDF（.nc）形式のファイルが生成されます。

NetCDF（Network Common Data Form） は気象・海洋・衛星データなどの多次元配列データを格納するための標準フォーマットです。xarrayが最も得意とするフォーマットでもあります。

プログラムにはひとつ重要な工夫があります。

# Windowsで日本語パスへの保存を避けるため、ホームディレクトリに保存
NC_PATH = os.path.join(os.path.expanduser("~"), "s2-amami.nc")

これは netCDF4のC言語ライブラリがWindowsで日本語パスを処理できないというバグへの対処です。os.path.expanduser("~") で英字のみで構成されたホームディレクトリ（例：C:\Users\username）のパスを取得しています。

さらに、ダウンロード時に resample_spatial(resolution=10, projection="EPSG:32652") を挟んでいます。これにはふたつの意図があります。

  1. 解像度の統一：Sentinel-2の各バンドは10m・20m・60mと解像度が異なります。10mに統一することでバンド間のピクセルが一致します。

  2. CRSの明示指定：サーバー側のバグで一部タイルのCRS情報が壊れていることへの回避策です。EPSG:32652 はUTM座標系（Zone 52N）で、日本付近の奄美大島に適した投影座標系です。

xarrayによるデータ操作

ダウンロードしたNetCDFを xarray.open_dataset() で読み込むと、Dataset オブジェクトが得られます。

このDatasetの構造は以下のとおりです。

Dimensions: (t: N, y: M, x: L)

Coordinates:
  * t    (t)    datetime64[ns]   ← 日時（複数シーン分）
  * x    (x)    float64          ← X座標（投影座標）
  * y    (y)    float64          ← Y座標（投影座標）

Data variables:
    crs   → 座標参照系の情報（配列ではない）
    B04   (t, y, x) float32   ← 赤バンド
    B03   (t, y, x) float32   ← 緑バンド
    B02   (t, y, x) float32   ← 青バンド
    SCL   (t, y, x) float32   ← Scene Classification Layer

※ノートブックのマークダウン説明に x: latitude / y: longitude と記載がありますが、投影座標系（UTM）を使用しているため、実際にはx/yは緯度経度ではなくメートル単位の投影座標です。

時間次元でのスライス

ds.sel(t="2024-03-11") のように、日時を文字列で指定することでその日のデータだけを取り出すことができます。これはxarrayがPandasのインデックス操作に似た形でラベルベースの選択をサポートしているためです。

RGBデータへの変換

複数のバンドをひとつの DataArray にまとめるには to_array(dim="bands") を使います。これにより (bands, t, y, x) という4次元配列になります。あとは plot.imshow() で直接描画できます。

表示の際の vmin=0, vmax=2000 は、Sentinel-2 Level-2AのDN値（デジタル数値）のスケールに対応します。生の値は0〜10000程度で、0〜2000の範囲を使うと暗め〜適正露出くらいの描画になります。

SCL（Scene Classification Layer）とは

SCL は各ピクセルが何であるかを分類した補助レイヤーです。プログラムでは ds["SCL"].plot.imshow(col='t') で時系列ごとに並べて表示しています。

雲関連は値3・8・9・10です。後の章で行う解析では、このSCLを使って雲ピクセルをマスクする処理が登場します。

rioxarrayでGeoTIFFに保存

rioxarray は xarray に空間参照情報（CRS・変換行列）を扱う機能を追加した拡張ライブラリです。rxr.open_rasterio() でNetCDFを開き、.rio.to_raster() でGeoTIFF形式に出力できます。

注意点として、open_rasterio() には decode_times=False を指定しています。これはNetCDFの時間次元のデコード方式がxarrayとrioxarrayで異なることへの対処で、decode_times=False で読んだ後に ds["t"] を上書きする形でxarrayで正しくデコードした時間情報を設定しています。

ルート2：STACによるデータ取得

STACとは

STAC（SpatioTemporal Asset Catalog） は衛星データのカタログ情報を標準化したJSONベースの仕様です。データ提供者がSTAC準拠のカタログを公開することで、利用者は統一されたAPIでデータを検索・取得できます。

プログラムでは https://earth-search.aws.element84.com/v1 というAWS上のSTACカタログを使用しています。これはCloud-Optimized GeoTIFF（COG）形式のSentinel-2データをクラウドで提供しています。

検索と結果の整理

pystac_client.Client.open() でカタログに接続したあと、client.search() で以下の条件で検索します。

  • collections：sentinel-2-l2a（Sentinel-2 Level-2A）

  • bbox（Bounding Box）：検索する地理範囲

  • datetime：検索する日付範囲

  • query：{"eo:cloud_cover": {"lt": 30}} → 雲量30%未満

検索結果は GeoDataFrame.from_features() でGeoDataFrameに変換されます。これにより雲量（eo:cloud_cover）などのメタデータを列として扱え、sort_values() で雲の少ない順に並び替えることが容易です。

AOI（Area of Interest）の座標変換

プログラムでは負の経度値（-220°など）が出てきますが、これは360°を超えた表現を正規化する処理が入っています。

# -220は140（東経140°）に変換される
if AREA[i][0] >= 0:
    AREA[i][0] = AREA[i][0] % 360
else:
    AREA[i][0] = -(abs(AREA[i][0]) % 360) + 360

このような処理が入っているのは、データソースによって経度の表現方法が異なる場合があるためです。

Cloud-Optimized GeoTIFF（COG）とは

STAC経由で取得できるデータの多くは COG（Cloud-Optimized GeoTIFF） 形式です。通常のGeoTIFFと違い、COGはファイル全体をダウンロードしなくても必要な部分だけをHTTPリクエストで取得できるように最適化されています。これにより、大きなファイルでも検索対象の小さな領域だけ高速に取得できます。

プログラムで selected_item[0][band].href で取得しているURLは、クラウドストレージ（AWS S3）上のCOGファイルのURLです。rasterioはCOGに対応しており、URLを直接渡すことでクラウド上のファイルをストリーミングで読み込めます。

rasterioでGeoTIFF作成・マスク処理

GeoTIFFの作成

B02（青）・B03（緑）・B04（赤）を別々にダウンロードして1つのGeoTIFFにまとめます。rasterioで新規ファイルを作る際は、以下のメタデータを明示的に指定します。

  • driver：'Gtiff'（GeoTIFF形式）

  • width / height：画像サイズ（B04と同じ）

  • count：バンド数（ここでは3）

  • crs：投影座標系（B04と同じ、epsg:32654）

  • transform：アフィン変換行列（ピクセル座標と地理座標の対応）

  • dtype：データ型（uint16：0〜65535の整数）

CRSが epsg:32654（UTM Zone 54N）を使用しているのは、この検索領域が日本本土付近（関東周辺）のためです。

AOIでのマスク処理

rasterio.mask.mask() で指定したポリゴン範囲だけを切り出します。このとき、画像のCRSとポリゴンのCRSが一致していないとエラーになります。プログラムでは GeoDataFrame に .to_crs(crs='epsg:32654') を適用して座標系を合わせています。

GDALによる8bit変換

Sentinel-2のDN値はuint16（16ビット）ですが、通常のディスプレイ表示には8bit（0〜255）が必要です。GDALの Translate() 関数で変換しています。

-scale 0 255 0 25

このスケールパラメータは「元の0〜255を出力の0〜25に変換する」ことを意味しています。実際のSentinel-2のDN値は0〜10000程度なので、これは非常に明るい設定になっています。用途に応じてスケールを調整する必要があります（例：-scale 0 3000 0 255 など）。

各ルートの比較

本記事のまとめ

次のステップ

このノートブックをマスターしたあとは、続く章で以下を学びます。

  • ch3.2：座標系の扱い（CRSの変換、投影変換）

  • ch3.3：GDALの詳細な使い方

  • ch4：バンド演算（NDVI等）、植生・道路・農地・沿岸の解析

  • ch5：機械学習（線形回帰、SVM）による衛星データ解析

  • ch6：教師なし分類（クラスタリング）

衛星データ取得という「データの入口」をしっかり理解したうえで、次の解析のステップに進むと全体像がつかみやすくなります。

本記事は Apache License 2.0 のもとで公開されているサンプルプログラム（https://github.com/tamanome/satelliteBook）の技術的解説を目的として作成しています。


いいなと思ったら応援しよう！