とらりもん - Landsatによるアマゾン熱帯林の変化の推定 Diff

2016/10/18 筑波大学生命環境系 奈佐原(西田)顕郎
!Landsatによるアマゾン熱帯林の変化の推定

　焼畑や伐採などによる熱帯雨林の破壊は、大きな環境問題のひとつであり、これを衛星リモートセンシングで同定することが、盛んに行われている。我々も、それをためしてみよう。
!作業ディレクトリを作る:

mkdir ~/work_Amazon
cd ~/work_Amazon

!Landsatデータの取得
　例によって、メリーランド大学の無料データベースを利用し、以下のデータを取得し、解凍する：

wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p230/r062/p230r62_5t19860803.TM-EarthSat-Orthorectified/*nn[1-4]*tif.gz
wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p230/r062/p230r62_5t19860803.TM-EarthSat-Orthorectified/*.met
wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p230/r062/p230r062_7x20010804.ETM-EarthSat-Orthorectified/*n[1-4]?.tif.gz
wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p230/r062/p230r062_7x20010804.ETM-EarthSat-Orthorectified/*.met
for i in *.gz; do gunzip $i; done

!GRASSの起動と設定
取得したデータのヘッダーデータp230r062_7x20010804.metなどの中をみると、

cat p230r062_7x20010804.met

(以下, 表示)
GROUP = PROJECTION_PARAMETERS  
REFERENCE_DATUM = "WGS84"
REFERENCE_ELLIPSOID = "WGS84"
GRID_CELL_ORIGIN = "Center"
UL_GRID_LINE_NUMBER = 1
....
RESAMPLING_OPTION = "NN"
MAP_PROJECTION = "UTM"
END_GROUP = PROJECTION_PARAMETERS
GROUP = UTM_PARAMETERS
ZONE_NUMBER = +21

という記述がある。従って、このデータは、測地系(datum)はWGS84, 投影法はUTM、そのzoneは21であることがわかる。そこで、新たに、これに対応するようなLOCATIONを、UTM21という名前で立ち上げよう。現在、ほかのLOCATIONに入ってGRASSを使っている人は、いったんexitコマンドでGRASSから抜け出して、改めてGRASSを起動してUTM21というLOCATIONを作り直す。

grass -text
cd ~/work_Amazon

!データの読み込み
では、データをGRASSに読み込もう：

r.in.gdal input=p230r62_5t19860803_nn1.tif output=1986B1 --o
r.in.gdal input=p230r62_5t19860803_nn2.tif output=1986B2 --o
r.in.gdal input=p230r62_5t19860803_nn3.tif output=1986B3 --o
r.in.gdal input=p230r62_5t19860803_nn4.tif output=1986B4 --o
r.in.gdal input=p230r062_7t20010804_z21_nn10.tif output=2001B1 --o
r.in.gdal input=p230r062_7t20010804_z21_nn20.tif output=2001B2 --o
r.in.gdal input=p230r062_7t20010804_z21_nn30.tif output=2001B3 --o
r.in.gdal input=p230r062_7t20010804_z21_nn40.tif output=2001B4 --o

ヘッダーデータを参照して、これらを輝度変換する:

g.region rast=1986B1
r.mapcalc "1986B1L=1986B1*0.6024314-1.5200000""1986B1L=1986B1*0.6024314-1.5200000" --o
r.mapcalc "1986B2L=1986B2*1.1750981-2.8399999""1986B2L=1986B2*1.1750981-2.8399999" --o
r.mapcalc "1986B3L=1986B3*0.8057647-1.1700000""1986B3L=1986B3*0.8057647-1.1700000" --o
r.mapcalc "1986B4L=1986B4*0.8145490-1.5100000""1986B4L=1986B4*0.8145490-1.5100000" --o
g.region rast=2001B1
r.mapcalc "2001B1L=(2001B1-1)*(191.6+6.2)/(255.0-1.0)-6.2""2001B1L=(2001B1-1)*(191.6+6.2)/(255.0-1.0)-6.2" --o
r.mapcalc "2001B2L=(2001B2-1)*(196.5+6.4)/(255.0-1.0)-6.4""2001B2L=(2001B2-1)*(196.5+6.4)/(255.0-1.0)-6.4" --o
r.mapcalc "2001B3L=(2001B3-1)*(152.9+5.0)/(255.0-1.0)-5.0""2001B3L=(2001B3-1)*(152.9+5.0)/(255.0-1.0)-5.0" --o
r.mapcalc "2001B4L=(2001B4-1)*(241.1+5.1)/(255.0-1.0)-5.1""2001B4L=(2001B4-1)*(241.1+5.1)/(255.0-1.0)-5.1" --o

白黒のカラーテーブルをテキストファイルとして作成:

cat << EOF >color_BW.txt
-1000 black
0 black
160 white
255 white
1000 white
nv black
EOF

このカラーテーブルを各mapに適用する。
(GRASS6.xの場合)
for i in 1 2 3 4; do
r.colors map=1986B${i}L color=rules < color_BW.txt
r.colors map=2001B${i}L color=rules < color_BW.txt
done

(GRASS7.xの場合)
for i in 1 2 3 4; do
r.colors map=1986B${i}L rules=color_BW.txt
r.colors map=2001B${i}L rules=color_BW.txt
done

表示画面を立ち上げる:
(GRASS6.xの場合)
d.mon x0
(GRASS7.xの場合)
d.mon wx0

表示してみる:
d.rgb blue=1986B1L green=1986B2L red=1986B3L
d.rgb blue=2001B1L green=2001B2L red=2001B3L

{{attach_view(1986_Amazon_Landsat_RGB.png)}}

!!!!!1986年

{{attach_view(2001_Amazon_Landsat_RGB.png)}}

!!!!!2001年

雲のかかっていない、中心から左寄りの領域をズームしよう：
d.erase
d.rast 1986B1L

(GRASS6.xの場合)
d.zoom
--> マウスで好きな範囲を選ぶ

(GRASS7.xの場合)
表示画面の上部のメニューアイコンから虫眼鏡(+)を選んで,
--> マウスで好きな範囲を選ぶ

ズームができたら再表示:
d.rgb blue=1986B1L green=1986B2L red=1986B3L
d.rgb blue=2001B1L green=2001B2L red=2001B3L

{{attach_view(1986_Amazon_Landsat_RGB_zoom.png)}}

!!!!!1986年

{{attach_view(2001_Amazon_Landsat_RGB_zoom.png)}}

!!!!!2001年

この状態でも、新たに道ができたりしている様子がよくわかるが、もうちょっと定量的に変化を見てみよう。

まず、各時期のNDVI画像を求めて表示してみる：

r.mapcalc "1986NDVI=(1986B4L-1986B3L)/(1986B4L+1986B3L)""1986NDVI=(1986B4L-1986B3L)/(1986B4L+1986B3L)" --o
d.rast 1986NDVI
r.mapcalc "2001NDVI=(2001B4L-2001B3L)/(2001B4L+2001B3L)""2001NDVI=(2001B4L-2001B3L)/(2001B4L+2001B3L)" --o
d.rast 2001NDVI

{{attach_view(1986_Amazon_NDVI.png)}}

!!!!!1986年

{{attach_view(2001_Amazon_NDVI.png)}}

!!!!!2001年

d.what.rastコマンドでそれぞれのNDVI画像の値を調べると、だいたい0.4くらいを境にして、裸地が同定できるようである。そこで、

両年とも植生なら0
両年とも裸地なら1
1986年に植生で2001年に裸地なら2
1986年に裸地で2001年に植生なら3

となるようなマップを作ってみよう：

(GRASS6.xの場合)
r.mapcalc "change=if(1986NDVI<0.4 && 2001NDVI<0.4)+if(0.4<=1986NDVI && 2001NDVI<0.4)*2+if(1986NDVI<0.4 && 0.4<=2001NDVI)*3"0.4<=2001NDVI)*3" --o

(GRASS7.xの場合)
r.mapcalc expression="change=if(1986NDVI<0.4 && 2001NDVI<0.4)+if(0.4<=1986NDVI && 2001NDVI<0.4)*2+if(1986NDVI<0.4 && 0.4<=2001NDVI)*3"0.4<=2001NDVI)*3" --o

d.rast change

{{attach_view(Amazon_NDVI_change.png)}}

黄色が0 (植生のまま)
水色が1 (裸地のまま)
紫色が2 (植生が裸地に)
赤色が3 (裸地が植生に)

{{attach_view(Amazon_NDVI_change_zoom.png)}}

!!!!!一部をズームしたところ

　これを見ると、新しい道路などで裸地が拡大する一方で、裸地に植生が戻っているようにも見える部分がある。では、それぞれの面積を出してみよう:

!!!!!課題： r.statsコマンドを用いて、植生面積の変化を検討せよ。

注意： 実際は、植生変動の解析はこのように単純ではなく、たとえばNDVIのしきい値を0.4とするのがどれだけ妥当か、ということや、各種のノイズ除去などの処理が必要だし、最終的な結果は、現地の何ヶ所かで、実際に検証しなければならない。