NDWI anomalies calculation in Google Earth Engine

Question

I am a Google Earth Engine beginning level user and I am trying to detect NDWI anomalies using Sentinel-2 & Landsat Harmonized data and this formula

I have watched and read some tutorials and this code seems to be "semi-working" except for the fact the NDWI Anomaly layer appears to be a plain grey image. I am struggling to understand what my mistakes are here. My area of interest here is a geometry polygon I drew myself which is the administrative border of the city of Berlin, so it is not included in the code as it is way too big

var V_display = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
var ndwi_palette = ['red','white', 'blue'] ;
function addndwi (input) {
  var ndwi = input.normalizedDifference(['B8', 'B11']).rename('ndwi');
  return input.addBands(ndwi);
}

// 2017 ndwi - REF
//
// 2018.06 ndwi - W1
//
// 2017 Referrence period  NDWI (June - September)    
var REF = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(geometry)
  .filterDate('2017-06-01', '2017-09-30')
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 40)
  .map(addndwi);    
print(REF);    
var REF_m = REF.median().clip(geometry);    
Map.addLayer(REF_m, V_display, "Referrence 2017 median",false);    
var ndwi_2017_REF = REF.select('ndwi').median().clip(geometry);    
Map.addLayer(ndwi_2017_REF, {min: -1, max: 1, palette: ndwi_palette}, 'NDWI Ref Year 2017');

// Calculating the Standard deviation for the reference period    
var REF_SD = REF.reduce(ee.Reducer.stdDev()).clip(geometry);
print (REF_SD);
var W_display = {bands: ['B4_stdDev', 'B3_stdDev', 'B2_stdDev'], min: -10, max: 3000};    
Map.addLayer(REF_SD, W_display, "Referrence 2017 Standard Deviation",false);


// 2018 June NDWI    
var W1 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(geometry)
  .filterDate('2018-06-01', '2018-06-30')
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20)
  .map(addndwi);    
print(W1);    
var W1_m = W1.median().clip(geometry);    
Map.addLayer(W1_m, V_display, "June 2018 median",false);    
var ndwi_2018_06 = W1.select('ndwi').median().clip(geometry);    
Map.addLayer(ndwi_2018_06, {min: -1, max: 1, palette: ndwi_palette}, 'NDWI June 2018');    
Map.centerObject(geometry);

// Calculating the NDWI Anomalies for June 2018    
var dif1 = ndwi_2018_06.subtract(ndwi_2017_REF);    
var Andwi_2018_06 = dif1.divide(REF_SD);      
print(Andwi_2018_06);    
Map.addLayer(Andwi_2018_06,
             {bands: ['B4_stdDev', 'B3_stdDev', 'B2_stdDev'], min: -32, max: 32}, 'NDWI Anomaly 2018 June');

Answer 1

My area of interest here is a geometry polygon as I drew myself, so it is no included in the code

Tip: You can copy the code for a polygon, and include it in your script, by clicking on the blue icon next to “Imports”:

Even if the polygon is too complex to include, it's useful to tell us a latitude and longitude so we can try out your code in the same geographical area you are.

I have watched and read some tutorials and this code seems to be "semi-working" except for the fact the NDWI Anomaly layer appears to be a plain grey image.

I ran your code and looked at some sample values with the Inspector tab:

NDWI Anomaly 2018 June: Image (3 bands)
   B4_stdDev: 0.000028036509884596032
   B3_stdDev: 0.000027888740315101327
   B2_stdDev: 0.00002702838865872788

These are very small values, compared to your chosen min and max of ±32. This is often the answer any time your image is gray (or black or white) — your visualization parameters are out of range or the wrong scale.

Try something more like min: -0.001, max: 0.001 to see more detail.

You can also use the layer controls (the gear icon on hover in the layer list here) to change the visualization parameters, including an automatic stretch based on the range of values in the visible area of the image.