Extract landsat pixel values over plots with measurement dates

Question

I am trying to get landsat reflectance values corresponding to field plots spread all over Finland in Google Earth Engine. Each plot is small, hence is represented by a single point (lat, lon). It also has a measurement date associated with it (eg: "2015-08-05", "2017-07-19"). For each such plot, I want to choose landsat images over it that were acquired within ± 6 months of the plot measurement date. For example, if the plot measurement date is 2015-08-15, only landsat images (over the plot) with an acquisition date between 2015-05-15 and 2015-11-15 (both dates inclusive) should be considered. And the pixel values over those plots should be returned (along with the plot ID, Landsat image ID, etc).

I have written some code for this (see below). Here 'myplots' is a featureCollection of features (points) each with two properties: plot_id (string) and m_date (string, measurement date, see above). My problem is that the code runs quite slow. That is, when I have ~60K plots spread over six years, it take more than an hour of processing time. Can you please suggest optimizations for this code? I feel that there is a better way than adding a new property to each plot feature.

My code:

var plot_list = [
  ee.Feature(ee.Geometry.Point(25.26,67.03), {plot_id: 'plot_1', m_date: '2015-06-14'}),
  ee.Feature(ee.Geometry.Point(26.29,65.07), {plot_id: 'plot_2', m_date: '2015-08-28'}),
  ee.Feature(ee.Geometry.Point(28.39,62.55), {plot_id: 'plot_3', m_date: '2016-07-01'})
];

var myplots = ee.FeatureCollection(plot_list);

// For any given image, intersect it with the plot set.
// The intersection is based on both the location of the plot
// and the measurement date associated with the plot.
// That is, 'location' means that the plot should fall in the
// image footprint. And 'measurement date' means that the measurement
// date of the plot should be 'near' the aquisition time of the image.
// In this context, 'near' could mean ± 1 month, ± 1 year, etc.
var getValsFromImage = function(img, fc){
  var inift = ee.FeatureCollection(fc);
  // gets the values for the points in the current img
  var fc2 = img.reduceRegions(myplots, ee.Reducer.first(), 30);
  // Discards null elements
  fc2 = fc2.filterMetadata('pixel_qa', 'not_equals', null);
  var i_Date = ee.Date(img.get('system:time_start'));
  var fc3 = fc2.filter(ee.Filter.dateRangeContains('daterange', i_Date));
  return inift.merge(fc3);
}

// ########### BLOCK: SELECTION OF SUITABLE LANDSAT IMAGES (BEGIN) ###################
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

// First, get a shapefile of Finland (in roi)
var allCountries = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017")
//Choose country using GEE Feature Collection 
var roi = allCountries.filterMetadata('country_na', 'equals', 'Finland')
// Now, get all suitable landsat images
var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterBounds(roi)
                  .filterMetadata('CLOUD_COVER_LAND', 'less_than', 5)
                  .filterMetadata('IMAGE_QUALITY_OLI', 'equals', 9)
                  .filterDate('2015-01-01', '2016-12-31')
                  .filter(ee.Filter.calendarRange(6, 8, 'month'))
                  .map(maskL8sr);
print (dataset.size())

// ########### BLOCK: SELECTION OF SUITABLE LANDSAT IMAGES (END) ###################

// Specification of time window to use.
// One of 'year', 'month' 'week', 'day', 'hour', 'minute', or 'second'
var timeWindow = 6;
var timeWindowUnits = 'month';

// Add a new property 'daterange' to each plot element.
myplots = myplots.map(function(fc){
  var fc_date = ee.Date(fc.get('m_date'));
  return fc.set('daterange',ee.DateRange(fc_date.advance(-1*timeWindow,timeWindowUnits), fc_date.advance(timeWindow,timeWindowUnits)))
});

var makeFC = function(){
  // Define an empty Collection for the iterator to fill
  var empty_fc = ee.FeatureCollection(ee.List([]));
  // Iterate, and fill the (empty) collection
  var newfc = ee.FeatureCollection(dataset.iterate(getValsFromImage, empty_fc));
  return(newfc)
}

// Export the FeatureCollection to a CSV file.
Export.table.toDrive({
  collection: ee.FeatureCollection(makeFC()),
  folder: 'GEE_Exports',
  description:'pts_intWithImgs',
  fileFormat: 'CSV'
});

Answer 1

The following script should do the same thing as you've done, except it relies on .map() instead of .iterate(), which improves parallelization. However, it uses reduceRegion instead of reduceRegions so it may be a wash in terms of time. I was reluctant to use reduceRegions since it looks like you have points that span years and months that you extract band values for only to remove many points because they are outside of the current target time window being worked on - potentially lots of reprocessing of the same points. Not sure if the method below is any better, but might give you a new perspective on the problem. Also note that I removed the filter on cloud cover since you are masking clouds.

// Define point features.
var plotList = [
  ee.Feature(ee.Geometry.Point(25.26,67.03), {plot_id: 'plot_1', m_date: '2015-06-14'}),
  ee.Feature(ee.Geometry.Point(26.29,65.07), {plot_id: 'plot_2', m_date: '2015-08-28'}),
  ee.Feature(ee.Geometry.Point(28.39,62.55), {plot_id: 'plot_3', m_date: '2016-07-01'})
];

// Aggregate features to collection.
var plotsFc = ee.FeatureCollection(plotList);

// Define a Landsat 8 cloud masking function to be mapped over an
// ImageCollection in the following getValsFromImage function. 
function maskL8sr(image) {
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  var qa = image.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

// Define a function to extract Landsat 8 band values per point. Images
// are filtered by intersection with point and nearness to point date.
// Returns a FeatureCollection that is a table where each row (feature)
// is an image observation for a given point and columns (properties) include
// image band values, image ID, image date, and any existing point properties.
var getValsFromImage = function(point){
  
  // Get dates to filter images by.
  var middleDate = ee.Date(point.get('m_date'));
  var startDate =  middleDate.advance(-1, 'month');
  var endDate =  middleDate.advance(1, 'month');
  
  // Get images that intersect point and date range (define above).
  var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                    .filterBounds(point.geometry())
                    .filterDate(startDate, endDate)
                    .map(maskL8sr);
  
  // Map over ImageCollection to extract band values for each image intersecting
  // the given point.
  var pointPerImg = dataset.map(function(img) {
    var vals = img.select(ee.List.sequence(0, 8))
        .reduceRegion({
          reducer: ee.Reducer.first(),
          geometry: point.geometry(),
          scale: 30
        });
    
    return point.set(vals).set({img_id: img.id(), img_date: img.date().format('YYYY-MM-dd')});
  });
  
  // Return the image samples as a FeatureCollection.
  return ee.FeatureCollection(pointPerImg);
};


// Flatten the collection of collections; result is a FeatureCollection. Filter
// out masked (null) observations.
var pointsFcImgBands = plotsFc.map(getValsFromImage).flatten()
  .filter(ee.Filter.notNull(['B1']));

// Print the results - if large set of points see following not on batch task.
print(pointsFcImgBands);


// It's quite likely that for point feature collections that
// cover a large area or contain many points, you'll need to complete this
// operation as a batch task by either exporting the final feature collection as
// an asset or as a CSV/Shapefile/GeoJSON to Google Drive. If your browser
// times out when printing it, definitely try exporting the results.
// Export the FeatureCollection to a CSV file.
Export.table.toDrive({
  collection: pointsFcImgBands,
  folder: 'GEE_Exports',
  description:'pts_intWithImgs',
  fileFormat: 'CSV'
});

Code Editor script