海老ノート

Google Earth Engine 苦闘の記録

SARの勉強１

SAR Sentinel-1

とっつきにくいイメージのSAR

前からいじってみたかったんだけど、よくわかんないので敬遠していたSAR。教科書や参考になるサイトを見ても数式やギリシア文字が登場し、ふーっとため息をつくばかり。成長を拒否するような自分には困ったものです。いろいろ扱いに注意がいるようだし、写真のように見えるLandsatやSentinel-2と比べると敷居があるなと。

Google Earth EngineでのSAR画像

Google Earth Engine にはSentinel-1とPALSAR-2/PALSARの全球モザイクが入っています。Sentinel-1はGRDは入っているけどSLCはなし。なので干渉SARを使うことはできないようです（GEEでは虚数が使えないからと言ってた）。

教材

Sentinel-1の扱いについてはチュートリアルのほか、Youtubeでのイントロ的な説明もあるのでここから勉強を始めることに。

www.youtube.com

その他、チュートリアルと動画で紹介のあった教科書やサイトなど勉強資料はこちら。

Layman's Interpretation Guide to L-band and C-band Synthetic Aperture Radar Data. Version 2.0

servirglobal.net

おなじみ宙畑先生も。

sorabatake.jp

こちらもおなじみとらりもん先生。

pen.envr.tsukuba.ac.jp

千葉大近藤先生

www.cr.chiba-u.jp

学んだことまとめ

分析に使うのは強度（InSARでないとき）

Sentinel-1でもPALSARでも、衛星からマイクロ波を照射→マイクロ波が地表で散乱→衛星に帰ってきて受信という流れで、帰ってきた信号の強度を測っている。で、後方散乱係数とか、それをデシベル単位にしたものがデータとして使える。この強度は地表の凸凹さツルツルさによって変わる。なぜかというとツルツルの表面（水面とか裸地）では斜め上方から照射されるマイクロ波は反射されて衛星方向に戻ってこないけど、樹木のように散乱する場合は一部が衛星方向に帰ってくるから。

〇バンドで、適した観察対象が変わる。

Sentinel-1はＣバンド、PALSAR/PALSAR-2はＬバンド。両者は周波数が違うから、得られるデータが違う。一般に、波長と同じくらいの大きさもしくはそれ以上の大きさの「物」にぶつかるとマイクロ波は反射/拡散されるが、それより小さいと通り抜けるらしい1。PALSAR-2の波長は23.5cmなので木の幹くらいの大きさ、一方Sentinel-1の波長は5.6cmで木の葉や小枝で、場合によっては草の葉でも検出できる。成熟した森林の場合バイオマスの多くは高木の幹に蓄積されるので、PALSAR-2は熱帯林のバイオマス推定に威力を発揮しているとのこと。

ノイズとエラー各種

光学センサーと異なる仕組みからノイズの種類も違っている。ここの理解はマイクロ波衛星データ活用のポイントなんだろうなと思った。

スペックルノイズ

SAR画像を見るとブツブツとごま塩状になっているけど、これがSpeckle。発生する理由は、いろんな方向に散乱したマイクロ波がある確率で強め合う/弱めあうことがあるから。GEEチュートリアルにも詳しい説明あり。これに対処するために各種フィルタが提唱されていて、GEEで使えるコードがGitHubで公開されている。

www.mdpi.com

単純な、例えば5 x 5の窓での平均のフィルタでは画像がぼやける。ので、ノイズは少なく、だけどエッジはくっきりとすべく考えられた各種フィルタがある。主に空間（n x n窓での統計量を使って畳み込んだ値と元の値を重みづけして加算とか）、場合によってはそれに時間方向（前後のデータからn x n窓平均に対する差を見て補正）のフィルタ組み合わせて使うよう。また、GEEを使った論文では複数のSAR画像をImageCollectionとしてmedian()やmean()をとったり、圃場などある面積の平均値をとったりすることで対処しているのもあった。

地形の影響

SARデータは衛星から斜め下方向を見ているので、山があると陰側のデータは得られないし、衛星との向きと斜面のなす角度によってマイクロ波の入射角が変わることで、計測値が影響をうける。また、同じ理由から、画像にすると山の近い側の斜面は実際よりも近くにあるようになっちゃう。

