中分解能衛星Terra/MODIS

筑波大学農林工学系 西田顕郎

衛星データの処理の基本を学ぶ。題材として、米国・NASAのMODISという衛星センサーのデータを用いる。MODISは、空間分解能250m〜1000m、時間分解能1日の、広域環境観測センサーであり、無料でデータを入手できる。入手方法・処理方法についてはこちらを参考にしてください。

1. grassを起動し、

location: latlon
mapset: PERMANENT

に入る。上のlocationやmapsetが無い場合は、こちらを参考にして作ること。

2. 以下のコマンドで、MODISデータ(この実習用に処理済)をダウンロード・解凍・展開する：

wget -c http://pen.agbi.tsukuba.ac.jp/~torarimon/data4download/MODIS_Tokyo.tar.gz
tar zxvf MODIS_Tokyo.tar.gz

3. データの概要を参照する：

cat README

個々のデータにはMODIS_Tokyo_2000_121_b01のようなファイル名がついているが、2000は年号、121はJulian day (1月1日を起点として通算日数で表現した日付。121は4月30日)、b01はバンド名(36本あるMODISの観測バンドのうちどれか)を表している。

このデータは、MODISが観測した可視・近赤外データを関東平野を切り出したもので、2000年春から年末まで、8日ごとのコンポジットである。コンポジットとは、雲のかかっていない部分をつぎはぎして合成したデータのことである。

READMEを参考にして、どれかひとつの時期のファイルを4つのバンド全てについて読み込む。これらのデータは、raw binary、すなわち、ピクセル値が単純に並んだバイナリファイルであり、これまでのGTOPO30やSRTMと同様なのだが、それらの有名なデータとは違い、自動的に読み込んでくれる便利なコマンドは無い。そこで、raw binaryデータを、全てのパラメータを明示的に指定して読み込むためのコマンド：r.in.binを使う：

r.in.bin input=MODIS_Tokyo_2000_121_b01 output=MODIS_Tokyo_2000_121_b01 west=138.916667 east=140.916667 north=36.75 south=35.0 r=420 c=480 byte=2 anull=-9999 -b -s
r.in.bin input=MODIS_Tokyo_2000_121_b02 output=MODIS_Tokyo_2000_121_b02 west=138.916667 east=140.916667 north=36.75 south=35.0 r=420 c=480 byte=2 anull=-9999 -b -s
r.in.bin input=MODIS_Tokyo_2000_121_b03 output=MODIS_Tokyo_2000_121_b03 west=138.916667 east=140.916667 north=36.75 south=35.0 r=420 c=480 byte=2 anull=-9999 -b -s
r.in.bin input=MODIS_Tokyo_2000_121_b04 output=MODIS_Tokyo_2000_121_b04 west=138.916667 east=140.916667 north=36.75 south=35.0 r=420 c=480 byte=2 anull=-9999 -b -s

ここで、r=420 c=480というのは、縦のピクセル数(row)と横のピクセル数(column)、 byte=2というのは、各ピクセルが２バイトの整数値であるという指定、anull=-9999というのは、null値が便宜上-9999という値になっていること、-bというのはエンディアンの変換指定である。（エンディアンについて忘れてしまった人は、こちらを参照）

さて、こうやって読み込んだデータを試しに表示してみよう：

d.mon x0
g.region rast=MODIS_Tokyo_2000_121_b01
d.rast MODIS_Tokyo_2000_121_b01

衛星センサーは人間の目と同様に、様々な波長の光を見分けているので、それらを適当に光の3原色(青、緑、赤)にわりあててディスプレーに表示すれば、目で見たかんじの画像になる。MODISはちょうど、バンド1が赤、バンド3が青、バンド4が緑に対応するので、そのように色を割り当てて合成表示してみよう。

まず、上で読んだラスターデータのカラーテーブルはデフォルトでは虹色のグラデーションになっているが、カラー合成をするには個々のラスターデータは白黒でなくては困るので、カラーテーブルを変換しよう：

r.colors map=MODIS_Tokyo_2000_121_b01 color=grey
r.colors map=MODIS_Tokyo_2000_121_b02 color=grey
r.colors map=MODIS_Tokyo_2000_121_b03 color=grey
r.colors map=MODIS_Tokyo_2000_121_b04 color=grey

で、カラー合成は以下のようにすればよい：

d.rgb red=MODIS_Tokyo_2000_121_b01 blue=MODIS_Tokyo_2000_121_b03 green=MODIS_Tokyo_2000_121_b04

(人の目と同様の色割り当て。トルーカラー画像と呼ばれる。)

d.rgb red=MODIS_Tokyo_2000_121_b02 blue=MODIS_Tokyo_2000_121_b04 green=MODIS_Tokyo_2000_121_b01

(近赤外を赤に割り当てるやりかた。植生が赤く強調される。フォールスカラー画像と呼ばれる。)

トルーカラー画像 フォールスカラー画像

可視・近赤外画像からは、植生の有無や活性を観察することができる。その最も簡単で一般的な方法は、NDVIという、次式で定義される指標を計算することである：

	NDVI = (NIR - Red) / (NIR + Red)

ここでNIRは近赤外反射率(MODISならバンド2など)、Redは赤の反射率(MODISならバンド1)である。NDVIに関する詳細は、「リモートセンシング論」などの講義や関連図書を参考にされたい。

それでは、MODISデータを用いてGRASS上でNDVIを計算してみよう：

r.mapcalc "MODIS_Tokyo_2000_121_NDVI = 100*(MODIS_Tokyo_2000_121_b02 - MODIS_Tokyo_2000_121_b01) / (MODIS_Tokyo_2000_121_b02 + MODIS_Tokyo_2000_121_b01)"

本来はNDVIは-1から1までの値をとるが、GRASSがカラー表示するには、カラーテーブルに対応させるために、100倍して-100〜100の範囲に拡大したほうが便利である。

実際はゼロ以下のNDVI値は植生ではなくて水や雪、雲などの特殊な条件なので、この場合、NDVIは0〜100と考えて良い。

そこで、以下のようなカラーテーブルを当てておこう：

vi color_NDVI.txt

-30000 black
0 black
100 white
101 black
30000 black
nv black 


r.colors map=MODIS_Tokyo_2000_121_NDVI color=rules < color_NDVI.txt
d.rast MODIS_Tokyo_2000_121_NDVI

上のように、対象とする画像が4枚程度の1セットなら、いちいちコマンドをキーボードからうちこむのも悪くは無いが、今回は、ほぼ1年間にわたって同一箇所の44時期のデータがあり、各時期には4つのバンド画像があるので、計156枚の画像を扱わねばならない。もともと生態環境を目的にした解析では、このように時系列データが膨大になることは珍しくないが、これらのファイル全てをいちいちキーボードでコマンドを打って読んで解析するのは大変である。たとえば、これら44時期のNDVIの時系列変化を解析する、というだけで大変な作業になってしまう。

そのような場合、UNIXのシェルスクリプト機能が非常に有用である。シェルスクリプトは、同じようなコマンドの繰り返しを自動で行ってくれる。シェルスクリプトを自由に書けるようになると、作業効率は劇的に高くなる。本サイトのシェルスクリプトのページなどを参考にして学習されたい。

今回ダウンロードした156枚の画像を読んでカラーテーブルを変更するシェルスクリプト(ファイル名はMODIS_Tokyo_2000_import.shとする)は下記のようになる：

vi MODIS_Tokyo_2000_import.sh

    #!/bin/sh for i in MODIS_Tokyo_2000_???_b0? do r.in.bin input=$i output=$i west=138.916667 east=140.916667 north=36.75 south=35.0 r=420 c=480 byte=2 anull=-9999 -b r.colors map=$i color=grey done 

解説： #!/bin/sh ... スクリプトの実行シェル環境の指定。shはbourneシェル。シェルによって文法が違うので、現在使用中のシェルと区別するためにこのような指定が必要。

for i in ... doとdoneで囲まれたコマンドを繰り返し実行する。その際、シェル変数iに次々と違う値を入れる。

MODIS_Tokyo_2000_???_b0? ... カレントディレクトリにあるファイルのファイル名を検索し、MODIS_Tokyo_2000_で始まって次に任意の3文字、次に_b0、次に任意の1文字が来るようなファイル名を抽出。これを満たすファイル名が、for i inのシェル変数$iに渡される。

シェルスクリプトができたら、まず実行権限を与え(chmod)、次にそれを実行する：

chmod +x MODIS_Tokyo_2000_import.sh
./MODIS_Tokyo_2000_import.sh

すると、あとはシェルが自動的に全てのファイルをGRASSに読み込んでカラーテーブルの変更をしてくれる。

その結果を確認しておこう：

g.list rast

このg.listコマンドは、現在使用中のLOCATIONにどのようなデータが登録されているかを表示するコマンドである。156個のラスターデータが登録されているはずである。

上のスクリプトは、カレントディレクトリを検索するシェル機能を利用したのでとてもシンプルに書けた。しかし、GRASSではラスターデータは必ずしもカレントディレクトリに格納されているわけではないので、スクリプトでそれらを扱うには、いったんラスターデータのリストを取得し、編集しなければならない。

そうした処理の例として、こんどはMODIS_Tokyo_2000_NDVI.shという名のシェルスクリプトを作って、NDVIの計算を全ての時期について行ってみよう。

vi MODIS_Tokyo_2000_NDVI.sh

    #!/bin/sh for i in MODIS_Tokyo*b01 do r.mapcalc "${i%b01}NDVI = 100*(${i%b01}b02 - ${i}) / (${i%b01}b02 + ${i})" r.colors map=${i%b01}NDVI color=rules < color_NDVI.txt done 

chmod +x MODIS_Tokyo_2000_NDVI.sh
./MODIS_Tokyo_2000_NDVI.sh

これで例えば、春・夏のNDVIの違いなどを見ることができる：

d.rast MODIS_Tokyo_2000_057_NDVI
d.rast MODIS_Tokyo_2000_193_NDVI

上の2つの時期の違いは、カラー合成として観察することもできる：

d.rgb blue=MODIS_Tokyo_2000_057_NDVI red=MODIS_Tokyo_2000_193_NDVI green=MODIS_Tokyo_2000_057_NDVI

この画像で、赤いところは夏しか植生がないところ、白いところは年中ずっと植生があるところ、黒いところは年中ずっと植生が無いところである(衛星データから推測する限り)。最後に、特定の場所のピクセル値を時系列で抜きだしてグラフにしてみよう：

vi MODIS_Tsukuba_2000_NDVI.sh

   #!/bin/sh for i in MODIS_Tokyo_*b01 do echo "140.0303 36.2924 ${i%_b01}" | r.what input=${i%b01}NDVI done 

./MODIS_Tsukuba_2000_NDVI.sh | awk 'BEGIN{FS="|"}{print $3,$4}' | sed 's/MODIS_Tokyo_2000_//' | sed 's/_//' > MODIS_Tsukuba_2000_NDVI.txt
gnuplot
gnuplot> plot "MODIS_Tsukuba_2000_NDVI.txt" using 1:2 w l

つくばにおけるNDVIの季節変化