MODIS前処理：調和解析（Harmonic model）への道３

ほかに方法はないものか

前回のチャレンジでは空白域ができちゃったもんで、なんか別の方法はないかなーと考えてみる。まあ、簡単なのは、QAのデータの判定を甘くする。Good dataに加えてMarginal dataをいれちゃえばいい。でも、「useful but look at detailed QA for more information」の部分がやや面倒。もう一つはある種の数値のフィルタをかけちゃう方法。具体的なやり方はいろいろありそうだけど、まあ、考え方は似てるはず。っつうわけで、数値を使ったフィルタリング方法を考えてみる。

最も単純な移動平均

基本は移動平均だよね、ということでこれは簡単に書ける。元ファイルに移動平均のバンドが追加されるはず。

var filteredModis = terra
                       .merge(aqua)
                       .filterDate('2003-01-01','2020-12-31')
                       .sort("system:time_start")
                       .select('NDVI');

var movingWindowDays = 2*8; // MODISが8日間隔なので、ここでは片方2回両側で5回の平均

var movingAverage = function(image){
  var dayEnd = image.date().advance(movingWindowDays, 'day') ;
  var dayStart = image.date().advance(-movingWindowDays, 'day') ;
  var selectedImageCollection = filteredModis.filterDate(dayStart, dayEnd);
  
  return image.addBands(selectedImageCollection.mean().rename('average')); 
};

var aveModis = filteredModis.map(movingAverage);

グラフはこんな感じ。場所は富士樹海だけど、ピンを落としなおしたから前回とは微妙に違うところです。やっぱデータが暴れてる。

5回の移動平均

TIMESATを参考にする

前に触ったことのある、リモセン時系列パッケージTIMESATっぽいことできないかな、と論文とウェブにあるマニュアルを読んでみる。まあ、いろんなテクが載っているわけですが、簡単にできそうなところから真似てみる。

TIMESAT流 外れ値を除去する前処理

マニュアルのほうの3.3に外れ値の前処理についての記述がある。３つの方法がある中で、簡単（でデフォルト）な処理方法をやってみる¹。これは標準偏差と前後の値を基準として、それを逸脱するデータを外れ値として扱うもの。具体的には

• （対象とする値の前後でとる）Moving window（日本語で移動窓？でいい？）の中央値を基準とし、その値±カットオフ値を逸脱するもの、
• 対象とする値の前後両隣の平均－カットオフ値よりも小さい、または
• 対象とする値の前後両隣の最大値＋カットオフ値よりも大きい
• ここで「カットオフ値」は全時系列データの標準偏差を基準として、ユーザーが〇倍と決めて設定

の手順で、外れ値除去をします。マニュアルに何も書いてなかったのでカットオフ値の設定に使う標準偏差は超単純に計算しましてみましたが、1回の階差の標準偏差としたほうがいい気も。

// Moving window の大きさを決める
var preWindowSize = 2;  // for one side
var preMovingWindowDays = preWindowSize*8+1;

// cutoff値を決める
var sdNdvi = filteredModis.select('NDVI')
                    .reduce(ee.Reducer.stdDev()); //単純な方法

var cutoff = 2; // とりあえず標準偏差の2倍
var cutoffImage = sdNdvi.multiply(ee.Image.constant(cutoff));

