1

I am working with Sentinel-2 time series in Google Earth Engine (GEE). I would like to calculate the area under the curve (AUC) defined by a smoothed time series using a Savizky-Golay filter.

The method should compute the integral on the interval defined by the start and the end of the time series [x=a, x=b] (see Image) or other method to compute an approximation of the AUC. Any ideas on how to approach this? Or, how can I get the function which defines the curve?

enter image description here

This is my code to get the smoothed time series. Source: https://courses.spatialthoughts.com/end-to-end-gee.html

// Aplying Savitzky-Golay Filter on a NDVI Time-Series
// This script uses the OEEL library to apply a 
// Savitzky-Golay filter on a imagecollection

// We require a regularly-spaced time-series without
// any masked pixels. So this script applies
// linear interpolation to created regularly spaced images
// from the original time-series

// Step-1: Prepare a NDVI Time-Series
// Step-2: Create an empty Time-Series with images at n days
// Step-3: Use Joins to find before/after images
// Step-4: Apply linear interpolation to fill each image
// Step-5: Apply Savitzky-Golay filter
// Step-6: Visualize the results

//##############################################################
// Step-1: Prepare a NDVI time-series
//##############################################################

var s2 = ee.ImageCollection("COPERNICUS/S2"); 
var geometry = ee.Geometry.Polygon([[
  [82.60642647743225, 27.16350437805251],
  [82.60984897613525, 27.1618529901377],
  [82.61088967323303, 27.163695288375266],
  [82.60757446289062, 27.16517483230927]
]]);
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)

var startDate = ee.Date('2017-01-01')
var endDate = startDate.advance(1, 'year')

var filtered = s2
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
  .filter(ee.Filter.bounds(geometry))

// Write a function for Cloud masking
function maskS2clouds(image) {
  var qa = image.select('QA60')
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0))
  return image.updateMask(mask)//.divide(10000)
      .select("B.*")
      .copyProperties(image, ["system:time_start"])
}

var filtered = filtered.map(maskS2clouds)


// Write a function that computes NDVI for an image and adds it as a band
function addNDVI(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).toFloat().rename('ndvi');
  return image.addBands(ndvi);
}


// Map the function over the collection
var withNdvi = filtered.map(addNDVI);

// Select 'ndvi' band
var ndviCol = withNdvi.select('ndvi')
print('Original Collection', ndviCol)
//##############################################################
// Step-2: Create an empty Time-Series with images at n days
//##############################################################

// Select the interval. We will have 1 image every n days
var n = 5;
var totalDays = endDate.difference(startDate, 'day');
var daysToInterpolate = ee.List.sequence(1, totalDays, n)

var initImages = daysToInterpolate.map(function(day) {
  var image = ee.Image().rename('ndvi').toFloat().set({
    'system:index': ee.Number(day).format('%d'),
    'system:time_start': startDate.advance(day, 'day').millis(),
    // Set a property so we can identify interpolated images
    'type': 'interpolated'
  })
  return image
})

var initCol = ee.ImageCollection.fromImages(initImages)
print('Empty Collection', initCol)

//##############################################################
// Step-3: Use Joins to find before/after images
//##############################################################

// Merge empty collection with the original collection so we can
// find images to interpolate from
var mergedCol = ndviCol.merge(initCol)

var mergedCol = mergedCol.map(function(image) {
  var timeImage = image.metadata('system:time_start').rename('timestamp')
  var timeImageMasked = timeImage.updateMask(image.mask().select(0))
  return image.addBands(timeImageMasked)
})

// Specify the time-window
// Set it so that we have at least 1 non-cloudy image in the period
var days = 60
var millis = ee.Number(days).multiply(1000*60*60*24)

var maxDiffFilter = ee.Filter.maxDifference({
  difference: millis,
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})

var lessEqFilter = ee.Filter.lessThanOrEquals({
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})


var greaterEqFilter = ee.Filter.greaterThanOrEquals({
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})


var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)
var join1 = ee.Join.saveAll({
  matchesKey: 'after',
  ordering: 'system:time_start',
  ascending: false})
  
var join1Result = join1.apply({
  primary: mergedCol,
  secondary: mergedCol,
  condition: filter1
})

var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)

var join2 = ee.Join.saveAll({
  matchesKey: 'before',
  ordering: 'system:time_start',
  ascending: true})
  
var join2Result = join2.apply({
  primary: join1Result,
  secondary: join1Result,
  condition: filter2
})

//##############################################################
// Step-4: Apply linear interpolation to fill each image
//##############################################################

// Once the joins are done, we don't need original NDVI images
// We keep only the blank images which now have matching NDVI images
// as properties
var filtered = join2Result.filter(ee.Filter.eq('type', 'interpolated'))

// Interpolatinon function
function interpolateImages(image) {
  var image = ee.Image(image)

  var beforeImages = ee.List(image.get('before'))
  var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
  var afterImages = ee.List(image.get('after'))
  var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()

  var t1 = beforeMosaic.select('timestamp').rename('t1')
  var t2 = afterMosaic.select('timestamp').rename('t2')

  var t = image.metadata('system:time_start').rename('t')

  var timeImage = ee.Image.cat([t1, t2, t])

  var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
    't': timeImage.select('t'),
    't1': timeImage.select('t1'),
    't2': timeImage.select('t2'),
  })

  var interpolated = beforeMosaic
    .add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
  var result = image.unmask(interpolated)
  return result.copyProperties(image, ['system:time_start'])
}

var interpolatedCol = ee.ImageCollection(
  filtered.map(interpolateImages)).select('ndvi')
print('Interpolated Collection', interpolatedCol)


//##############################################################
// Step-5: Apply Savitzky-Golay filter
//##############################################################


var oeel=require('users/OEEL/lib:loadAll');
// https://www.open-geocomputing.org/OpenEarthEngineLibrary/#.ImageCollection.SavatskyGolayFilter

// Use the same maxDiffFilter we used earlier
var maxDiffFilter = ee.Filter.maxDifference({
  difference: millis,
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})

// Use the default distanceFunction
var distanceFunction = function(infromedImage, estimationImage) {
  return ee.Image.constant(
      ee.Number(infromedImage.get('system:time_start'))
      .subtract(
        ee.Number(estimationImage.get('system:time_start')))
        );
  }

// Apply smoothing of order=3
var order = 3;
var smoothed = oeel.ImageCollection.SavatskyGolayFilter(
  interpolatedCol, 
  maxDiffFilter,
  distanceFunction,
  order)

// Select the d_0_ndvi band and rename it
var smoothed = smoothed.select(['d_0_ndvi'], ['smoothed'])

//##############################################################
// Step-6: Visualize the results
//##############################################################

// Chart the time-series at a single location
var title = 'Savitsky-Golay smoothing' +
  '(order = '+ order + ', window_size = ' + days + ')'

// Plot the original and fitted NDVI time-series
var chart = ui.Chart.image.series({
  imageCollection: ndviCol.merge(smoothed),
  region: geometry,
  reducer: ee.Reducer.mean(),
  scale: 20
}).setOptions({
      lineWidth: 1,
      title: title,
      interpolateNulls: true,
      vAxis: {title: 'NDVI'},
      hAxis: {title: '', format: 'YYYY-MMM'},
      series: {
        0: {color: 'blue', lineWidth: 1, 
          lineDashStyle: [1, 1], pointSize: 1,
          }, // Original NDVI
        1: {color: 'red', lineWidth: 2 }, // Smoothed NDVI
      },

    })
print(chart)

1 Answer 1

1

I think the most expedite way to do that is by using Simpson method with smoothed values obtained of downloaded CSV chart. It can be accessed with click in the button inside red square of following picture.

enter image description here

and the download CSV button appears as follows:

enter image description here

Smoothed values (73; as number of elements in smoothed Image Collection) from CSV file were converted in a list variable at the beginning of following code. As these values have a spacing of 5 days then, the value of h is 5 in following lines.

//simpson method

var smoothed = [0.3,0.317,0.379,0.452,0.519,0.573,0.61,0.629,0.628,0.613,0.587,0.556,0.521,
                0.485,0.448,0.411,0.363,0.317,0.273,0.227,0.185,0.148,0.12,0.103,0.093,0.086,
                0.082,0.079,0.076,0.071,0.066,0.059,0.053,0.051,0.052,0.057,0.067,0.081,0.103,
                0.129,0.165,0.209,0.252,0.292,0.328,0.363,0.398,0.431,0.459,0.479,0.491,0.494,
                0.488,0.475,0.457,0.43,0.397,0.357,0.316,0.285,0.254,0.227,0.206,0.189,0.181,
                0.171,0.157,0.151,0.162,0.192,0.239,0.302,0.397];

var n = ee.List(smoothed).size().subtract(1);
print("number of values", n.add(1));

var h = 5;

var serverList = ee.List.sequence(1, n, 2);

print(serverList);

serverList = serverList.map(function(n){

  n = ee.Number(n).int();
  var n_sub = n.subtract(1);
  var n_add = n.add(1);

  var sum = ee.Number.expression(
      '(h/3)*(v1 + 4*v2 + v3)',{
      'v1': ee.List(smoothed).get(n_sub),
      'v2': ee.List(smoothed).get(n),
      'v3': ee.List(smoothed).get(n_add),
      'h': h
      }
    );
  
  return sum;
  
});

var area = ee.List(serverList).reduce(ee.Reducer.sum());

print("area under the curve (AUC)", ee.Number(area).format('%.2f'));

After running above code in GEE code editor, result is printed in Console Tab as:

number of values
73
List (36 elements)
area under the curve (AUC)
105.02

I corroborated that this code works as expected by using the values from here (article in Spanish) where result (h = 1) corresponds an area of 72 (parabola X2 between [0,6] ).

4
  • This is a big step toward the optimal solution. I will work on applying this method automatically (no need to download the .CSV). The optimal solution is to be able to apply this process (interpolation, S-G smoothing, and AUC) to each pixel of the image. Any help is welcome
    – sermomon
    Aug 1 at 12:48
  • 1
    The optimal solution is yours, not mine ( applying this method automatically without no needs to download the .CSV). You didn't know about numerical methods of integration for obtaining area under the curve of Sentinel-2 time-serie; the main subject of your question. Remember that in this site is one question by question.
    – xunilk
    Aug 1 at 14:13
  • 1
    I edited my question for deleting code with client-side for-loop. Now, it uses server-side mapping. Execution is much faster than for-loop version.
    – xunilk
    Aug 1 at 16:55
  • This solution fails when the number of smoothed values is an even number. How can I should fix this?
    – sermomon
    Sep 22 at 8:47

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.