Inaccurate cloud shadow masks when using Earth Engine's script for Sentinel 2

Question

I am currently evaluating my cloud shadow masking options for Sentinel-2 data. One of the options I am exploring is using an example script published by Google that projects cloud shadow locations based on the s2cloudless mask in Earth Engine, dark pixels, and the MEAN_SOLAR_AZIMUTH_ANGLE. Google has provided an example script here: https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless. When I use the code in the link for my data, the shadow masks are very inaccurate however. And I would really like to understand how I could improve this. Let's take the following image as an example:

This image is in EPSG:32632 in Earth Engine with transform [10, 0, 300000, 0, -10, 5900040]. I add the following bands to the image using the code in the link: 'probability' (cloud probabilities from s2cloudless), 'clouds' (pixels with a threshold over 30 in the 'probability' mask), 'dark_pixels' (dark pixels that are considered potential shadows based on threshold of 1500 in NIR band), 'cloud_transform' (the transform to project shadows from clouds), 'shadows' (the pixels considered shadows: overlap of 'dark_pixels' and cloud transform), 'cloudmask' (a combination of shadows and clouds as mask). The issue is that the shadows seem inaccurate. For the above image, I plotted some of the aforementioned bands, which look as follows:

where the first image are the dark pixels, the second is the cloud mask (cloud probabilities over 30), the third picture is the "cloud transform" band, and the last picture is the "shadow" band (which is the overlap of the dark pixels with the cloud transform). As you can see, the dark_pixels and cloud_mask look pretty good. The projection and final result however are completely off. When I look at the plots, it seems something is not right with the projection (third picture) as the shadow of the cloud (which is the blob in the top right) in the projection is rotated in a way that does not correspond with the truth at all. It almost looks mirrored in shape and is located too far to the right. So my guess is either something is wrong with how "shadow_azimuth" in the code below is calculated, or how "cld_proj" is defined, or both. I've been changing some of the values to see if it helps, but I do not completely understand yet what I can do here. I hope there is someone that has experience with this and can help me, or explain to me where things might be going wrong.

Interesting is that in EE's code editor and map display, this does not happen (see answers below). But when I copy the code into Python (and make it Pythonic) the cloud transform gets twisted. So there seems to be a difference between EE's code editor output, and what I get from for example ee's getDownloadUrl() function. See link here for the Python code.

