lazymerge datafusion backend¶

Same mosaic as the conventions demo, but using the datafusion=True merge path. Source arrays are organized in a zarr-datafusion-search store layout: a /meta group with columnar metadata arrays, and source groups containing band subgroups with multiscale overviews.

Instead of building an in-memory ScanIndex upfront, merge(datafusion=True) issues per-chunk spatial SQL queries against the /meta group to find intersecting sources. This scales to stores with many sources.

The store parameter in merge() accepts any zarr-compatible store (zarr Store, Icechunk Session) or an obstore object store (LocalStore, S3Store). When datafusion=True, the store must be obstore-compatible since zarr-datafusion-search queries via obstore.

In [1]:

import tempfile

import numpy as np
import xarray as xr
import zarr
from affine import Affine
from obstore.store import LocalStore
from rasterix import RasterIndex

from lazymerge.conventions import (
    ProjAttrs,
    SpatialAttrs,
    read_multiscales,
    read_spatial,
    write_proj,
    write_spatial,
)
from lazymerge.merge import merge
from lazymerge.sources import select_overview

import tempfile import numpy as np import xarray as xr import zarr from affine import Affine from obstore.store import LocalStore from rasterix import RasterIndex from lazymerge.conventions import ( ProjAttrs, SpatialAttrs, read_multiscales, read_spatial, write_proj, write_spatial, ) from lazymerge.merge import merge from lazymerge.sources import select_overview

1. Create synthetic source arrays with DataFusion store layout¶

Four tiles arranged in a 2x2 grid. The left two tiles are in UTM 18N (EPSG:32618), and the right two are in UTM 17N (EPSG:32617) — covering adjacent geographic areas but stored in different coordinate reference systems.

Each source is a group containing a data/ band subgroup with base (10m) and overview (20m) arrays, plus the zarr multiscales convention.

We write to a local directory store so we can pass an obstore LocalStore to merge(). DataFusion queries require an obstore-compatible store.

store_path/
  meta/               # columnar metadata for DataFusion queries
  utm18n_tile_a/
    data/
      multiscales      # convention attr
      0/               # base (10m)
      1/               # overview (20m)
  ...

In [2]:

tmpdir = tempfile.mkdtemp()
store = zarr.storage.LocalStore(tmpdir)
root = zarr.open_group(store, mode="w")

sources = [
    # UTM 18N tiles (left side)
    {
        "name": "utm18n_tile_a",
        "crs": "EPSG:32618",
        "epsg": 32618,
        "bbox": (500000.0, 5990000.0, 505000.0, 5995000.0),
        "resolution": 10.0,
        "fill": 100.0,
    },
    {
        "name": "utm18n_tile_c",
        "crs": "EPSG:32618",
        "epsg": 32618,
        "bbox": (500000.0, 5995000.0, 505000.0, 6000000.0),
        "resolution": 10.0,
        "fill": 300.0,
    },
    # UTM 17N tiles (right side — same geographic area as UTM 18N
    # x=[503000, 508000], but expressed in UTM zone 17N coordinates)
    {
        "name": "utm17n_tile_b",
        "crs": "EPSG:32617",
        "epsg": 32617,
        "bbox": (895094.0, 6006920.0, 900511.0, 6012337.0),
        "resolution": 10.0,
        "fill": 200.0,
    },
    {
        "name": "utm17n_tile_d",
        "crs": "EPSG:32617",
        "epsg": 32617,
        "bbox": (894669.0, 6011912.0, 900086.0, 6017329.0),
        "resolution": 10.0,
        "fill": 400.0,
    },
]

for src in sources:
    xmin, ymin, xmax, ymax = src["bbox"]
    res = src["resolution"]
    width = int((xmax - xmin) / res)
    height = int((ymax - ymin) / res)

    scene_group = root.create_group(src["name"])
    band_group = scene_group.create_group("data")

    # Base resolution array (level 0)
    base = band_group.create_array("0", shape=(height, width), dtype="f4", chunks=(256, 256))
    base[:] = src["fill"]
    write_spatial(
        base,
        SpatialAttrs(
            dimensions=["y", "x"],
            transform=(res, 0.0, xmin, 0.0, -res, ymax),
            bbox=(xmin, ymin, xmax, ymax),
            shape=(height, width),
        ),
    )
    write_proj(base, ProjAttrs(code=src["crs"]))

    # 2x overview array (level 1)
    ovr_width = width // 2
    ovr_height = height // 2
    ovr_res = res * 2
    ovr = band_group.create_array("1", shape=(ovr_height, ovr_width), dtype="f4", chunks=(256, 256))
    ovr[:] = src["fill"]
    write_spatial(
        ovr,
        SpatialAttrs(
            dimensions=["y", "x"],
            transform=(ovr_res, 0.0, xmin, 0.0, -ovr_res, ymax),
            bbox=(xmin, ymin, xmax, ymax),
            shape=(ovr_height, ovr_width),
        ),
    )
    write_proj(ovr, ProjAttrs(code=src["crs"]))

    # Set multiscales convention on band group
    band_group.attrs["multiscales"] = {
        "layout": [
            {"asset": "0", "transform": {"scale": [1.0, 1.0], "translation": [0.0, 0.0]}},
            {
                "asset": "1",
                "derived_from": "0",
                "transform": {"scale": [2.0, 2.0], "translation": [0.5, 0.5]},
            },
        ],
    }


# Create obstore LocalStore pointing at the same directory for DataFusion queries
ob_store = LocalStore(tmpdir)

tmpdir = tempfile.mkdtemp() store = zarr.storage.LocalStore(tmpdir) root = zarr.open_group(store, mode="w") sources = [ # UTM 18N tiles (left side) { "name": "utm18n_tile_a", "crs": "EPSG:32618", "epsg": 32618, "bbox": (500000.0, 5990000.0, 505000.0, 5995000.0), "resolution": 10.0, "fill": 100.0, }, { "name": "utm18n_tile_c", "crs": "EPSG:32618", "epsg": 32618, "bbox": (500000.0, 5995000.0, 505000.0, 6000000.0), "resolution": 10.0, "fill": 300.0, }, # UTM 17N tiles (right side — same geographic area as UTM 18N # x=[503000, 508000], but expressed in UTM zone 17N coordinates) { "name": "utm17n_tile_b", "crs": "EPSG:32617", "epsg": 32617, "bbox": (895094.0, 6006920.0, 900511.0, 6012337.0), "resolution": 10.0, "fill": 200.0, }, { "name": "utm17n_tile_d", "crs": "EPSG:32617", "epsg": 32617, "bbox": (894669.0, 6011912.0, 900086.0, 6017329.0), "resolution": 10.0, "fill": 400.0, }, ] for src in sources: xmin, ymin, xmax, ymax = src["bbox"] res = src["resolution"] width = int((xmax - xmin) / res) height = int((ymax - ymin) / res) scene_group = root.create_group(src["name"]) band_group = scene_group.create_group("data") # Base resolution array (level 0) base = band_group.create_array("0", shape=(height, width), dtype="f4", chunks=(256, 256)) base[:] = src["fill"] write_spatial( base, SpatialAttrs( dimensions=["y", "x"], transform=(res, 0.0, xmin, 0.0, -res, ymax), bbox=(xmin, ymin, xmax, ymax), shape=(height, width), ), ) write_proj(base, ProjAttrs(code=src["crs"])) # 2x overview array (level 1) ovr_width = width // 2 ovr_height = height // 2 ovr_res = res * 2 ovr = band_group.create_array("1", shape=(ovr_height, ovr_width), dtype="f4", chunks=(256, 256)) ovr[:] = src["fill"] write_spatial( ovr, SpatialAttrs( dimensions=["y", "x"], transform=(ovr_res, 0.0, xmin, 0.0, -ovr_res, ymax), bbox=(xmin, ymin, xmax, ymax), shape=(ovr_height, ovr_width), ), ) write_proj(ovr, ProjAttrs(code=src["crs"])) # Set multiscales convention on band group band_group.attrs["multiscales"] = { "layout": [ {"asset": "0", "transform": {"scale": [1.0, 1.0], "translation": [0.0, 0.0]}}, { "asset": "1", "derived_from": "0", "transform": {"scale": [2.0, 2.0], "translation": [0.5, 0.5]}, }, ], } # Create obstore LocalStore pointing at the same directory for DataFusion queries ob_store = LocalStore(tmpdir)

2. Populate the /meta group¶

The /meta group holds columnar arrays that DataFusion queries against. Each row corresponds to one source scene. By convention, bbox values are stored in EPSG:4326 so queries work in a single CRS.

In [3]:

import struct

from pyproj import Transformer


def bbox_to_wkb(xmin, ymin, xmax, ymax):
    """Encode a bounding box as a WKB Polygon (little-endian)."""
    points = [(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax), (xmin, ymin)]
    wkb = struct.pack("<BII", 1, 3, 1)  # LE, Polygon, 1 ring
    wkb += struct.pack("<I", 5)  # 5 points
    for x, y in points:
        wkb += struct.pack("<dd", x, y)
    return wkb


meta = root.create_group("meta")

ids = np.array([s["name"] for s in sources])
epsgs = np.array([s["epsg"] for s in sources], dtype="int64")

# Affine transform components: (a=res, b=0, c=xmin, d=0, e=-res, f=ymax)
t0 = np.array([s["resolution"] for s in sources])
t1 = np.zeros(len(sources))
t2 = np.array([s["bbox"][0] for s in sources])  # xmin
t3 = np.zeros(len(sources))
t4 = np.array([-s["resolution"] for s in sources])
t5 = np.array([s["bbox"][3] for s in sources])  # ymax

shape_x = np.array(
    [int((s["bbox"][2] - s["bbox"][0]) / s["resolution"]) for s in sources], dtype="int64"
)
shape_y = np.array(
    [int((s["bbox"][3] - s["bbox"][1]) / s["resolution"]) for s in sources], dtype="int64"
)

# Compute bboxes in EPSG:4326 and encode as WKB (DataFusion convention)
bbox_wkb = []
for s in sources:
    xmin, ymin, xmax, ymax = s["bbox"]
    transformer = Transformer.from_crs(s["crs"], "EPSG:4326", always_xy=True)
    xs = [xmin, xmin, xmax, xmax]
    ys = [ymin, ymax, ymin, ymax]
    tx, ty = transformer.transform(xs, ys)
    lon_min, lon_max = min(tx), max(tx)
    lat_min, lat_max = min(ty), max(ty)
    bbox_wkb.append(bbox_to_wkb(lon_min, lat_min, lon_max, lat_max))

bbox_wkb_arr = np.array(bbox_wkb, dtype=object)

# Use dtype='string' for variable-length UTF-8 (zarr VariableLengthUTF8)
id_arr = meta.create_array("id", shape=ids.shape, dtype="string")
id_arr[:] = ids
# Store bbox as variable-length bytes (WKB polygons in EPSG:4326)
bbox_arr = meta.create_array("bbox", shape=bbox_wkb_arr.shape, dtype="variable_length_bytes")
bbox_arr[:] = bbox_wkb_arr
meta.create_array("proj:epsg", data=epsgs)
meta.create_array("transform_0", data=t0)
meta.create_array("transform_1", data=t1)
meta.create_array("transform_2", data=t2)
meta.create_array("transform_3", data=t3)
meta.create_array("transform_4", data=t4)
meta.create_array("transform_5", data=t5)
meta.create_array("shape_x", data=shape_x)
meta.create_array("shape_y", data=shape_y)

import struct from pyproj import Transformer def bbox_to_wkb(xmin, ymin, xmax, ymax): """Encode a bounding box as a WKB Polygon (little-endian).""" points = [(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax), (xmin, ymin)] wkb = struct.pack("<BII", 1, 3, 1) # LE, Polygon, 1 ring wkb += struct.pack("<I", 5) # 5 points for x, y in points: wkb += struct.pack("<dd", x, y) return wkb meta = root.create_group("meta") ids = np.array([s["name"] for s in sources]) epsgs = np.array([s["epsg"] for s in sources], dtype="int64") # Affine transform components: (a=res, b=0, c=xmin, d=0, e=-res, f=ymax) t0 = np.array([s["resolution"] for s in sources]) t1 = np.zeros(len(sources)) t2 = np.array([s["bbox"][0] for s in sources]) # xmin t3 = np.zeros(len(sources)) t4 = np.array([-s["resolution"] for s in sources]) t5 = np.array([s["bbox"][3] for s in sources]) # ymax shape_x = np.array( [int((s["bbox"][2] - s["bbox"][0]) / s["resolution"]) for s in sources], dtype="int64" ) shape_y = np.array( [int((s["bbox"][3] - s["bbox"][1]) / s["resolution"]) for s in sources], dtype="int64" ) # Compute bboxes in EPSG:4326 and encode as WKB (DataFusion convention) bbox_wkb = [] for s in sources: xmin, ymin, xmax, ymax = s["bbox"] transformer = Transformer.from_crs(s["crs"], "EPSG:4326", always_xy=True) xs = [xmin, xmin, xmax, xmax] ys = [ymin, ymax, ymin, ymax] tx, ty = transformer.transform(xs, ys) lon_min, lon_max = min(tx), max(tx) lat_min, lat_max = min(ty), max(ty) bbox_wkb.append(bbox_to_wkb(lon_min, lat_min, lon_max, lat_max)) bbox_wkb_arr = np.array(bbox_wkb, dtype=object) # Use dtype='string' for variable-length UTF-8 (zarr VariableLengthUTF8) id_arr = meta.create_array("id", shape=ids.shape, dtype="string") id_arr[:] = ids # Store bbox as variable-length bytes (WKB polygons in EPSG:4326) bbox_arr = meta.create_array("bbox", shape=bbox_wkb_arr.shape, dtype="variable_length_bytes") bbox_arr[:] = bbox_wkb_arr meta.create_array("proj:epsg", data=epsgs) meta.create_array("transform_0", data=t0) meta.create_array("transform_1", data=t1) meta.create_array("transform_2", data=t2) meta.create_array("transform_3", data=t3) meta.create_array("transform_4", data=t4) meta.create_array("transform_5", data=t5) meta.create_array("shape_x", data=shape_x) meta.create_array("shape_y", data=shape_y)

/Users/seanharkins/projects/lazymerge/docs/examples/.venv/lib/python3.13/site-packages/zarr/core/dtype/npy/bytes.py:1143: UnstableSpecificationWarning: The data type (VariableLengthBytes()) does not have a Zarr V3 specification. That means that the representation of arrays saved with this data type may change without warning in a future version of Zarr Python. Arrays stored with this data type may be unreadable by other Zarr libraries. Use this data type at your own risk! Check https://github.com/zarr-developers/zarr-extensions/tree/main/data-types for the status of data type specifications for Zarr V3.
  v3_unstable_dtype_warning(self)

Out[3]:

<Array file:///var/folders/mp/33cxt8xj36bbxj0jqdwbfwyc0000gn/T/tmp900t4_0d/meta/shape_y shape=(4,) dtype=int64>

3. Overview selection¶

Before running the full merge, let's demonstrate how overview selection works. This part uses only lazymerge (no DataFusion dependency).

In [4]:

band_group = root["utm18n_tile_a/data"]
overviews = read_multiscales(band_group)
base_spatial = read_spatial(band_group["0"])
native_res = abs(base_spatial.transform[0])

for _ov in overviews:
    pass

# At 10m target, no overview selected (use base)
sel_10 = select_overview(overviews, target_res=10.0, native_res=native_res)

# At 20m target, overview 1 is selected
sel_20 = select_overview(overviews, target_res=20.0, native_res=native_res)

band_group = root["utm18n_tile_a/data"] overviews = read_multiscales(band_group) base_spatial = read_spatial(band_group["0"]) native_res = abs(base_spatial.transform[0]) for _ov in overviews: pass # At 10m target, no overview selected (use base) sel_10 = select_overview(overviews, target_res=10.0, native_res=native_res) # At 20m target, overview 1 is selected sel_20 = select_overview(overviews, target_res=20.0, native_res=native_res)

4. Merge with DataFusion¶

With datafusion=True, merge() does not need a pre-built ScanIndex. Instead, each target chunk issues a DataFusion SQL query against the /meta group to find intersecting sources by bbox.

The store must be obstore-compatible (e.g. LocalStore, S3Store) when using the DataFusion path. For the convention path, any zarr-compatible store works (zarr Store, Icechunk Session, or a string path).

The sortby parameter controls the compositing order — sources returned earlier take priority when filling NaN pixels.

# No scan_store() call needed!
result_arr, result_spatial, result_proj = merge(
    store=ob_store,          # obstore LocalStore for DataFusion queries
    crs="EPSG:32618",
    bbox=(500000.0, 5990000.0, 508000.0, 6000000.0),
    resolution=10.0,
    chunk_size=(256, 256),
    bands="data",            # navigate into data/ band group
    datafusion=True,         # use DataFusion for source discovery
    sortby="id",             # order sources by id
)

3. Lazy merge¶

In [5]:

# Requires zarr-datafusion-search and obstore
result_arr, result_spatial, result_proj, time_coords = merge(
    store=ob_store,
    crs="EPSG:32618",
    bbox=(500000.0, 5990000.0, 508000.0, 6000000.0),
    resolution=10.0,
    chunk_size=(256, 256),
    bands="data",
    datafusion=True,
    sortby="id",
)

# Requires zarr-datafusion-search and obstore result_arr, result_spatial, result_proj, time_coords = merge( store=ob_store, crs="EPSG:32618", bbox=(500000.0, 5990000.0, 508000.0, 6000000.0), resolution=10.0, chunk_size=(256, 256), bands="data", datafusion=True, sortby="id", )

4. Open the lazy Cubed array with xarray¶

In [6]:

da = xr.DataArray(result_arr, dims=["y", "x"])

# Build georeferenced coordinates using rasterix
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 = xr.DataArray(result_arr, dims=["y", "x"]) # Build georeferenced coordinates using rasterix 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

Out[6]:

<xarray.DataArray 'array-003' (y: 1000, x: 800)> Size: 3MB
cubed.Array<array-003, shape=(1000, 800), dtype=float32, chunks=((256, 256, 256, 232), (256, 256, 256, 32))>
Coordinates:
  * y            (y) float64 8kB 6e+06 6e+06 6e+06 ... 5.99e+06 5.99e+06
  * x            (x) float64 6kB 5e+05 5e+05 5e+05 ... 5.08e+05 5.08e+05
  * spatial_ref  int64 8B 0
Indexes:
  ┌ x            RasterIndex (crs=EPSG:32618)
  └ y
    spatial_ref  CRSIndex (crs=EPSG:32618)

xarray.DataArray

'array-003'

• y: 1000
• x: 800

• cubed.Array<chunksize=(256, 256)>

  Array Chunk Bytes 3.2 MB 262.1 KB Shape (1000, 800) (256, 256) Count 3 arrays in Plan 16 Chunks Type float32 np.ndarray

• Coordinates: (3)
  • y
    (y)
    float64
    6e+06 6e+06 ... 5.99e+06 5.99e+06

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    [1000 values with dtype=float64]

  • x
    (x)
    float64
    5e+05 5e+05 ... 5.08e+05 5.08e+05

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

    [800 values with dtype=float64]

  • spatial_ref
    ()
    int64
    0

    array(0)

• Indexes: (2)
  • x
    y
    RasterIndex (crs=EPSG:32618)

    RasterIndex(crs=EPSG:32618)
        AxisAffineTransformIndex(AxisAffineTransform(a=10, b=0, c=5e+05, d=0, e=-10, f=6e+06, axis=X, dim='x'))
        AxisAffineTransformIndex(AxisAffineTransform(a=10, b=0, c=5e+05, d=0, e=-10, f=6e+06, axis=Y, dim='y'))

  • spatial_ref
    CRSIndex (crs=EPSG:32618)

    CRSIndex
    <Projected CRS: EPSG:32618>
    Name: WGS 84 / UTM zone 18N
    Axis Info [cartesian]:
    - E[east]: Easting (metre)
    - N[north]: Northing (metre)
    Area of Use:
    - name: Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.
    - bounds: (-78.0, 0.0, -72.0, 84.0)
    Coordinate Operation:
    - name: UTM zone 18N
    - method: Transverse Mercator
    Datum: World Geodetic System 1984 ensemble
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

5. Compute and visualize¶

Compute the result and plot it.

In [7]:

da = da.compute()
da

da = da.compute() da

/Users/seanharkins/projects/lazymerge/lazymerge/sources.py:183: RuntimeWarning: Successfully reconstructed a store defined in another Python module. Connection pooling will not be shared across store instances.
  zarr_table = await ZarrTable.from_obstore(store, "/meta")

Out[7]:

<xarray.DataArray 'array-003' (y: 1000, x: 800)> Size: 3MB
array([[300., 300., 300., ..., 400., 400., 400.],
       [300., 300., 300., ..., 400., 400., 400.],
       [300., 300., 300., ..., 400., 400., 400.],
       ...,
       [100., 100., 100., ..., 200., 200., 200.],
       [100., 100., 100., ..., 200., 200.,  nan],
       [100., 100., 100., ..., 200., 200.,  nan]],
      shape=(1000, 800), dtype=float32)
Coordinates:
  * y            (y) float64 8kB 6e+06 6e+06 6e+06 ... 5.99e+06 5.99e+06
  * x            (x) float64 6kB 5e+05 5e+05 5e+05 ... 5.08e+05 5.08e+05
  * spatial_ref  int64 8B 0
Indexes:
  ┌ x            RasterIndex (crs=EPSG:32618)
  └ y
    spatial_ref  CRSIndex (crs=EPSG:32618)

xarray.DataArray

'array-003'

• y: 1000
• x: 800

• 300.0 300.0 300.0 300.0 300.0 300.0 ... 200.0 200.0 200.0 200.0 nan

  array([[300., 300., 300., ..., 400., 400., 400.],
         [300., 300., 300., ..., 400., 400., 400.],
         [300., 300., 300., ..., 400., 400., 400.],
         ...,
         [100., 100., 100., ..., 200., 200., 200.],
         [100., 100., 100., ..., 200., 200.,  nan],
         [100., 100., 100., ..., 200., 200.,  nan]],
        shape=(1000, 800), dtype=float32)

• Coordinates: (3)
  • y
    (y)
    float64
    6e+06 6e+06 ... 5.99e+06 5.99e+06

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([5999995., 5999985., 5999975., 5999965., 5999955., 5999945.,
           5999935., 5999925., 5999915., 5999905., 5999895., 5999885.,
           5999875., 5999865., 5999855., 5999845., 5999835., 5999825.,
           5999815., 5999805., 5999795., 5999785., 5999775., 5999765.,
           5999755., 5999745., 5999735., 5999725., 5999715., 5999705.,
           5999695., 5999685., 5999675., 5999665., 5999655., 5999645.,
           5999635., 5999625., 5999615., 5999605., 5999595., 5999585.,
           5999575., 5999565., 5999555., 5999545., 5999535., 5999525.,
           5999515., 5999505., 5999495., 5999485., 5999475., 5999465.,
           5999455., 5999445., 5999435., 5999425., 5999415., 5999405.,
           5999395., 5999385., 5999375., 5999365., 5999355., 5999345.,
           5999335., 5999325., 5999315., 5999305., 5999295., 5999285.,
           5999275., 5999265., 5999255., 5999245., 5999235., 5999225.,
           5999215., 5999205., 5999195., 5999185., 5999175., 5999165.,
           5999155., 5999145., 5999135., 5999125., 5999115., 5999105.,
           5999095., 5999085., 5999075., 5999065., 5999055., 5999045.,
           5999035., 5999025., 5999015., 5999005., 5998995., 5998985.,
           5998975., 5998965., 5998955., 5998945., 5998935., 5998925.,
           5998915., 5998905., 5998895., 5998885., 5998875., 5998865.,
           5998855., 5998845., 5998835., 5998825., 5998815., 5998805.,
    ...
           5991175., 5991165., 5991155., 5991145., 5991135., 5991125.,
           5991115., 5991105., 5991095., 5991085., 5991075., 5991065.,
           5991055., 5991045., 5991035., 5991025., 5991015., 5991005.,
           5990995., 5990985., 5990975., 5990965., 5990955., 5990945.,
           5990935., 5990925., 5990915., 5990905., 5990895., 5990885.,
           5990875., 5990865., 5990855., 5990845., 5990835., 5990825.,
           5990815., 5990805., 5990795., 5990785., 5990775., 5990765.,
           5990755., 5990745., 5990735., 5990725., 5990715., 5990705.,
           5990695., 5990685., 5990675., 5990665., 5990655., 5990645.,
           5990635., 5990625., 5990615., 5990605., 5990595., 5990585.,
           5990575., 5990565., 5990555., 5990545., 5990535., 5990525.,
           5990515., 5990505., 5990495., 5990485., 5990475., 5990465.,
           5990455., 5990445., 5990435., 5990425., 5990415., 5990405.,
           5990395., 5990385., 5990375., 5990365., 5990355., 5990345.,
           5990335., 5990325., 5990315., 5990305., 5990295., 5990285.,
           5990275., 5990265., 5990255., 5990245., 5990235., 5990225.,
           5990215., 5990205., 5990195., 5990185., 5990175., 5990165.,
           5990155., 5990145., 5990135., 5990125., 5990115., 5990105.,
           5990095., 5990085., 5990075., 5990065., 5990055., 5990045.,
           5990035., 5990025., 5990015., 5990005.])

  • x
    (x)
    float64
    5e+05 5e+05 ... 5.08e+05 5.08e+05

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([500005., 500015., 500025., 500035., 500045., 500055., 500065.,
           500075., 500085., 500095., 500105., 500115., 500125., 500135.,
           500145., 500155., 500165., 500175., 500185., 500195., 500205.,
           500215., 500225., 500235., 500245., 500255., 500265., 500275.,
           500285., 500295., 500305., 500315., 500325., 500335., 500345.,
           500355., 500365., 500375., 500385., 500395., 500405., 500415.,
           500425., 500435., 500445., 500455., 500465., 500475., 500485.,
           500495., 500505., 500515., 500525., 500535., 500545., 500555.,
           500565., 500575., 500585., 500595., 500605., 500615., 500625.,
           500635., 500645., 500655., 500665., 500675., 500685., 500695.,
           500705., 500715., 500725., 500735., 500745., 500755., 500765.,
           500775., 500785., 500795., 500805., 500815., 500825., 500835.,
           500845., 500855., 500865., 500875., 500885., 500895., 500905.,
           500915., 500925., 500935., 500945., 500955., 500965., 500975.,
           500985., 500995., 501005., 501015., 501025., 501035., 501045.,
           501055., 501065., 501075., 501085., 501095., 501105., 501115.,
           501125., 501135., 501145., 501155., 501165., 501175., 501185.,
           501195., 501205., 501215., 501225., 501235., 501245., 501255.,
           501265., 501275., 501285., 501295., 501305., 501315., 501325.,
           501335., 501345., 501355., 501365., 501375., 501385., 501395.,
    ...
           506655., 506665., 506675., 506685., 506695., 506705., 506715.,
           506725., 506735., 506745., 506755., 506765., 506775., 506785.,
           506795., 506805., 506815., 506825., 506835., 506845., 506855.,
           506865., 506875., 506885., 506895., 506905., 506915., 506925.,
           506935., 506945., 506955., 506965., 506975., 506985., 506995.,
           507005., 507015., 507025., 507035., 507045., 507055., 507065.,
           507075., 507085., 507095., 507105., 507115., 507125., 507135.,
           507145., 507155., 507165., 507175., 507185., 507195., 507205.,
           507215., 507225., 507235., 507245., 507255., 507265., 507275.,
           507285., 507295., 507305., 507315., 507325., 507335., 507345.,
           507355., 507365., 507375., 507385., 507395., 507405., 507415.,
           507425., 507435., 507445., 507455., 507465., 507475., 507485.,
           507495., 507505., 507515., 507525., 507535., 507545., 507555.,
           507565., 507575., 507585., 507595., 507605., 507615., 507625.,
           507635., 507645., 507655., 507665., 507675., 507685., 507695.,
           507705., 507715., 507725., 507735., 507745., 507755., 507765.,
           507775., 507785., 507795., 507805., 507815., 507825., 507835.,
           507845., 507855., 507865., 507875., 507885., 507895., 507905.,
           507915., 507925., 507935., 507945., 507955., 507965., 507975.,
           507985., 507995.])

  • spatial_ref
    ()
    int64
    0

    array(0)

• Indexes: (2)
  • x
    y
    RasterIndex (crs=EPSG:32618)

    RasterIndex(crs=EPSG:32618)
        AxisAffineTransformIndex(AxisAffineTransform(a=10, b=0, c=5e+05, d=0, e=-10, f=6e+06, axis=X, dim='x'))
        AxisAffineTransformIndex(AxisAffineTransform(a=10, b=0, c=5e+05, d=0, e=-10, f=6e+06, axis=Y, dim='y'))

  • spatial_ref
    CRSIndex (crs=EPSG:32618)

    CRSIndex
    <Projected CRS: EPSG:32618>
    Name: WGS 84 / UTM zone 18N
    Axis Info [cartesian]:
    - E[east]: Easting (metre)
    - N[north]: Northing (metre)
    Area of Use:
    - name: Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.
    - bounds: (-78.0, 0.0, -72.0, 84.0)
    Coordinate Operation:
    - name: UTM zone 18N
    - method: Transverse Mercator
    Datum: World Geodetic System 1984 ensemble
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

In [8]:

da.plot(figsize=(12, 6), cmap="viridis", add_colorbar=True);

da.plot(figsize=(12, 6), cmap="viridis", add_colorbar=True);

6. Overview selection at 20m¶

Merge at 20m to trigger the 2x overview. merge() automatically selects the coarsest overview whose resolution is still <= the target.

In [11]:

result_20m, result_spatial_20m, result_proj_20m, time_coords = merge(
    store=ob_store,
    crs="EPSG:32618",
    bbox=(500000.0, 5990000.0, 508000.0, 6000000.0),
    resolution=20.0,
    chunk_size=(128, 128),
    bands="data",
    datafusion=True,
    sortby="id",
)


affine_20m = Affine(*result_spatial_20m.transform)
da_20m = xr.DataArray(result_20m, dims=["y", "x"])

raster_idx_20m = RasterIndex.from_transform(
    affine=affine_20m,
    width=da_20m.sizes["x"],
    height=da_20m.sizes["y"],
    x_dim="x",
    y_dim="y",
    crs=result_proj_20m.code,
)
coords_20m = xr.Coordinates.from_xindex(raster_idx_20m)
da_20m = da_20m.assign_coords(coords_20m)
da_20m = da_20m.proj.assign_crs(spatial_ref=result_proj_20m.code, allow_override=True)

da_20m

result_20m, result_spatial_20m, result_proj_20m, time_coords = merge( store=ob_store, crs="EPSG:32618", bbox=(500000.0, 5990000.0, 508000.0, 6000000.0), resolution=20.0, chunk_size=(128, 128), bands="data", datafusion=True, sortby="id", ) affine_20m = Affine(*result_spatial_20m.transform) da_20m = xr.DataArray(result_20m, dims=["y", "x"]) raster_idx_20m = RasterIndex.from_transform( affine=affine_20m, width=da_20m.sizes["x"], height=da_20m.sizes["y"], x_dim="x", y_dim="y", crs=result_proj_20m.code, ) coords_20m = xr.Coordinates.from_xindex(raster_idx_20m) da_20m = da_20m.assign_coords(coords_20m) da_20m = da_20m.proj.assign_crs(spatial_ref=result_proj_20m.code, allow_override=True) da_20m

Out[11]:

<xarray.DataArray 'array-009' (y: 500, x: 400)> Size: 800kB
cubed.Array<array-009, shape=(500, 400), dtype=float32, chunks=((128, 128, 128, 116), (128, 128, 128, 16))>
Coordinates:
  * y            (y) float64 4kB 6e+06 6e+06 6e+06 ... 5.99e+06 5.99e+06
  * x            (x) float64 3kB 5e+05 5e+05 5e+05 ... 5.08e+05 5.08e+05
  * spatial_ref  int64 8B 0
Indexes:
  ┌ x            RasterIndex (crs=EPSG:32618)
  └ y
    spatial_ref  CRSIndex (crs=EPSG:32618)

xarray.DataArray

'array-009'

• y: 500
• x: 400

• cubed.Array<chunksize=(128, 128)>

  Array Chunk Bytes 800.0 KB 65.5 KB Shape (500, 400) (128, 128) Count 3 arrays in Plan 16 Chunks Type float32 np.ndarray

• Coordinates: (3)
  • y
    (y)
    float64
    6e+06 6e+06 ... 5.99e+06 5.99e+06

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    [500 values with dtype=float64]

  • x
    (x)
    float64
    5e+05 5e+05 ... 5.08e+05 5.08e+05

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

    [400 values with dtype=float64]

  • spatial_ref
    ()
    int64
    0

    array(0)

• Indexes: (2)
  • x
    y
    RasterIndex (crs=EPSG:32618)

    RasterIndex(crs=EPSG:32618)
        AxisAffineTransformIndex(AxisAffineTransform(a=20, b=0, c=5e+05, d=0, e=-20, f=6e+06, axis=X, dim='x'))
        AxisAffineTransformIndex(AxisAffineTransform(a=20, b=0, c=5e+05, d=0, e=-20, f=6e+06, axis=Y, dim='y'))

  • spatial_ref
    CRSIndex (crs=EPSG:32618)

    CRSIndex
    <Projected CRS: EPSG:32618>
    Name: WGS 84 / UTM zone 18N
    Axis Info [cartesian]:
    - E[east]: Easting (metre)
    - N[north]: Northing (metre)
    Area of Use:
    - name: Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.
    - bounds: (-78.0, 0.0, -72.0, 84.0)
    Coordinate Operation:
    - name: UTM zone 18N
    - method: Transverse Mercator
    Datum: World Geodetic System 1984 ensemble
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

In [12]:

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5))

