はじめに

これまでに、SAR衛星や画像の原理や基礎知識を解説する記事を掲載してきました。

本記事では、Google Earth Engine（GEE）及び衛星データの標準規格であるSTAC APIを利用して、Sentinel-1 GRD画像を取得し、実際に可視化する手順を解説します。GRD（Ground Range Detected）画像とは、SAR画像を地図上の座標系（地上距離）に補正したもので、通常の地図と同じように扱えます。

衛星データやSARデータに触れたことがない方も、ぜひこの機会に試してみてください。

これからSAR画像（Sentinel-1 GRD画像）の解析を始めようという方の参考になれば幸いです。

なお、SAR画像を取得できるプラットフォームは、株式会社Tellusが運営するTellusをはじめ多く存在します。各プラットフォームは提供するデータや機能、料金体系などが異なるため、取得したい画像の種類や処理の目的、予算などに応じて最適なものを選ぶとよいでしょう。

今回ご紹介するGEEは、数多くの衛星データを取り揃えており、地理空間データを簡易に処理できる関数が豊富に用意されているのが大きな特徴です。特に、大量データを扱った時系列解析や、広範囲をカバーするデータセットを必要とするプロジェクトでは、GEEが有益です。ただし、商用利用の場合は有償となりますので、注意が必要です。

一方、STAC APIは無料で利用できるという点が大きな利点です。個人利用はもちろん、商用利用であってもコストを抑えながら衛星データを取り扱えるため、事業規模の大小問わず活用しやすいでしょう。STAC APIを使うと同じインターフェースで様々な衛星データカタログを検索取得できるため目的に合ったデータをスムーズに収集できます。

では、GEEの具体的な紹介に移ります。

Google Earth Engine（GEE）経由での取得からデータの可視化

GEEの概要

Google Earth Engine（GEE）とは、衛星画像を中心とした地理空間情報の解析・可視化ができるオンラインクラウドプラットフォームです。その名称から推測できるように Google が運営・管理しています。

従来、衛星画像の解析には、個別サイトから必要なデータを入手する必要があったり、データ容量が大きいためにハイスペックなパソコンを用意する必要があったりと様々な課題がありました。

そうした課題をまとめて解決したのが GEE となります。

GEE の登録

GEE を利用するには、まず公式サイトでアカウントを作成しましょう。

GEE のトップページ右上にある「Get Started」ボタン（下図）をクリックし、画面の指示に従って登録を進めてください。

登録時には、利用目的（非営利か商用か）やGoogle Cloudプロジェクトを所有しているかどうかなどを尋ねられます。途中でGoogle Cloudの無料トライアル登録画面が表示されますが、その画面の規約に同意すれば、無料トライアルの登録をスキップすることも可能です。GEEの利用にあたってGoogle Cloudの登録は不要なので、スキップしても構いません。

なお、非営利・研究・メディア・個人利用の場合は無料で使えますが、商用利用の場合は有償となる点にご注意ください。

詳しくはGEEウェブサイトにある「Noncomercial」ページ、「Commercial」ページ
を参照してください。

すべての登録情報を入力・確認すると、GEE を利用できるようになります。

Python での解析環境

GEEは当初、Webブラウザ上のみで使用できるサービスとして提供され、プログラミング言語もJavaScriptベースに限られていました。しかし、利用者の増加に伴い、Earth Engine Python APIがリリースされ、PythonでもGEEを扱えるようになりました。

さらに、それ以降、Qiusheng Wu氏によって、Python上で簡単にインタラクティブな解析や可視化を行うためのパッケージ「geemap」が開発され、利便性が一段と向上しています。

Geemap 公式サイト

Pythonを使ってGEEを動かす場合、特に指定された環境はありませんが、以下のようなインタラクティブ解析・可視化ができる環境が主な選択肢となります。