これに対処するためにDEMを使ってTerrain Flatteningとよばれる処理をする。GEEではスペックルノイズで紹介したコードでできる。

水分の影響

土壌や植物に含まれる水分によってマイクロ波の減衰・拡散されやすさが変わる。必要であれば、同じ時期の画像をつかったり、雨の多い時期の画像は避けるとか対応する。

入射角

衛星が斜め方向から下を見るとき一枚の画像で近い側と通り側ができる。それにより入射角がかわるので計測値も影響をうける。一定程度はTerrain Flatteningで緩和できる。また、VH/VVの比を取るなどすると影響が低減できるらしい。

計測値と植被率やバイオマスの関係が飽和したりする

これは時の通りで、ある一定の植被率を超えると計測値との関係が飽和することが多い。また、植被率が低すぎると計測値の誤差が大きくなり（ノイズの影響が大きくなるから？）推定値がブレブレになることが多いらしい。

というわけで、頑張って学習しました。ノイズの性格がよくわかっていないので、そこら辺をもうちょっと勉強したいと思います。

  1. 「透過」という表現がされている場合があってX線的なイメージでちょっと混乱した。たぶん、波長が長いと波の回折で障害物を回り込んで奥まで届くということなんじゃないかと思う。↩

石川チャレンジ５：まとめとおまけ

分類 まとめ

そんな簡単にはいかない

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

分析に時間がかかる

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

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

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

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

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

土地利用区分の偏り

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

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

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

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

www.mdpi.com

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

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

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

もっと良くするには？

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

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

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

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

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

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

石川チャレンジ４：GoogleEarthEngineで地図化

分類

やっとGEEに到達

あーだこーだとやって、これで行きましょうと行くまでの体感が長かった。実際はそんなにかかってないんだけど。 というわけで、GEEで分類します。見本があるのでそれに倣って書けば大丈夫。

developers.google.com

実際はこう。

//ポイントデータからラスタ情報を取得して教師データを作成
var trainingData = trainingLayers.sampleRegions({
                                collection: pointsInt,  // 教師データのポイントのFeatureCollection（シェープファイル）
                                properties: [label],   // 分類するクラスの列名
                                scale: 10, // 縮尺10m単位で取得
                                tileScale: 16 // 数を大きくしておくと（最大16）時間切れになる可能性がへる.
                                });
                                
// 分類機の作成
var classifier = ee.Classifier.smileRandomForest({
                                  numberOfTrees: 500, // ntree
                                  variablesPerSplit:5,    // Rの結果からのmtry
                                  minLeafPopulation: 1,
                                  bagFraction: 0.632, // RのrandamForestパッケージのデフォルトはこの数値らしい
                                }).train(trainingData, label, bands_all);
                                
// 分類機とラスタのレイヤを使って全域を分類
var classifiedRF = trainingLayers.classify(classifier);

結果をみる

石川県全体

あまり違いがわからない。

凡例は濃緑：高木林、薄緑：低木林、青緑：湿地、黄色：草地、オレンジ：耕作地、肌色：水田、薄赤：裸地、灰色：人工物、青：水域。

f:id:camarao:20210815102015p:plainf:id:camarao:20210815102025p:plainf:id:camarao:20210815102028p:plainf:id:camarao:20210815102033p:plain
左からB2-4, 8,8A, 11-12、+NDVI/NDWI、+夜間光と地形、+Sentinel-1

金沢の兼六園とか美術館とかのあたり

お手本はこちら↓

f:id:camarao:20210815102428p:plainf:id:camarao:20210815102429p:plain
左: Google Earth 右:環境省植生図からつくった図
結果↓
f:id:camarao:20210815102730p:plainf:id:camarao:20210815102735p:plain
左: B2-4, 8,8A, 11-12、右: +NDVI/NDWI、
f:id:camarao:20210815102741p:plainf:id:camarao:20210815102745p:plain
左: +夜間光と地形、右: +Sentinel-1

思ったよりいい感じ。変数（レイヤの種類）が増えるにつれてまとまっていく感じがある。耕作地と予測された場所（オレンジ色）が縁取りのようにあるのが気になる。

凡例は濃緑：高木林、薄緑：低木林、青緑：湿地、黄色：草地、オレンジ：耕作地、肌色：水田、薄赤：裸地、灰色：人工物、青：水域。

白山のメインのあたり

お手本

