4. 特徴量の計算

機械の気持ちになる作業

各種特徴量を並べたものを特徴ベクトルといいます。どんな特徴ベクトルを分類に使うかはつねに悩みどころです（リモセン屋の存在意義のひとつかもしれません）。特徴ベクトルの次元（使う特徴量の数）は、少なすぎれば土地被覆を分離しきれませんが、多すぎても次元の呪い [1] のせいで精度が落ちたり大量の教師が必要になったりします。というわけで、各土地被覆のスペクトル特性・空間的な特性をよく考えて、有効な特徴ベクトルを作成する必要があります。

たとえば、植生は近赤外光をよく反射するいっぽう、可視光をあまり反射しない（吸収する）という特性があります。そこで、光学衛星画像で植生とそれ以外をわけるのには、両者の反射率の差をみてやるのがいい、と考えられます。実際、「両者の反射率の差」を正規化した、正規化植生指数（NDVI）というのがよく使われます。

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

いろいろ考えた結果、以下の特徴量を計算することにしました。

Landsat8 2015/03/31の画像から、「band1〜7の反射率、NDVI、NDWI、NDVIのテクスチャ４種類（asm, contrast, entropy, correlation）」
DEMデータから、各ピクセルの「標高、傾斜、方位角」
Landsat8 2015/10/25と2015/03/31の差分（NDVIとNDWIのみ）

NDWI（正規化水指数）は水を強調し[2]、テクスチャは近傍のピクセルの不均一性を表します（同じ森林でも不均一性によって広葉樹/針葉樹/竹林などをわけられるかもしれない）。またLandsat8の2時期の差分をとることで、畑と水田の区別、常緑林と落葉林の区別などがつくかもしれません（水田には水が入る5月あたりの画像があると良かったかも、落葉林には冬の画像もあるとなお良かったかも）。

特徴量は合計18個、つまり特徴ベクトルの次元は18です。

ではLandsat8 2015/03/31の処理から始めましょう。

テクスチャ計算に必要なグレーレベル同時生起行列（GLCM）を求める際、数値を離散化をする必要があるので、途中NDVIを0～255の離散データにスケール変更しています。詳細は省きますが、上記例では「5x5のウィンドウの中で、あるピクセルと、隣接（距離1）するピクセルのグレースケールがどういう組み合わせか」を集計しマトリクスにまとめたものがGLCMになります。隣接するピクセルの方向もオプションで指定できますが、今回はデフォルトである「すべての方向（0°, 45°, 90°, 135°）についての平均値を求める」計算になっています。GLCMをもとにそれぞれ4種類の統計量を計算し、これをウィンドウの中心ピクセルの値として割り当てることで、NDVIの空間不均一性を示す4種類の特徴画像（ASM, Contrast, Entropy, Correlation）ができあがります。

つづいて傾斜と方位角を計算します。regionの縁がAW3D30シーンの分割境界にちょうどかぶっているため、まずちょっと広めのregionに再設定し、シーンモザイクしてから処理しましょう。

slopeが傾斜（°）、aspectが方位角（東が0°で時計回りに増加）です。ちなみに蛇足ですが、srtmカラースケールにしたDSMとaspect画像があると、次のようにして「それっぽい地形図」を描くことができます。

カラースケールの単位はメートルです。このようにすると筑波山の地形が見やすくなりますし、霞ヶ浦に注ぐ河川網（恋瀬川、桜川）なども見てとれます。

とやると、当該地域内の最大標高は872 mと出ます。実測の筑波山頂（女体山）が標高877 mなので、まあそういう感じです。

AW3D30では霞ヶ浦などの水域マスクデータも一緒に提供されていますが、今回はあえて使わずに、水域も含めて機械学習で分類していきましょう。

最後に、2時期の差分画像を計算します。

2時期の差分画像は、全特徴量について作ってもいいのですが、とりあえずNDVIとNDWIだけについて計算しました。あまり次元が増えすぎるのもどうかと思いますので・・・

次元の数が気になる人は適当にいくつか減らしてみたり、PCAとかやってもいいんじゃないでしょうか。

以上で必要な特徴画像が揃いました。機械学習でありがちなことですが、データ集めと前処理だけで既にそれなりに消耗している自分がいます。

教師をつくる作業

では教師データを作りましょう。大事なのは現地写真からつくった検証情報です。とりあえず中身を見てみます。

カテゴリ名、空間代表性（m）、緯度（°）、経度（°）の順に並んだデータです。例えば一行目だと、緯度経度（36.1115000372, 140.206746672）の点を中心に直径30 mくらいの範囲は竹林でした、という意味です。緯度経度の測位はGPSカメラ（単独測位）によっているので、精度は10 m程度、つまりざっくり言うと小数点以下6ケタめあたり以降の数字は意味ないのですが、とくに害もないのでこのまま読み込むことにします。ただし、カテゴリ名は整数値のほうが何かと扱いやすいので、最初に変換してしまいましょう。

じゃあなんであらかじめ変換したデータを用意しておかなかったのかと思いますよね。僕も思います

ではGRASSに読み込んでみましょう。点群の色付けのために、何でもいいのでエディタで以下のようなテキストファイル（reference_color.txt）を作ってカレントディレクトリに保存しておいてください。

1 110:54:42 #DF

2 44:82:60 #EF

3 133:122:60 #bamboo

4 209:112:71 #barren

5 255:255:0 #cropland

6 0:225:0 #grassland

7 255:0:0 #orchard

8 40:175:148 #rice

9 165:165:165 #urban

10 0:0:255 #water

予想通りですがかなり点が偏ってますね。霞ヶ浦の真ん中の点は遊覧船に乗ったのでしょうか。legendを描くときの注意ですが、GRASSはvectorデータに対する色凡例表示コマンドがないので（2018.12現在）、rasterデータに変換後凡例を表示しています。では、機械学習用のcsvファイルを作りましょう。これがトレーニング/検証データになります。けっきょく全てはこのデータ次第です。

r.whatはpointsで与えた点（ベクタ）におけるmapのラスタ値を抽出してくれます。出力の区切りはseparatorで指定します。CAT（土地被覆カテゴリ）のデータは「ラスタ化したreference_LULC」から引っ張ってくることでコードがシンプルになっています。また、欠測データは"*"で出力され、これがあるとpythonの機械学習時にエラーを起こすので、事前に除去しています。

さて、次回はいよいよこのデータを使って機械学習をやっていきます。grassのターミナル上で、pythonを起動しておいてください。

次章: Pythonによる機械学習

[1] Bishop, 2006. Pattern recognition and machine learning. 2006, pp. 33-38.

[2] 正確には改良正規化水指数（MNDWI）です（Xu, 2006）が、表記の煩雑さを避けるためここではNDWIで統一しています。

Google Sites

Report abuse