Microsoft Planetary ComputerのAPIを利用して、衛星データを扱う機会がありましたが、日本語の情報がなく英語を読まなければならなかったため、苦労しました。今後、衛星データを扱う人々が気軽にMicrosoft Planetary Computer APIを利用できるよう、日本語の情報があると便利だと感じ、この記事を書くことにしました。
また、衛星データによる分析は個人的に楽しい経験であり、この分析方法についても共有できればと思い、記事に取り上げました。

衛星データの利点

衛星データを手に入れる手順を説明する前に、衛星データがデータ分析においてなぜ魅力的な対象なのかを説明します。データサイエンスにおいて、分析対象に対するデータが不足していることは最大の問題の一つです。そのため、データ取得コストが安く、大量のデータを取得できることが望ましいです。現在、人工衛星は数多く打ち上げられ、低コストで入手できるだけでなく、リアルタイムにデータを取得できるため、農業や都市開発、海洋観測などに活用されています。このため、衛星データはデータ分析において魅力的な対象となっています。

Microsoft Planetary ComputerのAPIを利用してデータを取得する方法

本題に入ります。まず、Microsoft Planetary Computerからデータを取得するためには、アカウントを作成する必要があります。そのためには、以下のサイトにアクセスしてMicrosoft Planetary Computerのアカウントを作成しましょう。
Microsoft Planetary Computer

このような画面が見えたら右上にあるRequest accessをクリックしましょう。

このリクエストフォームに必要事項を入力し、送信してください。すると、1日か2日程度でPlanetary Computer Previewからメールが届き、アカウントリクエストが承認されたかどうかの通知が送られてきます。
アカウントリクエストが承認されたら、以下のサイトに飛んでログインしてください
Microsoft Planetary Computer APIs
手順に従ってサインインした後、プロファイルを確認してください。画像内のPrimary KeyとSecondary Keyが、あなたのAPIキーとなります。

Pythonで衛星データをダウンロードする手順を解説する前に今回使うデータのドメイン知識について解説しようと思います。

今回使用するデータはSAR（合成開口）データとなります。通常、衛星画像データは多くの観測日において雲に覆われており、地表の状態を確認することが困難です。しかし、SARデータは衛星から電波を送信し、地表面からの反射波を受信するという観測方法を用いるため、天候や昼夜の時間帯に関係なく高精度なデータを取得することができます。
SARデータは、様々なレーダー波長に対する物体の応答を記録したデータです。主に以下の種類の波長が用いられています（カッコ内はそれぞれの波長を利用している人工衛星の名前です）。

• P-band: ~69.0 cm (BIOMASS)
• L-band: ~23.5 cm (ALOS -2 PALSAR-2, SAOCOM-1, NISAR-L)
• S-band: ~9.4 cm (NovaSAR, NISAR-S)
• C-band: ~5.6 cm (Sentinel-1, Radarsat-2, RCM)
• X-band: ~3.1 cm (TerraSAR-X, TanDEM-X, COSMO-SkyMed)

X-bandのような比較的波長が短いマイクロ波は草地や樹冠、樹木の葉表面で散乱が発生するような物理応答を示し、L-bandのような比較的波長が長いマイクロ波は草地や樹冠を透過し、樹幹や地表面で散乱が発生するといった物理応答を示します。この応答の違いを用いることで、衛星データは意味のある数値を持つこととなります。詳しくはデモ分析で解説いたします。
また、レーダーは、水平偏波（進行方向に対して水平な面を持つ：H）と垂直偏波（振幅が進行方向に対して垂直な面を持つ：V）の2つに分類され、それぞれ送受信の組み合わせでHH（水平-水平）、HV（水平-垂直）、VV（垂直-垂直）、VH（垂直-水平）と分類されます。

さらに詳しい情報はこちらをご覧ください
Laymans_SAR_Interpretation_Guide_2.0a
国土技術政策総合研究所合成開口レーダ(SAR)に関する基本事項

次に、衛星データのダウンロード方法について解説します。手元のPythonでもダウンロードは可能ですが、今回はMicrosoft Planetary ComputerのHubを利用してみましょう。下の画像の上部付近にあるHubをクリックしてください。

Planetary Computer Hubは、Google Colabのように必要なライブラリがあらかじめインストールされている実行環境です。今回は、CPUを実行環境として選択します。

以下の画像のような画面になったら成功です。

