0

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. enter image description here

2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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