「ESAのセンチネルAPIに深く潜入」

「ESAのセンチネルAPIに潜入！美容・ファッション専門家が解説！」

Pythonを使用して衛星画像を取得、分析、視覚化する方法

この記事のすべての画像は著者によって作成されました。

欧州宇宙機関は、Copernicusを支援するためのSentinelミッションを実施しており、陸地、海洋、大気のモニタリングに使用されるレーダーと多光谱イメージング機器など、さまざまなリモートセンシングデータを提供しています。

現在、実行中のSentinelミッションは6つあり、そのうち3つはPython APIを介して簡単にアクセスできます。それらは、公式ソースによると次のとおりです：

Sentinel-1は、陸地と海洋のサービスのための全天候型の偏極軌道レーダーイメージングミッションです。Sentinel-1Aは2014年4月3日に、Sentinel-1Bは2016年4月25日に打ち上げられました。両者ともフランス領ギアナのヨーロッパ宇宙基地からソユーズロケットで軌道に乗せられました。Sentinel-1Bのミッションは2022年に終了し、Sentinel-1Cの導入計画が進行中です。

• 合成時系列データ生成としてのLLM
• 「Pythonを使用した最も近いバーを見つけるための近接解析」
• 戦略的データ分析（パート3）：診断的な問い

Sentinel-2は、陸地のモニタリングのための偏極軌道の多光谱高解像度イメージングミッションであり、植生、土壌および水域の被覆、内陸水路および沿岸地域のイメージなどを提供します。Sentinel-2は緊急サービスにも情報を提供できます。Sentinel-2Aは2015年6月23日に、Sentinel-2Bは2017年3月7日に打ち上げられました。

Sentinel-3は、高い精度と信頼性で海洋や陸地の表面の測定を行うマルチインストゥルメントミッションです。このミッションは、海洋の予測システムおよび環境・気候のモニタリングをサポートしています。Sentinel-3Aは2016年2月16日に、Sentinel-3Bは2018年4月25日に打ち上げられました。

さらに調査すると、Sentinel-1のデータは数メートルほどの空間分解能にまで達します。そして、視覚データの最高解像度はSentinel-2の10メートルであり、Sentinel-3はセンサータイプによっては100kmの規模で運用されます。

では、衛星データを取得する場所がわかりましたし、ソース（センサー）と空間解像度の豊富な選択肢もあります。このようなさまざまなタイプの衛星データはどのように使用されるのでしょうか？まずはじめに、以下は50以上のユースケースの一部です。

一般的な経験則では、ユースケースや問題の具体的な要素、ターゲット地域の地理特性や地形など、すべてがデータソースを選ぶ際に重要な要素です。ただし、日常の実践では、以下が主な要因となります：

• 価格（無料で利用可能なSentinelが好ましい）
• 数メートル程度の空間解像度で、さらに小さな都市構造物も捉えることができる
• 可視光線および近赤外線など、少なくともいくつかのバンドを持つ
• 時間的な頻度

これらの特徴により、Sentinel-2は地理空間データコミュニティで最も広く使用される衛星データソースであると言えます。この記事では、これらの要素を基に、Sentinelデータの取得方法とダウンロード時の予想される内容について詳しく説明します。また、画像レコードとそれに内包される情報の時間的な変化と異なる可能性についても詳しく探究します。

この記事では、EUの法律によってCopernicus Sentinel Data and Service Informationへの無料アクセスが認められているため、Copernicus Sentinelデータ2023を使用しました。

1. データ取得

まず、公式ドキュメントとサンプルコードに従ってAPI接続を設定します。さらに、画像をダウンロードする対象エリアが必要です。デバッグしやすいように、自分の故郷であるブダペストを選びました。そのために、OSMNxを使用して行政区画をダウンロードします。

import osmnx as ox # version: 1.0.1import matplotlib.pyplot as plt # version: 3.7.1city = 'Budapest'admin = ox.geocode_to_gdf(city)admin.plot()

次にSentinel APIにアクセスします：

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt # version 0.14# アカウントを取得するには、こちらでサインアップしてください：https://apihub.copernicus.eu/apihubuser = <ユーザー名を追加 password = <パスワードを追加 >api = SentinelAPI(user, password, 'https://apihub.copernicus.eu/apihub')  

クエリを実行するためには、場所を指定するスムーズなポリゴンを用意する方が良いです。そのために、ブダペストの行政区画の凸包を作成しました：

