Sentinel-1(SAR)とrasterioで洪水による被害範囲がどのように見えるのかを見てみた
本記事では、SARの観測データを用いると洪水がどのように見えるのか、衛星データ解析の初学者に向けてPythonによる分析コードを含めてご紹介します。
1. SARの復習
合成開口レーダー（SAR）は雲を突き抜けるマイクロ波によって地表の凹凸を観測します。この雲に関係なく観測できるという性質は、光学衛星にはない優位性です。
SAR（合成開口レーダ）のキホン～事例、分かること、センサ、衛星、波長～
天候に関わらず地表の様子が分かるということは、荒天に伴う災害時に力を発揮します。利用事例としては、例えばフィンランドの宇宙ベンチャー ICEYE がSARの観測データを洪水保険の開発に用いています。
SAR衛星事業者ICEYE、ガーナの洪水保険の開発プロジェクトに参画。被害状況の把握を担う【宇宙ビジネスニュース】
本記事では、SARの観測データを用いると洪水がどのように見えるのか、衛星データ解析の初学者に向けてPythonによる分析コードを含めてご紹介します。本記事で紹介している分析コードはGighubにも掲載していますので、ぜひ活用いただけますと幸いです。
2. SARと洪水の水面
なぜSARで洪水による被害範囲がわかるのでしょうか？少し考えてみましょう。
SARは衛星が発した電波の跳ね返り（後方散乱）を観測し、その跳ね返り具合によって地表の状態を判別できます。建物など凹凸が多い構造物や金属がむき出しの構造物では強い電波が返ってきます。建物ほどではなくても、植物なども程々の電波を返します。一方で水面のように滑らかで平らな面では衛星へと跳ね返る電波は非常に弱くなります。
洪水が発生すると、反射面が植物や建物などではなく、水面になります。すると返ってくる電波は、平常時と比べて非常に弱くなります。
つまり以前の観測と比べてSAR衛星が観測した電波の強度が下がっていれば、そこに新たな水面が出来ている、洪水が起きている可能性があると判断できます。もちろん、田んぼに水が張られただけかもしれませんので、最終的には他の情報も加味して考えます。
3. SARと洪水の文献紹介
近年の日本での大規模な洪水の1つに平成30年の西日本豪雨災害があります。その時に洪水に見舞われた岡山県倉敷市真備町の状況をJAXAの衛星ALOS-2のPALSAR-2で観測した論文があるので紹介します。
H. Gokon, F. Endo, and S. Koshimura (2023) Detecting Urban Floods with Small and Large Scale Analysis of ALOS-2/PALSAR-2 Data, Remote Sensing, 15(2), 532; https://doi.org/10.3390/rs15020532
論文では、洪水時の2018年7月7日と洪水前の2016年7月9日について、土地の勾配の影響を加味した後方散乱係数 σ0 (sigma nought、sigma naught とも。英語と米語）の変化を比べています。
田畑と思われるやや暗く映っている区域の中に、2016年は白く映っていたのに2018年の災害時は黒く消えている箇所があります。比較的わかりやすい洪水の痕跡です。しかし元々明るく映っている建物の密集地ではというとその判断は中々に難しいです。
この論文では観測された後方散乱の値を基に多数の領域（segment）に区分し、機械学習により各領域の浸水を判定し、建物の密集地も含めて、正解データである国土地理院提供の洪水範囲（赤線枠）を再現しています。学習とその評価については交差検証が行われています。
この論文では領域の区分する際に、２つのスケールにおいて領域成長法を用いています。ひとつは建物の密集地など地域的な特徴を捉える大きなスケールで、もうひとつは個々の建物等の特徴を捉える小さなスケールです。この大小2種類の領域区分を用いて、それぞれのスケールにおいて領域内の後方散乱の特徴や変化を数値化し、機械学習の入力としています。
具体的な数値や手法についてはここでは割愛します。
論文の図14を見ると、大小2種類のスケールの領域を組み合わせることで、機械学習により再現された浸水範囲が概ね正解の浸水範囲と重なっていることがわかります。
このように、論文では、2種類のスケールの領域を学習に入れることが重要であると強調していました。
同時に、今回は機械学習の入力データの前処理に多数の手作業が入っている点と、及び災害時と災害前の観測の間の物理的な変化（筆者注：例として建物が建つ/撤去される）や観測角度の違いがあった場合への対処が未考慮な点を今後の課題と挙げていました。
4. 簡単な分析
では、実際の洪水がどのように見えるのか、実際に手を動かして確認してみましょう。先に紹介した論文のように、建物が多いと洪水の判断は難しいので、広大な農地での例を見てみましょう。
今年、世界的なニュースとなった洪水のひとつに2022年末から2023年春にかけて米国カリフォルニア州の広範囲における事例があります。その内2023年3月に関してはICEYEのレポートが出ています。今回はこの一連の洪水の内、2023年1月上旬のサクラメント市周辺での洪水を取り上げます。
2023年1月12日の洪水の状況がYouTubeにありました。場所はサクラメント郊外、サクラメント空港の西、川を挟んだ西側です（次の図の左、OpenStreerMapで見る）。平時のこの辺りは地平線まで広がる畑です（20年以上前の記憶を呼び起こしながら書いてます）。
一方で、可視光では雲により地上の様子は判別できません。
結果を先にお見せすると、Sentinel-1（VV※）では次のように見えました。
※VVとは地面に対して垂直の向き（偏波）で発した電波の後方散乱を同じ垂直偏波で観測したということです。当記事末尾の用語集の偏波もご参考ください。
西経121.65度辺りのサクラメント市郊外の帯状の地域において大きくシグナルが下がり、緑色から青色に変化していることがわかります。また、先の動画が撮られた地点（洪水を横切る2本の高速道路の内、北の方の近く）周辺もシグナルが下がっていることがわかります。
データの取得と作図
それでは、この洪水の図（道路は除く）の作成方法を紹介します。
今回の分析ではSARの観測データに対し、緯度経度情報の簡単な付与（georeferencing）と切り出しのみを行います。地形による効果の補正や緯度経度に即したグリッドへの焼き直し（geocoding）は行いません。画角による歪み等の補正は最後の出力の際に matplotlib の機能を用いて、見た目のみ簡単な補正に留めます。
今回使用するSentinel-1のGround Range Detected (GRD) データは、ある程度のノイズ処理が済んでいる処理レベル1.5のデータです。その座標は衛星の観測方向に即した座標系（slant range coordinate）であり、緯度経度には沿っていません（geocodingは行われていません）。データはおおよそ 10m×10m の正方形のグリッドに焼き直されています。よって簡単に見るだけなら今回行う手順のいくつかは省略可能です（文末、余談をご参照ください）。またGoogle Earth Engineからはgeocodingまで済ませた処理レベル2のデータを取得できます。今回のようにただ見るだけならこちらの方が便利です。
ESAが提供する解析ライブラリのSNAPPYにより補正を行うこともできますが、Pythonのバージョン管理上の問題があるので、本記事ではよりシンプルなライブラリのみ使用してPythonで解析します。
手順０ 環境の構築（Anaconda）
手順１ SARの観測データの取得
手順２ 前処理（切り出し＋保存）
手順３ 描画
手順０ 環境の構築（Anaconda）
今回はできるだけ構築が簡単な環境で行います。
Anaconda でローカルに Python 3.9 の環境を作り、rasterio、scipy、matplotlib をインストールします。numpy と GDAL も自動的に含まれるはずです。Google Colaboratoryで構築する場合は !pip install で上記のライブラリをインストールすれば、以後のコードは同じです。
rasterio は GeoTIFF等のフォーマットの衛星データの分析に便利なライブラリです。rasterioの使い方については次の記事も参考にしてください。ただしこの記事での作図は配列のインデックス、つまり衛星が見た方向（slant range）に沿ったもので、緯度経度に合わせたグリッドに焼き直すgeoreferencingはしていません。
【宙畑】rasterioを使った分析
植物の分布から田畑の利用状況を推測する衛星データの利活用方法
手順１ SARの観測データの取得
今回は Alaska Satellite Facility Vertex から Sentinel-1 のデータを取得します。データの取得にはNASAのEarthObservationのアカウント登録（無料）が必要です。次の図を見ながらSentinel-1のデータ取得の手順を追っていきます。
まずは①Sentinel-1が選択されていることを確認し、②③欲しい領域を図示します。④⑤⑥次に時間を指定し、⑦検索します。⑧出てきたリストから良さそうなシーンを選び、⑨ダウンロードのカートから、⑩今回はGRDH（GRDの高解像度版）を選びダウンロードします。
また12日前（前の回帰）に同じ画角の観測があるはずです。④時間指定からの手順を繰り返し、同じ画角のGRDHデータをダウンロードします。
手順２ 前処理（切り出し+保存）
今回の作業で一番煩雑で難しいパートです。（コマンドから変換する方法もあります。文末の後日談１をご参照ください。）
特定領域を切り出すためには、衛星データの配列のインデックスと緯度経度を対応させなければなりません。インターネットには Sentinel-2 での rasterio の実装事例は多々見られますが、今回使用する Sentinel-1 とは座標系の保存方法が異なるために直接流用することはできません。Sentinel-2では衛星データの４角の緯度経度で制御されています（インスタンス変数 crs に格納されています）。一方でSentinel-1では衛星データの内側も含め多数のピクセルを基準点（ground control points, GCPs）として、それらの緯度経度を掲載しています。Sentinel-2は衛星の直下を見るため観測の歪みが少ないですが、Sentinel-1は斜め下を見るため観測に歪みが生じてしまうためかと思います。
まずはGeoTIFFファイルを rasterio で読み込みます。
import rasterio
import numpy as np
fn_zip =
fn_tiff =
data = rasterio.open("zip+file://"+fn_zip+"!"+fn_tiff)fn_zip には zipファイルのアドレスを、fn_tiff には zipファイルの中のTIFFファイルのアドレスを指定します。なおPython上で zipファイルの中のTIFFファイルのアドレスを探すには zipfileライブラリが便利です。
次に観測範囲の地理空間情報を取得します。観測範囲内のいくつかの代表点の地理空間情報から、観測範囲の各ピクセルの緯度・経度・高度を計算します。ここで紹介する手法の他にも手法がありますが、それについては本記事の末尾の後日談2をご覧ください。
読み込んだオブジェクト data のインスタンス変数 gcps に観測値の配列のいくつかの代表点（Ground Controll Points）の地理空間情報が入っています。
gcps, gcp_crs = data.gcps変数 gcp_crs には "WGS 84" のように文字列で座標系の情報が入っています。変数 gcps には基準点（GCPs）の位置情報（観測値の配列のインデックスと緯度経度と高度）が入っています。
変数 gcps は1次元の配列なので、2次元に直します。観測値の配列上において代表点は一定間隔に並んでいると仮定しています。不規則に並んでいた場合、このアルゴリズムはつかえません。その場合は gcpxx 等の要素数と gcps の要素数が異なりますので、必要なら確認してください。
gcp_idx = []
gcp_idy = []
for gg in gcps:
if gg.row not in gcp_idx:
gcp_idx.append(gg.row)
if gg.col not in gcp_idy:
gcp_idy.append(gg.col)
gcp_idx.sort()
gcp_idy.sort()
gcp_shp = (len(gcp_idx),len(gcp_idy))
print(gcp_shp)
gcpxx = np.zeros(gcp_shp)
gcpyy = np.zeros(gcp_shp)
gcpzz = np.zeros(gcp_shp)
for gg in gcps:
idx = gcp_idx.index(gg.row)
idy = gcp_idy.index(gg.col)
gcpxx[idx][idy] = gg.x
gcpyy[idx][idy] = gg.y
gcpzz[idx][idy] = gg.z次に代表点の座標を使った３次のスプライン補間を行い、配列の全要素の位置情報を構築します。高度情報は今回の表示では使いませんが、将来的な使用を考え、計算しておきます。
from scipy.interpolate import RectBivariateSpline
gcpInterpX = RectBivariateSpline(gcp_idx, gcp_idy, gcpxx)
gcpInterpY = RectBivariateSpline(gcp_idx, gcp_idy, gcpyy)
gcpInterpZ = RectBivariateSpline(gcp_idx, gcp_idy, gcpzz)
idx = np.arange(band.shape[0])
idy = np.arange(band.shape[1])
xx = gcpInterpX(idx, idy)
yy = gcpInterpY(idx, idy)
zz = gcpInterpZ(idx, idy)最後に観測値を取得し、切り出します。観測値は2次元の配列で入っています。シンプルに緯度経度で四隅を与えて、その範囲を含む最小の長方形の配列を切り出します。
band = data.read(1) # 観測値を取得
lons = [-122.15, -121.25]
lats = [38.2, 38.8]
msk = (lons[0] <= xx) * (xx <= lons[1]) * (lats[0] <= yy) * (yy <= lats[1])
trim_msk_x = np.any(msk, axis=0)
trim_msk_y = np.any(msk, axis=1)
band_trim = band[trim_msk_y,:][:,trim_msk_x]
xx_trim = xx[trim_msk_y,:][:,trim_msk_x]
yy_trim = yy[trim_msk_y,:][:,trim_msk_x]
zz_trim = zz[trim_msk_y,:][:,trim_msk_x]作図や解析の度にこの切り出し作業を行うのは手間なので、次のように fn_save で指定したnpzファイルに保存しておくと良いでしょう。
np.savez(fn_save, band=band_trim, lon=xx_trim, lat=yy_trim, alt=zz_trim)保存したデータは data = np.load(fn_save) で読み込みます。保存した各変数については band_trim = data["band"] のようにそれぞれ呼び出せます。
手順３ 描画
観測値は線形で保存されていますが、非常に大きくなるのでデシベル（dB、常用対数の10倍の値）に変換します。今回は vmin=8, vmax=3, vn=11 として、カラーマップ "gist_earth" を用いると良い具合の水の表現となります。
band_db = 10 * np.log10(band_trim)
levels = np.linspace(vmin, vmax, vn, endpoint=True)
fig, ax = plt.subplots(1,1, figsize=(8,5))
im = ax.contourf(xx_trim, yy_trim, band_db, levels, cmap="gist_earth")
cbar = fig.colorbar(im)
cbar.ax.set_title('dB')
ax.axis('equal')
plt.show()matplotlib の描画範囲の設定は各自で行ってください。切り抜いた大きさによっては描画に時間がかかる場合がありますが、その場合は配列の範囲を [::dd,::dd] として、データを間引いて表示してみてください。変数 dd は表示するピクセルの間隔です。また先述のように道路を追加すると見比べが楽になります。
余談 適切な座標設定を行わなかった場合
この座標の設定はデータによっては非常に重要になります。座標を設定しない場合の結果は、描画の contourf の最初の２つの引数 x_ql と y_ql を省略することで出ます。
今回使用した処理レベル1.5の GRDデータ ではあまり有難味がありません。衛星の観測方向の座標系ではありますが、グリッドが概ね正方形に整えられているためです。結果としては、緯度経度に合わせて回転や反転を行う程度です。
一方、今回は使用しなかった、より生のデータに近い処理レベル1の Single Look Complex (SLC) では大きな変化がもたらされます。これは Sentinel-1 の Interferometric Wide Swath (IWS) 観測モードではピクセルが 5 m × 20 mであり、正方形では無いからです。
また、例えば Ascending なら南から北、Descending なら北から南と言った具合に、衛星が観測した向きによって格納されているデータの順番と方位の関係が異なります。衛星の左右どちら側かによっても格納の順番と方位の関係が異なります。つまり、配列のインデックスのままだと上下左右反転した状態で表示される可能性があります。土地勘のある地域の観測データなら、閲覧の度に反転させることもできるかもしれませんが、得策ではないでしょう。
まとめ
今回行った処理は georeferencing と呼ばれる衛星データに緯度経度の地理空間情報を与える作業です。これを行うことで他の地理空間情報と組み合わせた解析が可能となります。分析の冒頭でお見せした道路情報と合わせた出力がその一例です。
光学観測では雲が邪魔で分からなかった地上の様子が、SARの観測では天候に関わらず観測できていることは重要です。洪水や土砂災害など大雨に起因する災害においては、往々にして現地の天候不順により観測が困難となります。SAR衛星はこの様な状況でも活用できます。
後日談１
SNAPPYやEoReaderだとSNAPに依存するので、よりお手軽に試せるように rasteria だけで処理するようにしました。しかし、ここで行った georeferencing を含む geocoding はGDALに含まれるソフトの1つ、gdalwarp でできるようです。筆者は動作を確認しきれていませんので、確認作業は読者にお任せ致します。
https://asf.alaska.edu/how-to/data-recipes/geocode-sentinel-1-with-gdal/
GDALについては次の宙畑の記事も参考にしてください。
GDALのおすすめコマンド5選と使用例【地理空間処理のいろは】
後日談２
Ground Control Pointsの地理空間情報から各ピクセルに座標を振り分けるgeoreferencingのさいに scipy.interpolate の RectBivariateSpline を使用しましたが、同ライブラリの griddata を使用することもできます。その場合、変数 gcps を2次元配列に変換する必要はなく、コードは短くなります。しかし、補間精度が3次から1次に下がり、手元環境での実行時間が5倍以上に伸びました。
from scipy.interpolate import griddata
gcpcords = np.array([[gg.row, gg.col] for gg in gcps])
gcpx, gcpy, gcpz = np.array([[gg.x, gg.y, gg.z] for gg in gcps]).T
idx = np.arange(band.shape[0])
idy = np.arange(band.shape[1])
idyy, idxx = np.meshgrid(idy, idx)
print(np.shape(idxx))
xx = griddata(gcpcords, gcpx, (idxx, idyy), method="linear")
yy = griddata(gcpcords, gcpy, (idxx, idyy), method="linear")
zz = griddata(gcpcords, gcpz, (idxx, idyy), method="linear")SAR 語彙録
geocoding：地理空間座標系に即したグリッドに焼き直すこと。同じ語句がGIS界隈では建物の住所などに緯度経度の情報を付与することを指すので注意。
georeferencing：衛星データに緯度経度などの地理空間座標情報を付与すること。
ground control point (GCP)：地上基準点。緯度経度が判明している地点のこと。
Ground Range Detected (GRD)：multi-looking処理を行い、地球の楕円体モデルに投影するgeoreferencingを行うこと。処理後のデータは後方散乱の強度のみ持つ（位相情報は失われる）。データのグリッドの形は概ね正方形となる。レベル1.5とも呼ばれる。
multi-looking処理：空間平均（spatial averaging）や間引き（decimation）によりスペックルノイズの緩和を行うこと。GRDの生成の際に使用される。
Single Look Complex (SLC)：レーダー信号を1枚の絵となるように加工すること。処理後のデータは後方散乱の強度と位相の情報を持ち、衛星の軌道と姿勢の情報から推察される地理空間情報が適用される。データのグリッドの形は衛星の観測モードにより異なる。レベル1とも呼ばれる。
slant range coordinate：衛星の観測方向に準じた座標系。SLCやGRDのデータの座標系。
terrain correction：衛星データを地球の地形の凹凸に合わせ、信号強度の標準化を行うこと。
後方散乱：SAR衛星が発信した電波が地上の物体に当たり、SAR衛星に返ってくること。
偏波：電波（電磁波）の電界の向き。詳しくは宙畑記事「SAR（合成開口レーダ）のキホン～事例、分かること、センサ、衛星、波長～」をご覧ください。
【図解】衛星データの前処理とは~概要、レベル別の処理内容と解説~
・Rüdiger Gens “Terrain correction of SAR imagery,” presentation slides