f:id:camarao:20210815103226p:plainf:id:camarao:20210815103228p:plainf:id:camarao:20210815103231p:plain
左: Google Earth、中：国土地理院地形図（地理院タイル）、右:環境省植生図からつくった図
結果↓
左からB2-4, 8,8A, 11-12、+NDVI/NDWI、+夜間光と地形、+Sentinel-1
f:id:camarao:20210815103443p:plainf:id:camarao:20210815103446p:plainf:id:camarao:20210815103455p:plainf:id:camarao:20210815103503p:plain

どれが正しいのかよくわかんないから、なんともいえない。一番左の図は低木（＝ハイマツとか）がほとんど描写されていない。

小松市粟津町周辺

お手本

f:id:camarao:20210815103926p:plainf:id:camarao:20210815103929p:plain
左: Google Earth 右:環境省植生図からつくった図

結果↓

f:id:camarao:20210815104007p:plainf:id:camarao:20210815104016p:plain
左: B2-4, 8,8A, 11-12、右: +NDVI/NDWI
f:id:camarao:20210815104024p:plainf:id:camarao:20210815104028p:plain
左: +夜間光と地形、右: +Sentinel-1

感想は金沢に同じ。

数値的な出来の判定
GoogleEarthEngineで別に発生させた8000点を使った評価

正解率は86.9%だった。

開放水域 耕作地 高木林 湿地 人工物 水田 草地 低木林 裸地
開放水域 1891 1 2 0 0 0 1 1
耕作地 2 52 39 0 24 45 10 4 9
高木林 20 65 4026 0 35 124 99 118 21
湿地 1 0 0 0 1 1 0 0 0
人工物 3 25 20 0 353 94 5 0 6
水田 5 18 25 3 27 526 5 2 2
草地 1 1 20 0 3 4 48 2 2
低木林 0 0 38 0 0 0 6 27 1
裸地 20 4 37 1 22 14 8 5 26
Rで別に発生させた7448点を使った評価

正解率は87.3％ (CI: 86.48 - 88.01%)だった。表は省略。前回作ったグラフと見比べるとGGEの結果は低い気がする。大体あってるけど、同じ結果と言っていいのかな？は自信なし。（続く）

石川チャレンジ３：Rで実験

分類

教師データはいくつ？

とりあえず10000のランダムなポイントでとしてはじめたんでした。で、サンプリングしてRに持ってきてみたんですけど、データの偏りがすごい。

開放水域 耕作地 高木林 湿地 人工物 水田 草地 低木林 裸地
2723 225 5949 8 627 1083 236 231 82

いやー、どう考えても8点（湿地）だけではどうしようもない。なので100,000くらいに増やしてみたんです。が、そしたら後々のランダムフォレストの分類機を作る手順でデータ多すぎですエラーがでた。ちゃんと説明で、そういうことは止めてねと書いてあるので、するべきではありませんでした。

これを踏まえて、元のデータに湿地や裸地等のデータの少ないカテゴリ補うようポイントを追加してみた。これだと教師データが多くなりすぎないはず（計14,600）↓

開放水域 耕作地 高木林 湿地 人工物 水田 草地 低木林 裸地
2723 885 5949 89 1286 1083 896 891 798

で、何点くらいあったらそれなりの結果がでるのか見てみる。テストデータは植生図から別途取得したランダムな点。合計14,600のデータから適当な数を抽出して、Rでランダムフォレストを回してみる。

f:id:camarao:20210814134429p:plainf:id:camarao:20210814134448p:plain
正解率（左）とカッパ係数（右）
f:id:camarao:20210814135942j:plain
各クラスのBalanced Accuracy

右上がりのグラフだけど、8,000点くらい1あればそれなりになりそう。それ以上は努力に見合わなそうというざっくりとした見立てで、この14,600→8,000点をランダムに抽出して教師データとすることに。内訳はこのとおり。振り返ると、耕作地、湿地、低木林はもっと多めに取ったほうがよかったっぽい。

開放水域 耕作地 高木林 湿地 人工物 水田 草地 低木林 裸地
1483 487 3238 56 693 591 494 488 462

Rで実験

今回の目的である、いろんなデータをレイヤに入れるとどのくらいうまく分類できるかな？というのをRで実験。 普通にRのrandomForestパッケージを使いました。ntreeは500、それといちおうtuneRFでmtryを求めて入れ込んでます（スクリプトは省略）。

バンド B2-4, 8,8A, 11-12のメディアン +NDVIとNDWIのメディアン +NDVIとNDWIの10,90分位 +夜間光 +地形データ +Sentinel1
正解率 0.854 0.854 0.859 0.864 0.868 0.871
カッパ係数 0.767 0.768 0.777 0.783 0.792 0.796

なんか変化が地味。ガーっとよくなるのかなと期待してたんだけど。 全部入りから作ったRで出せる変数の重要度のよくあるグラフはこうなった↓ f:id:camarao:20210814141649p:plain Sentinel1はそれほど役に立ってなさそう。B8やらNDVIが大事っていうのは納得（続く）。

  1. 試したところ10,000点くらいまではGEE上での分析が可能だった。↩

石川チャレンジ２：なかなか始まらない

分類

そう簡単にはいかない

いや、わかってましたよ。そんなサクサクはいかないだろうなーと。でも思った以上に考えることやらケースが多くて、勉強・経験不足が露呈します。どこらへんがなのか、というのは追々として、これまでのあらすじを書いてみましょうか。

作業手順

  1. QGISで植生図のカテゴリをリクラス
  2. QGISでランダムにポイント発生させて、ポイントに植生クラスを付与、シェープファイルに書き出し
  3. シェープファイルをGEEに読み込む
  4. GEEで特徴量として使う衛星データを整理して、読み込んだシェープファイルにそのデータをくっつける
  5. シェープファイルの各点のデータをcsvで書き出し。ダウンロードしてRでランダムフォレストの設定を調整
  6. 調整した数値などを使って、GEE上でランダムフォレスト。完成。

というわけで、結構手数があるわけです。で、5まで行ったと思ったら違ってたから1に戻るとかですね、そういうのがたびたび起きるわけで、あちゃーとなるわけです。モチベーションは下がるけど、一度通った道は早く感じるというのか、スクリプトはあるし思ったより楽なんだけど、モチベーションは下がる。慣れていれば落とし穴やらひっかけ問題にかからず、サクッと行くんでしょうが…。

GEEでの衛星データの整理

QGIS周りの作業は大胆に割愛。シェープファイルのアップロード法はこちらの説明の通り。

Sentinel2の前処理

GEEのサイトを見たらpythonの雲除去コードがあったので、こちらを借用。これをjavascriptに書き換えます（っつってもdefをvarに変えるとかそれくらい）。最後に雲を除去する関数を実行して、そこからNDVIとNDWI、各バンドのmedian、10分位、90分位を計算するようにしています。日付を微妙に設定しているのは、冬は雪で真っ白かもしれないから、そうじゃない時期だけにしたいから↓

https://code.earthengine.google.com/eccc6a6b863ab6ea8c3522cdc085f4dc

地形と夜間光も使ってみたい

地形データはJAXAのALOSから、最近流行ってるっぽい夜間光データはNOAAのものを使いました。JAXAのALOSは使い方にひと手間いるようですが、説明があるのでそのとおりやればできた。

// dem
var dsm = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").select('DSM').filterBounds(aoi);
var elevation = dsm.mosaic() // 標高

var proj = dsm.first().select(0).projection();
var slope = ee.Terrain.slope(dsm.mosaic()
                             .setDefaultProjection(proj))
                             .rename('slope'); // 傾斜

var aspect = ee.Terrain.aspect(dsm.mosaic()
                             .setDefaultProjection(proj))
                             .rename('aspect');　// 斜面方位

var mtpi = ee.Image('CSP/ERGo/1_0/Global/ALOS_mTPI').select('AVE').rename('mtpi');

var dem = elevation
             .addBands(slope)
             .addBands(aspect)
             .addBands(mtpi);
Sentinel1はよくわからないのでそのまま

Sentinel1の地形補正ですが、こちらはSRTMのDEMを使って入るらしい、というところまでわかりましたが、あとはちょっと…タイルごとに書き出す部分を削除したのを除いてVenter先生のコードをほぼそのまま使ってみます。出てくるのはmasterstackというマルチバンドのラスタで、ここに「vh、vv、dpol」「メディアン、標準偏差」、「軌道上向き、下向き」の3×2×2=12通りのデータが入っています。

