Note

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

Landsat-7 cloud/shadow-free composite

This example aims to create a Landsat-7 cloud/shadow free composite on, or as close as possible to 22-23 November 2016. The area of interest covers a range of natural, agricultural and urban areas around Stellenbosch, South Africa.

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 ee
import geedim as gd
import geemap.foliumap as geemap

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

Create and search a Landsat-7 collection

[3]:

# geojson search polygon
region = {
    'type': 'Polygon', 'coordinates': [[
        (19.075, -34.115), (19.075, -33.731), (18.723, -33.731),
        (18.723, -34.115), (19.075, -34.115)
    ]]
}

# create and search the Landsat-7 collection
coll = gd.MaskedCollection.from_name('LANDSAT/LE07/C02/T1_L2')
filt_coll = coll.search(
    '2016-11-01', '2016-12-19', region, cloudless_portion=40
)

# print the search results
print(filt_coll.schema_table, end='\n\n')
print(filt_coll.properties_table)

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 (%)
GRMSE      GEOMETRIC_RMSE_MODEL  Orthorectification RMSE (m)
SAA        SUN_AZIMUTH           Solar azimuth angle (deg)
SEA        SUN_ELEVATION         Solar elevation angle (deg)

ID                                          DATE              FILL CLOUDLESS GRMSE   SAA   SEA
------------------------------------------- ---------------- ----- --------- ----- ----- -----
LANDSAT/LE07/C02/T1_L2/LE07_175083_20161116 2016-11-16 08:37 65.43     80.93  4.89 67.06 61.31
LANDSAT/LE07/C02/T1_L2/LE07_175084_20161116 2016-11-16 08:38 67.79     72.44  6.74 65.06 60.50
LANDSAT/LE07/C02/T1_L2/LE07_175083_20161202 2016-12-02 08:37 65.36     99.99  4.90 73.95 62.07
LANDSAT/LE07/C02/T1_L2/LE07_175084_20161202 2016-12-02 08:38 67.02     99.98  7.48 71.70 61.41
LANDSAT/LE07/C02/T1_L2/LE07_175083_20161218 2016-12-18 08:37 66.27    100.00  4.95 78.16 61.24
LANDSAT/LE07/C02/T1_L2/LE07_175084_20161218 2016-12-18 08:38 67.05     99.98  7.13 75.91 60.67

Notes on search results

The 2016-11-16 images are closest to the target dates, but have some cloud in them.
No single image has full coverage (FILL=100) of the search area. In part, this is just due to the footprint of the images, but is also a result of the Landsat-7 SLC failure.

Find composite images

Here, we find cloud/shadow-free mosaic, and q-mosaic composite images, prioritising images closest to 2016-11-22 by specifying the date parameter.

[4]:

mosaic_im = filt_coll.composite(
    method=gd.CompositeMethod.mosaic, date='2016-11-22'
)
q_mosaic_im = filt_coll.composite(
    method=gd.CompositeMethod.q_mosaic, date='2016-11-22'
)

Visualise search result and composite images

You can select which images to show/hide by clicking the layer button on the top right.

Note

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

[5]:

l7_vis_params = dict(min=7300, max=13000, bands=[
                     'SR_B3', 'SR_B2', 'SR_B1'], gamma=1.5)
map = geemap.Map()

map.centerObject(ee.Geometry(region), 11)
for im_id in filt_coll.properties.keys():
    im = gd.MaskedImage.from_id(im_id, mask=False)
    map.addLayer(im.ee_image.clip(region), l7_vis_params, im_id[-20:])

map.addLayer(
    mosaic_im.ee_image.clip(region), l7_vis_params, 'Mosaic composite'
)
map.addLayer(
    q_mosaic_im.ee_image.clip(region), l7_vis_params, 'Q-mosaic composite'
)

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

Visualisation notes

The mosaic method image contains some artefacts due to remnant cloud in the masked input images.
By prioritising pixels with the highest distance to cloud, the q-mosaic method is robust to imperfect cloud/shadow masking, and produces a composite free of cloud artefacts.

Download

Lastly, we download the q-mosaic composite. crs, scale and region parameters must be specified, as the composite has no fixed projection.

[6]:

# download the q_mosaic composite image, specifying crs, scale and region as
# it has no fixed projection
q_mosaic_im.download(
    'l7_q_mosaic_im.tif', crs='EPSG:3857', scale=30, region=region,
    dtype='uint16', overwrite=True
)