2

I am trying to pan-sharpen a Landsat-8 image with an 'artificial' panchromatic band derived from calculating the average of the Sentinel-2 bands with a spatial resolution of 10 meters (i.e., Red, Green, Blue, NIR).

I am doing it in GEE with the code provided in the GEE documentation on spectral transformations. The process works well and I get a Landsat image with a resolution of 10 meters as a result, but only the RGB bands are retained in the final product, as the methodology provided in the documentation involves using the 'rgbToHsv()' and the 'hsvToRgb()' methods. How can I also pan-sharpen the infrared bands so that I end up getting a pan-sharpened product with the same bands as the original product instead of just the RGB bands?

This is my code:

var Spain = 
    /* color: #0b4a8b */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[-3.72860936293182, 40.36209427277739],
          [-3.72860936293182, 40.236409074308916],
          [-3.46493748793182, 40.236409074308916],
          [-3.46493748793182, 40.36209427277739]]], null, false),
    geometry2 = /* color: #d63000 */ee.Geometry.Point([-3.6054709910758365, 40.28987996402006]);



Map.centerObject(Spain, 10); // Polygon in Spain

var START_DATE = '2018-01-01'
var END_DATE = '2020-12-31'


////////////////////////////////////////////LANDSAT_IMAGE///////////////////////////////////


// Load a Landsat 8 top-of-atmosphere reflectance image.
var L8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA')
            .filterDate(START_DATE, END_DATE)
            .filterBounds(Spain)
            .sort('CLOUD_COVER')
            .first()
            .clip(Spain)

print(L8)

Map.addLayer(
    L8,
    {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.25, gamma: [1.1, 1.1, 1]},
    'rgb');


////////////////////////////////////////SENTINEL-2 IMAGE////////////////////////////////////


// get one S2 image
var S2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED") 
           .filterDate('2019-08-01', '2019-08-19')
           .filterBounds(geometry2) // Point within the polygon 'Spain'
           .sort('CLOUD_COVER')
           .first()
           .clip(Spain)
           
S2 = S2.select('B4', 'B3', 'B2', 'B8') //Only keep bands with a 10-meter resolution


//////////////////////////////////// NORMALIZE S2 IMAGE ////////////////////////////////////

// The method "hsvToRgb()", which will be used later to pan-sharpen the image, requires three normalized bands.
// The first two will be normalized later by the method "rgbToHsv()". As for the third, the code
// below normalizes the Sentinel-2 bands, from wihichit will be calculated, to the range [0-1], .

// calculate the min and max value of an image
var minMax = S2.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: S2.geometry(),
  scale: 10,
  maxPixels: 10e9,
  // tileScale: 16
}); 

print(minMax)

// use unit scale to normalize the pixel values
var unitScale = ee.ImageCollection.fromImages(
    S2.bandNames().map(function(name){
    name = ee.String(name);
    var band = S2.select(name);
    return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))))
                // multiply by 100 if you want to get range 0-100
                //.multiply(100);
})).toBands().rename(S2.bandNames());
  
// add to the map
Map.addLayer(S2, {min: 0, max: 3500, bands: ['B4', 'B3', 'B2']}, 'original')
Map.addLayer(unitScale, {min: 0, max: 0.25, bands: ['B4', 'B3', 'B2']}, 'unitscaled')


////////////////////////// CALCULATE ARTIFICIAL PANCHROMATIC BAND //////////////////////////

// Calculate 'artificial panchromatic band' from Sentinel-2 bands with a 10-meter resolution.
var S2_image = unitScale.expression("((B4+B3+B2+B8)/4)",
  {
    B4 : unitScale.select('B4'),
    B3: unitScale.select('B3'),
    B2: unitScale.select('B2'),
    B8: unitScale.select('B8')
  }).rename('Panchromatic');

Map.addLayer(
    S2_image.select('Panchromatic'),
    {min: 0, max: 0.25},
    'Panchromatic');  


////////////////////////////////////// PAN-SHARPENING //////////////////////////////////////


// Convert the RGB bands to the HSV color space.
var hsv = L8.select(['B4', 'B3', 'B2']).rgbToHsv(); // The resulting bands are normalized

// Swap in the panchromatic band and convert back to RGB.
var sharpened = ee.Image.cat([
  hsv.select('hue'), hsv.select('saturation'), S2_image.select('Panchromatic') //Here, I include the 'artificial' Sentinel-2 normalized panchromatic band.
]).hsvToRgb();

// Display the pan-sharpened result (10-meter resolution -> same as for the panchromatic band)
Map.addLayer(sharpened,
             {min: 0, max: 0.25, bands:['red','green','blue']},
             'pan-sharpened');

print(sharpened)

2 Answers 2

1

Here's one approach based of feedback from Noel in a thread on the Google Earth Engine Developers group. I think it should be OK to do like this. This assumes the band names are the same in L8 and S2.

var sharpened = sharpen(L8, S2)

function sharpen(toSharpen, image) {
  var low = image.reproject(image.projection().scale(3, 3))  // Force a fixed projection of 30m.
  var offset = image.subtract(low).reproject(image.projection())
  return toSharpen
    .add(offset)
}

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

UPDATE: Corrected the rescaling and removed the unused kernel variable. Thanks for spotting it Diego!

3
  • Thanks for the answer. The only thing, it seems that you are reprojecting "image" to a 20-meter resolution. To force a fixed projection of 30 meters you need to scale the projection by 3: "image.projection().scale(3, 3)" Dec 12, 2022 at 12:56
  • One question, what is the function of the "kernel" variable within the "sharpen()" function? The variable does not seem to be used later in the function and the pan-sharpening works well without it Dec 12, 2022 at 13:14
  • Thanks Diego - I've updated the answer. I was too quick in copy/pasting from the other script. Dec 14, 2022 at 12:56
1
var S2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED")
  .filterDate('2019-08-01', '2019-08-19')
  .filterBounds(table)
  .sort('CLOUDY_PIXEL_PERCENTAGE')
  .first()
  .select( ['B4', 'B3', 'B2', 'B8'] );
  
  print (S2)

I found this mistake in the .sort case you need to use CLOUDY_PIXEL_PERCENTAGE to Sentinel instead of CLOUD_COVER cause this is to LANDSAT COLLECTION.

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.