Note

This page was generated from a Jupyter notebook. To run and interact with it, you can download it here.

Sentinel-2 cloud/shadow-free composite

This example creates and downloads a large cloud/shadow free composite from a year of Sentinel-2 SR imagery. The area of interest is over Accra, Ghana. It demonstrates approaches for staying within Earth Engine computational and memory limits, and for compositing heavily clouded areas.

CLI commands equivalent to the API code snippets are given in the comments where possible.

Setup

geemap is required to run the notebook. You can uncomment the cell below to install it, if it isn’t installed already.

[1]:

# geemap should be installed if it isn't already.
# import sys
# if 'conda' in sys.prefix:
#     # install into the conda environment the notebook is being run from
#     !conda install --yes --prefix {sys.prefix} -c conda-forge geemap
# else:
#     # install into the python environment the notebook is being run from
#     !{sys.executable} -m pip install geemap

[2]:

import logging
import ee
import geedim as gd
import geemap.foliumap as geemap

# initialise earth engine with the high-volume endpoint
gd.Initialize()

Create and search a Sentinel-2 SR collection

In this step we search for images to composite, over the area and time period of interest.

Searching and compositing Sentinel-2 collections over large areas or time spans can be time-consuming. This is in part due to the resolution of the imagery, but also due to the computational burden of geedim cloud/shadow detection algorithms, where those are used. To reduce the number of images requiring cloud/shadow detection in the search, we specify a custom_filter that places an upper limit on the already existing CLOUDY_PIXEL_PERCENTAGE image granule property. This custom_filter is fast to apply and helps speed up the search.

While the custom_filter applies to properties of the image granule, the fill_portion and cloudless_portion filters apply to geedim calculated portions of the search region. Here they are used to limit the search results to relatively large and cloud/shadow-free image portions.

For consitency with the composite creation, we specify a cloud probability threshold of 40% (prob=40), and a cloud/shadow dilation distance of 100m (buffer=100).

[3]:

# geojson search polygon
region = {
    'type': 'Polygon', 'coordinates': [[
    [-0.37, 5.75], [-0.37, 5.5], [0.0, 5.5], [0.0, 5.75], [-0.37, 5.75]
  ]],
}

# create and search the Sentinel-2 TOA collection
coll = gd.MaskedCollection.from_name('COPERNICUS/S2_SR_HARMONIZED')
filt_coll = coll.search(
    start_date='2019-01-01', end_date='2020-01-01', region=region,
    cloudless_portion=75, fill_portion=30,
    custom_filter='CLOUDY_PIXEL_PERCENTAGE<25', prob=40, buffer=100,
)

# print the search results
print(f'Found {len(filt_coll.properties)} images.')
print(f'Image property descriptions:\n\n{filt_coll.schema_table}\n')
print(f'Search Results:\n\n{filt_coll.properties_table}')

# !geedim config --prob 40 --buffer 100 search -c s2-sr-hm -s 2019-01-01 -e 2020-01-01 --bbox -0.37 5.5 -0.0 5.75 -cf "CLOUDY_PIXEL_PERCENTAGE<25" -cp 75 -fp 30

Found 26 images.
Image property descriptions:

ABBREV     NAME                             DESCRIPTION
---------  -------------------------------  ----------------------------------------------
ID         system:id                        Earth Engine image id
DATE       system:time_start                Image capture date/time (UTC)
FILL       FILL_PORTION                     Portion of region pixels that are valid (%)
CLOUDLESS  CLOUDLESS_PORTION                Portion of filled pixels that are cloud/shadow
                                            free (%)
RADQ       RADIOMETRIC_QUALITY              Radiometric quality check
GEOMQ      GEOMETRIC_QUALITY                Geometric quality check
SAA        MEAN_SOLAR_AZIMUTH_ANGLE         Solar azimuth angle (deg)
SZA        MEAN_SOLAR_ZENITH_ANGLE          Solar zenith angle (deg)
VAA        MEAN_INCIDENCE_AZIMUTH_ANGLE_B1  View (B1) azimuth angle (deg)
VZA        MEAN_INCIDENCE_ZENITH_ANGLE_B1   View (B1) zenith angle (deg)

Search Results:

ID                                                                 DATE              FILL CLOUDLESS RADQ   GEOMQ     SAA   SZA    VAA   VZA
------------------------------------------------------------------ ---------------- ----- --------- ------ ------ ------ ----- ------ -----
COPERNICUS/S2_SR_HARMONIZED/20190109T101409_20190109T102253_T30NYM 2019-01-09 10:29 41.15     85.23 PASSED PASSED 139.50 37.20 110.79 11.21
COPERNICUS/S2_SR_HARMONIZED/20190122T102329_20190122T103839_T30NZM 2019-01-22 10:39 35.83     98.36 PASSED PASSED 139.22 33.99 291.77 11.04
COPERNICUS/S2_SR_HARMONIZED/20190122T102329_20190122T103839_T30NYM 2019-01-22 10:39 44.84     96.52 PASSED PASSED 138.15 34.57 284.86  7.01
COPERNICUS/S2_SR_HARMONIZED/20190129T101259_20190129T102648_T30NZM 2019-01-29 10:29 78.98     94.86 PASSED PASSED 133.66 34.76 103.87  7.45
COPERNICUS/S2_SR_HARMONIZED/20190129T101259_20190129T102648_T30NYM 2019-01-29 10:29 40.75     86.70 PASSED PASSED 132.71 35.41 110.77 11.22
COPERNICUS/S2_SR_HARMONIZED/20190203T101221_20190203T102328_T30NZM 2019-02-03 10:29 78.98     98.36 PASSED PASSED 131.54 34.00 102.95  7.46
COPERNICUS/S2_SR_HARMONIZED/20190203T101221_20190203T102328_T30NYM 2019-02-03 10:29 40.12     97.48 PASSED PASSED 130.60 34.67 110.22 11.22
COPERNICUS/S2_SR_HARMONIZED/20190211T102149_20190211T103444_T30NZM 2019-02-11 10:39 36.62     90.95 PASSED PASSED 130.61 30.63 291.60 11.03
COPERNICUS/S2_SR_HARMONIZED/20190211T102149_20190211T103444_T30NYM 2019-02-11 10:39 44.84     95.34 PASSED PASSED 129.57 31.32 284.89  7.01
COPERNICUS/S2_SR_HARMONIZED/20190221T102039_20190221T103744_T30NZM 2019-02-21 10:39 37.43     96.71 PASSED PASSED 125.14 28.46 291.63 11.03
COPERNICUS/S2_SR_HARMONIZED/20190221T102039_20190221T103744_T30NYM 2019-02-21 10:39 44.84     99.21 PASSED PASSED 124.14 29.19 284.91  6.97
COPERNICUS/S2_SR_HARMONIZED/20190226T102021_20190226T102750_T30NZM 2019-02-26 10:39 37.51     99.97 PASSED PASSED 122.02 27.33 292.19 11.05
COPERNICUS/S2_SR_HARMONIZED/20190226T102021_20190226T102750_T30NYM 2019-02-26 10:39 44.84     99.98 PASSED PASSED 121.05 28.09 285.84  6.97
COPERNICUS/S2_SR_HARMONIZED/20190303T102019_20190303T103501_T30NZM 2019-03-03 10:39 38.43     94.30 PASSED PASSED 118.64 26.18 291.45 11.01
COPERNICUS/S2_SR_HARMONIZED/20190305T101021_20190305T101942_T30NZM 2019-03-05 10:29 78.98     97.65 PASSED PASSED 114.79 27.96 102.84  7.49
COPERNICUS/S2_SR_HARMONIZED/20190305T101021_20190305T101942_T30NYM 2019-03-05 10:29 38.75     91.50 PASSED PASSED 114.03 28.77 110.18 11.24
COPERNICUS/S2_SR_HARMONIZED/20190422T102029_20190422T102825_T30NZM 2019-04-22 10:39 34.57     84.97 PASSED PASSED  70.43 20.33 291.93 11.07
COPERNICUS/S2_SR_HARMONIZED/20190422T102029_20190422T102825_T30NYM 2019-04-22 10:39 44.82     82.72 PASSED PASSED  71.12 21.16 284.90  7.05
COPERNICUS/S2_SR_HARMONIZED/20190502T102029_20190502T103323_T30NYM 2019-05-02 10:39 44.84     92.48 PASSED PASSED  62.60 21.88 285.00  7.07
COPERNICUS/S2_SR_HARMONIZED/20191029T102039_20191029T103707_T30NZM 2019-10-29 10:39 37.60     92.46 PASSED PASSED 140.78 24.88 291.61 11.05
COPERNICUS/S2_SR_HARMONIZED/20191029T102039_20191029T103707_T30NYM 2019-10-29 10:39 44.84     94.07 PASSED PASSED 139.25 25.45 284.87  6.98
COPERNICUS/S2_SR_HARMONIZED/20191203T102401_20191203T103859_T30NZM 2019-12-03 10:39 36.76     87.69 PASSED PASSED 149.20 32.68 292.18 11.04
COPERNICUS/S2_SR_HARMONIZED/20191203T102401_20191203T103859_T30NYM 2019-12-03 10:39 44.84     85.52 PASSED PASSED 147.94 33.15 285.92  6.99
COPERNICUS/S2_SR_HARMONIZED/20191223T102431_20191223T103653_T30NYM 2019-12-23 10:39 44.84     76.97 PASSED PASSED 146.26 35.50 285.85  6.93
COPERNICUS/S2_SR_HARMONIZED/20191225T101329_20191225T102600_T30NZM 2019-12-25 10:29 78.98     99.93 PASSED PASSED 143.96 36.51 103.79  7.49
COPERNICUS/S2_SR_HARMONIZED/20191225T101329_20191225T102600_T30NYM 2019-12-25 10:29 39.09     99.83 PASSED PASSED 142.92 37.05 110.76 11.25

Create composite image

Next we create a medoid composite image from the search result images. Again we specify a cloud probability threshold of 40% (prob=40), and a cloud/shadow dilation distance of 100m (buffer=100), to aggressively mask cloudy pixels. These parameters result in some masking of valid pixels (false positives), but with numerous component images, there are enough unmasked pixels to produce a full coverage composite.

Of the available compositing algorithms, the medoid method is the most memory and computation intensive. Computation scales roughly quadratically with the number of component images, so it is important to keep those to a necessary minimum. The number of images returned by the search (i.e. 26) is low enough to produce a reasonably fast composite in this case. Limiting the components to mostly cloud/shadow-free images, and aggressively masking these, helps improve the composite quality.

[4]:

medoid_im = filt_coll.composite('medoid', prob=40, buffer=100)

Visualise component and composite images

Now we use geemap to display the first 4 component images, their cloud/shadow masks, and the medoid composite. You can use the layer button in the top right hand corner to turn on/off the image layers. The composite will take a few minutes to load.

Note

You need to download and run the notebook to view the images.

[5]:

# geemap visualisation parameters
s2_vis_params = dict(min=0, max=4000, bands=['B4', 'B3', 'B2'], gamma=1.5)
mask_vis_params = dict(min=0, max=1)

map = geemap.Map()
map.centerObject(ee.Geometry(region), 11)

# iterate over the first 4 component image ID's
for im_id in list(filt_coll.properties.keys())[:4]:
    # create image and mask, and add to map
    im = gd.MaskedImage.from_id(im_id, prob=40, buffer=100)
    vis_image = im.ee_image.clip(region)
    cloudless_mask = vis_image.select('CLOUDLESS_MASK')

    map.addLayer(vis_image, s2_vis_params, f'Image {im_id[-38:]}')
    map.addLayer(cloudless_mask, mask_vis_params, f'Cl mask {im_id[-38:]}')

# add the composite to the map
map.addLayer(
    medoid_im.ee_image.clip(region), s2_vis_params, 'Medoid composite',
    # shown=False
)

region_im = ee.Image().byte().paint(
    featureCollection=ee.Geometry(region), width=2, color=1
)
map.addLayer(region_im, dict(palette=['FF0000']), 'Region')

map

[5]:

Make this Notebook Trusted to load map: File -> Trust Notebook

Download

Here we download the composite image directly. crs, scale and region parameters must be specified, as the composite has no fixed projection. Computation of the medoid composite is memory intensive, so we specify a max_tile_size of 16MB. Without this, the download would exceed the user memory limit.

[6]:

medoid_im.download(
    's2_medoid_im.tif', region=region, crs='EPSG:32735', scale=10,
    dtype='uint16', max_tile_size=16, overwrite=True,
)

# !geedim config --prob 40 --buffer 100 search -c s2-sr-hm -s 2019-01-01 -e 2020-01-01 --bbox -0.37 5.5 -0.0 5.75 -cf "CLOUDY_PIXEL_PERCENTAGE<25" -cp 75 -fp 30 composite -cm medoid download --crs EPSG:32735 --scale 10 --dtype uint16 --max-tile-size 16 -o

Export

Lastly we demonstrate an alternative to the direct download approach i.e. we export the composite image to an Earth Engine asset, and then download the computed asset image. Exporting is not subject to the user memory limit, and once computed, the asset image can be downloaded without the need for reducing max_tile_size. Bear in mind that exporting is usually slower than downloading, and exported GeoTIFFs can not be populated with STAC metadata.

[8]:

# 'geedim' is the name of a valid Earth Engine asset project
medoid_asset_id = f'projects/geedim/assets/s2_medoid_im'

# export to earth engine asset
medoid_task = medoid_im.export(
    medoid_asset_id, type='asset', region=region, crs='EPSG:32735',
    scale=10, dtype='uint16', wait=True,
)

# create a MaskedImage from the earth engine asset
medoid_asset_im = gd.MaskedImage.from_id(medoid_asset_id)
# download asset image
medoid_asset_im.download('s2_medoid_im.tif', overwrite=True)

# !geedim config --prob 40 --buffer 100 search -c s2-sr-hm -s 2019-01-01 -e 2020-01-01 --bbox -0.37 5.5 -0.0 5.75 -cf "CLOUDY_PIXEL_PERCENTAGE<25" -cp 75 -fp 30 composite -cm medoid export -t asset -f geedim --crs EPSG:32735 --scale 10 --dtype uint16 download -o

Preparing projects-geedim-assets-s2_medoid_im: done

There is no STAC entry for: projects/geedim/assets/s2_medoid_im

[ ]: