# Area under curve of Sentinel-2 time-series using Google Earth Engine

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?

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.centerObject(geometry)

var startDate = ee.Date('2017-01-01')

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
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
.select("B.*")
.copyProperties(image, ["system:time_start"])
}

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

// Map the function over the collection

// 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'),
// 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) {
})

// 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 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
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
//##############################################################

// 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)
``````

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.

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

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 sum = ee.Number.expression(
'(h/3)*(v1 + 4*v2 + v3)',{
'v1': ee.List(smoothed).get(n_sub),
'v2': ee.List(smoothed).get(n),
'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] ).

• 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. Aug 1, 2022 at 14:13
• 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. Aug 1, 2022 at 16:55
• This solution fails when the number of smoothed values is an even number. How can I should fix this? Sep 22, 2022 at 8:47
• @sermomon Sorry, just now I watch your comment. First, calculate for n-1 points (odd number of points). Afterward, add the last area as a single trapezoid. You should modify slightly the code for that. Nov 24, 2022 at 17:38