Landsat画像を用いたパンシャープン画像の作成

2011/10/26 筑波大学 田中 健太郎

衛星画像の解像度は、衛星センサの性能によって決っている。例えばLandsatに搭載されているETM+センサは、可視域〜中間赤外波長域までの観測域を持ち、それぞれ14.25m~60mの空間解像度を持っている。

ETMセンサ主要諸元

ここで空間解像度がもっとも良いのは、パンクロマチックバンド (バンド8) の15mである。しかしながら得られる画像は、モノクロである。 一方で可視波長域 (赤 : バンド3, 青 : バンド1, 緑 : バンド2) における空間解像度は28.5mであるため、トゥルーカラー画像もしくはフォルスカラー画像を作成した際の空間分解能は28.5mとなっている。

素晴らしいことに、パンクロマチックバンド (バンド8) と 可視波長域 (赤 : バンド3, 青 : バンド1, 緑 : バンド2) の画像を組み合わせることによって、高分解能化したカラー画像を得ることが出来る。こうして得られた合成画像のことをパンシャープン画像とよぶ。パンシャープン画像の作成アルゴリズムについては、様々な研究が行われているが、本章では RGB<->HSI変換によるパンシャープ処理を行う。

RGB画像は、Red画像、GREEN画像、Blue画像の3枚の画像を合成したカラー画像であり、馴染み深い画像であるかと思う。

HSI画像は、カラー画像をRGBの3成分によって表現するのではなく、Hue (色相)、Saturation (彩度), Intensity (明度)の3成分で表現した画像のことである。

RGB画像をHSI画像に変換する処理のことを、HSI変換と呼ぶ。

Landsat画像の取得

高分解能衛星Landsat: 植生指標NDVI, 温度画像の章と同様の場所から衛星画像をダウンロードする。

関東地方周辺の、2001年9月23日のLandsat-ETM+画像であり、ダウンロードするのはバンド1,2,3,4,8 の画像とする。

※バンド1,2,3,4のデータを既に入手している場合は下記1行目は必要無い。

wget -c ftp://ftp.glcf.umiacs.umd.edu/glcf/Landsat/WRS2/p107/r035/p107r035_7x20010924.ETM-EarthSat-Orthorectified/*nn[1-4]0.tif.gz
wget -c ftp://ftp.glcf.umiacs.umd.edu/glcf/Landsat/WRS2/p107/r035/p107r035_7x20010924.ETM-EarthSat-Orthorectified/p107r035_7p20010924_z54_nn80.tif.gz
 

もしくは、ここからダウンロードする。

解凍方法

for i in *tif.gz; do gunzip $i; done

GRASSへの読み込み

MAPSETは、高分解能衛星Landsat: 植生指標NDVI, 温度画像と同様のものを使用する。

grassへ衛星画像を取り込む場合、r.in.gdalコマンドを実行する。

r.in.gdal input=p107r035_7t20010924_z54_nn10.tif output=B1
r.in.gdal input=p107r035_7t20010924_z54_nn20.tif output=B2
r.in.gdal input=p107r035_7t20010924_z54_nn30.tif output=B3
r.in.gdal input=p107r035_7t20010924_z54_nn40.tif output=B4
r.in.gdal input=p107r035_7p20010924_z54_nn80.tif output=B8

これでGRASSへ衛星画像を取り組む作業は完了である。

GRASSに読み込んだ画像の空間解像とを確認してみよう。

例) B1
g.region rast=B1 -g
n=4103757.75
s=3876042.75
w=290714.25
e=546957.75
nsres=28.5
ewres=28.5
rows=7990
cols=8991
cells=71838090

例) B8
g.region rast=B8 -g
n=4103757.75
s=3876042.75
w=290714.25
e=546957.75
nsres=14.25
ewres=14.25
rows=15980
cols=17982
cells=287352360

空間解像度を示す値 (nsres, ewres) を見てみると、B1~B4の空間解像度が28.5m、B8の空間解像度が14.25mであることがわかる。

データの準備が整ったので、早速パンシャープン処理を行う。

パンシャープン処理には、GRASS内で実装されているi.rgb.his コマンドを使用する。

まずB1,B2,B3の画像をそれぞれ青, 緑, 赤に割り当て、his変換を行う。

$ g.region rast=B1
$ i.rgb.his red=B3 g=B2 b=B1 h=hue i=intensity s=saturation

※GRASS7.0.3では"i.rgb.his red=B3 green=B2 blue=B1 hue=hue i=intensity s=saturation"（追記：2018/04/10 菊島）

作成された画像(hue, intensity, saturation)を表示してみよう。

$ d.mon x0
$ d.rast hue
$ d.rast intensity
$ d.rast saturation

パンシャープン画像は、作成されたintensity画像とB8画像を入れ替え、hue画像とsaturation画像およびB8画像から、RGB画像を合成することによって作成される。

HSI画像からRGB画像への変換には、i.his.rgbコマンドを使用する。※最初に現在のRegionの空間解像度を空間解像度14.25に変更する。

g.region res=14.25
i.his.rgb  h=hue i=B8 s=saturation red=B3_pan g=B2_pan b=B1_pan

以上で、パンシャープン処理は完了である。

作成された画像を早速表示してみよう。

for i in 1 2 3 ; do r.colors B${i}_pan color=grey ; done
d.rgb  r=B3_pan g=B2_pan b=B1_pan

対象範囲を狭くし、元のRGB画像と比較してみよう。

g.region n=3972087.75 s=3966858 w=424578.75 e=430649.25
d.erase
d.rgb  r=B3_pan g=B2_pan b=B1_pan
d.rgb r=B3 g=B2 b=B1

左図 パンシャープン処理後のトゥルーカラー画像, 右図 パンシャープン処理前のトゥルーカラー画像

パンシャープン処理前の画像では水田の輪郭がぼんやりとしているのに対し、パンシャープン処理後では輪郭がはっきりしている事がわかる。

本章では、パンシャープン処理の方法としてHSI変換法を用いた。本章の最初にも述べたが、パンシャープン処理のアルゴリズムは他にもあり、様々な研究が行われている。 興味のある人は、色々と調べてみるとよい。最後にパンシャープン画像の色は擬似的なものであり、スペクトル特性が元の画像から変化しているという欠点を持つ。したがって定量的な解析を行う際は、十分注意する必要がある。