・Google Colaboratory (https://colab.google/)
・Jupyter Notebook（ローカル環境も含む）（https://jupyter.org/）
・QGIS with Python（https://qgis.org/ja/site/）

本記事では、Google Colaboratoryを使用します。

SAR 画像とGISデータの可視化

目標

東京都千代田区の行政区域データを使い、2019年1月3日に観測されたSentinel-1 GRD画像を取得して、千代田区の範囲を可視化することを目標にします。

画像取得条件は、以下の通りです。
-偏波：VV　※偏波とは何かは、前回の記事をご覧ください。
-観測モード：Interferometric Wide (IW)
-衛星進行方向：Descending（北から南への軌道）

GEEにおけるSentinel-1 GRD画像の概要

GEEには、Sentinel-1 GRD画像を含め多様な衛星データが蓄積されており、Earth Engine Data Catalogで確認できます。

カタログ情報によると、GEE 上のSentinel-1 GRD 画像（IDはCOPERNICUS/S1_GRD）には以下の前処理が施されているようです。
1.熱ノイズ除去
2.後方散乱係数への変換
3.オルソ補正（SRTM 30mDEMもしくはASTER DEMを使用）
4.デシベル単位（dB）への変換

今回はデシベル単位変換後の画像を使用しますが、デジタル単位変換前の画像（IDはCOPERNICUS/S1_GRD_FLOAT）も保管されています。

SARデータを解析する際は、手元のデータにどのような前処理が施されているかを把握しておく必要があります。前処理の有無や内容によって、解析結果が変わる可能性があるからです。解析を始める前に、必ず確認しておきましょう。

行政区域データの取得から可視化まで

行政区域データの取得から可視化までの流れをご紹介します。

今回は、東京都千代田区の行政区域データを使います。国土交通省のサイトから東京都の行政区域データをダウンロードしてください。

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html

ダウンロードした Zip ファイルを解凍すると、GeoJSON ファイルとシェープファイルの2種類が含まれています。今回は GeoJSON ファイルを使用します。

GEE で GeoJSON ファイルを扱う場合は、ee.〇〇（以下、「ee形式」）に変換する必要があります。以下のコードでは、geojson_to_ee 関数を使って変換を行っています。

PATH_GEOJSON = "N03-23_13_230101.geojson"


with open(PATH_GEOJSON, encoding='utf-8') as f:
    json_data = json.load(f)


ee_chiyoda = geojson_to_ee(json_data["features"][0])
print(type(ee_chiyoda))

このコードを実行すると、GeoJSON ファイルが ee.geometry.Geometry という形式に変換されます。次に、以下のコードで地図上に可視化してみましょう。

Map = geemap.Map()
Map.centerObject(ee_chiyoda, 13)
Map.addLayer(ee_chiyoda, {}, "東京千代田区")
Map

無事に可視化できました。

今回、GeoJSONファイルを利用したのでgeojson_to_ee関数を使ってee形式へと変換しました。シェープファイルを使う場合は、shp_to_ee関数でee形式に変換できます。

Sentinel-1 GRD画像の取得から可視化まで

次に、2019年1月3日に撮影された Sentinel-1 GRD 画像を取得し、千代田区を対象に可視化してみます。条件は以下の通りです。
-偏波：VV
-観測モード：Interferometric Wide (IW)
-衛星進行方向：Descending（地球の北側から南側に移動する軌道）

以下のコードを実行すると、2019年1月1日から12月31日までに該当条件を満たす GRD 画像（計 31 枚）が取得されます。

sentinel1 = ee.ImageCollection("COPERNICUS/S1_GRD")
# dB 変換前の場合は以下を使用: ee.ImageCollection("COPERNICUS/S1_GRD_FLOAT")


sentinel1_filtered = (
    sentinel1
    .filterDate("2019-01-01", "2019-12-31")
    .filterBounds(ee_chiyoda)
.filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VV"))
    .filter(ee.Filter.eq("instrumentMode", "IW"))
    .filter(ee.Filter.eq("orbitProperties_pass", "DESCENDING"))
)


print(sentinel1_filtered.size().getInfo()) # 31

実行結果として「31」が表示されます。これは、2019年1月から12月31日までのGRD画像が31枚あることを示しています。

また、少し複雑なコードとなりますが、下記のようにすることで任意の●枚目の画像の日付を取得できます。今回の例では1枚目の画像を取得しており、日付が2019年1月3日であることがわかります。

# 1枚目（.get() の引数が0）の画像を取得
s1_list = sentinel1_filtered.toList(sentinel1_filtered.size())
s1_first_image = ee.Image(s1_list.get(0))


# "system:time_start" プロパティから ee.Date オブジェクトを作成
date = ee.Date(s1_first_image.get('system:time_start'))


# 日付を "YYYY-MM-dd" 形式の文字列に変換
date_string = date.format('yyyy-MM-dd').getInfo()


# 結果を表示
print(date_string)

では、東京都千代田区内に限定したSentinel-1 GRD画像からVV バンドを抽出し、値のヒストグラムの下位2%と上位98%の値をパラメータとして設定して可視化します。

今回は下位2%と上位98%の値を求めるためにee.Reducer.percentile()関数を使いました。「Reducer」 という単語はGEE独特のものですが、対象領域のピクセルをまとめて単一の値（ここではパーセンタイル）に変換する機能だと考えると分かりやすいでしょう。

# Sentinel-1の1枚目の画像をクリップ
s1_first_clipped = s1_first_image.clip(ee_chiyoda)


# "VV" バンドの下位2%と上位98%の値を計算
percentiles = s1_first_clipped.select("VV").reduceRegion(
    reducer=ee.Reducer.percentile([2, 98]),
    geometry=ee_chiyoda,
    scale=10,
    maxPixels=1e9
)


# 計算結果から各パーセンタイル値を取得
min_value = percentiles.get("VV_p2")
max_value = percentiles.get("VV_p98")


# 可視化パラメータにパーセンタイル値を設定
vis_params = {"min": min_value, "max": max_value}


# レイヤーの追加
Map.addLayer(s1_first_clipped.select("VV"), vis_params, "Sentinel-1 First Image Clip")
Map.centerObject(ee_chiyoda, 13)
Map

可視化した画像をダウンロードしたい場合は、.ee_export_image_to_drive関数を使うことでGoolgeDriveにダウンロードすることが可能です。

# GoogleDrive上で画像をダウンロード
geemap.ee_export_image_to_drive(
    image=s1_first_clipped.select("VV"),
    description="sentinel1_first_clip",
    region=ee_chiyoda,
    scale=10)

追加情報になりますが、今回の画像には、VV と VH の両偏波情報が含まれています。一般的に、VV 偏波は建物など都市部で特に強い反射を示し、VH 偏波は森林や草地など自然環境で強い反射を示す特長があります。これら 2 種類の偏波情報を組み合わせることで、下図のような鮮やかな疑似カラー画像が作成できます。

なお、RGB では赤を VV、緑を VH、青を VV/VH に割り当てているため、各偏波の強度に応じてそれぞれの色が際立つ仕組みになっています。

疑似カラー画像は、単なるモノクロ画像とは違い、都市部と自然環境のコントラストを明確に示し、土地利用の状況をひと目で把握できる点が特長です。市街地の拡大や森林の分布などをわかりやすく可視化できるため、1枚の画像から各種情報を抽出したい時に役立ちます。

統計量（平均値）の算出

最後に、統計量（平均値）の算出方法をご紹介します。下記のコードで平均値を求めることができます。

# 統計量（平均値）の算出
mean = s1_first_clipped.select("VV").reduceRegion(
    reducer=ee.Reducer.mean(),
    geometry=ee_chiyoda)


# 平均値の表示
print(mean.getInfo()) # -2.345530128502706

今回は平均値を求めるために、ee.Reducer.mean() 関数を使いました。もし中央値や標準偏差といった他の統計量を計算したい場合は、ee.Reducer.median()やee.Reducer.std()などを使用できます。

STAC AP I経由での取得

STAC API の概要

STAC とは、衛星データ等の地理空間情報をまとめたカタログのことです。正式名称は、Spatio Temporal Asset Catalogs（STAC） です。

STAC では、取得日時や座標範囲、センサー情報などのメタデータを JSON や GeoJSON 形式にまとめることで、データ提供機関ごとに異なりがちな形式の差を解消し、検索やアクセス性を高めることを目的としています。詳しくは、こちらをご覧ください。

STAC API を利用する際は、どの API を使うかを選び、設定する必要があります。例えば、APIの一つにEarth Searchがあります。API の一覧はこちらで確認できます。

基本的に、Public APIs や Public Static Catalogs に分類される API は無償で利用できます。ただし、一部の API では事前登録が必要な場合があります。また、Protected APIs and Catalogs に分類される API の多くは、登録が必要または有償利用となります。

今回は、無償で多くのデータを利用できる Microsoft Planetary Computer を使用します。

Sentinel-1 GRD画像の取得から可視化まで

目標

先ほどと同様に、2019年1月3日に東京都千代田区で取得されたSentinel-1のGRD画像を使い、取得から可視化までを行うことを目標とします。条件は以下の通りです。
-偏波：VV
-観測モード：Interferometric Wide (IW)
-衛星進行方向：Descending

STAC API（Planetary Computer）におけるSentinel-1 GRD画像の概要

Planetary Computerの解説によると、Sentinel-1 GRD画像には以下の前処理が行われています。

　　1.マルチルック処理　※詳細は以前の記事をご覧ください。
　　2.測地系「WGS84」に基づくグラウンドレンジ（地上距離）への変換

一見すると、先ほどご紹介したGEE上のSentinel-1 GRD画像の前処理とは異なっているように思えますが、解説にESAのドキュメントへの言及があることから、基本的にはESAが標準的に行っている前処理が適用されていると推測できます。GEEの前処理も同じくESAに準拠しているため、最終的に Planetary Computer と GEE で取得される画像は同様の前処理を経ていると考えられます。とはいえ、Sentinel-1 の場合、前処理内容が明確に示されていないことも少なくないので、実際に利用・解析する際には、別のプラットフォームの画像と比較したり、目視で確認するなどしてチェックすることをおすすめします。

なお、放射強度の補正やPlanetDEMを用いたオルソ補正が施された画像が必要な場合は、Sentinel-1 Radiometrically Terrain Corrected（RTC）をご利用ください。

行政区域データの取得から可視化まで

先ほど国土交通省のウェブサイトからダウンロードした東京都の行政区域データ（GeoJSON ファイル）を使い、千代田区を抽出して可視化します。なお、後ほど活用するために、バウンディングボックス（bounding box）も作成しておきます。

# GeoJSONファイルのパスを指定
PATH_GEOJSON = "/content/drive/MyDrive/sorabatake_gee/N03-23_13_230101.geojson"


# GeoJSONファイルを読み込み、GeoDataFrameに変換
gdf = gpd.read_file(PATH_GEOJSON)


# "N03_004" 列が "千代田区" であるレコードをフィルタリング
chiyoda = gdf[gdf["N03_004"] == "千代田区"]


# 千代田区の領域のバウンディングボックス（最小限の外接矩形）を計算
chiyoda_bbox = chiyoda.total_bounds


# 千代田区の領域をプロット
chiyoda.plot();

Sentinel-1 GRD画像の取得から可視化まで

次に、東京都千代田区における 2019 年 1 月 3 日の Sentinel-1 GRD 画像を取得し、その可視化を行います。繰り返しになりますが、偏波は VV、観測モードは Interferometric Wide (IW)、衛星の進行方向（軌道）は Descendingです。

以下のコードを実行すると、2019年1月1日から12月31日までのGRD画像（31枚）を取得できます。

# # Microsoft Planetary Computerに接続
client = Client.open(
    url="https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace
)  


# 検索条件の設定
collection = "sentinel-1-grd"
datetime_range = "2019-01-01/2019-12-31"
bbox = chiyoda_bbox


# クエリ条件の設定
query_params = {
    "sar:instrument_mode": {"eq": "IW"},
    "sat:orbit_state": {"eq": "descending"}
}


# STAC APIを利用して、条件に合致するアイテムを検索
search = client.search(
    collections=collection,
    bbox=bbox,
    datetime=datetime_range,
    query=query_params
)


# 検索結果のアイテム群を取得
items = search.item_collection()


# 取得したアイテムの数を表示
print("取得したアイテム数:", len(items))

GEEの時と同様、31枚のGRD画像が取得できます。31枚のIDは下記コードで確認できます。

# 各アイテムのIDを順番に表示
for i, item in enumerate(items):
    print(f"{i}: {item.id}")

IDが逆順になっていますが、31枚目（IDは30）が 2019年1月3日の画像であることを確認できます。

早速、取得した画像を可視化したいところですが、Planetary Computerからのダウンロードには約10分ほどかかることがあります。

以下のコードでダウンロードできます。地理座標系としてEPSG:4326 を付与します。

# 30番目（2019年1月3日取得）のアイテムを選択し、
item = planetary_computer.sign(items[30])


# 選択したアイテムからVVバンドのアセットを取り出し、そのダウンロード用URLを取得
image_url = item.assets["vv"].href


# 書き出すファイル名を指定
OUTPUT_FILE_NEAREST = "/content/drive/MyDrive/sorabatake_gee_stac/image_georeferenced.tiff"


# EPSG:4326座標系に最近傍補間 (nearest)でリサンプリングを行って書き出し
gdal.Warp(OUTPUT_FILE_NEAREST, image_url, dstSRS="EPSG:4326", resampleAlg="nearest")

今回は、画像に位置情報を付加する際の補間手法として nearest neighbor という方式を使いましたが、bilinear など他の方式も選べるので、目的に応じて変更してください。

補間手法についてはこちらの記事でも解説しているので参考にしてください。

画像の書き出し（image_geocoded_nearest.tif）が無事に終われば（4分程度かかります）、フリーのGISソフト”QGIS”などで地図上に可視化できます。

今回のサンプルコードは、どなたでも気軽に再現できるよう、Planetary ComputerのAPI Keyを使わない形での紹介をしました。しかし、Planetary ComputerのAPI Keyを取得する形とすれば、特定の領域だけを効率よくダウンロードできるなど、柔軟なデータ取得や解析が可能になります。今回の解説を通じてSTAC APIの基本的な使い方を理解できた方は、ぜひAPI Keyを使った方法にもチャレンジしてみてください。

行政区域データを用いたマスク処理

先ほど書き出したデータは、東京都を含む Sentinel-1 のシーン全体をカバーしていますが、今回は千代田区のみを対象にしています。そこで、千代田区の行政区域データを用いて先ほど書き出した画像をマスク処理し、千代田区の範囲だけを切り出します。以下のコードでは、このマスク処理とマスク後の画像を書き出す手順を示しています。

# 書き出すファイル名を指定
OUTPUT_FILE_MASKED = "/content/drive/MyDrive/sorabatake_gee_stac/image_georeferenced_masked.tiff"


# 千代田区のジオメトリを取得
chiyoda_geom = [chiyoda.geometry.iloc[0]]


# 画像の読込とマスク処理
with rio.open(OUTPUT_FILE_NEAREST) as src:
  out_image, out_transform = rio.mask.mask(src, chiyoda_geom, crop=True)
  out_meta = src.meta.copy()


out_meta.update({
    "driver": "GTiff",
    "height": out_image.shape[1],
    "width": out_image.shape[2],
    "transform": out_transform
    })


# 画像をファイルに書き出し
with rasterio.open(OUTPUT_FILE_MASKED, "w", **out_meta) as dest:
    dest.write(out_image)

最後に、書き出したマスク済みファイルの可視化を行い、平均値などの統計情報を算出しましょう。

# 書き出した画像を読込
with rio.open(OUTPUT_FILE_MASKED) as src:
  img = src.read(1).astype("float32")


# 画像を可視化
p2, p98 = np.nanpercentile(img, [2, 98])
plt.imshow(img, cmap="gray", vmin=p2, vmax=p98);


# 画像の平均値を算出
img_mean = np.nanmean(img)
print(f"平均値: {img_mean}") # 385.673980712890

最後に

本記事では、GEE と STAC API を活用したSAR画像の取得から可視化までを取り扱いました。SAR画像を解析する際にはご参考いただけると幸いです。

繰り返しにはなりますが、SAR画像を解析する場合、自分が扱うSAR画像において何の前処理が施されているのかを確認・理解しておくことが必要です。

今回、GEE と STAC API（Planetary Computer）が提供するSentinel-1 GRD画像においてほぼ同じ前処理が施されていることを紹介しました。

この記事以降は、SAR画像を活用した解析記事をいくつか掲載していく予定ですので、引き続きフォローしていただけると嬉しいです。

本記事のコード（Google Colab）

本記事で使用したコードとGISデータは、以下のリンクからご覧いただけます。
本記事で使用したコードとGISデータ（Google Drive）