https://code.earthengine.google.com/8b7f12488967850eedbcb2930c439a82

見た感じ大体OKでしたが、白山の周りではちょっとしんどそうな谷間がありました。あと、海上の船？みたいな点が出てくるのも気になります。Sentinel1を弄り回すのはそのうちやってみたい。

sampleRegionsでデータを取得

公式分類ガイドと似たような感じでデータをとって、それをエクスポート。

// Overlay the points on the imagery to get training.
var sampledPoints = imageSentinel.select(bands).sampleRegions({  // imageSentinelというラスタのbandsというバンドから取得
  collection: points, // 点データのシェープファイルの指定
  scale: 10
});

// Google Driveに書き出し
Export.table.toDrive({
  collection: sampledPoints,
  description: 'sentinel2', // 書き出すファイル名
  fileFormat: 'CSV'
});

ということで、ここまでで手順5の半分まできた。こうしてようやく分析に入れるわけです。（続く）

石川チャレンジ１：土地被覆分類への道（予告）

分類 方向性

リモセンと言ったら植生分類

GEEのチュートリアルその他で土地被覆分類の紹介は分かりやすくされていて、その通りに手軽にできます。分類のアルゴリズムも決定木、SVM、ランダムフォレスト…と各種揃っています。分類はQGISだとかRだとか無料のソフトを使ってもできるけど、GEEだとダウンロードや前処理の手間がいらないのがメリット。とはいえ、まああえて探求しようというほどには興味がなかったんだけど、これを見ていたらやってみたくなった↓。 www.mdpi.com

GEEならば広範囲の分類もパワーで押し切ってできちゃうし、いろんなデータを手軽にくっつけたりできるのかと。というわけで、今回のテーマは、

  1. 地理的な解像度が高い分類を広範囲にできるのかな？っていうのと、

  2. Sentinel1の画像その他を組み合わせるとどうなるんだろうかな？ という２点です。元論文だと、8個のカテゴリーを90％くらいの正答率で分類できたそうです。すげー。

始める前に

山も海も里もある石川にしてみる

どこがいいかなと思ったんですが、広すぎない県、海から高山まである県ということで思いついた石川にしてみました。ただ、舳倉島その他沖合の島や岩礁は除いて、本州近辺だけ。

教師データはどうするか

教師データがないので、どうやって確保するかですが、とりあえず環境省植生図のシェープファイルを使って、ポリゴンからポイントを生成することにしました。描写スケールの問題で、正しい教師データが得られないポイントもあるけど、そこは目をつぶります。元の論文では教師データの数を5,000から50,000に増やしてみても分類精度は3％しか上がらなかったといっているから、とりあえず10,000ポイントくらいを作ってみます。

カテゴリー分け

環境省植生図は非常に細かい分類がなされています。非常に素晴らしい仕事です。それをリモセンで分類となると、ちょっと難しそうなので凡例を眺めて適当に分類しなおします。生態学的な妥当性というより、リモセンでできそうな方に歩み寄ったと。元論文では耕作地と一つにされているけど、そこは水田とそれ例外を耕作地に分類するという9分類にしてみましょうか。果樹園は森林と耕作地のどっち？とかとか細かく言いはじめるとよくわかんないんですが。この辺は様子をみて変えるかも。

というわけで、しばらく様子見ながらやってみることにします。

回帰分析の結果はどうなってるか

時系列 回帰

Google Earth Engineで使える線形回帰は４種類

GEEのガイドでは次の4種類のメソッドを使った線形回帰が紹介されている。これまでも使ったけど、係数以外の返り値がどうなってるのか詳しく見たい、似たように見えるけど、どこが違うのか見たいというのがやりたいこと。

  1. ee.Reducer.linearFit()
  2. ee.Reducer.linearRegression()
  3. ee.Reducer.robustLinearRegression()
  4. ee.Reducer.ridgeRegression()

簡単な違いの説明

単回帰分析 ee.Reducer.linearFit()

単回帰分析専門。引数はX（説明変数のバンド）とY（目的変数バンド）を指定すればOK。返り値は'scale'（傾き）と'offset'（切片）の２つのバンドを持つImage。

重回帰分析もできる ee.Reducer.linearRegression()

