はじめに
衛星データ分析シリーズの続きです。今回は Earth Big Data(EBD)SARトレーニングワークショップのChapter 3、変化検出(Change Detection)を実装しました。
元の教材である、THE SAR HANDBOOK Chapter 3 Training Module Part3 - Change Detection with Time Series Metrics and Log Ratio Methodは、手元に置いたローカルのSentinel-1 VRTファイルを読み込む形式で書かれています。しかし私たちの環境では大量のデータをローカルに置くのが難しく、またGoogle Earth Engine(GEE)を使えばSentinel-1データにクラウド上で直接アクセスできるので、GEEを使った実装に書き換えました。
SARを使う利点
今回扱うのは可視光画像ではなくSAR(合成開口レーダー)です。SARの最大の強みは雲の影響を受けないこと。可視光センサーでは雲に覆われると地表が見えなくなりますが、SARはマイクロ波を使うため雲を貫通して観測できます。熱帯や梅雨の時期でも安定してデータが取れるのが大きな魅力です。
さらに、SARのバックスキャッター(後方散乱)値は植生の体積散乱や地表の粗さに敏感で、森林の伐採・洪水・農地変化などを捉えるのに向いています。
利用環境・データ
- 環境: ローカルJupyter Notebook(Python 3.12、仮想環境)
- データ: Sentinel-1 GRD、VVポラリゼーション
- 対象地域: 西アフリカ バイオマスサイト(ROI: 経度2.5〜2.6°、緯度10.5〜10.6°)
- 期間: 2015年3月〜2015年11月(19シーン)
実装した変化検出手法は3つです:
- 時系列メトリクス - パーセンタイル差分閾値法
- 時系列変動係数(Coefficient of Variation)閾値法
- 対数比(Log Ratio)法
Pythonコードと実行結果
セル1: ライブラリのインポート
import os, sys
import datetime
from osgeo import gdal
%matplotlib inline
import matplotlib.pylab as plt
import matplotlib.patches as patches
from skimage import exposure
import numpy as np
import pandas as pd
import subprocess
import sys
try:
import geemap
except ImportError:
subprocess.check_call([sys.executable, "-m", "pip", "install", "-U", "geemap"])
import geemap
解説: 画像処理に osgeo.gdal、数値計算に numpy、可視化に matplotlib、ヒストグラム平坦化に skimage.exposure を使います。元の教材では import gdal と書かれていましたが、現在のGDALは from osgeo import gdal が正しいインポート方法です。
セル2: Google Earth Engineの初期化
import ee
try:
ee.Initialize()
except Exception as e:
ee.Authenticate()
ee.Initialize(project='<プロジェクト名>')
解説: GEEへの接続は ee.Initialize() で行います。初回は認証が必要で、ee.Authenticate() を呼ぶとブラウザが開きGoogleアカウントでのログインが求められます。project= には自分のGEEプロジェクトIDを指定します。
ハマりポイント: ローカル環境での認証は一度通れば
~/.config/earthengine/にキャッシュされますが、キャッシュが壊れると再認証が必要です。Colabで実行するのが一番手軽です。
セル3: ROIの定義とSentinel-1コレクションの取得
roi = ee.Geometry.Rectangle([2.5, 10.5, 2.6, 10.6])
s1_collection = (ee.ImageCollection('COPERNICUS/S1_GRD')
.filterBounds(roi)
.filterDate('2015-03-22', '2015-11-30')
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.select('VV'))
解説: GEEの COPERNICUS/S1_GRD コレクションからSentinel-1データを取得します。フィルタリングのポイントは以下の通り:
-
filterBounds(roi): 対象地域のデータだけに絞る -
filterDate(...): 取得期間を指定 -
instrumentMode='IW': Interferometric Wide Swath(干渉広域)モード。陸域観測の標準モード -
select('VV'): VVポラリゼーションのみ選択。VVは垂直送信・垂直受信で、森林変化に感度が高い
元の教材ではVVとVH(クロスポラリゼーション)の両方を使うオプションもありましたが、今回はVVに絞りました。
セル4: 出力用ROIの定義
center_point = [2.595, 10.595]
roi = ee.Geometry.Point(center_point).buffer(10000) # 半径10kmの円
# プレビュー/単一シーン処理用に1シーン目をROIでクリップ
image_clipped = s1_collection.first().clip(roi)
解説: 解析範囲(ROI)を定義し、1シーン目をクリップします。GEEはクラウド上でデータを処理できるため、ローカルへのファイルダウンロードは不要です。出力GeoTIFFに付ける投影情報も、後述のセル15でGEEのメタデータから直接取得します。
セル5: 1シーン目のプレビュー(GEEから直接)
preview = geemap.ee_to_numpy(image_clipped) # dBスケール
preview2d = preview[:, :, 0] if preview.ndim == 3 else preview
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(preview2d, cmap='gray')
ax.set_title('Sentinel-1 VV (dB) - first scene (GEE direct)')
plt.colorbar(im, ax=ax, label='backscatter (dB)')
plt.show()
print(f"shape: {preview2d.shape}")
print(f"range: {np.nanmin(preview2d):.2f} ~ {np.nanmax(preview2d):.2f} dB")
解説: 1シーン目のVVバンドを geemap.ee_to_numpy() でNumPy配列に変換し、グレースケールで表示します。ファイルを介さずGEEから直接取得しているのがポイントです。
値が -26 〜 4 付近に分布する点に注目してください。これはdB(デシベル)スケールです。GEEのSentinel-1 GRDはデフォルトでdBスケールで出力されます。元の教材ではDN値(デジタルナンバー)から振幅を読んでいましたが、GEE版では最初からdBで来るため変換式が異なります。
セル6: 時系列スタックの作成とdB→パワー変換
dates = s1_collection.aggregate_array('system:time_start').getInfo()
dates = [datetime.datetime.fromtimestamp(d/1000).date() for d in dates]
tindex = pd.DatetimeIndex(dates)
print(f"Number of images in collection: {len(tindex)}")
print(f"Date range: {tindex[0].date()} to {tindex[-1].date()}")
# スタック化(サーバーサイドで効率的に処理)
stacked_image = s1_collection.toBands()
stacked_image_clipped = stacked_image.clipToBoundsAndScale(
geometry=roi, maxDimension=512
)
# バンド数を取得
num_bands = stacked_image.bandNames().length().getInfo()
# NumPy配列に変換
rasterDB = geemap.ee_to_numpy(stacked_image_clipped)
# dB → パワー変換
rasterPwr = 10.**(rasterDB / 10.)
# GEEは (height, width, bands) 形式で返すので転置
if rasterPwr.shape[2] == num_bands:
rasterPwr = np.transpose(rasterPwr, (2, 0, 1))
print(f"Final array shape: {rasterPwr.shape}")
解説: ここが元の教材との最大の違いです。
元の教材(ローカルVRTファイル)の場合:
# DN値 → パワー(振幅の2乗)
CF = np.power(10., -8.3)
rasterPwr = np.power(rasterDN, 2.) * CF
GEE版の場合:
# dB → パワー(標準的なdB逆変換)
rasterPwr = 10.**(rasterDB / 10.)
GEEのSentinel-1はすでに 10 * log10(パワー) の計算が済んだdB値で出力されているため、逆変換式 10^(dB/10) でパワー値に戻します。数学的には正しい変換です。
また、s1_collection.toBands() でコレクション全体をマルチバンド画像に変換し、clipToBoundsAndScale(maxDimension=512) で512×512ピクセルにリスケールしています(GEEのリクエストサイズ制限50MBに収めるため)。
実行結果:
Number of images in collection: 19
Date range: 2015-03-23 to 2015-11-30
Final array shape: (19, 512, 512)
19時期分のデータが揃いました。
セル7: 多時期カラーコンポジット表示
rgb_bands = (0, int(num_bands/2), num_bands-1) # 最初・中間・最後
rgb_idx = np.array(rgb_bands)
rgb = np.dstack((rasterPwr[rgb_idx[0]], rasterPwr[rgb_idx[1]], rasterPwr[rgb_idx[2]]))
rgb_dates = (tindex[rgb_idx[0]].date(),
tindex[rgb_idx[1]].date(),
tindex[rgb_idx[2]].date())
rgb_stretched = rgb.copy()
for i in range(rgb_stretched.shape[2]):
band = rgb_stretched[:,:,i]
valid_mask = band != 0
if valid_mask.any():
v_min, v_max = band[valid_mask].min(), band[valid_mask].max()
band_norm = (band - v_min) / (v_max - v_min) if v_max > v_min else band
rgb_stretched[:,:,i] = exposure.equalize_hist(band_norm)
fig, ax = plt.subplots(1, 2, figsize=(16, 8))
fig.suptitle('Multi-temporal Sentinel-1 backscatter image R:{} G:{} B:{}'
.format(rgb_dates[0], rgb_dates[1], rgb_dates[2]))
ax[0].imshow(rgb)
ax[0].set_title('Unstreched')
ax[1].imshow(rgb_stretched)
ax[1].set_title('Histogram Equalized')
解説: 時系列の「最初・中間・最後」の3時期をRGB合成して表示します。変化があった場所は色がついて見え、変化のない場所はグレーになります。ヒストグラム平坦化(equalize_hist)を適用することでコントラストが改善されます。元の教材と比べ、0値マスクの処理を丁寧に行っている点が改良点です。
セル8: 時系列メトリクスの計算
def timeseries_metrics(raster, ndv=0):
if not raster.dtype.name.find('float') > -1:
raster = raster.astype(np.float32)
if ndv != np.nan:
raster[np.equal(raster, ndv)] = np.nan
tsmetrics = {}
rperc = np.nanpercentile(raster, [5, 50, 95], axis=0)
tsmetrics['mean'] = np.nanmean(raster, axis=0)
tsmetrics['max'] = np.nanmax(raster, axis=0)
tsmetrics['min'] = np.nanmin(raster, axis=0)
tsmetrics['range'] = tsmetrics['max'] - tsmetrics['min']
tsmetrics['median'] = rperc[1]
tsmetrics['p5'] = rperc[0]
tsmetrics['p95'] = rperc[2]
tsmetrics['prange'] = rperc[2] - rperc[0]
tsmetrics['var'] = np.nanvar(raster, axis=0)
tsmetrics['cov'] = tsmetrics['var'] / tsmetrics['mean']
return tsmetrics
metrics = timeseries_metrics(rasterPwr, ndv=0)
解説: 各ピクセルについて19時期分の時系列統計量を計算します。計算されるメトリクスは10種類:
| メトリクス | 意味 |
|---|---|
mean |
時系列平均バックスキャッター |
median |
中央値(外れ値に頑健) |
max/min |
最大・最小値 |
range |
最大値 − 最小値(変化幅) |
p5/p95 |
5・95パーセンタイル |
prange |
p95 − p5(外れ値に頑健な変化幅) |
var |
分散 |
cov |
変動係数(= 分散/平均) |
ハマりポイント(解決済み): 元の教材では
rasterPwr.filled(np.nan)を使っていましたが、これはNumPyのマスク配列(MaskedArray)専用のメソッドです。GEE版ではgeemap.ee_to_numpy()が通常のNumPy配列を返すのでfilled()は使えません。timeseries_metrics(rasterPwr, ndv=0)として呼び出すことで、0値をNaNとして扱い正しく計算できます。
実行結果(一部):
dict_keys(['mean', 'max', 'min', 'range', 'median', 'p5', 'p95', 'prange', 'var', 'cov'])
セル9: メトリクスの可視化
metric_keys = ['mean', 'median', 'max', 'min',
'p95', 'p5', 'range', 'prange', 'var', 'cov']
fig = plt.figure(figsize=(16, 40))
idx = 1
for i in metric_keys:
ax = fig.add_subplot(5, 2, idx)
data = metrics[i].flatten()
data_valid = data[~np.isnan(data)]
vmin = np.percentile(data_valid, 2)
vmax = np.percentile(data_valid, 98)
ax.imshow(metrics[i], vmin=vmin, vmax=vmax, cmap='gray')
ax.set_title(f'{i.upper()} (min={vmin:.2e}, max={vmax:.2e})')
ax.axis('off')
idx += 1
解説: 10種類のメトリクスを一覧で可視化します。元の教材では固定のスケール(vmin=0.0001, vmax=0.3)を使っていましたが、GEE版では各メトリクスのデータから自動的に2〜98パーセンタイルを計算してスケールを決めています。データに合わせた表示になるため見やすくなっています。
セル10: 変化検出① パーセンタイル差分閾値法
plt.hist(metrics['prange'].flatten(), bins=100, range=(0, 0.3))
_=plt.axvline(0.27, color='red')
thres = 0.27
mask = metrics['prange'] < thres
maskpdiff = ~mask
plt.figure(figsize=(8, 8))
plt.imshow(mask, cmap='gray')
plt.legend(['$p_{95} - p_5 > 0.27$'], loc='center right')
_=plt.title('Threshold Classifier on Percentile Difference ($P_{95} - P_5 > 0.27$)')
解説: 時系列のパーセンタイル差分(prange = P95 − P5)をヒストグラムで確認し、閾値を設定します。変化が大きいピクセルは prange が大きくなるため、閾値を超えたピクセルを変化点として検出します。最大−最小(range)と違い、上下5%の外れ値を除いた差分なので、ノイズや単発のスパイクに頑健な変化幅が得られます。
セル11: 変化検出② 変動係数閾値法
plt.hist(metrics['cov'].flatten(), bins=100, range=(0, 0.05))
_=plt.axvline(0.025, color='red')
thres = 0.025
mask = metrics['cov'] < thres
maskcv = ~mask
plt.figure(figsize=(8, 8))
plt.imshow(mask, cmap='gray')
_=plt.title('Threshold Classifier on Time Series Coefficient of Variation')
解説: 変動係数(COV = 分散/平均)を閾値で分類します。COVが高いピクセルは時系列のばらつきが大きく、何らかの変化が起きている可能性が高いです。COVは平均値で正規化されているため、バックスキャッターの絶対値に左右されにくいのが特徴です。
セル12: 時系列の確認とログ比法の日付選定
tsmean = 10 * np.log10(np.nanmean(rasterPwr, axis=(1, 2)))
ts = pd.Series(tsmean, index=tindex)
for i in range(len(ts)):
print(i, ts.index[i].date(), ts[i])
実行結果:
0 2015-03-23 -11.497
1 2015-04-04 -11.706
...
9 2015-07-21 -8.340
...
18 2015-11-30 -11.543
解説: np.nanmean(rasterPwr, axis=(1,2)) で各時期の画像全体の平均パワー値を計算し、dBに変換して時系列プロットします。axis=(1,2) はピクセル方向(height・width)で集計するという意味で、結果は (19,) の1次元配列になります。
乾季(3〜4月・11月)はバックスキャッターが低く(-11〜-12 dB)、雨季(7〜9月)は高い(-8〜-7 dB)傾向が読み取れます。これは西アフリカの植生の季節変動を捉えています。
セル13: 変化検出③ 対数比(Log Ratio)法
Xr = rasterPwr[1] # 参照画像: 2015-04-04 (-11.7 dB)
Xi = rasterPwr[15] # 比較画像: 2015-10-01 (-8.1 dB)
r = np.log10(Xi / Xr)
fig, ax = plt.subplots(2, 1, figsize=(8, 16))
ax[0].imshow(np.abs(r), vmin=0, vmax=0.3, cmap='gray')
_=ax[1].hist(r.flatten(), bins=100, range=(-0.4, 0.4))
解説: 2つの時期の対数比を計算します。
$$r = \log_{10}\left(\frac{X_i}{X_r}\right)$$
$r > 0$: バックスキャッター増加(植生成長・湿潤化)、$r < 0$: 減少(伐採・乾燥化)。元の教材では同じ季節の異なる年を比較(2015年 vs 2017年)していましたが、GEE版では2015年内の乾季(4月)と雨季(10月)を比較しています。これは2015年のみのデータを取得しているためで、目的が「変化検出」から「季節変動の検出」に変わっている点に注意が必要です。
セル14: 対数比の変化ピクセル抽出
stddev = np.std(r)
thres = 3 * stddev
mask = np.logical_and(r > -1*thres, r < thres)
masklr = ~mask
fig, ax = plt.subplots(figsize=(8, 16))
ax.imshow(mask, cmap='gray')
_=ax.set_title('Log Ratio Classifier of the 2015-04-04 - 2015-10-01 Log Ratio Images')
解説: 対数比の平均から ±3σ(標準偏差の3倍)を超えたピクセルを変化点として分類します。統計的に「通常の変動範囲」を超えた変化を検出するシンプルで頑健な方法です。
$$cp_{x,y} = (r_{x,y} < \bar{r} - 3\sigma_r) \text{ or } (r_{x,y} > \bar{r} + 3\sigma_r)$$
セル15: 投影情報の取得(GEEメタデータから)
from osgeo import osr
n_bands, height, width = rasterPwr.shape
proj_info = stacked_image_clipped.projection().getInfo()
crs_code = proj_info.get('crs', 'EPSG:4326')
t = proj_info.get('transform', None)
# CRSコード → WKT
srs = osr.SpatialReference()
srs.SetFromUserInput(crs_code)
proj = srs.ExportToWkt()
def _looks_degenerate(tr):
return (tr is None or len(tr) != 6 or
(abs(tr[0]) == 1 and abs(tr[4]) == 1 and tr[2] == 0 and tr[5] == 0))
if not _looks_degenerate(t):
# GEE transform [a,b,c,d,e,f] → GDAL GeoTransform (xoff,xscale,xshear,yoff,yshear,yscale)
geotrans = (t[2], t[0], t[1], t[5], t[3], t[4])
else:
# フォールバック: ROIの外接矩形を配列サイズ(512×512)へ割り当てて逆算
coords = ee.Geometry(roi).bounds().transform(crs_code, 1).coordinates().getInfo()[0]
xs = [c[0] for c in coords]; ys = [c[1] for c in coords]
geotrans = (min(xs), (max(xs)-min(xs))/width, 0.0,
max(ys), 0.0, -(max(ys)-min(ys))/height)
print(f"CRS: {crs_code}")
print(f"GeoTransform: {geotrans}")
解説: GeoTIFF出力に必要な投影座標系(CRS)と地理変換(GeoTransform)をGEEのメタデータから直接組み立てます。
ポイントは、解析に使った stacked_image_clipped(512×512にリスケール済み)の投影情報をそのまま使うことです。projection().getInfo() が具体的なtransformを返せばそれを採用し、退化した値(単位行列)しか返さない場合はROIの外接矩形を実際の配列サイズ(512×512)で割ってピクセルサイズを逆算します。こうすることで、配列と地理座標が必ず一致します。
セル16: 各メトリクスをGeoTIFFとして出力
def CreateGeoTiff(Name, Array, DataType, NDV, bandnames=None,
ref_image=None, GeoT=None, Projection=None):
if len(Array.shape) == 2:
Array = np.array([Array])
if ref_image != None:
refimg = gdal.Open(ref_image)
GeoT = refimg.GetGeoTransform()
Projection = refimg.GetProjection()
driver = gdal.GetDriverByName('GTIFF')
Array[np.isnan(Array)] = NDV
DataSet = driver.Create(Name, Array.shape[2], Array.shape[1],
Array.shape[0], DataType)
DataSet.SetGeoTransform(GeoT)
DataSet.SetProjection(Projection)
for i, image in enumerate(Array, 1):
DataSet.GetRasterBand(i).WriteArray(image)
DataSet.GetRasterBand(i).SetNoDataValue(NDV)
DataSet.FlushCache()
return Name
output_dir = 'sentinel1_geotiff_output'
metrics_dir = os.path.join(output_dir, 'metrics')
os.makedirs(metrics_dir, exist_ok=True) # 親のoutput_dirも自動作成される
output_files = []
for metric_name, metric_data in metrics.items():
output_file = os.path.join(metrics_dir, f'metric_{metric_name}.tif')
CreateGeoTiff(
Name=output_file, Array=metric_data,
DataType=gdal.GDT_Float32, NDV=np.nan,
bandnames=[metric_name],
GeoT=geotrans, Projection=proj # GEE由来の投影情報を直接渡す
)
output_files.append(output_file)
print(f"✓ {metric_name:12s} → {output_file}")
解説: 10種類のメトリクスをそれぞれ独立したGeoTIFFとして書き出します。
CreateGeoTiff() 関数の処理の流れは以下の通りです:
-
len(Array.shape) == 2のときnp.array([Array])で先頭に1次元追加し、2次元配列も3次元(バンド, 高さ, 幅)として統一的に扱う -
driver.Create(Name, 幅, 高さ, バンド数, DataType)でGeoTIFFの器を作成。引数の順序がArray.shape[2](幅)→Array.shape[1](高さ)であることに注意 -
SetGeoTransform()/SetProjection()で地理参照を埋め込む。ここに渡しているのがセル15でGEEから取得したgeotrans/projです -
Array[np.isnan(Array)] = NDVでNaNをNoData値に置換し、各バンドへWriteArray()で書き込み
呼び出し側では GeoT=geotrans, Projection=proj を直接渡しているため、GEEのメタデータだけでGeoTIFFを書き出せます(gdal は書き出し処理にのみ使用)。ループで metrics の全エントリを回し、metric_<名前>.tif として10ファイルを出力します。
実行結果:
✓ mean → sentinel1_geotiff_output/metrics/metric_mean.tif
✓ max → sentinel1_geotiff_output/metrics/metric_max.tif
...(10ファイル出力)
セル17: VRT(Virtual Raster Table)の作成
vrt_file = os.path.join(output_dir, 'metrics.vrt')
cmd = ['gdalbuildvrt', '-separate', '-overwrite', '-vrtnodata', 'nan',
vrt_file] + output_files
result = subprocess.run(cmd, capture_output=True, text=True)
if result.returncode == 0:
print(f"VRT file: {vrt_file}")
for i, metric_name in enumerate(metrics.keys(), 1):
print(f" Band {i}: {metric_name}")
解説: 10個のGeoTIFFを1つのVRTファイルにまとめます。VRTは元ファイルへの参照を持つだけの軽量ファイル(約6KB)で、QGISで開くとバンドを切り替えながら全メトリクスを確認できます。
実行結果:
Band 1: mean
Band 2: max
...
Band 10: cov
QGISでこのVRTを開き、OpenStreetMapと重ね合わせることで変化検出結果を地図上で確認できます。
まとめ
EBDのSARハンドブックChapter 3をGoogle Earth Engine(Sentinel-1 GRD)で再実装しました。
うまくいったこと
- GEEの
COPERNICUS/S1_GRDから19シーンを取得し、3つの変化検出手法(パーセンタイル差分、変動係数、対数比)を全て実装できた -
toBands()+clipToBoundsAndScale()によるサーバーサイドスタック処理でデータ取得が効率化できた - dB→パワー変換(
10**(dB/10))とtimeseries_metrics(rasterPwr, ndv=0)の組み合わせが正しく動作することを確認した - 出力GeoTIFFの投影情報をGEEメタデータから直接取得することで、中間ファイルを介さずに解析配列(512×512)と地理座標を一致させられた
- GeoTIFF出力とVRT作成まで完走し、QGISで可視化できた
課題・気づき
- 比較対象が1年内のみ: 元の教材は異なる年の同季節を比較するが、GEE版は2015年内の乾季↔雨季比較になった。変化検出というより季節変動の可視化になっている
- GEE認証の不安定さ: ローカルでの認証キャッシュが壊れることがあった。Colabで実行するのが現状最も安定している(これについては別途調査中)
- VHバンドの活用: 元の教材はVHとVVの比較も示唆していた。VHを加えることでよりリッチな変化検出ができるはず
次のステップは、ROIを日本国内(鳥取砂丘や釧路湿原など)に変えて、日本の植生や土地被覆変化に対してどう機能するかを試してみたいと思います。
参考:EBD SAR Training Workshop for Forest Applications, Josef Kellndorfer Ph.D., Earth Big Data LLC, 2018







