Create a Sentinel-2 RGB image

from shapely import box
from tilebox.datasets import Client, field

# west, south, east, north
area = box(14.42, -23.43, 14.58, -23.25)  # Sandwich Harbour, Namibia

collection = Client().dataset("open_data.aws_earth.sentinel2").collection("L2A")
scenes = collection.query(
    temporal_extent=("2024-06-17", "2024-06-18"),
    spatial_extent=area,
    filter=field("cloud_cover") < 1,
)

datapoint = scenes.sortby("cloud_cover").isel(time=0)
print(datapoint.stac_id.item(), datapoint.cloud_cover.item())

from tilebox.datasets.assets import AssetCollection

assets = AssetCollection.from_datapoint(datapoint)

import asyncio

import numpy as np
from tilebox.storage.aio import Client as StorageClient
from tilebox.storage.geotiff import window_from_bounds

async def read_rgb():
    storage = StorageClient()

    async def read_band(key):
        asset = assets[key]
        geotiff = await storage.open_geotiff(asset)
        window = window_from_bounds(geotiff, area.bounds, crs="EPSG:4326")
        raster = await geotiff.read(window=window)
        pixels = raster.data[0].astype(np.float32)
        return pixels * asset.raster.scale + asset.raster.offset

    red, green, blue = await asyncio.gather(
        read_band("red"),
        read_band("green"),
        read_band("blue"),
    )
    return np.stack((red, green, blue), axis=-1)

rgb = asyncio.run(read_rgb())

from PIL import Image, ImageOps

display_rgb = np.power(np.clip(rgb / 0.3, 0, 1), 1 / 2.2)
image = Image.fromarray((display_rgb * 255).astype(np.uint8))
image = ImageOps.autocontrast(image, cutoff=0.5)
image.save("sentinel2-sandwich-harbour.png")