山梨チャレンジ３：調和解析（Harmonic model）への道6
目的はなんだっけ？
前回までで、このまま進めてもダメだなと思って気づいたんです。目的はなんだったっけ？そう、調和解析やってみるってことだった。そして、それは達成済みだった。
で、なんで今行き詰ってるのかなと思うと、山梨のデータをきれいにしたいという、手段が目的になりつつあって、そこんところがこんがらがってるのねと。まとめると、調和解析はできた、けど山梨のデータは普通にHarmonic modelをあてはめるのは難しい。ということでここからは目的をアップデートして、「山梨のデータをどうやったらきれいに使えるか？」にしようと思います。
新テーマ「山梨の植物季節をMODISから読み取る」
で思ったんです。何年かに一回くらいは晴れの日があるから、数年分のおなじ撮影日のデータで最大値なり、90％分位なりをとれば、雲のない画像ができるんじゃないかと。で、それを平滑化したりすれば行けそう。
データを見てみよう
横軸をDOY（１年〇日目）にして、過去６年分のNDVI値を表すとこうなる。
まずは扱いやすい勝沼のNDVI
問題の多い鳳凰三山のNDVI
勝沼はよし。鳳凰三山だけど、こうしてみると、やっぱり晴れている年はある。正月付近の高いNDVI値は謎だけど、うまいこと季節変化をくみ取れそう。
90分位点（90 percentile）でreduce
最大値だと上側に振れるエラーもありそうなんで、とりあえず90パーセンタイルで処理してみましょう。さっきの見本は過去6年分だったけど、90パーセンタイルにする切りのいいところで過去10年（2011-2020）にしてみる。
Medianフィルターをかける
鳳凰三山の９月とか、スパイクというほどではない凸凹が残っている。やんなくてもよさそうだけど、窓サイズ5のメディアンフィルタをかけてみる。
Harmonic modelのあてはめ
これだけ滑らかだったらHarmonic modelあてはめなくても値の計算で出してもいいと思うけど、あてはめます。なぜなら波の大きさとかピーク時期とかの計算が楽だから。するとこうなる↓
まとめ
感想
というわけで、こういう方法だったらフェノロジーが出せることがわかりました。ただ、数年分まとめたデータが必要なので「今年のデータがどうしてもみたい（農業とか）」「経年変化がみたい」なんかだとこれじゃあ無理です。雲が少ない場所だったら前回までの方法でできるかもしれないけど。あと、この方法である場所の雲がない（と思われる）データをまとめてベースラインとして、それとの比較で何か言える場合もあるかもしれない。
今回のコードはこちらに：https://code.earthengine.google.com/22094521a2dbc7fa9626e27b14f34d18
地図にするとおしゃれ
Nick先生の動画にもあるんですが、地図にするとおしゃれです。というか、これがやってみたかったんだよ。
ピーク日を色の環で、波の高さを色の濃さ（濃い：波が高い、薄い：波が低い）で表すと素敵な地図に。
地図の解釈
山梨では山の中腹で色が薄い帯がみえる。これは亜高山帯に分布する常緑の針葉樹林と多分一致してて、こういう林ではNDVIの季節変化が小さいんだと思う（樹海みたいに）。あとは都市部でも色は薄い（納得）。高山は色が黄色っぽい（=ピークが遅い）。まあ、そうなんだ、と思いました。
あと、関東近郊を見ると色が緑の埼玉県深谷市付近と、水色の三浦半島の南側が目立つ。GoogleEarth画像でみると畑のようなので、それぞれネギと大根のせいだと思う。どちらも冬の野菜だから、夏場は畑に何も生えてなくて秋から冬にかけて成長のピーク（=NDVIのピーク）がやってくるんだろうね。こういうことがわかるとは思っていなかった。
何に役立つのか？
思いつくことを書いてみる；
- 常緑樹林と落葉樹林の分類するときにバンドとして入れると分類精度とかが上がるといいなと思う
- たとえば10年ごととかに区切って比較すれば長期間で生物季節（展葉や落葉の時期、作物の栽培・収穫時期）の変化がわかるのかも（温暖化とか）
- 上の項目と重なるけど、今年は例年よりあれだねーとか、リアルタイム的に見るときの基準として使う
- 気候による生物地理的な区分（バイオリージョンみたいな）の元データに使える気がする。そういう縮尺で作っていないんだろうけど、もうちょっと細かく作ってあったらと思うシェープファイルとかあるから、その辺の改良？につかえるのかも。まあ、この辺は既にいろいろ使われていそう。
今回の学び
1週間ほどで多くを学びました。最初に書いたのを今見るとダメなので、今後みるともっとダメに見えてくるでしょう。忘れないうちにメモを。
よく使う関数 updateMask
データ抽出やエラー値の除去に欠かせないImage.updateMask()。updateMaskのカッコの中にマスクしたい場所を0/1で表した画像をいれると、1の部分だけ残って0の部分はマスクされる。マスクしたい場所を示すためには.eq()、lt()、gt()などを使えば画像が作れる。
var maskQA = function(image) { return image.updateMask(image.select("SummaryQA").eq(0)); // バンド"SummaryQA"の値が0の部分を残して、それ以外はマスクする };
よく使う関数 .bind()からの.map()
ImageCollectionを構成する各Imageに対して処理を行うとき、Google Earth Engineでは.map()を使うことが多い1。ただ、.mapする関数は引数を一つしか取れない（ImageCollection）ので、複数の値を与えなければいけない関数の時は困る。
その時は.bind()を使えばうまいこと行く。そのほかにも関数を入れ子にして作る方法がある（こっちが主流っぽい）けど、頭がこんがらがったので.bind()を使った。便利。
// 引数が２つある関数 var funcDiffImage = function(centreImage, eachImage){ //最後にくる引数がImageCollectionの各Image（.bind()ではnullのところにあたる） var eachDate = eachImage.date(); var dateDiff = centreImage.date().difference(eachDate, 'day').toFloat(); return eachImage.set('system:dayDiff', dateDiff); }; // 関数名.bind(null, 指定する引数) var add_funcDiffImage = funcDiffImage.bind(null, image); // ImageCollection.map(bindした関数) var addDayDiff = filteredImageCollection.map(add_funcDiffImage);