# クエリを簡素化するために、入力ポリゴンの凸包を抽出しますadmin_polygon = admin.convex_hull.geometry.to_list()[0]admin_polygon

出力：

指定した場所で特定の期間内、選択したプラットフォームの衛星画像を検索することができます。また、クラウドカバレッジに基づいてフィルタリングもできます。つまり、画像があまりにも曇っている場合はすぐに破棄します。

# ここで、場所（ポリゴンに基づいて）# 時間枠# スペースプローブ# 雲カバー率のレベルを指定できますproducts = api.query(admin_polygon,                     date=('20150623', '20231006'),                     platformname='Sentinel-2',                     cloudcoverpercentage=(0, 100))len(products)

これらのセルの出力を見ると、Sentinelのドキュメントに従って、2015年6月23日（ミッションの開始）から2023年10月6日（この記事を書いていた時点）までに、合計3876枚の衛星画像がブダペストの行政区画と重なって記録されました。クラウドカバー率を100に設定したため、クラウドカバレッジに基づいてフィルタリングは適用されませんでした。したがって、過去8年間のすべての画像識別子が揃っているはずです。

また、製品リストには衛星画像の識別子とメタデータが含まれていますが、画像そのものは含まれていません。また、同じ操作をSentinel-3で行うと、解像度ははるかに低いですが、結果として約2万件の画像レコードを取得することができます。

2. メタデータの探索

製品リストをPandas DataFrameに変換して分析を開始しましょう！

import pandas as pd # version: 1.4.2products_gdf = api.to_geodataframe(products)products_gdf = products_gdf.sort_values(['beginposition'], ascending=[True])print(products_gdf.keys())print(len(products_gdf.keys()))products_gdf.head(3)

このブロックの結果：

衛星画像識別子によってインデックス付けされたテーブルのキーの数をカウントすると、それだけで41 の特徴量のカラムを含むデータの豊かさがわかります。

これらのフィールドには多くの技術情報が含まれていますが、いくつかの特徴に注目してみましょう。一方、空間と時間の次元は、生成日時と開始位置（日時情報として）、ジオメトリ（ポリゴン、GIS、データ型として）でエンコードされています。他方、画像に基づいて土地被覆タイプを特徴付ける興味深いメトリクスがいくつかあります。すでにクエリで目にしたcloudcoverpercentage、vegetationpercentage、waterpercentage、snowicepercentage です。これらの環境指標は、異なる材料のスペクトル特性に基づいて派生されています。注意：これらの値は、タイル全体の総平均を捉えた集約スコアです。このトピックについて詳しくはこちらをご覧ください。

3. 空間の次元

ジオメトリの次元があるので、地図上の見え方を見てみましょう！以下のコードでは、ランダムなタイルのセットを可視化しています。これは何回かの実行の後、完全に代表的なものです。可視化には、FoliumとCartoDB Dark_Matterのベースマップを使用しました。

import foliumimport geopandas as gpdx, y = admin_polygon.centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=8, tiles='CartoDB Dark_Matter')# 一部のランダムなタイルを可視化するための設定polygon_style = { 'fillColor': '#39FF14', 'color': 'black',  'weight': 3, 'opacity': 0}geojson_data = products_gdf[['geometry']].sample(10).to_json()folium.GeoJson(    geojson_data,    style_function=lambda feature: polygon_style).add_to(m)# 管理境界を上に追加admin_style = {'fillColor': '#00FFFF',  'color': 'black','weight': 3, 'opacity': 100.0  }admin_geojson_data = admin[['geometry']].to_json()folium.GeoJson(    admin_geojson_data,    style_function=lambda feature: admin_style).add_to(m)# 地図を表示m

このコードブロックの出力：

このビジュアルを見ると、いくつかのタイルが繰り返し表示されることが明らかになります。また、何らかの理由で、これらのタイルのいくつかは都市の行政境界を2つに分割しています。これは、データを完全にターゲットエリアと重なり合わせて分析したい場合に不都合な状況になる可能性があります。解決策の1つは、所望の行政エリアと完全にオーバーラップしないタイルをフィルタリングすることです：

def compute_overlapping_area(tile, admin):    return tile.intersection(admin_polygon).area / admin_polygon.areaproducts_gdf['overlapping_area_fraction'] = products_gdf.geometry.apply(lambda x: compute_overlapping_area(x, admin_polygon))products_gdf_f = products_gdf[products_gdf.overlapping_area_fraction==1]print(len(products_gdf))print(len(products_gdf_f))products_gdf_f.head(3)

このセルの結果：