今回は試しにSentinel 1 Radiometrically Terrain Corrected (RTC)をダウンロードします。Sentinel-1はC-bandの波長のレーダーデータを提供する衛星です。
まずはデータをダウンロードする関数を組みましょう

get_data.py

import warnings
warnings.filterwarnings('ignore')
import ipyleaflet
import matplotlib.pyplot as plt
from IPython.display import Image
import numpy as np
import pandas as pd
#Planetary Computer Tools
import pystac
import pystac_client
import odc
from pystac_client import Client
from pystac.extensions.eo import EOExtension as eo
from odc.stac import stac_load
import planetary_computer as pc
import xarray as xr
pc.settings.set_subscription_key('ここに上で確認したPrimary KeyかSecondary Keyを打ち込んでください')

def get_data(longitude,latitude,time,assets,strides):
    #取ってくる範囲を指定
    bbox = (longitude, latitude, longitude+strides, latitude+strides)
    time = time #興味のある時間範囲を指定
    catalog = pystac_client.Client.open(
        "https://planetarycomputer.microsoft.com/api/stac/v1"
    )
    search = catalog.search(
        collections=["sentinel-1-rtc"], bbox=bbox, datetime=time
    )
    items = list(search.get_all_items())
    resolution = 10  #1ピクセル何メートルか
    scale = resolution / 111320.0 #EPSG:4326に合うよう変換
    data = stac_load(items, patch_url=pc.sign, bbox=bbox, crs="EPSG:4326",bands = assets,
    resolution=scale)
    return data

上のコードの引数の説明をします。　

longitude, latitude: float型、興味のある緯度経度の四隅のうちの一つです。
time: str型、衛星データのうち、興味のある時間範囲を指定します。分析に応じて決定します。
assets: list(str)型、興味のある衛星データを、str型を要素とするリストで指定します。
strides: float型、取得する緯度経度の範囲を指定する引数です。

上のコードの

get_data.py

search = catalog.search(
        collections=["sentinel-1-rtc"], bbox=bbox, datetime=time
    )

の"sentinel-1-rtc"の部分を他のMicrosoft Planetary Computer内のデータセットに変えれば、他のデータの情報を得ることもできます。

pystacに関しては公式のドキュメントをご覧ください
PySTAC

デモ分析

ここまででAPIを利用して衛星からのデータをダウンロードする方法は説明しました。しかし、このデータが何の役に立つのかわからないままでは、単なるデータの収集に終わってしまい、面白みもありません。簡単な探索的データ分析を行い、衛星データの面白さを伝えたいと思います。

上空から稲の栽培地域のVV、VH、VH/VVの特徴を確認しよう

そこで、今回は新潟県内における米の栽培地域と非栽培地域のVV、VH、VH/VVなどを時系列でプロットしてみたいと思います。VV、VHの値は衛星からのレーダーに対する物理的応答を反映しているため、水、葉、木の幹、コンクリート、砂浜など、それぞれの物理的応答に応じて値が変化することが知られています。では、稲はどうでしょうか？簡単なサーベイによると、稲は5月頃に種を蒔き、8月には穂が出て、9月には穂が垂れ下がり、収穫を迎えるとのです。このことから、5月頃はまだ種が水田の底に沈んでいるため、水と同様のVH、VVの応答を示すと考えられます。そして、8月頃には葉が出てくるため、葉と同様のVV、VHの応答が見られるはずです。また、今回使用するSentinel-1はCバンドで波長が比較的短いため、穂が出てきた時に物理的応答が顕著になると考えられます。

データを取得する前に、仮説を確かめるために、Google Earthで実際に新潟県を見てみることで仮説の信頼性を高めましょう。

新潟県の水田、6月と8月の様子


田植え直後の6月に比べて、8月には緑がより生い茂っている様子がわかります。
新潟県の森、6月と8月の様子


常緑樹は季節によって緑の量に大きな変化は見られません。　
新潟県の湖、6月と8月の様子


湖も季節によって変化はなさそうです。

ここまで確認できたところで、データをダウンロードしてみましょう。

#"2022-05-01/2022-09-30"の間のデータを取る
#strides引数は地図をみて決定した。

rice = get_data(138.482403,37.230511,"2022-05-01/2022-09-30",["vv", "vh"],0.0001)
forest = get_data(138.501666,37.212620,"2022-05-01/2022-09-30",["vv", "vh"],0.001)
water = get_data(138.502852,37.210894,"2022-05-01/2022-09-30",["vv", "vh"],0.0001)
#取得した期間は変数にしてしまう
time = rice["time"].compute().data
rice

