はじめて衛星画像を触ってみたので、ノートを兼ねて分析手順を書きます。

画像のソース

ALOS-2（だいち2号）のサンプル画像を使います。今回は1月に配布された能登半島のL2.1(GeoTIFF)データを使って分析します。L1.1のデータを使うと、バイナリファイルの処理といくつかの補正処理などが必要で、初心者向けではないので、一定程度の処理済みのL2.1を使います。

参考するソースコード

今回はTellusの記事を参考して、特に「2時期のSAR強度画像を用いたRGB疑似カラー合成画像による浸水範囲の推定」の部分の手順を再現してみます。

ファイルの中身を確認

ダウンロードしたzipは画像データのtifといくつか付随するファイルがあります。特にテキストファイルの中に重要なサマリー情報が入っているので、そちらで座標範囲・撮影日時・スペックなども手軽に確認できます。

サイトの情報と照り合わせるとよりわかりやすいと思います。今回比較したい画像は2024/01/02（地震後）と2023/06/06（地震前）で、

方向は「Desc」（降交）なので、画像の並び方は下から上に
方向がわかったので、つまり「Frame」の部分は下から2820 -> 2830 -> 2840の順に
地図を参考して、今回比較したいのは輪島市で、Frameとしては2830となる
「Beam」はU2で、つまり高分解能観測モードで、実際は3m分解能で、無料で取得できる衛星画像の場合だとかなり貴重かな
画像データはHH偏波の1枚だけ、仕様資料を読むと他の種類もあるけど、今回提供されるのはHH偏波だけ

画像データを読んで確認

地理空間処理用のgdalを使ってtifを読み込む

import numpy as np
from osgeo import gdal, gdalconst, gdal_array

src_366 = gdal.Open('./IMG-HH-ALOS2487932830-230606-UBSL2.1GUD.tif', gdalconst.GA_ReadOnly)
src_412 = gdal.Open('./IMG-HH-ALOS2518982830-240102-UBSL2.1GUD.tif', gdalconst.GA_ReadOnly)

それぞれの属性取得とArrayに変換

RasterXSize：水平方向ピクセル数
RasterYSize：鉛直方向ピクセル数
RasterCount：バンド数

print(src_366.RasterXSize, src_366.RasterYSize) # 29100 32800
print(src_412.RasterXSize, src_412.RasterYSize) # 29080 32800
print(src_366.RasterCount, src_412.RasterCount) # 1

a_366 = src_366.GetRasterBand(1).ReadAsArray()
a_412 = src_412.GetRasterBand(1).ReadAsArray()

同じFrameの衛星画像だけど、サイズがほんの少し違うので、処理する時に気を付かないといけません。

画像を描画してみる

plt.imshow(a_412,vmax=10000,cmap='gray')
plt.show

地図と照り合わせると位置が大体把握できると思います。衛星画像は必ず長方形の中に納まるではなくて、斜めの画像範囲が多いので、切り抜きの時に注意が必要です。

座標を指定して画像を切り抜き

raw画像のサイズが大きくて、分析しにくいので、輪島市だけを切り抜けます。
座標は約37.37, 136.88から37.42, 136.93までとして、1枚目の切り抜きのソースコードはこちらになります。
ソースコードはGPTで生成されて少し編集を入れました。わりと難易度が高い処理でも生成できました。

from osgeo import gdal, osr

src_srs = osr.SpatialReference()
src_srs.ImportFromWkt(src_366.GetProjectionRef())

tgt_srs = osr.SpatialReference()
tgt_srs.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(src_srs, tgt_srs)

geotransform = src_366.GetGeoTransform()

width = src_366.RasterXSize
height = src_366.RasterYSize

ulx, uly = geotransform[0], geotransform[3]
lrx = ulx + width * geotransform[1]
lry = uly + height * geotransform[5]

# Transform corners to lat/lon
ul_lon, ul_lat, _ = transform.TransformPoint(ulx, uly)
lr_lon, lr_lat, _ = transform.TransformPoint(lrx, lry)

desired_ul_lon, desired_ul_lat = 136.88, 37.42
desired_lr_lon, desired_lr_lat = 136.93, 37.37

inverse_transform = osr.CoordinateTransformation(tgt_srs, src_srs)

desired_ulx, desired_uly, _ = inverse_transform.TransformPoint(desired_ul_lat, desired_ul_lon)
desired_lrx, desired_lry, _ = inverse_transform.TransformPoint(desired_lr_lat, desired_lr_lon)

output_path = 'cropped_src_366.tif'
gdal.Translate(output_path, src_366, projWin=[desired_ulx, desired_uly, desired_lrx, desired_lry])

切り抜けた画像を確認

c_src_366 = gdal.Open('./cropped_src_366.tif', gdalconst.GA_ReadOnly)
c_a_366 = c_src_366.GetRasterBand(1).ReadAsArray()
plt.imshow(c_a_366,vmax=10000,cmap='gray')

画像を見ると、輪島市とその周辺です。

print(c_src_366.RasterXSize, c_src_366.RasterYSize) # 1815 2183
print(c_src_412.RasterXSize, c_src_412.RasterYSize) # 1815 2183

もう1回サイズを確認すると、サイズは一致しています。
ただしまだ少し大きいので、さらに切り抜けを行います。

サイズを指定して画像を切り抜き

gdal.Translate(
  'tran_cropped_src_366.tif',
  c_src_366,
  outputType=gdal.GDT_Byte,
  scaleParams=[[0,10000]],
  srcWin=[250,250,1000,1000] # minX,minY,deltaX,deltaY
)
gdal.Translate(
  'tran_cropped_src_412.tif',
  c_src_412,
  outputType=gdal.GDT_Byte,
  scaleParams=[[0,10000]],
  srcWin=[250,250,1000,1000] # minX,minY,deltaX,deltaY
)

再度切り抜けた画像を確認します。

t_c_src_366 = gdal.Open('./tran_cropped_src_366.tif', gdalconst.GA_ReadOnly)
t_c_a_366 = t_c_src_366.GetRasterBand(1).ReadAsArray()
plt.imshow(t_c_a_366,cmap='gray')

これで大体の輪島市の部分を切り出しました。

疑似カラー合成画像で比較

最後は画像の描画と比較の部分です。

上記の「2時期のSAR強度画像を用いたRGB疑似カラー合成画像による浸水範囲の推定」記載してあるソースコードを参照して、地震前・地震後・比較の3画像を横に並んで出します。

import matplotlib.image as mpimg

img_before_8bit = gdal.Open('tran_cropped_src_366.tif', gdalconst.GA_ReadOnly)
array_before_8bit = img_before_8bit.ReadAsArray()
Xsize=img_before_8bit.RasterXSize
Ysize=img_before_8bit.RasterYSize
dtype=gdal.GDT_Byte
band=1
outfile = 'RGBcomposite.tif'
outRGB= gdal.GetDriverByName('GTiff').Create(outfile, Xsize, Ysize, band, dtype)
outRGB.SetProjection(img_before_8bit.GetProjection())
outRGB.SetGeoTransform(img_before_8bit.GetGeoTransform())
outRGB.GetRasterBand(1).WriteArray(array_before_8bit)
outRGB.FlushCache() 
image1 = mpimg.imread(outfile)

img_after_8bit = gdal.Open('tran_cropped_src_412.tif', gdalconst.GA_ReadOnly)
array_after_8bit = img_after_8bit.ReadAsArray()
Xsize=img_after_8bit.RasterXSize
Ysize=img_after_8bit.RasterYSize
dtype=gdal.GDT_Byte
band=1
outfile = 'RGBcomposite.tif'
outRGB= gdal.GetDriverByName('GTiff').Create(outfile, Xsize, Ysize, band, dtype)
outRGB.SetProjection(img_after_8bit.GetProjection())
outRGB.SetGeoTransform(img_after_8bit.GetGeoTransform())
outRGB.GetRasterBand(1).WriteArray(array_after_8bit)
outRGB.FlushCache() 
image2 = mpimg.imread(outfile)

img_before_8bit = gdal.Open('tran_cropped_src_366.tif', gdalconst.GA_ReadOnly)
img_after_8bit = gdal.Open('tran_cropped_src_412.tif', gdalconst.GA_ReadOnly)
array_before_8bit = img_before_8bit.ReadAsArray()
array_after_8bit = img_after_8bit.ReadAsArray()
Xsize=img_before_8bit.RasterXSize
Ysize=img_before_8bit.RasterYSize
dtype=gdal.GDT_Byte
band=3
outfile = 'RGBcomposite.tif'
outRGB= gdal.GetDriverByName('GTiff').Create(outfile, Xsize, Ysize, band, dtype)
outRGB.SetProjection(img_before_8bit.GetProjection())
outRGB.SetGeoTransform(img_before_8bit.GetGeoTransform())
outRGB.GetRasterBand(1).WriteArray(array_before_8bit)
outRGB.GetRasterBand(2).WriteArray(array_after_8bit)
outRGB.GetRasterBand(3).WriteArray(array_after_8bit)
outRGB.FlushCache() 
image3 = mpimg.imread(outfile)

fig, axes = plt.subplots(1, 3, figsize=(21, 7))
axes[0].imshow(image1)#, cmap="gray")
axes[0].set_title("before")
axes[1].imshow(image2)#, cmap="gray")
axes[1].set_title("after")
axes[2].imshow(image3)#, cmap="gray")
axes[2].set_title("compare")
plt.show()

比較結果

シアン色と赤色の部分を注目してください。シアン色は散乱体の出現で、赤色は散乱体の消失となります。

一番目立つのは、左真ん中のシアン色とやや真ん中下の方のシアン色で、前者は海岸なので、おそらく地震で上昇した地盤と見られます。後者は市内の学校の校庭なので、おそらく漂流物が堆積した結果と見られます。

赤色の部分が少ないけど、右の方の海岸の部分少し見られます。推測ですが、こちらは浸水によって強度が小さくなると考えられます。