イントロ
今回やってみたのは光学衛星画像の分析となります。
偶々衛星画像分析の記事を書いていますので、ぜひ見てみてください。
ALOS-2 / PALSAR-2のSAR画像で災害状況を分析
https://qiita.com/leolui2013/items/c64957ce85d8c23bb12e
雲画像予測をチャレンジしてみた
https://qiita.com/leolui2013/items/dd4c4110813df6c1a620
概要
衛星画像の種類が様々ありますが、人間にとって一番わかりやすいのはやはり光学画像です。
ただし、無料しかも解像度が高い光学画像は中々手に入れることが難しいので、今回は2011年に運用終了のALOSに搭載されたAVNIR-2の画像を使って分析します。
分析内容としては、2008年と2011年に、ちょうどお台場の開発が進んでいたので、時系列の比較ができます。
分析プロセス
データソース
衛星画像は、TellusのAPIから入手しました。
1B1データなので、オルソ補正は済みだけど、幾何補正・放射量較正などの前処理が必要です。こちらの処理を割愛します。
分析範囲
今回はお台場に絞って分析したくて、座標はかなり狭くまで指定しています。
min_lon, max_lon = 139.76, 139.80
min_lat, max_lat = 35.6, 35.65
上記の範囲を絞って、トリミングした画像はこちらになります。
tif_path = 'data/AVNIR-2_2008/processed_tif_2008.tif'
with rasterio.open(tif_path) as src:
window = Window(0, 0, 445, 490)
subset = src.read(window=window)
rgb = np.dstack([subset[2], subset[1], subset[0]])
plt.figure(figsize=(8,8))
plt.imshow(rgb)
plt.title("RGB Composite")
plt.axis('off')
plt.show()
2008年の画像
tif_path = 'data/AVNIR-2_2011/processed_tif_2011.tif'
with rasterio.open(tif_path) as src:
window = Window(0, 0, 445, 490)
subset = src.read(window=window)
rgb = np.dstack([subset[2], subset[1], subset[0]])
plt.figure(figsize=(8,8))
plt.imshow(rgb)
plt.title("RGB Composite")
plt.axis('off')
plt.show()
2011年の画像
特に右上と真ん中の部分は開発が進んだことがわかりますね。
変化ベクトル分析（Change Vector Analysis）
ここから時系列の比較をします。
まずはCVA、変化ベクトル分析（Change Vector Analysis）、を使ってみます。
ベクトルのノルムはnp.linalg.norm()で取得します。
tif_2011 = 'data/AVNIR-2_2011/processed_tif_2011.tif'
with rasterio.open(tif_2011) as src:
window = Window(0, 0, 445, 490)
subset_2011 = src.read(window=window)
tif_2008 = 'data/AVNIR-2_2008/processed_tif_2008.tif'
with rasterio.open(tif_2008) as src:
window = Window(0, 0, 445, 490)
subset_2008 = src.read(window=window)
subset_2008 = subset_2008[:, :subset_2011.shape[1], :]
change_magnitude = np.linalg.norm(subset_2011 - subset_2008, axis=0)
結果をヒートマップで表示させます。
plt.figure(figsize=(8, 6))
plt.imshow(change_magnitude, cmap='viridis', interpolation='nearest')
plt.colorbar(label='Change Magnitude')
plt.title('Change Magnitude Heatmap')
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.show()
これで結果が非常にわかりやすいですね。海での船を除いて、緑の部分はこの3年間で開発した部分となります。主に右上の部分は一番目立つとなっています。
主成分分析（PCA）
次に、別の手法も使って比較してみます。
主成分分析は、都市部での建物の高さや密度の変化が抽出できます。
こちらはsklearnのライブラリを使います。
from sklearn.decomposition import PCA
stacked_data = np.stack([subset_2008, subset_2011], axis=0)
reshaped_data = stacked_data.reshape(stacked_data.shape[0] * stacked_data.shape[1], -1).T
pca = PCA(n_components=3)
pca_result = pca.fit_transform(reshaped_data)
rows, cols = subset_2008.shape[1:]
pca_image = pca_result.reshape(rows, cols, -1)
結果をカラーバーグラフで表示させます。
for i in range(3):
plt.figure(figsize=(8, 6))
plt.imshow(pca_image[:, :, i], cmap='viridis')
plt.colorbar(label=f'Principal Component {i+1}')
plt.title(f'PCA Component {i+1}')
plt.axis('off')
plt.show()
今回は適当に3つの主成分に圧縮しましたが、結果はCVAと比べて、より真ん中の部分の開発状況も確認できましたね。
偽色合成画像（NIR-R-G False Color Composite）
最後に、比較ではないけど、都市開発にあたって、お台場での植生も確認していきたいです。
NIR-R-G（近赤外-赤-緑）の偽色合成画像で、地表面の特性を強調することができます。
rgb_urban = np.dstack([subset_2011[3], subset_2011[2], subset_2011[1]])
rgb_normalized = rgb_urban / np.max(rgb_urban)
plt.figure(figsize=(8, 8))
plt.imshow(rgb_normalized)
plt.title("NIR-Red-Green Composite")
plt.axis('off')
plt.show()
結果はこちらになります。
やや赤いのは、植生または植物の部分です。特に潮風公園の部分は一番目立つですね。
灰色の部分は都市部・建物・道路となって、お台場のほとんどもこの色になっています。
最後暗い黒色や青色は水域です。