最近、衛星データを扱う機会があった為、忘備録も兼ねてデータアクセスから可視化までをご紹介したいと思います。
また、本記事では利用しているTellus APIを含むTellusプラットフォームの紹介と衛星データに対する簡単な解析も行います。
Tellus（テルース）は、政府衛星データを利用した新たなビジネスマーケットプレイスを創出することを目的とした、日本初のオープン＆フリーな衛星データプラットフォームです。複数のデータをかけ合わせ、新たなビジネス創出を促進するためのあらゆるファンクションを提供します。
https://www.tellusxdp.com/ja/about
今回利用するTellus APIはTellusプラットフォームの内の一つです。
これまで衛星データを利用するとなると、特殊なデータフォーマットを扱い、高い処理能力を有したマシンが必要とされる面から産業に於いては限定的な利用状況でした。
これらの課題を解決する為に分析・アプリケーション開発などを行うクラウド環境を提供しているプラットフォームがTellusプラットフォームになります。
一口に衛星データといっても様々なセンサーデータが存在しています。
代表的なデータを紹介します。
光学衛星データは、可視光および近赤外線を使用して地表の画像を撮影します。特徴は以下の通りです。
代表的な光学衛星：Landsatシリーズ、Sentinel-2。
合成開口レーダー（SAR）データは、マイクロ波を利用して観測します。特徴は以下の通りです。
代表的なSAR衛星：Sentinel-1、ALOS-2。
ハイパースペクトルデータは、広範囲の波長にわたる詳細なスペクトル情報を提供します。特徴は以下の通りです。
代表的なハイパースペクトル衛星：Hyperion、EnMAP。
マルチスペクトルデータは、複数の波長帯域にわたるデータを取得します。特徴は以下の通りです。
代表的なマルチスペクトル衛星：Landsatシリーズ、Sentinel-2。
今回は光学衛星データを利用して可視化まで行います。
まずは、Tellus APIを利用する為にTellus Travelerに登録して、API tokenを発行する必要があります。
以下のURLにアクセスして新規アカウント作成を行います。
基本的には画面の案内に従って登録を行えば問題ないかと思います。
必要なものとしては、登録用のメールアドレスと2要素認証用の携帯電話番号が必要です。
アカウント作成が完了したら、以下のURLにアクセスするとアカウント設定画面にアクセスできます。
https://www.tellusxdp.com/account/setting/
「APIトークン」を押下するとAPIトークンの管理画面に移ります。
「トークンの発行」を押下してトークン名を決めるだけでAPIトークンが発行されます。
これで衛星データを利用する準備が整いました。
いよいよ衛星データを取得して可視化してみます。
今回は、Python(Jupyter notebook)を利用して行っていきます。
必要なライブラリをインストールしたいと思います。
APIを利用するのですが、今回はPythonクライアントを利用してデータの取得を行います。
また、衛星データを利用する際に地理データを参照する必要もある為、GISデータを扱えるGeoPandasもインストールします。
以下のコマンドにてインストールします。
python -m pip install geopandas tellus-traveler rioxarray matplotlibAPIコール時に緯度経度を利用して検索をかける為、GISデータを扱う必要がある為インストールします。
Tellus APIを利用する際のPythonクライアント
衛星データはGeoTiff形式で提供される為、画像処理する際に必要
最終的に可視化する為に利用
まずは、取得したAPIトークンをクライアントに設定します。
import tellus_traveler
# Tellus アクセストークンを取得
tellus_traveler.api_token = "[Your API Token]"Tellus APIでは、衛星データはデータセットと呼ばれる衛星センサー毎にまとめられた単位で管理されており、その中の特定エリアを観測した結果(光学センサなら画像)をシーンと呼びます。
今回は任意の場所のシーンを検索して表示するところまで行います。
まずは、どのようなデータセットが存在するのかデータセット一覧を取得します。
# 利用可能なデータセットを取得
datasets = tellus_traveler.datasets()
len(datasets)以下のような出力が得られました。
23どうやら23件のデータセットが利用可能なようです。
続いて取得したデータセット情報を見ていきましょう。
datasets[0]以下得られた出力です。
{'id': '1a41a4b1-4594-431f-95fb-82f9bdc35d6b',
'provider': {'name': 'テルース', 'description': 'テルースが提供する公式データです'},
'tags': [],
'published_at': '2021-10-08T14:24:12.960797+09:00',
'can_order_access_right': True,
'can_order_cut_data': False,
'is_order_required': False,
'minimum_purchase_square_kilometer': 0,
'square_kilometer_per_price': 0,
'related_site': 'https://www.eorc.jaxa.jp/ALOS/a/jp/index_j.htm',
'copyright': '©JAXA',
'properties': ['sat:orbit_state',
'sar:observation_direction',
'view:off_nadir',
'sar:polarizations',
'sar:frequency_band',
'sat:relative_orbit',
'tellus:sat_frame',
'sar:instrument_mode',
'processing:level',
'sar:product_type',
'gsd',
'palsar2:beam'],
'prices': [],
'name': '【Tellus公式】PALSAR-2_L1.1',
'description': 'JAXAが開発したPALSAR-2というSARセンサのデータです。',
'terms_of_use': '/api/traveler/v1/datasets/1a41a4b1-4594-431f-95fb-82f9bdc35d6b/terms-of-use-url/',
'manual': '/api/traveler/v1/datasets/1a41a4b1-4594-431f-95fb-82f9bdc35d6b/manual-url/',
'permission': {'allow_network_type': 'tellus'}}データセットの情報はJSON形式で返ってきます。
上記のデータセットは、JAXAが開発したPALSAR-2というSARセンサのデータセットの様です。
こちら調べると、JAXAが開発した「だいち2号（ALOS-2）」に搭載されたセンサのデータセットでした。
今回は、光学データを利用したいのでJAXAが開発した「だいち（ALOS）」に搭載されたAVNIR-2という光学センサで観測したデータに絞り込んでいきたいと思います。
# AVNIR-2のデータセットを抽出
avnir2_dataset = next(dataset for dataset in datasets if "AVNIR-2" in dataset["name"])
avnir2_dataset出力は以下の通りです。
{'id': 'ea71ef6e-9569-49fc-be16-ba98d876fb73',
'provider': {'name': 'テルース', 'description': 'テルースが提供する公式データです'},
'tags': [],
'published_at': '2021-10-08T14:26:06.089400+09:00',
'can_order_access_right': True,
'can_order_cut_data': False,
'is_order_required': False,
'minimum_purchase_square_kilometer': 0,
'square_kilometer_per_price': 0,
'related_site': 'https://www.eorc.jaxa.jp/ALOS/jp/alos/a1_about_j.htm',
'copyright': '©JAXA',
'properties': ['sat:orbit_state',
'tellus:pointing_angle',
'tellus:bands',
'sat:relative_orbit',
'tellus:sat_frame',
'processing:level',
'eo:cloud_cover',
'gsd'],
'prices': [],
'name': '【Tellus公式】AVNIR-2_1B1',
'description': '解像度10mの広域撮影を目的とした光学カラー画像です。\nJAXAのAVNIR-2センサデータから生成されています。',
'terms_of_use': '/api/traveler/v1/datasets/ea71ef6e-9569-49fc-be16-ba98d876fb73/terms-of-use-url/',
'manual': '/api/traveler/v1/datasets/ea71ef6e-9569-49fc-be16-ba98d876fb73/manual-url/',
'permission': {'allow_network_type': 'global'}}これでデータセットを特定できました。
つづいて、可視化したいシーンを検索していきます。
試しに東京都中野区のシーンを可視化したいので、まずは中野区の緯度経度を入手します。
国土地理院の行政区域データを
こちらから入手します。
ダウンロードが終わりましたら、zipファイルを解凍して.geojsonファイルを読み込んで中野区を矩形で囲える緯度経度を算出します。
# 国土地理院のデータから中野区の緯度経度を検索
import geopandas as gpd
tokyo_gdf = gpd.read_file('data/N03-20230101_13_GML/N03-23_13_230101.geojson')
nakano_ku_gdf = tokyo_gdf[tokyo_gdf['N03_004'] == '中野区']
nakano_ku_bbox = nakano_ku_gdf.total_bounds
print(nakano_ku_bbox)出力は以下の通りです。
[139.62432764 35.67634792 139.69433114 35.73538167]つづいて入手した緯度経度から指定範囲に該当するシーンを検索します。
検索後は、何件ヒットしたか確認します。
# tellus_travelerを使い検索
search = tellus_traveler.search(
datasets=[avnir2_dataset['id']],
bbox=nakano_ku_bbox,
start_datetime="2011-01-01T00:00:00Z",
end_datetime="2012-01-01T00:00:00Z",
)
search.total()出力は以下の通りです。
33つのシーンが検索できたようです。
では、実際にどの程度の範囲カバーできるシーンなのか確認してみます。
# 各シーンをリストに格納
scenes = list(search.scenes())
# 地図上での入手可能なシーンの位置を確認
# シーンの情報を元にGeoPandasデータフレームを作成
search_results_gdf = gpd.GeoDataFrame.from_features(scenes)
# 地図上での入手可能なシーンの位置を確認
# グラフで可視化
search_results_gdf.set_crs(epsg=4326).explore("tellus:name")出力は以下の通りです。
重なっていて分かりづらいですが、どのシーンでも中野区はカバー出来ている様です。
なので、最新のシーンを選びたいと思います。
# シーンを絞り込み
# 撮影時期が最新のシーンを検索
scenes_dates = [[x.__geo_interface__['properties']['end_datetime'], x['tellus:name']] for x in scenes]
scenes_dates = sorted(scenes_dates, key=lambda x:x[0], reverse=True)
scene = next(scene for scene in scenes if scene['tellus:name'] == scenes_dates[0][1])
scene.propertiesこれで可視化したいシーンを検索できました。
最後に得られたシーンをダウンロードして、指定の緯度経度で矩形に切り取り出力させてみたいと思います。
まずは、シーンのGeoTiffファイルをダウンロードします。
今回利用するGeoTiffはCOG(Cloud Optimized GeoTiff)と呼ばれるクラウドに特化したGeoTiffファイルをダウンロードします。
files = scene.files()
target_file = next(file for file in files if "webcog" in file["name"])
download_dir = 'data'
path = target_file.download(download_dir)続いてダウンロードしたファイルをrioxarrayを利用して読み込み、指定の緯度経度で画像を抜き出して表示します。
import matplotlib.pyplot as plt
import rioxarray
# COGを読み込み
data = rioxarray.open_rasterio(path, masked=True)
# 画像データを指定座標で切り抜き
clipped_data = data.rio.clip_box(*nakano_ku_bbox)
# Matplotlibで可視化
fig, ax = plt.subplots(figsize=(8, 8))
clipped_data.sel(band=[1, 2, 3]).astype("uint8").plot.imshow(ax=ax)
nakano_ku_gdf.plot(ax=ax, color="none")
ax.set_title(f"Nakano-Ku True Color\n{scene['tellus:name']}")
print(avnir2_dataset["copyright"])出力は以下の通りです。
これで、衛星データの可視化ができました。
今回利用した「だいち（ALOS）」に搭載されたAVNIR-2という光学センサですが、観測波長が可視光以外にも近赤外線をとらえることができております。
この近赤外線ですが、植物の葉に強く反射される為、植物の位置の可視化ができます。
先ほどの可視化例では、Bandの各波長に合わせてRGBを設定しましたが、Rの部分をBand4に設定することで植物を赤く可視化できます。
import matplotlib.pyplot as plt
import rioxarray
# COGを読み込み
data = rioxarray.open_rasterio(path, masked=True)
# 画像データを指定座標で切り抜き
clipped_data = data.rio.clip_box(*nakano_ku_bbox)
# Matplotlibで可視化
fig, ax = plt.subplots(figsize=(8, 8))
clipped_data.sel(band=[4, 1, 2]).astype("uint8").plot.imshow(ax=ax)
nakano_ku_gdf.plot(ax=ax, color="none")
ax.set_title(f"Nakano-Ku True Color\n{scene['tellus:name']}")
print(avnir2_dataset["copyright"])出力は以下の通りです。
赤みがかった部分が植物が存在していると考えられる場所です。
今回はTellus APIを利用して衛星データの取得から可視化までの方法を紹介しました。
衛星データもこんなに簡単に扱える時代になっていて驚きました。
次は光学センサ以外のデータも扱ってみたいと思います。
https://tellus-traveler.readthedocs.io/en/latest/