2

I am trying to make an image composite from Landsat 8 collection 2 based on Geometric Median (GM) algorithm. This is another option (like median or medoid) for creating a composite for time series analysis (more detail on GM can be found in this paper).

There is a reducer in GEE to compute GM, ee.Reducer.geometricMedian, and I tried it on my image collection. But there is some discrepancies and problematic patches in the composite which I don't know how to handle.

For comparison, I make 2 composite, one with GM and the other with classical median. enter image description here

// Area of Interest
var aoi = ee.Geometry.MultiPoint([
[46.72759630878907,38.171211790897694],
[48.38103869160157,38.171211790897694],
[48.38103869160157,39.25559070083476],
[46.72759630878907,39.25559070083476],
[46.72759630878907,38.171211790897694]
  ]).bounds()
  
Map.centerObject(aoi, 8)
Map.addLayer(aoi, {}, 'Area of Interest')

// Defining the period for time series 
var startDate = '2018-01-01'
var endDate = '2019-12-31'

// Collection the image collection
var l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
var l8Filt = l8.filterBounds(aoi).filterDate(startDate, endDate)

// Function for masking low-quality pixels 
function prepSrL8(image) {
  // Develop masks for unwanted pixels (fill, cloud, cloud shadow).
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);


  // Apply the scaling factors to the appropriate bands.
  var getFactorImg = function(factorNames) {
    var factorList = image.toDictionary().select(factorNames).values();
    return ee.Image.constant(factorList);
  };
  var scaleImg = getFactorImg([
    'REFLECTANCE_MULT_BAND_.|TEMPERATURE_MULT_BAND_ST_B10']);
  var offsetImg = getFactorImg([
    'REFLECTANCE_ADD_BAND_.|TEMPERATURE_ADD_BAND_ST_B10']);
  var scaled = image.select('SR_B.|ST_B10').multiply(scaleImg).add(offsetImg);

  // Replace original bands with scaled bands and apply masks.
  return image.addBands(scaled, null, true)
    .updateMask(qaMask).updateMask(saturationMask);
}

// Median Composite 
var composite_Median = l8Filt.map(prepSrL8).median().clip(aoi);
Map.addLayer(composite_Median, {
              bands: ['SR_B4','SR_B3','SR_B2'], min: 0, max: 0.2
              },'Median Composite')
              
// Geometric Median Composite 
var composite_GM = l8Filt.map(prepSrL8).reduce(ee.Reducer.geometricMedian(19)).clip(aoi)
print ('Geometricla Median', composite_GM)
Map.addLayer(composite_GM, 
            {bands: ['geometric_median_3','geometric_median_2','geometric_median_1'], 
            min: 0, max: 0.2},
            'GM Composite')

code snapshot

1 Answer 1

3

Geometric median is based on distance in all bands, and you're including all 19 bands in the distance calculation. That includes the distance to the bit-packed saturation and QA bands, as well as the distance to cloud, etc. (So, that's going to result in garbage).

Your bands also aren't normalized; the optical bands are 0 to 1 but the temperature band is in Kelvin. So the distance in the thermal band is going to dominate the computation. (So that's also going to result in garbage)

And finally, the thermal band is pretty messy in this region for the specified time. It's not even fully populated, so the geometric median can't even be computed in places.

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

Just use the optical bands.

1
  • Many thanks for your advice on this Noel. It perfectly works now.
    – rez
    Nov 9, 2022 at 23:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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