STAC + VirtualiZarr¶
End-to-end example: query a STAC API for Landsat 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.
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 usgs-landsat 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
import zarr
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
import zarr
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 0x125ab6510>
1. Configure cloud storage¶
Landsat Collection 2 Level-2 data lives in the usgs-landsat 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 = "usgs-landsat"
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 = "usgs-landsat"
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 [7]:
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-29T18:57:02.297127Z 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¶
Query Earth Search for Landsat Collection 2 Level-2 scenes around Lake Tahoe. The bbox (-121.5, 38.5, -118.5, 40.5) spans the UTM 10N / 11N boundary at -120° longitude, so we get scenes in both zones.
We request 6 scenes and ingest the red, green, and blue asset hrefs.
In [8]:
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,
collections=["landsat-c2-l2"],
bbox=(-121.5, 38.5, -118.5, 40.5),
max_items=6,
datetime="2025-10-01/2025-10-15",
asset_hrefs=["red", "green", "blue"],
)
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,
collections=["landsat-c2-l2"],
bbox=(-121.5, 38.5, -118.5, 40.5),
max_items=6,
datetime="2025-10-01/2025-10-15",
asset_hrefs=["red", "green", "blue"],
)
rows
Out[8]:
6
4. Query metadata with DataFusion¶
Use DataFusion SQL to inspect the ingested metadata and verify we have scenes in both UTM zones.
In [9]:
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[9]:
| asset_blue | asset_green | asset_red | bbox | created | datetime | description | earthsearch:payload_id | eo:cloud_cover | grid:code | gsd | id | landsat:cloud_cover_land | landsat:collection_category | landsat:collection_number | landsat:correction | landsat:product_generated | landsat:scene_id | landsat:wrs_path | landsat:wrs_row | landsat:wrs_type | platform | proj:epsg | sci:doi | shape_x | shape_y | transform_0 | transform_1 | transform_2 | transform_3 | transform_4 | transform_5 | updated | view:off_nadir | view:sun_azimuth | view:sun_elevation |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 30.47 | WRS2-043032 | 30 | 30.47 | T1 | 02 | L2SP | 2025-10-12T16:19:40Z | LC90430322025284LGN00 | 043 | 032 | 2 | landsat-9 | 32611 | 10.5066/P9OGBGM6 | 7901 | 8011 | 30.0 | 0.0 | 150885.0 | 0.0 | -30.0 | 4587915.0 | 0 | 158.43982425 | 39.96946877 | ||||||||||
| 9.05 | WRS2-043033 | 30 | 9.05 | T1 | 02 | L2SP | 2025-10-12T16:32:40Z | LC90430332025284LGN00 | 043 | 033 | 2 | landsat-9 | 32610 | 10.5066/P9OGBGM6 | 7581 | 7701 | 30.0 | 0.0 | 632385.0 | 0.0 | -30.0 | 4425015.0 | 0 | 157.54660497 | 41.19354606 | ||||||||||
| 3.34 | WRS2-043034 | 30 | 3.35 | T1 | 02 | L2SP | 2025-10-12T16:22:19Z | LC90430342025284LGN00 | 043 | 034 | 2 | landsat-9 | 32610 | 10.5066/P9OGBGM6 | 7591 | 7721 | 30.0 | 0.0 | 598785.0 | 0.0 | -30.0 | 4265415.0 | 0 | 156.62566494 | 42.408338 | ||||||||||
| 0.42 | WRS2-042032 | 30 | 0.42 | T1 | 02 | L2SP | 2025-11-18T02:08:11Z | LC80420322025285LGN00 | 042 | 032 | 2 | landsat-8 | 32611 | 10.5066/P9OGBGM6 | 7811 | 7931 | 30.0 | 0.0 | 284685.0 | 0.0 | -30.0 | 4583715.0 | 0 | 158.64254999 | 39.62665708 | ||||||||||
| 0.12 | WRS2-042033 | 30 | 0.12 | T1 | 02 | L2SP | 2025-11-18T12:16:19Z | LC80420332025285LGN00 | 042 | 033 | 2 | landsat-8 | 32611 | 10.5066/P9OGBGM6 | 7811 | 7941 | 30.0 | 0.0 | 243885.0 | 0.0 | -30.0 | 4425915.0 | 0 | 157.76121163 | 40.85308421 | ||||||||||
| 0.57 | WRS2-042034 | 30 | 0.57 | T1 | 02 | L2SP | 2025-11-18T02:22:30Z | LC80420342025285LGN00 | 042 | 034 | 2 | landsat-8 | 32611 | 10.5066/P9OGBGM6 | 7821 | 7941 | 30.0 | 0.0 | 202485.0 | 0.0 | -30.0 | 4268115.0 | 0 | 156.85262094 | 42.07091531 |
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 [10]:
Copied!
virtualizer = default_virtualizer(registry, target_resolution=30.0)
virtualizer = default_virtualizer(registry, target_resolution=30.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 [11]:
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=30.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
nodata=0,
virtualize=virtualizer,
)
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=30.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
nodata=0,
virtualize=virtualizer,
)
6b. Explain plan (dry-run)¶
Before running compute, use explain() to see exactly what source regions
each target chunk would read — without fetching any pixel data. This helps
diagnose performance and verify overview selection.
In [11]:
Copied!
from lazymerge.explain import explain
plan = explain(
store=ic_session.store,
crs="EPSG:32611",
bbox=(200000.0, 4260000.0, 340000.0, 4490000.0),
resolution=30.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
nodata=0,
)
from lazymerge.explain import explain
plan = explain(
store=ic_session.store,
crs="EPSG:32611",
bbox=(200000.0, 4260000.0, 340000.0, 4490000.0),
resolution=30.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
nodata=0,
)
In [12]:
Copied!
plan.to_dataframe()
plan.to_dataframe()
Out[12]:
| chunk_row | chunk_col | chunk_height | chunk_width | n_source_reads | source_path | source_crs | overview_path | native_resolution | overview_resolution | region_row_start | region_row_end | region_col_start | region_col_end | region_height | region_width | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 1024 | 1024 | 3 | LC09_L2SP_043032_20251011_02_T1/red/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 3263 | 4288 | 1637 | 2662 | 1025 | 1025 |
| 1 | 0 | 0 | 1024 | 1024 | 3 | LC09_L2SP_043032_20251011_02_T1/green/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 3263 | 4288 | 1637 | 2662 | 1025 | 1025 |
| 2 | 0 | 0 | 1024 | 1024 | 3 | LC09_L2SP_043032_20251011_02_T1/blue/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 3263 | 4288 | 1637 | 2662 | 1025 | 1025 |
| 3 | 0 | 1 | 1024 | 1024 | 3 | LC09_L2SP_043032_20251011_02_T1/red/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 3263 | 4288 | 2661 | 3686 | 1025 | 1025 |
| 4 | 0 | 1 | 1024 | 1024 | 3 | LC09_L2SP_043032_20251011_02_T1/green/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 3263 | 4288 | 2661 | 3686 | 1025 | 1025 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 304 | 7 | 4 | 499 | 571 | 9 | LC08_L2SP_042033_20251012_02_T1/green/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 5031 | 5531 | 2633 | 3205 | 500 | 572 |
| 305 | 7 | 4 | 499 | 571 | 9 | LC08_L2SP_042033_20251012_02_T1/blue/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 5031 | 5531 | 2633 | 3205 | 500 | 572 |
| 306 | 7 | 4 | 499 | 571 | 9 | LC08_L2SP_042034_20251012_02_T1/red/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 0 | 271 | 4013 | 4585 | 271 | 572 |
| 307 | 7 | 4 | 499 | 571 | 9 | LC08_L2SP_042034_20251012_02_T1/green/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 0 | 271 | 4013 | 4585 | 271 | 572 |
| 308 | 7 | 4 | 499 | 571 | 9 | LC08_L2SP_042034_20251012_02_T1/blue/0 | EPSG:32611 | 0 | 30.0 | 30.0 | 0 | 271 | 4013 | 4585 | 271 | 572 |
309 rows × 16 columns
7. Open the lazy Cubed array with xarray¶
In [12]:
Copied!
da = xr.DataArray(result_arr, dims=["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({"bands": bands})
da
da = xr.DataArray(result_arr, dims=["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({"bands": bands})
da
Out[12]:
<xarray.DataArray 'array-004' (bands: 3, y: 7667, x: 4667)> Size: 429MB
cubed.Array<array-004, shape=(3, 7667, 4667), dtype=float32, chunks=((1, 1, 1), (1024, 1024, 1024, 1024, 1024, 1024, 1024, 499), (1024, 1024, 1024, 1024, 571))>
Coordinates:
* bands (bands) <U5 60B 'red' 'green' 'blue'
* y (y) float64 61kB 4.49e+06 4.49e+06 ... 4.26e+06 4.26e+06
* x (x) float64 37kB 2e+05 2e+05 2.001e+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 — each chunk task virtualizes any sources it needs before reading pixel data.
In [13]:
Copied!
import logging
logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
import logging
logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
In [14]:
Copied!
da = da.compute()
da
da = da.compute()
da
INFO:lazymerge.virtualize:Virtualized LC09_L2SP_043032_20251011_02_T1 (3 bands, levels [0], native_res=30.0) INFO:lazymerge.virtualize:Virtualized LC08_L2SP_042032_20251012_02_T1 (3 bands, levels [0], native_res=30.0) INFO:lazymerge.virtualize:Virtualized LC08_L2SP_042033_20251012_02_T1 (3 bands, levels [0], native_res=30.0) INFO:lazymerge.virtualize:Virtualized LC09_L2SP_043033_20251011_02_T1 (3 bands, levels [0], native_res=30.0) INFO:lazymerge.virtualize:Virtualized LC09_L2SP_043034_20251011_02_T1 (3 bands, levels [0], native_res=30.0) INFO:lazymerge.virtualize:Virtualized LC08_L2SP_042034_20251012_02_T1 (3 bands, levels [0], native_res=30.0)
Out[14]:
<xarray.DataArray 'array-004' (bands: 3, y: 7667, x: 4667)> Size: 429MB
array([[[29472., 28274., 26576., ..., 12613., 12388., 12296.],
[24593., 23967., 20958., ..., 12663., 13021., 12416.],
[16734., 17361., 15448., ..., 12481., 13763., 13227.],
...,
[ 7672., 8198., 8366., ..., 11542., 11750., 11939.],
[ 7704., 8113., 8495., ..., 11527., 11862., 12243.],
[ 7887., 8010., 8582., ..., 11931., 11721., 12248.]],
[[28647., 27264., 25150., ..., 11548., 11376., 11334.],
[22637., 22411., 19189., ..., 11544., 11798., 11368.],
[14974., 15456., 13774., ..., 11404., 12367., 11998.],
...,
[ 7831., 8421., 8566., ..., 10844., 11005., 11142.],
[ 7855., 8397., 8793., ..., 10819., 11078., 11359.],
[ 7885., 8181., 8748., ..., 11144., 10981., 11372.]],
[[29795., 28273., 26883., ..., 10193., 10105., 10051.],
[26258., 23727., 20988., ..., 10325., 10346., 10143.],
[17353., 17270., 14127., ..., 10134., 10825., 10439.],
...,
[ 7531., 7805., 7914., ..., 9563., 9661., 9907.],
[ 7359., 7617., 7910., ..., 9689., 9828., 10043.],
[ 7508., 7582., 7984., ..., 9860., 9812., 10051.]]],
shape=(3, 7667, 4667), dtype=float32)
Coordinates:
* bands (bands) <U5 60B 'red' 'green' 'blue'
* y (y) float64 61kB 4.49e+06 4.49e+06 ... 4.26e+06 4.26e+06
* x (x) float64 37kB 2e+05 2e+05 2.001e+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)We can see that the Virtualizarr arrays have been created for all the sources required by the chunks in the Cubed array.
In [15]:
Copied!
root = zarr.open(ic_session.store)
root.tree()
root = zarr.open(ic_session.store)
root.tree()
Out[15]:
/
├── LC08_L2SP_042032_20251012_02_T1
│ ├── blue
│ │ └── 0
│ │ └── 0 (7931, 7811) uint16
│ ├── green
│ │ └── 0
│ │ └── 0 (7931, 7811) uint16
│ └── red
│ └── 0
│ └── 0 (7931, 7811) uint16
├── LC08_L2SP_042033_20251012_02_T1
│ ├── blue
│ │ └── 0
│ │ └── 0 (7941, 7811) uint16
│ ├── green
│ │ └── 0
│ │ └── 0 (7941, 7811) uint16
│ └── red
│ └── 0
│ └── 0 (7941, 7811) uint16
├── LC08_L2SP_042034_20251012_02_T1
│ ├── blue
│ │ └── 0
│ │ └── 0 (7941, 7821) uint16
│ ├── green
│ │ └── 0
│ │ └── 0 (7941, 7821) uint16
│ └── red
│ └── 0
│ └── 0 (7941, 7821) uint16
├── LC09_L2SP_043032_20251011_02_T1
│ ├── blue
│ │ └── 0
│ │ └── 0 (8011, 7901) uint16
│ ├── green
│ │ └── 0
│ │ └── 0 (8011, 7901) uint16
│ └── red
│ └── 0
│ └── 0 (8011, 7901) uint16
├── LC09_L2SP_043033_20251011_02_T1
│ ├── blue
│ │ └── 0
│ │ └── 0 (7701, 7581) uint16
│ ├── green
│ │ └── 0
│ │ └── 0 (7701, 7581) uint16
│ └── red
│ └── 0
│ └── 0 (7701, 7581) uint16
├── LC09_L2SP_043034_20251011_02_T1
│ ├── blue
│ │ └── 0
│ │ └── 0 (7721, 7591) uint16
│ ├── green
│ │ └── 0
│ │ └── 0 (7721, 7591) uint16
│ └── red
│ └── 0
│ └── 0 (7721, 7591) uint16
└── meta
├── asset_blue (6,) StringDType()
├── asset_green (6,) StringDType()
├── asset_red (6,) StringDType()
├── bbox (6,) object
├── created (6,) datetime64[ms]
├── datetime (6,) datetime64[ms]
├── description (6,) StringDType()
├── earthsearch:payload_id (6,) StringDType()
├── eo:cloud_cover (6,) float64
├── grid:code (6,) StringDType()
├── gsd (6,) int64
├── id (6,) StringDType()
├── landsat:cloud_cover_land (6,) float64
├── landsat:collection_category (6,) StringDType()
├── landsat:collection_number (6,) StringDType()
├── landsat:correction (6,) StringDType()
├── landsat:product_generated (6,) StringDType()
├── landsat:scene_id (6,) StringDType()
├── landsat:wrs_path (6,) StringDType()
├── landsat:wrs_row (6,) StringDType()
├── landsat:wrs_type (6,) StringDType()
├── platform (6,) StringDType()
├── proj:epsg (6,) int64
├── sci:doi (6,) StringDType()
├── shape_x (6,) int64
├── shape_y (6,) int64
├── transform_0 (6,) float64
├── transform_1 (6,) float64
├── transform_2 (6,) float64
├── transform_3 (6,) float64
├── transform_4 (6,) float64
├── transform_5 (6,) float64
├── updated (6,) datetime64[ms]
├── view:off_nadir (6,) int64
├── view:sun_azimuth (6,) float64
└── view:sun_elevation (6,) float64
In [16]:
Copied!
width = 8
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.title("Landsat rgb — merged from UTM 10N + 11N sources")
plt.tight_layout()
width = 8
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.title("Landsat rgb — merged from UTM 10N + 11N sources")
plt.tight_layout()
Slicing support¶
In [17]:
Copied!
mid_y = (da.y.values[0] + da.y.values[-1]) / 2
da = da.sel(y=slice(mid_y, da.y.values[-1]))
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.title("Landsat rgb — merged from UTM 10N + 11N sources")
plt.tight_layout()
mid_y = (da.y.values[0] + da.y.values[-1]) / 2
da = da.sel(y=slice(mid_y, da.y.values[-1]))
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.title("Landsat rgb — merged from UTM 10N + 11N sources")
plt.tight_layout()
9. Use additional SQL filters¶
We can also add additional SQL filters to narrow the criteria for source granules
In [18]:
Copied!
result_arr, result_spatial, result_proj, time_coord = merge(
store=ic_session.store,
crs="EPSG:32611",
bbox=(200000.0, 4260000.0, 340000.0, 4490000.0),
resolution=30.0,
chunk_size=(1024, 1024),
bands=["red", "green", "blue"],
datafusion=True,
sortby="datetime",
nodata=0,
sql_filter='"eo:cloud_cover" < 30',
virtualize=virtualizer,
)
result_arr, result_spatial, result_proj, time_coord = merge(
store=ic_session.store,
crs="EPSG:32611",
bbox=(200000.0, 4260000.0, 340000.0, 4490000.0),
resolution=30.0,
chunk_size=(1024, 1024),
bands=["red", "green", "blue"],
datafusion=True,
sortby="datetime",
nodata=0,
sql_filter='"eo:cloud_cover" < 30',
virtualize=virtualizer,
)
In [19]:
Copied!
da = xr.DataArray(result_arr, dims=["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)
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.title("Landsat rgb — cloud_cover < 30")
plt.tight_layout()
da = xr.DataArray(result_arr, dims=["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)
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.title("Landsat rgb — cloud_cover < 30")
plt.tight_layout()
10. Temporal grouping¶
Merge also supports a temporal_grouping parameter that allows aggregation into temporal bins and returns a coordinate array that can be used when constructing an Xarray DataArray.
In [20]:
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=30.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
nodata=0,
temporal_grouping="P1D",
virtualize=virtualizer,
)
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=30.0,
chunk_size=(1024, 1024),
bands=bands,
datafusion=True,
sortby="datetime",
nodata=0,
temporal_grouping="P1D",
virtualize=virtualizer,
)
In [21]:
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[21]:
<xarray.DataArray 'array-012' (time: 2, bands: 3, y: 7667, x: 4667)> Size: 859MB
cubed.Array<array-012, shape=(2, 3, 7667, 4667), dtype=float32, chunks=((1, 1), (1, 1, 1), (1024, 1024, 1024, 1024, 1024, 1024, 1024, 499), (1024, 1024, 1024, 1024, 571))>
Coordinates:
* time (time) datetime64[s] 16B 2025-10-11 2025-10-12
* bands (bands) <U5 60B 'red' 'green' 'blue'
* y (y) float64 61kB 4.49e+06 4.49e+06 ... 4.26e+06 4.26e+06
* x (x) float64 37kB 2e+05 2e+05 2.001e+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 [22]:
Copied!
da = da.sel(time="2025-10-12")
da
da = da.sel(time="2025-10-12")
da
Out[22]:
<xarray.DataArray 'array-012' (bands: 3, y: 7667, x: 4667)> Size: 429MB
cubed.Array<array-014, shape=(3, 7667, 4667), dtype=float32, chunks=((1, 1, 1), (1024, 1024, 1024, 1024, 1024, 1024, 1024, 499), (1024, 1024, 1024, 1024, 571))>
Coordinates:
* bands (bands) <U5 60B 'red' 'green' 'blue'
* y (y) float64 61kB 4.49e+06 4.49e+06 ... 4.26e+06 4.26e+06
* x (x) float64 37kB 2e+05 2e+05 2.001e+05 ... 3.4e+05 3.4e+05
* spatial_ref int64 8B 0
time datetime64[s] 8B 2025-10-12
Indexes:
┌ x RasterIndex (crs=EPSG:32611)
└ y
spatial_ref CRSIndex (crs=EPSG:32611)In [23]:
Copied!
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.tight_layout()
height = width * da.sizes["y"] / da.sizes["x"]
da.plot.imshow(rgb="bands", robust=True, figsize=(width, height))
plt.tight_layout()
