GEEでOtsuの二値化
できるはずだけど面倒くさそうと思っていた
画像処理でおなじみ大津の二値化。考え方はシンプルだし、GoogleEarthEngineでできるんだろうけどコード書くの面倒だなと。以前八郎潟の湛水の有無を見るのにSARデータ使った時にも、二値化したいと思っていたけど適当にごまかした。が、先日検索してみると随分前にコードが紹介されてました、しかも公式Nick先生。というわけで、こちらをお手本に試します。
お手本コードのメモ
このコードのポイントは ee.Reducer.histogram()を使うところ。お手本では2~256クラスに分けるよう引数で指定して、ヒストグラムの各階級ごとに平均とデータ数をまとめてBSSと呼んでいる分割の良さを表す指標を計算している。
これまで面倒と思ってたのは長いリストを作って処理しないといけないところ。0~255のDNくらいだとそこまで大変ではないだろうけど、NDWIを0.001刻みでとなると処理に時間かかりそうだなと。そこのところ、このお手本ではhistogramの関数をうまく使うことでざっくりと最大256段階に分けている。
お手本へのリンクhttps://code.earthengine.google.com/6f3abb73b6642198485e36000024827a
もっと詳細にみたい
なるほどと思いつつ、手を加えたくなっていじってみる。Nick先生もざっくりだよ、とおっしゃっているし:
We’re not going to search exhaustively. We’re just going to search over the thresholds that are represented by the bins in a histogram. The advantage of that approach is that it only requires a single pass over the data.
それなら、ざっくりから細かくしてみようかと。蛇足だったけど。
力ずくのコード
単純にお手本コードのee.Reducer.histogram()の一つ目の引数を大きくすれば閾値のレベルは細かくなる。ただ、それじゃあやった気がしないというわけで、お手本で求められた閾値の前後を細かく区切ってBSSがどうなるかを調べることにします。リストにした閾値候補ごとにmaskを作って、地道にBSSを計算してみる。大体のつくりはお手本コードをお借りしてます。nClassを小さくして、計算を数回繰り返せば早くて細かい計算ができるはず。
var nClass = 30 // 閾値候補の数。30くらいなら我慢のいらない速さで回る。 var bucketWidth = ee.Dictionary(histogram.get('B5_histogram')).get('bucketWidth'); var thresholdList = ee.List.sequence({start: threshold.subtract(bucketWidth), end: threshold.add(bucketWidth), count: nClass}); //閾値候補のリスト print(thresholdList); //今回メインの関数 var detailedOtsu = function(image){ var globalMean = image.reduceRegion({ reducer: ee.Reducer.mean(), geometry: polygon, scale: 30, bestEffort: true }).values().getNumber(0) // 全体の平均 var funcBss = thresholdList.map(function(i){ var greaterThanMask = image.gt(ee.Image.constant(i)) var greaterThan = image.mask(greaterThanMask) //閾値候補超のピクセル抽出 var lessThanEq = image.mask(greaterThanMask.eq(0)) //閾値候補以下のピクセル抽出 var funcStats = function(image){ var intraMean = image.reduceRegion({ reducer:ee.Reducer.mean(), geometry:polygon, scale:30, bestEffort: true }).values().getNumber(0) //マスクしたピクセルの平均 var intraCount = image.reduceRegion({ reducer: ee.Reducer.count(), geometry: polygon, scale: 30, bestEffort: true }).values().getNumber(0) //マスクしたピクセル数 return intraCount.multiply(intraMean.subtract(globalMean).pow(2))} return funcStats(greaterThan).add(funcStats(lessThanEq))} ) print(ui.Chart.array.values(ee.Array(funcBss), 0, thresholdList)) //BSSのグラフ作成 return ee.Array(thresholdList).sort(funcBss).get([-1]) //BSS最大の閾値候補の選択 }
どっちでもいい感じ
お手本コードでの閾値は8127.50。一方細かく見た方では8184.88。グラフだとこういう感じ。
閾値が60違うとどう変わるのか?
あまり違わない。画像で見てみるとほぼ同じ。
結局水と陸が混ざっているピクセルの扱いの問題みたいから、もっと解像度が高い画像か現地データがないと精度の評価は難しそう。
k-meansでもできるのでは?
GEEでも使えるお手軽分類、k-meansを2クラスでやってもできるんじゃないかと思ってやってみたら、これもほぼ同じ結果になった。
k-meansがいいのはGEEで書くコードが断然短くて楽。悪いのはどちらが水でどちらが陸のクラスなのかはランダムに決まるから、時系列で多くの画像から水域を抽出する処理をするような場合はひと手間が必要になるはずというところか。
まとめ:どこまで頑張ってどこで妥協するか
そういうことです。Landsat8の画像でざっと見る限りではどっちでも、というところ。値を取得する範囲となるポリゴンの位置によっても結構変わるだろうし。より詳細に突き詰めれば細かくなるんだろうけど、それがどの程度意味あるかっていうのはわからない、正解データがあれば数値で評価できるんだろうけど。バランスは大事です。
今回のコードはこちら: https://code.earthengine.google.com/b95c53a32bca944425f0f7d3b08e2f5d
GEEで分布推定モデル2:MaxEnt
SDM界の iris こと Bradypus
とりあえず動かせるかみるために、分布推定モデルでお題になるミツユビナマケモノ Bradypus のデータを使ってやってみる。 お手本はRのライブラリ dismo を使ったSDM各種を解説したページから。
データの準備
ミツユビナマケモノの分布地の座標は元祖のページからcsvがダウンロードできる(ほかでもいろいろある)。それをGoogle Earth EngineのAssetsにアップロードする。
環境データはお手本にそっておなじみBioClimから8レイヤと生物区系(バイオーム)。生物区系のデータはGEEのデータカタログにあったEco Regionを使ったけど、これはお手本データと微妙に違うかもしれない。次の手順は環境データ間の相関をみるだけど省略。お手本データがそれをやっているはずということで。1 これらをまとめて一つのee.Imageにすれば準備は一区切り。
// 生物区系 var ecoreg = ee.FeatureCollection("RESOLVE/ECOREGIONS/2017") .filter(ee.Filter.notNull(['BIOME_NUM'])) .reduceToImage({ // ベクタからラスタに変換 properties: ['BIOME_NUM'], reducer: ee.Reducer.first()}) .rename('ecoreg'); // bioClim から使うバンドを抽出して生物区系とスタック。 // そのあと使う部分をクリップして、ラスタをESPG4326、解像度1㎞にそろえる。 var envLayer = ee.Image("WORLDCLIM/V1/BIO").select(['bio01', 'bio05', 'bio06', 'bio07', 'bio08', 'bio12', 'bio16', 'bio17']) .addBands(ecoreg) .clip(AOI) .reproject({crs:'EPSG:4326', scale: 1000}) print(envLayer)
ラスタを揃えるのだけでもローカルでやると結構時間がかかる(ダウンロードも入れるともっと)から、この作業をサクッとできるだけでもGEEは便利。
Pseudo Absenceの点の発生
MaxEntでおなじみのPseudo Absence(訳は偽不在データ?)をee.Image.sample()を使って発生させる。発生させる範囲はお手本にそって調査範囲+その周り2の陸域とした。在地点の周囲○kmで…とかもっと複雑な方法もあるけど、それはお手本コードの方を参照してほしい。
発生させる点の数を決めるnumPixelsを4500としたのは、海の上の点(無効になる)が入っちゃうからか、使える点の数を1000位にするにはこのくらいが必要だったため。
// mask作り。在地点のピクセルと水域。 var mask = presencePoints.reduceToImage({ properties: ['presence'], reducer: ee.Reducer.first() }) .reproject('EPSG:4326', null, ee.Number(1000)); var watermask = ee.Image("CGIAR/SRTM90_V4").select('elevation').gte(0) // AOIとその周りに偽不在のポイントを発生。 var AreaForPA = mask.mask().updateMask(watermask).clip(AOI.buffer(800000)); Map.addLayer(AreaForPA, {},'Area to create pseudo-absences', 0); var pseudoAbsPoints = AreaForPA.sample({region: null, scale: 1000, numPixels: 4500, seed: 1, geometries: true, tileScale: 16}); pseudoAbsPoints = pseudoAbsPoints.map(function(feat){ return feat.set('presence', 0); });
で、このポイントデータを在地点のものとくっつけた後、ee.Image.sampleRegions()を使ってラスタの情報を取得すれば準備完了。
結果
後は前回書いた通りにMaxEntを走らせればOK。.explain()を見れば結果が出てくる。 各環境要因の重要度はContributionsに、テストデータのAUCはTest AUCに入っている。重要度はお手本と結構違う。AUCは似たような値(0.837 vs 0.834)。 地図での出力はこんな感じ。 だいたい似てる。
二値化してみる
これも色々説はあるみたいけど、お手本のとおり感度と特異度の和が最大になる点を閾値として2値化してみる。この値は.explain()で出てくるなかのThresholds -> Minimum test sensitivity plus specificityに入っている。で、これが入れ子になったDictionaryに入っていて取り出し方がわからず面倒だった。この閾値の値はお手本と結構違っていた(0.32 vs 0.17)。
まとめ
GEEでいいのは手軽にできること。
GEEでできないのは各環境要因のresponceのグラフが出てこない、AUCのグラフも出てこない、重要度のグラフが出てこない。とにかくグラフは出てこない。 responceのグラフを作るのに必要な数値は.explain()で出てくるLambdasに入っているのかもしれない(よくわかってない)。AUCは頑張れば似たようなグラフをGEEでも作れる(やってみた)。重要度のグラフも書けなくはない。
GEEとローカルを組み合わせるのがいいのかも
全体としてローカルでやった方が親切な気がした。パラメータの設定とか結局いろいろいじるだろうから、その辺は時間がかかってもローカルでやった方がいいんじゃないかなと思う。
例えば、対象種が1種の場合はデータの準備をGEEでやってラスタをエクスポートして、それを使ってローカルで分析を進めるのがいいのかなと思った。今回は端折ったけどデータの準備で環境要因間の相関やらVIFやらの部分がGEEでは手間がかかるけど、Rだと簡単なのでそういう分業がよさそう。
また、対象種が1種でも複数の環境データに対して予測する場合;
- 例えば年間の降雪量/日のような変動が大きい要因を扱う場合
- 森林伐採で森林面積が変わっている場所での分析など
そういう時は訓練1回で予測が複数回となる。こういうときもモデルの吟味はローカルでやってこれでいけるという設定に至ってからGEE上で走らせればいいんのかなと思った。
コードはこちら
https://code.earthengine.google.com/dad2d00114a4b4e99bdc945907590e5b
一部使わせてもらって大いに参考にしたコードはこちら smithsonian.github.io
-
多いに参考にしたSupplemental Materials - Appendix 1 - Implementation of species distribution models in Google Earth Engineでは、ランダムに点を発生させてGEE上で相関を分析してます。すごい。↩
-
正確には1.25倍だそうですが、目分量で800kmのバッファとした。調査範囲も目分量で合わせたからちょっとずれてる。↩
GEEで分布推定モデル1:できるらしい
できるんだろうなと思っていたら
この間のランダムフォレスト回帰で思ったのが、分布推定モデル(SDM)できそう。で、検索したらありました。
スクリプトとその解説も公開されています。 データクリーニングから始まるかなり長くてがっちりしたスクリプトです。
で、この論文はランダムフォレストで分布推定してます。すごいなあと思ったところはたくさんあるんだけど、特にGoogleEarthEngineでは手間な交差検証をスクリプトの中で仕上げてるのがすごい。調査範囲に配置したグリッドを10グループに分割して交差検証してます。
分布推定モデル(SDM)がGoogleEarthEngineでできたらいいとは思う
GEE上でSDMが動かせるメリットっていくつも思いつく。Crego et al. (2022)でも書いてある通り、NDVIとかの衛星由来のデータをモデルに簡単にぶち込めるのは大きい。あとは複数年や多くの種など沢山のモデルを走らせて、それらの出力結果をつかって考察するような仕事(分布適地の経年変化とか、潜在的なホットスポット解析とか)だといい。
なんだかんだのMaxEnt
生態学分野のSDMといえば、MaxEnt。乱用慎むべし的な話も聞くけれども、それでも何年も多くの人が使っているってことはそれなりに納得感のある結果が得られてるってことなんでしょう、たぶん。で、GoogleEarthEngineのドキュメントを見てたら、ClassifierにそのMaxEntがあった。
正直よくわかっていない
MaxEnt触ったことがある方はご存じの通り、調整できるパラメータが大量にある。
ee.Classifier.amnhMaxent(categoricalNames, outputFormat, autoFeature, linear, quadratic, product, threshold, hinge, hingeThreshold, l2lqThreshold, lq2lqptThreshold, addSamplesToBackground, addAllSamplesToBackground, betaMultiplier, betaHinge, betaLqp, betaCategorical, betaThreshold, extrapolate, doClamp, writeClampGrid, randomTestPoints, seed)
で、わかんないから、とりあえずオートでやってみてAUCみてこんなもんなのかなと(ということをやるから、MaxEntについていろいろ言われるのは知ってる)。で不安で人に聞く/任せるまでが流れ。
というように理解が浅くていので、ここでは動くスクリプトの書き方を記録するにとどめます。レビューも解説論文も出てるので気になる方はそちらを見ていただいたうえ、ご自身の責任でどうぞ。
肝の部分はシンプル
コードの書き方は他のee.Classifilerと似ている。autoFeatureをtrueにしておけば、さらっと動かせる。しかも結構早い。
// モデル作成 // 教師データ(FeatureCollection)をtraining、環境データ(マルチバンドのImage)をenvLayerとした場合 var classifier = ee.Classifier.amnhMaxent({ categoricalNames: ['ecoreg'], // カテゴリカルデータのバンド名を指定 autoFeature: true, // 設定はオートで randomTestPoints: 20 // 分布データの20%をランダムに選んでテストデータにする }) .train({ features: training, classProperty: 'presence', inputProperties: envLayer.bandNames() }); // MaxEnt で推定 var imageClassified = envLayer.classify(classifier); // 予測結果がImageで出力される print(classifier.explain()) //モデルの内容をコンソールに出力
というわけで、次回以降忘備録をつけ足します。
GEEでランダムフォレスト2:関東の人口を予測
人口データの準備
前回のとおり、e-statからデータをダウンロード。具体的には、関東地方7都県の4次メッシュの人口データとshpファイルをダウンロードして、Rを使って各メッシュの人口と中心座標をまとめた1枚のcsvファイルを作りました。
そこから、まずは海上のメッシュを除外。その後メッシュ当たり人口1,000人ごとのクラスに分けて層別サンプリングをして、約1,000メッシュを抽出し、shpファイルに書き出しました。
衛星データの準備
衛星データは深く吟味せず、データカタログをみて使えそうと思ったのを入れてみます。全部で16種類。
500mメッシュということで、MODISの500mのを基本に気温や標高、夜間光(これは解像度約1kmだった)を入れました。MODISデータのSurface Reflectanceは無雪期だけを使うことにして、NDVI、NDWI(水面の方)、NDBI(建物の指標らしい)を計算して追加。気温は月別だったので、年平均、最高、最低の3つにまとめました。
ランダムフォレストをする
人口のshpをassetへアップロードした後ee.FeatureCollection()にして、各種データの数値をsampleRegions()を使って取得、という流れは分類の場合と同じ。今回は8:2で訓練:検証データに分けてやってみました。訓練データ800っていうのは、Rで見てみたらそのくらい以上あればRMSEがあまり変わらなそうだったから。ntreeとかmtryもRでよさげだった値をGEEに設定。
ランダムフォレスト回帰のコードの書き方は前回の通り、.setOutputMode('REGRESSION')をつけるだけでOK。ドキュメントを見ると、他のee.Classifierでもできるものがあるよう。SVMもできるのか?
結果と評価
Google Earth Engineの結果
チュートリアルのコードの通り実行すると、変数の重要度のグラフとか予測と実測の散布図までできちゃう。素晴らしい。.explain()を使えばモデルの内容が出てくる。で、GEEでの結果はこういう感じ。avg_radというのは夜間光。夜間光、NDWI、MODISのバンド5が大事といっています。
Rの結果
randomForest()を使用。
似ているけど微妙に違う
ざっと見た感じ似てます。が、よく見ると重要な変数とかが違っている。両者で使っているデータが同じではない(GEE、Rそれぞれで訓練データ・検証データに分割しちゃったから、失敗)のと、まあランダムなのでということなのかな?RMSEはGEEでは983.6、Rで996.9とこれまた微妙な違い。大体似ているので、いいことにしたいと思います。
地図
それなりに、それっぽい地図ができました。適当に作ったにしては思ったよりもそれっぽくできている。
富士山の上のほう(関東地方外だけど)で人口が多めに出ているのは植物は生えていないか雪で色が他と違うからだろう。あと、成田空港も人口密集になってるけど、夜明るいからなんだろうと思う。
もっとよくできるという点はいくつもあるけど、今回はこれで終わり;
- メッシュとラスタの位置が対応していない
- sampleRegions()でのデータ取得は点からバッファ500mの平均(按分)でやった方が妥当だと思う
- 予測に使う変数を整理。ほかにも探してみて追加
- 人口0の訓練データをもうちょっと増やしてもいいのかも。
今回のコードはこちら
https://code.earthengine.google.com/f9b1398016c10d332e34d44832407085
- 人口データはe-statからダウンロードしました 出典:政府統計の総合窓口(e-Stat)(https://www.e-stat.go.jp/)
- 大変参考になったチュートリアルはこちら(再掲)
GEEでランダムフォレスト1:できるらしいと分かって
できたらいいな、できるはずと思っていました
GoolgeEarthEngineでランダムフォレストの分類ができるのだから、回帰もできるんじゃないかとは思っていたんですがググってもやり方が出てこなくってあれー?と思っていました。線形回帰はできるっていうのは知っていたけど、データがややこしくなってくるとそれじゃあ足りない。Rで分類器を作ってその中をほじくってGEEで条件分けとかするしかないのかなと、思っていたんですが。
超親切なチュートリアル
久しぶりにGEEをいじってみるかと思ってランダムフォレスト回帰をググったところ、教材がありました↓ 下記リンクの「04.2 DigitalSoilMappinginEarthEngine_RandomForestRegression.docx」に方法が説明されています。
このチュートリアルは土壌中の炭素の量?マッピングという課題に対して、NDVIなどの衛星データやらDEMやらを使って答えるという内容になっています。ランダムフォレスト「回帰」にするための変更は、.setOutputMode('REGRESSION')を追記するだけと超シンプル。前に石川県の植生分類を比較したときのコードと今回のを比べてみても、ほぼ同じ。
// 分類の場合 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 classifier = ee.Classifier.smileRandomForest({ numberOfTrees: 500, // ntree variablesPerSplit: 10, // mtry minLeafPopulation: 5, // Rule of Thumb: 1 for classification, 5 for regression bagFraction: 0.632, maxNodes: null, seed: 0}) .setOutputMode('REGRESSION') .train({features: training, classProperty: 'population', inputProperties: bands });
チュートリアルを試したところGEEでもランダムフォレストはサクッと動いて、これは便利そうだと。興味のある方は是非どうぞ。
ちょうどいいネタが思いつかない
で、自分でもでもやってみたくなったわけ(といいつつほぼコピペ)ですが、面白そうなネタで使えそうなデータがない。コメの収量とか味とか、あるいはワインやコーヒーの味とかを気象データや土壌データ、NDVIなどから予測とかできると楽しそう!と思ったんですが、座標のくっついたそういうデータはネット上に見つからず。考えてみれば農家の経営にかかわる情報だから、そりゃあ勝手に公開はできないよなと反省しました。酒蔵座標×日本酒の辛さとかならできるかも、とかも思いましたが、趣味でやるにはデータ準備が面倒なのでお手軽な公開データに頼ることにしました。
とりあえず人口
というわけで、簡単に使える人口データを使ってやってみます。今回はe-statにある4次メッシュ(約500m)の人口データの関東地方分を使うことにして、GEE上にある適当なデータを使って予測してみることに。Rでも同じようなことをやってみて、結果を見てみます。
SARの勉強5:データを見てみる3 植林地、落葉樹林、草地、農地
森林の種類で違うかはわからない
何か所で見てみたけど、いまいちわからない。 というのは、斜面の方位とか傾斜の大きさでVVもVHも相当影響を受けていて、森林の種類の違いなのか分からなかった。 グラフでも画像でもいまいち。切り取り方がへたくそなのか。
こんな感じのコード: https://code.earthengine.google.com/673838da213a372c9d9e373d2cd49856
森林と草原は違う
これは違いそう。場所は信州の霧ヶ峰。森林>草原です。草丈とか密度によっても変わるのかもしれないけど。
画像で見てもわかる。
上のとほぼ同じコード: https://code.earthengine.google.com/0e957a3e58cf855195a9618f06965839
時系列として扱える:農地
前にMODISのNDVIでやったような、Harmonic Modelへの当てはめもできる。三角関数をフィットさせて位相と振幅(=山のピークの位置と高さ)を求めればよい。
春は何も植わってなくて、秋に収穫してるのかな?と推測できる。
これをうまく使えば、畑の作物マッピングにつかえるかもしれない。畑といえばなんとなく宮崎のイメージだったので、Google Earthでみて畑が多そうだったところを選んで分析。以前の先生のスクリプトとほぼ同じ感じで書き出してみる。
続いてそのあたりの筆ポリゴンをダウンロードして畑だけを抽出して、畑単位の.reduceRegions()もしてみる。建造物とかが隣にあると影響を受けそうだったから、平均値ではなく.median()を取った。 現地データがあれば、作物の分類に使えそう。面白い。
コードはこちら(処理にかなりの時間がかかる): https://code.earthengine.google.com/7be06b74e747af389a62cd85e47c9d9c
SARの勉強5:データを見てみる2 八郎潟
八郎潟のSAR画像が面白かった
前回見た八郎潟が面白かったので、詳しくみる。 何が面白かったって、あるタイミングで黒い部分が一気に増えたこと。これは田んぼに水を張ったからだと思う。
時系列グラフを見る
八郎潟の田んぼを3か所適当に選んでポリゴン→Featureにしてから、グラフを書く。 平らだからDESENDINGとAECSENDINGを一緒にしたら、ガタガタのグラフになった。DESとAECで値が違っているみたい。
ならば、DESCENDINGだけ選択してグラフを書き直し。
水を張る時期は分かる。稲刈りは不明
VVもVHも5月に急に0近い値に落ち込む時期がある。時期的にも、そして土とかよりも水は後方散乱係数が低いという性質からも、水を張った時期だなと推測できる。その後少しずつ値が上がるのは稲が育って水面が隠されてということだろう。
で、稲刈りは?よくわかんない。積雪の影響もよくわからない。
水張った田んぼ地図を作る
ちょっとコードを書き足したら簡単な地図が作れる。 田んぼのポリゴンが公開されているので、このポリゴンをもとにSAR画像を.reduceRegions()するだけ。こういった手間のかかる情報が公開されているのは感謝です。
大潟村のポリゴンのうち、「田」だけを選択して使った。
それっぽいのができる。
色塗りの閾値は適当に決めたけど、ちゃんとやったらちゃんとした地図になるはず。後は各回の画像を組み合わせて条件抽出すれば一枚になった水張り時期地図ができるはず。 データの取得がもっと細かいほうが面白そうだけど、Sentinel-1は12日だからしょうがない。ASCENDINGとDESCENDINGを組み合わせても、この地域は両者の通過時の差が1日だけだったからそんなに変わらない。何に使えるかというと、
まとめ
ちょうどいい撮影間隔の場所であれば、田んぼに水を入れる時期がいい感じでわかる(だろう)。
ちょうどいい撮影間隔の場所で、十分な現場データがあれば、稲の生長がわかるのかも。ASCとDESの差も何らかのスムージングをするか、校正する回帰式で整理できたらいいかも。
水を張っていない田んぼの検出はできそう。実際に耕作されている田んぼだけ抽出というのはできそう。ただ、山間の田んぼや一枚の面積が小さい田んぼだとしんどい場合もありそう。
コードはこちら
https://code.earthengine.google.com/390219941d904d9e56687f75fae1b036