STAC + VirtualiZarr (Sentinel-2)¶
End-to-end example: query a STAC API for Sentinel COGs in California spanning
two adjacent UTM zones (10N and 11N), then merge them into a single grid using
merge(datafusion=True, virtualize=...).
Virtualization happens lazily — COGs are virtualized on-the-fly during
compute() via the default_virtualizer callback, so only sources that
intersect the target grid are processed and only when they don't already exist.
The UTM zone boundary runs roughly along -120° longitude through the Sierra Nevada. We pick a bbox around Lake Tahoe that straddles both zones.
AWS credentials: This notebook reads from the public e84-earth-search-sentinel-data bucket
(requester-pays) in us-west-2. You need valid AWS credentials configured.
In [1]:
Copied!
import warnings
import matplotlib.pyplot as plt
import xarray as xr
from affine import Affine
from rasterix import RasterIndex
from lazymerge.merge import merge
from lazymerge.virtualize import default_virtualizer
warnings.filterwarnings("ignore")
xr.set_options(display_expand_indexes=True)
import warnings
import matplotlib.pyplot as plt
import xarray as xr
from affine import Affine
from rasterix import RasterIndex
from lazymerge.merge import merge
from lazymerge.virtualize import default_virtualizer
warnings.filterwarnings("ignore")
xr.set_options(display_expand_indexes=True)
Out[1]:
<xarray.core.options.set_options at 0x11a5b6900>
1. Configure cloud storage¶
Sentinel-2 data lives in the e84-earth-search-sentinel-data S3 bucket
(requester-pays, us-west-2).
In [2]:
Copied!
import boto3
from obspec_utils.registry import ObjectStoreRegistry
from obstore.auth.boto3 import Boto3CredentialProvider
from obstore.store import S3Store
scheme = "s3://"
bucket = "e84-earth-search-sentinel-data"
region = "us-west-2"
session = boto3.Session()
object_store = S3Store(
bucket=bucket,
region=region,
request_payer=True,
credential_provider=Boto3CredentialProvider(session=session),
)
registry = ObjectStoreRegistry({f"{scheme}{bucket}": object_store})
import boto3
from obspec_utils.registry import ObjectStoreRegistry
from obstore.auth.boto3 import Boto3CredentialProvider
from obstore.store import S3Store
scheme = "s3://"
bucket = "e84-earth-search-sentinel-data"
region = "us-west-2"
session = boto3.Session()
object_store = S3Store(
bucket=bucket,
region=region,
request_payer=True,
credential_provider=Boto3CredentialProvider(session=session),
)
registry = ObjectStoreRegistry({f"{scheme}{bucket}": object_store})
2. Set up Icechunk repository¶
Icechunk provides a versioned zarr store that supports virtual chunk references back to the original COG files on S3.
In [3]:
Copied!
import icechunk
location = "data/icechunk_california"
storage = icechunk.local_filesystem_storage(location)
config = icechunk.RepositoryConfig.default()
config.caching = icechunk.CachingConfig(
num_bytes_chunks=500_000_000, # 500 MB chunk data cache
)
s3_chunk_store = icechunk.s3_store(
region=region,
requester_pays=True,
)
config.set_virtual_chunk_container(
icechunk.VirtualChunkContainer(f"{scheme}{bucket}/", s3_chunk_store),
)
credentials = icechunk.containers_credentials(
{f"{scheme}{bucket}/": icechunk.s3_credentials(from_env=True)},
)
repo = icechunk.Repository.open_or_create(
storage=storage,
config=config,
authorize_virtual_chunk_access=credentials,
)
ic_session = repo.writable_session("main")
import icechunk
location = "data/icechunk_california"
storage = icechunk.local_filesystem_storage(location)
config = icechunk.RepositoryConfig.default()
config.caching = icechunk.CachingConfig(
num_bytes_chunks=500_000_000, # 500 MB chunk data cache
)
s3_chunk_store = icechunk.s3_store(
region=region,
requester_pays=True,
)
config.set_virtual_chunk_container(
icechunk.VirtualChunkContainer(f"{scheme}{bucket}/", s3_chunk_store),
)
credentials = icechunk.containers_credentials(
{f"{scheme}{bucket}/": icechunk.s3_credentials(from_env=True)},
)
repo = icechunk.Repository.open_or_create(
storage=storage,
config=config,
authorize_virtual_chunk_access=credentials,
)
ic_session = repo.writable_session("main")
2026-05-29T19:00:35.363936Z WARN icechunk::storage::object_store: The LocalFileSystem storage is not safe for concurrent commits. If more than one thread/process will attempt to commit at the same time, prefer using object stores. at icechunk/src/storage/object_store.rs:92
3. Ingest STAC search results¶
Search Sentinel scenes and ingest the metadata with the red, green, and blue asset hrefs.
In [4]:
Copied!
from zarr_datafusion_search import ingest_stac_search
rows = await ingest_stac_search(
url="https://earth-search.aws.element84.com/v1",
session=ic_session,
bbox=(-121.5, 38.5, -118.5, 40.5),
datetime="2025-10-01/2025-10-15",
collections=["sentinel-2-c1-l2a"],
max_items=10,
asset_hrefs=["red", "green", "blue"],
query={"eo:cloud_cover": {"lt": 10}},
)
rows
from zarr_datafusion_search import ingest_stac_search
rows = await ingest_stac_search(
url="https://earth-search.aws.element84.com/v1",
session=ic_session,
bbox=(-121.5, 38.5, -118.5, 40.5),
datetime="2025-10-01/2025-10-15",
collections=["sentinel-2-c1-l2a"],
max_items=10,
asset_hrefs=["red", "green", "blue"],
query={"eo:cloud_cover": {"lt": 10}},
)
rows
Out[4]:
10
4. Query metadata with DataFusion¶
Use DataFusion SQL to inspect the ingested metadata and verify we have scenes in both UTM zones.
In [5]:
Copied!
from datafusion import SessionContext
from geodatafusion import register_all
from zarr_datafusion_search import ZarrTable
ic_session = repo.writable_session("main")
zarr_table = await ZarrTable.from_icechunk(session=ic_session, group_path="/meta")
ctx = SessionContext()
register_all(ctx)
ctx.register_table_provider("meta", zarr_table)
df = ctx.sql(
"SELECT * FROM meta ORDER BY datetime",
)
df
from datafusion import SessionContext
from geodatafusion import register_all
from zarr_datafusion_search import ZarrTable
ic_session = repo.writable_session("main")
zarr_table = await ZarrTable.from_icechunk(session=ic_session, group_path="/meta")
ctx = SessionContext()
register_all(ctx)
ctx.register_table_provider("meta", zarr_table)
df = ctx.sql(
"SELECT * FROM meta ORDER BY datetime",
)
df
Out[5]:
| asset_blue | asset_green | asset_red | bbox | constellation | created | datetime | earthsearch:payload_id | eo:cloud_cover | grid:code | id | mgrs:grid_square | mgrs:latitude_band | mgrs:utm_zone | platform | proj:epsg | s2:cloud_shadow_percentage | s2:datastrip_id | s2:datatake_id | s2:datatake_type | s2:degraded_msi_data_percentage | s2:generation_time | s2:high_proba_clouds_percentage | s2:medium_proba_clouds_percentage | s2:nodata_pixel_percentage | s2:not_vegetated_percentage | s2:processing_baseline | s2:product_type | s2:product_uri | s2:reflectance_conversion_factor | s2:saturated_defective_pixel_percentage | s2:snow_ice_percentage | s2:thin_cirrus_percentage | s2:tile_id | s2:unclassified_percentage | s2:vegetation_percentage | s2:water_percentage | storage:platform | storage:region | storage:requester_pays | updated | view:azimuth | view:incidence_angle | view:sun_azimuth | view:sun_elevation |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| sentinel-2 | 0.234831 | MGRS-11TKE | KE | T | 11 | sentinel-2c | 32611 | 0.003685 | INS-NOBS | 0.0134 | 0.049996 | 0.106056 | 14.551552 | 90.295559 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 3.5e-05 | 0.078779 | 0.044552 | 4.145623 | 4.957866 | AWS | us-west-2 | False | 103.39706481599862 | 8.226388589336729 | 162.817826061657 | 40.6701397768753 | ||||||||||||||
| sentinel-2 | 0.279838 | MGRS-10TGK | GK | T | 10 | sentinel-2c | 32610 | 0.003699 | INS-NOBS | 0.0143 | 0.059549 | 0.129535 | 25.285274 | 89.158958 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 3.6e-05 | 0.090754 | 0.04839 | 4.358395 | 5.647048 | AWS | us-west-2 | False | 104.27783984543423 | 8.624376857205261 | 162.649924254137 | 40.6431223664525 | ||||||||||||||
| sentinel-2 | 0.060511 | MGRS-11SLD | LD | S | 11 | sentinel-2c | 32611 | 3e-06 | INS-NOBS | 0.0144 | 0.002943 | 0.034141 | 0.0 | 97.190738 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.021695 | 0.023427 | 0.043334 | 1.820538 | 0.748289 | AWS | us-west-2 | False | 200.10979481301962 | 2.892353807621798 | 164.195726992237 | 41.7824106469429 | ||||||||||||||
| sentinel-2 | 0.015442 | MGRS-11SKD | KD | S | 11 | sentinel-2c | 32611 | 0.0 | INS-NOBS | 0.0211 | 0.000392 | 0.004303 | 0.21976 | 63.256186 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.00286 | 0.010747 | 0.017371 | 31.539431 | 4.918656 | AWS | us-west-2 | False | 104.90327192418415 | 6.89134039053623 | 162.710280334676 | 41.5515224355808 | ||||||||||||||
| sentinel-2 | 0.016039 | MGRS-10SGJ | GJ | S | 10 | sentinel-2c | 32610 | 0.0 | INS-NOBS | 0.0155 | 0.000293 | 0.003261 | 7.203414 | 52.777588 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.002195 | 0.012485 | 0.02053 | 41.350394 | 5.318559 | AWS | us-west-2 | False | 103.21261451408677 | 7.85596235915169 | 162.444003721581 | 41.5063913718383 | ||||||||||||||
| sentinel-2 | 0.012961 | MGRS-10SFJ | FJ | S | 10 | sentinel-2c | 32610 | 0.0 | INS-NOBS | 0.0 | 0.0 | 0.00036 | 96.313667 | 16.721104 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.0 | 0.012601 | 0.00864 | 79.38875 | 0.330134 | AWS | us-west-2 | False | 108.65371859143146 | 11.50286086589293 | 160.982084293749 | 41.2051032455647 | ||||||||||||||
| sentinel-2 | 0.079038 | MGRS-11SLC | LC | S | 11 | sentinel-2c | 32611 | 0.000574 | INS-NOBS | 0.0107 | 0.005474 | 0.04783 | 0.0 | 94.780922 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.01857 | 0.025733 | 0.048261 | 2.165311 | 2.67311 | AWS | us-west-2 | False | 243.40666102784363 | 3.8419365247390145 | 164.078898419319 | 42.6663965710168 | ||||||||||||||
| sentinel-2 | 0.011456 | MGRS-11SKC | KC | S | 11 | sentinel-2c | 32611 | 0.00065 | INS-NOBS | 0.0214 | 0.002538 | 0.006576 | 0.0 | 51.878589 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.111386 | 0.002342 | 0.003105 | 45.780191 | 0.745844 | AWS | us-west-2 | False | 102.1427886624837 | 4.93699741381131 | 162.593209564776 | 42.4331142678946 | ||||||||||||||
| sentinel-2 | 0.011798 | MGRS-10SGH | GH | S | 10 | sentinel-2c | 32610 | 0.000109 | INS-NOBS | 0.0211 | 0.000345 | 0.00288 | 0.0 | 41.815534 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.033925 | 0.008573 | 0.008457 | 55.799222 | 0.983388 | AWS | us-west-2 | False | 105.22571248878553 | 6.72451497425614 | 162.230325800081 | 42.3698680625101 | ||||||||||||||
| sentinel-2 | 0.005109 | MGRS-10SFH | FH | S | 10 | sentinel-2c | 32610 | 0.0 | INS-NOBS | 0.0 | 0.000334 | 0.001241 | 79.153657 | 55.602968 | 05.11 | S2MSI2A | 1.00270392241554 | 0 | 0.0 | 0.003533 | 0.008913 | 42.907405 | 1.171952 | AWS | us-west-2 | False | 107.31339183787854 | 10.846582766031384 | 160.769769734024 | 42.0662180365237 |
5. Create virtualizer¶
Create a default_virtualizer callback that will lazily virtualize COGs
during merge. Pass target_resolution to automatically select the best COG
overview level for the output grid — only level 0 (full resolution) and the
optimal overview are virtualized, avoiding unnecessary work.
Virtualization is deferred — nothing happens until compute() is called,
and only sources that intersect each target chunk are virtualized.
This virtualization step only happens once so subsequent calls will only need to read data and won't need to virtualize sources and will be significantly faster.
In [6]:
Copied!
virtualizer = default_virtualizer(registry, target_resolution=20.0)
virtualizer = default_virtualizer(registry, target_resolution=20.0)
6. Lazy merge¶
Define a target grid in UTM 11N (EPSG:32611) covering the Lake Tahoe area. Sources in UTM 10N will be reprojected on the fly during merge.
The virtualize callback triggers lazy VirtualiZarr creation — each source
COG is virtualized into the Icechunk store the first time a chunk needs it.
Coordination ensures each source is virtualized exactly once, even when
multiple chunks reference the same source concurrently.
In [7]:
Copied!
bands = ["red", "green", "blue"]
result_arr, result_spatial, result_proj, time_coords = merge(
store=ic_session.store,
crs="EPSG:32611",
bbox=(200000.0, 4260000.0, 340000.0, 4490000.0),
resolution=20.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
temporal_grouping="P1D",
nodata=0,
virtualize=virtualizer,
)
result_arr
bands = ["red", "green", "blue"]
result_arr, result_spatial, result_proj, time_coords = merge(
store=ic_session.store,
crs="EPSG:32611",
bbox=(200000.0, 4260000.0, 340000.0, 4490000.0),
resolution=20.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
temporal_grouping="P1D",
nodata=0,
virtualize=virtualizer,
)
result_arr
Out[7]:
|
7. Open the lazy Cubed array with xarray¶
In [9]:
Copied!
da = xr.DataArray(result_arr, dims=["time", "bands", "y", "x"])
affine = Affine(*result_spatial.transform)
raster_idx = RasterIndex.from_transform(
affine=affine,
width=da.sizes["x"],
height=da.sizes["y"],
x_dim="x",
y_dim="y",
crs=result_proj.code,
)
coords = xr.Coordinates.from_xindex(raster_idx)
da = da.assign_coords(coords)
da = da.proj.assign_crs(spatial_ref=result_proj.code, allow_override=True)
da = da.assign_coords({"time": time_coords, "bands": bands})
da
da = xr.DataArray(result_arr, dims=["time", "bands", "y", "x"])
affine = Affine(*result_spatial.transform)
raster_idx = RasterIndex.from_transform(
affine=affine,
width=da.sizes["x"],
height=da.sizes["y"],
x_dim="x",
y_dim="y",
crs=result_proj.code,
)
coords = xr.Coordinates.from_xindex(raster_idx)
da = da.assign_coords(coords)
da = da.proj.assign_crs(spatial_ref=result_proj.code, allow_override=True)
da = da.assign_coords({"time": time_coords, "bands": bands})
da
Out[9]:
<xarray.DataArray 'array-004' (time: 1, bands: 3, y: 11500, x: 7000)> Size: 966MB
cubed.Array<array-004, shape=(1, 3, 11500, 7000), dtype=float32, chunks=((1,), (1, 1, 1), (1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 236), (1024, 1024, 1024, 1024, 1024, 1024, 856))>
Coordinates:
* time (time) datetime64[s] 8B 2025-10-12
* bands (bands) <U5 60B 'red' 'green' 'blue'
* y (y) float64 92kB 4.49e+06 4.49e+06 ... 4.26e+06 4.26e+06
* x (x) float64 56kB 2e+05 2e+05 2e+05 ... 3.4e+05 3.4e+05 3.4e+05
* spatial_ref int64 8B 0
Indexes:
┌ x RasterIndex (crs=EPSG:32611)
└ y
spatial_ref CRSIndex (crs=EPSG:32611)8. Compute and visualize¶
Compute to materialize the array and then plot it. This is where virtualization actually happens — each chunk task virtualizes any sources it needs before reading pixel data.
In [ ]:
Copied!
import logging
logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
import logging
logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
In [10]:
Copied!
da = da.compute()
da
da = da.compute()
da
INFO:lazymerge.virtualize:Virtualized S2C_T11TKE_20251012T184845_L2A (3 bands, levels [0], native_res=10.1) INFO:lazymerge.virtualize:Virtualized S2C_T10TGK_20251012T184845_L2A (3 bands, levels [0, 1], native_res=9.2) INFO:lazymerge.virtualize:Virtualized S2C_T11SLD_20251012T184845_L2A (3 bands, levels [0], native_res=10.4) INFO:lazymerge.virtualize:Virtualized S2C_T11SKD_20251012T184845_L2A (3 bands, levels [0], native_res=10.6) INFO:lazymerge.virtualize:Virtualized S2C_T10SGJ_20251012T184845_L2A (3 bands, levels [0], native_res=10.7) INFO:lazymerge.virtualize:Virtualized S2C_T11SKC_20251012T184845_L2A (3 bands, levels [0], native_res=10.6) INFO:lazymerge.virtualize:Virtualized S2C_T11SLC_20251012T184845_L2A (3 bands, levels [0], native_res=10.4) INFO:lazymerge.virtualize:Virtualized S2C_T10SGH_20251012T184845_L2A (3 bands, levels [0], native_res=10.6)
Out[10]:
<xarray.DataArray 'array-004' (time: 1, bands: 3, y: 11500, x: 7000)> Size: 966MB
array([[[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[1271., 1211., 1138., ..., 2332., 2362., 2062.],
[1505., 1561., 1194., ..., 1794., 1809., 1902.],
[1490., 1418., 1238., ..., 2128., 2092., 1834.]],
[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[1393., 1301., 1160., ..., 2246., 2180., 1974.],
[1512., 1540., 1304., ..., 1832., 1746., 1852.],
[1504., 1386., 1297., ..., 2018., 1966., 1787.]],
[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[1112., 1102., 1045., ..., 1905., 1935., 1750.],
[1225., 1232., 1081., ..., 1567., 1631., 1644.],
[1244., 1186., 1086., ..., 1740., 1700., 1578.]]]],
shape=(1, 3, 11500, 7000), dtype=float32)
Coordinates:
* time (time) datetime64[s] 8B 2025-10-12
* bands (bands) <U5 60B 'red' 'green' 'blue'
* y (y) float64 92kB 4.49e+06 4.49e+06 ... 4.26e+06 4.26e+06
* x (x) float64 56kB 2e+05 2e+05 2e+05 ... 3.4e+05 3.4e+05 3.4e+05
* spatial_ref int64 8B 0
Indexes:
┌ x RasterIndex (crs=EPSG:32611)
└ y
spatial_ref CRSIndex (crs=EPSG:32611)In [11]:
Copied!
da = da.sel(time="2025-10-12")
width = 8
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.tight_layout()
da = da.sel(time="2025-10-12")
width = 8
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.tight_layout()
