Harmonized Landsat-7 NDVI time series values are all the same when calculated with Google Earth Engine

Question

I wanted to calculate the NDVI time series of multiple points using Google Earth Engine. I'm using Collection 2 Level 2 Landsat-7 and Landsat-8 imagery to calculate the NDVI time series. However, when I try to harmonize the Landsat-7 ETM values to Landsat-8 OLI values using the tutorial, I get the same value of 0,74 for each NDVI value of the Landsat-7 images. I use .map to apply the harmınization function to all my images but can not solve the problem. How can I get the real/calculated NDVI values?

Here is my code:

 // Import the necessary libraries
var geometry = ee.FeatureCollection('users/maggie06/Kirsehir-5S');

//Image Collection Filter
var colFilter = ee.Filter.and(
    ee.Filter.bounds(geometry), ee.Filter.date('2014-09-22', '2017-07-02'),
    ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
    ee.Filter.or(
        ee.Filter.eq('IMAGE_QUALITY', 9),
        ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));

// Filter the images to only include images with less than 50% cloud coverage
var filtered_l7 = l7.filter(colFilter);
var filtered_l8 = l8.filter(colFilter);

// Applies scaling factors.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
//  var thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true);
//              .addBands(thermalBand, null, true);
}

var scaled_l7 = filtered_l7.map(applyScaleFactors);
var scaled_l8 = filtered_l8.map(applyScaleFactors);

// Rename the Landsat 7 bands in the collection
var renamed_l7 = scaled_l7.map(function(image) {
  return image.select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7', 'QA_PIXEL'],
                       ['blue', 'green', 'red', 'NIR', 'swir1', 'swir2', 'pixel_qa']);
});

// Rename the Landsat 8 bands in the collection
var renamed_l8 = scaled_l8.map(function(image) {
  return image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7', 'QA_PIXEL'],
                       ['blue', 'green', 'red', 'NIR', 'swir1', 'swir2', 'pixel_qa']);
});

// Mask clouds, cloud shadows, and snow
// Define the mask function using the pixel_qa band
var mask = function(img) {
  var qa = img.select(['pixel_qa']);
  // Clouds
  var cloud = qa.bitwiseAnd(1 << 3).neq(0);
  // Cloud shadows
  var shadow = qa.bitwiseAnd(1 << 4).neq(0);
  // Snow
  //var snow = qa.bitwiseAnd(1 << 5).neq(0);
  
  var mask = cloud.or(shadow);
  //.or(snow);
  return img.updateMask(mask.not());
};

var masked_l7 = renamed_l7.map(mask);
var masked_l8 = renamed_l8.map(mask);

//Harmonization Formula
var coefficients = {
  itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172]).multiply(10000),
  slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};

// Define function to apply harmonization transformation.
var harmonised_l7 = masked_l7.map(function (image) {
  return image.select(['blue', 'green', 'red', 'NIR', 'swir1', 'swir2'])
    .multiply(coefficients.slopes)
    .add(coefficients.itcps)
    .round()
    .toShort()
    .addBands(image.select('pixel_qa')
  );
});

// Apply the Transform function to the L7 ETM collection
//var harmonised_l7 = masked_l7.map(etm2oli);

// Merge the two collections
var merged_collections = harmonised_l7.merge(masked_l8);

// Define the NDVI function
var calcNDVI = function(img) {
  var ndvi = img.normalizedDifference(['NIR','red']).rename('NDVI');
  return img.addBands(ndvi);
};

// Apply the NDVI function to the collection
var ndvi_collection = merged_collections.map(calcNDVI);

// Function to map over the FeatureCollection - Saving NDVI values to GoogleDrive:
var extractNDVI = function(feat) {
  // get feature geometry
  var geom = feat.geometry();
  // function to iterate over the yearly ImageCollection
  // the initial object for the iteration is the feature
  var addProp = function(img, f) {
    // cast Feature
    var newf = ee.Feature(f);
    // get date as string
    //var date = img.date().format();
    var date = ee.String(img.get('system:index'));
    // extract the value (first) of 'waterClass' in the feature
    var value = img.reduceRegion(ee.Reducer.first(), geom, 30).get('NDVI');
    // if the value is not null, set the values as a property of the feature. The name of the property will be the date
    return ee.Feature(ee.Algorithms.If(value,
                                       newf.set(date, ee.String(value)),
                                       newf.set(date, ee.String('No data'))));
  };
  var ndvi_values = ee.Feature(ndvi_collection.iterate(addProp, feat));
  return ndvi_values;
};

var ndvi_values = geometry.map(extractNDVI);

// Export the NDVI values and dates to an Excel file
Export.table.toDrive({
  collection: ndvi_values,
  description: 'Kirsehir_5s_NDVI_harmonized',
  fileFormat: 'CSV'
});

Answer 1

I don't think that is the intended use of ee.ImageCollection.iterate, but there is a simpler way of retrieving time series over multiple points for an image collection: using ee.Image.reduceRegions.

Here's a full minimal example consisting of: (1) defining a small feature collection (2 points), (2) retrieving a simple NDVI image collection (from Landsat 9), (3) retrieving the time series as a FeatureCollection, (4) making a chart showing the time series over the two points. You can easily integrate step 3 with your feature colleciton and your NDVI image collection.

var fc = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(-123.75, 41.26), {name:"point1"}),
  ee.Feature(ee.Geometry.Point(-110.28, 38.44), {name:"point2"})]
)
Map.addLayer(fc)

var ndvi = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
.filterDate("2023-01-01", "2023-05-30")
.filterBounds(fc)
.map(function(img){
  return img.addBands(
    img.normalizedDifference(["SR_B5", "SR_B4"]).rename("NDVI")
    )
}).select("NDVI")

var timeSeries = ndvi.map(function(img){
return img.select("NDVI")
  .reduceRegions({
    collection: fc,
    reducer:ee.Reducer.first(),
    scale: 100  // see https://developers.google.com/earth-engine/guides/scale#scale-of-analysis
    })   // Use reduceRegions in img, to get the first() value for each feature (point) in the collection.
  .filter(ee.Filter.notNull(["first"])) // Don't keep the features in which the value is null 
  // (because each image doesn't cover all points in the feature collection)
  .select(["first", "name"],["NDVI", "name"]) // rename the property "first" (from the ee.Reducer.first() to NDVI)
  .map(function(f){
    return f.set({date:img.date(), 
    landsatIndex:img.get("system:index")})
  }) // Assign the "date" property to each feature.. and optionally other metadata from each Landsat image. 
}).flatten() // Map the above function over all the images, and flatten the resulting FeatureCollection
.sort("geometry")// Sort the features by date

print(timeSeries)

var chart = ui.Chart.feature.groups(timeSeries, "date", "NDVI","name")

print(
  chart.setOptions({"title": "NDVI time series over two points"})
  )

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