0

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'
});

2
  • 1
    .iterate() could be slowing this down. You should rearrange so that you .map() over your feature collection of points. Essentially, for each point, filter the Landsat collection by bounds and nearest date to point date, apply cloud mask (any other image processing), reduce the image to get values for point intersection, return the current feature (point) with the image values .set() as properties. If you can't get it going with this tip or someone else does not answer, I'll answer in a few days. Jun 11 '20 at 18:53
  • Thanks! I am not entirely sure how to do this but will try (on monday) and let you know. Jun 12 '20 at 13:23
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

1
  • Hi Justin, Thanks so much for that. I saw performance improvements with your code in the toy example. That is, your code executed in 11s, while mine had taken 15s. I am pretty sure it should improve the execution time with my 60K plots, and six year. I have not found time to adapt your code to what I have (without breaking some down-stream stuff). But hopefully, I will get to it soon, and will then post some performance improvement numbers with this bigger set, too. Thanks a lot! Jun 26 '20 at 10:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.