山梨チャレンジ４：フェノロジーの指標を出してみる

NDVIの波から直接いけるんじゃない？

前回悟ったように、そしてまとめたように数年分のデータをまとめれば波の形が見えてくる。Harmonicモデルをあてはめてピークの高さや時期を抽出したけど、波がみえてるんだったら、そこから直接取ればいいんじゃないと。

目標

TIMESATのトップページには9つの指標が載っている。季節のはじめ・終わり・中間・長さ。最大値、低い値、波の高さ、それとNDVIの積分値2種類。この辺を再現したい。

で、日にちについては現在８日毎データで、これを使ってもあまり差が出なそうだから、毎日のデータを作ってやっている。

まずは8日間隔→毎日に補間

こういうことをやるのが妥当か、やるとしたらどういう方法がいいのかという話は置いといて¹、できる方法を考えます。

局所回帰やSavitzky-Golayフィルタで平滑化するときに毎日データに対して回帰する。
局所回帰したアウトプット（8日毎）を等間隔で補間して毎日にする。

で、両方やってみたんですが、一つ目はかなり時間がかかったので、二つ目がよさそうです。ただ、単純な局所回帰だといいけど、SGフィルタで平滑化したものに使うのは如何なものかと。

// 毎日データへと穴埋めする関数
var fillingDay = function(imageCollection, image){
  var nextImage = imageCollection.filterDate(image.date().advance(1,'day'),image.date().advance(10,'day')).first(); //8日後にデータないと困るから（そんなことはないはず）、とりあえず10日分
  var changePerDay = nextImage.select('fitted').subtract(image.select('fitted')).divide(8).rename('fitted');　// 平均した1日当たりの変化量

  var dayN = function(image, day){
    return image.add(changePerDay.multiply(ee.Image.constant(day))) // 穴埋めNDVIデータ
                .set('system:time_start', image.date().advance(day,'day').millis())　// メタデータを与える
                .set('system:time_end', image.date().advance(day,'day').millis());
  }; 
  
 var days = ee.List([1,2,3,4,5,6,7]);
 var output = ee.ImageCollection(days.map(dayN.bind(null, image)));
 return ee.ImageCollection([image]).merge(output);
};

ここまでのリンク：https://code.earthengine.google.com/74b5b053b0f7551e960ed0d8dcf80c01

指標を取り出す

季節のはじめ・終わり・中間・長さ

季節のタイミングは一連のものなので、似たようにして取り出せるはず。やり方が簡単に思いついた最大・最小値²とその時期を求めて、そこから進めていく。平滑化したImageCollectionのImageに日付のバンドを付け加えて、こうする↓

// 最大値
var maximum = phenometricsIC.reduce(ee.Reducer.max(2)) // reducerは一つ目のバンドにかけて、それに対応するその他のバンドも値を出力する
                        .rename('fitted','date');
// 最小値（左側：時間の早い側）
var leftBase = phenometricsIC.map(function(image){return image.updateMask(image.select('date').lt(maximum.select('date')))}) //最大値の日より早いデータ限定
                             .reduce(ee.Reducer.min(2))
                             .rename('fitted','date');

TIMESATではシーズンのはじめ・終わりを「山の高さ＊〇倍」としていて、〇倍の部分は各自で入力することになっている。とりあえず20%(0.2倍)としてみる。

var threshold = function(baseImage, maxImage, fraction){
  var difference = maxImage.select('fitted').subtract(baseImage.select('fitted'));
  return difference.multiply(fraction).add(baseImage.select('fitted')); 
} ;

//　基準となるNDVIの計算
var leftThreshold = threshold(leftBase, maximum, 0.20).rename('threshold');　// 0.2倍

var startOfSeason = phenometricsIC.map(function(image){
                                            var thresholdMask = image.select('fitted').gte(leftThreshold);　//NDVI基準値より大きい
                                            var leftBaseDateMask = image.select('date').gt(leftBase.select('date'));　//左側の最小値よりも遅い日
                                            var maxDateMask = image.select('date').lt(maximum.select('date'));　//最大値の日よりも早い日
                                            return image.updateMask(thresholdMask.and(leftBaseDateMask.and(maxDateMask)))})
                                      .filterDate(startGrowthYear, endGrowthYear) 
                                      .select('date')
                                      .min(); //フィルタリングされた中で一番早い日

同じような方法で季節の終わりも求められる。そうすると単純なラスタのバンド間演算で季節の中間・長さも求められます。