def add_shadow_bands(img, NIR_DRK_THRESH=0.15, CLD_PRJ_DIST=1.5):
    # Identify water pixels from the SCL band.
    not_water = img.select('SCL').neq(6)
    # Identify dark NIR pixels that are not water (potential cloud shadow pixels).
    SR_BAND_SCALE = 1e4
    dark_pixels = img.select('B8').lt(NIR_DRK_THRESH*SR_BAND_SCALE).multiply(not_water).rename('dark_pixels')
    # Determine the direction to project cloud shadow from clouds (assumes UTM projection).
    shadow_azimuth = ee.Number(90).subtract(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')));
    # Project shadows from clouds for the distance specified by the CLD_PRJ_DIST input.
    cld_proj = (img.select('clouds').directionalDistanceTransform(shadow_azimuth, CLD_PRJ_DIST*10)
        .reproject(**{'crs': img.select(1).projection(), 'scale': 20})
        .select('distance')
        .mask()
        .rename('cloud_transform'))
    # Identify the intersection of dark pixels with cloud shadow projection.
    shadows = cld_proj.multiply(dark_pixels).rename('shadows')
    # Add dark pixels, cloud projection, and identified shadows as image bands.
    return img.addBands(ee.Image([dark_pixels, cld_proj, shadows]))

In case it helps, the image is part of the satellite image in Earth Engine with the following id: COPERNICUS/S2_SR/20210405T105021_20210405T105323_T32ULD.

Answer 1

Here's a minimal, complete example following the same tutorial over the location you show for that image.

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

var geometry = ee.Geometry.Point([6.263813807653249, 52.989728804949266]);
Map.centerObject(geometry, 14)
var v = {
  min:0,
  max:2000,
  bands:["B4","B3","B2"]
}
var CLD_PRB_THRESH = 30
var NIR_DRK_THRESH = 0.15
var CLD_PRJ_DIST = 1.5
var BUFFER = 50

var s2prob = ee.Image('COPERNICUS/S2_CLOUD_PROBABILITY/20210405T105021_20210405T105323_T32ULD')
var img = ee.Image("COPERNICUS/S2_SR/20210405T105021_20210405T105323_T32ULD")
.addBands(s2prob) 
.addBands(s2prob.gt(CLD_PRB_THRESH).rename("clouds"))

function add_shadow_bands(img){
  var not_water = img.select("SCL").neq(6)
  var SR_BAND_SCALE=1e4
  var dark_pixels = img.select("B8").lt(NIR_DRK_THRESH*SR_BAND_SCALE).multiply(not_water).rename('dark_pixels')
  var shadow_azimuth = ee.Number(90).subtract(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')));
  var cld_proj = (img.select('clouds').directionalDistanceTransform(shadow_azimuth, CLD_PRJ_DIST*10)
    .reproject({'crs': img.select(0).projection(), 'scale': 100})
    .select('distance')
    .mask()
    .rename('cloud_transform'))

  var shadows = cld_proj.multiply(dark_pixels).rename('shadows')
  return img.addBands(ee.Image([dark_pixels, cld_proj, shadows]))
}

function add_cld_shdw_mask(img){
    var is_cld_shdw = img.select('clouds').add(img.select('shadows')).gt(0)

    is_cld_shdw = (is_cld_shdw.focalMin(2).focalMax(BUFFER*2/20)
        .reproject({'crs': img.select([0]).projection(), 'scale': 20})
        .rename('cloudmask'))

    return img.addBands(is_cld_shdw)
}
img = add_shadow_bands(img)


Map.addLayer(img,v, "Image")
Map.addLayer(img.select("dark_pixels").selfMask(),{min:0,max:1, palette:["red"]}, "Dark Pixels", false)
Map.addLayer(img.select("clouds").selfMask(),{min:0,max:1, palette:["purple"]}, "Clouds")
Map.addLayer(img.select("cloud_transform").selfMask(),{min:0,max:1, palette:["pink"]}, "Cloud Transform", false)
Map.addLayer(img.select("shadows").selfMask(),{min:0,max:1, palette:["blue"]}, "shadows")

and here's the result showing clouds in purple and shadow in blue. As you can see, most of the shadow is identified.. of course you shouldn't expect a perfect match, especially for small clouds.

Here's the cloud_transform band in semi-transparency:

As you can see, I was unable to reproduce the bad cloud_transform that you show, which is pointing in the wrong direction. Would you mind sharing the complete code you used?

Answer 2

I took a look at the notebook you provided in the comment to my previous answer. It seems that your get_ee_cloudprobs function that you use to retrieve the data as a numpy array might not be working as you intend (I suspect there is a misrepresentation of the orientation in one or both of the x,y coordinates). I suggest you use other methods to do this, e.g. you could use the geemap package which includes a geemap.ee_to_numpy, or you could just export the image to drive as a tif file and then read it in python, etc..

Here's a screenshot of your notebook working as intended, using geemap (which now comes pre-installed) to visualize the results. Again, the results look fine as you can see below:

Note: I only had to modify one line of your code:

json_bounds = "{'type': 'FeatureCollection', 'features': [{'id': '0', 'type': 'Feature', 'properties': {}, 'geometry': {'type': 'Polygon', 'coordinates': [[[7.347377176243512, 52.22621087043558], [7.347377176243512, 53.23647271090539], [6.004760925069468, 53.23647271090539], [6.004760925069468, 52.22621087043558], [7.347377176243512, 52.22621087043558]]]}}]}"
featureCollection = ee.FeatureCollection(json_bounds)

to this:

json_bounds = {'type': 'FeatureCollection', 'features': [{'id': '0', 'type': 'Feature', 'properties': {}, 'geometry': {'type': 'Polygon', 'coordinates': [[[7.347377176243512, 52.22621087043558], [7.347377176243512, 53.23647271090539], [6.004760925069468, 53.23647271090539], [6.004760925069468, 52.22621087043558], [7.347377176243512, 52.22621087043558]]]}}]}
featureCollection = ee.FeatureCollection(json_bounds)

i.e. remove the double quotes from json_bounds.

For reference, this is how the cloud_transform looks with your get_ee_cloudprobs function: