海老ノート

Google Earth Engine 苦闘の記録

石川チャレンジ5:まとめとおまけ

そんな簡単にはいかない

簡単そう、お手軽そうと思って始めましたが、予想より3倍くらいは大変でした。なんで大変だったんかなと思いだしてみます。

分析に時間がかかる

サクサク進まないと感じた理由が時間がかかること。GGEで8000点を使って石川県サイズのランダムフォレストをして、書き出すとおよそ30分かかる。で、結果を都度見ながらやりたいなとMap.addLayerを使うとタイムアウトでできない。だからお手本論文ではRでいろいろ実験する流れになってたんだなと思った。

Rでもランダムフォレストにそれなりに時間はかかるけど、随分と時間の節約にはなる。必要なサンプル数の計算とかループで回してしまえば待っているだけでできるし。

正解はどれだ?分析に適したデータとは。

現地データに基づいたこれだけの立派な植生図が全国的に整備されていることはすばらしい(だからリモセンは不要ともいえる)です。その上勝手にデータを使ってる身分であれなんですが、10m単位の分類をするために適したデータとは言えないのかもしれない。

そもそも環境省植生図はそこまで小面積を区分する仕様になっていないのでケチをつけているわけではなくって、今回はその使用がカバーしていない作業での利用なので限界はありそうだなという認識です。

土地利用区分の偏り

湿地の分布面積少なすぎ。どう頑張っても教師データを十分に確保できない。そうなると、教師データの点と点(さらにはテストデータとの点とも)が近すぎたり、データの質も問題が多くなりそう。テストデータにしたって湿地は少なくって困ります。

今回は少ないデータを勝手に増やしたけど、適当にやっただけ。少ないデータを増やすと分布(?)は描けるようになるんだろうから、特徴量のデータの分布がすごい離れている場合はその方がいいはず。だけど、重なりが大きい場合は閾値が実際の世界で一番よい点よりも少ない側データに寄る気がする。自信ないけど。

もう一点、不均衡データでは正解率が指標として頼りない。今回の場合分布面積の多い高木林やら開放水域にしておけば正解率は上がるんだけど、そんな地図があっても役に立たない。ので、もっと他の指標(F値とか?)も検討しないといけない。

ということで、不均衡データどうするか問題は不勉強だったことがわかりました。RではできるけれどもGEEでは特に記述がなさそう、と思ったら対処したという論文(SIにスクリプトがあった)があった。今回の宿題です。

www.mdpi.com

Sentinel-1あまり役に立たなかった説

期待していたSentinel-1ですが、正解率の変化とランダムフォレストの出力見る限り、劇的に役に立ったわけではなさそう。

なんでかな?と考えるに、Sentinel-1は森林と草地、あるいは水域を区分するのに使えるとされているけど、石川県の場合は標高とか傾斜でも随分と区分できちゃったってことなのかなと推測した。平地は市街地・田畑。傾斜があれば森林みたいな。お手本論文のようにヨーロッパの平地が多いところであれば、役立つ場面が多いのかもしれない。まだ使いどころがわからないので、Sentinel-1どんな感じなのか見るっていうのも宿題。

もっと良くするには?

教師データの質を上げる工夫

お手本論文でも教師データの質を上げる工夫がいくつかされています。今回のデータからできそうなことと言ったら、たとえばポリゴンの端の方のデータは除外するとか、事前にNDVI/NDWIなどの指標で外れ値を示す地点をふるいにかけるとか、できるのかもしれない。

あとは、湿地のような少なすぎるクラスのデータを得るべく、富山や福井もちょっと含んでみたりとか。今回さぼった不均衡データの対処については既に書いた通りです。

いくつかの小区域に区分?

効果があるかは不明だけど、例えば、加賀平地・山地、能登とか3区分くらいに小分けにするとよさそうな気もした。計算も早くなるだろうし。

SVM

ランダムフォレストじゃないアルゴリズムを使ってみたらいいんじゃない?と、いうわけです。SVMをちょっとやってみたんだけど、実際正解率は少し高い(RF 87.4% vs SVM 87.7%)。ランダムフォレストだとよくあった絶対にないところに出てくる耕作地がないっていう違いがあったけど、あとはどうなんだろう?↓

f:id:camarao:20210815115507p:plainf:id:camarao:20210815115515p:plainf:id:camarao:20210815115523p:plain
"左: 金沢、中: 白山、右: 粟津"
GGEでSVMを使う前処理としてデータの標準化が必要なので、その分のコードを書き足さないといけないこと、時間がかかること(RF 30分 vs SVM 2時間+)がちょっと面倒。標準化の作業も結構重そうだったので、各バンドの平均と標準偏差は別途計算して書き出し、別のスクリプトでそれを読み込んで処理することにした。

// レイヤの準備
var trainingLayers = s1.addBands(s2).addBands(nightLight).addBands(dem).select(bands_all);

// 平均と標準偏差のReducerをくっつける
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

// くっつけたReducerでバンドの平均と標準偏差を計算.
var stats = trainingLayers.reduceRegion({
  reducer: reducers,
  scale: 10,
  tileScale: 16,
  bestEffort: true,
});

// 出てくるのがDictionaryなので、地理情報を持たないFeatureに(FeatureCollectionじゃないとcsvに書き出せない)
var feature = ee.Feature(null, stats);

// Wrap the Feature in a FeatureCollection for export.
var featureCollection = ee.FeatureCollection([feature]);

Export.table.toDrive(featureCollection,'LayerMeanSd');
var LayerMeanSd =  ee.FeatureCollection('アップロードしたcsvファイル').first().toDictionary();

var standardise = function(band){
  var stringBand = ee.String(band)
  var mean = LayerMeanSd.get(stringBand.cat('_mean'));
  var sd = LayerMeanSd.get(stringBand.cat('_stdDev'));
  var img = trainingLayers.select(stringBand)
  
  return ((img.subtract(ee.Image.constant(mean))).divide(ee.Image.constant(sd)));
};

var standardisedLayer = ee.ImageCollection.fromImages(
                              trainingLayers.bandNames().map(standardise))
                          .toBands()
                          .rename(trainingLayers.bandNames());

// 教師データ取得
var trainingData = standardisedLayer.sampleRegions({
                                collection: pointsInt,
                                properties: [label],
                                scale: 10,
                                tileScale: 16 // this is to save memory. increase number to save.
                                });
                                
// 分類機作成
var classifier = ee.Classifier.libsvm({kernelType: 'RBF', //ガウジアンカーネル
                                        gamma: 0.03571 ,    //Rで求めたガンマ
                                        cost: 1}        //Rで求めたコスト
                                        ).train(trainingData, label, bands_all);

// 分類
var classifiedSVM = standardisedLayer.classify(classifier);

まとめ

勉強になりました。やはり手を動かして気づいたこと、学ぶことは多いです。 今回のコードはこちら。長くなったので複数に分けてexportsとrequireを使って参照してます。

Sentinel-1 前処理: https://code.earthengine.google.com/5d2b7597da09a617359b2185285bbfad

Sentinel-2 前処理: https://code.earthengine.google.com/dd330eb74dd956224d2e839d4597a526

地形データ前処理: https://code.earthengine.google.com/0365f239fab9967f09476a383d6a42ae

Rで作業するためのポイントのデータ取得https://code.earthengine.google.com/418e0471c4936e3697b4a62c38366ff9

メインのコード: https://code.earthengine.google.com/52c3323b874235dfa2132e05b3ffc7a0

あと、プロの作った土地利用図はこちら↓

www.eorc.jaxa.jp