//メインの関数
var preMask = function(image){
  var imageDate  = image.date();
  var dayEnd   = imageDate.advance(preMovingWindowDays, 'day') ;
  var dayStart = imageDate.advance(-preMovingWindowDays, 'day') ;

  // 移動窓の範囲の画像を選択
  var filteredImageCollection = filteredModis.filterDate(dayStart, dayEnd)
                                       .select('NDVI'); 
  
  // 中央値を算出して、±SDの画像をつくって比較演算でマスクの条件を決める
  var medNDVI = filteredImageCollection.median();
  var medAddSd = medNDVI.add(cutoffImage);
  var medSubSd = medNDVI.subtract(cutoffImage);
  var medMask = image.select('NDVI').gt(medSubSd).and(image.select('NDVI').lt(medAddSd));
  
  // 日数の差を数えて、dayDiffというメタデータをつくって書き込む関数
  var funcDiffImage = function(eachImage, centreImage){
    var eachDate = eachImage.date();
    var dateDiff = centreImage.date().difference(eachDate,'day').toFloat();
    return centreImage.set('system:dayDiff', dateDiff);
  };
  
  var add_funcDiffImage = funcDiffImage.bind(null, image);
  var addDayDiff = filteredImageCollection.map(add_funcDiffImage);
  
  // ひとつ前と後ろの画像を選ぶ
  var imgBefore = addDayDiff.filterMetadata('system:dayDiff',"less_than", 0)
                            .sort('system:time_start', false)
                            .first();
  var imgAhead = addDayDiff.filterMetadata('system:dayDiff',"greater_than", 0)
                           .sort('system:time_start')
                           .first();

  // カットオフ値を足したり引いたりして、比較からマスクの条件をつくる
  var meanBeforeAhead =  ee.ImageCollection.fromImages([imgBefore, imgAhead]).select('NDVI').mean().subtract(cutoffImage);
  var maxBeforeAhead  =  ee.ImageCollection.fromImages([imgBefore, imgAhead]).select('NDVI').max().add(cutoffImage);
  var sdMask = image.select('NDVI').gt(meanBeforeAhead).and(image.select('NDVI').lt(maxBeforeAhead));

  //マスクする
  var medFiltered = image.updateMask(medMask)
                         .select('NDVI');

  var sdFiltered = image.updateMask(sdMask)
                        .select('NDVI');

  var bothFiltered = image.updateMask(sdMask.and(medMask))
                        .select('NDVI');


  return image.addBands(ee.Image(sdFiltered).rename('sd'))
              .addBands(ee.Image(medFiltered).rename('median')
              .addBands(ee.Image(bothFiltered).rename('both')));
};

var preProcessed = filteredModis.map(preMask);

長くなったけど、短く上手な人は短く書くんだろうなきっと。例によって富士樹海ですが、明らかにダメそうなデータは除去されている。

ノイズ除去処理後

Savitzky-Golay フィルタ²

えー、これはかなり苦闘しました。大作になったので入りきりません。アカウント持っていたらリンクから見られるはず。移動窓の数値に二次関数をあてはめるってことで、なだらかになってしまいがちな移動平均よりピークを捉えやすいとのこと。まあ、富士樹海じゃああまり威力発揮してませんね。

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

SGフィルタ

結局どれがいいのか

で、わかった、てか知ってた。結局万能な方法はなくデータと相談しながら妥当な方法を探っていく。

ノイズらしきデータをどんどん削っていけばいいのかもしれないけど、そうすると歩留まり悪くなるしね、ピークやら底やらを見逃すかもしれない。移動平均やSGフィルタは欠測が多くなりすぎると困る³。移動窓を大きくとるとなだらかになるけど、ピークがはっきりしなくなる。

というわけで、方法を組み合わせたりしてみるのがよさそう。あとは、狭い範囲であれば矛盾少なく対処できそうだけど、広範囲になればなるほどひずみが出そう。

---

1. ちなみに後の二つ簡単でない方法は要STL decompositionでちょっと込み入りそうなので、取り組みません。↩

2. TIMESATのはピーク周辺で移動窓サイズを狭めたり、雲の影響を抑えるため全体的に上に持ち上げるなど高等テクが付け加わった「Adaptive Savitzky-Golay filtering」で、ずっと複雑です（俺には再現できなそう）。ここでは最も単純な移動窓サイズ固定の二次式あてはめだけをやってます。↩

3. 特にSGフィルタは二次式をあてはめるので、欠測で変な値が残ってるととんでもない値が出現することあり。↩