0

I need a time series of Sentinel-2 images freed from clouds and with the masked pixels of each image replaced by interpolated values (linear over time). Thus, first, I applied the cloud mask as reported in Sentinel-2 Cloud Masking with s2cloudless. Subsequently, I tried to translate (from JavaScript to Python) the code reported in Temporal Gap-Filling with Linear Interpolation in GEE (first example) in which a linear interpolation is applied on a time series in order to fill the cloud gaps. The cloud masking works well: the output is an image collection formed by images in which the cloudy pixels are NaN. However, I wasn't able to apply the second step, namely the linear interpolation so that null pixels (clouds) were replaced by time-interpolated values (from the previous and subsequent images of the time series). Below the code:

'''CLOUD MASK CODE'''
startyear = 2019
endyear = 2020
years = range(startyear,endyear)
startdate = ee.Date.fromYMD(startyear,6,1)
enddate = ee.Date.fromYMD(endyear,5,31)

bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12']

CLOUD_FILTER = 80 
CLD_PRB_THRESH = 10 
NIR_DRK_THRESH = 0.15 
CLD_PRJ_DIST = 0.5   
BUFFER = 50

def S2_collection_cloud(aoi, startdate, enddate):
    s2 = ee.ImageCollection('COPERNICUS/S2_SR').filterBounds(aoi).filterDate(startdate, enddate).filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', CLOUD_FILTER)).map(lambda image: image.clip(area))
    s2_cloudless = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY').filterBounds(aoi).filterDate(startdate, enddate).map(lambda image: image.clip(area))
    return ee.ImageCollection(ee.Join.saveFirst('s2cloudless').apply(**{ 
        'primary': s2, 'secondary': s2_cloudless, 'condition': ee.Filter.equals(**{
            'leftField': 'system:index', 'rightField': 'system:index'})}))

def add_cloud_bands(img):
    cloud_prb = ee.Image(img.get('s2cloudless')).select('probability')
    is_cloud = cloud_prb.gt(CLD_PRB_THRESH).rename('clouds')
    return img.addBands(ee.Image([cloud_prb, is_cloud]))

def add_shadow_bands(img):
    not_water = img.select('SCL').neq(6)
    SR_BAND_SCALE = 1e4  # =10000
    dark_pixels = img.select('B8').lt(NIR_DRK_THRESH*SR_BAND_SCALE).multiply(not_water).rename('dark_pixels')
    shadow_azimuth = ee.Number(90).subtract(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')));
    cloud_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'))
    shadows = cloud_proj.multiply(dark_pixels).rename('shadows')
    return img.addBands(ee.Image([dark_pixels, cloud_proj, shadows]))

def add_cld_shdw_mask(img):
    img_cloud = add_cloud_bands(img)
    img_cloud_shadow = add_shadow_bands(img_cloud)
    is_cld_shdw = img_cloud_shadow.select('clouds').add(img_cloud_shadow.select('shadows')).gt(0)
    is_cld_shdw = (is_cld_shdw.focal_min(2).focal_max(BUFFER*2/20)
        .reproject(**{'crs': img.select([0]).projection(), 'scale': 20})
        .rename('cloudmask'))
    return img_cloud_shadow.addBands(is_cld_shdw)

def apply_cld_shdw_mask(img):
    not_cld_shdw = img.select('cloudmask').Not()
    return img.select('B.*').updateMask(not_cld_shdw)

S2_col_cloudmask = S2_collection_cloud(area, startdate, enddate) #the area was set previously
S2_col_cloudmask_disp = S2_col_cloudmask.map(add_cld_shdw_mask)
S2_col_nocloud= S2_col_cloudmask_disp.map(apply_cld_shdw_mask) 


'''START THE INTERPOLATION TEMPORAL INTERPOLATION CODE'''

def add_time_filter_band(image):
    timeImage = image.metadata('system:time_start').rename('timestamp')
    timeImageMasked = timeImage.updateMask(image.mask().select(0))
    return image.addBands(timeImageMasked)

S2_col_nocloud = S2_col_nocloud.map(add_time_filter_band)

days = 30
millis = ee.Number(days).multiply(1000*60*60*24)

maxDiffFilter = ee.Filter.maxDifference(difference= days,
                                         leftField= 'system:time_start',
                                         rightField= 'system:time_start')

lessEqFilter = ee.Filter.lessThanOrEquals(leftField= 'system:time_start',
                                           rightField= 'system:time_start')

greaterEqFilter = ee.Filter.greaterThanOrEquals(leftField= 'system:time_start',
                                                 rightField= 'system:time_start')

filter1 = ee.Filter.And(maxDiffFilter, lessEqFilter)

join1 = ee.Join.saveAll(matchesKey= 'after',
                        ordering= 'system:time_start',
                        ascending=False)

join1Result = join1.apply(primary= collections_list[0],
                          secondary= collections_list[0],
                          condition= filter1)

filter2 = ee.Filter.And(maxDiffFilter, greaterEqFilter)

join2 = ee.Join.saveAll(matchesKey= 'before',
                        ordering= 'system:time_start',
                        ascending = True)

join2Result = join2.apply(primary= join1Result,
                          secondary= join1Result,
                          condition= filter2)

def interpolation(image):
    image = ee.Image(image)
    beforeImages = ee.List(image.get('before'))
    beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
    afterImages = ee.List(image.get('after'))
    afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
    t1 = beforeMosaic.select('timestamp').rename('t1')
    t2 = afterMosaic.select('timestamp').rename('t2')
    t = image.metadata('system:time_start').rename('t')
    timeImage = ee.Image.cat([t1, t2, t])
    timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {'t': timeImage.select('t'),
                                                              't1': timeImage.select('t1'),
                                                              't2': timeImage.select('t2')})
    interpolated = beforeMosaic.add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
    result = image.unmask(interpolated)
    return image.unmask(interpolated).copyProperties(image, ['system:time_start'])

interpolatedCol = ee.ImageCollection(join2Result.map(interpolation))

I would have expected a final image collection (interpolatedCol) with all the images whose masked pixels (clouds) were filled with interpolated values, but this is not the case! Please help! Any other solution that allows me time-interpolation after using the S2cloudless product is welcome!

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy