Reindex via indexer

docs/how_to/examples.md

@@ -15,6 +15,10 @@ This notebook accompanies the [_Old format, no problem!: Cloud-optimizing the GO
 
 1. [Ingesting the GOES-16 archive](https://github.com/zarr-developers/VirtualiZarr/blob/main/examples/V2/goes-16-ingest.ipynb) - End-to-end notebook building the virtual store from the netCDF archive.
 
+### ITS_LIVE glacier-velocity mosaic
+
+1. [Mosaicking ITS_LIVE granules into a virtual cube](https://github.com/zarr-developers/VirtualiZarr/blob/main/examples/V2/its_live.ipynb) - Aligns granules that share a global grid but cover different regions using native xarray `concat(..., join="outer")`, stacks them along `time`, and writes the sparse virtual cube to Icechunk.
+
 ## V1 Examples
 
 !!! note

docs/how_to/validation.md

@@ -231,7 +231,6 @@ VirtualiZarr provides some utilities for creating virtual datasets with such spa
 
 The `ManifestArray.with_fill_value_only` method returns a new ManifestArray with the same schema (shape, chunks, codecs, dimension names, attributes) as a given ManifestArray, but with an empty chunk manifest and the given `fill_value`.
 This is useful for filling in missing files/variables, by creating virtual datasets containing manifestarrays which have no chunk references.
-These empty manifestarrays can then be concatenated with the manifestarrays containing chunk references, to create a virtual dataset which a logical grid spanning regions with both present and missing data.
 
 ```python exec="on" session="usage" source="material-block" result="code"
 # a virtual variable from a file where the data is present
@@ -255,6 +254,16 @@ print(combined)
 The combined variable is still a single `ManifestArray` of shape `(2, 3, 4)`, but only the first timestep references real chunks.
 Reading data from the second timestep would return the array's `fill_value`.
 
+In the previously mentioned case when a dataset is "inherently sparse at a global level, whilst still being dense at a regional level", when users want to concatenate these sparse datasets they can use normal xarray concatenation patterns with outer join.
+
+```
+xr.concat([ds_1, ds_2], dim="time", join="outer")
+```
+This will internally align the sparse datasets on a larger logical grid spanning both regions with null chunks for the sparse regions where no source chunks exist prior to concatenating them.
+
+![Aligning two regionally-dense but globally-sparse datasets onto a shared logical grid before concatenation along time](grid_alignment.png)
+
+
 ### Inhomogeneous Codecs
 
 Zarr arrays are by definition made up of a grid of chunks with identical codecs (e.g. for compression).

examples/V2/README.md

@@ -20,6 +20,10 @@ Virtualizes the same GOES-16 file using `CachingReadableStore` and `SplittingRea
 uv run --script examples/V2/goes_with_caching_stores.py
 ```
 
+### [its_live.ipynb](its_live.ipynb)
+
+Mosaics two [ITS_LIVE](https://its-live.jpl.nasa.gov/) glacier-velocity granules that share a global grid but each cover only part of it. Uses native xarray `concat(..., join="outer")` to align them onto the shared grid (sparse, with fill where a granule has no data) and stack them along `time`, keeping the data variables virtual (`ManifestArray`) throughout, then writes the result to Icechunk without copying any pixel data.
+
 ## Performance Comparison
 
 Run both scripts and compare the "Virtualization time" printed at the end of each:

examples/V2/its_live.ipynb

@@ -0,0 +1,288 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "its-live-intro",
+   "metadata": {},
+   "source": [
+    "# Mosaicking ITS\\_LIVE granules into a virtual cube\n",
+    "\n",
+    "This notebook virtualizes two [ITS\\_LIVE](https://its-live.jpl.nasa.gov/) glacier-velocity\n",
+    "granules that share a global grid but each cover only part of it, then uses native xarray\n",
+    "`concat(..., join=\"outer\")` to align them onto the shared grid (sparse, with fill where a\n",
+    "granule has no data) and stack them along `time`. The data variables stay virtual\n",
+    "(`ManifestArray`) throughout — only the coordinates are loaded — and the result is written to\n",
+    "Icechunk without copying any pixel data.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "adc1e1cd-7d0f-4345-a6b0-54a70ab044de",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "\n",
+    "import numpy as np\n",
+    "import obstore\n",
+    "import xarray as xr\n",
+    "from obspec_utils.registry import ObjectStoreRegistry\n",
+    "\n",
+    "import virtualizarr as vz\n",
+    "from virtualizarr.parsers import HDFParser"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f5cf7b30-577a-46c3-ac48-a6f62ba1ea86",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bucket = \"s3://its-live-data\"\n",
+    "key = \"test-space/virtual-cubes\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "356d1624-041a-4928-a45c-ae0070ce0b56",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "granule_1 = \"LC08_L1GT_020121_20231013_20231102_02_T2_X_LC09_L1GT_020121_20231106_20231106_02_T2_G0120V02_P084.nc\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "42bca045-986d-4d17-834b-7526ac8afabb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "store = obstore.store.from_url(bucket, region=\"us-west-2\", skip_signature=True)\n",
+    "registry = ObjectStoreRegistry({bucket: store})\n",
+    "parser = HDFParser(drop_variables=[\"mapping\", \"img_pair_info\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ce09e8ef-8959-459f-8c1f-4676cb21a735",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vds_1 = vz.open_virtual_dataset(\n",
+    "    url=os.path.join(bucket, key, granule_1),\n",
+    "    parser=parser,\n",
+    "    registry=registry,\n",
+    "    # load x/y/time as real, indexed coordinates so xarray has an index on every\n",
+    "    # dimension for alignment. Data variables stay virtual (ManifestArray).\n",
+    "    loadable_variables=[\"time\", \"y\", \"x\"],\n",
+    "    decode_times=True,\n",
+    ")\n",
+    "vds_1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "884dffd7-f125-4c1f-add0-25075fc7ce14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "granule_2 = \"LC08_L1GT_020120_20201121_20210315_02_T2_X_LC08_L1GT_020120_20210124_20210305_02_T2_G0120V02_P051.nc\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2198e8ee-62fc-431d-986b-b4f59a9a4667",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vds_2 = vz.open_virtual_dataset(\n",
+    "    url=os.path.join(bucket, key, granule_2),\n",
+    "    parser=parser,\n",
+    "    registry=registry,\n",
+    "    loadable_variables=[\"time\", \"y\", \"x\"],\n",
+    "    decode_times=True,\n",
+    ")\n",
+    "vds_2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9e7d4099-bbfc-416b-b0e4-d8eb1bfb0b57",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The ITS_LIVE data uses a descending y coordinate.  Internally, xarray's alignment machinery requires ascending\n",
+    "# ordering so the y coordinate needs to be flipped prior to concatenation.\n",
+    "flip = lambda ds: ds.assign_coords(y=-ds.y)  # descending y  ->  ascending -y\n",
+    "flip_1 = flip(vds_1)\n",
+    "flip_2 = flip(vds_2)\n",
+    "cube = xr.concat(\n",
+    "    [flip_1, flip_2], dim=\"time\", join=\"outer\", combine_attrs=\"drop_conflicts\"\n",
+    ")\n",
+    "cube = flip(cube)\n",
+    "cube"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1d2b78bc-9560-4a0a-b55c-21933e11ff21",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import shutil\n",
+    "\n",
+    "import icechunk as ic\n",
+    "\n",
+    "# the referenced chunk data lives on this public, anonymous S3 bucket\n",
+    "url_prefix = \"s3://its-live-data/\"\n",
+    "store_path = \"its_live_cube.icechunk\"\n",
+    "\n",
+    "# start fresh so this cell is re-runnable\n",
+    "shutil.rmtree(store_path, ignore_errors=True)\n",
+    "\n",
+    "# register a virtual chunk container so icechunk knows how to read the s3 refs\n",
+    "config = ic.RepositoryConfig.default()\n",
+    "config.set_virtual_chunk_container(\n",
+    "    ic.VirtualChunkContainer(\n",
+    "        url_prefix, ic.s3_store(region=\"us-west-2\", anonymous=True)\n",
+    "    )\n",
+    ")\n",
+    "repo = ic.Repository.create(\n",
+    "    storage=ic.local_filesystem_storage(store_path),\n",
+    "    config=config,\n",
+    "    authorize_virtual_chunk_access=ic.containers_credentials(\n",
+    "        {url_prefix: ic.s3_credentials(anonymous=True)}\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "\n",
+    "# icechunk metadata is strict JSON and cannot represent NaN/inf, but ITS_LIVE\n",
+    "# stores NaN in some attributes (e.g. 'stable_shift_stationary'); drop those.\n",
+    "def _drop_nonfinite_attrs(ds):\n",
+    "    ds = ds.copy()\n",
+    "\n",
+    "    def clean(attrs):\n",
+    "        keep = {}\n",
+    "        for k, v in attrs.items():\n",
+    "            arr = np.asarray(v)\n",
+    "            if arr.dtype.kind == \"f\" and not np.isfinite(arr).all():\n",
+    "                continue\n",
+    "            keep[k] = v\n",
+    "        return keep\n",
+    "\n",
+    "    ds.attrs = clean(ds.attrs)\n",
+    "    for var in ds.variables.values():\n",
+    "        var.attrs = clean(var.attrs)\n",
+    "    return ds\n",
+    "\n",
+    "\n",
+    "# write the cube: only the virtual references are stored, no pixel data is copied\n",
+    "session = repo.writable_session(\"main\")\n",
+    "_drop_nonfinite_attrs(cube).vz.to_icechunk(session.store)\n",
+    "snapshot_id = session.commit(\n",
+    "    \"its_live virtual cube: 2 granules mosaicked + stacked on time\"\n",
+    ")\n",
+    "print(\"committed snapshot\", snapshot_id)\n",
+    "\n",
+    "# reopen from the committed store\n",
+    "cube_roundtrip = xr.open_zarr(\n",
+    "    repo.readonly_session(\"main\").store, consolidated=False, zarr_format=3\n",
+    ")\n",
+    "cube_roundtrip"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6aa95554-5a0c-481c-bc2a-945576737e00",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "v = cube_roundtrip[\"v\"]\n",
+    "dv = v.isel(time=1) - v.isel(time=0)\n",
+    "dv"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fd10354b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# surface speed `v` at each time step, on a shared color scale\n",
+    "v = cube_roundtrip[\"v\"].load()\n",
+    "vmax = float(np.nanpercentile(v, 98))\n",
+    "\n",
+    "nt = v.sizes[\"time\"]\n",
+    "fig, axes = plt.subplots(1, nt, figsize=(7 * nt, 6), constrained_layout=True)\n",
+    "axes = np.atleast_1d(axes)\n",
+    "for ax, t in zip(axes, range(nt)):\n",
+    "    im = v.isel(time=t).plot.imshow(\n",
+    "        ax=ax, cmap=\"viridis\", vmin=0, vmax=vmax, add_colorbar=False\n",
+    "    )\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.set_title(str(v[\"time\"].values[t])[:10])\n",
+    "fig.colorbar(im, ax=list(axes), label=\"speed (m/yr)\", shrink=0.8)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1cc2f9b4-ce0b-4614-8b91-e2416390434e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "lim = np.nanpercentile(np.abs(dv), 98)\n",
+    "\n",
+    "# plot\n",
+    "fig, ax = plt.subplots(figsize=(9, 7))\n",
+    "dv.plot.imshow(\n",
+    "    ax=ax,\n",
+    "    cmap=\"RdBu_r\",  # red = speedup, blue = slowdown\n",
+    "    vmin=-lim,\n",
+    "    vmax=lim,  # symmetric, centered on zero\n",
+    "    cbar_kwargs={\"label\": \"Δ speed (m/yr)\"},\n",
+    ")\n",
+    "ax.set_aspect(\"equal\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}