とらりもんHOME  Index  Search  Changes  Login

Landsatによるアマゾン熱帯林の変化の推定

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" --o
r.mapcalc "1986B2L=1986B2*1.1750981-2.8399999" --o
r.mapcalc "1986B3L=1986B3*0.8057647-1.1700000" --o
r.mapcalc "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" --o
r.mapcalc "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" --o
r.mapcalc "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 
1986_Amazon_Landsat_RGB.png
1986年
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
1986_Amazon_Landsat_RGB_zoom.png
1986年
2001_Amazon_Landsat_RGB_zoom.png
2001年

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

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

r.mapcalc "1986NDVI=(1986B4L-1986B3L)/(1986B4L+1986B3L)" --o
d.rast 1986NDVI
r.mapcalc "2001NDVI=(2001B4L-2001B3L)/(2001B4L+2001B3L)" --o
d.rast 2001NDVI
1986_Amazon_NDVI.png
1986年
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" --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" --o
d.rast change
Amazon_NDVI_change.png
黄色が0 (植生のまま) 
水色が1 (裸地のまま)
紫色が2 (植生が裸地に)
赤色が3 (裸地が植生に)
Amazon_NDVI_change_zoom.png
一部をズームしたところ

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

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

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

Last modified:2016/10/18 10:46:19
Keyword(s):
References:[GIS入門]