0

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:

Image with cloud

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:

cloud_and_shadow_masks

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.

2
  • Thank you for your reply. Here are the bounding box coordinates for this image: POLYGON ((315182.01234129263 5873848.984715646, 317007.7388380545 5873848.984715646, 317007.7388380545 5875589.032606219, 315182.01234129263 5875589.032606219, 315182.01234129263 5873848.984715646)). They're in EPSG:32632. You can read them as a shapely geometry and then use: meadow_geom_ee = ee.Geometry.Polygon(mapping(meadow_envelope_geom)['coordinates'], proj='epsg:32632', evenOdd=False) where the meadow_envelope_geom is the shapely object. Oct 19, 2023 at 11:02
  • I deleted my comment, because I found the location before you replied -- see the answer below. Would you mind sharing the complete code (perhaps through a Colab notebook, since you are using python). Oct 19, 2023 at 11:07

2 Answers 2

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.

enter image description here

Here's the cloud_transform band in semi-transparency:

enter image description here

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?

1
  • Thanks again. Here is my code: link. I compared them and they look pretty much the same. I am not sure why my projection and result are so strange. I had to copy my code and change some variables to hard-code, and I couldnt test it because my google refuses to accept google colab notebooks in authentication. If something is missing or not working please let me know! Oct 19, 2023 at 13:21
1
+50

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:

enter image description here

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:

enter image description here

4
  • Thank you so much! it is strange that this is happening, cause my function is essentially theee_to_numpy but with the X and Y coordinates added. I need those so I can recreate the projection in Python and plot the results on top of my other rasters (for that I need the top leftcoordinates). I will try to get those in a different way. I will accept answer and reputation soon. May I be so inpolite to ask if you could add a small explanation of what we can tweak in the code to improve the shadow mask, bsides the obvious shadow threshold. I will accept anyway btw when I can in some hours Oct 19, 2023 at 14:07
  • So in EE coder all of this works. But in Python, something different happens. I tested it with GEE's ee_to_numpy function and EE's downloadUrl function. In both cases, the results are still like my original post. I Exactly copied the code from Earth Engine. I now also included the ee.downloadUrl variation in the notebook here link. So something different happens in Python vs EE coder. I have no clue why this could be, I am testing the projections but see nothing suspicious besides the output. Any ideas? Oct 20, 2023 at 11:38
  • I tried your notebook again, locally. It works as expected. The reason your image looks flipped is because you are plotting it without the actual coordinates (you only read the data as array using gdal, but not the coordinates). You are using imshow without any spatial coordinates, so it's showing you the data how it was read, but no spatial context. I suggest you either read the coordinates with gdal (e.g. see gis.stackexchange.com/questions/57834/… ) or use another method. Oct 22, 2023 at 5:21
  • I recommend using xarray, since it will read everything you need (coordinates and data) then when you plot it you will see the correct orientation of the image. In summary, there is nothing wrong in the ee code, the actual tif file is ok (you can also check in other gis software e.g. qgis). Oct 22, 2023 at 5:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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