山梨チャレンジ2：調和解析（Harmonic model）への道５

どうして、うまくいかないのか

問題のある前回の続きですが、この雲でダメなデータを何とかする必要があるわけです。考えた対処方針は３つ。

• 時間を延ばす：重くなるから5年分でやったけど、もっと長期間にすれば晴れてる日がたくさんあるかも。
• 前処理方法を変えてみる：ちょっと書いたけど、全体の標準偏差とか雑じゃない？勝沼みたいに二山型の分布になるだろうところもあるし。差分でやったらどうか？
• 何らかの閾値を設ける：低いNDVIをバッサリカットする。

上から順に穏当な方法でしょうか。最後のはもう、わがままというかなんというか、のような気がします。では、やってみましょう。前処理の設定は、移動窓サイズ：５（前後各２回）、閾値：２×標準偏差です。

時間を延ばす

filterDateの日付を変えただけ。10年分に伸ばしてみたけど、だめでした。結局毎年同じ時期に雲の画像があるみたい。

前処理後（10年に伸ばす）

Harmonic model（10年に伸ばす）

気持ちピークが早くなっているけど、これが正しくないという自信あり。

前処理方法を変えてみる

差分のSDのほうがいいんじゃないか

外れ値の検出に全体をザザっと計算して出した標準偏差をつかっているんだけど、これでいいのだろうかと不安でした。もしかしたら、論文を読み落としているかもしれなくてそれも不安¹。

まあ、いいんです、樹海のように年間を通して変化の少ない地点では↓

樹海のNDVIのヒストグラム

赤線は正規分布の曲線で、あってはいないけど、まあね²。 ただ、これが勝沼になるとこうなる。

勝沼のNDVIのヒストグラム

というわけで、二山型になる地点では標準偏差が大きくなるから、一山の地点より外れ値が甘くなる可能性が大。

で、差分を取るとこうなる。

樹海の差分のヒストグラム

勝沼の差分のヒストグラム

少なくとも勝沼は元よりはいい気がするけど、それほど当てはまりはよくない。考えてみると、季節差分も考慮すると…と（不可能ではなさそうだけど）大変険しい道に続いていることに気が付きました。なので、これがいいとは思えないけど、一旦打ち止めにして、差分でやってみる。

// 次の回に取得された画像との差分をとる関数
var addNextNdvi = function(imageCollection, image){
  var date = image.date();
  var nextNDVI = imageCollection.filterDate(date.advance(1, 'day'), date.advance(40, 'day')).first().select('NDVI'); //対象の画像の撮影日の翌日から40日後の間に撮影された画像を抽出し、その中で撮影日が早いものを選ぶ
  var differenceNDVI = image.select('NDVI').subtract(nextNDVI) ;　//対象の画像と上の行で選ばれた画像の差分
  return image.addBands(nextNDVI.rename('nextNDVI')).addBands(differenceNDVI.rename('differenceNDVI'));　//バンド「differenceNDVI」を付け加える
};

// ImageCollectionの最後の方のエラーを避けるため
var shortenedModis = filteredNDVI.filterDate(filteredModis.first().date(),filteredModis.sort("system:time_start",false).first().date().advance(-40, 'day'));
var nextAddedModis = shortenedModis.map(addNextNdvi.bind(null, filteredModis));

// 準備
var sdDiffNdvi = nextAddedModis.select('differenceNDVI') .reduce(ee.Reducer.stdDev()); //差分の標準偏差
var sdDiffCutoffImage = sdDiffNdvi.multiply(ee.Image.constant(cutoff)); // 標準偏差のcutoff（ここでは２）倍
var bindDiffFunc = preMask.bind(null, sdDiffCutoffImage);　// preMaskは以前作った関数

// 以下、以前のもの参照。

結局だめそう

TIMESATの前処理よりはいいのかもしれないけど（ひいき目）、だめですね。結局データに問題があるところは標準偏差も大きくなるので、フィルタリングが難しい場所が出てきてしまうのかも。

差分での前処理後

差分で前処理したHarmonic model

何らかの閾値を設ける

最後の手段です。強引にハードルを設けます。「全体のNDVI中央値－２SD」は無効！と決めます。平均でなく中央値としたのは、平均だと雲画像（NDVIが低い）に引きずられるから。2SDは根拠なしです。こういうスクリプトを付け加える↓

var cutoffTotal = filteredModis.select('NDVI').reduce(ee.Reducer.stdDev()).multiply(ee.Image.constant(cutoff));
var totalMask = function(image){return image.updateMask(image.gt(filteredModis.select('NDVI').median().subtract(cutoffTotal)))};

var totalMaskedNDVI = filteredModis.map(totalMask);

これと差分の前処理を組み合わせるとなんとか。

前処理後（閾値＋差分）

Harmonic model（閾値＋差分）

まとめ

一応できたけど、ちょっと無理がある。このデータに合わせて作っているから、ほかの地域への応用性というか汎用性も低そうだし。

このまま突き進んでもいいことなさそう、と、作業してるうちに別のアプローチを思いついたので、次はそちらへ。

---

1. TIMESATの元のプログラムを見たらいいのかもしれない、そうは思いますがやり方がわからないし、手一杯なところに加えてFORTRANを勉強するのも…↩

2. このノートで統計的な吟味は基本的にスルー。ちゃんとやるときはちゃんとやらないとダメですが、ここの目的はGoogle Earth Engineでのスクリプトの書き方の勉強なので。↩