0

I have been working on a project requiring me to work with Sentinel-2 data on cloud and shadow masking. I have used the following code for masking clouds and cloud shadows using Google Earth Engine. The code is written in Python but I am comfortable with a JavaScript answer too.

COUNTRY_NAME = 'Malawi'
AOI == (ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017").
       filter(ee.Filter.eq('country_na', COUNTRY_NAME)).geometry()) 

START_DATE = '2020-01-01'
END_DATE = '2020-12-31'
CLOUD_FILTER = 30
CLD_PRB_THRESH = 50
NIR_DRK_THRESH = 0.15
CLD_PRJ_DIST = 1
BUFFER = 50`

Here we are selecting and building Sentinel-2 collection -

def get_s2_sr_cld_col(aoi, start_date, end_date):
   # Import and filter S2 SR.
   s2_sr_col = (ee.ImageCollection('COPERNICUS/S2_SR')
       .filterBounds(aoi)
       .filterDate(start_date, end_date)
       .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', CLOUD_FILTER)))

   # Import and filter s2cloudless.
   s2_cloudless_col = (ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
             .filterBounds(aoi)
             .filterDate(start_date, end_date))

   # Join the filtered s2cloudless collection to the SR collection by the 'system:index' property.
   return ee.ImageCollection(ee.Join.saveFirst('s2cloudless').apply(**{
       'primary': s2_sr_col,
       'secondary': s2_cloudless_col,
       'condition': ee.Filter.equals(**{
            'leftField': 'system:index',
            'rightField': 'system:index'
        })
    }))

Next, we are adding the cloud components -

def add_cloud_bands(img):
    # Get s2cloudless image, subset the probability band.
    cld_prb = ee.Image(img.get('s2cloudless')).select('probability')

    # Condition s2cloudless by the probability threshold value.
    is_cloud = cld_prb.gt(CLD_PRB_THRESH).rename('clouds')

    # Add the cloud probability layer and cloud mask as image bands.
    return img.addBands(ee.Image([cld_prb, is_cloud]))

Next, we are adding cloud shadow components -

def add_shadow_bands(img):
    # 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(0).projection(), 'scale': 100})
    .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]))    

Next, we are defining a function to assemble all of the cloud and cloud shadow components and produce the final mask -

def add_cld_shdw_mask(img):
    # Add cloud component bands.
    img_cloud = add_cloud_bands(img)

    # Add cloud shadow component bands.
    img_cloud_shadow = add_shadow_bands(img_cloud)

    # Combine cloud and shadow mask, set cloud and shadow as value 1, else 0.
    is_cld_shdw = img_cloud_shadow.select('clouds').add(img_cloud_shadow.select('shadows')).gt(0)

    # Remove small cloud-shadow patches and dilate remaining pixels by BUFFER input.
    # 20 m scale is for speed, and assumes clouds don't require 10 m precision.
    is_cld_shdw = (is_cld_shdw.focalMin(2).focalMax(BUFFER*2/20)
    .reproject(**{'crs': img.select([0]).projection(), 'scale': 20})
    .rename('cloudmask'))

    # Add the final cloud-shadow mask to the image.
    return img_cloud_shadow.addBands(is_cld_shdw)

Defining and applying cloud mask to each image in the collection

def apply_cld_shdw_mask(img):
    # Subset the cloudmask baand and invert it so clouds/shadow are 0, else 1.
    not_cld_shdw = img.select('cloudmask').Not()

    # Subset reflectance bands and update their masks, return the result.
    return img.select('B.*').updateMask(not_cld_shdw)

Then, making a function to split all images in a list for each month starting from January and ending in December.

Now, calling all the functions.

s2_sr_cld_col = get_s2_sr_cld_col(AOI, START_DATE, END_DATE)

s2_sr_cld_mask_add = (s2_sr_cld_col.map(add_cld_shdw_mask))

s2_sr_cld_col_list = split_image_collection_by_month(s2_sr_cld_mask_add)

s2_median_list = []
for image in s2_sr_cld_col_list:
s2_median_list.append(image.map(apply_cld_shdw_mask)
                         .median())

Now, I want to do an average superimpose or mosaicing (I am not sure with the correct word). So, I want to do something like this, so, for example, if I select February from the list I have which will be from s2_median_list. Now, for the places where there is a cloud mask has been applied, I want to do an average superimpose where there are no pixels. So, for the month of February, I need to do mosaicing from the month of January and March and put the value of their average pixels in the places where February has no pixel values.

I just started using Google Earth Engine this week and have very little idea on how to proceed with it.

2
  • 1
    You need to chose either python or javascript, since it is required to ask for one solution per question.
    – Padmanabha
    Commented Jun 7, 2023 at 19:45
  • You might want to take a look to the where function. Commented Jun 7, 2023 at 21:58

1 Answer 1

0

If you ensure you've got your monthly composites in an image collection and each have a timestamp, you can use a join to get images for the adjacent months.

var gapFilledCollection = ee.ImageCollection(
  ee.Join.saveAll('images')
    .apply({
      primary: monthlyCollection, 
      secondary: monthlyCollection, 
      condition: ee.Filter.maxDifference({
        difference: 31 * 24 * 60 * 60 * 1000,
        leftField: 'system:time_start',
        rightField: 'system:time_start'
      })
    })
    .map(function (image) {
      var fill = ee.ImageCollection(ee.List(
        image.get('images')
      )).mean()
      return ee.Image(image).unmask(fill)
    })
)

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

2
  • But, I have three images, one from January, the second from February, and the third one from March, and all of them are without the time stamp. Now, I want to do a composite for the month of February where the null pixels in that month can be superimposed by the average value of from the month of January and February if they exists for those null areas. If there are no value for January, then only the value of March will be taken for the remaining area and vice-versa.
    – Raj Bajaj
    Commented Jun 8, 2023 at 17:17
  • You can set a timestamp on your images. I assume you mistyped your comment and if February is masked out, you want the mean of January and March. That's what the script in my answer does. Commented Jun 8, 2023 at 19:39

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.