Is there a way to use fusion algorithms like STARFM in Google Earth Engine?

Question

Landsat and MODIS products both have their advantages. Landsat with its high spatial resolution and MODIS with its high temporal resolution. I've read a lot about downloading the files and fuse them locally with algorithms like STARFM in Python for example. Is there a way to fuse both collections directly in Google Earth Engine to save computational time?

Answer 1

There is Google earth engine code in GitHub. The link is below:

https://github.com/MacGallagher/GEE_STARFM_FUSION

The code I copied is below:

/**** Start of imports. If edited, may not auto-convert in the playground. ****/
var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR"),
    modis = ee.ImageCollection("MODIS/006/MOD09GQ"),
    geometry4 = /* color: #d63000 */ee.Geometry.Polygon(
        [[[-116.004638671875, 42.809506838324204],
          [-115.59539794921875, 42.78129125156276],
          [-115.499267578125, 43.004647127794435],
          [-115.499267578125, 43.28920196020127],
          [-116.00189208984375, 43.35114690203121],
          [-116.19415283203125, 43.07691312608711]]]),
    geometry = /* color: #98ff00 */ee.Geometry.Polygon(
        [[[-115.93391418457031, 43.09346209534859],
          [-115.80276489257812, 43.095718426590174],
          [-115.80276489257812, 43.161366298103566],
          [-115.9332275390625, 43.16086543618915]]]);
/***** End of imports. If edited, may not auto-convert in the playground. *****/
// based on script by Faye Peters, University of Louisville
// modified by Megan Gallagher, Boise State University
// for code inquiries, or adjustments please email [email protected]
// code exists for: 
//    landsat 5,7, and 8 NDVI merging with MODIS Terra daily
//    MODIS Terra and Sentinel-2
//    Landsat pre tier and modis daily terra - works with base R code
//    Landsat 8 and Sentinel-2
////////////////////////////////////////Variables//////////////////////////////////////////////
// import Landsat 8 and MODIS daily terra surface reflectance 250m
var region = geometry; //region you want to export 
var bounds = geometry4 //outer bounds of image to catch boundary effect, make larger than region
var date_begin = '2016-03-02' //start date of data collection, must be a landsat image date
var date_end ='2016-10-13' // end date of data collection, must be the day AFTER last landsat image

var csv_title = '2016_Dates_BOP' //title of csv output
var ls_title = '2016_landsat'   // title of landsat tif file
var mod_title = '2016_mod'      // title of modis tif file

/////////////////////////////////////////Code////////////////////////////////////////////////////////////////////////

var starfm = function(modis, l8){

//Preliminary filtering of MODIS and Landsat
var filt_mod = modis.filterDate(date_begin, date_end);
var filt_l8 = l8.filterDate(date_begin, date_end)
              .filterBounds(bounds);

//add bounds to subset area
var subset_bounds = bounds.transform('EPSG:4326', 30).bounds().getInfo();


// get rid of bad pixels but keep metadata date information
var removeBadObservations = function(image){
    var valid_data_mask = ee.Image(image).select('pixel_qa').bitwiseAnd(2).neq(0);
var numberBandsHaveData =
  image.mask().reduce(ee.Reducer.sum());
var allOrNoBandsHaveData =
 numberBandsHaveData.eq(0).or(numberBandsHaveData.gte(9));
 var allBandsHaveData = allOrNoBandsHaveData;

 //Make sure no band is just under zero
 var allBandsGT = image.reduce(ee.Reducer.min()).gt(-0.001);
 var result = ee.Image(image).mask(image.mask().and(valid_data_mask).and(allBandsHaveData).and(allBandsGT));
 return result.copyProperties(ee.Image(image),['system:time_start']);
};
//NDVI functions for Landsat and MODIS
var getNDVI_mod = function(image){
return image
    .addBands(image.normalizedDifference(['sur_refl_b02','sur_refl_b01']).multiply(10000).rename('NDVI'));}; 

var getNDVI_l8= function(image){
var ndvi = ee.Image(image).normalizedDifference(['B5','B4']);
return ndvi.copyProperties(ee.Image(image),['system:time_start']);
};


var filt2_l8= filt_l8.filterBounds(subset_bounds).aside(print).map(removeBadObservations).map(getNDVI_l8);
// update mask for different areas
var subset_mask = ee.Image().byte().paint(geometry4, "id").add(1);

var filtered_modis = filt_mod.filterBounds(subset_bounds).aside(print).map(getNDVI_mod);

// Pull out the date:
var extract_modis_date = function(row) {
var d = ee.Date(row.get('system:time_start'));
var d2 = ee.Date.fromYMD(d.get('year'), d.get('month'), d.get('day'));
var result = ee.Feature(null, {'date': d2});
result = result.set({'date': d2});
return result;
};

var getQABits = function(image, start, end, newName) {
 //// Compute the bits we need to extract.
var pattern = 0;
  for (var i = start; i <= end; i++) {
 pattern += Math.pow(2, i);
 }
 //// Return a single band image of the extracted QA bits, giving the band a new name.
return image.select([0], [newName])
 .bitwiseAnd(pattern)
 .rightShift(start); };
 
print(filtered_modis);

filtered_modis = filtered_modis.map(function(image){
var quality = getQABits(image.select(2), 4, 5, 'QAMask');
quality = quality.eq(3).not();
return image.clip(subset_bounds).mask(image.mask().multiply(subset_mask).multiply(quality));
 });

filtered_modis = filtered_modis.select('NDVI'); //Take only NDVI band
 
 var modis_multiband = ee.Image(filtered_modis.filterDate(date_begin, date_end).iterate( function(x, modis_multiband) {
 return ee.Image(modis_multiband).addBands(ee.Image(x));
 }, filtered_modis.first()));

var dates_modis = filtered_modis.map(extract_modis_date);
 print(dates_modis.getInfo());
 
var filt2_l8_ndvi = filt2_l8.map(function(image) {
return ee.Image(image)
 .clip(subset_bounds)
.mask(ee.Image(image).mask().multiply(subset_mask));
 });
 
//use this to choose the range of days, and check for overlap
var day_expand = 1;

var reduceLandsatNDVI = function(MODISdate) {
 MODISdate = ee.Date(MODISdate.get('date'));
 
var ndvi_subset = ee.ImageCollection(filt2_l8_ndvi).filterDate( MODISdate,
 MODISdate.advance(day_expand, 'day') );
 
 ndvi_subset = ndvi_subset.map(function (image) {
  var diff = MODISdate.difference(ee.Date(ee.Image(image).get('system:time_start')), 'day').abs();
  return ee.Image(image).set('diff', diff);
 });
 
ndvi_subset = ndvi_subset.sort('diff');
 var ndvi_first = ndvi_subset.reduce('first');
 var ndvi_mean = ndvi_subset.reduce('mean');

return ee.Algorithms.If(
 ndvi_first.bandNames(),
 ndvi_first.eq(0).multiply(ndvi_mean).add(ndvi_first),
 ee.Image(0)
 );
 };
 
 var extract_landsat_date = function(MODISdate) {
 MODISdate = ee.Date(MODISdate.get('date'));
 
var ndvi_subset =
 ee.ImageCollection(filt2_l8_ndvi).filterDate( MODISdate,
 MODISdate.advance(day_expand, 'day') );

ndvi_subset = ndvi_subset.map(function (image) {
 var diff = MODISdate.difference(ee.Date(ee.Image(image).get('system:time_start')), 'day').abs();
return ee.Image(image).set('diff', diff);
 });

 ndvi_subset = ndvi_subset.sort('diff');
 var d = ndvi_subset.aggregate_first('system:time_start');
 var count = ndvi_subset.aggregate_count('system:time_start');
 
 d = ee.Algorithms.If(
 ee.Number(count).gt(0),
 ee.Date(d),
 ee.Date('1971-01-01')
 );
 d = ee.Date(d);

 var d2 = ee.Date.fromYMD(d.get('year'), d.get('month'),
 d.get('day'));

 var result = ee.Feature(null, {'LSdate': d2, 'MODISdate':
 MODISdate, 'CountLSScenes': count});
 result = result.set({'LSdate': d2, 'MODISdate': MODISdate,
 'CountLSScenes': count});
 return result;
 };

var ls_collection = dates_modis.map(reduceLandsatNDVI);
print(ls_collection,'ls collection');


var dates_landsat = dates_modis.map(extract_landsat_date);


Export.table(dates_landsat, csv_title);

var ls_multiband = ls_collection.iterate( function(x,
 ls_multiband) 
{ return ee.Image(ls_multiband).addBands(ee.Image(x));
}, ls_collection.first());

ls_multiband = ee.Image(ls_multiband).multiply(10000).int16(); 
ls_multiband = ls_multiband.mask(ls_multiband.mask().multiply(ls_multiband.neq(0))); 


//remove repeated first layer
var ls_multiband2=ee.Image(ls_multiband.slice(1))
var modis_multiband2=ee.Image(modis_multiband.slice(1))


Export.image.toDrive({
  image: modis_multiband2,
  description: mod_title,
  crs:'EPSG:4326 ',
 region:region,
 scale:30
});


Export.image.toDrive({
  image: ls_multiband2, 
  description: ls_title, 
  crs:'EPSG:4326',
  region:region,
  scale:30
});

return 'Done!'
}

var running = starfm(modis,l8)
print(running)
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////