GEEで分布推定モデル２：MaxEnt

SDM界の iris こと Bradypus

とりあえず動かせるかみるために、分布推定モデルでお題になるミツユビナマケモノ Bradypus のデータを使ってやってみる。お手本はRのライブラリ dismo を使ったSDM各種を解説したページから。

データの準備

ミツユビナマケモノの分布地の座標は元祖のページからcsvがダウンロードできる（ほかでもいろいろある）。それをGoogle Earth EngineのAssetsにアップロードする。

環境データはお手本にそっておなじみBioClimから８レイヤと生物区系（バイオーム）。生物区系のデータはGEEのデータカタログにあったEco Regionを使ったけど、これはお手本データと微妙に違うかもしれない。次の手順は環境データ間の相関をみるだけど省略。お手本データがそれをやっているはずということで。¹ これらをまとめて一つの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()を使って発生させる。発生させる範囲はお手本にそって調査範囲＋その周り²の陸域とした。在地点の周囲○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）。地図での出力はこんな感じ。

f:id:camarao:20220308182607p:plain — Google Earth Engineで出力したMaxEnt

だいたい似てる。

二値化してみる

これも色々説はあるみたいけど、お手本のとおり感度と特異度の和が最大になる点を閾値として2値化してみる。この値は.explain()で出てくるなかのThresholds -> Minimum test sensitivity plus specificityに入っている。で、これが入れ子になったDictionaryに入っていて取り出し方がわからず面倒だった。この閾値の値はお手本と結構違っていた（0.32 vs 0.17）。

f:id:camarao:20220308182614p:plain — 二値化して在不在にした図。だいたいあってる。

まとめ

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のバッファとした。調査範囲も目分量で合わせたからちょっとずれてる。↩