MODIS_TFA_Daily_Mean_Export_2_Assets

// This code uses Temporal Fourier Analysis to build 365 images, for each day of year, 
// of the climatological tempratures based on the daily mean temperatures.

// Change the path to your own assets folder
var Assets_path = 'users/shilosh/MODIS_TFA_Daily_'

var Temperature_Band = 'LST_Daily_Mean';
var Day_Temperature_Band = 'LST_Day_1km';
var Night_Temperature_Band = 'LST_Night_1km';
var collection = 'MODIS/006/MYD11A1';
// var geometry = ee.Geometry.Rectangle([33.2,29.0,36.6,34.0]); //Israel
var geometry = ee.Geometry.Rectangle([20,50,25,54]); //Poland
// var geometry = ee.Geometry.Rectangle([43,0,48,4]); //France
// var geometry = ee.Geometry.Rectangle([49,7,54,11]); //Germany
// var geometry = ee.Geometry.Rectangle([-76,36.0,-79,40]); //Maryland
// var geometry = ee.Geometry.Rectangle([-90,41,-95,44.0]);//Iowa
// var geometry = ee.Geometry.Rectangle([-60,-6,-55,-2]); //Brazil
// var geometry = ee.Geometry.Rectangle([-50,-20,-45,-16]); //SouthBrazil
// var geometry = ee.Geometry.Rectangle([-69,-34,-64,-30]); //Argentina
// var geometry = ee.Geometry.Rectangle([12,24,17,28]); //Libya (Sahara)
// var geometry = ee.Geometry.Rectangle([105,30,110,34]); //China
// var geometry = ee.Geometry.Rectangle([145,-37,150,-33]); //SouthAustralia
// var geometry = ee.Geometry.Rectangle([135,-25,140,-20]); //MidAustralia
// var geometry = ee.Geometry.Rectangle([-130,58,-125,62]); //Canada (Mid latitudes)


// var TFA_File_Name = 'MODIS_Daily_TFA_Israel'
var TFA_File_Name = 'MODIS_Daily_TFA_Poland'
// var TFA_File_Name = 'MODIS_Daily_TFA_France'
// var TFA_File_Name = 'MODIS_Daily_TFA_Germany'
// var TFA_File_Name = 'MODIS_Daily_TFA_Maryland'
// var TFA_File_Name = 'MODIS_Daily_TFA_Iowa'
// var TFA_File_Name = 'MODIS_Daily_TFA_Brazil'
// var TFA_File_Name = 'MODIS_Daily_TFA_SouthBrazil'
// var TFA_File_Name = 'MODIS_Daily_TFA_Argentina'
// var TFA_File_Name = 'MODIS_Daily_TFA_Libya'
// var TFA_File_Name = 'MODIS_Daily_TFA_China'
// var TFA_File_Name = 'MODIS_Daily_TFA_SouthAustralia'
// var TFA_File_Name = 'MODIS_Daily_TFA_MidAustralia'
// var TFA_File_Name = 'MODIS_Daily_TFA_Canada'

var geometry_json = geometry.toGeoJSON();
// The number of cycles per year to model.
var harmonics = 3;

// Make a list of harmonic frequencies to model.  
// These also serve as band name suffixes.
var harmonicFrequencies = ee.List.sequence(1, harmonics);
// Function to get a sequence of band names for harmonic terms.
var getNames = function(base, list) {
  return ee.List(list).map(function(i) { 
    return ee.String(base).cat(ee.Number(i).int());
  });
};

// Construct lists of names for the harmonic terms.
var cosNames = getNames('cos_', harmonicFrequencies);
var sinNames = getNames('sin_', harmonicFrequencies);

//convert Kelvin to Celsius
var k2celsius = function(image) {
      return image.updateMask(image.select(Day_Temperature_Band))
                  .updateMask(image.select(Night_Temperature_Band))
                  .reduce( ee.Reducer.mean()).rename(Temperature_Band)
                  .multiply(ee.Image(0.02))
                  .subtract(ee.Image(273.15))   //convert Kelvin to Celsius
                  .set('doy', image.get('doy'));
};
// Add a band with doy to the colection.
function createDoyBand(img) {
  var d = ee.Number(img.get('doy'))
  return img.addBands(ee.Image(d).rename('t'));
}

// Function to compute the specified number of harmonics
// and add them as bands. 
    // for p = 0, Nharmonics-1 do begin
    //   a[p] = a[p] + 2*TS_t*cos(2*!pi*(p+1)*(t+1)/N)
    //   b[p] = b[p] + 2*TS_t*sin(2*!pi*(p+1)*(t+1)/N)
var addHarmonics = function(freqs) {
  return function(image) {
    // Make an image of frequencies.
    var frequencies = ee.Image.constant(freqs);
    // This band should represent lst.
    var lst = ee.Image(image).select(Temperature_Band);
    // This band should represent days from start.
    var time = ee.Image(image).select('t');
    // Get the cosine terms.
    var cosines = lst.multiply(2.0).multiply((frequencies.multiply(2.0 * Math.PI).multiply(time).divide(365.0)).cos())
      .rename(cosNames);
    // Get the sin terms.
    var sines = lst.multiply(2.0).multiply((frequencies.multiply(2.0 * Math.PI).multiply(time).divide(365.0)).sin())
      .rename(sinNames);
    return image.addBands(cosines).addBands(sines);
  };
};

var LST_Day = ee.ImageCollection(ee.List.sequence(1, 365).map(function (doy){
  return  ee.ImageCollection(collection)
                        .select(Day_Temperature_Band)
                        .filter(ee.Filter.calendarRange(doy, doy, 'day_of_year'))
                        .mean()
                        .set({doy: doy})
                        //.addBands(ee.Image(ee.Number(doy)).rename('t'))
}))
var LST_Night = ee.ImageCollection(ee.List.sequence(1, 365).map(function (doy){
  return  ee.ImageCollection(collection)
                        .select(Night_Temperature_Band)
                        .filter(ee.Filter.calendarRange(doy, doy, 'day_of_year'))
                        .mean()
                        .set({doy: doy})
                        //.addBands(ee.Image(ee.Number(doy)).rename('t'))
}))
print('LST_Day = ', LST_Day)
// // For each day cloudy pixels with no data, insert yearly mean values.
// var LST_Mean
// function BlendWithYearlyMean(img) {
//   return img.blend(LST_Mean);
// }
// LST_Mean = LST_Day.mean()
// LST_Day = LST_Day.map(BlendWithYearlyMean)
// LST_Mean = LST_Night.mean()
// LST_Night = LST_Night.map(BlendWithYearlyMean)

// Use an equals filter to specify how the collections match.
var Filter = ee.Filter.equals({
  leftField: 'doy',
  rightField: 'doy'
});
  // Define the join.
  var innerJoin = ee.Join.inner('primary', 'secondary');

// Join LST_Day with LST_Night by DOY
  // Apply the join.
  var CFSV2_JoinInner = innerJoin.apply(LST_Day, LST_Night, Filter);
  
// Add band
var LST = CFSV2_JoinInner.map(function(f) {
  var Night = ee.Image(f.get('secondary'));
  var Day = ee.Image(f.get('primary'));
  
  return Day.addBands(Night);
})

LST = ee.ImageCollection(LST)
            .map(k2celsius)
            .map(createDoyBand)
print('LST = ', LST)
LST = LST.map(addHarmonics(harmonicFrequencies));

  // a = a / N_TS
  // b = b / N_TS
var num_days = 365;//ee.Date(lastDay).difference(ee.Date(firstDay), 'day');
var a_coef_1 = LST.select('cos_1').sum().divide(num_days);
var a_coef_2 = LST.select('cos_2').sum().divide(num_days);
var a_coef_3 = LST.select('cos_3').sum().divide(num_days);
var b_coef_1 = LST.select('sin_1').sum().divide(num_days);
var b_coef_2 = LST.select('sin_2').sum().divide(num_days);
var b_coef_3 = LST.select('sin_3').sum().divide(num_days);

  // ; ========================================
  // ; prepare the first 3 harmonics of the TFA
  // ; ========================================
  // H0        = total(TS) / N_TS                                ;mean
  // FTA       = MAKE_ARRAY(N_TS, /FLOAT, VALUE=0.)
  // H         = MAKE_ARRAY(Nharmonics, N_TS, /FLOAT, VALUE=0.)
  // omega     = MAKE_ARRAY(Nharmonics, /FLOAT, VALUE=0.)
  // Phase     = MAKE_ARRAY(Nharmonics, /FLOAT, VALUE=0.)
  // Amplitude = MAKE_ARRAY(Nharmonics, /FLOAT, VALUE=0.)
  // for p = 0, Nharmonics-1 do begin
  //   omega[p]     = 2. * !pi * (p+1) / N
  //   Phase[p]     = atan( -b[p], a[p] )
  //   Amplitude[p] = sqrt( a[p]^2 + b[p]^2 )

