Cloud masking of Sentinel 2 using Google Earth Engine

Table of Content

Cloud masking of Sentinel 2 using Google Earth Engine

Source: https://github.com/samsammurphy/cloud-masking-sentinel2

Cloud masking of Sentinel 2 using Google Earth Engine

Cloud, and shadow, masking of Sentinel 2 images using Google Earth Engine Python API. Refactored from javascipt posted in the forum in this thread: Sentinel 2 cloud masking.

Dependencies

Usage

Authenticate Earth Engine (if necessary)

earthengine authenticate

Run the jupyter notebook

git clone https://github.com/samsammurphy/cloud-masking-sentinel2.git
cd cloud-masking-sentinel2
jupyter-notebook cloud-masking-sentinel2.ipynb

Cloud masking with Sentinel 2

Cloud and cloud shadow masking of Sentinel 2 images in Python. Refactored from javascipt taken from this thread: Sentinel 2 cloud masking

** Wavebands used for cloud and cloud shadow masking **
Band Use Wavelength (nm) Resolution (m)
B1 Aerosols 443 60
B2 Blue 490 10
B3 Green 560 10
B4 Red 665 10
B6 Red Edge 2 740 20
B8a Red Edge 4 865 20
B9 Water vapor 940 60
B10 Cirrus 1375 60
B11 SWIR 1 1610 20
B12 SWIR 2 2190 20
QA60 ESA Cloud n/a 60

NB: 8a (red4) used instead of 8 (nir) for narrower spectral width (i.e. 20 vs 115 nm)

In [ ]:
from IPython.display import display, Image
import math
import ee
ee.Initialize()
In [ ]:
def rescale(img, thresholds):
    """
    Linear stretch of image between two threshold values.
    """
    return img.subtract(thresholds[0]).divide(thresholds[1] - thresholds[0])
In [ ]:
def sentinelCloudScore(img):
    """
    Computes spectral indices of cloudyness and take the minimum of them.
    
    Each spectral index is fairly lenient because the group minimum 
    is a somewhat stringent comparison policy. side note -> this seems like a job for machine learning :)
    
    originally written by Matt Hancher for Landsat imagery
    adapted to Sentinel by Chris Hewig and Ian Housman
    """
    
    # cloud until proven otherwise
    score = ee.Image(1)

    # clouds are reasonably bright
    score = score.min(rescale(img.select(['blue']), [0.1, 0.5]))
    score = score.min(rescale(img.select(['aerosol']), [0.1, 0.3]))
    score = score.min(rescale(img.select(['aerosol']).add(img.select(['cirrus'])), [0.15, 0.2]))
    score = score.min(rescale(img.select(['red']).add(img.select(['green'])).add(img.select('blue')), [0.2, 0.8]))

    # clouds are moist
    ndmi = img.normalizedDifference(['red4','swir1'])
    score=score.min(rescale(ndmi, [-0.1, 0.1]))

    # clouds are not snow.
    ndsi = img.normalizedDifference(['green', 'swir1'])
    score=score.min(rescale(ndsi, [0.8, 0.6])).rename(['cloudScore'])
    
    return img.addBands(score)
In [ ]:
def ESAcloudMask(img):
    """
    European Space Agency (ESA) clouds from 'QA60', i.e. Quality Assessment band at 60m
     
    parsed by Nick Clinton
    """

    qa = img.select('QA60')

    # bits 10 and 11 are clouds and cirrus
    cloudBitMask = int(2**10)
    cirrusBitMask = int(2**11)

    # both flags set to zero indicates clear conditions.
    clear = qa.bitwiseAnd(cloudBitMask).eq(0).And(\
           qa.bitwiseAnd(cirrusBitMask).eq(0))
    
    # clouds is not clear
    cloud = clear.Not().rename(['ESA_clouds'])

    # return the masked and scaled data.
    return img.addBands(cloud)  
In [ ]:
def shadowMask(img,cloudMaskType):
    """
    Finds cloud shadows in images
    
    Originally by Gennadii Donchyts, adapted by Ian Housman
    """
    
    def potentialShadow(cloudHeight):
        """
        Finds potential shadow areas from array of cloud heights
        
        returns an image stack (i.e. list of images) 
        """
        cloudHeight = ee.Number(cloudHeight)
        
        # shadow vector length
        shadowVector = zenith.tan().multiply(cloudHeight)
        
        # x and y components of shadow vector length
        x = azimuth.cos().multiply(shadowVector).divide(nominalScale).round()
        y = azimuth.sin().multiply(shadowVector).divide(nominalScale).round()
        
        # affine translation of clouds
        cloudShift = cloudMask.changeProj(cloudMask.projection(), cloudMask.projection().translate(x, y)) # could incorporate shadow stretch?
        
        return cloudShift
  
    # select a cloud mask
    cloudMask = img.select(cloudMaskType)
    
    # make sure it is binary (i.e. apply threshold to cloud score)
    cloudScoreThreshold = 0.5
    cloudMask = cloudMask.gt(cloudScoreThreshold)

    # solar geometry (radians)
    azimuth = ee.Number(img.get('solar_azimuth')).multiply(math.pi).divide(180.0).add(ee.Number(0.5).multiply(math.pi))
    zenith  = ee.Number(0.5).multiply(math.pi ).subtract(ee.Number(img.get('solar_zenith')).multiply(math.pi).divide(180.0))

    # find potential shadow areas based on cloud and solar geometry
    nominalScale = cloudMask.projection().nominalScale()
    cloudHeights = ee.List.sequence(500,4000,500)        
    potentialShadowStack = cloudHeights.map(potentialShadow)
    potentialShadow = ee.ImageCollection.fromImages(potentialShadowStack).max()

    # shadows are not clouds
    potentialShadow = potentialShadow.And(cloudMask.Not())

    # (modified) dark pixel detection 
    darkPixels = toa.normalizedDifference(['green', 'swir2']).gt(0.25)

    # shadows are dark
    shadows = potentialShadow.And(darkPixels).rename(['shadows'])
    
    # might be scope for one last check here. Dark surfaces (e.g. water, basalt, etc.) cause shadow commission errors.
    # perhaps using a NDWI (e.g. green and nir)

    return img.addBands(shadows)
In [ ]:
def quicklook(bandNames, mn, mx, region, gamma=False, title=False):
    """
    Displays images in notebook
    """
    
    if title:
        print('\n',title)
 
    if not gamma:
        gamma = 1
        
    visual = Image(url=toa.select(bandNames).getThumbUrl({
                'region':region,
                'min':mn,
                'max':mx,
                'gamma':gamma,
                'title':title
                }))
    
    display(visual)
In [ ]:
# region of interest
geom = ee.Geometry.Point(-155.0844, 19.7189)

# start and end of time series
startDate = ee.Date('1980-01-01')
stopDate  = ee.Date('2020-01-01')
In [ ]:
# image collection
ic = ee.ImageCollection('COPERNICUS/S2')\
    .filterBounds(geom)\
    .filterDate(startDate,stopDate)
    
# single image
img = ee.Image(ic.first())
In [ ]:
# top of atmosphere reflectance
toa = img.select(['B1','B2','B3','B4','B6','B8A','B9','B10', 'B11','B12'],\
                 ['aerosol', 'blue', 'green', 'red', 'red2','red4','h2o', 'cirrus','swir1', 'swir2'])\
                 .divide(10000).addBands(img.select('QA60'))\
                 .set('solar_azimuth',img.get('MEAN_SOLAR_AZIMUTH_ANGLE'))\
                 .set('solar_zenith',img.get('MEAN_SOLAR_ZENITH_ANGLE'))
In [ ]:
# clouds
toa = sentinelCloudScore(toa)
toa = ESAcloudMask(toa)
In [ ]:
# cloud shadow
toa = shadowMask(toa,'cloudScore')
In [ ]:
# display region
region = geom.buffer(10000).bounds().getInfo()['coordinates']
In [ ]:
# quicklooks
quicklook(['red','green','blue'], 0, 0.25, region, gamma=1.5, title='RGB')
quicklook('cloudScore', 0, 1, region, title='Cloud Score')
quicklook('ESA_clouds', 0, 1, region, title = 'ESA Clouds (QA60)')
quicklook('shadows', 0, 1, region, title = 'Shadow mask')

0 thoughts on “Cloud masking of Sentinel 2 using Google Earth Engine”

  1. Hi, thank you for your example.. I was wondering how will the code work if I want a composite of a couple of months?

Leave a Reply

Your email address will not be published.