Composite Multi-Product Imagery

Composite imagery from two data sources and display as a single image.

from descarteslabs.catalog import Image, properties as p
from descarteslabs.utils import display
import numpy as np

# Define a bounding box around Taos in a GeoJSON

taos = {
    "type": "Polygon",
    "coordinates": [
        [
            [-105.71868896484375, 36.33725319397006],
            [-105.2105712890625, 36.33725319397006],
            [-105.2105712890625, 36.73668306473141],
            [-105.71868896484375, 36.73668306473141],
            [-105.71868896484375, 36.33725319397006],
        ]
    ],
}

# Create an ImageCollection
search = (
    Image.search()
    .intersects(taos)
    .filter(
        p.product_id.any_of(["usgs:landsat:oli-tirs:c2:l2:v0", "esa:sentinel-2:l2a:v1"])
    )
    .filter("2018-05-01" <= p.acquired < "2018-06-01")
    .filter(p.cloud_fraction < 0.2)
    .sort("acquired")
    .limit(15)
)
images = search.collect()

See which images we have, and how many per product:

print(images)

ImageCollection of 7 images
  * Dates: May 10, 2018 to May 29, 2018
  * Products: usgs:landsat:oli-tirs:c2:l2:v0: 3, esa:sentinel-2:l2a:v1: 4

And if you’re curious, which image IDs:

print(images.each.id)

'usgs:landsat:oli-tirs:c2:l2:v0:LC08_L2SP_034034_20180510_20200901_02_T1'
'esa:sentinel-2:l2a:v1:S2A_MSIL2A_20180514T173911_N0001_R098_T13SDA_20200312T042108'
'usgs:landsat:oli-tirs:c2:l2:v0:LC08_L2SP_033035_20180519_20200901_02_T1'
'esa:sentinel-2:l2a:v1:S2A_MSIL2A_20180524T174051_N0001_R098_T13SDA_20200311T040700'
'usgs:landsat:oli-tirs:c2:l2:v0:LC08_L2SP_032035_20180528_20200901_02_T1'
'esa:sentinel-2:l2a:v1:S2B_MSIL2A_20180529T173859_N0001_R098_T13SDA_20200312T074340'
'esa:sentinel-2:l2a:v1:S2B_MSIL2A_20180529T173859_N0500_R098_T13SDA_20230830T051704'

Make a median composite of the images.

# Request a stack of all the images using the same GeoContext with lower resolution
arr_stack = images.stack("red green blue", resolution=60, data_type="Float64")

# Composite the images based on the median pixel value
composite = np.ma.median(arr_stack, axis=0)
display(composite, title="Taos Composite", size=2)