0

I am trying to derive 5 statistical values: 1) maximum; 2) minimum; 3) standard deviation; 4) medium; and 5) 1090 interval mean from the time series data. And I know there are existed functions in GEE which support the above calculations:

ee.Reducer.min()
ee.Reducer.max()
ee.Reducer.stdDev()
ee.Reducer.median()
ee.Reducer.intervalMean(10, 90)

In the following code, the user can specify multiple points (var points) and a period (var globOptions.startDate, var globOptions.endDate) to collect the data from Landsat images. I specified 5 points and the period from Jan 1st 2014 to Dec 31st 2014 in this example. By running this code, we can get two CSV files, they are the calculation results of NDVI and NDWI. Each table has 5 rows, which correspond to the 5 points I have specified. And each column represent a Landsat image in the year of 2014. That is, we can calculate a NDVI value and a NDWI value at any single specified point on any single Landsat image in 2014. In case no value can be calculated at a specified point on a specified Landsat image, the null value is assigned as -99999.

var points = /* color: #d63000 */ee.FeatureCollection(
        [ee.Feature(
            ee.Geometry.Point([-81.28391315624998, 25.966529677689568]),
            {
              "system:index": "0"
            }),
        ee.Feature(
            ee.Geometry.Point([-81.03122760937498, 25.917133026546388]),
            {
              "system:index": "1"
            }),
        ee.Feature(
            ee.Geometry.Point([-80.93235065624998, 25.590597209051275]),
            {
              "system:index": "2"
            }),
        ee.Feature(
            ee.Geometry.Point([-81.00376178906248, 25.268135383731057]),
            {
              "system:index": "3"
            }),
        ee.Feature(
            ee.Geometry.Point([-81.42124225781248, 25.650034070641205]),
            {
              "system:index": "4"
            })]);

var site = ee.Geometry.Polygon([[-80.79, 24.98], [-80.08, 26.43], [-80.71, 26.16],[-81.07, 26.20]], null, false);
var globOptions = {
  versionID: '_SR',
  outFolder: 'SR',
  startDate: '2014-01-01',
  endDate: '2014-12-31',
  bandSelect: ['blue', 'green', 'red', 'nir', 'swir1', 'swir2'],
  bands8: ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
  bands7: ['B1', 'B2', 'B3', 'B4', 'B5', 'B7'], 
  maskAltitude: 100,  
  maskDepth: -100, 
  maskDistance: 50000,
  maskApplySRTM: false,
  parallelScale: 8,
  trainingValidationRatio: 0.0001,
  nTrees: 10, 
  outScale: 30, 
  conPixels: 100
};

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);
  };
  
  // A function to mask out cloudy shadow pixels.
var cloud_shadows = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
      // Return an image masking out cloudy areas.
  };
  
  // A function to mask out cloud pixels.
var clouds = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 5,5, 'Cloud').eq(0);
  // Return an image masking out cloudy areas.
};

var maskClouds = function(image) {
  var cs = cloud_shadows(image);
  var c = clouds(image);
  image = image.updateMask(cs);
  return image.updateMask(c);
  };
//// Import the Landsat 8 TOA image collection.
//var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');

  var L4collection = ee.ImageCollection('LANDSAT/LT04/C01/T1_SR')
      .filterDate(globOptions.startDate, globOptions.endDate)
      .map(maskClouds)
      .select(globOptions.bands7, globOptions.bandSelect);
  var L5collection = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
      .filterDate(globOptions.startDate, globOptions.endDate)
      .map(maskClouds)
      .select(globOptions.bands7, globOptions.bandSelect);
  var L7collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
      .filterDate(globOptions.startDate,globOptions.endDate)
      .map(maskClouds)
      .select(globOptions.bands7, globOptions.bandSelect);
  var L8collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
      .filterDate(globOptions.startDate, globOptions.endDate)
      .map(maskClouds)
      .select(globOptions.bands8, globOptions.bandSelect);
  var collectionFull = ee.ImageCollection(L4collection
  .merge(L5collection)
  .merge(L7collection)
  .merge(L8collection))


var getImage = function(id) {
  return ee.Image(collection.filter(ee.Filter.eq('system:index', id)).first())
}

var points = points.map(function(feature) {
  return ee.Feature(feature.geometry(), {'id': feature.id()})
})

//NDVI

function addNDVI(image) {
  var ndvi = image.normalizedDifference(['nir','red']).rename('ndvi')
  return image.addBands([ndvi])
}

var collection = collectionFull
    .map(addNDVI)
    .filter(ee.Filter.bounds(points))
    
var triplets = collection.map(function(image) {
  return image.select('ndvi').reduceRegions({
    collection: points, 
    reducer: ee.Reducer.first().setOutputs(['ndvi']), 
    scale: 10,
  })// reduceRegion doesn't return any output if the image doesn't intersect
    // with the point or if the image is masked out due to cloud
    // If there was no ndvi value found, we set the ndvi to a NoData value -99999
    .map(function(feature) {
    var ndvi = ee.List([feature.get('ndvi'), -99999])
      .reduce(ee.Reducer.firstNonNull())
    return feature.set({'ndvi': ndvi, 'imageID': image.id()})
    })
  }).flatten();
  
  var format = function(table, rowId, colId) {
  var rows = table.distinct(rowId); 
  var joined = ee.Join.saveAll('matches').apply({
    primary: rows, 
    secondary: table, 
    condition: ee.Filter.equals({
      leftField: rowId, 
      rightField: rowId
    })
  });
        
  return joined.map(function(row) {
      var values = ee.List(row.get('matches'))
        .map(function(feature) {
          feature = ee.Feature(feature);
          return [feature.get(colId), feature.get('ndvi')];
        });
      return row.select([rowId]).set(ee.Dictionary(values.flatten()));
    });
};

var LandsatResults = format(triplets, 'id', 'imageID');

Export.table.toDrive({
    collection: LandsatResults,
    description: 'NDVI_time_series',
    folder: 'earthengine',
    fileNamePrefix: 'ndvi_time_series',
    fileFormat: 'CSV'
})

//NDWI

function addNDWI(image) {
  var ndwi = image.normalizedDifference(['green','nir']).rename('ndwi')
  return image.addBands([ndwi])
}

var collection2 = collectionFull
    .map(addNDWI)
    .filter(ee.Filter.bounds(points))
    
var triplets2 = collection2.map(function(image) {
  return image.select('ndwi').reduceRegions({
    collection: points, 
    reducer: ee.Reducer.first().setOutputs(['ndwi']), 
    scale: 10,
  })
    .map(function(feature) {
    var ndwi = ee.List([feature.get('ndwi'), -99999])
      .reduce(ee.Reducer.firstNonNull())
    return feature.set({'ndwi': ndwi, 'imageID': image.id()})
    })
  }).flatten();
  
  var format2 = function(table, rowId, colId) {
  var rows = table.distinct(rowId); 
  var joined = ee.Join.saveAll('matches').apply({
    primary: rows, 
    secondary: table, 
    condition: ee.Filter.equals({
      leftField: rowId, 
      rightField: rowId
    })
  });
        
  return joined.map(function(row) {
      var values = ee.List(row.get('matches'))
        .map(function(feature) {
          feature = ee.Feature(feature);
          return [feature.get(colId), feature.get('ndwi')];
        });
      return row.select([rowId]).set(ee.Dictionary(values.flatten()));
    });
};

var LandsatResults2 = format2(triplets2, 'id', 'imageID');

Export.table.toDrive({
    collection: LandsatResults2,
    description: 'NDWI_time_series',
    folder: 'earthengine',
    fileNamePrefix: 'ndwi_time_series',
    fileFormat: 'CSV'
})

In other words, this code gives us the original data. And what I want is, by inserting the 5 existed GEE functions into my code (or any other alternative way is welcomed), we can export a table of statistical values based on the original data. The export table is expected to have five rows (represent five specified points) and ten columns (2 indices * 5 statistical values). There is another tricky thing: before calculating the statistical values, we should firstly discard those -99999 values, because they are actually null value. Any hint or suggestion is appreciated.

2

EE works best with images, so in general, try to wait converting your data into features as long as possible. I'd do things in this order:

  1. Create time-series image collection (add NDVI, NDWI, mask clouds)
  2. Calculate the statistics on that image collection, by calling ee.ImageCollection.reduce()
  3. Extract data for your points by calling reduceRegions() on the resulting image.
  4. Export the feature collection.

https://code.earthengine.google.com/ec58b1afd837f6460e9e34cf320b65fa

function preProcess(image) {
  image = addNDWI(image)
  image = addNDVI(image)     
  image = maskClouds(image)
  return image.select(['ndvi', 'ndwi'])
}

var stats = collectionFull
  .filterBounds(points)
  .map(preProcess)
  .reduce(
    ee.Reducer.minMax()
      .combine(ee.Reducer.stdDev(), '', true)
      .combine(ee.Reducer.median(), '', true)    
      .combine(ee.Reducer.intervalMean(10, 90), '', true)    
  )
  .reduceRegions({
    collection: points, 
    reducer: ee.Reducer.firstNonNull(), 
    scale: 30
  })

Export.table.toDrive({
    collection: stats,
    description: 'NDWI_time_series',
    folder: 'earthengine',
    fileNamePrefix: 'ndwi_time_series',
    fileFormat: 'CSV'
})
1
  • It works great! Thank you! – Tsui Raymond Sep 16 '20 at 11:07

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.