var H0 = LST.select(Temperature_Band).mean();
    // Get the omega terms.
var omega_1 = ee.Image(1.0).multiply(2.0 * Math.PI).divide(365.0);
var omega_2 = ee.Image(2.0).multiply(2.0 * Math.PI).divide(365.0);
var omega_3 = ee.Image(3.0).multiply(2.0 * Math.PI).divide(365.0);

var phase_1 = a_coef_1.atan2(b_coef_1.multiply(ee.Image(-1.0)));
var phase_2 = a_coef_2.atan2(b_coef_2.multiply(ee.Image(-1.0)));
var phase_3 = a_coef_3.atan2(b_coef_3.multiply(ee.Image(-1.0)));

var amplitude_1 = (a_coef_1.pow(2).add(b_coef_1.pow(2))).sqrt();
var amplitude_2 = (a_coef_2.pow(2).add(b_coef_2.pow(2))).sqrt();
var amplitude_3 = (a_coef_3.pow(2).add(b_coef_3.pow(2))).sqrt();

  //   for t = 0, N_TS-1 do begin
  //     H[p,t] = Amplitude[p] * cos(omega[p]*(t+1) + Phase[p])    
  
// Function to add a H band.
var addH = function(image) {
  var H_1 = amplitude_1.multiply((omega_1.multiply(image.select('t')).add(phase_1)).cos()).rename('H_1');
  var H_2 = amplitude_2.multiply((omega_2.multiply(image.select('t')).add(phase_2)).cos()).rename('H_2');
  var H_3 = amplitude_3.multiply((omega_3.multiply(image.select('t')).add(phase_3)).cos()).rename('H_3');
  // FTA[t] = H0 + H[0,t] + H[1,t] + H[2,t]
  var TFA = H0.add(H_1).add(H_2).add(H_3).rename('TFA');

  var Anomalies = TFA.subtract(image.select(Temperature_Band)).rename('Anomalies');
  return image.addBands(TFA).addBands(Anomalies);
};

var LST_TFA = LST.map(addH);

// Define the visualization parameters.
var vizParams = {
  bands: ['mean', 'phase', 'amplitude'],
  min:   [285.0,2.0,5.0],
  max:   [295.0,3.0,10.0],
  gamma: [1.0, 1.0, 1.0]
};

//display the image.
var img = ee.Image(H0.addBands(phase_1).addBands(amplitude_1).rename('mean', 'phase', 'amplitude'));
// Map.addLayer(img, vizParams, 'RGB');
            
// var modelTFA = LST_TFA.select(['TFA', Temperature_Band, 'Anomalies']);
// var modisTFA = LST_TFA.select(['TFA', Temperature_Band]);

// print(LST_TFA)
//Iterating over the image collection using this function....
// var geometry = ee.FeatureCollection('ft:19Un8ogYFYTnR9p_VsMQ1tlzbR-bVwZoj7ipKy-bC');
var TFA_Images = LST_TFA.select('TFA')
  .sort('t')
  .iterate(function(img, all) {
      return ee.Image(all).addBands(img);
  }, ee.Image().select());

// print('TFA_Images = ', TFA_Images );
// Map.addLayer(ee.Image(TFA_Images))

// Export one single file with all the bands
Export.image.toAsset({
  image: ee.Image(TFA_Images),
  assetId: Assets_path + TFA_File_Name,
  region: geometry_json,
  crs: ee.Image('MODIS/006/MYD11A1/2018_01_01').projection().crs().getInfo(),
  crsTransform: [926.625433056,0,-20015109.354,0,-926.625433055,10007554.677],
  description: 'exportToAsset-MODIS-TFA-' + TFA_File_Name,
  pyramidingPolicy: {'.default': 'sample'}
});