3. データ準備と読み込み

データ準備

ではさっそく、茨城県つくば市界隈(140.0E-140.5E, 36.0N-36.5N)の土地被覆分類を公開衛星データを使ってやってみましょう。USGSのLandsat8地表面反射率プロダクト、およびJAXAの数値地表モデルであるAW3D30 [1]を入力データとして使ってみます。

ここから以下のものをDLし、tarボールは展開してください。

landsat8_20150331.tar.gz #2015/03/31のband1〜7地表面反射率

landsat8_20151025.tar.gz #2015/10/25のband1〜7地表面反射率

aw3d30.tar.gz #N35-36, E139-140のDSMデータ

tsukuba_photos_reference.csv #教師&検証データ

展開

Landsat8データは事前にGoogle Earth EngineからDLしたものですが、USGSの本家にいけば他にもいろいろな入手方法が指示されていると思います。Google Earth Engineの使い方を知りたい方はこちらのコンテンツを参考にしてください。AW3D30はJAXAのサイトにユーザ登録後DLできます。

Landsat8, AW3D30ともに投影法は正距円筒図法(lat-lon)です。教師&検証データは、地上検証データセットであるSACLAJをもとに作成したものです。対象範囲内で撮影された空間代表性30 m以上の地上写真(合計3839点)を、今回の土地被覆分類用に次の10カテゴリに整理しなおしました:

DF(落葉林), EF(常緑林), bamboo, barren, croplands, grasslands, orchard, rice, urban, water


さて、まだ解析が始まってもいないですが、この段階で早くも疑問が湧いてきます。。。

  1. 土地被覆のカテゴリはこの10カテゴリでいいの?

  2. 写真撮影地点には空間的に偏りがあると思うけどいいの?

  3. カテゴリ間でデータサイズが違ってるけどいいの?

しかし今回は練習ですのでスルーしましょう。いちおう少しだけ考察らしきものを挙げておきますが早く実践をやりたい人はデータ読み込みまでスキップしてください。

  1. カテゴリ数について、一般には多くなるほど分類精度が低くなる傾向にありますが、少なすぎても特徴がうまく分離できず良くありません。地域ごとに特徴的な土地被覆カテゴリを追加する必要もあるでしょう。有名なのはIGBPの分類定義で、常緑針葉樹・常緑広葉樹・落葉針葉樹・落葉広葉樹・混交林・密な灌木地・疎な灌木地・木本サバンナ・草本サバンナ・草原・湿地・耕作地・都市・耕作地と自然植生の混ざった土地・雪氷域・裸地・水域の17カテゴリにわけます。今回はこのカテゴリ分けを参考にしつつも、よしなに分けてみました。畑と水田を分けるのは重要な気がします。

  2. よくありません。理想をいえば空間的にランダムにとるべきです。現地調査の場合、地形とか交通とか地権の都合で、完全に空間的にランダムに調査するのは困難です。こういう場合、現地写真を一次情報として参考にしつつ、Google Earthなどの高解像度衛星画像を判読して、ランダムに教師点を増やす、ということをやります。ただその場合にも、空間的なランダム性を重視するあまり、ミクセル(ひとつの画素内に複数の土地被覆が混在したもの)などの質の悪い教師をとってしまえば、それはそれで精度の低下につながりうるので、まあこれもよしなに、としか言えません。

  3. 機械学習的には、カテゴリ間の偏りを無くすようにデータオーギュメンテーション(足りないカテゴリの水増し)とかアンダーサンプリング(多すぎるカテゴリの間引き)をやれ、みたいなことが言われますが、実感から言うとそれで精度が上がるとは限らない気がします。データを間引く場合、そもそも貴重な教師のサイズが減ってしまいます(かさ増しはどうやるのがいいのかイマイチわからん)。次善の策として、既存の土地被覆図などを利用した層別ランダムサンプリングで各カテゴリ同サイズのサンプルを確保、という道もあります。しかし、カテゴリの偏りはそもそもの土地被覆の現れやすさを反映しているわけで、この特徴を破壊することの問題もありそうです―カテゴリの偏りはベイズの識別規則における事前確率みたいな効果をもつような雰囲気があります。というわけでこれもよしなに。けっきょく過程(仮定?)はどうあれ、汎用性と精度こそ正義という筋肉な世界です。


ひとしきり瞑想が終わったところで、データを読み込むとしましょう。まずはGRASSを起動し、前回作っておいたEPSG:4326 (wgs84, latlon) の環境に入ります。

続いてデータを読み込みます。以下、> マークはGRASSコンソールであることを表します。

--oマークは「上書き」のオプションです。GRASS7.0以降では、これを明示的につけないと既存マップがある場合に上書きをしてくれなくなったので(警告が出てスキップになる)、癖としてつけるようにしています。今回は初回読み込みなので、なくても問題ありません。

regionの設定をしましょう。分解能は1arcsec (約30 m) にします。

では、とりあえず春のtrue color画像をみてみましょう。

筑波山と霞ヶ浦あたりを範囲にいれておきました。ぐりぐりズームしたりして遊びましょう。

d.monでコンソールから出力モニタを立ち上げてますが、GUI操作が好きな人はデフォルトで立ち上がっているGUIから表示してもいいと思います。2015/10/25のほうが見たければ、DATE=20151025に変えてください。

地表面反射率画像は、デフォルトの0~1.0のグレースケールで表示するとかなり暗くなってしまうため、0~0.2の範囲でストレッチをかけたカラースケールを作成しています。ちなみに、こういった手動でのカラースケール作成が面倒な場合は、次のコマンドも使えます:

指定したバンドごとに自動でストレッチをかけたカラースケールを割り当ててくれます。デフォルトでは2%quantile~98%quantileが[0,1]に対応するよう線形にストレッチがかかります。

次回は読み込んだデータからピクセルごとに特徴量を計算してみましょう。GRASS大活躍です。

次章: 特徴量の計算