da.plot(ax=ax1, cmap="viridis", vmin=100, vmax=400, add_colorbar=False)
ax1.set_title(f"10m — base ({da.sizes['y']}x{da.sizes['x']})")
ax1.set_aspect("equal")

da_20m.plot(ax=ax2, cmap="viridis", vmin=100, vmax=400, add_colorbar=False)
ax2.set_title(f"20m — overview ({da_20m.sizes['y']}x{da_20m.sizes['x']})")
ax2.set_aspect("equal")

fig.suptitle("Overview selection: 10m uses base, 20m uses 2x overview", fontsize=14)
fig.tight_layout();

import matplotlib.pyplot as plt fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5)) da.plot(ax=ax1, cmap="viridis", vmin=100, vmax=400, add_colorbar=False) ax1.set_title(f"10m — base ({da.sizes['y']}x{da.sizes['x']})") ax1.set_aspect("equal") da_20m.plot(ax=ax2, cmap="viridis", vmin=100, vmax=400, add_colorbar=False) ax2.set_title(f"20m — overview ({da_20m.sizes['y']}x{da_20m.sizes['x']})") ax2.set_aspect("equal") fig.suptitle("Overview selection: 10m uses base, 20m uses 2x overview", fontsize=14) fig.tight_layout();

/Users/seanharkins/projects/lazymerge/lazymerge/sources.py:183: RuntimeWarning: Successfully reconstructed a store defined in another Python module. Connection pooling will not be shared across store instances.
  zarr_table = await ZarrTable.from_obstore(store, "/meta")