一般線形モデル。引数はX（説明変数のバンド数）とY（目的変数のバンド数）の二つ。デフォルトでXに定数項が入っていないので含めるのを忘れないこと。返り値は2バンドのArray Imageで、1バンド目の ['coefficients'] に推定された切片と係数が、2バンド目の ['residuals'] に目的変数のRMSEが入っている。

決定係数や検定の結果は入っていない。決定係数はRMSEと分散から計算できる↓

// R squared from RMSE
var rSquared = function(residualImage, coefficientsImage, imageCollection){
  var filteredIC = imageCollection.select('celsius');

  var n = filteredIC.count();
  var p = ee.Image.constant(coefficientsImage.bandNames().length().subtract(1));
  var one = ee.Image.constant(1);

  var variance = filteredIC.reduce(ee.Reducer.variance()); // .sampleVariance()だと不偏標本分散のようで値が異なる
  var rmseSq = residualImage.multiply(residualImage);

  var rSquared = one.subtract(rmseSq.divide(variance)); // 1 - rmseの二乗/分散
  var adjRSquared = ((n.subtract(one)).multiply(rSquared).subtract(p)).divide (n.subtract(p).subtract(one)); // ややこしい回帰のときはこれじゃあ力不足かも。

  return residualImage.addBands(rSquared)
                        .addBands(adjRSquared)
                        .rename(['RMSE','R-sqared','adjusted R-squared']);
};

これでRのlmでのsummaryの出力と一致。決定係数と自由度調整済み決定係数の計算式はwikipediaと他のブログに教えてもらいました。

Rのlmの決定係数R^2が実数範囲で負をとりえない話 - 備忘ログ

Rのsummary関数を使わずに、決定係数を求めてみる | while(isプログラマ)

Xに説明変数をn乗した項を含めて多項式回帰にしたり、その他三角関数や対数関数の回帰も可能。 気温データを使ったここまでの練習コードはこちら https://code.earthengine.google.com/1c414c3696fa60bf1e620715f2ea8faf

ロバスト回帰分析 ee.Reducer.robustLinearRegression()

外れ値の影響を受けにくいロバスト回帰分析。ほとんどee.Reducer.linearRegression()と同じだけど、引数にベータが加わるらしい。詳しくは論文を読まないとわからなそう（読んでもわからないかも）だったので、調べていません。

リッジ回帰分析 ee.Reducer.ridgeRegression()

リッジ回帰の引数はX（説明変数のバンド数）、Y（目的変数のバンド数）、ラムダ（正則化項のパラメータ、デフォルト値0.1）の３つ。.linearRegression()と違い、Xに定数項を含める必要はない。返り値は3バンドのArray Imageで、1バンド目の ['coefficients'] に推定された切片と係数が、2バンド目の ['residuals'] にはRMSEではなく残差平方和の平方根1が入っている。3バンド目の['pValue'] は偏回帰係数をt検定したp値2が入っている。

ラムダを0にすれば、普通の線形回帰と同じ。なので、t検定の結果が欲しい場合は.linearRegression()ではなく.ridgeRegression()を使えば出てくる。residulasの出力がよくわかんないけど、残差平方和の平方根からも分散を使えば決定係数は計算できる。だいたいは.linearRegression()と同じなので省略、こちらをどうぞ↓

ridgeRegression()でλを0にした場合：https://code.earthengine.google.com/6f256e27b04e5dcc252f8ab9363ac90d

その他の回帰的なもの

チュートリアルにノンパラメトリックな方法が紹介されている。ランダムフォレスト回帰とかできるのかもしれないけど、探せなかった。調べたらできるのかも。 developers.google.com

回帰の分散分析（F検定）はやりかたがわからない

調べてもわからなかったのが、回帰式自体が有意かどうかのF検定の方法。F値と自由度は求められるけど、その値がF分布でどうなのよという部分がわからない。自作するとしたらベータ分布やらガンマ分布やら出てきてちょっと…という感じ。先人の知恵があるのか。そのうち調べたい3。

  1. リファレンスを読むと"Additional outputs are a vector of the root mean square of the residuals of each dependent variable …“と書いてあるけど、Rで確かめたところmeanにはなっていないみたいで謎。↩

  2. 数学的な言い回しがわかんないんだけど、係数のt検定ということを言いたい。↩

  3. さっき見たら ee.Image.erf() というのがあったんだけど、これでいけるのかな？わかんない。↩