このフィルタを適用することで、タイルの約半分が削除されました。さて、マップを見てみましょう。すべてのタイルが市の境界と完全に重なり、都市が半分に分割されていないか確認してみましょう：

import foliumimport geopandas as gpdx, y = admin_polygon.centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=8, tiles='CartoDB Dark_Matter')# 一部のランダムなタイルを可視化するための設定polygon_style = { 'fillColor': '#39FF14', 'color': 'black',  'weight': 3, 'opacity': 0}geojson_data = products_gdf_f[['geometry']].sample(10).to_json()folium.GeoJson(    geojson_data,    style_function=lambda feature: polygon_style).add_to(m)# 管理境界を上に追加admin_style = {'fillColor': '#00FFFF',  'color': 'black','weight': 3, 'opacity': 100.0  }admin_geojson_data = admin[['geometry']].to_json()folium.GeoJson(    admin_geojson_data,    style_function=lambda feature: admin_style).add_to(m)# 地図を表示m

4. データセットの時間的次元

まずは、ブダペストをカバーする1日、1週間、1ヶ月あたりの画像枚数を確認しましょう。時間を測るために、beginpositionフィールドに頼ります。

# 前提として、'beginposition' はあなたのGeoDataFrame内のタイムスタンプの列だと仮定します
# 'beginposition' をDateTimeインデックスに変換します
products_gdf_f_cntr = products_gdf_f.copy()
products_gdf_f_cntr['beginposition'] = pd.to_datetime(products_gdf_f_cntr['beginposition'])
products_gdf_f_cntr.set_index('beginposition', inplace=True)

# データを日単位、週単位、月単位でリサンプリングして、行数をカウントします
daily_counts = products_gdf_f_cntr.resample('D').count()
weekly_counts = products_gdf_f_cntr.resample('W').count()
monthly_counts = products_gdf_f_cntr.resample('M').count()

# サブプロットを作成して、データを表示します
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
for idx, (count_name, count_val) in enumerate([('毎日の枚数', daily_counts), ('毎週の枚数', weekly_counts), ('毎月の枚数', monthly_counts), ]):
    ax[idx].plot(count_val.index[0:250], count_val['geometry'].to_list()[0:250])
    ax[idx].set_xlabel('日付')
    ax[idx].set_ylabel('枚数')
    ax[idx].set_title(count_name)

plt.tight_layout()
plt.suptitle('異なる時間枠で撮影された衛星画像の枚数', fontsize=20, y=1.15)
plt.show()

これらの図は、Sentinel-2プローブの最初の250日、最初の250週、および最初の250ヶ月（全期間）を示しています。最初の図は、毎日のスナップショットが隔日で撮られていることを示しています。2番目のグラフは、前のグラフの週間平均を取っており、最初の2年間は週に1〜2回、2017年から2018年にかけては週に5〜6回ブダペストが撮影されていることを示しています。最後のグラフは、全期間を表示しており、これらの25枚の画像が月平均の標準レベルになったことを示しています。

5. 陸地被覆変数の時間次元

次に、vegetationpercentage、waterpercentage、snowicepercentage、およびcloudcoverpercentageの時間的な変化を見てみましょう。先ほどの画像が示すように、初期の年には異なる、おそらくノイズの結果が表示される場合があるため、注意が必要です。ここでは、データのこれらの年を削除しないでおくことにします。合計して8年間のデータがあり、そのうち3年をカットするのは情報の損失が大きすぎる可能性があるためです。まずは、週単位の集計で時間の経過とともに生の値を見てみましょう。

import pandas as pd
import matplotlib.pyplot as plt

# 前提として、'beginposition' はあなたのGeoDataFrame内のタイムスタンプの列だと仮定します
# 'beginposition' をDateTimeインデックスに変換します
products_gdf_f_cntr = products_gdf_f.copy()
products_gdf_f_cntr['beginposition'] = pd.to_datetime(products_gdf_f_cntr['beginposition'])
products_gdf_f_cntr.set_index('beginposition', inplace=True)

# データを週単位で平均値を計算してリサンプリングします
weekly_averages = products_gdf_f_cntr[['vegetationpercentage', 'waterpercentage', 'snowicepercentage', 'cloudcoverpercentage']].resample('W').mean()

# 4つのサブプロットを持つ複数プロットの図を作成します
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 15))

# 'vegetationpercentage' を緑の線でプロットします
ax1.plot(weekly_averages.index, weekly_averages['vegetationpercentage'], color='green', label='週平均植生割合')
ax1.set_xlabel('日付')
ax1.set_ylabel('割合')
ax1.set_title('週平均植生割合')
ax1.legend()

# 'waterpercentage' を青の線でプロットします
ax2.plot(weekly_averages.index, weekly_averages['waterpercentage'], color='blue', label='週平均水割合')
ax2.set_xlabel('日付')
ax2.set_ylabel('割合')
ax2.set_title('週平均水割合')
ax2.legend()

# 'snowicepercentage' をシアンの線でプロットします
ax3.plot(weekly_averages.index, weekly_averages['snowicepercentage'], color='cyan', label='週平均雪/氷割合')
ax3.set_xlabel('日付')
ax3.set_ylabel('割合')
ax3.set_title('週平均雪/氷割合')
ax3.legend()

# 'cloudcoverpercentage' を灰色の線でプロットします
ax4.plot(weekly_averages.index, weekly_averages['cloudcoverpercentage'], color='gray', label='週平均雲カバー割合')
ax4.set_xlabel('日付')
ax4.set_ylabel('割合')
ax4.set_title('週平均雲カバー割合')
ax4.legend()

plt.tight_layout()
plt.show()

同様に、月間集計と共に:

これらの時系列データはいくつかの興味深い点を示しています:

• 植生の割合は明らかに季節変化を示しており、毎年の春にはすべてが緑に変わり、秋にはこの緑が薄れて50〜60％からほぼゼロに減少します。
• それに比べて、水の割合は年間通して0.8％前後で変動しており、観測期間全体にわたって一貫しています。これは研究エリアの表面水の量が非常に少ないためです。それでも、水の減少は冬季により頻繁に起こるようであり、一部の淡水部分が凍結することを示しています。
• 雪に関しては、最も顕著なピークは冬に発生します。それでも、これらの冬の大部分にいた私の経験から言えることは、私たちにはあまり雪がなかったということです。したがって、特に冬季に1〜2％の値が測定されることは、雲のノイズや誤分類の原因となる可能性があります。
• 雲については、植生に大きく関連しており、季節パターンに従っています。

これらの観察結果は、これらの指標の相関でも見ることができます:

products_gdf_f_cntr[['vegetationpercentage', 'waterpercentage', 'snowicepercentage', 'cloudcoverpercentage']].corr()

6. 衛星画像のダウンロード

まず、Sentinel-2とSentinel 3の両方を同じ週の8月にクエリし、できるだけ雲カバーを制限します。利用可能なスナップショットの数を確認してください:

# タイル製品IDのクエリproducts_sent = api.query(admin_polygon, date=('20230806', '20230813'), platformname='Sentinel-2', cloudcoverpercentage=(0, 1))products_sent = api.to_geodataframe(products_sent)f, ax = plt.subplots(1,1,figsize=(6,4))admin.plot(ax=ax, color = 'none', edgecolor = 'k')ax.set_title('Sentinel-2, タイルの数 = ' + str(len(products_sent)))products_sent.plot(ax=ax, alpha = 0.3)# Budapestと完全に重ならないタイルをフィルタリングするproducts_sent['overlapping_area_fraction'] = products_sent.geometry.apply(lambda x: compute_overlapping_area(x, admin_polygon))products_sent = products_sent[products_sent.overlapping_area_fraction==1]f, ax = plt.subplots(1,1,figsize=(6,4))admin.plot(ax=ax, color = 'none', edgecolor = 'k')ax.set_title('Sentinel-2, タイルの数 = ' + str(len(products_sent)))products_sent.plot(ax=ax, alpha = 0.3)len(products_sent)

このコードブロックの結果:

次に、イメージデータをダウンロードします:

# 最初のタイルを衛星画像としてダウンロードproduct_ids = products_sent.index.to_list()for prod in product_ids:    api.download(prod)

注記— ここではこの通知を受け取る場合もあります。それの場合は数時間待ってからダウンローダーを再実行してください。

製品a3c61497-d77d-48da-9a4d-394986d2fe1dはオンラインではありません。長期アーカイブからの回収をトリガーします。

7. ダウンロードした画像を開いて可視化する

ここでは、フォルダ構造に関する非常に良いビジュアルとデータ形式の詳細な説明を見つけることができます。画像ディレクトリを開くと、異なるバンドが見つかるでしょう。13つのバンドのそれぞれの意味と空間分解能は、 この 記事によってうまくまとめられています。13つのバンドの空間分解能は10〜60メートルです。いくつかのハイライト:

