[GIS] Pansharpening Sentinel-2 imagery in Google Earth Engine

Question

I am attempting to pansharpen Sentinel-2 imagery in GEE in order to determine the amount of land covered in water at my field site locations (probably will calculate MNDWI on sharpened image). I do not have any experience in GEE or with pansharpening, so I am not sure whether my issue is a coding issue or a mathematical issue. I have been following this code < https://leclab.wixsite.com/spatial/post/pansharpening-sentinel-2-imagery-in-google-earth-engine >, which pansharpens the Sentinel imagery to all 10m bands using a high-pass filter (HPF) resolution merge. I have been following the code using my own Sentinel-2 scene, but am receiving an "internal error" when I try to print the result of the MSstdev calculation in Section 5c. I can print all images and numbers before this step, and nothing after this step will print (although the entire code will run without any errors).

I have included a snippet of the code below where the error occurs, if anyone might know what the issue is.

Answer 1

Thanks for sharing the article. I haven't seen this before, and I'm sure I'll get some use of this at one point or another. Here's my stab at getting the code running, and cleaning up a bit:

var image = ee.Image('COPERNICUS/S2/20190630T100031_20190630T100212_T32TQM')
var sharpened = panSharpen({
  image: image,
  bestEffort: true
})

print('sharpened', sharpened)
var visParams = {bands: 'B5,B8A,B12', min: 600, max: 4000}
Map.addLayer(image, visParams, 'image')
Map.addLayer(sharpened, visParams, 'sharpened')
Map.centerObject(ee.Geometry.Point([12.49248, 41.89015]), 14)


/**
 * Pansharpens a Sentinel 2 image.
 * 
 * Arguments:
 * 
 * params - a client-side object containing:
 *  
 *    image (Image, required) 
 *        The image to pansharpen
 * 
 *    geometry (Geometry, default: image.geometry()) 
 *        The region to pansharpen
 * 
 *    crs (Projection, default: projection of image's first band)
 *        The projection to work in.
 * 
 *    maxPixels (Long, default: 10000000)
 *        The maximum number of pixels to reduce.
 * 
 *    bestEffort (Boolean, default: false)
 *        If the geometry would contain more pixels than maxPixels, 
 *        compute and use a larger scale which would allow the operation to succeed.
 * 
 *    tileScale (Float, default: 1)
 *        A scaling factor used to reduce aggregation tile size; 
 *        using a larger tileScale (e.g. 2 or 4) may enable computations 
 *        that run out of memory with the default.
 */
function panSharpen(params) {
  if (params && !(params.image instanceof ee.Image))
    throw Error('panSharpen(params): You must provide an params object with an image key.')

  var image = params.image
  var geometry = params.geometry || image.geometry()
  var crs = params.crs || image.select(0).projection()
  var maxPixels = params.maxPixels
  var bestEffort = params.bestEffort || false
  var tileScale = params.tileScale || 1

  image = image.clip(geometry)

  var bands10m = ['B2', 'B3', 'B4', 'B8']
  var bands20m = ['B5', 'B6', 'B7', 'B8A', 'B11', 'B12']
  var panchromatic = image
    .select(bands10m)
    .reduce(ee.Reducer.mean())
  var image20m = image.select(bands20m)
  var image20mResampled = image20m.resample('bilinear')

  var stats20m = image20m
    .reduceRegion({
      reducer: ee.Reducer.stdDev().combine(
        ee.Reducer.mean(), null, true
      ),
      geometry: geometry,
      scale: 20,
      crs: crs, 
      bestEffort: bestEffort, 
      maxPixels: maxPixels, 
      tileScale: tileScale
    })
    .toImage()

  var mean20m = stats20m
    .select('.*_mean')
    .regexpRename('(.*)_mean', '$1')

  var stdDev20m = stats20m
    .select('.*_stdDev')
    .regexpRename('(.*)_stdDev', '$1')

  var kernel = ee.Kernel.fixed({
    width: 5,
    height: 5, 
    weights: [
      [-1, -1, -1, -1, -1],
      [-1, -1, -1, -1, -1],
      [-1, -1, 24, -1, -1],
      [-1, -1, -1, -1, -1],
      [-1, -1, -1, -1, -1]
    ], 
    x: -3, 
    y: -3, 
    normalize: false
  })

  var highPassFilter = panchromatic
    .convolve(kernel)
    .rename('highPassFilter')

  var stdDevHighPassFilter = highPassFilter
    .reduceRegion({
      reducer: ee.Reducer.stdDev(),
      geometry: geometry,
      scale: 10,
      crs: crs, 
      bestEffort: bestEffort, 
      maxPixels: maxPixels, 
      tileScale: tileScale
    })
    .getNumber('highPassFilter')

  function calculateOutput(bandName) {
    bandName = ee.String(bandName)
    var W = ee.Image().expression(
      'stdDev20m / stdDevHighPassFilter * modulatingFactor', {
        stdDev20m: stdDev20m.select(bandName),
        stdDevHighPassFilter: stdDevHighPassFilter,
        modulatingFactor: 0.25
      }
    )
    return ee.Image()
      .expression(
        'image20mResampled + (HPF * W)', {
          image20mResampled: image20mResampled.select(bandName),
          HPF: highPassFilter,
          W: W
      }
    )
    .uint16()
  }

  var output = ee.ImageCollection(
      bands20m.map(calculateOutput)
    )
    .toBands()
    .regexpRename('.*_(.*)', '$1')

  var statsOutput = output
    .reduceRegion({
      reducer: ee.Reducer.stdDev().combine(
        ee.Reducer.mean(), null, true
      ),
      geometry: geometry,
      scale: 10,
      crs: crs, 
      bestEffort: bestEffort, 
      maxPixels: maxPixels, 
      tileScale: tileScale
    })
    .toImage()

  var meanOutput = statsOutput
    .select('.*_mean')
    .regexpRename('(.*)_mean', '$1')

  var stdDevOutput = statsOutput
    .select('.*_stdDev')
    .regexpRename('(.*)_stdDev', '$1')

  var sharpened = ee.Image()
    .expression(
      '(output - meanOutput) / stdDevOutput * stdDev20m + mean20m', {
        output: output,
        meanOutput: meanOutput,
        stdDevOutput: stdDevOutput,
        stdDev20m: stdDev20m,
        mean20m: mean20m
      }
    )
    .uint16() 

  return image
    .addBands(sharpened, null, true)
    .select(image.bandNames())
}  

https://code.earthengine.google.com/86ccd4ff98f3806eb7290ef0cf926782