Pythonで衛星画像処理 ~MajorTOM~

この記事では、ESAのΦ Labが提供する衛星画像データセット用フレームワーク「Major TOM」の紹介と、pythonを使った簡単な衛星画像処理をしたいと思います。(実行環境はGoogle Colabratoryです。)

衛星画像処理をしようと思うと、データを入手するのに手間がかかったり、専用のソフトウェアの使い方を覚えなければならなかったりと、何かとハードルが高いものです。

しかし、このMajor TOMを使うと簡単に衛星画像データにアクセスできます。
また、この研究では、MajorTOM-Coreという地球の大部分をカバーするオープンアクセスなデータセットを提供しています。これを用いれば、アカウント登録などをせずに衛星データを取得できます。

Major TOMのリンク
Git Hub
Hugging Face
arXiv

対象読者

• pythonで簡単に衛星画像処理をしてみたい方
• 複数の衛星データを結合して利用するためのフレームワークを知りたい方

Major TOMとは

Major TOM (Terrestrial Observation Metaset)は、ESA(Europian Space Agency)のΦ Labが提供する大規模な衛星画像データセット用のフレームワークです。

従来、衛星データはそれぞれ特有のデータ構造をしており、複数のデータセットを組み合わせて使うことが難しいという問題がありました。そこで、複数のデータセットを組み合わせてアクセスできる共通の枠組みとしてこのMajor TOMが提案されました。将来的に、この枠組みに基づいて、様々なデータセットを組み合わせて使うことが容易になっていくと考えられます。

衛星画像処理

早速ですが、pythonで衛星からのデータの処理をしていきたいと思います。

使用するデータセット　「MajorTOM-Core」

MajorTOM-CoreはESAによって提供されている大規模でオープンな衛星データセットです。
地表の大部分(各国の軍事施設などを除く)をカバーしており、誰でも簡単にアクセスして使えるデータセットです。

現時点(2024年3月)でSentinel 2が観測したデータを公開しており、「Level 2A」、「 Level 1C」の２種類のデータが提供されています。
(この辺りは、自分は詳しくないのですが、このLevelというのは事前処理のレベルで、Level 1Cにシーン分類や気象に関する補正を加えたものがLevel 2Aということだと思われます:参考)

衛星データのバンド等の詳細は、「Sentinel Hub」から確認できます。

データの読み込み

公開している衛星データをすべて読み込もうとすると、そのデータ量が膨大であるため大変です。
そこであらかじめ欲しい領域のデータを経度、緯度の範囲で指定してデータを取得します。

その際便利なのが、MajorTOM-Core-ViewerというHugging FaceのSpaceです。
地図上でデータを取得したい地点まで移動して、「FIND SAMPLE」のボタンを押すとその地点を含む衛星画像の中心座標(経度、緯度)とその地点のデータがどのフォルダに含まれているか表示されます。(時々うまく行かないことがあります。その際はページをリフレッシュするとうまくいきます。)

コード

今回は試しにLondonの中心部のデータを処理していきたいと思います。

データの読み込み

データの読み込みまではこちらの公式のColabを参考にしてください。

「🔥 PyTorch Dataset with a local copy」のセルまで実行してください。

ちなみに私はLondonを含む領域を

「1. 📅 Filtering based on location, time, and cloud cover」のセクション２つ目のセルで

london = box(-0.1,51.0,0.0,51.9)

(東経-0.1~東経0、北緯51~北緯51.9の四角形)として実行しています

その直下のセルで region = london　としてください。

filtered_df = filter_metadata(gdf,
                              cloud_cover = (0,10), 
                              region=london
                              daterange=('2020-01-01', '2025-01-01'), 
                              nodata=(0.0,0.0) 
                              )

指標の計算と可視化

次に衛星データが格納されているフォルダを指定すると、いくつかの指標(NDVI,NDWIなど)を計算して可視化する関数を作成します。
指標の計算については、

• 宙畑のサイト、Sentinel Hub
• Sentinel HubのページNDVI、NDWI

を参考にしました。