• 青色(B2)、緑色(B3)、赤色(B4)、および近赤外線 (NIR) (B8) チャンネルは10メートルの解像度です。
• 次に、植生の赤外線(Red Edge) (B5)、近赤外線 (NIR) (B6、B7、およびB8A)、および短波長赤外線 (SWIR) (B11およびB12) チャンネルは10メートルの解像度です。
• 最後に、沿岸エアロゾル (B1) と SWIRシラスバンド (B10) は60メートルのピクセルサイズです。

こうなります：

# ダウンロードしたフォルダを解凍した後:import osimage_path = 'S2B_MSIL1C_20230810T094549_N0509_R079_T34TCT_20230810T124346.SAFE/GRANULE/L1C_T34TCT_A033567_20230810T095651/IMG_DATA'sorted(os.listdir(image_path))

これが、rasterio を使用して、B4、赤色バンドを可視化した場合のタイル全体の見た目です：

import rasteriofrom rasterio.plot import showimage_file = 'T34TCT_20230810T094549_B04.jp2' with rasterio.open(image_path + '/' + image_file) as src:        image = src.read(1)  # 必要に応じてバンドインデックスを変更    plt.figure(figsize=(10, 10))    plt.imshow(image, cmap='Reds')  # カラーマップを変更できます    plt.title(image_file)    plt.colorbar()    plt.show()

出力:

次に、ブダペストに焦点を当て、都市の行政境界によってフルラスタータイルをマスクし、R、G、およびBバンドを個別に処理します：

from rasterio import maskf, ax = plt.subplots(1,3,figsize=(15,5))for idx, (band_name, band_num, color_map) in enumerate([('ブルー', 'B02', 'Blues'), ('グリーン', 'B03', 'Greens'), ('レッド', 'B04', 'Reds')]):       raster_path = image_path + '/T34TCT_20230810T094549_' + band_num + '.jp2'    with rasterio.open(raster_path) as src:        polygons = admin.copy().to_crs(src.crs)        geom = polygons.geometry.iloc[0]        masked_image, _ = mask.mask(src, [geom], crop=True)    ax[idx].imshow(masked_image[0], cmap=color_map)    ax[idx].set_title('ブダペスト センチネル 2 - ' + band_name + ' バンド')

最後に、これらを組み合わせてブダペストのRGBイメージを作成しましょう。まず、フルタイルをRGBイメージに組み立て、それを読み込んで公式指示に従ってヒストグラム均等化を行い、最終的なプロットに至ります。

# バンドの場所を取得band_blue = '/T34TCT_20230810T094549_B02.jp2'band_green = '/T34TCT_20230810T094549_B03.jp2'band_red = '/T34TCT_20230810T094549_B04.jp2'# バンドを読み込んで完全なRGBタイルを作成b2   = rasterio.open(image_path + '/' + band_blue)b3   = rasterio.open(image_path + '/' + band_green)b4   = rasterio.open(image_path + '/' + band_red)# フルタイルをtifファイルとしてエクスポートするmeta = b4.metameta.update({"count": 3})prefire_rgb_path = 'budapest_rgb.tif'with rasterio.open(prefire_rgb_path, 'w', **meta) as dest:    dest.write(b2.read(1),1)    dest.write(b3.read(1),2)    dest.write(b4.read(1),3)        # 切り抜かれたバージョンを読み込んで表示from skimage import exposureimg = rasterio.open('budapest_rgb_cropped.tif')image = np.array([img.read(3), img.read(2), img.read(1)])image = image.transpose(1,2,0)# ヒストグラム均等化を行うp2, p98 = np.percentile(image, (2,98))image = exposure.rescale_intensity(image, in_range=(p2, p98)) / 100000f, ax = plt.subplots(1,1,figsize=(15,15))rasterio.plot.show(image.transpose(2,0,1), transform=img.transform, ax = ax)ax.axis('off')plt.savefig('budapest_rgb_cropped.png', dpi = 700, bbox_inches = 'tight')

出力:

まとめ

まとめると、この記事で起こったことは以下の通りです:

• センチネル衛星プラットフォームの概要
• 画像識別子のクエリ方法とそのメタデータの取得例
• タイルの集約情報を基にした時間的な分析の簡単な方法
• 単一の画像のダウンロード、保存、可視化方法

これらのステップは、都市計画から環境モニタリングや農業まで、さまざまなユースケースにおいて衛星画像処理と分析を日常的な地理空間データサイエンスツールに追加することを目指しています。

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!