上のコードを実行すると、出力は以下の通りになります。

get_data関数は、xarray型で取得したデータを返す関数です。xarrayは、VVやVHなどの情報の他に、どの緯度経度で観測されたものか、どの時間で観測されたものかという情報も含むデータ型です。PandasのDataFrame型に時間と座標が追加されたようなものだと考えてください。
非常に便利なライブラリなので衛星画像を扱う際は以下のサイトを参考にすると良いです。
xarray公式ドキュメント
 xarray を用いたデータ解析

次に、観測範囲の平均を取りましょう。

rice_vv = np.mean(np.mean(rice["vv"].compute().data,axis = 1),axis = 1)
rice_vh = np.mean(np.mean(rice["vh"].compute().data,axis = 1),axis = 1)
forest_vv = np.mean(np.mean(forest["vv"].compute().data,axis = 1),axis = 1)
forest_vh = np.mean(np.mean(forest["vh"].compute().data,axis = 1),axis = 1)
water_vv = np.mean(np.mean(water["vv"].compute().data,axis = 1),axis = 1)
water_vh = np.mean(np.mean(water["vh"].compute().data,axis = 1),axis = 1)

さらに、これらのVH/VVも計算します。

rice_vh_by_vv = rice_vh/rice_vv
forest_vh_by_vv = forest_vh/forest_vv
water_vh_by_vv = water_vh/water_vv

これらを描画したものが以下の通りになります。

グラフを見ると、3つのグラフ全てで、田植え期に当たる5月は、WaterとRiceの値が近いことがわかります。これは、田植え期間はWaterと近い物理的応答を示すという仮説に合致していることを示唆しています。
また、3つのグラフ全てでWaterとRiceは10月にはまた近い値になっていることがわかります。これは稲が刈り終わって、また水と同じ物理的特性を示すようになったからであると考えられます。

3つのグラフの特徴を見ていきましょう

VV

全体的にForestが高い値を示しているが、期間を通して大きな変化は見られない。
Waterは期間中に変化が大きく見られる。
Riceは全体的に低い値を示しており、期間を通して変化が小さい。

VH

全体的にForestの値が高いが、変化が小さい。
WaterとRiceは9月に入るまで似たような値の傾向を示している。
9月初め、ForestとWaterが低下している中、Riceのみが上昇している。

VH/VV

全体的にRiceが高く、また分散も大きい
ForestはRiceとWaterに挟まれた位置にある傾向があり、分散が小さい

以上のグラフを見ると、VV、VH、VH/VVの値が物体によって異なることが分かります。衛星データが提供する空中からの観測情報が、ここまでの情報を提供できることは、非常に興味深い点です。
[2023/4/3 15:30追記] 「湖も季節によって変化はなさそうです。」と述べましたが、グラフを見ると明らかに時系列で変動していることがわかります。例えば夏の間は、蒸発などにより水位が下がったり、月の満ち欠けにより水位が変動していることが考えられます。そのため、衛星データはダムの水位予測にも応用できるかもしれません。

より具体的な分類モデルの構築のために

各変数の時系列データから、分散や最大値、最小値などの統計量を抽出して、それらを特徴量とするテーブルデータを作成し、農地か否かの判別に使用する。これらの統計量は、時系列データの情報を要約することができるため、分類モデルの学習に適した形式となります。
一次元畳み込み層などのディープニューラルネットワークを用いて、時系列データを時系列データのまま処理する

このような問題に対して二次元の畳み込みを用いた論文があります。
Paddy Rice Mapping in Thailand Using Time-Series Sentinel-1 Data and Deep Learning Model

Sentinel-2-LAのデータを活用して新たに植生変数を作る

植生変数とは衛星画像を用いた植生の解析のために用いられる指標のことです。この変数を用いて稲の収量を予測した論文もあります。
Predicting rice yield at pixel scale through synthetic use of crop and deep learning models with satellite data in South and North Korea

まとめ

衛星データの解析は、思っていたよりも非常に興味深かったため、今回記事を書かせていただきました。衛星画像に興味を持たれた方は、ぜひ自分なりのアプローチを試してみてください。

*この記事はChatGPTに日本語を添削してもらいながら書きました。もし間違いがあればお知らせください。

Microsoft Planetary ComputerのAPIを利用して衛星からのデータをダウンロードしてデータ分析をする

目次