詳しい解説は省きますが、この関数で求める指標は
NDVI (Normalized difference vegetation index) :植生を知るために使う指標
NDWI (Normalized Difference Water Index) :　水圏を知るために使う指標
です。
また、Sentinel-2の「band01」は20nmの波長の強度で、エアロゾルの検知(水蒸気等の影響もあるため補正が必要)に用いられるため可視化しています。

import os
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from matplotlib import gridspec

def Basic_Indices(path):
  Pathlist = os.listdir(path)
  Pathlist = sorted(Pathlist)

  band1 = plt.imread(path +'/'+ Pathlist[0]).astype(np.float32)
  band2 = plt.imread(path +'/'+ Pathlist[1]).astype(np.float32)
  band3 = plt.imread(path +'/'+ Pathlist[2]).astype(np.float32)
  band4 = plt.imread(path +'/'+ Pathlist[3]).astype(np.float32)
  band5 = plt.imread(path +'/'+ Pathlist[4]).astype(np.float32)
  band6 = plt.imread(path +'/'+ Pathlist[5]).astype(np.float32)
  band7 = plt.imread(path +'/'+ Pathlist[6]).astype(np.float32)
  band8 = plt.imread(path +'/'+ Pathlist[7]).astype(np.float32)
  band9 = plt.imread(path +'/'+ Pathlist[8]).astype(np.float32)
  band10 = plt.imread(path +'/'+ Pathlist[9]).astype(np.float32)
  band11 = plt.imread(path +'/'+ Pathlist[10]).astype(np.float32)
  RGB = Image.open(path +'/'+ Pathlist[13])
  RGB = np.array(RGB)

  NDVI = (band8-band4)/(band8+band4)
  NDWI = (band3 - band8)/(band3+band8)

  

  fig = plt.figure(figsize=(20,5))
  spec = gridspec.GridSpec(ncols=4, nrows=1, width_ratios=[1, 1, 1, 1],height_ratios=[1])
  ax0 = fig.add_subplot(spec[0])
  
  ax1 = fig.add_subplot(1,4,1)
  ax1.set_aspect('equal')
  ax2 = fig.add_subplot(1,4,2)
  ax2.set_aspect('equal')
  ax3 = fig.add_subplot(1,4,3)
  ax3.set_aspect('equal')
  ax4 = fig.add_subplot(1,4,4)
  ax4.set_aspect('equal')

  ax1.imshow(RGB)
  NDVImap = sns.heatmap(NDVI,vmin =-1.0 ,vmax = 1.0,cbar = False,cmap='YlGn',ax=ax2)
  NDWImap = sns.heatmap(NDWI,vmax = 1,cbar = False,cmap='GnBu',ax=ax3)

  band1map = sns.heatmap(band1,cbar = False,cmap = 'CMRmap',ax=ax4)
  
  
  band1map.set_title("443nm")
  NDVImap.set_title('NDVI')
  NDWImap.set_title('NDWI')

  ax0.axis("off")
  ax1.axis("off")
  ax2.axis("off")
  ax3.axis("off")
  ax4.axis("off")
  plt.show()

  return

Londonのデータが入っているフォルダを指定して実行(直下に各バンドのデータが入っているフォルダを指定してください。)

Basic_Indices('/content/data/L2A/573U/573U_1L/S2B_MSIL2A_20171123T111349_N0206_R137_T30UXC_20171123T132843')

実行結果

Londonの実行結果
※443nmの図の上のタイトル20nmとなっていますがそれはバンドの幅であり波長の中心は443nmです

NDVIを可視化した画像では、公園などで植物が生えているところが濃い緑で表示されているのがわかります。
また、NDWIを可視化した画像では、テムズ川が青く表示されているのがわかります。

最後に

記事を最後まで読んでいただきありがとうございました。間違っている点、不明な点があれば、遠慮なくコメント欄にてご指摘をお願いします。

参考文献

宙畑

Sentinel Hub

Major TOMのリンク(再掲)
Git Hub
Hugging Face
arXiv