Google Earth Engine – Calculate Stdev While Reducing Resolution

google-earth-enginereducers

I have a MODIS image that I'd like to resample to 12km resolution. For each 12km pixel, I'd like two pieces of information: (1) the mean of the MODIS pixels within the 12km pixel and (2) the standard deviation of the MODIS pixels within the 12km pixel.

Consider the following example:

var imageCollection = ee.ImageCollection("UMT/NTSG/v2/MODIS/NPP");    
var data = imageCollection.first().select("annualNPP")
    
var data_mean =data .reproject({ crs: 'EPSG:4326', scale: 12000})
                                        .reduceResolution({
                                          reducer: ee.Reducer.mean(), 
                                          maxPixels: 65535  }).rename("mean")

var data_stdev =data .reproject({ crs: 'EPSG:4326', scale: 12000})
                                        .reduceResolution({
                                          reducer: ee.Reducer.stdDev(), 
                                          maxPixels: 65535  }).rename("stdev")
Map.addLayer(data_mean,{},"mean")
Map.addLayer(data_stdev,{},"stdev")

This gives me an image with the mean of the MODIS pixels within the 12km pixel but the standard deviation at each 12km pixel is zero.

I thought maybe the standard deviation was small, so I tried multiplying the data_stdev image by 1,000, but that did not solve the problem.

Perhaps there is another way to approach this?

Best Answer

The reproject goes after the reduceResolution.

Related Solutions

Google Earth Engine – Calculating Gaussian Standardized Variable

I did a couple of changes:

Increased the CLOUDY_PIXEL_PERCENTAGE limit to 10, to make sure you get imagery for your date of interest.
Switched out the cloud masking to use the COPERNICUS/S2_CLOUD_PROBABILITY collection. It's slower but better.
Limited the date range used to calculate the statistics, not to include the landslide. When the landslide is part of the statistics, the pixel values afterwards wouldn't be an anomaly.

The below doesn't replicate the results exactly, but at least highlights the landslide.

var geometry = ee.Geometry.Polygon([
  [16.53793865386438,40.509788101097925],
  [16.552594243921753,40.509788101097925],
  [16.552594243921753,40.51740638148722],
  [16.53793865386438,40.51740638148722],
  [16.53793865386438,40.509788101097925]
])
var anomalyThreshold = -2
var cloudThreshold = 30

var filter = ee.Filter.and(
  ee.Filter.bounds(geometry),
  ee.Filter.date('2017-01-01', '2020-01-01')
)

// Uses slower but better cloud masking
var ndviCollection = ee.ImageCollection(
    ee.Join.saveFirst('cloudProbability').apply({
        primary: ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
          .filter(filter)
          .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)), // Increased the limit, to include imagery for your date
        secondary: ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
          .filter(filter),
        condition: ee.Filter.equals({leftField: 'system:index', rightField: 'system:index'})
    })
).map(toNdvi)

var stats = ndviCollection 
  .filterDate('2017-01-01', '2019-02-01') // Maybe limit this date range to some stable reference period?
  .reduce(
    ee.Reducer.mean()
        .combine(ee.Reducer.stdDev(), null, true)
  )
var gaussianCollection = ndviCollection.map(toGaussian)
var anomalyCount = gaussianCollection
  .filterDate('2019-02-01', '2019-12-01')
  .map(function (image) {
    return image.updateMask(image.lte(anomalyThreshold))
  })
  .reduce(ee.Reducer.count())
  .unmask(0)

var image = gaussianCollection
  .filterDate(ee.Date('2019-02-09').getRange('day'))
  .mosaic()

Map.addLayer(image, {min: -5, max: 5, palette: '#67001f, #b2182b, #d6604d, #f4a582, #fddbc7, #f7f7f7, #d1e5f0, #92c5de, #4393c3, #2166ac, #053061'}, 'image', true)
Map.addLayer(anomalyCount, {min: 0, max: 52, palette: '#042333, #2c3395, #744992, #b15f82, #eb7958, #fbb43d, #e8fa5b'}, 'anomaly count')
Map.centerObject(geometry)


function toGaussian(ndvi) {
  return ee.Image()
    .expression('(ndvi - ndvi_mean) / ndvi_stdDev', {
      ndvi: ndvi,
      ndvi_mean: stats.select('ndvi_mean'),
      ndvi_stdDev: stats.select('ndvi_stdDev')
    })
    .copyProperties(ndvi, ndvi.propertyNames())
}

function toNdvi(image) {
  var cloudFree = ee.Image(image.get('cloudProbability')).lt(cloudThreshold)
  return image.updateMask(cloudFree)
    .normalizedDifference(['B8', 'B4'])
    .float()
    .rename('ndvi')
    .copyProperties(image, image.propertyNames())
}

function maskS2clouds(image) {
  var qa = image.select('QA60')
  var cloudBitMask = 1 << 10
  var cirrusBitMask = 1 << 11
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
    .and(qa.bitwiseAnd(cirrusBitMask).eq(0))
  return image.updateMask(mask)
}

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

Related Question