Do an average mosaic on an image in Google Earth Engine

Question

So 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 less idea on how to proceed with it.

Thanks.

Answer 1

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