【问题标题】:Is there a way to use fusion algorithms like STARFM in Google Earth Engine?有没有办法在谷歌地球引擎中使用像 STARFM 这样的融合算法?
【发布时间】:2020-11-06 11:48:39
【问题描述】:

Landsat 和 MODIS 产品各有优势。具有高空间分辨率的 Landsat 和具有高时间分辨率的 MODIS。我已经阅读了很多关于下载文件并将它们与 Python 中的 STARFM 等算法在本地融合的内容。有没有办法直接在 Google 地球引擎中融合这两个集合以节省计算时间?

【问题讨论】:

    标签: algorithm resolution google-earth-engine satellite-image landsat


    【解决方案1】:

    GitHub 中有 Google 地球引擎代码。链接如下:

    https://github.com/MacGallagher/GEE_STARFM_FUSION

    我复制的代码如下:

    /**** 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 megangallagher@u.boisestate.edu
    // 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)
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    

    【讨论】:

      猜你喜欢
      • 2013-08-02
      • 1970-01-01
      • 1970-01-01
      • 2014-02-11
      • 1970-01-01
      • 2020-05-10
      • 2019-01-05
      • 1970-01-01
      • 2011-09-08
      相关资源
      最近更新 更多