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.
where
function.