AWSのオープンデータとPythonを使って、宇宙から地球を覗いてみよう！

初めに

みなさんこんにちは！

宇宙ビジネスってご存知でしょうか？

「宇宙ビジネス」とだけ言われてもあまり想像がつかないかもしれませんが、GPSやロケットの製造や運用、打ち上げや制御のためのシステム・衛星テレビ配信サービスやその運用・それらを受診するためのアンテナや端末など、「宇宙」にまつわるものは案外たくさんあります。

リモートセンシングと呼ばれる衛星画像を取り扱った「地球観測」もその一つですね。

そんな宇宙ビジネスですが、こちらの記事によると2019年での市場規模は40兆円にまで到達しているそうです。

（宙畑より（https://sorabatake.jp/216/ ））

クラウド市場は2021年の段階で1800億ドル弱（当時のドル/円は115円程度だったので、約20兆円）のようなので、およそ2倍！
すごい市場規模ですね。

（クラソルより（https://business.ntt-east.co.jp/content/cloudsolution/column-374.html#section-3 ））

AWSでも色々やってる

AWSでもAWS Ground Stationのような宇宙関連サービスなどは展開されておりますし、先日行われたAWS re:InventでもAWS SageMakerがGeospatial MLに対応するなど、対応が増えてきているところです。

で、僕のようなWeb × GIS（地理空間情報）エンジニアとして気になるのは宇宙に関する「データ」の部分なんですが、衛星データというのは基本的に「とっても高い」です！！！（数百万円とかそんなレベル）

というのも当然で、衛星は打ち上げるだけでも大変（最近は比較的安価に打ち上げられるようになっている）ですが、打ち上げる以上は「宇宙と地球とで」何かしらのデータのやり取りをするのが普通です。

リモートセンシングの分野では衛星から可視光や近赤外線などの様々なデータを受け取り、それを農作物の育成状況の確認などに利用しますが、この際に、衛星データを受け取り、適切に解析するための大きな基地局を地上に作る必要があります。

ただ衛星からのデータのみを利用したい一般の事業者にとってはとても高いハードルとなることがよくわかりますね。

このため課題を解決するのがAWS Ground Stationだったりはするんですが、ただのエンジニアにはなんにせよハードルが高いぐぬぬ。。。

というあなたに朗報です！

なんと！AWSではそんな衛星データをオープンデータとして無料で公開しております！！！！
今回はそのデータとPython(Jupyter notebook)を利用して衛星データの世界をのぞいていきます！

衛星データを触ってみる

Sentinel-2って？

Sentinel-2は欧州宇宙機関と欧州連合が運用を行なっている衛星で、2015年に初号機（2A）が打ち上げられて以来、2022年現在まで4機打ち上げられてます。

