Access Sentinel 2 Data from AWS

https://registry.opendata.aws/sentinel-2-l2a-cogs/

import odc.ui
import yaml
from IPython.display import display
from odc.algo import to_rgba
from pystac_client import Client

from odc.stac import stac2ds, stac_load

cfg = """---
"*":
  warnings: ignore # Disable warnings about duplicate common names
sentinel-s2-l2a-cogs:
  assets:
    '*':
      data_type: uint16
      nodata: 0
      unit: '1'
    SCL:
      data_type: uint8
      nodata: 0
      unit: '1'
    visual:
      data_type: uint8
      nodata: 0
      unit: '1'
  aliases:  # Alias -> Canonical Name
    red: B04
    green: B03
    blue: B02
"""
cfg = yaml.load(cfg, Loader=yaml.SafeLoader)

catalog = Client.open("https://earth-search.aws.element84.com/v0")

Find STAC Items to Load

km2deg = 1.0 / 111
x, y = (113.887, -25.843)  # Center point of a query
r = 100 * km2deg

query = catalog.search(
    collections=["sentinel-s2-l2a-cogs"],
    datetime="2021-09-16",
    limit=10,
    bbox=(x - r, y - r, x + r, y + r),
)

items = list(query.get_items())
print(f"Found: {len(items):d} datasets")

Found: 9 datasets

Construct Dask Dataset

Note that even though there are 9 STAC Items on input, there is only one timeslice on output. This is because of groupy="solar_day". With that setting stac_load will place all items that occured on the same day (as adjusted for the timezone) into one image plane.

# Since we will plot it on a map we need to use `EPSG:3857` projection
crs = "epsg:3857"
zoom = 2 ** 5  # overview level 5

xx = stac_load(
    items,
    bands=("red", "green", "blue"),
    crs=crs,
    resolution=10 * zoom,
    chunks={},  # <-- use Dask
    groupby="solar_day",
    stac_cfg=cfg,
)
display(xx)

<xarray.Dataset>
Dimensions:      (time: 1, y: 1098, x: 833)
Coordinates:
  * time         (time) datetime64[ns] 2021-09-16T02:34:44
  * y            (y) float64 -2.797e+06 -2.798e+06 ... -3.148e+06 -3.148e+06
  * x            (x) float64 1.255e+07 1.255e+07 ... 1.282e+07 1.282e+07
    spatial_ref  int32 3857
Data variables:
    red          (time, y, x) uint16 dask.array<chunksize=(1, 1098, 833), meta=np.ndarray>
    green        (time, y, x) uint16 dask.array<chunksize=(1, 1098, 833), meta=np.ndarray>
    blue         (time, y, x) uint16 dask.array<chunksize=(1, 1098, 833), meta=np.ndarray>
Attributes:
    crs:           epsg:3857
    grid_mapping:  spatial_ref

Load data and convert to RGBA

%%time
rgba = to_rgba(xx, clamp=(1, 3000))
_rgba = rgba.compute()

CPU times: user 2.43 s, sys: 1.96 s, total: 4.39 s
Wall time: 25.7 s

Display Image on a map

dss = list(stac2ds(items, cfg))
_map = odc.ui.show_datasets(dss, style={"fillOpacity": 0.1}, scroll_wheel_zoom=True)
ovr = odc.ui.mk_image_overlay(_rgba)
_map.add_layer(ovr)
display(_map)

---