Landsatによるアマゾン熱帯林の変化の推定
2016/10/18 筑波大学生命環境系 奈佐原(西田)顕郎
Landsatによるアマゾン熱帯林の変化の推定
焼畑や伐採などによる熱帯雨林の破壊は、大きな環境問題のひとつであり、これを衛星リモートセンシングで同定することが、盛んに行われている。我々も、それをためしてみよう。
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年
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年
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年
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
黄色が0 (植生のまま) 水色が1 (裸地のまま) 紫色が2 (植生が裸地に) 赤色が3 (裸地が植生に)
一部をズームしたところ
これを見ると、新しい道路などで裸地が拡大する一方で、裸地に植生が戻っているようにも見える部分がある。では、それぞれの面積を出してみよう:
課題: r.statsコマンドを用いて、植生面積の変化を検討せよ。
注意: 実際は、植生変動の解析はこのように単純ではなく、たとえばNDVIのしきい値を0.4とするのがどれだけ妥当か、ということや、各種のノイズ除去などの処理が必要だし、最終的な結果は、現地の何ヶ所かで、実際に検証しなければならない。
Keyword(s):
References:[GIS入門]