(Sentinel-2の外観図(https://www.restec.or.jp/satellite/sentinel-2-a-2-b.html ))

地球一周（というと語弊があるかもしれませんが…）約10日かかるらしく、4機で回っているので、2.5日に一回は最新の衛星データが取得できるわけですね。すごい！

「高速道路の渋滞情報」のようなリアルタイム性があるわけではありませんが、地球全域のデータを2.5日に1度更新できるというのはまぁほぼリアルタイムと言っても過言ではないでしょう。

Sentinel-2の他に無料で利用できる有名な衛星データとしては、Landsat-8/9なんかもあって、微調整は必要ですが、これらも合わせることでさらに高頻度でデータを取得することが可能になりますが、今回はSentinel-2のみ触っていきます。

データを眺めてみる

EO BrowserというWebサイトではSentinelやLandsatの衛星画像を地理的範囲や撮影時期を絞って検索することができます。

今回はこれを利用して自分の欲しい範囲のデータを取得してみましょう！

EO Browserでデータを検索してみる

サイトに移動したら、取得したい範囲に画面を移動しましょう。

今回は札幌市周辺のデータが欲しかったので、その辺りまで地図を移動させました。

画面左側にはメニューが表示されていますが、これをSentinel-2のみチェックを入れ、Advanced searchをオンに変更、L2Aを選択（L1のデータは取り扱いが難しいので、他の衛星を検索する時もL2にするのをお勧めします）します。

雲が被っている画像はあまり利用したくないので、「Max. cloud coverage」は20%くらいまで下げちゃいましょう。

メニュー下部に移動すると、「Time range」というのが選べるので今回は2020年の3月のデータを取得してみましょう！

この状態で「Search」のボタンを押すと検索結果が出てきます！

ちょうど札幌市のデータはありませんでしたが、その周辺のデータで良いとしましょう。

表示された水色の台形をクリックすると以下のようにデータが表示されます。

この状態でリンクのアイコンをクリックすると「AWS path」なるものが表示され、s3://sentinel-s2-l2a/tiles/54/T/VN/2020/3/31/0/のようなS3 URIが表示されたかと思います。

ここに対してAWS CLIなどを通じてGetObjectのAPIにリクエストをすることで、データが無料で取得できる、というわけです。

S3のURIが表示されていることからも分かるとおり、Sentinel-2のデータはオープンデータとして、AWSのS3にホスティングされているんですね。

AWSにあるデータを検索してみる

AWSオープンデータレジストリ

ちょっと話は変わりまして、AWSのオープンデータというのはRegistry of Open Data on AWSというところで管理されています。
https://registry.opendata.aws/

機械学習用のデータがメインストリームな気はしますが、多種多様なデータが置いてあり、全て検索することが可能です。

こちらからSentinel-2のデータを検索できます。

以下のようにテキストボックスに「sentinel-2」と入力するだけで候補がサジェストされます。

先ほどEO Browserで検索したのは一番上の「Sentinel-2」のデータセットの中の衛星画像なんですが、今度は「Sentinel-2 Cloud-Optimized GeoTIFFs」の方を選んでみましょう。

GeoTIFFとは衛星データなどの地理情報が付与されたTIFF画像のことですが、S3などのクラウドサーバーに配置し、HTTP　rangeリクエストで画像の一部分だけ取得することができるようになり、クラウドに最適化されたGeoTIFFのことをCloud-Optimized GeoTIFFと呼びます。

衛星画像はGB、TB級になることはよくあり、一部のデータはPB級になるとかなんとか。

なのでいちいち全てのデータをローカルマシンにダウンロードしてから解析する、というわけにはいきません。

それを解決したのが「Cloud-Optimized GeoTIFF」というわけです。

ちなみに、「COG（シーオージー、コグ）」のように略されることが多いです。

詳しくはこちらを参照してください。
https://www.cogeo.org/in-depth.html

今回は、そもそもデータの範囲を絞って検索したりするので、普通のGeoTIFFでも良いんですが、なんかかっこいいのでこっちを触っていきたいですが、「Registry of Open Data on AWS」では地理的検索や、期間の絞り込みなどはできません。

Earth SearchのAPIを利用する

Earth Searchというサイトがあります。

こちらはAWSオープンデータレジストリに格納されているSentinel-2やLandsat-8などを検索できるAPIになります。

このAPIはSTACと呼ばれるCOGの検索を手軽に行えるようにするためのカタログ仕様で実装されたSTAC APIというもので、仕様が決まっているので、Pythonからシームレスに検索するためのライブラリが用意されていたりします。

（STACの詳しい説明はこちらから）

まずはこちらのインストールを行なっていきましょう。

Pythonで衛星データを扱ってみる

必要なライブラリをインストールする

衛星画像などの地理空間情報の解析にはAnacondaをおすすめしています。

インストールの方法なんかはMacでGISデータ分析を始めるためにサクッとAnacondaとjupyter labをインストールしてみるに記載していたりしますので、参考にしてみてください。

Anacondaのインストールが終わり、仮想環境を立ち上げアクティベートが完了している想定で話を進めていきます。

以下のようにライブラリをインストールしていきましょう。

$ conda install -c conda-forge gdal pystac pandas geopandas boto3 numpy matplotlib shapely
$ python -m pip install pystac_client

地理空間情報界隈では、GDALのインストールが特に鬼門で、Win/Mac/Linuxでそれぞれハマりポイントがあり、誰でも一度は泣いたことがあると思いますが、Anacondaを使えば比較的用意にインストールできるようになっています。

Earth SearchのAPIを触ってみる

それでは実際にPythonでEarth SearchのAPIを利用しデータを検索してみます。

これ以降はjupyterの利用を推奨している他、jupyterで動作するコードを書いていきますが、.pyファイルなどの生のPythonでも大体は動作すると思います。

まずはpystac_clientにAPIを読み込ませていきましょう。

stac_catalog_url = "https://earth-search.aws.element84.com/v0"

# 今回利用するのは、STAC API
# STAC APIを開くときは、Clientを使う
# 静的なSTAC Catalogの場合は、Catalogを使う
# どちらにしろ、STACの入り口である最上位の「カタログ」を取得できる
catalog = Client.open(stac_catalog_url)
print(catalog)

このカタログ内にどんなコレクション（データの集合みたいなもの）が含まれているか見てみましょう。

# 関連アイテムをグループ化し、アイテムの集約やメタデータの提供などを行う「コレクション」
# カタログはトップレベルのコレクション内に、カタログやコレクションなどのネストされたレイヤーを持てる
collections = [(c.id, c.title) for c in catalog.get_collections()]
collections

そうすると以下のような出力が得られるはずです。

[('sentinel-s2-l2a', 'Sentinel 2 L2A'),
 ('sentinel-s2-l1c', 'Sentinel 2 L1C'),
 ('sentinel-s2-l2a-cogs', 'Sentinel 2 L2A COGs'),
 ('landsat-8-l1-c1', 'Landsat-8 L1 Collection-1')]

このAPIからはLandsat-8、Sentinel-2のL1とL2、そしてL2のCOGのデータが取得できることが分かります。

今回はCOGを利用したいので、sentinel-s2-l2a-cogsを指定していきましょう。

# idを指定して、コレクションを取得する
# sentinel-s2-l2aはjp2kファイルで格納されているが、こっちはCOGで配信されている
collection_id = 'sentinel-s2-l2a-cogs'

collection = catalog.get_collection(collection_id)
collection_df = pd.DataFrame.from_dict(collection.to_dict()['item_assets'], orient='index')
collection_df

取得したデータは視認性向上のためにpandasのDataFrameに変換していきます。

以下のようにBandが12個ある（B01~B12）のがわかりました。

今度はデータを検索していきましょう。

検索に最低限必要な情報はEO Browserで検索した時とほとんど同じで、以下になります。

• collectionのid
• 地理的範囲
• 期間
• 雲の量

これらの中で、ちょっとだけ設定が難しい「地理的範囲」の設定をしていきましょう。

geojson.ioでAOIを作る

衛星データを地理的な範囲で絞る場合には、対象エリアを指し示すArea Of Interest(AOI)というものが利用されます。

これを簡単に作成できるgeojson.ioというサイトを利用してみましょう。
（GeoJSONとは、緯度・経度などの地理空間情報を格納するために仕様が拡張されたJSONファイルです。）

サイトにアクセスすると、いい感じの地球儀が表示されると思います。

自分の好きな場所にズームしたり、移動したりしましょう。

今回は札幌市の範囲にズームしてみました。

画面右側のバーに表示されている四角いアイコンをクリックしてみてください。

そうすると地図に四角形が書き込めるようになりますので、指定したい領域を指定します。

左上の端点と、右下の端点をクリックすることで以下のように四角形が作成できます。

この状態で、テキストエディターを開き、新たなファイルを作成、ファイル名を「aoi.geojson」として、作業中のjupyter notebookファイル（.ipynb）と同様のディレクトリに配置してください。

ファイルの中には、geojson.ioの右側に表示されているJSONを書き込みます。

これでAOIの作成が完了しました。

作成したAOIなどを利用して、検索してみる

以下のようにして検索します。

# カタログの基本的な構成である「アイテム」を検索する
# カタログもコレクションもアイテムを持てる
# アイテムは1つの時空間的なリソース（衛星シーンなど）を持てる
aoi = gpd.read_file('aoi.geojson')
geom = aoi.geometry[0]

time_range = "2020-03-01/2020-03-31"

search = catalog.search(
    collections = [collection_id],
    intersects = geom,
    datetime = time_range,
    query = ["eo:cloud_cover<10"],
    limit = 100
)

print(f"{search.matched()} items found")

1~2行目では地理空間情報を扱うためにpandasを拡張した、geopandasを使ってGeoJSONを読み取り、最初の範囲（今回は範囲が1つしかないのでこれでOK）を抽出しています。

取り出した範囲を利用して、地理的な検索を行いました。

結果は以下のようになっており、8個のデータが取得できたことがわかります。

8 items found

検索結果はSTAC Itemと呼ばれるGeoJSON形式で整理されたデータが返ってきます。

このデータをパースし、geopandasのGeoDataFrameに変換していきます。

items_dict = search.get_all_items_as_dict()["features"]
items_df = pd.DataFrame(items_dict)
items_df["raw_geom"] = items_df["geometry"]

# geometryカラムをshaperlyのPolygonに形式に変換して、GeoDataFrameにする
items_df['geometry'] = items_df['geometry'].apply(lambda x: Polygon(x['coordinates'][0]))
items_gdf = gpd.GeoDataFrame(items_df, geometry='geometry', crs='EPSG:4326')

この中から不要なカラムを除去して、GeoJSONで出力しましょう。

# idとgeometryだけを抽出
items = items_gdf[["id", "geometry"]]
# gdfをgeojsonで出力
items.to_file('items.geojson', driver='GeoJSON')

出力されたitem.geojsonを見てみましょう。

VSCodeをお使いの方はGeo Data Viewerという拡張機能を利用することでVSCode上でも閲覧できますし、QGISと呼ばれる世界で最も利用されているオープンソースのGISソフトウェアを利用することでも閲覧可能です。

複数のポリゴンが同じ場所に重なっているのでとても見づらいですが、8個のポリゴンが表示されているのがわかります。

ポリゴン数は異なりますが、先ほどEO Browserで検索した時と同じ形状のポリゴンが表示されいますね。

衛星画像を取得する

検索したSTAC Itemには当然、実データ（GeoTIFFなど）のS3 URIなども一緒に格納されています。

いよいよこのデータを取得していきましょう。

以下のようにすると8個のアイテムのURLが取得できます。

たくさんのバンドがありますが、普通っぽい画像を使うのはやめて、今回はRed Edge（band5）を利用していきましょう！

# S3のURIが格納されている
# 全てのアイテムのURLを取得する
# 今回はband5のみを取得する
urls = []
for i in range(len(items_df)):
    urls.append(items_df['assets'][i]['B05']['href'])
urls

8個のURLが取得できました。

['https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/VN/2020/3/S2B_54TVN_20200331_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/WN/2020/3/S2A_54TWN_20200330_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/XN/2020/3/S2A_54TXN_20200330_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/VN/2020/3/S2A_54TVN_20200326_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/WN/2020/3/S2A_54TWN_20200326_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/WN/2020/3/S2B_54TWN_20200325_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/WN/2020/3/S2B_54TWN_20200321_0_L2A/B05.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/54/T/WN/2020/3/S2B_54TWN_20200315_0_L2A/B05.tif']

URLの形式になっていますが、直接ダウンロードすることはできず、データの取得にはAWSのアカウントが必要です。

今回は以下の設定は済んでいるとして話を進めていきますのでご了承ください。

• AWSのアカウントがある
• S3へのアクセス権限が設定されたユーザーのアクセスキーとシークレットがマシンに設定されている
  • もしくはセッショントークンが取得され有効期限内である、インスタンスにロールが設定されているなど

とにかく、boto3からS3へアクセスできるような状態になっていればOKです！

私はSSOを利用しているので、boto3のセッションを作成しました。

session = boto3.Session(profile_name='PowerUserAccess-xxxxxxxxxxxx')

URLの形式になっているので、S3 URIの形式に変換しつつ、いよいよデータをダウンロードしてみましょう。

# boto3のセッションを使って、S3のファイルをダウンロードする
# 最初の画像だけを取得
target_url = urls[0]

# URLを置換して、S3のURIに変換
s3_uri = target_url.replace('https://sentinel-cogs.s3.us-west-2.amazonaws.com', 's3://sentinel-cogs')

print(s3_uri)

bucket_name = s3_uri.split('/')[2]
key = '/'.join(s3_uri.split('/')[3:])
local_path = os.path.join(os.getcwd(), key.split('/')[-1])

session.resource('s3').Bucket(bucket_name).download_file(key, local_path)

gdalを利用してダウンロードしたデータを開き、matplotlibで可視化してみます！

# ファイルを開く
band5_image = gdal.Open(local_path)

# 画像をmatplotlibで表示する
plt.imshow(band5_image.ReadAsArray())

おぉー！それっぽい画像が表示されました！
ちゃんとマシン内に「B05.tif」というデータも取得できているかと思います。

これで好きな場所・好きな時間の地球の様子が確認でき、解析もできる環境が整いましたね！

前提知識が多すぎて大変だったかと思いますが、なんとかデータのダウンロードと可視化ができましたね…お疲れ様でした…

衛星データを編集してみる

…ということで、ここから先はオマケ程度ですが、この画像をちょろっとだけいじっていきたいと思います。

今回のデータは符号なし16bitの画像なので、値の範囲は0~65535になりますが、実際に格納されている値はどんなもんなんでしょうか？

ヒストグラムを出力してデータの範囲を見てみましょう。

# 横に長く表示する設定
plt.figure(figsize=(20, 3))
# ヒストグラムを表示する
plt.hist(band5_image.ReadAsArray().flatten(), bins=100, range=(0, 65535))

ほとんどの値が0付近に集中しているのがわかりますね。

なので、このデータを0~10000の範囲のみ抽出し、更にRed Edgeっぽく（？）表示させてみましょう！

# 表示する最小値と、最大値を指定して、画像を表示する
# Red Edgeなので、植生が多いほど、値が高くなる
# カラーの指定も行う
plt.imshow(band5_image.ReadAsArray(), vmin=0, vmax=10000, cmap='OrRd')
# カラーバーも表示
plt.colorbar()
# タイトルも表示
plt.title('Band5')

いい感じですね！笑

最後に、表示範囲だと海の領域が多いので、地表面のみ切り抜いてみましょう！

先ほどからちょこちょこ名前が上がっているGDALというツールは「地理情報に特化した画像データ処理プログラム」です。

これを利用して地理的に画像を切り抜くようなことも可能なんですが、それはまたの機会ということで、シンプルにピクセルの範囲を指定して切り抜きます！

北海道には羊蹄山という大きな山があるので、こちらを切り抜いていきましょう

# gdal_translateで切り出し位置と幅を指定して画像を切り出す
# 上の画像を見ながら、羊蹄山付近が切り抜けるようにピクセルを決定した
min_x = 4000
min_y = 2600
delta_x = 500
delta_y = 500

cliped_file_path = "B5_cliped.tif"

ds = gdal.Translate(
    cliped_file_path, band5_image, srcWin=[min_x, min_y, delta_x, delta_y]
)
ds = None

これを表示させてみましょう。

# 切り抜いた画像を表示する
cliped_image = gdal.Open(cliped_file_path)
plt.imshow(cliped_image.ReadAsArray(), vmin=0, vmax=10000, cmap='OrRd')
plt.colorbar()

だいぶカッコいい画像が作成できましたね！！

終わりに

みなさんいかがだったでしょうか！

必要な知識が多くちょっと大変だったと思いますが、これでも昔に比べればすっっっっっごく楽になったそうです！

実際、普段の業務でAWSを使っている方であればデータのダウンロード自体はEO BrowserとAWS CLIだけで出来てしまいます。

これを機に、みなさんもAWSにあるオープンデータを利用して自由に宇宙から地球を覗いてみてはいかがでしょうか！