From 249ecd6ac260d1ea91ba59ac7325fa81f8ad19bf Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 00:44:11 +0000 Subject: [PATCH 01/17] phase 1 of issue #38 --- src/zagg/viz/__init__.py | 65 +++++++ src/zagg/viz/shardmap.py | 375 +++++++++++++++++++++++++++++++++++++++ tests/test_viz.py | 227 ++++++++++++++++++++++++ 3 files changed, 667 insertions(+) create mode 100644 src/zagg/viz/__init__.py create mode 100644 src/zagg/viz/shardmap.py create mode 100644 tests/test_viz.py diff --git a/src/zagg/viz/__init__.py b/src/zagg/viz/__init__.py new file mode 100644 index 0000000..251dd2d --- /dev/null +++ b/src/zagg/viz/__init__.py @@ -0,0 +1,65 @@ +"""Shard-map viewer (issue #38). + +Two layers: + +- :mod:`zagg.viz.shardmap` -- the **headless render core** (phase 1). Pure + Python, no browser or ipyleaflet import: it turns a :class:`ShardMap` (and an + optional :class:`~zagg.catalog.sources.Catalog`) into GeoJSON + ``FeatureCollection`` dicts -- shard/chunk outlines, granule footprints, and + viewport-clipped cell outlines at the shard order. Fully unit-testable with + no widget stack installed. +- :func:`show_shardmap` -- the **ipyleaflet wrapper** (phase 2). Builds an + interactive map from a saved ShardMap (basemap + shard layer + a toggleable + footprint layer + a zoom-thresholded grid layer). ``ipyleaflet`` is imported + lazily inside the widget functions so phase-1 core and the test suite never + require it; it lives in the optional ``viz`` extra (``pip install zagg[viz]``). + +Both inputs are reused directly off the existing surface -- ``ShardMap`` / +``Catalog`` round-trips and per-grid ``shard_footprint`` / ``signature`` -- so +there is no viewer-specific file type or second tessellation (issue #38). +""" +from __future__ import annotations + +from zagg.viz.shardmap import ( + granule_footprints, + grid_from_signature, + render_shardmap, + shard_outlines, + viewport_cells, +) + + +def show_shardmap(shardmap_path, catalog=None, **kwargs): + """Build an interactive ipyleaflet map for a saved ShardMap (phase 2). + + Thin lazy passthrough to :func:`zagg.viz.leaflet.show_shardmap` so importing + :mod:`zagg.viz` never pulls in ``ipyleaflet`` -- only calling this does. + Install the widget stack with ``pip install zagg[viz]``. + + Parameters + ---------- + shardmap_path : str + Path to a ``ShardMap`` JSON file (``ShardMap.to_json``). + catalog : str or Catalog, optional + A geoparquet path or a loaded ``Catalog`` for the granule-footprint + layer. When omitted, only the shard layer is drawn. + **kwargs + Forwarded to :func:`zagg.viz.leaflet.show_shardmap`. + + Returns + ------- + ipyleaflet.Map + """ + from zagg.viz.leaflet import show_shardmap as _show + + return _show(shardmap_path, catalog=catalog, **kwargs) + + +__all__ = [ + "render_shardmap", + "shard_outlines", + "granule_footprints", + "viewport_cells", + "grid_from_signature", + "show_shardmap", +] diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py new file mode 100644 index 0000000..97f2fbd --- /dev/null +++ b/src/zagg/viz/shardmap.py @@ -0,0 +1,375 @@ +"""Headless render core for the shard-map viewer (issue #38, phase 1). + +Pure Python: turns a :class:`~zagg.catalog.shardmap.ShardMap` (and an optional +:class:`~zagg.catalog.sources.Catalog`) into GeoJSON ``FeatureCollection`` dicts +in WGS84. No browser, no ipyleaflet -- everything here is unit-testable with +just the core deps (``shapely`` + the grid backends). + +Three layers are produced: + +- :func:`shard_outlines` -- one polygon feature per shard, straight off + ``grid.shard_footprint(key)``. The grid is reconstructed from the map's own + ``grid_signature`` (:func:`grid_from_signature`) so no second grid spec is + needed. +- :func:`granule_footprints` -- one polygon feature per granule footprint, + decoded from a ``Catalog`` (``granule_records``). +- :func:`viewport_cells` -- shard-order cell outlines clipped to a viewport + bbox, emitted **only** when ``<= max_shards`` shards intersect the viewport + (the "grid-on-zoom" gate -- never a global graticule, issue #38). + +Antimeridian handling +--------------------- +HEALPix shard polygons near +-180 deg come back from mortie's ``mort2polygon`` +with longitudes that, read as a flat ring, span more than a hemisphere (e.g. +180 -> -178) and would render as a band wrapping the whole globe. The mortie +path already normalizes vertices that merely *touch* the antimeridian; for the +ones that genuinely *cross* it, :func:`_split_antimeridian` cuts the ring at ++-180 into a ``MultiPolygon`` so GeoJSON consumers draw it correctly. +""" +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np + +# Longitude span (deg) above which a ring is treated as antimeridian-crossing. +_ANTIMERIDIAN_SPAN = 180.0 + + +def grid_from_signature(signature: dict): + """Reconstruct an output grid from a ``ShardMap.grid_signature``. + + The viewer only needs ``shard_footprint`` / ``children`` off the grid, both + of which are fully determined by the signature -- so the map is + self-describing and no separate config is required. + + Parameters + ---------- + signature : dict + A grid ``signature()`` dict (``type`` is ``"healpix"`` or + ``"rectilinear"``). + + Returns + ------- + OutputGrid + + Raises + ------ + ValueError + If ``signature['type']`` is unknown. + """ + gtype = signature.get("type") + if gtype == "healpix": + from zagg.grids import HealpixGrid + + return HealpixGrid( + parent_order=signature["parent_order"], + child_order=signature["child_order"], + layout=signature.get("layout", "fullsphere"), + ) + if gtype == "rectilinear": + from zagg.grids import RectilinearGrid + + a, _b, c, _d, e, f = signature["affine"] + height, width = signature["shape"] + res_x, res_y = a, -e + xmin, ymax = c, f + xmax = c + a * width + ymin = f + e * height + return RectilinearGrid( + crs=signature["crs"], + resolution=(res_x, res_y), + bounds=[xmin, ymin, xmax, ymax], + chunk_shape=tuple(signature["chunk_shape"]), + ) + raise ValueError(f"unknown grid signature type: {gtype!r}") + + +# ── GeoJSON geometry helpers ───────────────────────────────────────────────── + +def _ring_coords(geom) -> list[list[float]]: + """Exterior-ring ``[[lon, lat], ...]`` for a shapely Polygon.""" + x, y = geom.exterior.coords.xy + return [[float(lon), float(lat)] for lon, lat in zip(x, y)] + + +def _split_antimeridian(geom): + """Split a Polygon that crosses +-180 deg into hemisphere-local parts. + + Returns a GeoJSON ``geometry`` dict -- a ``Polygon`` when the ring stays + within a hemisphere, or a ``MultiPolygon`` cut at the antimeridian when it + crosses. The cut unwraps the ring (western vertices shifted +360), clips + against ``[-180, 180]`` and ``[180, 540]`` half-planes, then rewraps the + eastern part back into ``[-180, 180]`` -- shapely-only, no extra deps. + """ + from shapely.geometry import Polygon, box + + ring = _ring_coords(geom) + lons = np.array([pt[0] for pt in ring]) + if lons.max() - lons.min() <= _ANTIMERIDIAN_SPAN: + return {"type": "Polygon", "coordinates": [ring]} + + # Unwrap: lift western-hemisphere vertices by +360 so the ring is monotone + # across the seam (e.g. 180, -178 -> 180, 182). + unwrapped = [[lon + 360.0 if lon < 0 else lon, lat] for lon, lat in ring] + poly = Polygon(unwrapped) + if not poly.is_valid: + poly = poly.buffer(0) + + west = poly.intersection(box(-180.0, -90.0, 180.0, 90.0)) + east = poly.intersection(box(180.0, -90.0, 540.0, 90.0)) + + parts: list = [] + for part, shift in ((west, 0.0), (east, -360.0)): + for sub in getattr(part, "geoms", [part]): + if sub.is_empty or sub.geom_type != "Polygon": + continue + x, y = sub.exterior.coords.xy + parts.append([[float(lon + shift), float(lat)] for lon, lat in zip(x, y)]) + + if not parts: + return {"type": "Polygon", "coordinates": [ring]} + return {"type": "MultiPolygon", "coordinates": [[p] for p in parts]} + + +def _polygon_geometry(geom) -> dict: + """GeoJSON geometry for a shapely Polygon/MultiPolygon, antimeridian-safe. + + Coordinates are plain ``list`` (not shapely's tuples) so the result is + canonical, ``json``-stable GeoJSON. + """ + if geom.geom_type == "MultiPolygon": + coords: list = [] + for sub in geom.geoms: + g = _split_antimeridian(sub) + if g["type"] == "Polygon": + coords.append(g["coordinates"]) + else: + coords.extend(g["coordinates"]) + return {"type": "MultiPolygon", "coordinates": coords} + if geom.geom_type == "Polygon": + return _split_antimeridian(geom) + raise ValueError(f"unsupported geometry type for GeoJSON: {geom.geom_type}") + + +def _polygonal(geom): + """Reduce an intersection result to Polygon/MultiPolygon, or None. + + A clip can degenerate to a point/line on a shared edge or yield a + ``GeometryCollection`` mixing dimensions; keep only the polygonal part. + """ + if geom.is_empty: + return None + if geom.geom_type in ("Polygon", "MultiPolygon"): + return geom + if geom.geom_type == "GeometryCollection": + from shapely.geometry import MultiPolygon + + polys = [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")] + if not polys: + return None + merged = MultiPolygon( + [p for g in polys for p in getattr(g, "geoms", [g])] + ) + return merged + return None + + +def _feature(geometry: dict, properties: dict) -> dict: + return {"type": "Feature", "geometry": geometry, "properties": properties} + + +def _collection(features: list[dict]) -> dict: + return {"type": "FeatureCollection", "features": features} + + +# ── layers ─────────────────────────────────────────────────────────────────── + +def shard_outlines(shardmap) -> dict: + """Shard/chunk outlines as a GeoJSON ``FeatureCollection``. + + One feature per shard key in the map, with the shard's footprint polygon and + its granule count under ``properties``. The grid is reconstructed from the + map's ``grid_signature`` -- no separate grid spec needed. + + Parameters + ---------- + shardmap : ShardMap + A built or loaded shard map. + + Returns + ------- + dict + GeoJSON ``FeatureCollection`` (WGS84). + """ + grid = grid_from_signature(shardmap.grid_signature) + features = [] + for key, granules in zip(shardmap.shard_keys, shardmap.granules): + geom = grid.shard_footprint(key) + features.append( + _feature( + _polygon_geometry(geom), + {"shard_key": _jsonable(key), "n_granules": len(granules)}, + ) + ) + return _collection(features) + + +def granule_footprints(catalog) -> dict: + """Granule footprints as a GeoJSON ``FeatureCollection``. + + One polygon feature per granule in the catalog (its footprint + exterior ring), with the granule id under ``properties``. + + Parameters + ---------- + catalog : Catalog + A loaded catalog (provides ``granule_records``). + + Returns + ------- + dict + GeoJSON ``FeatureCollection`` (WGS84). + """ + from shapely.geometry import Polygon + + features = [] + for rec in catalog.granule_records(): + ring = list(zip(np.asarray(rec["lons"]), np.asarray(rec["lats"]))) + if len(ring) < 4: + continue + features.append( + _feature(_polygon_geometry(Polygon(ring)), {"id": rec["id"]}) + ) + return _collection(features) + + +def viewport_cells(shardmap, bbox, *, max_shards: int = 4) -> dict: + """Shard-order cell outlines clipped to a viewport, gated on visible shards. + + Implements the "grid-on-zoom" behavior (issue #38): cell outlines **at the + shard order** are drawn only when ``<= max_shards`` shards intersect + ``bbox``. When more shards are visible the viewport is too zoomed-out for a + useful grid and an empty collection is returned -- never a global graticule. + + Parameters + ---------- + shardmap : ShardMap + A built or loaded shard map. + bbox : tuple of float + Viewport ``(lon_min, lat_min, lon_max, lat_max)`` in WGS84. + max_shards : int + Maximum number of intersecting shards for the grid to render. + + Returns + ------- + dict + GeoJSON ``FeatureCollection`` of shard-cell outlines clipped to the + viewport, or an empty collection when the gate is not met. + """ + from shapely.geometry import box + + grid = grid_from_signature(shardmap.grid_signature) + view = box(bbox[0], bbox[1], bbox[2], bbox[3]) + + visible = [ + key + for key in shardmap.shard_keys + if grid.shard_footprint(key).intersects(view) + ] + if not visible or len(visible) > max_shards: + return _collection([]) + + features = [] + for key in visible: + geom = grid.shard_footprint(key) + clipped = _polygonal(geom.intersection(view)) + if clipped is None: + continue + features.append( + _feature(_polygon_geometry(clipped), {"shard_key": _jsonable(key)}) + ) + return _collection(features) + + +# ── top-level assembly ─────────────────────────────────────────────────────── + +def render_shardmap(shardmap, catalog=None, *, bbox=None, max_shards: int = 4) -> dict: + """Assemble all viewer layers for a shard map into one dict of collections. + + Parameters + ---------- + shardmap : ShardMap or str + A ``ShardMap`` or a path to a ShardMap JSON file. + catalog : Catalog or str, optional + A ``Catalog`` or a geoparquet path. When given, the granule-footprint + layer is included. + bbox : tuple of float, optional + Viewport ``(lon_min, lat_min, lon_max, lat_max)``. When given, the + viewport-clipped cell-outline layer is included. + max_shards : int + Visible-shard gate for the viewport cell layer. + + Returns + ------- + dict + ``{"shards": FC, "granules": FC | None, "cells": FC | None}`` where each + value is a GeoJSON ``FeatureCollection`` (or ``None`` when its input was + not provided). + """ + shardmap = _load_shardmap(shardmap) + out = {"shards": shard_outlines(shardmap), "granules": None, "cells": None} + if catalog is not None: + out["granules"] = granule_footprints(_load_catalog(catalog)) + if bbox is not None: + out["cells"] = viewport_cells(shardmap, bbox, max_shards=max_shards) + return out + + +# ── small loaders / utils ──────────────────────────────────────────────────── + +def _load_shardmap(shardmap): + """Accept a ShardMap or a JSON path; return a ShardMap.""" + if isinstance(shardmap, (str, Path)): + from zagg.catalog.shardmap import ShardMap + + return ShardMap.from_json(str(shardmap)) + return shardmap + + +def _load_catalog(catalog): + """Accept a Catalog or a geoparquet path; return a Catalog.""" + if isinstance(catalog, (str, Path)): + from zagg.catalog.sources import Catalog + + return Catalog.from_geoparquet(str(catalog)) + return catalog + + +def _jsonable(key): + """Render a shard key (int or tuple) as a JSON-safe value.""" + if isinstance(key, (list, tuple)): + return [int(k) for k in key] + try: + return int(key) + except (TypeError, ValueError): + return key + + +def _is_geojson(obj) -> bool: + """True if ``obj`` round-trips as JSON and is a FeatureCollection.""" + return ( + isinstance(obj, dict) + and obj.get("type") == "FeatureCollection" + and json.loads(json.dumps(obj)) == obj + ) + + +__all__ = [ + "grid_from_signature", + "shard_outlines", + "granule_footprints", + "viewport_cells", + "render_shardmap", +] diff --git a/tests/test_viz.py b/tests/test_viz.py new file mode 100644 index 0000000..37dfa3d --- /dev/null +++ b/tests/test_viz.py @@ -0,0 +1,227 @@ +"""Tests for the headless shard-map render core (issue #38, phase 1). + +Pure Python -- no ipyleaflet / widget stack. Exercises GeoJSON emission off a +small saved ShardMap fixture, the catalog footprint layer, the viewport +grid-on-zoom gate, and antimeridian splitting. +""" + +import json + +import numpy as np +import pyarrow as pa +import pytest +import stac_geoparquet.arrow as sga + +from zagg.catalog.shardmap import ShardMap +from zagg.catalog.sources import Catalog +from zagg.grids import HealpixGrid, RectilinearGrid +from zagg.viz import ( + granule_footprints, + grid_from_signature, + render_shardmap, + shard_outlines, + viewport_cells, +) +from zagg.viz.shardmap import _is_geojson, _split_antimeridian + +# ── fixtures ───────────────────────────────────────────────────────────────── + +def _item(gid, lon0, lon1, lat0=38.85, lat1=38.93): + ring = [[lon0, lat0], [lon1, lat0], [lon1, lat1], [lon0, lat1], [lon0, lat0]] + return { + "type": "Feature", "stac_version": "1.0.0", "id": gid, + "geometry": {"type": "Polygon", "coordinates": [ring]}, + "bbox": [lon0, lat0, lon1, lat1], + "properties": {"datetime": "2025-06-01T00:00:00Z"}, + "collection": "TEST", "stac_extensions": [], "links": [], + "assets": { + "data": {"href": f"https://h/{gid}.h5", "roles": ["data"]}, + "data_s3": {"href": f"s3://b/{gid}.h5", "roles": ["data"]}, + }, + } + + +@pytest.fixture +def rect_grid(): + return RectilinearGrid( + "EPSG:32618", 10, [359400, 4300740, 369400, 4310740], [250, 250] + ) + + +@pytest.fixture +def shardmap(rect_grid): + """A tiny hand-built ShardMap over a 4x4 chunk grid (no fetch needed).""" + keys = [0, 1, 5] + granules = [ + [{"id": "Ga", "s3": "s3://b/a.h5", "https": "https://h/a.h5"}], + [{"id": "Gb", "s3": "s3://b/b.h5", "https": "https://h/b.h5"}], + [ + {"id": "Gb", "s3": "s3://b/b.h5", "https": "https://h/b.h5"}, + {"id": "Gc", "s3": "s3://b/c.h5", "https": "https://h/c.h5"}, + ], + ] + return ShardMap(rect_grid.signature(), keys, granules, {"backend": "test"}) + + +@pytest.fixture +def catalog(): + items = [_item("Ga", -76.62, -76.57), _item("Gb", -76.55, -76.50)] + return Catalog( + pa.table(sga.parse_stac_items_to_arrow(items)), + {"collection": "TEST"}, + ) + + +# ── grid_from_signature ────────────────────────────────────────────────────── + +class TestGridFromSignature: + def test_rectilinear_round_trip(self, rect_grid): + g = grid_from_signature(rect_grid.signature()) + assert g.signature() == rect_grid.signature() + assert g.shard_footprint(0).equals(rect_grid.shard_footprint(0)) + + def test_healpix_round_trip(self): + hp = HealpixGrid(3, 7, layout="fullsphere") + g = grid_from_signature(hp.signature()) + assert g.signature() == hp.signature() + + def test_unknown_type_raises(self): + with pytest.raises(ValueError, match="unknown grid signature type"): + grid_from_signature({"type": "mystery"}) + + +# ── shard_outlines ─────────────────────────────────────────────────────────── + +class TestShardOutlines: + def test_one_feature_per_shard(self, shardmap): + fc = shard_outlines(shardmap) + assert _is_geojson(fc) + assert len(fc["features"]) == len(shardmap.shard_keys) + + def test_properties_and_geometry(self, shardmap): + fc = shard_outlines(shardmap) + feat = fc["features"][0] + assert feat["properties"]["shard_key"] == 0 + assert feat["properties"]["n_granules"] == 1 + assert feat["geometry"]["type"] in ("Polygon", "MultiPolygon") + + def test_wgs84_coords(self, shardmap): + # UTM 18N grid near 38.9N, -76.6E -> lons in [-77, -76], lats ~[38.8, 39]. + fc = shard_outlines(shardmap) + ring = fc["features"][0]["geometry"]["coordinates"][0] + lons = [pt[0] for pt in ring] + lats = [pt[1] for pt in ring] + assert all(-78 < lon < -75 for lon in lons) + assert all(38 < lat < 40 for lat in lats) + + def test_valid_json(self, shardmap): + fc = shard_outlines(shardmap) + assert json.loads(json.dumps(fc)) == fc + + +# ── granule_footprints ─────────────────────────────────────────────────────── + +class TestGranuleFootprints: + def test_one_feature_per_granule(self, catalog): + fc = granule_footprints(catalog) + assert _is_geojson(fc) + assert len(fc["features"]) == 2 + ids = {f["properties"]["id"] for f in fc["features"]} + assert ids == {"Ga", "Gb"} + + def test_geometry_is_polygon(self, catalog): + fc = granule_footprints(catalog) + assert fc["features"][0]["geometry"]["type"] == "Polygon" + + +# ── viewport_cells (grid-on-zoom gate) ─────────────────────────────────────── + +class TestViewportCells: + def test_gate_open_few_shards(self, shardmap): + # Tight viewport over one chunk -> <= max_shards visible, grid drawn. + fp = grid_from_signature(shardmap.grid_signature).shard_footprint(0) + lon0, lat0, lon1, lat1 = fp.bounds + fc = viewport_cells(shardmap, (lon0, lat0, lon1, lat1), max_shards=4) + assert _is_geojson(fc) + assert len(fc["features"]) >= 1 + + def test_gate_closed_too_many_shards(self, shardmap): + # Global viewport over all shards but max_shards=1 -> gated empty. + fc = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=1) + assert fc["features"] == [] + + def test_gate_closed_no_shards(self, shardmap): + # Viewport far from the grid -> nothing visible -> empty. + fc = viewport_cells(shardmap, (10, 10, 11, 11), max_shards=4) + assert fc["features"] == [] + + def test_clipped_to_viewport(self, shardmap): + grid = grid_from_signature(shardmap.grid_signature) + fp = grid.shard_footprint(0) + lon0, lat0, lon1, lat1 = fp.bounds + # Half-width viewport -> clipped cell stays within the viewport bbox. + view = (lon0, lat0, (lon0 + lon1) / 2, lat1) + fc = viewport_cells(shardmap, view, max_shards=4) + for feat in fc["features"]: + for ring in feat["geometry"]["coordinates"]: + for lon, lat in ring: + assert view[0] - 1e-6 <= lon <= view[2] + 1e-6 + + +# ── antimeridian splitting ─────────────────────────────────────────────────── + +class TestAntimeridian: + def test_healpix_shard_near_antimeridian_splits(self): + # A HEALPix parent cell straddling +-180 -> MultiPolygon, each part in + # one hemisphere (no globe-spanning band). + grid = HealpixGrid(2, 6, layout="fullsphere") + key = int( + grid.coverage( + [(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))] + )[0] + ) + sm = ShardMap( + grid.signature(), [key], + [[{"id": "G", "s3": "s", "https": "h"}]], {}, + ) + fc = shard_outlines(sm) + geom = fc["features"][0]["geometry"] + assert geom["type"] == "MultiPolygon" + for poly in geom["coordinates"]: + lons = [pt[0] for pt in poly[0]] + assert max(lons) - min(lons) <= 180.0 + + def test_non_crossing_polygon_stays_polygon(self): + from shapely.geometry import Polygon + + poly = Polygon([(10, 0), (12, 0), (12, 2), (10, 2), (10, 0)]) + geom = _split_antimeridian(poly) + assert geom["type"] == "Polygon" + + +# ── render_shardmap assembly ───────────────────────────────────────────────── + +class TestRenderShardmap: + def test_shards_only(self, shardmap): + out = render_shardmap(shardmap) + assert _is_geojson(out["shards"]) + assert out["granules"] is None + assert out["cells"] is None + + def test_with_catalog_and_bbox(self, shardmap, catalog): + out = render_shardmap(shardmap, catalog, bbox=(-180, -90, 180, 90)) + assert _is_geojson(out["shards"]) + assert _is_geojson(out["granules"]) + assert _is_geojson(out["cells"]) + + def test_from_json_path(self, shardmap, tmp_path): + path = tmp_path / "sm.json" + shardmap.to_json(str(path)) + out = render_shardmap(str(path)) + assert len(out["shards"]["features"]) == len(shardmap.shard_keys) + + def test_from_geoparquet_path(self, shardmap, catalog, tmp_path): + cat_path = tmp_path / "cat.parquet" + catalog.to_geoparquet(str(cat_path)) + out = render_shardmap(shardmap, str(cat_path)) + assert len(out["granules"]["features"]) == 2 From ae8159590a28811086143219c3ec2bb2631e9fe8 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 00:50:05 +0000 Subject: [PATCH 02/17] fold review: antimeridian seam-jump detection, keep holes (#38) --- src/zagg/viz/shardmap.py | 72 ++++++++++++++++++++++++++-------------- tests/test_viz.py | 26 +++++++++++++++ 2 files changed, 74 insertions(+), 24 deletions(-) diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py index 97f2fbd..ec18c52 100644 --- a/src/zagg/viz/shardmap.py +++ b/src/zagg/viz/shardmap.py @@ -88,32 +88,51 @@ def grid_from_signature(signature: dict): # ── GeoJSON geometry helpers ───────────────────────────────────────────────── +def _ring_list(ring) -> list[list[float]]: + """A shapely ring's ``[[lon, lat], ...]`` as plain floats.""" + x, y = ring.coords.xy + return [[float(lon), float(lat)] for lon, lat in zip(x, y)] + + def _ring_coords(geom) -> list[list[float]]: """Exterior-ring ``[[lon, lat], ...]`` for a shapely Polygon.""" - x, y = geom.exterior.coords.xy - return [[float(lon), float(lat)] for lon, lat in zip(x, y)] + return _ring_list(geom.exterior) + + +def _crosses_antimeridian(geom) -> bool: + """True if the polygon's exterior ring jumps the +-180 seam. + + A *jump* -- two consecutive vertices > 180 deg apart in longitude -- means + the ring crosses the antimeridian. This is distinct from a merely wide + polygon (e.g. a swath from lon -170 to +170 across 0 deg), whose vertices + step continuously and never jump, so it is left intact (review of #38 + phase 1). + """ + lons = np.array([pt[0] for pt in _ring_coords(geom)]) + return bool(np.any(np.abs(np.diff(lons)) > _ANTIMERIDIAN_SPAN)) def _split_antimeridian(geom): """Split a Polygon that crosses +-180 deg into hemisphere-local parts. - Returns a GeoJSON ``geometry`` dict -- a ``Polygon`` when the ring stays - within a hemisphere, or a ``MultiPolygon`` cut at the antimeridian when it - crosses. The cut unwraps the ring (western vertices shifted +360), clips - against ``[-180, 180]`` and ``[180, 540]`` half-planes, then rewraps the - eastern part back into ``[-180, 180]`` -- shapely-only, no extra deps. + Returns a GeoJSON ``geometry`` dict -- a ``Polygon`` (interior rings kept) + when the ring does not cross the seam, or a ``MultiPolygon`` cut at the + antimeridian when it does. The cut unwraps the polygon (western vertices + shifted +360), clips against the ``[-180, 180]`` and ``[180, 540]`` + half-planes, then rewraps the eastern part back into ``[-180, 180]`` -- + shapely-only, no extra deps. Holes are carried through the unwrap/clip. """ from shapely.geometry import Polygon, box - ring = _ring_coords(geom) - lons = np.array([pt[0] for pt in ring]) - if lons.max() - lons.min() <= _ANTIMERIDIAN_SPAN: - return {"type": "Polygon", "coordinates": [ring]} + if not _crosses_antimeridian(geom): + return {"type": "Polygon", "coordinates": _polygon_rings(geom)} # Unwrap: lift western-hemisphere vertices by +360 so the ring is monotone - # across the seam (e.g. 180, -178 -> 180, 182). - unwrapped = [[lon + 360.0 if lon < 0 else lon, lat] for lon, lat in ring] - poly = Polygon(unwrapped) + # across the seam (e.g. 180, -178 -> 180, 182). Interiors come along. + def _unwrap(ring): + return [[lon + 360.0 if lon < 0 else lon, lat] for lon, lat in _ring_list(ring)] + + poly = Polygon(_unwrap(geom.exterior), [_unwrap(r) for r in geom.interiors]) if not poly.is_valid: poly = poly.buffer(0) @@ -125,12 +144,17 @@ def _split_antimeridian(geom): for sub in getattr(part, "geoms", [part]): if sub.is_empty or sub.geom_type != "Polygon": continue - x, y = sub.exterior.coords.xy - parts.append([[float(lon + shift), float(lat)] for lon, lat in zip(x, y)]) + parts.append(_polygon_rings(sub, shift=shift)) if not parts: - return {"type": "Polygon", "coordinates": [ring]} - return {"type": "MultiPolygon", "coordinates": [[p] for p in parts]} + return {"type": "Polygon", "coordinates": _polygon_rings(geom)} + return {"type": "MultiPolygon", "coordinates": parts} + + +def _polygon_rings(geom, *, shift: float = 0.0) -> list[list[list[float]]]: + """GeoJSON ring list ``[exterior, *interiors]`` for a Polygon, lon-shifted.""" + rings = [geom.exterior, *geom.interiors] + return [[[lon + shift, lat] for lon, lat in _ring_list(r)] for r in rings] def _polygon_geometry(geom) -> dict: @@ -273,18 +297,18 @@ def viewport_cells(shardmap, bbox, *, max_shards: int = 4) -> dict: grid = grid_from_signature(shardmap.grid_signature) view = box(bbox[0], bbox[1], bbox[2], bbox[3]) + # Build each footprint once (reproject + densify is not free), then gate. visible = [ - key - for key in shardmap.shard_keys - if grid.shard_footprint(key).intersects(view) + (key, fp) + for key, fp in ((k, grid.shard_footprint(k)) for k in shardmap.shard_keys) + if fp.intersects(view) ] if not visible or len(visible) > max_shards: return _collection([]) features = [] - for key in visible: - geom = grid.shard_footprint(key) - clipped = _polygonal(geom.intersection(view)) + for key, fp in visible: + clipped = _polygonal(fp.intersection(view)) if clipped is None: continue features.append( diff --git a/tests/test_viz.py b/tests/test_viz.py index 37dfa3d..a288786 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -198,6 +198,32 @@ def test_non_crossing_polygon_stays_polygon(self): geom = _split_antimeridian(poly) assert geom["type"] == "Polygon" + def test_wide_non_crossing_polygon_kept_intact(self): + # A swath that steps continuously from -170 across 0 to +170 spans >180 + # in total but never jumps the seam between consecutive vertices; it + # must stay a single Polygon, not get split into ±180 slivers (review of + # #38 phase 1 -- the old total-span gate over-split this). + from shapely.geometry import Polygon + + ring = [ + (-170, -10), (-85, -10), (0, -10), (85, -10), (170, -10), + (170, 10), (85, 10), (0, 10), (-85, 10), (-170, 10), (-170, -10), + ] + geom = _split_antimeridian(Polygon(ring)) + assert geom["type"] == "Polygon" + lons = [pt[0] for pt in geom["coordinates"][0]] + assert min(lons) == -170 and max(lons) == 170 + + def test_holes_preserved(self): + # An exterior ring with an interior hole keeps the hole through GeoJSON. + from shapely.geometry import Polygon + + shell = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)] + hole = [(3, 3), (7, 3), (7, 7), (3, 7), (3, 3)] + geom = _split_antimeridian(Polygon(shell, [hole])) + assert geom["type"] == "Polygon" + assert len(geom["coordinates"]) == 2 # exterior + one interior ring + # ── render_shardmap assembly ───────────────────────────────────────────────── From e7fce374dd6297cf53fd482da02909b9349834ab Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 00:50:12 +0000 Subject: [PATCH 03/17] phase 2 of issue #38 --- pyproject.toml | 7 +++ src/zagg/viz/leaflet.py | 129 ++++++++++++++++++++++++++++++++++++++++ tests/test_viz.py | 38 ++++++++++++ 3 files changed, 174 insertions(+) create mode 100644 src/zagg/viz/leaflet.py diff --git a/pyproject.toml b/pyproject.toml index ca71d94..ef0cd12 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -61,6 +61,13 @@ lambda = [ catalog = [ "stac-geoparquet>=0.7.0", ] +# Interactive shard-map viewer (issue #38). Optional — kept out of core and out +# of `lambda` so the deployment layer stays lean. `pip install zagg[viz]`. +# ipyleaflet is a Jupyter widget; @espg approved adding it for the viewer +# (https://github.com/englacial/zagg/issues/38#issuecomment-4713639466). +viz = [ + "ipyleaflet>=0.19", +] analysis = [ "cubed>=0.24.0", "cubed-xarray>=0.0.9", diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py new file mode 100644 index 0000000..e85c275 --- /dev/null +++ b/src/zagg/viz/leaflet.py @@ -0,0 +1,129 @@ +"""ipyleaflet wrapper for the shard-map viewer (issue #38, phase 2). + +Builds an interactive map from a saved ShardMap: a basemap, the shard-outline +layer, an optional (toggleable) granule-footprint layer, and a grid layer that +draws shard-order cell outlines only when the viewport is zoomed in far enough +to show ``<= max_shards`` shards (the "grid-on-zoom" gate -- never a global +graticule, issue #38). + +All ``ipyleaflet`` imports are local to the functions here so importing +:mod:`zagg.viz` (and the phase-1 render core / test suite) never requires the +widget stack. Install it with ``pip install zagg[viz]``. + +The geometry is produced entirely by the headless core in +:mod:`zagg.viz.shardmap` -- this module only wires those GeoJSON collections +onto ipyleaflet layers and keeps the grid layer in sync with the viewport. +""" +from __future__ import annotations + +from zagg.viz.shardmap import ( + _load_catalog, + _load_shardmap, + granule_footprints, + shard_outlines, + viewport_cells, +) + +# Layer styles (kept terse; tweakable by callers via the returned Map). +_SHARD_STYLE = {"color": "#1f78b4", "weight": 1, "fillOpacity": 0.05} +_FOOTPRINT_STYLE = {"color": "#e31a1c", "weight": 1, "fillOpacity": 0.10} +_GRID_STYLE = {"color": "#333333", "weight": 1, "fillOpacity": 0.0} + + +def _center_zoom(fc: dict): + """Center ``(lat, lon)`` for a FeatureCollection's bbox (zoom is left default).""" + lons: list[float] = [] + lats: list[float] = [] + for feat in fc["features"]: + _walk_coords(feat["geometry"]["coordinates"], lons, lats) + if not lons: + return (0.0, 0.0) + return ((min(lats) + max(lats)) / 2, (min(lons) + max(lons)) / 2) + + +def _walk_coords(coords, lons, lats): + """Collect lon/lat from an arbitrarily nested GeoJSON coordinate array.""" + if coords and isinstance(coords[0], (int, float)): + lons.append(coords[0]) + lats.append(coords[1]) + return + for sub in coords: + _walk_coords(sub, lons, lats) + + +def show_shardmap( + shardmap_path, + catalog=None, + *, + max_shards: int = 4, + zoom: int = 3, + basemap=None, +): + """Build an interactive ipyleaflet map for a saved ShardMap. + + Parameters + ---------- + shardmap_path : str or ShardMap + Path to a ``ShardMap`` JSON file (or an in-memory ``ShardMap``). + catalog : str or Catalog, optional + A geoparquet path or a loaded ``Catalog``. When given, a toggleable + granule-footprint layer is added (off by default in the layer control). + max_shards : int + Visible-shard gate for the grid-on-zoom layer. + zoom : int + Initial map zoom. + basemap : ipyleaflet basemap, optional + Overrides the default OpenStreetMap basemap. + + Returns + ------- + ipyleaflet.Map + Map with a shard layer, optional footprint layer, a zoom-thresholded + grid layer, and a ``LayersControl`` for toggling layers. + """ + # Import first so a missing `viz` extra fails clearly, before any work. + from ipyleaflet import GeoJSON, LayersControl, Map, basemaps + + shardmap = _load_shardmap(shardmap_path) + shards_fc = shard_outlines(shardmap) + + center = _center_zoom(shards_fc) + m = Map(basemap=basemap or basemaps.OpenStreetMap.Mapnik, center=center, zoom=zoom) + + shard_layer = GeoJSON(data=shards_fc, style=_SHARD_STYLE, name="shards") + m.add(shard_layer) + + if catalog is not None: + footprint_fc = granule_footprints(_load_catalog(catalog)) + footprint_layer = GeoJSON( + data=footprint_fc, style=_FOOTPRINT_STYLE, name="granule footprints" + ) + m.add(footprint_layer) + + # Grid-on-zoom: an empty layer kept in sync with the viewport. Recomputed on + # every bounds change; the headless core's gate returns nothing when too + # many shards are visible, so this never becomes a global graticule. + grid_layer = GeoJSON( + data={"type": "FeatureCollection", "features": []}, + style=_GRID_STYLE, + name="grid (shard cells)", + ) + m.add(grid_layer) + + def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) + bounds = m.bounds + if not bounds: + return + (south, west), (north, east) = bounds + grid_layer.data = viewport_cells( + shardmap, (west, south, east, north), max_shards=max_shards + ) + + m.observe(_refresh_grid, names="bounds") + _refresh_grid() + + m.add(LayersControl(position="topright")) + return m + + +__all__ = ["show_shardmap"] diff --git a/tests/test_viz.py b/tests/test_viz.py index a288786..7f7531e 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -251,3 +251,41 @@ def test_from_geoparquet_path(self, shardmap, catalog, tmp_path): catalog.to_geoparquet(str(cat_path)) out = render_shardmap(shardmap, str(cat_path)) assert len(out["granules"]["features"]) == 2 + + +# ── phase 2: ipyleaflet wrapper (skips when the viz extra isn't installed) ──── + +class TestShowShardmap: + def test_import_core_without_ipyleaflet(self): + # The headless core and zagg.viz import must not require ipyleaflet. + import zagg.viz # noqa: F401 + + assert hasattr(zagg.viz, "shard_outlines") + + def test_build_map(self, shardmap, tmp_path): + pytest.importorskip("ipyleaflet") + from ipyleaflet import Map + + from zagg.viz import show_shardmap + + path = tmp_path / "sm.json" + shardmap.to_json(str(path)) + m = show_shardmap(str(path)) + assert isinstance(m, Map) + # shard layer + grid layer (+ basemap) present. + assert len(m.layers) >= 2 + + def test_build_map_with_catalog(self, shardmap, catalog, tmp_path): + pytest.importorskip("ipyleaflet") + from ipyleaflet import Map + + from zagg.viz import show_shardmap + + sm_path = tmp_path / "sm.json" + cat_path = tmp_path / "cat.parquet" + shardmap.to_json(str(sm_path)) + catalog.to_geoparquet(str(cat_path)) + m = show_shardmap(str(sm_path), catalog=str(cat_path)) + assert isinstance(m, Map) + names = {getattr(layer, "name", None) for layer in m.layers} + assert "granule footprints" in names From 69b53012b9f75fead0587ee8a4226029b0244e9f Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 10:15:51 +0000 Subject: [PATCH 04/17] phase 3 of issue #38 --- docs/quickstart.md | 7 + notebooks/shardmap_viewer.ipynb | 334 ++++++++++++++++++++++++++++++++ 2 files changed, 341 insertions(+) create mode 100644 notebooks/shardmap_viewer.ipynb diff --git a/docs/quickstart.md b/docs/quickstart.md index 0eb8c46..6ce0e13 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -61,6 +61,13 @@ This produces a JSON file (e.g., `shardmap_ATL06_2024-01-06_2024-04-07.json`) th parent morton cells to the S3 URLs of HDF5 granules containing data for those cells. The processing step consumes this file. +To inspect the chunking interactively -- shard outlines, granule footprints, +and a grid that appears on zoom -- use the shard-map viewer +(`pip install zagg[viz]`). See the +[shard-map viewer notebook](https://github.com/englacial/zagg/blob/main/notebooks/shardmap_viewer.ipynb), +which runs on a synthetic example (no network needed) and includes manual +in-browser verification instructions. + ## Local Processing The simplest path -- no AWS Lambda needed: diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb new file mode 100644 index 0000000..5f91503 --- /dev/null +++ b/notebooks/shardmap_viewer.ipynb @@ -0,0 +1,334 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Shard-map viewer\n", + "\n", + "This notebook demonstrates and gives **verification instructions** for the\n", + "`zagg.viz` shard-map viewer (issue\n", + "[#38](https://github.com/englacial/zagg/issues/38)).\n", + "\n", + "The viewer has two layers:\n", + "\n", + "- **Headless render core** (`zagg.viz.render_shardmap` and friends) -- pure\n", + " Python, no browser/widget stack. It turns a `ShardMap` (and an optional\n", + " `Catalog` / STAC-geoparquet file) into WGS84 GeoJSON `FeatureCollection`\n", + " dicts: shard outlines, granule footprints, and viewport-clipped grid cells.\n", + "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) -- builds an interactive\n", + " map from those collections: a basemap, the shard-outline layer, a toggleable\n", + " granule-footprint layer, and a grid layer that only draws when you zoom in\n", + " far enough (the \"grid-on-zoom\" gate -- never a global graticule).\n", + "\n", + "## Install\n", + "\n", + "The headless core needs only the zagg core deps (`shapely` + the grid\n", + "backends). The interactive map needs the optional `viz` extra:\n", + "\n", + "```bash\n", + "pip install zagg[viz] # or: uv sync --extra viz\n", + "```\n", + "\n", + "Everything below uses a tiny **synthetic** ShardMap and catalog, so the\n", + "notebook runs with no network, S3, or Earthdata credentials." + ], + "id": "cell-0" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Build a small synthetic ShardMap and Catalog\n", + "\n", + "We hand-build a `ShardMap` over a 4x4-chunk rectilinear grid (UTM 18N, near\n", + "Washington DC) with three populated shards, plus a two-granule STAC `Catalog`.\n", + "This mirrors the test fixtures in `tests/test_viz.py`. In real use you would\n", + "instead load a saved shard map and catalog produced by `python -m zagg.catalog`\n", + "(see [quickstart](https://github.com/englacial/zagg/blob/main/docs/quickstart.md)):\n", + "\n", + "```python\n", + "from zagg.catalog.shardmap import ShardMap\n", + "from zagg.catalog.sources import Catalog\n", + "sm = ShardMap.from_json(\"shardmap_ATL06_cycle22.json\")\n", + "cat = Catalog.from_geoparquet(\"catalog_ATL06_cycle22.parquet\")\n", + "```" + ], + "id": "cell-1" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pyarrow as pa\n", + "import stac_geoparquet.arrow as sga\n", + "\n", + "from zagg.catalog.shardmap import ShardMap\n", + "from zagg.catalog.sources import Catalog\n", + "from zagg.grids import HealpixGrid, RectilinearGrid\n", + "\n", + "# A 10 km x 10 km rectilinear grid at 10 m, chunked 250x250 -> one chunk == one\n", + "# shard. Three shards are populated with synthetic granule references.\n", + "grid = RectilinearGrid(\n", + " \"EPSG:32618\", 10, [359400, 4300740, 369400, 4310740], [250, 250]\n", + ")\n", + "shard_keys = [0, 1, 5]\n", + "granules = [\n", + " [{\"id\": \"Ga\", \"s3\": \"s3://b/a.h5\", \"https\": \"https://h/a.h5\"}],\n", + " [{\"id\": \"Gb\", \"s3\": \"s3://b/b.h5\", \"https\": \"https://h/b.h5\"}],\n", + " [\n", + " {\"id\": \"Gb\", \"s3\": \"s3://b/b.h5\", \"https\": \"https://h/b.h5\"},\n", + " {\"id\": \"Gc\", \"s3\": \"s3://b/c.h5\", \"https\": \"https://h/c.h5\"},\n", + " ],\n", + "]\n", + "shardmap = ShardMap(grid.signature(), shard_keys, granules, {\"backend\": \"demo\"})\n", + "shardmap.shard_keys" + ], + "id": "cell-2" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# A tiny STAC catalog of two granule footprints over the same area.\n", + "def _stac_item(gid, lon0, lon1, lat0=38.85, lat1=38.93):\n", + " ring = [[lon0, lat0], [lon1, lat0], [lon1, lat1], [lon0, lat1], [lon0, lat0]]\n", + " return {\n", + " \"type\": \"Feature\", \"stac_version\": \"1.0.0\", \"id\": gid,\n", + " \"geometry\": {\"type\": \"Polygon\", \"coordinates\": [ring]},\n", + " \"bbox\": [lon0, lat0, lon1, lat1],\n", + " \"properties\": {\"datetime\": \"2025-06-01T00:00:00Z\"},\n", + " \"collection\": \"DEMO\", \"stac_extensions\": [], \"links\": [],\n", + " \"assets\": {\n", + " \"data\": {\"href\": f\"https://h/{gid}.h5\", \"roles\": [\"data\"]},\n", + " \"data_s3\": {\"href\": f\"s3://b/{gid}.h5\", \"roles\": [\"data\"]},\n", + " },\n", + " }\n", + "\n", + "\n", + "items = [_stac_item(\"Ga\", -76.62, -76.57), _stac_item(\"Gb\", -76.55, -76.50)]\n", + "catalog = Catalog(\n", + " pa.table(sga.parse_stac_items_to_arrow(items)), {\"collection\": \"DEMO\"}\n", + ")\n", + "len(list(catalog.granule_records()))" + ], + "id": "cell-3" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Headless render core: GeoJSON FeatureCollections\n", + "\n", + "`render_shardmap` assembles every layer into one dict of GeoJSON\n", + "`FeatureCollection`s. Pass a `catalog` to include granule footprints and a\n", + "`bbox` to include the viewport-clipped grid cells. Each value is plain,\n", + "JSON-serializable GeoJSON -- no widgets required, so this is the part covered\n", + "by the unit tests." + ], + "id": "cell-4" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from zagg.viz import render_shardmap\n", + "\n", + "layers = render_shardmap(shardmap, catalog)\n", + "{name: (fc and len(fc[\"features\"])) for name, fc in layers.items()}" + ], + "id": "cell-5" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# One feature per shard, with the shard key and its granule count. The grid is\n", + "# rebuilt from the map's own `grid_signature` -- the map is self-describing.\n", + "shards_fc = layers[\"shards\"]\n", + "feat = shards_fc[\"features\"][0]\n", + "print(\"geometry type:\", feat[\"geometry\"][\"type\"])\n", + "print(\"properties:\", feat[\"properties\"])\n", + "[f[\"properties\"] for f in shards_fc[\"features\"]]" + ], + "id": "cell-6" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# One feature per granule footprint, decoded from the catalog.\n", + "[f[\"properties\"][\"id\"] for f in layers[\"granules\"][\"features\"]]" + ], + "id": "cell-7" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Grid-on-zoom gate\n", + "\n", + "`viewport_cells` draws shard-order cell outlines clipped to a viewport, but\n", + "**only** when `<= max_shards` shards intersect the viewport. Zoom in (a tight\n", + "bbox) and the grid appears; zoom out (a global bbox) and it returns an empty\n", + "collection rather than a globe-spanning graticule. The interactive map wires\n", + "this to the live `bounds`, so the grid blinks on as you zoom in." + ], + "id": "cell-8" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from zagg.viz import viewport_cells\n", + "from zagg.viz.shardmap import grid_from_signature\n", + "\n", + "# Tight viewport over a single shard -> gate open, cells drawn.\n", + "fp = grid_from_signature(shardmap.grid_signature).shard_footprint(0)\n", + "lon0, lat0, lon1, lat1 = fp.bounds\n", + "zoomed_in = viewport_cells(shardmap, (lon0, lat0, lon1, lat1), max_shards=4)\n", + "\n", + "# Global viewport -> too many shards visible -> gate closed, empty.\n", + "zoomed_out = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=1)\n", + "\n", + "print(\"zoomed-in cells:\", len(zoomed_in[\"features\"]))\n", + "print(\"zoomed-out cells (gated):\", len(zoomed_out[\"features\"]))" + ], + "id": "cell-9" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Antimeridian handling\n", + "\n", + "HEALPix shards near +-180 deg are split into a `MultiPolygon` so they don't\n", + "render as a band wrapping the whole globe. Here a HEALPix parent cell that\n", + "straddles the seam comes back as a `MultiPolygon`, each part confined to one\n", + "hemisphere." + ], + "id": "cell-10" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from zagg.viz import shard_outlines\n", + "\n", + "hp = HealpixGrid(2, 6, layout=\"fullsphere\")\n", + "key = int(\n", + " hp.coverage(\n", + " [(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))]\n", + " )[0]\n", + ")\n", + "sm_hp = ShardMap(hp.signature(), [key], [[{\"id\": \"G\", \"s3\": \"s\", \"https\": \"h\"}]], {})\n", + "geom = shard_outlines(sm_hp)[\"features\"][0][\"geometry\"]\n", + "geom[\"type\"]" + ], + "id": "cell-11" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Interactive map (ipyleaflet)\n", + "\n", + "`show_shardmap` builds an `ipyleaflet.Map` from a **saved** ShardMap (and an\n", + "optional catalog). We round-trip the synthetic objects to disk first, since the\n", + "viewer accepts either in-memory objects or file paths.\n", + "\n", + "**This cell needs `zagg[viz]` and is best run in JupyterLab** -- the widget\n", + "renders as a live, pannable/zoomable map there. Under headless `nbconvert` it\n", + "constructs the `Map` object (so the cell still executes) but won't show the\n", + "interactive tiles. See the verification checklist below for what to look for." + ], + "id": "cell-12" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "from pathlib import Path\n", + "\n", + "from zagg.viz import show_shardmap\n", + "\n", + "tmp = Path(tempfile.mkdtemp())\n", + "sm_path = tmp / \"shardmap.json\"\n", + "cat_path = tmp / \"catalog.parquet\"\n", + "shardmap.to_json(str(sm_path))\n", + "catalog.to_geoparquet(str(cat_path))\n", + "\n", + "# Construct the interactive map. In JupyterLab the returned Map renders inline.\n", + "m = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=12)\n", + "print(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\n", + "m" + ], + "id": "cell-13" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Manual in-browser verification\n", + "\n", + "Run this notebook in **JupyterLab** with the `viz` extra installed\n", + "(`pip install zagg[viz]`), execute the cells top-to-bottom, and confirm in the\n", + "map rendered by section 3:\n", + "\n", + "1. **Basemap** -- an OpenStreetMap context map renders and pans/zooms.\n", + "2. **Shard outlines** -- the blue `shards` polygons (one per shard key) sit\n", + " over the data area; hovering/clicking exposes their `shard_key` /\n", + " `n_granules` properties.\n", + "3. **Granule footprints toggle** -- open the layer control (top-right) and\n", + " toggle the red `granule footprints` layer on/off; the granule polygons\n", + " appear/disappear.\n", + "4. **Grid-on-zoom** -- when zoomed out, the `grid (shard cells)` layer is empty\n", + " (no global graticule). Zoom in until only a few shards fill the viewport and\n", + " the grey shard-cell outlines appear; zoom back out and they vanish. This is\n", + " the `max_shards` gate driven off the map's live `bounds`.\n", + "\n", + "To verify the **headless** core without a browser, run the test suite:\n", + "\n", + "```bash\n", + "uv run pytest tests/test_viz.py -v\n", + "```\n", + "\n", + "For a real dataset, replace the synthetic objects in section 1 with a shard map\n", + "and catalog built via `python -m zagg.catalog` (see the quickstart) and pass\n", + "their paths to `show_shardmap`." + ], + "id": "cell-14" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 9eabce811e361650db6d5f0b3485df89eee73a56 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 10:18:43 +0000 Subject: [PATCH 05/17] fold review: clarify saved-file round-trip in viewer notebook (#38) --- notebooks/shardmap_viewer.ipynb | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index 5f91503..dc01cf9 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -245,18 +245,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## 3. Interactive map (ipyleaflet)\n", - "\n", - "`show_shardmap` builds an `ipyleaflet.Map` from a **saved** ShardMap (and an\n", - "optional catalog). We round-trip the synthetic objects to disk first, since the\n", - "viewer accepts either in-memory objects or file paths.\n", - "\n", - "**This cell needs `zagg[viz]` and is best run in JupyterLab** -- the widget\n", - "renders as a live, pannable/zoomable map there. Under headless `nbconvert` it\n", - "constructs the `Map` object (so the cell still executes) but won't show the\n", - "interactive tiles. See the verification checklist below for what to look for." - ], + "source": "## 3. Interactive map (ipyleaflet)\n\n`show_shardmap` builds an `ipyleaflet.Map` from a ShardMap (and an optional\ncatalog), accepting either in-memory objects or file paths. Here we round-trip\nthe synthetic objects to disk to demonstrate the saved-file path you'd use in\npractice.\n\n**This cell needs `zagg[viz]` and is best run in JupyterLab** -- the widget\nrenders as a live, pannable/zoomable map there. Under headless `nbconvert` it\nconstructs the `Map` object (so the cell still executes) but won't show the\ninteractive tiles. See the verification checklist below for what to look for.", "id": "cell-12" }, { @@ -331,4 +320,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file From e9db108ef0adeb72852b5400a9684410868e7055 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 22:08:42 +0000 Subject: [PATCH 06/17] update notebook to real CMR data (issue #38) --- notebooks/shardmap_viewer.ipynb | 440 +++++++++++++++++++------------- 1 file changed, 264 insertions(+), 176 deletions(-) diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index dc01cf9..f99fee7 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -2,140 +2,220 @@ "cells": [ { "cell_type": "markdown", + "id": "cell-0", "metadata": {}, "source": [ "# Shard-map viewer\n", "\n", - "This notebook demonstrates and gives **verification instructions** for the\n", - "`zagg.viz` shard-map viewer (issue\n", - "[#38](https://github.com/englacial/zagg/issues/38)).\n", + "This notebook demonstrates the `zagg.viz` shard-map viewer (issue\n", + "[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\n", + "ATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\n", + "Login required — anonymous CMR queries are publicly accessible).\n", "\n", "The viewer has two layers:\n", "\n", - "- **Headless render core** (`zagg.viz.render_shardmap` and friends) -- pure\n", - " Python, no browser/widget stack. It turns a `ShardMap` (and an optional\n", - " `Catalog` / STAC-geoparquet file) into WGS84 GeoJSON `FeatureCollection`\n", - " dicts: shard outlines, granule footprints, and viewport-clipped grid cells.\n", - "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) -- builds an interactive\n", - " map from those collections: a basemap, the shard-outline layer, a toggleable\n", - " granule-footprint layer, and a grid layer that only draws when you zoom in\n", - " far enough (the \"grid-on-zoom\" gate -- never a global graticule).\n", + "- **Headless render core** (`zagg.viz.render_shardmap` and friends) — pure\n", + " Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n", + " `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines,\n", + " granule footprints, and viewport-clipped grid cells.\n", + "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) — builds an interactive\n", + " basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n", "\n", "## Install\n", "\n", - "The headless core needs only the zagg core deps (`shapely` + the grid\n", - "backends). The interactive map needs the optional `viz` extra:\n", - "\n", "```bash\n", - "pip install zagg[viz] # or: uv sync --extra viz\n", + "pip install zagg[viz] # core + ipyleaflet widget\n", + "pip install stac-geoparquet # catalog fetch dep (zagg catalog extra)\n", "```\n", "\n", - "Everything below uses a tiny **synthetic** ShardMap and catalog, so the\n", - "notebook runs with no network, S3, or Earthdata credentials." - ], - "id": "cell-0" + "**Requires an internet connection** for the anonymous CMR-STAC query in\n", + "section 1. No credentials are needed.\n", + "\n", + "## Areas of interest\n", + "\n", + "The notebook queries two AOIs:\n", + "\n", + "1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) — January\n", + " 2020. Dense ICESat-2 coverage near the pole; O(10–40) granules in two\n", + " weeks.\n", + "2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) — June\n", + " 2020. Second example showing the same pipeline on an Arctic AOI.\n", + "\n", + "Both inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\n", + "day one — see section 3 for the round-trip to disk and reload." + ] }, { "cell_type": "markdown", + "id": "cell-1", + "metadata": {}, + "source": [ + "## 1. Antarctic Peninsula — fetch real ATL06 granules from NASA CMR-STAC\n", + "\n", + "`CMRSource` speaks directly to NASA's CMR-STAC endpoint (`requests`).\n", + "No Earthdata Login or credentials are needed for anonymous granule-metadata\n", + "queries. The returned `Catalog` wraps a stac-geoparquet Arrow table (one row\n", + "per granule, both S3 and HTTPS asset hrefs preserved)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cell-2", "metadata": {}, + "outputs": [], "source": [ - "## 1. Build a small synthetic ShardMap and Catalog\n", + "from zagg.catalog.sources import CMRSource, Query\n", "\n", - "We hand-build a `ShardMap` over a 4x4-chunk rectilinear grid (UTM 18N, near\n", - "Washington DC) with three populated shards, plus a two-granule STAC `Catalog`.\n", - "This mirrors the test fixtures in `tests/test_viz.py`. In real use you would\n", - "instead load a saved shard map and catalog produced by `python -m zagg.catalog`\n", - "(see [quickstart](https://github.com/englacial/zagg/blob/main/docs/quickstart.md)):\n", + "# Antarctic Peninsula AOI — two weeks in January 2020.\n", + "AOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\n", + "START_DATE = \"2020-01-01\"\n", + "END_DATE = \"2020-01-15\"\n", "\n", - "```python\n", - "from zagg.catalog.shardmap import ShardMap\n", - "from zagg.catalog.sources import Catalog\n", - "sm = ShardMap.from_json(\"shardmap_ATL06_cycle22.json\")\n", - "cat = Catalog.from_geoparquet(\"catalog_ATL06_cycle22.parquet\")\n", - "```" - ], - "id": "cell-1" + "query = Query(\n", + " short_name=\"ATL06\",\n", + " version=\"006\",\n", + " start_date=START_DATE,\n", + " end_date=END_DATE,\n", + " region=AOI_BBOX,\n", + " provider=\"NSIDC_CPRD\",\n", + ")\n", + "\n", + "catalog = CMRSource().fetch(query)\n", + "print(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-3", "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "import pyarrow as pa\n", - "import stac_geoparquet.arrow as sga\n", + "# Inspect the catalog schema — each row is a STAC Item with WKB geometry.\n", + "print(\"Schema columns:\", catalog.table.schema.names[:10], \"...\")\n", + "print(\"Total rows:\", catalog.table.num_rows)\n", "\n", + "# Decode a few granule records to confirm footprints are present.\n", + "records = catalog.granule_records()\n", + "print(f\"\\n{len(records)} records with valid footprint geometry.\")\n", + "for rec in records[:3]:\n", + " print(f\" {rec['id']} | https: {str(rec['https'])[:70]}...\")" + ] + }, + { + "cell_type": "markdown", + "id": "cell-4", + "metadata": {}, + "source": [ + "## 2. Build a ShardMap on a HEALPix grid\n", + "\n", + "We use a HEALPix grid (`parent_order=6, child_order=12`) — the same\n", + "configuration as `src/zagg/configs/atl06.yaml`. The `mortie` backend is\n", + "used for HEALPix intersection and requires no extra install. If `spherely`\n", + "is installed, pass `backend='auto'` to prefer exact S2 intersections.\n", + "\n", + "`ShardMap.build` maps each grid shard (a parent-order HEALPix cell) to the\n", + "set of granules whose footprint intersects it. The manifest is self-contained\n", + "— it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", + "never needs the Catalog again at run time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cell-5", + "metadata": {}, + "outputs": [], + "source": [ "from zagg.catalog.shardmap import ShardMap\n", - "from zagg.catalog.sources import Catalog\n", - "from zagg.grids import HealpixGrid, RectilinearGrid\n", + "from zagg.grids import HealpixGrid\n", "\n", - "# A 10 km x 10 km rectilinear grid at 10 m, chunked 250x250 -> one chunk == one\n", - "# shard. Three shards are populated with synthetic granule references.\n", - "grid = RectilinearGrid(\n", - " \"EPSG:32618\", 10, [359400, 4300740, 369400, 4310740], [250, 250]\n", - ")\n", - "shard_keys = [0, 1, 5]\n", - "granules = [\n", - " [{\"id\": \"Ga\", \"s3\": \"s3://b/a.h5\", \"https\": \"https://h/a.h5\"}],\n", - " [{\"id\": \"Gb\", \"s3\": \"s3://b/b.h5\", \"https\": \"https://h/b.h5\"}],\n", - " [\n", - " {\"id\": \"Gb\", \"s3\": \"s3://b/b.h5\", \"https\": \"https://h/b.h5\"},\n", - " {\"id\": \"Gc\", \"s3\": \"s3://b/c.h5\", \"https\": \"https://h/c.h5\"},\n", - " ],\n", - "]\n", - "shardmap = ShardMap(grid.signature(), shard_keys, granules, {\"backend\": \"demo\"})\n", - "shardmap.shard_keys" - ], - "id": "cell-2" + "# HEALPix grid matching atl06.yaml (parent_order=6, child_order=12, fullsphere).\n", + "grid = HealpixGrid(parent_order=6, child_order=12, layout=\"fullsphere\")\n", + "\n", + "# Build: catalog footprints are intersected with the HEALPix shard cells\n", + "# that cover the AOI. Use backend='auto' to prefer spherely when available.\n", + "shardmap = ShardMap.build(catalog, grid, backend=\"mortie\")\n", + "\n", + "print(f\"ShardMap: {len(shardmap.shard_keys)} shards, \"\n", + " f\"{shardmap.metadata['total_pairs']} granule-shard pairs\")\n", + "print(f\"Build time: {shardmap.metadata['build_wall_s']:.2f}s \"\n", + " f\"backend: {shardmap.metadata['backend']}\")\n", + "print(f\"Grid signature: {shardmap.grid_signature}\")" + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-6", "metadata": {}, "outputs": [], "source": [ - "# A tiny STAC catalog of two granule footprints over the same area.\n", - "def _stac_item(gid, lon0, lon1, lat0=38.85, lat1=38.93):\n", - " ring = [[lon0, lat0], [lon1, lat0], [lon1, lat1], [lon0, lat1], [lon0, lat0]]\n", - " return {\n", - " \"type\": \"Feature\", \"stac_version\": \"1.0.0\", \"id\": gid,\n", - " \"geometry\": {\"type\": \"Polygon\", \"coordinates\": [ring]},\n", - " \"bbox\": [lon0, lat0, lon1, lat1],\n", - " \"properties\": {\"datetime\": \"2025-06-01T00:00:00Z\"},\n", - " \"collection\": \"DEMO\", \"stac_extensions\": [], \"links\": [],\n", - " \"assets\": {\n", - " \"data\": {\"href\": f\"https://h/{gid}.h5\", \"roles\": [\"data\"]},\n", - " \"data_s3\": {\"href\": f\"s3://b/{gid}.h5\", \"roles\": [\"data\"]},\n", - " },\n", - " }\n", - "\n", - "\n", - "items = [_stac_item(\"Ga\", -76.62, -76.57), _stac_item(\"Gb\", -76.55, -76.50)]\n", - "catalog = Catalog(\n", - " pa.table(sga.parse_stac_items_to_arrow(items)), {\"collection\": \"DEMO\"}\n", - ")\n", - "len(list(catalog.granule_records()))" - ], - "id": "cell-3" + "# Inspect a few shard assignments.\n", + "for key, granules in zip(shardmap.shard_keys[:4], shardmap.granules[:4]):\n", + " ids = [g[\"id\"] for g in granules]\n", + " print(f\" shard {key:6d}: {len(ids)} granule(s) — {ids[:2]}{'...' if len(ids) > 2 else ''}\")" + ] }, { "cell_type": "markdown", + "id": "cell-7", "metadata": {}, "source": [ - "## 2. Headless render core: GeoJSON FeatureCollections\n", + "## 3. Persist to disk and reload (round-trip)\n", + "\n", + "Both the ShardMap JSON and the STAC-geoparquet Catalog are supported as\n", + "file-path inputs to `show_shardmap`. Here we round-trip both to disk to\n", + "demonstrate the saved-file path you would use in practice (e.g. after\n", + "running `python -m zagg.catalog --config atl06.yaml ...`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cell-8", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "from pathlib import Path\n", + "\n", + "tmp = Path(tempfile.mkdtemp())\n", + "sm_path = tmp / \"shardmap_ATL06_peninsula_jan2020.json\"\n", + "cat_path = tmp / \"catalog_ATL06_peninsula_jan2020.parquet\"\n", + "\n", + "shardmap.to_json(str(sm_path))\n", + "catalog.to_geoparquet(str(cat_path))\n", + "\n", + "print(f\"ShardMap -> {sm_path} ({sm_path.stat().st_size / 1024:.1f} KB)\")\n", + "print(f\"Catalog -> {cat_path} ({cat_path.stat().st_size / 1024:.1f} KB)\")\n", + "\n", + "# Reload to verify round-trip.\n", + "from zagg.catalog.sources import Catalog\n", + "\n", + "sm_rt = ShardMap.from_json(str(sm_path))\n", + "cat_rt = Catalog.from_geoparquet(str(cat_path))\n", + "print(f\"\\nRound-trip OK: {len(sm_rt.shard_keys)} shards, {len(cat_rt)} granules\")" + ] + }, + { + "cell_type": "markdown", + "id": "cell-9", + "metadata": {}, + "source": [ + "## 4. Headless render core: GeoJSON FeatureCollections\n", "\n", "`render_shardmap` assembles every layer into one dict of GeoJSON\n", - "`FeatureCollection`s. Pass a `catalog` to include granule footprints and a\n", - "`bbox` to include the viewport-clipped grid cells. Each value is plain,\n", - "JSON-serializable GeoJSON -- no widgets required, so this is the part covered\n", - "by the unit tests." - ], - "id": "cell-4" + "`FeatureCollection`s. Pass a `catalog` to include granule footprints and\n", + "a `bbox` to include viewport-clipped grid cells. Each value is plain,\n", + "JSON-serializable GeoJSON — no widgets required." + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-10", "metadata": {}, "outputs": [], "source": [ @@ -143,168 +223,176 @@ "\n", "layers = render_shardmap(shardmap, catalog)\n", "{name: (fc and len(fc[\"features\"])) for name, fc in layers.items()}" - ], - "id": "cell-5" + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-11", "metadata": {}, "outputs": [], "source": [ - "# One feature per shard, with the shard key and its granule count. The grid is\n", - "# rebuilt from the map's own `grid_signature` -- the map is self-describing.\n", + "# One feature per shard, with the shard key and its granule count.\n", "shards_fc = layers[\"shards\"]\n", "feat = shards_fc[\"features\"][0]\n", "print(\"geometry type:\", feat[\"geometry\"][\"type\"])\n", "print(\"properties:\", feat[\"properties\"])\n", - "[f[\"properties\"] for f in shards_fc[\"features\"]]" - ], - "id": "cell-6" + "\n", + "# All populated shards in the Antarctic Peninsula window.\n", + "[f[\"properties\"] for f in shards_fc[\"features\"]][:5]" + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-12", "metadata": {}, "outputs": [], "source": [ - "# One feature per granule footprint, decoded from the catalog.\n", - "[f[\"properties\"][\"id\"] for f in layers[\"granules\"][\"features\"]]" - ], - "id": "cell-7" + "# One feature per granule footprint, decoded from the real Catalog.\n", + "print(f\"{len(layers['granules']['features'])} granule footprint(s)\")\n", + "[f[\"properties\"][\"id\"] for f in layers[\"granules\"][\"features\"]][:5]" + ] }, { "cell_type": "markdown", + "id": "cell-13", "metadata": {}, "source": [ "### Grid-on-zoom gate\n", "\n", - "`viewport_cells` draws shard-order cell outlines clipped to a viewport, but\n", - "**only** when `<= max_shards` shards intersect the viewport. Zoom in (a tight\n", - "bbox) and the grid appears; zoom out (a global bbox) and it returns an empty\n", - "collection rather than a globe-spanning graticule. The interactive map wires\n", - "this to the live `bounds`, so the grid blinks on as you zoom in." - ], - "id": "cell-8" + "`viewport_cells` draws shard-order cell outlines clipped to a viewport but\n", + "**only** when `<= max_shards` shards intersect the viewport. Zoom in and\n", + "the grid appears; zoom out and it returns an empty collection rather than a\n", + "globe-spanning graticule. The interactive map wires this to the live map\n", + "`bounds`, so the grid blinks on as you zoom in." + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-14", "metadata": {}, "outputs": [], "source": [ "from zagg.viz import viewport_cells\n", "from zagg.viz.shardmap import grid_from_signature\n", "\n", - "# Tight viewport over a single shard -> gate open, cells drawn.\n", - "fp = grid_from_signature(shardmap.grid_signature).shard_footprint(0)\n", + "# Tight viewport over one shard -> gate open, cells drawn.\n", + "first_shard = shardmap.shard_keys[0]\n", + "fp = grid_from_signature(shardmap.grid_signature).shard_footprint(first_shard)\n", "lon0, lat0, lon1, lat1 = fp.bounds\n", - "zoomed_in = viewport_cells(shardmap, (lon0, lat0, lon1, lat1), max_shards=4)\n", "\n", - "# Global viewport -> too many shards visible -> gate closed, empty.\n", + "zoomed_in = viewport_cells(shardmap, (lon0, lat0, lon1, lat1), max_shards=4)\n", "zoomed_out = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=1)\n", "\n", "print(\"zoomed-in cells:\", len(zoomed_in[\"features\"]))\n", "print(\"zoomed-out cells (gated):\", len(zoomed_out[\"features\"]))" - ], - "id": "cell-9" + ] }, { "cell_type": "markdown", + "id": "cell-15", "metadata": {}, "source": [ - "### Antimeridian handling\n", - "\n", - "HEALPix shards near +-180 deg are split into a `MultiPolygon` so they don't\n", - "render as a band wrapping the whole globe. Here a HEALPix parent cell that\n", - "straddles the seam comes back as a `MultiPolygon`, each part confined to one\n", - "hemisphere." - ], - "id": "cell-10" + "## 5. Interactive map — Antarctic Peninsula (ipyleaflet)\n", + "\n", + "`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\n", + "Catalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\n", + "as file paths or in-memory objects.\n", + "\n", + "**Run in JupyterLab** with `pip install zagg[viz]` to see the live map.\n", + "Under headless `nbconvert` the Map object is constructed (no error) but\n", + "tiles won't display.\n", + "\n", + "### Verification checklist\n", + "\n", + "1. **Basemap** — OpenStreetMap context map renders, pans, zooms.\n", + "2. **Shard outlines** — blue polygons over the Peninsula; click one for\n", + " its `shard_key` and `n_granules` properties.\n", + "3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n", + " track footprints appear/disappear.\n", + "4. **Grid-on-zoom** — zoom in until few shards fill the viewport to see\n", + " grey shard-cell outlines; zoom back out and they vanish." + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-16", "metadata": {}, "outputs": [], "source": [ - "from zagg.viz import shard_outlines\n", + "from zagg.viz import show_shardmap\n", "\n", - "hp = HealpixGrid(2, 6, layout=\"fullsphere\")\n", - "key = int(\n", - " hp.coverage(\n", - " [(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))]\n", - " )[0]\n", - ")\n", - "sm_hp = ShardMap(hp.signature(), [key], [[{\"id\": \"G\", \"s3\": \"s\", \"https\": \"h\"}]], {})\n", - "geom = shard_outlines(sm_hp)[\"features\"][0][\"geometry\"]\n", - "geom[\"type\"]" - ], - "id": "cell-11" + "# File-path interface — same as after `python -m zagg.catalog ...`.\n", + "m = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\n", + "print(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\n", + "m" + ] }, { "cell_type": "markdown", + "id": "cell-17", "metadata": {}, - "source": "## 3. Interactive map (ipyleaflet)\n\n`show_shardmap` builds an `ipyleaflet.Map` from a ShardMap (and an optional\ncatalog), accepting either in-memory objects or file paths. Here we round-trip\nthe synthetic objects to disk to demonstrate the saved-file path you'd use in\npractice.\n\n**This cell needs `zagg[viz]` and is best run in JupyterLab** -- the widget\nrenders as a live, pannable/zoomable map there. Under headless `nbconvert` it\nconstructs the `Map` object (so the cell still executes) but won't show the\ninteractive tiles. See the verification checklist below for what to look for.", - "id": "cell-12" + "source": [ + "## 6. Second AOI — Jakobshavn Glacier, West Greenland\n", + "\n", + "A second example using an Arctic AOI to show the same pipeline working\n", + "independently on a different region. Jakobshavn Glacier is one of\n", + "Greenland's fastest-moving glaciers and a standard ICESat-2 study region." + ] }, { "cell_type": "code", "execution_count": null, + "id": "cell-18", "metadata": {}, "outputs": [], "source": [ - "import tempfile\n", - "from pathlib import Path\n", + "# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\n", + "AOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\n", + "START2, END2 = \"2020-06-01\", \"2020-06-15\"\n", "\n", - "from zagg.viz import show_shardmap\n", + "query2 = Query(\n", + " short_name=\"ATL06\",\n", + " version=\"006\",\n", + " start_date=START2,\n", + " end_date=END2,\n", + " region=AOI2_BBOX,\n", + " provider=\"NSIDC_CPRD\",\n", + ")\n", "\n", - "tmp = Path(tempfile.mkdtemp())\n", - "sm_path = tmp / \"shardmap.json\"\n", - "cat_path = tmp / \"catalog.parquet\"\n", - "shardmap.to_json(str(sm_path))\n", - "catalog.to_geoparquet(str(cat_path))\n", + "catalog2 = CMRSource().fetch(query2)\n", + "shardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n", "\n", - "# Construct the interactive map. In JupyterLab the returned Map renders inline.\n", - "m = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=12)\n", - "print(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\n", - "m" - ], - "id": "cell-13" + "print(f\"Greenland AOI: {len(catalog2)} granules, \"\n", + " f\"{len(shardmap2.shard_keys)} shards\")\n", + "\n", + "# Headless render — confirm layers are populated.\n", + "layers2 = render_shardmap(shardmap2, catalog2)\n", + "print(\"Layer feature counts:\",\n", + " {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" + ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, + "id": "cell-19", "metadata": {}, + "outputs": [], "source": [ - "## 4. Manual in-browser verification\n", - "\n", - "Run this notebook in **JupyterLab** with the `viz` extra installed\n", - "(`pip install zagg[viz]`), execute the cells top-to-bottom, and confirm in the\n", - "map rendered by section 3:\n", - "\n", - "1. **Basemap** -- an OpenStreetMap context map renders and pans/zooms.\n", - "2. **Shard outlines** -- the blue `shards` polygons (one per shard key) sit\n", - " over the data area; hovering/clicking exposes their `shard_key` /\n", - " `n_granules` properties.\n", - "3. **Granule footprints toggle** -- open the layer control (top-right) and\n", - " toggle the red `granule footprints` layer on/off; the granule polygons\n", - " appear/disappear.\n", - "4. **Grid-on-zoom** -- when zoomed out, the `grid (shard cells)` layer is empty\n", - " (no global graticule). Zoom in until only a few shards fill the viewport and\n", - " the grey shard-cell outlines appear; zoom back out and they vanish. This is\n", - " the `max_shards` gate driven off the map's live `bounds`.\n", - "\n", - "To verify the **headless** core without a browser, run the test suite:\n", + "# Save and display interactive map for the Greenland AOI.\n", + "tmp2 = Path(tempfile.mkdtemp())\n", + "sm2_path = tmp2 / \"shardmap_ATL06_greenland_jun2020.json\"\n", + "cat2_path = tmp2 / \"catalog_ATL06_greenland_jun2020.parquet\"\n", "\n", - "```bash\n", - "uv run pytest tests/test_viz.py -v\n", - "```\n", + "shardmap2.to_json(str(sm2_path))\n", + "catalog2.to_geoparquet(str(cat2_path))\n", "\n", - "For a real dataset, replace the synthetic objects in section 1 with a shard map\n", - "and catalog built via `python -m zagg.catalog` (see the quickstart) and pass\n", - "their paths to `show_shardmap`." - ], - "id": "cell-14" + "m2 = show_shardmap(str(sm2_path), catalog=str(cat2_path), zoom=5)\n", + "m2" + ] } ], "metadata": { @@ -315,9 +403,9 @@ }, "language_info": { "name": "python", - "version": "3.12" + "version": "3.12.0" } }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} From b57f7c67ce3d205fcd3bb5042ce3cc596be82709 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 22:12:22 +0000 Subject: [PATCH 07/17] fold review: guard None https, fix loop var shadow, add re-run imports (#38) --- notebooks/shardmap_viewer.ipynb | 59 +++++++++++++++++---------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index f99fee7..9692d0e 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -10,22 +10,21 @@ "This notebook demonstrates the `zagg.viz` shard-map viewer (issue\n", "[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\n", "ATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\n", - "Login required — anonymous CMR queries are publicly accessible).\n", + "Login required \u2014 anonymous CMR queries are publicly accessible).\n", "\n", "The viewer has two layers:\n", "\n", - "- **Headless render core** (`zagg.viz.render_shardmap` and friends) — pure\n", + "- **Headless render core** (`zagg.viz.render_shardmap` and friends) \u2014 pure\n", " Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n", " `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines,\n", " granule footprints, and viewport-clipped grid cells.\n", - "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) — builds an interactive\n", + "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) \u2014 builds an interactive\n", " basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n", "\n", "## Install\n", "\n", "```bash\n", - "pip install zagg[viz] # core + ipyleaflet widget\n", - "pip install stac-geoparquet # catalog fetch dep (zagg catalog extra)\n", + "pip install \"zagg[viz,catalog]\" # core + ipyleaflet widget + stac-geoparquet\n", "```\n", "\n", "**Requires an internet connection** for the anonymous CMR-STAC query in\n", @@ -35,14 +34,14 @@ "\n", "The notebook queries two AOIs:\n", "\n", - "1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) — January\n", - " 2020. Dense ICESat-2 coverage near the pole; O(10–40) granules in two\n", + "1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) \u2014 January\n", + " 2020. Dense ICESat-2 coverage near the pole; O(10\u201340) granules in two\n", " weeks.\n", - "2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) — June\n", + "2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) \u2014 June\n", " 2020. Second example showing the same pipeline on an Arctic AOI.\n", "\n", "Both inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\n", - "day one — see section 3 for the round-trip to disk and reload." + "day one \u2014 see section 3 for the round-trip to disk and reload." ] }, { @@ -50,7 +49,7 @@ "id": "cell-1", "metadata": {}, "source": [ - "## 1. Antarctic Peninsula — fetch real ATL06 granules from NASA CMR-STAC\n", + "## 1. Antarctic Peninsula \u2014 fetch real ATL06 granules from NASA CMR-STAC\n", "\n", "`CMRSource` speaks directly to NASA's CMR-STAC endpoint (`requests`).\n", "No Earthdata Login or credentials are needed for anonymous granule-metadata\n", @@ -67,7 +66,7 @@ "source": [ "from zagg.catalog.sources import CMRSource, Query\n", "\n", - "# Antarctic Peninsula AOI — two weeks in January 2020.\n", + "# Antarctic Peninsula AOI \u2014 two weeks in January 2020.\n", "AOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\n", "START_DATE = \"2020-01-01\"\n", "END_DATE = \"2020-01-15\"\n", @@ -92,7 +91,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Inspect the catalog schema — each row is a STAC Item with WKB geometry.\n", + "# Inspect the catalog schema \u2014 each row is a STAC Item with WKB geometry.\n", "print(\"Schema columns:\", catalog.table.schema.names[:10], \"...\")\n", "print(\"Total rows:\", catalog.table.num_rows)\n", "\n", @@ -100,7 +99,7 @@ "records = catalog.granule_records()\n", "print(f\"\\n{len(records)} records with valid footprint geometry.\")\n", "for rec in records[:3]:\n", - " print(f\" {rec['id']} | https: {str(rec['https'])[:70]}...\")" + " print(f\" {rec['id']} | https: {(rec['https'] or 'None')[:70]}...\")" ] }, { @@ -110,14 +109,14 @@ "source": [ "## 2. Build a ShardMap on a HEALPix grid\n", "\n", - "We use a HEALPix grid (`parent_order=6, child_order=12`) — the same\n", + "We use a HEALPix grid (`parent_order=6, child_order=12`) \u2014 the same\n", "configuration as `src/zagg/configs/atl06.yaml`. The `mortie` backend is\n", "used for HEALPix intersection and requires no extra install. If `spherely`\n", "is installed, pass `backend='auto'` to prefer exact S2 intersections.\n", "\n", "`ShardMap.build` maps each grid shard (a parent-order HEALPix cell) to the\n", "set of granules whose footprint intersects it. The manifest is self-contained\n", - "— it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", + "\u2014 it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", "never needs the Catalog again at run time." ] }, @@ -153,9 +152,9 @@ "outputs": [], "source": [ "# Inspect a few shard assignments.\n", - "for key, granules in zip(shardmap.shard_keys[:4], shardmap.granules[:4]):\n", - " ids = [g[\"id\"] for g in granules]\n", - " print(f\" shard {key:6d}: {len(ids)} granule(s) — {ids[:2]}{'...' if len(ids) > 2 else ''}\")" + "for key, shard_granules in zip(shardmap.shard_keys[:4], shardmap.granules[:4]):\n", + " ids = [g[\"id\"] for g in shard_granules]\n", + " print(f\" shard {key:6d}: {len(ids)} granule(s) \u2014 {ids[:2]}{'...' if len(ids) > 2 else ''}\")" ] }, { @@ -209,7 +208,7 @@ "`render_shardmap` assembles every layer into one dict of GeoJSON\n", "`FeatureCollection`s. Pass a `catalog` to include granule footprints and\n", "a `bbox` to include viewport-clipped grid cells. Each value is plain,\n", - "JSON-serializable GeoJSON — no widgets required." + "JSON-serializable GeoJSON \u2014 no widgets required." ] }, { @@ -295,24 +294,24 @@ "id": "cell-15", "metadata": {}, "source": [ - "## 5. Interactive map — Antarctic Peninsula (ipyleaflet)\n", + "## 5. Interactive map \u2014 Antarctic Peninsula (ipyleaflet)\n", "\n", "`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\n", "Catalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\n", "as file paths or in-memory objects.\n", "\n", - "**Run in JupyterLab** with `pip install zagg[viz]` to see the live map.\n", + "**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\n", "Under headless `nbconvert` the Map object is constructed (no error) but\n", "tiles won't display.\n", "\n", "### Verification checklist\n", "\n", - "1. **Basemap** — OpenStreetMap context map renders, pans, zooms.\n", - "2. **Shard outlines** — blue polygons over the Peninsula; click one for\n", + "1. **Basemap** \u2014 OpenStreetMap context map renders, pans, zooms.\n", + "2. **Shard outlines** \u2014 blue polygons over the Peninsula; click one for\n", " its `shard_key` and `n_granules` properties.\n", - "3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n", + "3. **Granule footprints toggle** \u2014 layer-control (top-right); the ICESat-2\n", " track footprints appear/disappear.\n", - "4. **Grid-on-zoom** — zoom in until few shards fill the viewport to see\n", + "4. **Grid-on-zoom** \u2014 zoom in until few shards fill the viewport to see\n", " grey shard-cell outlines; zoom back out and they vanish." ] }, @@ -325,7 +324,7 @@ "source": [ "from zagg.viz import show_shardmap\n", "\n", - "# File-path interface — same as after `python -m zagg.catalog ...`.\n", + "# File-path interface \u2014 same as after `python -m zagg.catalog ...`.\n", "m = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\n", "print(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\n", "m" @@ -336,7 +335,7 @@ "id": "cell-17", "metadata": {}, "source": [ - "## 6. Second AOI — Jakobshavn Glacier, West Greenland\n", + "## 6. Second AOI \u2014 Jakobshavn Glacier, West Greenland\n", "\n", "A second example using an Arctic AOI to show the same pipeline working\n", "independently on a different region. Jakobshavn Glacier is one of\n", @@ -350,7 +349,9 @@ "metadata": {}, "outputs": [], "source": [ - "# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\n", + "# Imports also available from earlier cells; repeated here for standalone re-runs.\n", + "from zagg.viz import render_shardmap, show_shardmap\n", + "# Second AOI: Jakobshavn Glacier, West Greenland \u2014 June 2020.\n", "AOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\n", "START2, END2 = \"2020-06-01\", \"2020-06-15\"\n", "\n", @@ -369,7 +370,7 @@ "print(f\"Greenland AOI: {len(catalog2)} granules, \"\n", " f\"{len(shardmap2.shard_keys)} shards\")\n", "\n", - "# Headless render — confirm layers are populated.\n", + "# Headless render \u2014 confirm layers are populated.\n", "layers2 = render_shardmap(shardmap2, catalog2)\n", "print(\"Layer feature counts:\",\n", " {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" From 987c8e858e412ae7ab9a0a7324088c2cb9ee32c9 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 22 Jun 2026 19:24:48 +0000 Subject: [PATCH 08/17] phase C of issue #38 --- notebooks/shardmap_viewer.ipynb | 163 +++++------------------------- src/zagg/viz/crs.py | 174 ++++++++++++++++++++++++++++++++ src/zagg/viz/leaflet.py | 78 +++++++++++--- src/zagg/viz/shardmap.py | 54 ++++++++-- tests/test_viz.py | 130 ++++++++++++++++++++++++ 5 files changed, 443 insertions(+), 156 deletions(-) create mode 100644 src/zagg/viz/crs.py diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index 9692d0e..df2e72c 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -4,52 +4,14 @@ "cell_type": "markdown", "id": "cell-0", "metadata": {}, - "source": [ - "# Shard-map viewer\n", - "\n", - "This notebook demonstrates the `zagg.viz` shard-map viewer (issue\n", - "[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\n", - "ATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\n", - "Login required \u2014 anonymous CMR queries are publicly accessible).\n", - "\n", - "The viewer has two layers:\n", - "\n", - "- **Headless render core** (`zagg.viz.render_shardmap` and friends) \u2014 pure\n", - " Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n", - " `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines,\n", - " granule footprints, and viewport-clipped grid cells.\n", - "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) \u2014 builds an interactive\n", - " basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n", - "\n", - "## Install\n", - "\n", - "```bash\n", - "pip install \"zagg[viz,catalog]\" # core + ipyleaflet widget + stac-geoparquet\n", - "```\n", - "\n", - "**Requires an internet connection** for the anonymous CMR-STAC query in\n", - "section 1. No credentials are needed.\n", - "\n", - "## Areas of interest\n", - "\n", - "The notebook queries two AOIs:\n", - "\n", - "1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) \u2014 January\n", - " 2020. Dense ICESat-2 coverage near the pole; O(10\u201340) granules in two\n", - " weeks.\n", - "2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) \u2014 June\n", - " 2020. Second example showing the same pipeline on an Arctic AOI.\n", - "\n", - "Both inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\n", - "day one \u2014 see section 3 for the round-trip to disk and reload." - ] + "source": "# Shard-map viewer\n\nThis notebook demonstrates the `zagg.viz` shard-map viewer (issue\n[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\nATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\nLogin required — anonymous CMR queries are publicly accessible).\n\nThe viewer has two layers:\n\n- **Headless render core** (`zagg.viz.render_shardmap` and friends) — pure\n Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines,\n granule footprints, and viewport-clipped grid cells.\n- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) — builds an interactive\n basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n\n**Polar-aware projection.** `show_shardmap` picks the display CRS from the\nmap's extent: a polar AOI renders in NASA polar-stereographic (EPSG:3031\nAntarctic / EPSG:3413 Arctic) with a **GIBS** basemap — undistorted at the\npole — while mid-latitude AOIs stay on Web Mercator + OpenStreetMap.\n\n## Install\n\n```bash\npip install \"zagg[viz,catalog]\" # core + ipyleaflet widget + stac-geoparquet\n```\n\n**Requires an internet connection** for the anonymous CMR-STAC query in\nsection 1. No credentials are needed.\n\n## Areas of interest\n\nThe notebook queries two AOIs:\n\n1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) — January\n 2020. Dense ICESat-2 coverage near the pole; O(10–40) granules in two\n weeks. Renders in EPSG:3031 (Antarctic Polar Stereographic).\n2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) — June\n 2020. Second example showing the same pipeline on an Arctic AOI\n (EPSG:3413).\n\nBoth inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\nday one — see section 3 for the round-trip to disk and reload." }, { "cell_type": "markdown", "id": "cell-1", "metadata": {}, "source": [ - "## 1. Antarctic Peninsula \u2014 fetch real ATL06 granules from NASA CMR-STAC\n", + "## 1. Antarctic Peninsula — fetch real ATL06 granules from NASA CMR-STAC\n", "\n", "`CMRSource` speaks directly to NASA's CMR-STAC endpoint (`requests`).\n", "No Earthdata Login or credentials are needed for anonymous granule-metadata\n", @@ -63,26 +25,7 @@ "id": "cell-2", "metadata": {}, "outputs": [], - "source": [ - "from zagg.catalog.sources import CMRSource, Query\n", - "\n", - "# Antarctic Peninsula AOI \u2014 two weeks in January 2020.\n", - "AOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\n", - "START_DATE = \"2020-01-01\"\n", - "END_DATE = \"2020-01-15\"\n", - "\n", - "query = Query(\n", - " short_name=\"ATL06\",\n", - " version=\"006\",\n", - " start_date=START_DATE,\n", - " end_date=END_DATE,\n", - " region=AOI_BBOX,\n", - " provider=\"NSIDC_CPRD\",\n", - ")\n", - "\n", - "catalog = CMRSource().fetch(query)\n", - "print(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" - ] + "source": "from zagg.catalog.sources import CMRSource, Query\n\n# Antarctic Peninsula AOI — two weeks in January 2020.\nAOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\nSTART_DATE = \"2020-01-01\"\nEND_DATE = \"2020-01-15\"\n\nquery = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START_DATE,\n end_date=END_DATE,\n region=AOI_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog = CMRSource().fetch(query)\nprint(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" }, { "cell_type": "code", @@ -91,7 +34,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Inspect the catalog schema \u2014 each row is a STAC Item with WKB geometry.\n", + "# Inspect the catalog schema — each row is a STAC Item with WKB geometry.\n", "print(\"Schema columns:\", catalog.table.schema.names[:10], \"...\")\n", "print(\"Total rows:\", catalog.table.num_rows)\n", "\n", @@ -109,14 +52,14 @@ "source": [ "## 2. Build a ShardMap on a HEALPix grid\n", "\n", - "We use a HEALPix grid (`parent_order=6, child_order=12`) \u2014 the same\n", + "We use a HEALPix grid (`parent_order=6, child_order=12`) — the same\n", "configuration as `src/zagg/configs/atl06.yaml`. The `mortie` backend is\n", "used for HEALPix intersection and requires no extra install. If `spherely`\n", "is installed, pass `backend='auto'` to prefer exact S2 intersections.\n", "\n", "`ShardMap.build` maps each grid shard (a parent-order HEALPix cell) to the\n", "set of granules whose footprint intersects it. The manifest is self-contained\n", - "\u2014 it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", + "— it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", "never needs the Catalog again at run time." ] }, @@ -154,7 +97,7 @@ "# Inspect a few shard assignments.\n", "for key, shard_granules in zip(shardmap.shard_keys[:4], shardmap.granules[:4]):\n", " ids = [g[\"id\"] for g in shard_granules]\n", - " print(f\" shard {key:6d}: {len(ids)} granule(s) \u2014 {ids[:2]}{'...' if len(ids) > 2 else ''}\")" + " print(f\" shard {key:6d}: {len(ids)} granule(s) — {ids[:2]}{'...' if len(ids) > 2 else ''}\")" ] }, { @@ -208,7 +151,7 @@ "`render_shardmap` assembles every layer into one dict of GeoJSON\n", "`FeatureCollection`s. Pass a `catalog` to include granule footprints and\n", "a `bbox` to include viewport-clipped grid cells. Each value is plain,\n", - "JSON-serializable GeoJSON \u2014 no widgets required." + "JSON-serializable GeoJSON — no widgets required." ] }, { @@ -293,27 +236,7 @@ "cell_type": "markdown", "id": "cell-15", "metadata": {}, - "source": [ - "## 5. Interactive map \u2014 Antarctic Peninsula (ipyleaflet)\n", - "\n", - "`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\n", - "Catalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\n", - "as file paths or in-memory objects.\n", - "\n", - "**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\n", - "Under headless `nbconvert` the Map object is constructed (no error) but\n", - "tiles won't display.\n", - "\n", - "### Verification checklist\n", - "\n", - "1. **Basemap** \u2014 OpenStreetMap context map renders, pans, zooms.\n", - "2. **Shard outlines** \u2014 blue polygons over the Peninsula; click one for\n", - " its `shard_key` and `n_granules` properties.\n", - "3. **Granule footprints toggle** \u2014 layer-control (top-right); the ICESat-2\n", - " track footprints appear/disappear.\n", - "4. **Grid-on-zoom** \u2014 zoom in until few shards fill the viewport to see\n", - " grey shard-cell outlines; zoom back out and they vanish." - ] + "source": "## 5. Interactive map — Antarctic Peninsula (ipyleaflet, polar-stereographic)\n\n`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\nCatalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\nas file paths or in-memory objects.\n\n**Projection is auto-selected from the map's extent.** This Antarctic AOI is\nentirely south of −60°, so the viewer renders in **EPSG:3031** (Antarctic\nPolar Stereographic) with a NASA **GIBS** polar basemap instead of the\ndistorted Web Mercator default. An Arctic AOI (poleward of +60°) gets\n**EPSG:3413**; mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Pass\n`crs=\"EPSG:3857\"` to force Mercator, or `crs=\"EPSG:3413\"`/`\"EPSG:3031\"` to\nforce a specific pole. Vector layers stay WGS84 GeoJSON — proj4leaflet\nreprojects them client-side.\n\n**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\nUnder headless `nbconvert` the Map object is constructed (no error) but\ntiles won't display.\n\n### Verification checklist\n\n1. **Basemap** — the GIBS polar basemap renders, pans, zooms, and the\n continent is undistorted at the pole (no Web-Mercator stretching).\n2. **Shard outlines** — blue polygons over the Peninsula; click one for\n its `shard_key` and `n_granules` properties.\n3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n track footprints appear/disappear.\n4. **Grid-on-zoom** — see the next cell, which centers the map on a single\n shard so the grey shard-cell grid actually appears." }, { "cell_type": "code", @@ -321,21 +244,28 @@ "id": "cell-16", "metadata": {}, "outputs": [], - "source": [ - "from zagg.viz import show_shardmap\n", - "\n", - "# File-path interface \u2014 same as after `python -m zagg.catalog ...`.\n", - "m = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\n", - "print(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\n", - "m" - ] + "source": "from zagg.viz import show_shardmap\n\n# File-path interface — same as after `python -m zagg.catalog ...`.\n# CRS auto-selected from the AOI extent: this Antarctic map renders in EPSG:3031.\nm = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\nprint(\"display CRS:\", m.crs[\"name\"])\nprint(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\nm" + }, + { + "cell_type": "markdown", + "id": "a326b760", + "source": "### Zoom to one shard so the grid appears\n\nThe grid layer is gated on `<= max_shards` (default 4) shards intersecting the\nviewport — over the whole AOI far more than 4 shards are visible, so the gate\nstays shut and you only see the shard outlines. Centering the map on a single\nshard's bbox and zooming in puts `<= 4` shards in view, which fires the\n`Map.bounds` observer and makes `viewport_cells` emit the grey shard-cell grid.\nRun this cell, then watch the **grid (shard cells)** layer turn on.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "f497f2c3", + "source": "from zagg.viz.shardmap import grid_from_signature\n\n# Center on a single shard's bbox so <= max_shards fill the viewport.\none_shard = shardmap.shard_keys[0]\nslon0, slat0, slon1, slat1 = grid_from_signature(shardmap.grid_signature).shard_footprint(\n one_shard\n).bounds\nshard_center = ((slat0 + slat1) / 2, (slon0 + slon1) / 2) # (lat, lon)\n\nm.center = shard_center\nm.zoom = 9 # tight enough that <= 4 shards intersect -> grid gate opens\n\n# The grid layer is driven off m.bounds; setting center/zoom re-fires the\n# observer. In a live JupyterLab session the grey \"grid (shard cells)\" layer\n# now appears. (Under nbconvert m.bounds stays empty until the browser sizes\n# the map, so the grid only renders interactively.)\nprint(f\"centered on shard {one_shard} at {shard_center}, zoom={m.zoom}\")\nm", + "metadata": {}, + "execution_count": null, + "outputs": [] }, { "cell_type": "markdown", "id": "cell-17", "metadata": {}, "source": [ - "## 6. Second AOI \u2014 Jakobshavn Glacier, West Greenland\n", + "## 6. Second AOI — Jakobshavn Glacier, West Greenland\n", "\n", "A second example using an Arctic AOI to show the same pipeline working\n", "independently on a different region. Jakobshavn Glacier is one of\n", @@ -348,33 +278,7 @@ "id": "cell-18", "metadata": {}, "outputs": [], - "source": [ - "# Imports also available from earlier cells; repeated here for standalone re-runs.\n", - "from zagg.viz import render_shardmap, show_shardmap\n", - "# Second AOI: Jakobshavn Glacier, West Greenland \u2014 June 2020.\n", - "AOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\n", - "START2, END2 = \"2020-06-01\", \"2020-06-15\"\n", - "\n", - "query2 = Query(\n", - " short_name=\"ATL06\",\n", - " version=\"006\",\n", - " start_date=START2,\n", - " end_date=END2,\n", - " region=AOI2_BBOX,\n", - " provider=\"NSIDC_CPRD\",\n", - ")\n", - "\n", - "catalog2 = CMRSource().fetch(query2)\n", - "shardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n", - "\n", - "print(f\"Greenland AOI: {len(catalog2)} granules, \"\n", - " f\"{len(shardmap2.shard_keys)} shards\")\n", - "\n", - "# Headless render \u2014 confirm layers are populated.\n", - "layers2 = render_shardmap(shardmap2, catalog2)\n", - "print(\"Layer feature counts:\",\n", - " {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" - ] + "source": "# Imports also available from earlier cells; repeated here for standalone re-runs.\nfrom zagg.viz import render_shardmap, show_shardmap\n# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render — confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" }, { "cell_type": "code", @@ -382,18 +286,7 @@ "id": "cell-19", "metadata": {}, "outputs": [], - "source": [ - "# Save and display interactive map for the Greenland AOI.\n", - "tmp2 = Path(tempfile.mkdtemp())\n", - "sm2_path = tmp2 / \"shardmap_ATL06_greenland_jun2020.json\"\n", - "cat2_path = tmp2 / \"catalog_ATL06_greenland_jun2020.parquet\"\n", - "\n", - "shardmap2.to_json(str(sm2_path))\n", - "catalog2.to_geoparquet(str(cat2_path))\n", - "\n", - "m2 = show_shardmap(str(sm2_path), catalog=str(cat2_path), zoom=5)\n", - "m2" - ] + "source": "# Save and display interactive map for the Greenland AOI.\ntmp2 = Path(tempfile.mkdtemp())\nsm2_path = tmp2 / \"shardmap_ATL06_greenland_jun2020.json\"\ncat2_path = tmp2 / \"catalog_ATL06_greenland_jun2020.parquet\"\n\nshardmap2.to_json(str(sm2_path))\ncatalog2.to_geoparquet(str(cat2_path))\n\n# Arctic AOI -> auto-selects EPSG:3413 (NSIDC Sea Ice Polar Stereographic North).\nm2 = show_shardmap(str(sm2_path), catalog=str(cat2_path), zoom=5)\nprint(\"display CRS:\", m2.crs[\"name\"])\nm2" } ], "metadata": { @@ -409,4 +302,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/src/zagg/viz/crs.py b/src/zagg/viz/crs.py new file mode 100644 index 0000000..b043ce0 --- /dev/null +++ b/src/zagg/viz/crs.py @@ -0,0 +1,174 @@ +"""CRS selection for the shard-map viewer (issue #38, phase C). + +The default ipyleaflet basemap is EPSG:3857 (Web Mercator), which is unusable +above ~85 deg and badly distorts polar regions -- exactly the cryosphere AOIs +zagg exists for. This module picks a projection from a shard map's geographic +extent: NASA polar-stereographic for high-latitude AOIs, Web Mercator otherwise. + +- ``EPSG:3031`` -- Antarctic Polar Stereographic, when the map bbox is entirely + poleward of ~-60 deg latitude. +- ``EPSG:3413`` -- NSIDC Sea Ice Polar Stereographic North (Arctic), when the + bbox is entirely poleward of ~+60 deg. +- ``EPSG:3857`` -- Web Mercator, the mid-latitude / global fallback. + +For the two polar cases the matching proj4leaflet projection definition (proj4 +string, origin, resolutions, bounds) and a NASA **GIBS** WMTS basemap URL are +provided so :mod:`zagg.viz.leaflet` can wire them onto an ipyleaflet ``Map``. +The vector layers stay WGS84 GeoJSON -- proj4leaflet reprojects them +client-side -- so the headless render core is unchanged. + +Everything here is pure Python (no ipyleaflet): the bbox helper and the picker +are unit-testable with just the core deps. +""" +from __future__ import annotations + +from zagg.viz.shardmap import grid_from_signature + +# Latitude (deg) past which an AOI is treated as polar. A bbox whose nearer +# edge to the pole is beyond this cutoff selects a polar-stereographic CRS. +_POLAR_LAT = 60.0 + +# proj4leaflet projection definitions for the two NASA polar-stereographic +# grids. ``resolutions`` and ``bounds`` follow the GIBS tile matrix sets so the +# WMTS basemaps below line up (the published EPSG:3413 / EPSG:3031 GIBS tile +# pyramids, 8192 m down to 256 m per pixel over a +-4194304 m extent). +_PROJ_3413 = { + "name": "EPSG:3413", + "proj4def": ( + "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 " + "+a=6378137 +b=6356752.3142 +units=m +no_defs" + ), + "origin": [-4194304, 4194304], + "bounds": [[-4194304, -4194304], [4194304, 4194304]], + "resolutions": [8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0], +} + +_PROJ_3031 = { + "name": "EPSG:3031", + "proj4def": ( + "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 " + "+a=6378137 +b=6356752.3142 +units=m +no_defs" + ), + "origin": [-4194304, 4194304], + "bounds": [[-4194304, -4194304], [4194304, 4194304]], + "resolutions": [8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0], +} + +# NASA GIBS WMTS basemap (BlueMarble shaded relief, a static all-time layer) per +# polar tile matrix set. ``{z}/{y}/{x}`` is filled in by the leaflet TileLayer. +_GIBS_BASE = ( + "https://gibs.earthdata.nasa.gov/wmts/epsg{epsg}/best/" + "BlueMarble_ShadedRelief_Bathymetry/default/{matrix}/{{z}}/{{y}}/{{x}}.jpeg" +) + +_GIBS_3413 = { + "url": _GIBS_BASE.format(epsg="3413", matrix="500m"), + "attribution": "NASA EOSDIS GIBS", + "name": "GIBS BlueMarble (Arctic)", +} +_GIBS_3031 = { + "url": _GIBS_BASE.format(epsg="3031", matrix="500m"), + "attribution": "NASA EOSDIS GIBS", + "name": "GIBS BlueMarble (Antarctic)", +} + +# What each selectable CRS carries: the proj4leaflet definition (None for Web +# Mercator, which ipyleaflet ships natively) and the GIBS basemap (None -> the +# caller's Mercator default). +_CRS_INFO = { + "EPSG:3413": {"projection": _PROJ_3413, "basemap": _GIBS_3413}, + "EPSG:3031": {"projection": _PROJ_3031, "basemap": _GIBS_3031}, + "EPSG:3857": {"projection": None, "basemap": None}, +} + + +def shardmap_bbox(shardmap) -> tuple[float, float, float, float]: + """WGS84 ``(lon_min, lat_min, lon_max, lat_max)`` of a shard map's footprints. + + Computed from the union of the shard footprints (rebuilt from the map's own + ``grid_signature``), so it reflects exactly the area the viewer will draw. + + Raises + ------ + ValueError + If the map has no shards. + """ + grid = grid_from_signature(shardmap.grid_signature) + keys = list(shardmap.shard_keys) + if not keys: + raise ValueError("shard map has no shards") + lon_min = lat_min = float("inf") + lon_max = lat_max = float("-inf") + for key in keys: + x0, y0, x1, y1 = grid.shard_footprint(key).bounds + lon_min, lat_min = min(lon_min, x0), min(lat_min, y0) + lon_max, lat_max = max(lon_max, x1), max(lat_max, y1) + return (lon_min, lat_min, lon_max, lat_max) + + +def pick_crs(shardmap, override=None) -> str: + """Select an EPSG CRS for displaying ``shardmap``. + + Returns ``EPSG:3031`` when the map bbox is entirely poleward of ~-60 deg, + ``EPSG:3413`` when entirely poleward of ~+60 deg, else ``EPSG:3857`` (Web + Mercator). An explicit ``override`` (one of those three) is returned as-is. + + Parameters + ---------- + shardmap : ShardMap + A built or loaded shard map. + override : str, optional + Force a CRS. Must be one of ``"EPSG:3031"``, ``"EPSG:3413"``, + ``"EPSG:3857"``. + + Returns + ------- + str + One of the three EPSG codes above. + + Raises + ------ + ValueError + If ``override`` is given but not a supported CRS. + """ + if override is not None: + if override not in _CRS_INFO: + raise ValueError( + f"unsupported crs override {override!r}; " + f"expected one of {sorted(_CRS_INFO)}" + ) + return override + + _lon_min, lat_min, _lon_max, lat_max = shardmap_bbox(shardmap) + if lat_max <= -_POLAR_LAT: + return "EPSG:3031" + if lat_min >= _POLAR_LAT: + return "EPSG:3413" + return "EPSG:3857" + + +def crs_info(crs: str) -> dict: + """proj4leaflet ``projection`` and GIBS ``basemap`` dicts for an EPSG ``crs``. + + Returns ``{"projection": dict | None, "basemap": dict | None}``. ``None`` + values mean "use ipyleaflet's native Web Mercator / the caller's default + basemap" -- i.e. the EPSG:3857 case. + + Raises + ------ + ValueError + If ``crs`` is not a supported code. + """ + if crs not in _CRS_INFO: + raise ValueError( + f"unsupported crs {crs!r}; expected one of {sorted(_CRS_INFO)}" + ) + return _CRS_INFO[crs] + + +def is_polar(crs: str) -> bool: + """True for a polar-stereographic CRS (no +-180 antimeridian seam).""" + return crs in ("EPSG:3413", "EPSG:3031") + + +__all__ = ["shardmap_bbox", "pick_crs", "crs_info", "is_polar"] diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index e85c275..8faa510 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -1,4 +1,4 @@ -"""ipyleaflet wrapper for the shard-map viewer (issue #38, phase 2). +"""ipyleaflet wrapper for the shard-map viewer (issue #38, phases 2 & C). Builds an interactive map from a saved ShardMap: a basemap, the shard-outline layer, an optional (toggleable) granule-footprint layer, and a grid layer that @@ -6,16 +6,20 @@ to show ``<= max_shards`` shards (the "grid-on-zoom" gate -- never a global graticule, issue #38). +The CRS is chosen from the map's extent (:mod:`zagg.viz.crs`): polar AOIs get a +NASA polar-stereographic projection (EPSG:3413/3031) with a matching **GIBS** +WMTS basemap, mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Vector +layers stay WGS84 GeoJSON -- proj4leaflet reprojects them client-side -- so the +headless render core is unchanged; the only seam difference is that the +-180 +antimeridian split is skipped under a polar CRS (there is no such seam there). + All ``ipyleaflet`` imports are local to the functions here so importing :mod:`zagg.viz` (and the phase-1 render core / test suite) never requires the widget stack. Install it with ``pip install zagg[viz]``. - -The geometry is produced entirely by the headless core in -:mod:`zagg.viz.shardmap` -- this module only wires those GeoJSON collections -onto ipyleaflet layers and keeps the grid layer in sync with the viewport. """ from __future__ import annotations +from zagg.viz.crs import crs_info, is_polar, pick_crs from zagg.viz.shardmap import ( _load_catalog, _load_shardmap, @@ -51,6 +55,24 @@ def _walk_coords(coords, lons, lats): _walk_coords(sub, lons, lats) +def _leaflet_crs(projection: dict): + """A proj4leaflet ``ipyleaflet.projections`` CRS dict from a projection def. + + ipyleaflet's ``Map.crs`` accepts a dict matching proj4leaflet's + ``L.Proj.CRS`` signature (``name`` + ``proj4def`` + an ``options`` block with + ``origin``/``bounds``/``resolutions``). :mod:`zagg.viz.crs` carries those + values per polar EPSG so they line up with the GIBS tile matrix set. + """ + return { + "name": projection["name"], + "custom": True, + "proj4def": projection["proj4def"], + "origin": projection["origin"], + "bounds": projection["bounds"], + "resolutions": projection["resolutions"], + } + + def show_shardmap( shardmap_path, catalog=None, @@ -58,9 +80,16 @@ def show_shardmap( max_shards: int = 4, zoom: int = 3, basemap=None, + crs=None, ): """Build an interactive ipyleaflet map for a saved ShardMap. + The display CRS is auto-selected from the map's extent: a polar AOI gets a + NASA polar-stereographic projection (EPSG:3413 Arctic / EPSG:3031 Antarctic) + with a matching GIBS WMTS basemap; mid-latitude AOIs keep Web Mercator + + OpenStreetMap. Pass ``crs=`` to force one of ``"EPSG:3031"``, + ``"EPSG:3413"``, ``"EPSG:3857"``. + Parameters ---------- shardmap_path : str or ShardMap @@ -73,7 +102,9 @@ def show_shardmap( zoom : int Initial map zoom. basemap : ipyleaflet basemap, optional - Overrides the default OpenStreetMap basemap. + Overrides the default basemap (OSM for Web Mercator, GIBS for polar). + crs : str, optional + Force the display CRS instead of auto-picking from the map extent. Returns ------- @@ -82,19 +113,39 @@ def show_shardmap( grid layer, and a ``LayersControl`` for toggling layers. """ # Import first so a missing `viz` extra fails clearly, before any work. - from ipyleaflet import GeoJSON, LayersControl, Map, basemaps + from ipyleaflet import GeoJSON, LayersControl, Map, TileLayer, basemaps shardmap = _load_shardmap(shardmap_path) - shards_fc = shard_outlines(shardmap) - center = _center_zoom(shards_fc) - m = Map(basemap=basemap or basemaps.OpenStreetMap.Mapnik, center=center, zoom=zoom) + selected_crs = pick_crs(shardmap, override=crs) + polar = is_polar(selected_crs) + info = crs_info(selected_crs) + # Under a polar CRS the +-180 seam does not exist, so skip the split. + split_seam = not polar + + shards_fc = shard_outlines(shardmap, split_seam=split_seam) + + map_kwargs = {"center": _center_zoom(shards_fc), "zoom": zoom} + if info["projection"] is not None: + map_kwargs["crs"] = _leaflet_crs(info["projection"]) + + if basemap is not None: + map_kwargs["basemap"] = basemap + elif info["basemap"] is not None: + gibs = info["basemap"] + map_kwargs["basemap"] = TileLayer( + url=gibs["url"], attribution=gibs["attribution"], name=gibs["name"] + ) + else: + map_kwargs["basemap"] = basemaps.OpenStreetMap.Mapnik + + m = Map(**map_kwargs) shard_layer = GeoJSON(data=shards_fc, style=_SHARD_STYLE, name="shards") m.add(shard_layer) if catalog is not None: - footprint_fc = granule_footprints(_load_catalog(catalog)) + footprint_fc = granule_footprints(_load_catalog(catalog), split_seam=split_seam) footprint_layer = GeoJSON( data=footprint_fc, style=_FOOTPRINT_STYLE, name="granule footprints" ) @@ -116,7 +167,10 @@ def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) return (south, west), (north, east) = bounds grid_layer.data = viewport_cells( - shardmap, (west, south, east, north), max_shards=max_shards + shardmap, + (west, south, east, north), + max_shards=max_shards, + split_seam=split_seam, ) m.observe(_refresh_grid, names="bounds") diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py index ec18c52..1695e45 100644 --- a/src/zagg/viz/shardmap.py +++ b/src/zagg/viz/shardmap.py @@ -157,12 +157,33 @@ def _polygon_rings(geom, *, shift: float = 0.0) -> list[list[list[float]]]: return [[[lon + shift, lat] for lon, lat in _ring_list(r)] for r in rings] -def _polygon_geometry(geom) -> dict: - """GeoJSON geometry for a shapely Polygon/MultiPolygon, antimeridian-safe. +def _plain_geometry(geom) -> dict: + """GeoJSON geometry for a Polygon/MultiPolygon with no antimeridian split. - Coordinates are plain ``list`` (not shapely's tuples) so the result is + Used under a polar-stereographic CRS, where the +-180 seam is a Mercator/ + WGS84 artifact that does not exist in the projected plane (the singularity is + the opposite pole instead). Coordinates are plain ``list`` so the result is canonical, ``json``-stable GeoJSON. """ + if geom.geom_type == "MultiPolygon": + return { + "type": "MultiPolygon", + "coordinates": [_polygon_rings(sub) for sub in geom.geoms], + } + if geom.geom_type == "Polygon": + return {"type": "Polygon", "coordinates": _polygon_rings(geom)} + raise ValueError(f"unsupported geometry type for GeoJSON: {geom.geom_type}") + + +def _polygon_geometry(geom, *, split_seam: bool = True) -> dict: + """GeoJSON geometry for a shapely Polygon/MultiPolygon, antimeridian-safe. + + When ``split_seam`` is False (polar CRS), the +-180 split is skipped -- see + :func:`_plain_geometry`. Coordinates are plain ``list`` (not shapely's + tuples) so the result is canonical, ``json``-stable GeoJSON. + """ + if not split_seam: + return _plain_geometry(geom) if geom.geom_type == "MultiPolygon": coords: list = [] for sub in geom.geoms: @@ -210,7 +231,7 @@ def _collection(features: list[dict]) -> dict: # ── layers ─────────────────────────────────────────────────────────────────── -def shard_outlines(shardmap) -> dict: +def shard_outlines(shardmap, *, split_seam: bool = True) -> dict: """Shard/chunk outlines as a GeoJSON ``FeatureCollection``. One feature per shard key in the map, with the shard's footprint polygon and @@ -221,6 +242,9 @@ def shard_outlines(shardmap) -> dict: ---------- shardmap : ShardMap A built or loaded shard map. + split_seam : bool + Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a + polar-stereographic CRS, where there is no antimeridian seam. Returns ------- @@ -233,14 +257,14 @@ def shard_outlines(shardmap) -> dict: geom = grid.shard_footprint(key) features.append( _feature( - _polygon_geometry(geom), + _polygon_geometry(geom, split_seam=split_seam), {"shard_key": _jsonable(key), "n_granules": len(granules)}, ) ) return _collection(features) -def granule_footprints(catalog) -> dict: +def granule_footprints(catalog, *, split_seam: bool = True) -> dict: """Granule footprints as a GeoJSON ``FeatureCollection``. One polygon feature per granule in the catalog (its footprint @@ -250,6 +274,9 @@ def granule_footprints(catalog) -> dict: ---------- catalog : Catalog A loaded catalog (provides ``granule_records``). + split_seam : bool + Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a + polar-stereographic CRS, where there is no antimeridian seam. Returns ------- @@ -264,12 +291,15 @@ def granule_footprints(catalog) -> dict: if len(ring) < 4: continue features.append( - _feature(_polygon_geometry(Polygon(ring)), {"id": rec["id"]}) + _feature( + _polygon_geometry(Polygon(ring), split_seam=split_seam), + {"id": rec["id"]}, + ) ) return _collection(features) -def viewport_cells(shardmap, bbox, *, max_shards: int = 4) -> dict: +def viewport_cells(shardmap, bbox, *, max_shards: int = 4, split_seam: bool = True) -> dict: """Shard-order cell outlines clipped to a viewport, gated on visible shards. Implements the "grid-on-zoom" behavior (issue #38): cell outlines **at the @@ -285,6 +315,9 @@ def viewport_cells(shardmap, bbox, *, max_shards: int = 4) -> dict: Viewport ``(lon_min, lat_min, lon_max, lat_max)`` in WGS84. max_shards : int Maximum number of intersecting shards for the grid to render. + split_seam : bool + Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a + polar-stereographic CRS, where there is no antimeridian seam. Returns ------- @@ -312,7 +345,10 @@ def viewport_cells(shardmap, bbox, *, max_shards: int = 4) -> dict: if clipped is None: continue features.append( - _feature(_polygon_geometry(clipped), {"shard_key": _jsonable(key)}) + _feature( + _polygon_geometry(clipped, split_seam=split_seam), + {"shard_key": _jsonable(key)}, + ) ) return _collection(features) diff --git a/tests/test_viz.py b/tests/test_viz.py index 7f7531e..ac7f124 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -22,6 +22,7 @@ shard_outlines, viewport_cells, ) +from zagg.viz.crs import crs_info, is_polar, pick_crs, shardmap_bbox from zagg.viz.shardmap import _is_geojson, _split_antimeridian # ── fixtures ───────────────────────────────────────────────────────────────── @@ -72,6 +73,20 @@ def catalog(): ) +@pytest.fixture +def antarctic_shardmap(): + """ShardMap on an EPSG:3031 grid -> footprints entirely south of -60 deg.""" + g = RectilinearGrid("EPSG:3031", 100000, [-1000000, -1000000, 0, 0], [5, 5]) + return ShardMap(g.signature(), [0, 1], [[], []], {"backend": "test"}) + + +@pytest.fixture +def arctic_shardmap(): + """ShardMap on an EPSG:3413 grid -> footprints entirely north of +60 deg.""" + g = RectilinearGrid("EPSG:3413", 100000, [-1000000, -1000000, 0, 0], [5, 5]) + return ShardMap(g.signature(), [0, 1], [[], []], {"backend": "test"}) + + # ── grid_from_signature ────────────────────────────────────────────────────── class TestGridFromSignature: @@ -253,6 +268,88 @@ def test_from_geoparquet_path(self, shardmap, catalog, tmp_path): assert len(out["granules"]["features"]) == 2 +# ── phase C: CRS selection (headless, no browser) ──────────────────────────── + +class TestCrsSelection: + def test_midlatitude_is_web_mercator(self, shardmap): + # UTM 18N grid near 38.9N -> mid-latitude -> Web Mercator default. + assert pick_crs(shardmap) == "EPSG:3857" + + def test_antarctic_is_3031(self, antarctic_shardmap): + assert pick_crs(antarctic_shardmap) == "EPSG:3031" + + def test_arctic_is_3413(self, arctic_shardmap): + assert pick_crs(arctic_shardmap) == "EPSG:3413" + + def test_explicit_override_wins(self, shardmap): + # A mid-latitude map forced onto a polar CRS returns the override. + assert pick_crs(shardmap, override="EPSG:3031") == "EPSG:3031" + + def test_bad_override_raises(self, shardmap): + with pytest.raises(ValueError, match="unsupported crs override"): + pick_crs(shardmap, override="EPSG:9999") + + def test_bbox_from_footprints(self, antarctic_shardmap): + lon_min, lat_min, lon_max, lat_max = shardmap_bbox(antarctic_shardmap) + assert lat_max <= -60.0 + assert lon_min >= -180.0 and lon_max <= 180.0 + + def test_empty_map_bbox_raises(self, shardmap): + empty = ShardMap(shardmap.grid_signature, [], [], {}) + with pytest.raises(ValueError, match="no shards"): + shardmap_bbox(empty) + + +class TestCrsInfo: + def test_web_mercator_has_no_projection_or_basemap(self): + info = crs_info("EPSG:3857") + assert info["projection"] is None + assert info["basemap"] is None + assert not is_polar("EPSG:3857") + + @pytest.mark.parametrize("crs", ["EPSG:3031", "EPSG:3413"]) + def test_polar_has_projection_and_gibs_basemap(self, crs): + info = crs_info(crs) + assert is_polar(crs) + proj = info["projection"] + assert proj["name"] == crs + assert proj["proj4def"].startswith("+proj=stere") + assert proj["origin"] and proj["bounds"] and proj["resolutions"] + assert "gibs.earthdata.nasa.gov" in info["basemap"]["url"] + assert f"epsg{crs.split(':')[1]}" in info["basemap"]["url"] + + def test_bad_crs_raises(self): + with pytest.raises(ValueError, match="unsupported crs"): + crs_info("EPSG:9999") + + +# ── phase C: CRS-aware antimeridian seam ───────────────────────────────────── + +class TestSeamAwareLayers: + def test_polar_skips_antimeridian_split(self): + # A HEALPix cell straddling +-180 is a MultiPolygon under the Mercator + # split, but stays a single Polygon when split_seam=False (polar CRS). + grid = HealpixGrid(2, 6, layout="fullsphere") + key = int( + grid.coverage( + [(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))] + )[0] + ) + sm = ShardMap( + grid.signature(), [key], + [[{"id": "G", "s3": "s", "https": "h"}]], {}, + ) + split = shard_outlines(sm)["features"][0]["geometry"] + unsplit = shard_outlines(sm, split_seam=False)["features"][0]["geometry"] + assert split["type"] == "MultiPolygon" + assert unsplit["type"] == "Polygon" + + def test_granule_footprints_seam_flag(self, catalog): + fc = granule_footprints(catalog, split_seam=False) + assert _is_geojson(fc) + assert all(f["geometry"]["type"] == "Polygon" for f in fc["features"]) + + # ── phase 2: ipyleaflet wrapper (skips when the viz extra isn't installed) ──── class TestShowShardmap: @@ -289,3 +386,36 @@ def test_build_map_with_catalog(self, shardmap, catalog, tmp_path): assert isinstance(m, Map) names = {getattr(layer, "name", None) for layer in m.layers} assert "granule footprints" in names + + def test_polar_map_uses_3031_crs_and_gibs(self, antarctic_shardmap, tmp_path): + pytest.importorskip("ipyleaflet") + from zagg.viz import show_shardmap + + path = tmp_path / "sm.json" + antarctic_shardmap.to_json(str(path)) + m = show_shardmap(str(path)) + # proj4leaflet CRS for EPSG:3031 is wired onto the Map. + assert m.crs["name"] == "EPSG:3031" + # GIBS Antarctic tile basemap is present. + urls = [getattr(layer, "url", "") for layer in m.layers] + assert any("epsg3031" in u for u in urls) + + def test_midlatitude_map_stays_web_mercator(self, shardmap, tmp_path): + pytest.importorskip("ipyleaflet") + from zagg.viz import show_shardmap + + path = tmp_path / "sm.json" + shardmap.to_json(str(path)) + m = show_shardmap(str(path)) + # No custom proj4leaflet CRS -> ipyleaflet's default Web Mercator. + assert m.crs["name"] == "EPSG3857" + assert not m.crs.get("custom", False) + + def test_crs_override(self, shardmap, tmp_path): + pytest.importorskip("ipyleaflet") + from zagg.viz import show_shardmap + + path = tmp_path / "sm.json" + shardmap.to_json(str(path)) + m = show_shardmap(str(path), crs="EPSG:3413") + assert m.crs["name"] == "EPSG:3413" From e6f3e60adabf2b9f07066e0ad0b4d71dc8fa272a Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 22 Jun 2026 19:29:26 +0000 Subject: [PATCH 09/17] fold review: align polar CRS defs to GIBS, strengthen seam test (#38) --- src/zagg/viz/crs.py | 34 +++++++++++++++++++++++----------- src/zagg/viz/leaflet.py | 8 ++++---- tests/test_viz.py | 24 ++++++++++++++++++++---- 3 files changed, 47 insertions(+), 19 deletions(-) diff --git a/src/zagg/viz/crs.py b/src/zagg/viz/crs.py index b043ce0..f86ca24 100644 --- a/src/zagg/viz/crs.py +++ b/src/zagg/viz/crs.py @@ -29,29 +29,36 @@ _POLAR_LAT = 60.0 # proj4leaflet projection definitions for the two NASA polar-stereographic -# grids. ``resolutions`` and ``bounds`` follow the GIBS tile matrix sets so the -# WMTS basemaps below line up (the published EPSG:3413 / EPSG:3031 GIBS tile -# pyramids, 8192 m down to 256 m per pixel over a +-4194304 m extent). +# grids. These mirror ipyleaflet's own bundled ``projections.EPSG3413/3031 +# ["NASAGIBS"]`` (proj4 string, origin, bounds, and the full 7-level GIBS 500m +# tile-matrix-set resolution pyramid, 16384 m down to 256 m per pixel over a +# +-4194304 m extent) so the WMTS basemaps below line up zoom-for-zoom. Kept +# here as plain dicts so this module stays ipyleaflet-free and headlessly +# testable. +_GIBS_RESOLUTIONS = [16384.0, 8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0] +_GIBS_ORIGIN = [-4194304, 4194304] +_GIBS_BOUNDS = [[-4194304, -4194304], [4194304, 4194304]] + _PROJ_3413 = { "name": "EPSG:3413", "proj4def": ( "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 " - "+a=6378137 +b=6356752.3142 +units=m +no_defs" + "+ellps=WGS84 +datum=WGS84 +units=m +no_defs" ), - "origin": [-4194304, 4194304], - "bounds": [[-4194304, -4194304], [4194304, 4194304]], - "resolutions": [8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0], + "origin": _GIBS_ORIGIN, + "bounds": _GIBS_BOUNDS, + "resolutions": _GIBS_RESOLUTIONS, } _PROJ_3031 = { "name": "EPSG:3031", "proj4def": ( "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 " - "+a=6378137 +b=6356752.3142 +units=m +no_defs" + "+datum=WGS84 +units=m +no_defs" ), - "origin": [-4194304, 4194304], - "bounds": [[-4194304, -4194304], [4194304, 4194304]], - "resolutions": [8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0], + "origin": _GIBS_ORIGIN, + "bounds": _GIBS_BOUNDS, + "resolutions": _GIBS_RESOLUTIONS, } # NASA GIBS WMTS basemap (BlueMarble shaded relief, a static all-time layer) per @@ -88,6 +95,11 @@ def shardmap_bbox(shardmap) -> tuple[float, float, float, float]: Computed from the union of the shard footprints (rebuilt from the map's own ``grid_signature``), so it reflects exactly the area the viewer will draw. + Note: for a footprint that encloses a pole the longitude extent degenerates + to roughly the full ``[-180, 180]`` (the cell wraps the pole), so callers + should rely on the latitude extent for a pole-enclosing map. :func:`pick_crs` + keys on latitude only for exactly this reason. + Raises ------ ValueError diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index 8faa510..34235ab 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -58,10 +58,10 @@ def _walk_coords(coords, lons, lats): def _leaflet_crs(projection: dict): """A proj4leaflet ``ipyleaflet.projections`` CRS dict from a projection def. - ipyleaflet's ``Map.crs`` accepts a dict matching proj4leaflet's - ``L.Proj.CRS`` signature (``name`` + ``proj4def`` + an ``options`` block with - ``origin``/``bounds``/``resolutions``). :mod:`zagg.viz.crs` carries those - values per polar EPSG so they line up with the GIBS tile matrix set. + ipyleaflet's ``Map.crs`` accepts a flat dict (``name``, ``custom``, + ``proj4def``, ``origin``, ``bounds``, ``resolutions``) -- the same shape as + its bundled ``projections.EPSG3413["NASAGIBS"]``. :mod:`zagg.viz.crs` carries + those values per polar EPSG so they line up with the GIBS tile matrix set. """ return { "name": projection["name"], diff --git a/tests/test_viz.py b/tests/test_viz.py index ac7f124..9007557 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -344,10 +344,26 @@ def test_polar_skips_antimeridian_split(self): assert split["type"] == "MultiPolygon" assert unsplit["type"] == "Polygon" - def test_granule_footprints_seam_flag(self, catalog): - fc = granule_footprints(catalog, split_seam=False) - assert _is_geojson(fc) - assert all(f["geometry"]["type"] == "Polygon" for f in fc["features"]) + def test_granule_footprints_seam_flag(self): + # A footprint straddling +-180 (ring jumps 178 -> -178): split_seam=True + # cuts it into a MultiPolygon, but split_seam=False (polar CRS) keeps the + # single ring -- so the flag, not the data, decides the geometry type. + ring = [[178.0, 70.0], [-178.0, 70.0], [-178.0, 72.0], [178.0, 72.0], [178.0, 70.0]] + item = { + "type": "Feature", "stac_version": "1.0.0", "id": "Gx", + "geometry": {"type": "Polygon", "coordinates": [ring]}, + "bbox": [-180.0, 70.0, 180.0, 72.0], + "properties": {"datetime": "2025-06-01T00:00:00Z"}, + "collection": "TEST", "stac_extensions": [], "links": [], + "assets": {"data": {"href": "https://h/Gx.h5", "roles": ["data"]}}, + } + cat = Catalog( + pa.table(sga.parse_stac_items_to_arrow([item])), {"collection": "TEST"} + ) + split = granule_footprints(cat)["features"][0]["geometry"] + unsplit = granule_footprints(cat, split_seam=False)["features"][0]["geometry"] + assert split["type"] == "MultiPolygon" + assert unsplit["type"] == "Polygon" # ── phase 2: ipyleaflet wrapper (skips when the viz extra isn't installed) ──── From 09111597c3e518c1ecea08c1c93995cbfa4731c4 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 22 Jun 2026 20:13:43 +0000 Subject: [PATCH 10/17] fix grid-on-zoom hang via STRtree shard index + debounce (#38) --- src/zagg/viz/leaflet.py | 44 +++++++++++++++++++-- src/zagg/viz/shardmap.py | 85 ++++++++++++++++++++++++++++++++++++---- tests/test_viz.py | 75 +++++++++++++++++++++++++++++++++++ 3 files changed, 193 insertions(+), 11 deletions(-) diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index 34235ab..0b59205 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -19,11 +19,14 @@ """ from __future__ import annotations +import threading + from zagg.viz.crs import crs_info, is_polar, pick_crs from zagg.viz.shardmap import ( _load_catalog, _load_shardmap, granule_footprints, + shard_index, shard_outlines, viewport_cells, ) @@ -33,6 +36,35 @@ _FOOTPRINT_STYLE = {"color": "#e31a1c", "weight": 1, "fillOpacity": 0.10} _GRID_STYLE = {"color": "#333333", "weight": 1, "fillOpacity": 0.0} +# Coalesce a burst of pan/zoom ``bounds`` ticks into one grid refresh after the +# viewport settles (seconds). Without this the grid recomputes on every +# intermediate frame of a drag. +_GRID_DEBOUNCE_S = 0.2 + + +def _debounce(wait, func): + """Wrap ``func`` so rapid calls coalesce to one call ``wait`` s after the last. + + Each call cancels the pending :class:`threading.Timer` and reschedules, so a + burst of events fires ``func`` once when they stop. Returns the wrapper with a + ``cancel()`` to tear the timer down. + """ + timer: list = [None] + + def wrapper(*args, **kwargs): + if timer[0] is not None: + timer[0].cancel() + timer[0] = threading.Timer(wait, lambda: func(*args, **kwargs)) + timer[0].daemon = True + timer[0].start() + + def cancel(): + if timer[0] is not None: + timer[0].cancel() + + wrapper.cancel = cancel + return wrapper + def _center_zoom(fc: dict): """Center ``(lat, lon)`` for a FeatureCollection's bbox (zoom is left default).""" @@ -151,8 +183,9 @@ def show_shardmap( ) m.add(footprint_layer) - # Grid-on-zoom: an empty layer kept in sync with the viewport. Recomputed on - # every bounds change; the headless core's gate returns nothing when too + # Grid-on-zoom: an empty layer kept in sync with the viewport. The shard + # footprints are indexed once (STRtree) so each refresh is a cheap query, not + # a full footprint rebuild; the headless core's gate returns nothing when too # many shards are visible, so this never becomes a global graticule. grid_layer = GeoJSON( data={"type": "FeatureCollection", "features": []}, @@ -161,6 +194,8 @@ def show_shardmap( ) m.add(grid_layer) + index = shard_index(shardmap) + def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) bounds = m.bounds if not bounds: @@ -171,9 +206,12 @@ def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) (west, south, east, north), max_shards=max_shards, split_seam=split_seam, + index=index, ) - m.observe(_refresh_grid, names="bounds") + # Debounce so a burst of intermediate pan/zoom ticks coalesces into one + # refresh after movement settles. + m.observe(_debounce(_GRID_DEBOUNCE_S, _refresh_grid), names="bounds") _refresh_grid() m.add(LayersControl(position="topright")) diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py index 1695e45..8f1057d 100644 --- a/src/zagg/viz/shardmap.py +++ b/src/zagg/viz/shardmap.py @@ -299,7 +299,71 @@ def granule_footprints(catalog, *, split_seam: bool = True) -> dict: return _collection(features) -def viewport_cells(shardmap, bbox, *, max_shards: int = 4, split_seam: bool = True) -> dict: +class ShardIndex: + """One-time shard-footprint index for fast viewport queries (issue #38). + + Building a shard footprint is not free -- on HEALPix it is + ``mortie.tools.mort2polygon`` per key. The interactive viewer re-queries on + every pan/zoom, so regenerating all footprints per event freezes the kernel. + This builds the ``(key, footprint)`` list **once** and indexes the + footprints in a shapely ``STRtree``; :meth:`query` then answers a viewport + in ``O(log N + k)`` instead of ``O(N)`` footprint rebuilds. + """ + + def __init__(self, shardmap): + from shapely import STRtree + + grid = grid_from_signature(shardmap.grid_signature) + self.keys = list(shardmap.shard_keys) + self.footprints = [grid.shard_footprint(k) for k in self.keys] + self._tree = STRtree(self.footprints) + + def query(self, view): + """``[(key, footprint)]`` whose footprint truly intersects ``view``. + + The ``STRtree`` returns bbox-overlap candidates; each is refined with a + precise ``footprint.intersects(view)`` so the result matches the old + exhaustive scan exactly. + """ + out = [] + for i in self._tree.query(view): + fp = self.footprints[i] + if fp.intersects(view): + out.append((self.keys[i], fp)) + return out + + +# Cache one ShardIndex per shardmap object so repeated viewport queries (and the +# back-compat viewport_cells path) never rebuild footprints. Keyed by id() and +# guarded by a weak ref so a collected shardmap drops its entry. +_INDEX_CACHE: dict = {} + + +def shard_index(shardmap) -> ShardIndex: + """Return a cached :class:`ShardIndex` for ``shardmap`` (built once).""" + import weakref + + key = id(shardmap) + entry = _INDEX_CACHE.get(key) + if entry is not None: + return entry[0] + + index = ShardIndex(shardmap) + + def _drop(_ref, key=key): + _INDEX_CACHE.pop(key, None) + + try: + ref = weakref.ref(shardmap, _drop) + except TypeError: + ref = None # not weak-referenceable; entry simply lingers + _INDEX_CACHE[key] = (index, ref) + return index + + +def viewport_cells( + shardmap, bbox, *, max_shards: int = 4, split_seam: bool = True, index=None +) -> dict: """Shard-order cell outlines clipped to a viewport, gated on visible shards. Implements the "grid-on-zoom" behavior (issue #38): cell outlines **at the @@ -307,6 +371,10 @@ def viewport_cells(shardmap, bbox, *, max_shards: int = 4, split_seam: bool = Tr ``bbox``. When more shards are visible the viewport is too zoomed-out for a useful grid and an empty collection is returned -- never a global graticule. + Footprints are taken from a one-time :class:`ShardIndex` (built and cached + per shardmap, or passed in via ``index=``), so repeated calls -- e.g. the + interactive viewer's per-pan refresh -- never rebuild them. + Parameters ---------- shardmap : ShardMap @@ -318,6 +386,9 @@ def viewport_cells(shardmap, bbox, *, max_shards: int = 4, split_seam: bool = Tr split_seam : bool Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a polar-stereographic CRS, where there is no antimeridian seam. + index : ShardIndex, optional + A prebuilt index to query. When omitted, the cached index for + ``shardmap`` is used (:func:`shard_index`). Returns ------- @@ -327,15 +398,11 @@ def viewport_cells(shardmap, bbox, *, max_shards: int = 4, split_seam: bool = Tr """ from shapely.geometry import box - grid = grid_from_signature(shardmap.grid_signature) + if index is None: + index = shard_index(shardmap) view = box(bbox[0], bbox[1], bbox[2], bbox[3]) - # Build each footprint once (reproject + densify is not free), then gate. - visible = [ - (key, fp) - for key, fp in ((k, grid.shard_footprint(k)) for k in shardmap.shard_keys) - if fp.intersects(view) - ] + visible = index.query(view) if not visible or len(visible) > max_shards: return _collection([]) @@ -431,5 +498,7 @@ def _is_geojson(obj) -> bool: "shard_outlines", "granule_footprints", "viewport_cells", + "ShardIndex", + "shard_index", "render_shardmap", ] diff --git a/tests/test_viz.py b/tests/test_viz.py index 9007557..47d28fa 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -182,6 +182,61 @@ def test_clipped_to_viewport(self, shardmap): for lon, lat in ring: assert view[0] - 1e-6 <= lon <= view[2] + 1e-6 + def test_footprints_built_once_across_queries(self, shardmap, monkeypatch): + # Regression for the grid-on-zoom hang (PR #44): footprints must be + # generated once for the STRtree index and never again per query. Wrap + # the grid's shard_footprint with a call counter and run many viewports. + from zagg.grids.rectilinear import RectilinearGrid + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() # fresh index for this shardmap + calls = {"n": 0} + orig = RectilinearGrid.shard_footprint + + def counting(self, key): + calls["n"] += 1 + return orig(self, key) + + monkeypatch.setattr(RectilinearGrid, "shard_footprint", counting) + + full = (-180, -90, 180, 90) + for _ in range(10): + viewport_cells(shardmap, full, max_shards=10) + # Exactly one footprint build per shard for index construction; the ten + # queries add nothing. + assert calls["n"] == len(shardmap.shard_keys) + + def test_prebuilt_index_reused(self, shardmap, monkeypatch): + # Passing an explicit index never touches shard_footprint at query time. + from zagg.grids.rectilinear import RectilinearGrid + from zagg.viz.shardmap import ShardIndex + + index = ShardIndex(shardmap) + + def boom(self, key): # pragma: no cover - must not be called + raise AssertionError("shard_footprint rebuilt during query") + + monkeypatch.setattr(RectilinearGrid, "shard_footprint", boom) + fc = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=10, index=index) + assert _is_geojson(fc) + + def test_index_query_matches_exhaustive_scan(self, shardmap): + # The STRtree-backed visible set equals the old exhaustive intersects + # scan, for both a tight and a wide bbox. + from shapely.geometry import box + + from zagg.viz.shardmap import grid_from_signature, shard_index + + grid = grid_from_signature(shardmap.grid_signature) + idx = shard_index(shardmap) + for bbox in ((-180, -90, 180, 90), grid.shard_footprint(0).bounds): + view = box(*bbox) + indexed = {k for k, _ in idx.query(view)} + exhaustive = { + k for k in shardmap.shard_keys if grid.shard_footprint(k).intersects(view) + } + assert indexed == exhaustive + # ── antimeridian splitting ─────────────────────────────────────────────────── @@ -366,6 +421,26 @@ def test_granule_footprints_seam_flag(self): assert unsplit["type"] == "Polygon" +# ── debounce wrapper (pure stdlib; no widget stack) ────────────────────────── + +class TestDebounce: + def test_rapid_calls_coalesce_to_one(self): + # Importing leaflet pulls only stdlib + zagg.viz at module level (the + # ipyleaflet imports are local to the functions), so this is headless. + import time + + from zagg.viz.leaflet import _debounce + + calls = {"n": 0} + deb = _debounce(0.05, lambda: calls.__setitem__("n", calls["n"] + 1)) + for _ in range(20): + deb() # each cancels and reschedules the prior timer + assert calls["n"] == 0 # nothing fired yet (still within the window) + time.sleep(0.2) + assert calls["n"] == 1 # the burst coalesced into a single refresh + deb.cancel() + + # ── phase 2: ipyleaflet wrapper (skips when the viz extra isn't installed) ──── class TestShowShardmap: From 342e43a34df2b761448281fcc28c5ff69dff54ec Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 22 Jun 2026 20:17:18 +0000 Subject: [PATCH 11/17] fold review: expose debounce cancel hook (#38) --- src/zagg/viz/leaflet.py | 7 +++++-- tests/test_viz.py | 3 +++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index 0b59205..35825f3 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -210,8 +210,11 @@ def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) ) # Debounce so a burst of intermediate pan/zoom ticks coalesces into one - # refresh after movement settles. - m.observe(_debounce(_GRID_DEBOUNCE_S, _refresh_grid), names="bounds") + # refresh after movement settles. Keep a handle so the pending timer can be + # cancelled (exposed as ``m.cancel_grid_refresh``) for clean teardown. + debounced_refresh = _debounce(_GRID_DEBOUNCE_S, _refresh_grid) + m.observe(debounced_refresh, names="bounds") + m.cancel_grid_refresh = debounced_refresh.cancel _refresh_grid() m.add(LayersControl(position="topright")) diff --git a/tests/test_viz.py b/tests/test_viz.py index 47d28fa..613c1a9 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -462,6 +462,9 @@ def test_build_map(self, shardmap, tmp_path): assert isinstance(m, Map) # shard layer + grid layer (+ basemap) present. assert len(m.layers) >= 2 + # Debounced bounds observer is wired with a reachable cancel hook. + assert callable(getattr(m, "cancel_grid_refresh", None)) + m.cancel_grid_refresh() def test_build_map_with_catalog(self, shardmap, catalog, tmp_path): pytest.importorskip("ipyleaflet") From 93373edc8d0eff8dba52f671070bfbf49166b8f9 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 22 Jun 2026 21:32:37 +0000 Subject: [PATCH 12/17] render nested child-cell grid on zoom; main-thread debounce (#38) --- src/zagg/viz/leaflet.py | 48 +++++-- src/zagg/viz/shardmap.py | 129 ++++++++++++++++--- tests/test_viz.py | 263 ++++++++++++++++++++++++++++++++++----- 3 files changed, 384 insertions(+), 56 deletions(-) diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index 35825f3..90001d6 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -17,9 +17,10 @@ :mod:`zagg.viz` (and the phase-1 render core / test suite) never requires the widget stack. Install it with ``pip install zagg[viz]``. """ + from __future__ import annotations -import threading +import asyncio from zagg.viz.crs import crs_info, is_polar, pick_crs from zagg.viz.shardmap import ( @@ -45,22 +46,45 @@ def _debounce(wait, func): """Wrap ``func`` so rapid calls coalesce to one call ``wait`` s after the last. - Each call cancels the pending :class:`threading.Timer` and reschedules, so a - burst of events fires ``func`` once when they stop. Returns the wrapper with a - ``cancel()`` to tear the timer down. + The coalesced call is scheduled on the **kernel's own asyncio event loop** + (``loop.call_later``) -- the same main thread that owns the ipywidgets comm + -- not on a background ``threading.Timer``. That distinction is the fix for + the kernel-crash report (PR #44): the Jupyter widget comm channel is not + thread-safe, so mutating a widget traitlet (here ``grid_layer.data``) from a + timer thread can hang or corrupt the comm and crash the kernel. Scheduling on + the running loop keeps every refresh on the main thread. + + Each call cancels the pending :class:`asyncio.TimerHandle` and reschedules, + so a burst of events fires ``func`` once when they stop. Returns the wrapper + with a ``cancel()`` to tear the pending call down. + + When no event loop is running (e.g. a plain script / headless test) the call + is run synchronously -- there is no comm to protect and no loop to schedule + on, so coalescing is moot. """ - timer: list = [None] + handle: list = [None] + + def _fire(args, kwargs): + handle[0] = None + func(*args, **kwargs) def wrapper(*args, **kwargs): - if timer[0] is not None: - timer[0].cancel() - timer[0] = threading.Timer(wait, lambda: func(*args, **kwargs)) - timer[0].daemon = True - timer[0].start() + if handle[0] is not None: + handle[0].cancel() + handle[0] = None + try: + loop = asyncio.get_running_loop() + except RuntimeError: + loop = None + if loop is None: + func(*args, **kwargs) + return + handle[0] = loop.call_later(wait, _fire, args, kwargs) def cancel(): - if timer[0] is not None: - timer[0].cancel() + if handle[0] is not None: + handle[0].cancel() + handle[0] = None wrapper.cancel = cancel return wrapper diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py index 8f1057d..54c0031 100644 --- a/src/zagg/viz/shardmap.py +++ b/src/zagg/viz/shardmap.py @@ -13,9 +13,12 @@ needed. - :func:`granule_footprints` -- one polygon feature per granule footprint, decoded from a ``Catalog`` (``granule_records``). -- :func:`viewport_cells` -- shard-order cell outlines clipped to a viewport - bbox, emitted **only** when ``<= max_shards`` shards intersect the viewport - (the "grid-on-zoom" gate -- never a global graticule, issue #38). +- :func:`viewport_cells` -- the grid-on-zoom layer clipped to a viewport bbox, + emitted **only** when ``<= max_shards`` shards intersect the viewport (the + gate -- never a global graticule, issue #38). For HEALPix this is the finer + **child cells at ``child_order``** nesting inside the visible shards (generated + viewport-bounded via ``mortie.morton_coverage``, not by per-shard fan-out); + for other grids it is the shard footprint clipped to the viewport. Antimeridian handling --------------------- @@ -26,6 +29,7 @@ ones that genuinely *cross* it, :func:`_split_antimeridian` cuts the ring at +-180 into a ``MultiPolygon`` so GeoJSON consumers draw it correctly. """ + from __future__ import annotations import json @@ -88,6 +92,7 @@ def grid_from_signature(signature: dict): # ── GeoJSON geometry helpers ───────────────────────────────────────────────── + def _ring_list(ring) -> list[list[float]]: """A shapely ring's ``[[lon, lat], ...]`` as plain floats.""" x, y = ring.coords.xy @@ -214,9 +219,7 @@ def _polygonal(geom): polys = [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")] if not polys: return None - merged = MultiPolygon( - [p for g in polys for p in getattr(g, "geoms", [g])] - ) + merged = MultiPolygon([p for g in polys for p in getattr(g, "geoms", [g])]) return merged return None @@ -231,6 +234,7 @@ def _collection(features: list[dict]) -> dict: # ── layers ─────────────────────────────────────────────────────────────────── + def shard_outlines(shardmap, *, split_seam: bool = True) -> dict: """Shard/chunk outlines as a GeoJSON ``FeatureCollection``. @@ -361,15 +365,81 @@ def _drop(_ref, key=key): return index +# Edge samples per side for a child-cell polygon. Child cells are small +# near-squares, so a low step traces them faithfully while keeping each ring -- +# and the GeoJSON shipped over the widget comm -- light. (HealpixGrid.shard_ +# footprint uses step=32 for the much larger parent diamonds.) +_CHILD_CELL_STEP = 8 + + +def _child_cell_polygon(cell): + """Polygon (WGS84 lon/lat) for a single HEALPix child cell morton id.""" + from mortie.tools import mort2polygon + from shapely.geometry import Polygon + + verts = mort2polygon(int(cell), step=_CHILD_CELL_STEP) + lats = np.array([v[0] for v in verts]) + lons = np.array([v[1] for v in verts]) + return Polygon(zip(lons, lats)) + + +def _healpix_child_cells(grid, visible_keys, view, bbox, max_cells): + """Child cells at ``child_order`` nesting inside the visible HEALPix shards. + + The grid-on-zoom for HEALPix is the **finer child grid that subdivides each + shard** (issue #38) -- not the shard outline redrawn. Children are generated + in a **viewport-bounded** way: a single ``mortie.morton_coverage`` query over + the viewport bbox at ``child_order`` returns only the child cells overlapping + the view, so the count scales with the viewport, not with the + ``4^(child_order - parent_order)`` per-shard fan-out. The returned cells are + filtered to those whose parent (``clip2order(parent_order, ...)``) is one of + the visible shards, so the rendered grid lines up *inside* the shard outlines. + + Returns ``(cells, parents)`` morton-id arrays, or ``None`` when the in-view + child count exceeds ``max_cells`` (the analogous gate to ``max_shards`` -- + keeps a single dense-viewport refresh bounded rather than emitting tens of + thousands of polygons). + """ + from mortie import clip2order, morton_coverage + + lon_min, lat_min, lon_max, lat_max = bbox + view_lats = np.array([lat_min, lat_min, lat_max, lat_max, lat_min]) + view_lons = np.array([lon_min, lon_max, lon_max, lon_min, lon_min]) + cov = np.asarray(morton_coverage([view_lats], [view_lons], order=grid.child_order)) + if cov.size == 0: + return np.empty(0, dtype=np.int64), np.empty(0, dtype=np.int64) + parents = clip2order(grid.parent_order, cov) + keep = np.isin(parents, np.asarray(list(visible_keys))) + cells = cov[keep] + if cells.size > max_cells: + return None + return cells, parents[keep] + + def viewport_cells( - shardmap, bbox, *, max_shards: int = 4, split_seam: bool = True, index=None + shardmap, + bbox, + *, + max_shards: int = 4, + max_cells: int = 2000, + split_seam: bool = True, + index=None, ) -> dict: - """Shard-order cell outlines clipped to a viewport, gated on visible shards. - - Implements the "grid-on-zoom" behavior (issue #38): cell outlines **at the - shard order** are drawn only when ``<= max_shards`` shards intersect - ``bbox``. When more shards are visible the viewport is too zoomed-out for a - useful grid and an empty collection is returned -- never a global graticule. + """Grid-on-zoom cell outlines clipped to a viewport, gated on visible shards. + + Implements the "grid-on-zoom" behavior (issue #38): the finer grid is drawn + only when ``<= max_shards`` shards intersect ``bbox``. When more shards are + visible the viewport is too zoomed-out for a useful grid and an empty + collection is returned -- never a global graticule. + + For a **HEALPix** grid the emitted features are the **child cells at + ``child_order`` that nest inside the visible shards** -- the grid that + subdivides each shard 4-for-1 per order step -- not the shard outline + redrawn. They are generated viewport-bounded via a single + :func:`mortie.morton_coverage` query at ``child_order`` (so the cell count + scales with the viewport, not with the ``4^(child_order - parent_order)`` + per-shard fan-out), filtered to the visible shards, and clipped to ``bbox``. + For any other grid the shard footprint clipped to the viewport is emitted. Footprints are taken from a one-time :class:`ShardIndex` (built and cached per shardmap, or passed in via ``index=``), so repeated calls -- e.g. the @@ -383,6 +453,10 @@ def viewport_cells( Viewport ``(lon_min, lat_min, lon_max, lat_max)`` in WGS84. max_shards : int Maximum number of intersecting shards for the grid to render. + max_cells : int + Maximum number of in-view HEALPix child cells to render. Above this the + viewport is too dense (zoomed out relative to ``child_order``) and an + empty collection is returned, keeping a single refresh bounded. split_seam : bool Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a polar-stereographic CRS, where there is no antimeridian seam. @@ -393,8 +467,8 @@ def viewport_cells( Returns ------- dict - GeoJSON ``FeatureCollection`` of shard-cell outlines clipped to the - viewport, or an empty collection when the gate is not met. + GeoJSON ``FeatureCollection`` of grid outlines clipped to the viewport, + or an empty collection when a gate is not met. """ from shapely.geometry import box @@ -406,6 +480,29 @@ def viewport_cells( if not visible or len(visible) > max_shards: return _collection([]) + grid = grid_from_signature(shardmap.grid_signature) + + from zagg.grids import HealpixGrid + + if isinstance(grid, HealpixGrid): + visible_keys = {int(key) for key, _ in visible} + result = _healpix_child_cells(grid, visible_keys, view, bbox, max_cells) + if result is None: + return _collection([]) + cells, parents = result + features = [] + for cell, parent in zip(cells, parents): + clipped = _polygonal(_child_cell_polygon(cell).intersection(view)) + if clipped is None: + continue + features.append( + _feature( + _polygon_geometry(clipped, split_seam=split_seam), + {"cell": int(cell), "shard_key": int(parent)}, + ) + ) + return _collection(features) + features = [] for key, fp in visible: clipped = _polygonal(fp.intersection(view)) @@ -422,6 +519,7 @@ def viewport_cells( # ── top-level assembly ─────────────────────────────────────────────────────── + def render_shardmap(shardmap, catalog=None, *, bbox=None, max_shards: int = 4) -> dict: """Assemble all viewer layers for a shard map into one dict of collections. @@ -456,6 +554,7 @@ def render_shardmap(shardmap, catalog=None, *, bbox=None, max_shards: int = 4) - # ── small loaders / utils ──────────────────────────────────────────────────── + def _load_shardmap(shardmap): """Accept a ShardMap or a JSON path; return a ShardMap.""" if isinstance(shardmap, (str, Path)): diff --git a/tests/test_viz.py b/tests/test_viz.py index 613c1a9..86a4954 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -27,14 +27,19 @@ # ── fixtures ───────────────────────────────────────────────────────────────── + def _item(gid, lon0, lon1, lat0=38.85, lat1=38.93): ring = [[lon0, lat0], [lon1, lat0], [lon1, lat1], [lon0, lat1], [lon0, lat0]] return { - "type": "Feature", "stac_version": "1.0.0", "id": gid, + "type": "Feature", + "stac_version": "1.0.0", + "id": gid, "geometry": {"type": "Polygon", "coordinates": [ring]}, "bbox": [lon0, lat0, lon1, lat1], "properties": {"datetime": "2025-06-01T00:00:00Z"}, - "collection": "TEST", "stac_extensions": [], "links": [], + "collection": "TEST", + "stac_extensions": [], + "links": [], "assets": { "data": {"href": f"https://h/{gid}.h5", "roles": ["data"]}, "data_s3": {"href": f"s3://b/{gid}.h5", "roles": ["data"]}, @@ -44,9 +49,7 @@ def _item(gid, lon0, lon1, lat0=38.85, lat1=38.93): @pytest.fixture def rect_grid(): - return RectilinearGrid( - "EPSG:32618", 10, [359400, 4300740, 369400, 4310740], [250, 250] - ) + return RectilinearGrid("EPSG:32618", 10, [359400, 4300740, 369400, 4310740], [250, 250]) @pytest.fixture @@ -73,6 +76,22 @@ def catalog(): ) +@pytest.fixture +def healpix_shardmap(): + """A HEALPix shardmap (parent_order=6, child_order=12) over a tight AOI. + + Child cells subdivide each order-6 shard 4-for-1 per order step + (``4^(12-6)`` per shard), so the grid-on-zoom must show those nested child + cells -- not the shard outline redrawn. + """ + g = HealpixGrid(parent_order=6, child_order=12, layout="fullsphere") + lats = np.array([10.0, 10.5, 10.5, 10.0, 10.0]) + lons = np.array([20.0, 20.0, 20.5, 20.5, 20.0]) + keys = [int(k) for k in g.coverage([(lats, lons)])] + granules = [[{"id": "G", "s3": "s", "https": "h"}] for _ in keys] + return ShardMap(g.signature(), keys, granules, {"backend": "test"}) + + @pytest.fixture def antarctic_shardmap(): """ShardMap on an EPSG:3031 grid -> footprints entirely south of -60 deg.""" @@ -89,6 +108,7 @@ def arctic_shardmap(): # ── grid_from_signature ────────────────────────────────────────────────────── + class TestGridFromSignature: def test_rectilinear_round_trip(self, rect_grid): g = grid_from_signature(rect_grid.signature()) @@ -107,6 +127,7 @@ def test_unknown_type_raises(self): # ── shard_outlines ─────────────────────────────────────────────────────────── + class TestShardOutlines: def test_one_feature_per_shard(self, shardmap): fc = shard_outlines(shardmap) @@ -136,6 +157,7 @@ def test_valid_json(self, shardmap): # ── granule_footprints ─────────────────────────────────────────────────────── + class TestGranuleFootprints: def test_one_feature_per_granule(self, catalog): fc = granule_footprints(catalog) @@ -151,6 +173,7 @@ def test_geometry_is_polygon(self, catalog): # ── viewport_cells (grid-on-zoom gate) ─────────────────────────────────────── + class TestViewportCells: def test_gate_open_few_shards(self, shardmap): # Tight viewport over one chunk -> <= max_shards visible, grid drawn. @@ -238,21 +261,162 @@ def test_index_query_matches_exhaustive_scan(self, shardmap): assert indexed == exhaustive +# ── viewport_cells: HEALPix child-cell nesting + viewport bound ────────────── + + +class TestViewportCellsHealpix: + def _tight_view(self, grid, key, frac=0.3): + """A viewport centered in shard ``key``, ``frac`` of its bbox each side.""" + lon0, lat0, lon1, lat1 = grid.shard_footprint(key).bounds + cx, cy = (lon0 + lon1) / 2, (lat0 + lat1) / 2 + w, h = (lon1 - lon0) * frac / 2, (lat1 - lat0) * frac / 2 + return (cx - w, cy - h, cx + w, cy + h) + + def test_emits_child_order_cells_not_shard_outline(self, healpix_shardmap): + # The grid-on-zoom for HEALPix is the finer child-cell grid, so a single + # visible shard must yield *many* child cells -- not one shard-clip + # feature (the bug @espg reported: grid == shards, no nesting). + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + grid = grid_from_signature(healpix_shardmap.grid_signature) + key = healpix_shardmap.shard_keys[0] + view = self._tight_view(grid, key, frac=0.3) + fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=5000) + assert _is_geojson(fc) + assert len(fc["features"]) > 1 # not a lone shard-outline clip + + def test_cells_nest_within_parent_shard(self, healpix_shardmap): + # Every emitted cell's parent (clip2order at parent_order) is one of the + # visible shards, and the cell IDs are genuine order-child_order cells. + from mortie import clip2order, infer_order_from_morton + + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + grid = grid_from_signature(healpix_shardmap.grid_signature) + key = int(healpix_shardmap.shard_keys[0]) + view = self._tight_view(grid, key, frac=0.2) + fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=5000) + shard_keys = {int(k) for k in healpix_shardmap.shard_keys} + assert fc["features"] + for feat in fc["features"]: + cell = feat["properties"]["cell"] + parent = feat["properties"]["shard_key"] + assert int(infer_order_from_morton(cell)) == grid.child_order + assert int(clip2order(grid.parent_order, np.asarray([cell]))[0]) == parent + assert parent in shard_keys + + def test_cell_union_tiles_visible_shard(self, healpix_shardmap): + # A viewport covering exactly one whole shard: the union of the emitted + # child cells must reconstruct that shard footprint (4-for-1 nesting), + # confirming the grid lines up inside the shard outline. + from shapely.ops import unary_union + + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + # Use a smaller level diff so a whole-shard view stays cheap. + g = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere") + lats = np.array([10.0, 10.5, 10.5, 10.0, 10.0]) + lons = np.array([20.0, 20.0, 20.5, 20.5, 20.0]) + keys = [int(k) for k in g.coverage([(lats, lons)])] + sm = ShardMap(g.signature(), keys, [[] for _ in keys], {}) + key = keys[0] + shard_fp = g.shard_footprint(key) + # The shard's bbox overlaps its 4 neighbors (a HEALPix diamond's bbox is + # larger than the diamond), so allow them through the gate; cells are + # then filtered to this shard via the parent check below. + view = shard_fp.bounds + fc = viewport_cells(sm, view, max_shards=8, max_cells=5000) + cell_polys = [ + g.shard_footprint(f["properties"]["cell"]) + for f in fc["features"] + if f["properties"]["shard_key"] == key + ] + union = unary_union(cell_polys) + # Union of the shard's child cells reconstructs the shard footprint. + assert union.intersection(shard_fp).area / shard_fp.area > 0.98 + + def test_cell_count_scales_with_viewport_not_fan_out(self, healpix_shardmap): + # Zooming in (smaller viewport) must emit *fewer* cells. A naive + # 4^(child-parent) per-shard enumeration would emit a constant 4096/shard + # regardless of zoom; viewport-bounded coverage shrinks with the view. + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + grid = grid_from_signature(healpix_shardmap.grid_signature) + key = healpix_shardmap.shard_keys[0] + wide = viewport_cells( + healpix_shardmap, + self._tight_view(grid, key, 0.4), + max_shards=4, + max_cells=100000, + ) + tight = viewport_cells( + healpix_shardmap, + self._tight_view(grid, key, 0.05), + max_shards=4, + max_cells=100000, + ) + assert len(tight["features"]) < len(wide["features"]) + # And both are far below the full 4^(12-6)=4096 per-shard enumeration. + assert len(wide["features"]) < grid.n_children + + def test_no_full_enumeration_on_zoom(self, healpix_shardmap, monkeypatch): + # A zoomed-in query must not enumerate every child of a shard. Guard + # generate_morton_children: the viewport path uses morton_coverage, so a + # full child enumeration would be a regression. + import mortie + + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + + def boom(*a, **k): # pragma: no cover - must not be called + raise AssertionError("full child enumeration on a zoomed-in query") + + monkeypatch.setattr(mortie, "generate_morton_children", boom) + grid = grid_from_signature(healpix_shardmap.grid_signature) + key = healpix_shardmap.shard_keys[0] + fc = viewport_cells( + healpix_shardmap, + self._tight_view(grid, key, 0.1), + max_shards=4, + max_cells=5000, + ) + assert _is_geojson(fc) + + def test_max_cells_gate_returns_empty(self, healpix_shardmap): + # A whole-shard view at child_order 12 is tens of thousands of cells; + # the max_cells gate keeps a refresh bounded by returning empty rather + # than emitting them all (the dense-viewport bound). + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + grid = grid_from_signature(healpix_shardmap.grid_signature) + key = healpix_shardmap.shard_keys[0] + view = grid.shard_footprint(key).bounds + fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=200) + assert fc["features"] == [] + + # ── antimeridian splitting ─────────────────────────────────────────────────── + class TestAntimeridian: def test_healpix_shard_near_antimeridian_splits(self): # A HEALPix parent cell straddling +-180 -> MultiPolygon, each part in # one hemisphere (no globe-spanning band). grid = HealpixGrid(2, 6, layout="fullsphere") key = int( - grid.coverage( - [(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))] - )[0] + grid.coverage([(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))])[0] ) sm = ShardMap( - grid.signature(), [key], - [[{"id": "G", "s3": "s", "https": "h"}]], {}, + grid.signature(), + [key], + [[{"id": "G", "s3": "s", "https": "h"}]], + {}, ) fc = shard_outlines(sm) geom = fc["features"][0]["geometry"] @@ -276,8 +440,17 @@ def test_wide_non_crossing_polygon_kept_intact(self): from shapely.geometry import Polygon ring = [ - (-170, -10), (-85, -10), (0, -10), (85, -10), (170, -10), - (170, 10), (85, 10), (0, 10), (-85, 10), (-170, 10), (-170, -10), + (-170, -10), + (-85, -10), + (0, -10), + (85, -10), + (170, -10), + (170, 10), + (85, 10), + (0, 10), + (-85, 10), + (-170, 10), + (-170, -10), ] geom = _split_antimeridian(Polygon(ring)) assert geom["type"] == "Polygon" @@ -297,6 +470,7 @@ def test_holes_preserved(self): # ── render_shardmap assembly ───────────────────────────────────────────────── + class TestRenderShardmap: def test_shards_only(self, shardmap): out = render_shardmap(shardmap) @@ -325,6 +499,7 @@ def test_from_geoparquet_path(self, shardmap, catalog, tmp_path): # ── phase C: CRS selection (headless, no browser) ──────────────────────────── + class TestCrsSelection: def test_midlatitude_is_web_mercator(self, shardmap): # UTM 18N grid near 38.9N -> mid-latitude -> Web Mercator default. @@ -380,19 +555,20 @@ def test_bad_crs_raises(self): # ── phase C: CRS-aware antimeridian seam ───────────────────────────────────── + class TestSeamAwareLayers: def test_polar_skips_antimeridian_split(self): # A HEALPix cell straddling +-180 is a MultiPolygon under the Mercator # split, but stays a single Polygon when split_seam=False (polar CRS). grid = HealpixGrid(2, 6, layout="fullsphere") key = int( - grid.coverage( - [(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))] - )[0] + grid.coverage([(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))])[0] ) sm = ShardMap( - grid.signature(), [key], - [[{"id": "G", "s3": "s", "https": "h"}]], {}, + grid.signature(), + [key], + [[{"id": "G", "s3": "s", "https": "h"}]], + {}, ) split = shard_outlines(sm)["features"][0]["geometry"] unsplit = shard_outlines(sm, split_seam=False)["features"][0]["geometry"] @@ -405,16 +581,18 @@ def test_granule_footprints_seam_flag(self): # single ring -- so the flag, not the data, decides the geometry type. ring = [[178.0, 70.0], [-178.0, 70.0], [-178.0, 72.0], [178.0, 72.0], [178.0, 70.0]] item = { - "type": "Feature", "stac_version": "1.0.0", "id": "Gx", + "type": "Feature", + "stac_version": "1.0.0", + "id": "Gx", "geometry": {"type": "Polygon", "coordinates": [ring]}, "bbox": [-180.0, 70.0, 180.0, 72.0], "properties": {"datetime": "2025-06-01T00:00:00Z"}, - "collection": "TEST", "stac_extensions": [], "links": [], + "collection": "TEST", + "stac_extensions": [], + "links": [], "assets": {"data": {"href": "https://h/Gx.h5", "roles": ["data"]}}, } - cat = Catalog( - pa.table(sga.parse_stac_items_to_arrow([item])), {"collection": "TEST"} - ) + cat = Catalog(pa.table(sga.parse_stac_items_to_arrow([item])), {"collection": "TEST"}) split = granule_footprints(cat)["features"][0]["geometry"] unsplit = granule_footprints(cat, split_seam=False)["features"][0]["geometry"] assert split["type"] == "MultiPolygon" @@ -423,26 +601,53 @@ def test_granule_footprints_seam_flag(self): # ── debounce wrapper (pure stdlib; no widget stack) ────────────────────────── + class TestDebounce: - def test_rapid_calls_coalesce_to_one(self): + def test_rapid_calls_coalesce_to_one_on_loop(self): # Importing leaflet pulls only stdlib + zagg.viz at module level (the # ipyleaflet imports are local to the functions), so this is headless. - import time + # The debounce now schedules on the *running* event loop (the kernel's + # main thread) rather than a background threading.Timer -- that is the + # fix for the comm-thread crash (PR #44). Drive it on a loop and assert + # the burst coalesces to one call, fired on the loop thread. + import asyncio + import threading from zagg.viz.leaflet import _debounce + async def scenario(): + calls = {"n": 0, "tid": None} + + def hit(): + calls["n"] += 1 + calls["tid"] = threading.get_ident() + + deb = _debounce(0.05, hit) + for _ in range(20): + deb() # each cancels and reschedules the pending loop callback + assert calls["n"] == 0 # nothing fired yet (still within the window) + await asyncio.sleep(0.2) + return calls + + result = asyncio.run(scenario()) + assert result["n"] == 1 # the burst coalesced into a single refresh + assert result["tid"] == threading.get_ident() # ran on the loop's thread + + def test_no_loop_runs_synchronously(self): + # With no running event loop (plain script / headless), there is no comm + # to protect and no loop to schedule on, so the call runs inline. + from zagg.viz.leaflet import _debounce + calls = {"n": 0} deb = _debounce(0.05, lambda: calls.__setitem__("n", calls["n"] + 1)) - for _ in range(20): - deb() # each cancels and reschedules the prior timer - assert calls["n"] == 0 # nothing fired yet (still within the window) - time.sleep(0.2) - assert calls["n"] == 1 # the burst coalesced into a single refresh + deb() + assert calls["n"] == 1 deb.cancel() # ── phase 2: ipyleaflet wrapper (skips when the viz extra isn't installed) ──── + class TestShowShardmap: def test_import_core_without_ipyleaflet(self): # The headless core and zagg.viz import must not require ipyleaflet. From 888632183cd6cb8baf953256b135dda75f5a336f Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 22 Jun 2026 21:43:01 +0000 Subject: [PATCH 13/17] fold review: split antimeridian seam before clipping child cells (#38) --- src/zagg/viz/leaflet.py | 5 ++-- src/zagg/viz/shardmap.py | 57 +++++++++++++++++++++++++++++++++++++--- tests/test_viz.py | 51 +++++++++++++++++++++++++++++++++++ 3 files changed, 108 insertions(+), 5 deletions(-) diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index 90001d6..4938f01 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -234,8 +234,9 @@ def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) ) # Debounce so a burst of intermediate pan/zoom ticks coalesces into one - # refresh after movement settles. Keep a handle so the pending timer can be - # cancelled (exposed as ``m.cancel_grid_refresh``) for clean teardown. + # refresh after movement settles, scheduled on the kernel's event loop (main + # thread) -- never a background thread mutating the widget comm. Keep a handle + # so the pending callback can be cancelled (``m.cancel_grid_refresh``). debounced_refresh = _debounce(_GRID_DEBOUNCE_S, _refresh_grid) m.observe(debounced_refresh, names="bounds") m.cancel_grid_refresh = debounced_refresh.cancel diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py index 54c0031..be62fdc 100644 --- a/src/zagg/viz/shardmap.py +++ b/src/zagg/viz/shardmap.py @@ -156,6 +156,51 @@ def _unwrap(ring): return {"type": "MultiPolygon", "coordinates": parts} +def _seam_safe_polygon(geom): + """Cut a seam-crossing Polygon into hemisphere-local parts (shapely). + + Returns a shapely Polygon/MultiPolygon already in ``[-180, 180]`` -- the + same unwrap (+360) / clip / rewrap as :func:`_split_antimeridian`, but + returning geometry instead of a GeoJSON dict so callers can keep clipping or + measuring it. A polygon that doesn't cross the seam is returned unchanged. + + This must run **before** any ``box`` clip on a HEALPix child cell: a cell + straddling +-180 comes back from ``mort2polygon`` as a flat ring spanning + ~360 deg (e.g. ``-179.98 ... 180``), so intersecting that band with a + viewport box is meaningless -- the cell would vanish or mis-clip (PR #44). + """ + from shapely.geometry import MultiPolygon, Polygon, box + + if not _crosses_antimeridian(geom): + return geom + + def _unwrap(ring): + return [[lon + 360.0 if lon < 0 else lon, lat] for lon, lat in _ring_list(ring)] + + poly = Polygon(_unwrap(geom.exterior), [_unwrap(r) for r in geom.interiors]) + if not poly.is_valid: + poly = poly.buffer(0) + + parts: list = [] + for clip, shift in ( + (box(-180.0, -90.0, 180.0, 90.0), 0.0), + (box(180.0, -90.0, 540.0, 90.0), -360.0), + ): + piece = poly.intersection(clip) + for sub in getattr(piece, "geoms", [piece]): + if sub.is_empty or sub.geom_type != "Polygon": + continue + rings = [ + [[lon + shift, lat] for lon, lat in _ring_list(r)] + for r in (sub.exterior, *sub.interiors) + ] + parts.append(Polygon(rings[0], rings[1:])) + + if not parts: + return geom + return parts[0] if len(parts) == 1 else MultiPolygon(parts) + + def _polygon_rings(geom, *, shift: float = 0.0) -> list[list[list[float]]]: """GeoJSON ring list ``[exterior, *interiors]`` for a Polygon, lon-shifted.""" rings = [geom.exterior, *geom.interiors] @@ -383,7 +428,7 @@ def _child_cell_polygon(cell): return Polygon(zip(lons, lats)) -def _healpix_child_cells(grid, visible_keys, view, bbox, max_cells): +def _healpix_child_cells(grid, visible_keys, bbox, max_cells): """Child cells at ``child_order`` nesting inside the visible HEALPix shards. The grid-on-zoom for HEALPix is the **finer child grid that subdivides each @@ -486,13 +531,19 @@ def viewport_cells( if isinstance(grid, HealpixGrid): visible_keys = {int(key) for key, _ in visible} - result = _healpix_child_cells(grid, visible_keys, view, bbox, max_cells) + result = _healpix_child_cells(grid, visible_keys, bbox, max_cells) if result is None: return _collection([]) cells, parents = result features = [] for cell, parent in zip(cells, parents): - clipped = _polygonal(_child_cell_polygon(cell).intersection(view)) + poly = _child_cell_polygon(cell) + # Cut a seam-crossing cell at +-180 *before* clipping: its flat ring + # spans ~360 deg, so a box clip on the unsplit band drops/mis-clips + # it (PR #44). Skipped under a polar CRS (no seam there). + if split_seam: + poly = _seam_safe_polygon(poly) + clipped = _polygonal(poly.intersection(view)) if clipped is None: continue features.append( diff --git a/tests/test_viz.py b/tests/test_viz.py index 86a4954..d82f13c 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -400,6 +400,57 @@ def test_max_cells_gate_returns_empty(self, healpix_shardmap): fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=200) assert fc["features"] == [] + def test_child_cells_respect_antimeridian_seam(self): + # Child cells straddling +-180 split into hemisphere-local parts (no + # globe-spanning band) when split_seam=True, and stay single Polygons + # under a polar CRS (split_seam=False) -- same seam handling as shards. + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + g = HealpixGrid(parent_order=4, child_order=8, layout="fullsphere") + key = int( + g.coverage([(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))])[0] + ) + sm = ShardMap(g.signature(), [key], [[]], {}) + view = g.shard_footprint(key).bounds + split = viewport_cells(sm, view, max_shards=8, max_cells=5000, split_seam=True) + assert any(f["geometry"]["type"] == "MultiPolygon" for f in split["features"]) + for feat in split["features"]: + geom = feat["geometry"] + polys = geom["coordinates"] if geom["type"] == "MultiPolygon" else [geom["coordinates"]] + for poly in polys: + lons = [pt[0] for pt in poly[0]] + assert max(lons) - min(lons) <= 180.0 # no globe-spanning band + unsplit = viewport_cells(sm, view, max_shards=8, max_cells=5000, split_seam=False) + assert all(f["geometry"]["type"] == "Polygon" for f in unsplit["features"]) + + def test_seam_cell_clipped_to_true_sliver_not_global_band(self): + # Regression (PR #44 review): a child cell straddling +-180 must be split + # at the seam *before* clipping. The unsplit cell's flat ring spans + # ~360 deg; clipping that band to a seam-hugging viewport yields a wrong, + # oversized geometry (or drops the cell). After the fix, every emitted + # part is a small hemisphere-local sliver -- never a near-global band. + from shapely.geometry import shape + + from zagg.viz import shardmap as sm_mod + + sm_mod._INDEX_CACHE.clear() + g = HealpixGrid(parent_order=4, child_order=9, layout="fullsphere") + key = int( + g.coverage([(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))])[0] + ) + sm = ShardMap(g.signature(), [key], [[]], {}) + # A narrow viewport hugging the seam at the shard's latitude band. + lon0, lat0, lon1, lat1 = g.shard_footprint(key).bounds + view = (179.0, lat0, 180.0, (lat0 + lat1) / 2) + fc = viewport_cells(sm, view, max_shards=8, max_cells=20000, split_seam=True) + assert fc["features"] # the seam region is not silently emptied + # An order-9 cell footprint is ~0.009 deg^2. Clipping the unsplit flat + # band to the view instead fills the whole view strip (~0.19 deg^2, 20x + # bigger) -- so a small per-cell area bound pins the seam-first fix. + for feat in fc["features"]: + assert shape(feat["geometry"]).area < 0.05 # true cell, not a band + # ── antimeridian splitting ─────────────────────────────────────────────────── From 9a7bfd87cfa6674d8edeb1b52afe36c1078e4de3 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 23 Jun 2026 11:19:01 +0000 Subject: [PATCH 14/17] minimal viewer: drop cell-grid rendering (#38) --- notebooks/shardmap_viewer.ipynb | 141 ++++++------ src/zagg/viz/__init__.py | 25 ++- src/zagg/viz/leaflet.py | 108 +--------- src/zagg/viz/shardmap.py | 323 +--------------------------- tests/test_viz.py | 370 +------------------------------- 5 files changed, 119 insertions(+), 848 deletions(-) diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index df2e72c..29641ad 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -4,14 +4,58 @@ "cell_type": "markdown", "id": "cell-0", "metadata": {}, - "source": "# Shard-map viewer\n\nThis notebook demonstrates the `zagg.viz` shard-map viewer (issue\n[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\nATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\nLogin required — anonymous CMR queries are publicly accessible).\n\nThe viewer has two layers:\n\n- **Headless render core** (`zagg.viz.render_shardmap` and friends) — pure\n Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines,\n granule footprints, and viewport-clipped grid cells.\n- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) — builds an interactive\n basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n\n**Polar-aware projection.** `show_shardmap` picks the display CRS from the\nmap's extent: a polar AOI renders in NASA polar-stereographic (EPSG:3031\nAntarctic / EPSG:3413 Arctic) with a **GIBS** basemap — undistorted at the\npole — while mid-latitude AOIs stay on Web Mercator + OpenStreetMap.\n\n## Install\n\n```bash\npip install \"zagg[viz,catalog]\" # core + ipyleaflet widget + stac-geoparquet\n```\n\n**Requires an internet connection** for the anonymous CMR-STAC query in\nsection 1. No credentials are needed.\n\n## Areas of interest\n\nThe notebook queries two AOIs:\n\n1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) — January\n 2020. Dense ICESat-2 coverage near the pole; O(10–40) granules in two\n weeks. Renders in EPSG:3031 (Antarctic Polar Stereographic).\n2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) — June\n 2020. Second example showing the same pipeline on an Arctic AOI\n (EPSG:3413).\n\nBoth inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\nday one — see section 3 for the round-trip to disk and reload." + "source": [ + "# Shard-map viewer\n", + "\n", + "This notebook demonstrates the `zagg.viz` shard-map viewer (issue\n", + "[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\n", + "ATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\n", + "Login required \u2014 anonymous CMR queries are publicly accessible).\n", + "\n", + "The viewer has two layers:\n", + "\n", + "- **Headless render core** (`zagg.viz.render_shardmap` and friends) \u2014 pure\n", + " Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n", + " `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines\n", + " and granule footprints.\n", + "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) \u2014 builds an interactive\n", + " basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n", + "\n", + "**Polar-aware projection.** `show_shardmap` picks the display CRS from the\n", + "map's extent: a polar AOI renders in NASA polar-stereographic (EPSG:3031\n", + "Antarctic / EPSG:3413 Arctic) with a **GIBS** basemap \u2014 undistorted at the\n", + "pole \u2014 while mid-latitude AOIs stay on Web Mercator + OpenStreetMap.\n", + "\n", + "## Install\n", + "\n", + "```bash\n", + "pip install \"zagg[viz,catalog]\" # core + ipyleaflet widget + stac-geoparquet\n", + "```\n", + "\n", + "**Requires an internet connection** for the anonymous CMR-STAC query in\n", + "section 1. No credentials are needed.\n", + "\n", + "## Areas of interest\n", + "\n", + "The notebook queries two AOIs:\n", + "\n", + "1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) \u2014 January\n", + " 2020. Dense ICESat-2 coverage near the pole; O(10\u201340) granules in two\n", + " weeks. Renders in EPSG:3031 (Antarctic Polar Stereographic).\n", + "2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) \u2014 June\n", + " 2020. Second example showing the same pipeline on an Arctic AOI\n", + " (EPSG:3413).\n", + "\n", + "Both inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\n", + "day one \u2014 see section 3 for the round-trip to disk and reload." + ] }, { "cell_type": "markdown", "id": "cell-1", "metadata": {}, "source": [ - "## 1. Antarctic Peninsula — fetch real ATL06 granules from NASA CMR-STAC\n", + "## 1. Antarctic Peninsula \u2014 fetch real ATL06 granules from NASA CMR-STAC\n", "\n", "`CMRSource` speaks directly to NASA's CMR-STAC endpoint (`requests`).\n", "No Earthdata Login or credentials are needed for anonymous granule-metadata\n", @@ -25,7 +69,7 @@ "id": "cell-2", "metadata": {}, "outputs": [], - "source": "from zagg.catalog.sources import CMRSource, Query\n\n# Antarctic Peninsula AOI — two weeks in January 2020.\nAOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\nSTART_DATE = \"2020-01-01\"\nEND_DATE = \"2020-01-15\"\n\nquery = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START_DATE,\n end_date=END_DATE,\n region=AOI_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog = CMRSource().fetch(query)\nprint(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" + "source": "from zagg.catalog.sources import CMRSource, Query\n\n# Antarctic Peninsula AOI \u2014 two weeks in January 2020.\nAOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\nSTART_DATE = \"2020-01-01\"\nEND_DATE = \"2020-01-15\"\n\nquery = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START_DATE,\n end_date=END_DATE,\n region=AOI_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog = CMRSource().fetch(query)\nprint(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" }, { "cell_type": "code", @@ -34,7 +78,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Inspect the catalog schema — each row is a STAC Item with WKB geometry.\n", + "# Inspect the catalog schema \u2014 each row is a STAC Item with WKB geometry.\n", "print(\"Schema columns:\", catalog.table.schema.names[:10], \"...\")\n", "print(\"Total rows:\", catalog.table.num_rows)\n", "\n", @@ -52,14 +96,14 @@ "source": [ "## 2. Build a ShardMap on a HEALPix grid\n", "\n", - "We use a HEALPix grid (`parent_order=6, child_order=12`) — the same\n", + "We use a HEALPix grid (`parent_order=6, child_order=12`) \u2014 the same\n", "configuration as `src/zagg/configs/atl06.yaml`. The `mortie` backend is\n", "used for HEALPix intersection and requires no extra install. If `spherely`\n", "is installed, pass `backend='auto'` to prefer exact S2 intersections.\n", "\n", "`ShardMap.build` maps each grid shard (a parent-order HEALPix cell) to the\n", "set of granules whose footprint intersects it. The manifest is self-contained\n", - "— it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", + "\u2014 it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", "never needs the Catalog again at run time." ] }, @@ -97,7 +141,7 @@ "# Inspect a few shard assignments.\n", "for key, shard_granules in zip(shardmap.shard_keys[:4], shardmap.granules[:4]):\n", " ids = [g[\"id\"] for g in shard_granules]\n", - " print(f\" shard {key:6d}: {len(ids)} granule(s) — {ids[:2]}{'...' if len(ids) > 2 else ''}\")" + " print(f\" shard {key:6d}: {len(ids)} granule(s) \u2014 {ids[:2]}{'...' if len(ids) > 2 else ''}\")" ] }, { @@ -149,9 +193,8 @@ "## 4. Headless render core: GeoJSON FeatureCollections\n", "\n", "`render_shardmap` assembles every layer into one dict of GeoJSON\n", - "`FeatureCollection`s. Pass a `catalog` to include granule footprints and\n", - "a `bbox` to include viewport-clipped grid cells. Each value is plain,\n", - "JSON-serializable GeoJSON — no widgets required." + "`FeatureCollection`s. Pass a `catalog` to include granule footprints. Each\n", + "value is plain, JSON-serializable GeoJSON \u2014 no widgets required." ] }, { @@ -198,74 +241,52 @@ }, { "cell_type": "markdown", - "id": "cell-13", + "id": "cell-15", "metadata": {}, "source": [ - "### Grid-on-zoom gate\n", + "## 5. Interactive map \u2014 Antarctic Peninsula (ipyleaflet, polar-stereographic)\n", "\n", - "`viewport_cells` draws shard-order cell outlines clipped to a viewport but\n", - "**only** when `<= max_shards` shards intersect the viewport. Zoom in and\n", - "the grid appears; zoom out and it returns an empty collection rather than a\n", - "globe-spanning graticule. The interactive map wires this to the live map\n", - "`bounds`, so the grid blinks on as you zoom in." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cell-14", - "metadata": {}, - "outputs": [], - "source": [ - "from zagg.viz import viewport_cells\n", - "from zagg.viz.shardmap import grid_from_signature\n", + "`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\n", + "Catalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\n", + "as file paths or in-memory objects.\n", + "\n", + "**Projection is auto-selected from the map's extent.** This Antarctic AOI is\n", + "entirely south of \u221260\u00b0, so the viewer renders in **EPSG:3031** (Antarctic\n", + "Polar Stereographic) with a NASA **GIBS** polar basemap instead of the\n", + "distorted Web Mercator default. An Arctic AOI (poleward of +60\u00b0) gets\n", + "**EPSG:3413**; mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Pass\n", + "`crs=\"EPSG:3857\"` to force Mercator, or `crs=\"EPSG:3413\"`/`\"EPSG:3031\"` to\n", + "force a specific pole. Vector layers stay WGS84 GeoJSON \u2014 proj4leaflet\n", + "reprojects them client-side.\n", "\n", - "# Tight viewport over one shard -> gate open, cells drawn.\n", - "first_shard = shardmap.shard_keys[0]\n", - "fp = grid_from_signature(shardmap.grid_signature).shard_footprint(first_shard)\n", - "lon0, lat0, lon1, lat1 = fp.bounds\n", + "**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\n", + "Under headless `nbconvert` the Map object is constructed (no error) but\n", + "tiles won't display.\n", "\n", - "zoomed_in = viewport_cells(shardmap, (lon0, lat0, lon1, lat1), max_shards=4)\n", - "zoomed_out = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=1)\n", + "### Verification checklist\n", "\n", - "print(\"zoomed-in cells:\", len(zoomed_in[\"features\"]))\n", - "print(\"zoomed-out cells (gated):\", len(zoomed_out[\"features\"]))" + "1. **Basemap** \u2014 the GIBS polar basemap renders, pans, zooms, and the\n", + " continent is undistorted at the pole (no Web-Mercator stretching).\n", + "2. **Shard outlines** \u2014 blue polygons over the Peninsula; click one for\n", + " its `shard_key` and `n_granules` properties.\n", + "3. **Granule footprints toggle** \u2014 layer-control (top-right); the ICESat-2\n", + " track footprints appear/disappear." ] }, - { - "cell_type": "markdown", - "id": "cell-15", - "metadata": {}, - "source": "## 5. Interactive map — Antarctic Peninsula (ipyleaflet, polar-stereographic)\n\n`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\nCatalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\nas file paths or in-memory objects.\n\n**Projection is auto-selected from the map's extent.** This Antarctic AOI is\nentirely south of −60°, so the viewer renders in **EPSG:3031** (Antarctic\nPolar Stereographic) with a NASA **GIBS** polar basemap instead of the\ndistorted Web Mercator default. An Arctic AOI (poleward of +60°) gets\n**EPSG:3413**; mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Pass\n`crs=\"EPSG:3857\"` to force Mercator, or `crs=\"EPSG:3413\"`/`\"EPSG:3031\"` to\nforce a specific pole. Vector layers stay WGS84 GeoJSON — proj4leaflet\nreprojects them client-side.\n\n**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\nUnder headless `nbconvert` the Map object is constructed (no error) but\ntiles won't display.\n\n### Verification checklist\n\n1. **Basemap** — the GIBS polar basemap renders, pans, zooms, and the\n continent is undistorted at the pole (no Web-Mercator stretching).\n2. **Shard outlines** — blue polygons over the Peninsula; click one for\n its `shard_key` and `n_granules` properties.\n3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n track footprints appear/disappear.\n4. **Grid-on-zoom** — see the next cell, which centers the map on a single\n shard so the grey shard-cell grid actually appears." - }, { "cell_type": "code", "execution_count": null, "id": "cell-16", "metadata": {}, "outputs": [], - "source": "from zagg.viz import show_shardmap\n\n# File-path interface — same as after `python -m zagg.catalog ...`.\n# CRS auto-selected from the AOI extent: this Antarctic map renders in EPSG:3031.\nm = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\nprint(\"display CRS:\", m.crs[\"name\"])\nprint(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\nm" - }, - { - "cell_type": "markdown", - "id": "a326b760", - "source": "### Zoom to one shard so the grid appears\n\nThe grid layer is gated on `<= max_shards` (default 4) shards intersecting the\nviewport — over the whole AOI far more than 4 shards are visible, so the gate\nstays shut and you only see the shard outlines. Centering the map on a single\nshard's bbox and zooming in puts `<= 4` shards in view, which fires the\n`Map.bounds` observer and makes `viewport_cells` emit the grey shard-cell grid.\nRun this cell, then watch the **grid (shard cells)** layer turn on.", - "metadata": {} - }, - { - "cell_type": "code", - "id": "f497f2c3", - "source": "from zagg.viz.shardmap import grid_from_signature\n\n# Center on a single shard's bbox so <= max_shards fill the viewport.\none_shard = shardmap.shard_keys[0]\nslon0, slat0, slon1, slat1 = grid_from_signature(shardmap.grid_signature).shard_footprint(\n one_shard\n).bounds\nshard_center = ((slat0 + slat1) / 2, (slon0 + slon1) / 2) # (lat, lon)\n\nm.center = shard_center\nm.zoom = 9 # tight enough that <= 4 shards intersect -> grid gate opens\n\n# The grid layer is driven off m.bounds; setting center/zoom re-fires the\n# observer. In a live JupyterLab session the grey \"grid (shard cells)\" layer\n# now appears. (Under nbconvert m.bounds stays empty until the browser sizes\n# the map, so the grid only renders interactively.)\nprint(f\"centered on shard {one_shard} at {shard_center}, zoom={m.zoom}\")\nm", - "metadata": {}, - "execution_count": null, - "outputs": [] + "source": "from zagg.viz import show_shardmap\n\n# File-path interface \u2014 same as after `python -m zagg.catalog ...`.\n# CRS auto-selected from the AOI extent: this Antarctic map renders in EPSG:3031.\nm = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\nprint(\"display CRS:\", m.crs[\"name\"])\nprint(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\nm" }, { "cell_type": "markdown", "id": "cell-17", "metadata": {}, "source": [ - "## 6. Second AOI — Jakobshavn Glacier, West Greenland\n", + "## 6. Second AOI \u2014 Jakobshavn Glacier, West Greenland\n", "\n", "A second example using an Arctic AOI to show the same pipeline working\n", "independently on a different region. Jakobshavn Glacier is one of\n", @@ -278,7 +299,7 @@ "id": "cell-18", "metadata": {}, "outputs": [], - "source": "# Imports also available from earlier cells; repeated here for standalone re-runs.\nfrom zagg.viz import render_shardmap, show_shardmap\n# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render — confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" + "source": "# Imports also available from earlier cells; repeated here for standalone re-runs.\nfrom zagg.viz import render_shardmap, show_shardmap\n# Second AOI: Jakobshavn Glacier, West Greenland \u2014 June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render \u2014 confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" }, { "cell_type": "code", @@ -302,4 +323,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/src/zagg/viz/__init__.py b/src/zagg/viz/__init__.py index 251dd2d..c1f67b7 100644 --- a/src/zagg/viz/__init__.py +++ b/src/zagg/viz/__init__.py @@ -2,22 +2,23 @@ Two layers: -- :mod:`zagg.viz.shardmap` -- the **headless render core** (phase 1). Pure - Python, no browser or ipyleaflet import: it turns a :class:`ShardMap` (and an - optional :class:`~zagg.catalog.sources.Catalog`) into GeoJSON - ``FeatureCollection`` dicts -- shard/chunk outlines, granule footprints, and - viewport-clipped cell outlines at the shard order. Fully unit-testable with +- :mod:`zagg.viz.shardmap` -- the **headless render core**. Pure Python, no + browser or ipyleaflet import: it turns a :class:`ShardMap` (and an optional + :class:`~zagg.catalog.sources.Catalog`) into GeoJSON ``FeatureCollection`` + dicts -- shard/chunk outlines and granule footprints. Fully unit-testable with no widget stack installed. -- :func:`show_shardmap` -- the **ipyleaflet wrapper** (phase 2). Builds an - interactive map from a saved ShardMap (basemap + shard layer + a toggleable - footprint layer + a zoom-thresholded grid layer). ``ipyleaflet`` is imported - lazily inside the widget functions so phase-1 core and the test suite never - require it; it lives in the optional ``viz`` extra (``pip install zagg[viz]``). +- :func:`show_shardmap` -- the **ipyleaflet wrapper**. Builds an interactive map + from a saved ShardMap (context basemap with auto polar-projection switching + + shard layer + a toggleable granule-footprint layer). ``ipyleaflet`` is + imported lazily inside the widget functions so the render core and the test + suite never require it; it lives in the optional ``viz`` extra + (``pip install zagg[viz]``). Both inputs are reused directly off the existing surface -- ``ShardMap`` / ``Catalog`` round-trips and per-grid ``shard_footprint`` / ``signature`` -- so there is no viewer-specific file type or second tessellation (issue #38). """ + from __future__ import annotations from zagg.viz.shardmap import ( @@ -25,12 +26,11 @@ grid_from_signature, render_shardmap, shard_outlines, - viewport_cells, ) def show_shardmap(shardmap_path, catalog=None, **kwargs): - """Build an interactive ipyleaflet map for a saved ShardMap (phase 2). + """Build an interactive ipyleaflet map for a saved ShardMap. Thin lazy passthrough to :func:`zagg.viz.leaflet.show_shardmap` so importing :mod:`zagg.viz` never pulls in ``ipyleaflet`` -- only calling this does. @@ -59,7 +59,6 @@ def show_shardmap(shardmap_path, catalog=None, **kwargs): "render_shardmap", "shard_outlines", "granule_footprints", - "viewport_cells", "grid_from_signature", "show_shardmap", ] diff --git a/src/zagg/viz/leaflet.py b/src/zagg/viz/leaflet.py index 4938f01..3e2c021 100644 --- a/src/zagg/viz/leaflet.py +++ b/src/zagg/viz/leaflet.py @@ -1,10 +1,7 @@ -"""ipyleaflet wrapper for the shard-map viewer (issue #38, phases 2 & C). +"""ipyleaflet wrapper for the shard-map viewer (issue #38). Builds an interactive map from a saved ShardMap: a basemap, the shard-outline -layer, an optional (toggleable) granule-footprint layer, and a grid layer that -draws shard-order cell outlines only when the viewport is zoomed in far enough -to show ``<= max_shards`` shards (the "grid-on-zoom" gate -- never a global -graticule, issue #38). +layer, and an optional (toggleable) granule-footprint layer. The CRS is chosen from the map's extent (:mod:`zagg.viz.crs`): polar AOIs get a NASA polar-stereographic projection (EPSG:3413/3031) with a matching **GIBS** @@ -20,74 +17,17 @@ from __future__ import annotations -import asyncio - from zagg.viz.crs import crs_info, is_polar, pick_crs from zagg.viz.shardmap import ( _load_catalog, _load_shardmap, granule_footprints, - shard_index, shard_outlines, - viewport_cells, ) # Layer styles (kept terse; tweakable by callers via the returned Map). _SHARD_STYLE = {"color": "#1f78b4", "weight": 1, "fillOpacity": 0.05} _FOOTPRINT_STYLE = {"color": "#e31a1c", "weight": 1, "fillOpacity": 0.10} -_GRID_STYLE = {"color": "#333333", "weight": 1, "fillOpacity": 0.0} - -# Coalesce a burst of pan/zoom ``bounds`` ticks into one grid refresh after the -# viewport settles (seconds). Without this the grid recomputes on every -# intermediate frame of a drag. -_GRID_DEBOUNCE_S = 0.2 - - -def _debounce(wait, func): - """Wrap ``func`` so rapid calls coalesce to one call ``wait`` s after the last. - - The coalesced call is scheduled on the **kernel's own asyncio event loop** - (``loop.call_later``) -- the same main thread that owns the ipywidgets comm - -- not on a background ``threading.Timer``. That distinction is the fix for - the kernel-crash report (PR #44): the Jupyter widget comm channel is not - thread-safe, so mutating a widget traitlet (here ``grid_layer.data``) from a - timer thread can hang or corrupt the comm and crash the kernel. Scheduling on - the running loop keeps every refresh on the main thread. - - Each call cancels the pending :class:`asyncio.TimerHandle` and reschedules, - so a burst of events fires ``func`` once when they stop. Returns the wrapper - with a ``cancel()`` to tear the pending call down. - - When no event loop is running (e.g. a plain script / headless test) the call - is run synchronously -- there is no comm to protect and no loop to schedule - on, so coalescing is moot. - """ - handle: list = [None] - - def _fire(args, kwargs): - handle[0] = None - func(*args, **kwargs) - - def wrapper(*args, **kwargs): - if handle[0] is not None: - handle[0].cancel() - handle[0] = None - try: - loop = asyncio.get_running_loop() - except RuntimeError: - loop = None - if loop is None: - func(*args, **kwargs) - return - handle[0] = loop.call_later(wait, _fire, args, kwargs) - - def cancel(): - if handle[0] is not None: - handle[0].cancel() - handle[0] = None - - wrapper.cancel = cancel - return wrapper def _center_zoom(fc: dict): @@ -133,7 +73,6 @@ def show_shardmap( shardmap_path, catalog=None, *, - max_shards: int = 4, zoom: int = 3, basemap=None, crs=None, @@ -152,9 +91,7 @@ def show_shardmap( Path to a ``ShardMap`` JSON file (or an in-memory ``ShardMap``). catalog : str or Catalog, optional A geoparquet path or a loaded ``Catalog``. When given, a toggleable - granule-footprint layer is added (off by default in the layer control). - max_shards : int - Visible-shard gate for the grid-on-zoom layer. + granule-footprint layer is added. zoom : int Initial map zoom. basemap : ipyleaflet basemap, optional @@ -165,8 +102,8 @@ def show_shardmap( Returns ------- ipyleaflet.Map - Map with a shard layer, optional footprint layer, a zoom-thresholded - grid layer, and a ``LayersControl`` for toggling layers. + Map with a shard layer, an optional footprint layer, and a + ``LayersControl`` for toggling layers. """ # Import first so a missing `viz` extra fails clearly, before any work. from ipyleaflet import GeoJSON, LayersControl, Map, TileLayer, basemaps @@ -207,41 +144,6 @@ def show_shardmap( ) m.add(footprint_layer) - # Grid-on-zoom: an empty layer kept in sync with the viewport. The shard - # footprints are indexed once (STRtree) so each refresh is a cheap query, not - # a full footprint rebuild; the headless core's gate returns nothing when too - # many shards are visible, so this never becomes a global graticule. - grid_layer = GeoJSON( - data={"type": "FeatureCollection", "features": []}, - style=_GRID_STYLE, - name="grid (shard cells)", - ) - m.add(grid_layer) - - index = shard_index(shardmap) - - def _refresh_grid(event=None): # noqa: ARG001 (traitlets observe signature) - bounds = m.bounds - if not bounds: - return - (south, west), (north, east) = bounds - grid_layer.data = viewport_cells( - shardmap, - (west, south, east, north), - max_shards=max_shards, - split_seam=split_seam, - index=index, - ) - - # Debounce so a burst of intermediate pan/zoom ticks coalesces into one - # refresh after movement settles, scheduled on the kernel's event loop (main - # thread) -- never a background thread mutating the widget comm. Keep a handle - # so the pending callback can be cancelled (``m.cancel_grid_refresh``). - debounced_refresh = _debounce(_GRID_DEBOUNCE_S, _refresh_grid) - m.observe(debounced_refresh, names="bounds") - m.cancel_grid_refresh = debounced_refresh.cancel - _refresh_grid() - m.add(LayersControl(position="topright")) return m diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py index be62fdc..7419454 100644 --- a/src/zagg/viz/shardmap.py +++ b/src/zagg/viz/shardmap.py @@ -1,11 +1,11 @@ -"""Headless render core for the shard-map viewer (issue #38, phase 1). +"""Headless render core for the shard-map viewer (issue #38). Pure Python: turns a :class:`~zagg.catalog.shardmap.ShardMap` (and an optional :class:`~zagg.catalog.sources.Catalog`) into GeoJSON ``FeatureCollection`` dicts in WGS84. No browser, no ipyleaflet -- everything here is unit-testable with just the core deps (``shapely`` + the grid backends). -Three layers are produced: +Two layers are produced: - :func:`shard_outlines` -- one polygon feature per shard, straight off ``grid.shard_footprint(key)``. The grid is reconstructed from the map's own @@ -13,12 +13,6 @@ needed. - :func:`granule_footprints` -- one polygon feature per granule footprint, decoded from a ``Catalog`` (``granule_records``). -- :func:`viewport_cells` -- the grid-on-zoom layer clipped to a viewport bbox, - emitted **only** when ``<= max_shards`` shards intersect the viewport (the - gate -- never a global graticule, issue #38). For HEALPix this is the finer - **child cells at ``child_order``** nesting inside the visible shards (generated - viewport-bounded via ``mortie.morton_coverage``, not by per-shard fan-out); - for other grids it is the shard footprint clipped to the viewport. Antimeridian handling --------------------- @@ -44,9 +38,9 @@ def grid_from_signature(signature: dict): """Reconstruct an output grid from a ``ShardMap.grid_signature``. - The viewer only needs ``shard_footprint`` / ``children`` off the grid, both - of which are fully determined by the signature -- so the map is - self-describing and no separate config is required. + The viewer only needs ``shard_footprint`` off the grid, which is fully + determined by the signature -- so the map is self-describing and no separate + config is required. Parameters ---------- @@ -156,51 +150,6 @@ def _unwrap(ring): return {"type": "MultiPolygon", "coordinates": parts} -def _seam_safe_polygon(geom): - """Cut a seam-crossing Polygon into hemisphere-local parts (shapely). - - Returns a shapely Polygon/MultiPolygon already in ``[-180, 180]`` -- the - same unwrap (+360) / clip / rewrap as :func:`_split_antimeridian`, but - returning geometry instead of a GeoJSON dict so callers can keep clipping or - measuring it. A polygon that doesn't cross the seam is returned unchanged. - - This must run **before** any ``box`` clip on a HEALPix child cell: a cell - straddling +-180 comes back from ``mort2polygon`` as a flat ring spanning - ~360 deg (e.g. ``-179.98 ... 180``), so intersecting that band with a - viewport box is meaningless -- the cell would vanish or mis-clip (PR #44). - """ - from shapely.geometry import MultiPolygon, Polygon, box - - if not _crosses_antimeridian(geom): - return geom - - def _unwrap(ring): - return [[lon + 360.0 if lon < 0 else lon, lat] for lon, lat in _ring_list(ring)] - - poly = Polygon(_unwrap(geom.exterior), [_unwrap(r) for r in geom.interiors]) - if not poly.is_valid: - poly = poly.buffer(0) - - parts: list = [] - for clip, shift in ( - (box(-180.0, -90.0, 180.0, 90.0), 0.0), - (box(180.0, -90.0, 540.0, 90.0), -360.0), - ): - piece = poly.intersection(clip) - for sub in getattr(piece, "geoms", [piece]): - if sub.is_empty or sub.geom_type != "Polygon": - continue - rings = [ - [[lon + shift, lat] for lon, lat in _ring_list(r)] - for r in (sub.exterior, *sub.interiors) - ] - parts.append(Polygon(rings[0], rings[1:])) - - if not parts: - return geom - return parts[0] if len(parts) == 1 else MultiPolygon(parts) - - def _polygon_rings(geom, *, shift: float = 0.0) -> list[list[list[float]]]: """GeoJSON ring list ``[exterior, *interiors]`` for a Polygon, lon-shifted.""" rings = [geom.exterior, *geom.interiors] @@ -248,27 +197,6 @@ def _polygon_geometry(geom, *, split_seam: bool = True) -> dict: raise ValueError(f"unsupported geometry type for GeoJSON: {geom.geom_type}") -def _polygonal(geom): - """Reduce an intersection result to Polygon/MultiPolygon, or None. - - A clip can degenerate to a point/line on a shared edge or yield a - ``GeometryCollection`` mixing dimensions; keep only the polygonal part. - """ - if geom.is_empty: - return None - if geom.geom_type in ("Polygon", "MultiPolygon"): - return geom - if geom.geom_type == "GeometryCollection": - from shapely.geometry import MultiPolygon - - polys = [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")] - if not polys: - return None - merged = MultiPolygon([p for g in polys for p in getattr(g, "geoms", [g])]) - return merged - return None - - def _feature(geometry: dict, properties: dict) -> dict: return {"type": "Feature", "geometry": geometry, "properties": properties} @@ -348,231 +276,11 @@ def granule_footprints(catalog, *, split_seam: bool = True) -> dict: return _collection(features) -class ShardIndex: - """One-time shard-footprint index for fast viewport queries (issue #38). - - Building a shard footprint is not free -- on HEALPix it is - ``mortie.tools.mort2polygon`` per key. The interactive viewer re-queries on - every pan/zoom, so regenerating all footprints per event freezes the kernel. - This builds the ``(key, footprint)`` list **once** and indexes the - footprints in a shapely ``STRtree``; :meth:`query` then answers a viewport - in ``O(log N + k)`` instead of ``O(N)`` footprint rebuilds. - """ - - def __init__(self, shardmap): - from shapely import STRtree - - grid = grid_from_signature(shardmap.grid_signature) - self.keys = list(shardmap.shard_keys) - self.footprints = [grid.shard_footprint(k) for k in self.keys] - self._tree = STRtree(self.footprints) - - def query(self, view): - """``[(key, footprint)]`` whose footprint truly intersects ``view``. - - The ``STRtree`` returns bbox-overlap candidates; each is refined with a - precise ``footprint.intersects(view)`` so the result matches the old - exhaustive scan exactly. - """ - out = [] - for i in self._tree.query(view): - fp = self.footprints[i] - if fp.intersects(view): - out.append((self.keys[i], fp)) - return out - - -# Cache one ShardIndex per shardmap object so repeated viewport queries (and the -# back-compat viewport_cells path) never rebuild footprints. Keyed by id() and -# guarded by a weak ref so a collected shardmap drops its entry. -_INDEX_CACHE: dict = {} - - -def shard_index(shardmap) -> ShardIndex: - """Return a cached :class:`ShardIndex` for ``shardmap`` (built once).""" - import weakref - - key = id(shardmap) - entry = _INDEX_CACHE.get(key) - if entry is not None: - return entry[0] - - index = ShardIndex(shardmap) - - def _drop(_ref, key=key): - _INDEX_CACHE.pop(key, None) - - try: - ref = weakref.ref(shardmap, _drop) - except TypeError: - ref = None # not weak-referenceable; entry simply lingers - _INDEX_CACHE[key] = (index, ref) - return index - - -# Edge samples per side for a child-cell polygon. Child cells are small -# near-squares, so a low step traces them faithfully while keeping each ring -- -# and the GeoJSON shipped over the widget comm -- light. (HealpixGrid.shard_ -# footprint uses step=32 for the much larger parent diamonds.) -_CHILD_CELL_STEP = 8 - - -def _child_cell_polygon(cell): - """Polygon (WGS84 lon/lat) for a single HEALPix child cell morton id.""" - from mortie.tools import mort2polygon - from shapely.geometry import Polygon - - verts = mort2polygon(int(cell), step=_CHILD_CELL_STEP) - lats = np.array([v[0] for v in verts]) - lons = np.array([v[1] for v in verts]) - return Polygon(zip(lons, lats)) - - -def _healpix_child_cells(grid, visible_keys, bbox, max_cells): - """Child cells at ``child_order`` nesting inside the visible HEALPix shards. - - The grid-on-zoom for HEALPix is the **finer child grid that subdivides each - shard** (issue #38) -- not the shard outline redrawn. Children are generated - in a **viewport-bounded** way: a single ``mortie.morton_coverage`` query over - the viewport bbox at ``child_order`` returns only the child cells overlapping - the view, so the count scales with the viewport, not with the - ``4^(child_order - parent_order)`` per-shard fan-out. The returned cells are - filtered to those whose parent (``clip2order(parent_order, ...)``) is one of - the visible shards, so the rendered grid lines up *inside* the shard outlines. - - Returns ``(cells, parents)`` morton-id arrays, or ``None`` when the in-view - child count exceeds ``max_cells`` (the analogous gate to ``max_shards`` -- - keeps a single dense-viewport refresh bounded rather than emitting tens of - thousands of polygons). - """ - from mortie import clip2order, morton_coverage - - lon_min, lat_min, lon_max, lat_max = bbox - view_lats = np.array([lat_min, lat_min, lat_max, lat_max, lat_min]) - view_lons = np.array([lon_min, lon_max, lon_max, lon_min, lon_min]) - cov = np.asarray(morton_coverage([view_lats], [view_lons], order=grid.child_order)) - if cov.size == 0: - return np.empty(0, dtype=np.int64), np.empty(0, dtype=np.int64) - parents = clip2order(grid.parent_order, cov) - keep = np.isin(parents, np.asarray(list(visible_keys))) - cells = cov[keep] - if cells.size > max_cells: - return None - return cells, parents[keep] - - -def viewport_cells( - shardmap, - bbox, - *, - max_shards: int = 4, - max_cells: int = 2000, - split_seam: bool = True, - index=None, -) -> dict: - """Grid-on-zoom cell outlines clipped to a viewport, gated on visible shards. - - Implements the "grid-on-zoom" behavior (issue #38): the finer grid is drawn - only when ``<= max_shards`` shards intersect ``bbox``. When more shards are - visible the viewport is too zoomed-out for a useful grid and an empty - collection is returned -- never a global graticule. - - For a **HEALPix** grid the emitted features are the **child cells at - ``child_order`` that nest inside the visible shards** -- the grid that - subdivides each shard 4-for-1 per order step -- not the shard outline - redrawn. They are generated viewport-bounded via a single - :func:`mortie.morton_coverage` query at ``child_order`` (so the cell count - scales with the viewport, not with the ``4^(child_order - parent_order)`` - per-shard fan-out), filtered to the visible shards, and clipped to ``bbox``. - For any other grid the shard footprint clipped to the viewport is emitted. - - Footprints are taken from a one-time :class:`ShardIndex` (built and cached - per shardmap, or passed in via ``index=``), so repeated calls -- e.g. the - interactive viewer's per-pan refresh -- never rebuild them. - - Parameters - ---------- - shardmap : ShardMap - A built or loaded shard map. - bbox : tuple of float - Viewport ``(lon_min, lat_min, lon_max, lat_max)`` in WGS84. - max_shards : int - Maximum number of intersecting shards for the grid to render. - max_cells : int - Maximum number of in-view HEALPix child cells to render. Above this the - viewport is too dense (zoomed out relative to ``child_order``) and an - empty collection is returned, keeping a single refresh bounded. - split_seam : bool - Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a - polar-stereographic CRS, where there is no antimeridian seam. - index : ShardIndex, optional - A prebuilt index to query. When omitted, the cached index for - ``shardmap`` is used (:func:`shard_index`). - - Returns - ------- - dict - GeoJSON ``FeatureCollection`` of grid outlines clipped to the viewport, - or an empty collection when a gate is not met. - """ - from shapely.geometry import box - - if index is None: - index = shard_index(shardmap) - view = box(bbox[0], bbox[1], bbox[2], bbox[3]) - - visible = index.query(view) - if not visible or len(visible) > max_shards: - return _collection([]) - - grid = grid_from_signature(shardmap.grid_signature) - - from zagg.grids import HealpixGrid - - if isinstance(grid, HealpixGrid): - visible_keys = {int(key) for key, _ in visible} - result = _healpix_child_cells(grid, visible_keys, bbox, max_cells) - if result is None: - return _collection([]) - cells, parents = result - features = [] - for cell, parent in zip(cells, parents): - poly = _child_cell_polygon(cell) - # Cut a seam-crossing cell at +-180 *before* clipping: its flat ring - # spans ~360 deg, so a box clip on the unsplit band drops/mis-clips - # it (PR #44). Skipped under a polar CRS (no seam there). - if split_seam: - poly = _seam_safe_polygon(poly) - clipped = _polygonal(poly.intersection(view)) - if clipped is None: - continue - features.append( - _feature( - _polygon_geometry(clipped, split_seam=split_seam), - {"cell": int(cell), "shard_key": int(parent)}, - ) - ) - return _collection(features) - - features = [] - for key, fp in visible: - clipped = _polygonal(fp.intersection(view)) - if clipped is None: - continue - features.append( - _feature( - _polygon_geometry(clipped, split_seam=split_seam), - {"shard_key": _jsonable(key)}, - ) - ) - return _collection(features) - - # ── top-level assembly ─────────────────────────────────────────────────────── -def render_shardmap(shardmap, catalog=None, *, bbox=None, max_shards: int = 4) -> dict: - """Assemble all viewer layers for a shard map into one dict of collections. +def render_shardmap(shardmap, catalog=None) -> dict: + """Assemble the viewer layers for a shard map into one dict of collections. Parameters ---------- @@ -581,25 +289,17 @@ def render_shardmap(shardmap, catalog=None, *, bbox=None, max_shards: int = 4) - catalog : Catalog or str, optional A ``Catalog`` or a geoparquet path. When given, the granule-footprint layer is included. - bbox : tuple of float, optional - Viewport ``(lon_min, lat_min, lon_max, lat_max)``. When given, the - viewport-clipped cell-outline layer is included. - max_shards : int - Visible-shard gate for the viewport cell layer. Returns ------- dict - ``{"shards": FC, "granules": FC | None, "cells": FC | None}`` where each - value is a GeoJSON ``FeatureCollection`` (or ``None`` when its input was - not provided). + ``{"shards": FC, "granules": FC | None}`` where each value is a GeoJSON + ``FeatureCollection`` (or ``None`` when its input was not provided). """ shardmap = _load_shardmap(shardmap) - out = {"shards": shard_outlines(shardmap), "granules": None, "cells": None} + out = {"shards": shard_outlines(shardmap), "granules": None} if catalog is not None: out["granules"] = granule_footprints(_load_catalog(catalog)) - if bbox is not None: - out["cells"] = viewport_cells(shardmap, bbox, max_shards=max_shards) return out @@ -647,8 +347,5 @@ def _is_geojson(obj) -> bool: "grid_from_signature", "shard_outlines", "granule_footprints", - "viewport_cells", - "ShardIndex", - "shard_index", "render_shardmap", ] diff --git a/tests/test_viz.py b/tests/test_viz.py index d82f13c..a35c7c2 100644 --- a/tests/test_viz.py +++ b/tests/test_viz.py @@ -1,8 +1,9 @@ -"""Tests for the headless shard-map render core (issue #38, phase 1). +"""Tests for the shard-map viewer (issue #38). -Pure Python -- no ipyleaflet / widget stack. Exercises GeoJSON emission off a -small saved ShardMap fixture, the catalog footprint layer, the viewport -grid-on-zoom gate, and antimeridian splitting. +Pure Python -- no ipyleaflet / widget stack for the render core. Exercises +GeoJSON emission off a small saved ShardMap fixture, the catalog footprint +layer, antimeridian splitting, CRS picking, and the ipyleaflet wrapper +(skipped when the viz extra isn't installed). """ import json @@ -20,7 +21,6 @@ grid_from_signature, render_shardmap, shard_outlines, - viewport_cells, ) from zagg.viz.crs import crs_info, is_polar, pick_crs, shardmap_bbox from zagg.viz.shardmap import _is_geojson, _split_antimeridian @@ -76,22 +76,6 @@ def catalog(): ) -@pytest.fixture -def healpix_shardmap(): - """A HEALPix shardmap (parent_order=6, child_order=12) over a tight AOI. - - Child cells subdivide each order-6 shard 4-for-1 per order step - (``4^(12-6)`` per shard), so the grid-on-zoom must show those nested child - cells -- not the shard outline redrawn. - """ - g = HealpixGrid(parent_order=6, child_order=12, layout="fullsphere") - lats = np.array([10.0, 10.5, 10.5, 10.0, 10.0]) - lons = np.array([20.0, 20.0, 20.5, 20.5, 20.0]) - keys = [int(k) for k in g.coverage([(lats, lons)])] - granules = [[{"id": "G", "s3": "s", "https": "h"}] for _ in keys] - return ShardMap(g.signature(), keys, granules, {"backend": "test"}) - - @pytest.fixture def antarctic_shardmap(): """ShardMap on an EPSG:3031 grid -> footprints entirely south of -60 deg.""" @@ -171,287 +155,6 @@ def test_geometry_is_polygon(self, catalog): assert fc["features"][0]["geometry"]["type"] == "Polygon" -# ── viewport_cells (grid-on-zoom gate) ─────────────────────────────────────── - - -class TestViewportCells: - def test_gate_open_few_shards(self, shardmap): - # Tight viewport over one chunk -> <= max_shards visible, grid drawn. - fp = grid_from_signature(shardmap.grid_signature).shard_footprint(0) - lon0, lat0, lon1, lat1 = fp.bounds - fc = viewport_cells(shardmap, (lon0, lat0, lon1, lat1), max_shards=4) - assert _is_geojson(fc) - assert len(fc["features"]) >= 1 - - def test_gate_closed_too_many_shards(self, shardmap): - # Global viewport over all shards but max_shards=1 -> gated empty. - fc = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=1) - assert fc["features"] == [] - - def test_gate_closed_no_shards(self, shardmap): - # Viewport far from the grid -> nothing visible -> empty. - fc = viewport_cells(shardmap, (10, 10, 11, 11), max_shards=4) - assert fc["features"] == [] - - def test_clipped_to_viewport(self, shardmap): - grid = grid_from_signature(shardmap.grid_signature) - fp = grid.shard_footprint(0) - lon0, lat0, lon1, lat1 = fp.bounds - # Half-width viewport -> clipped cell stays within the viewport bbox. - view = (lon0, lat0, (lon0 + lon1) / 2, lat1) - fc = viewport_cells(shardmap, view, max_shards=4) - for feat in fc["features"]: - for ring in feat["geometry"]["coordinates"]: - for lon, lat in ring: - assert view[0] - 1e-6 <= lon <= view[2] + 1e-6 - - def test_footprints_built_once_across_queries(self, shardmap, monkeypatch): - # Regression for the grid-on-zoom hang (PR #44): footprints must be - # generated once for the STRtree index and never again per query. Wrap - # the grid's shard_footprint with a call counter and run many viewports. - from zagg.grids.rectilinear import RectilinearGrid - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() # fresh index for this shardmap - calls = {"n": 0} - orig = RectilinearGrid.shard_footprint - - def counting(self, key): - calls["n"] += 1 - return orig(self, key) - - monkeypatch.setattr(RectilinearGrid, "shard_footprint", counting) - - full = (-180, -90, 180, 90) - for _ in range(10): - viewport_cells(shardmap, full, max_shards=10) - # Exactly one footprint build per shard for index construction; the ten - # queries add nothing. - assert calls["n"] == len(shardmap.shard_keys) - - def test_prebuilt_index_reused(self, shardmap, monkeypatch): - # Passing an explicit index never touches shard_footprint at query time. - from zagg.grids.rectilinear import RectilinearGrid - from zagg.viz.shardmap import ShardIndex - - index = ShardIndex(shardmap) - - def boom(self, key): # pragma: no cover - must not be called - raise AssertionError("shard_footprint rebuilt during query") - - monkeypatch.setattr(RectilinearGrid, "shard_footprint", boom) - fc = viewport_cells(shardmap, (-180, -90, 180, 90), max_shards=10, index=index) - assert _is_geojson(fc) - - def test_index_query_matches_exhaustive_scan(self, shardmap): - # The STRtree-backed visible set equals the old exhaustive intersects - # scan, for both a tight and a wide bbox. - from shapely.geometry import box - - from zagg.viz.shardmap import grid_from_signature, shard_index - - grid = grid_from_signature(shardmap.grid_signature) - idx = shard_index(shardmap) - for bbox in ((-180, -90, 180, 90), grid.shard_footprint(0).bounds): - view = box(*bbox) - indexed = {k for k, _ in idx.query(view)} - exhaustive = { - k for k in shardmap.shard_keys if grid.shard_footprint(k).intersects(view) - } - assert indexed == exhaustive - - -# ── viewport_cells: HEALPix child-cell nesting + viewport bound ────────────── - - -class TestViewportCellsHealpix: - def _tight_view(self, grid, key, frac=0.3): - """A viewport centered in shard ``key``, ``frac`` of its bbox each side.""" - lon0, lat0, lon1, lat1 = grid.shard_footprint(key).bounds - cx, cy = (lon0 + lon1) / 2, (lat0 + lat1) / 2 - w, h = (lon1 - lon0) * frac / 2, (lat1 - lat0) * frac / 2 - return (cx - w, cy - h, cx + w, cy + h) - - def test_emits_child_order_cells_not_shard_outline(self, healpix_shardmap): - # The grid-on-zoom for HEALPix is the finer child-cell grid, so a single - # visible shard must yield *many* child cells -- not one shard-clip - # feature (the bug @espg reported: grid == shards, no nesting). - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - grid = grid_from_signature(healpix_shardmap.grid_signature) - key = healpix_shardmap.shard_keys[0] - view = self._tight_view(grid, key, frac=0.3) - fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=5000) - assert _is_geojson(fc) - assert len(fc["features"]) > 1 # not a lone shard-outline clip - - def test_cells_nest_within_parent_shard(self, healpix_shardmap): - # Every emitted cell's parent (clip2order at parent_order) is one of the - # visible shards, and the cell IDs are genuine order-child_order cells. - from mortie import clip2order, infer_order_from_morton - - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - grid = grid_from_signature(healpix_shardmap.grid_signature) - key = int(healpix_shardmap.shard_keys[0]) - view = self._tight_view(grid, key, frac=0.2) - fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=5000) - shard_keys = {int(k) for k in healpix_shardmap.shard_keys} - assert fc["features"] - for feat in fc["features"]: - cell = feat["properties"]["cell"] - parent = feat["properties"]["shard_key"] - assert int(infer_order_from_morton(cell)) == grid.child_order - assert int(clip2order(grid.parent_order, np.asarray([cell]))[0]) == parent - assert parent in shard_keys - - def test_cell_union_tiles_visible_shard(self, healpix_shardmap): - # A viewport covering exactly one whole shard: the union of the emitted - # child cells must reconstruct that shard footprint (4-for-1 nesting), - # confirming the grid lines up inside the shard outline. - from shapely.ops import unary_union - - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - # Use a smaller level diff so a whole-shard view stays cheap. - g = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere") - lats = np.array([10.0, 10.5, 10.5, 10.0, 10.0]) - lons = np.array([20.0, 20.0, 20.5, 20.5, 20.0]) - keys = [int(k) for k in g.coverage([(lats, lons)])] - sm = ShardMap(g.signature(), keys, [[] for _ in keys], {}) - key = keys[0] - shard_fp = g.shard_footprint(key) - # The shard's bbox overlaps its 4 neighbors (a HEALPix diamond's bbox is - # larger than the diamond), so allow them through the gate; cells are - # then filtered to this shard via the parent check below. - view = shard_fp.bounds - fc = viewport_cells(sm, view, max_shards=8, max_cells=5000) - cell_polys = [ - g.shard_footprint(f["properties"]["cell"]) - for f in fc["features"] - if f["properties"]["shard_key"] == key - ] - union = unary_union(cell_polys) - # Union of the shard's child cells reconstructs the shard footprint. - assert union.intersection(shard_fp).area / shard_fp.area > 0.98 - - def test_cell_count_scales_with_viewport_not_fan_out(self, healpix_shardmap): - # Zooming in (smaller viewport) must emit *fewer* cells. A naive - # 4^(child-parent) per-shard enumeration would emit a constant 4096/shard - # regardless of zoom; viewport-bounded coverage shrinks with the view. - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - grid = grid_from_signature(healpix_shardmap.grid_signature) - key = healpix_shardmap.shard_keys[0] - wide = viewport_cells( - healpix_shardmap, - self._tight_view(grid, key, 0.4), - max_shards=4, - max_cells=100000, - ) - tight = viewport_cells( - healpix_shardmap, - self._tight_view(grid, key, 0.05), - max_shards=4, - max_cells=100000, - ) - assert len(tight["features"]) < len(wide["features"]) - # And both are far below the full 4^(12-6)=4096 per-shard enumeration. - assert len(wide["features"]) < grid.n_children - - def test_no_full_enumeration_on_zoom(self, healpix_shardmap, monkeypatch): - # A zoomed-in query must not enumerate every child of a shard. Guard - # generate_morton_children: the viewport path uses morton_coverage, so a - # full child enumeration would be a regression. - import mortie - - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - - def boom(*a, **k): # pragma: no cover - must not be called - raise AssertionError("full child enumeration on a zoomed-in query") - - monkeypatch.setattr(mortie, "generate_morton_children", boom) - grid = grid_from_signature(healpix_shardmap.grid_signature) - key = healpix_shardmap.shard_keys[0] - fc = viewport_cells( - healpix_shardmap, - self._tight_view(grid, key, 0.1), - max_shards=4, - max_cells=5000, - ) - assert _is_geojson(fc) - - def test_max_cells_gate_returns_empty(self, healpix_shardmap): - # A whole-shard view at child_order 12 is tens of thousands of cells; - # the max_cells gate keeps a refresh bounded by returning empty rather - # than emitting them all (the dense-viewport bound). - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - grid = grid_from_signature(healpix_shardmap.grid_signature) - key = healpix_shardmap.shard_keys[0] - view = grid.shard_footprint(key).bounds - fc = viewport_cells(healpix_shardmap, view, max_shards=4, max_cells=200) - assert fc["features"] == [] - - def test_child_cells_respect_antimeridian_seam(self): - # Child cells straddling +-180 split into hemisphere-local parts (no - # globe-spanning band) when split_seam=True, and stay single Polygons - # under a polar CRS (split_seam=False) -- same seam handling as shards. - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - g = HealpixGrid(parent_order=4, child_order=8, layout="fullsphere") - key = int( - g.coverage([(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))])[0] - ) - sm = ShardMap(g.signature(), [key], [[]], {}) - view = g.shard_footprint(key).bounds - split = viewport_cells(sm, view, max_shards=8, max_cells=5000, split_seam=True) - assert any(f["geometry"]["type"] == "MultiPolygon" for f in split["features"]) - for feat in split["features"]: - geom = feat["geometry"] - polys = geom["coordinates"] if geom["type"] == "MultiPolygon" else [geom["coordinates"]] - for poly in polys: - lons = [pt[0] for pt in poly[0]] - assert max(lons) - min(lons) <= 180.0 # no globe-spanning band - unsplit = viewport_cells(sm, view, max_shards=8, max_cells=5000, split_seam=False) - assert all(f["geometry"]["type"] == "Polygon" for f in unsplit["features"]) - - def test_seam_cell_clipped_to_true_sliver_not_global_band(self): - # Regression (PR #44 review): a child cell straddling +-180 must be split - # at the seam *before* clipping. The unsplit cell's flat ring spans - # ~360 deg; clipping that band to a seam-hugging viewport yields a wrong, - # oversized geometry (or drops the cell). After the fix, every emitted - # part is a small hemisphere-local sliver -- never a near-global band. - from shapely.geometry import shape - - from zagg.viz import shardmap as sm_mod - - sm_mod._INDEX_CACHE.clear() - g = HealpixGrid(parent_order=4, child_order=9, layout="fullsphere") - key = int( - g.coverage([(np.array([0.0, 1, 1, 0, 0]), np.array([179.0, 179, 180, 180, 179]))])[0] - ) - sm = ShardMap(g.signature(), [key], [[]], {}) - # A narrow viewport hugging the seam at the shard's latitude band. - lon0, lat0, lon1, lat1 = g.shard_footprint(key).bounds - view = (179.0, lat0, 180.0, (lat0 + lat1) / 2) - fc = viewport_cells(sm, view, max_shards=8, max_cells=20000, split_seam=True) - assert fc["features"] # the seam region is not silently emptied - # An order-9 cell footprint is ~0.009 deg^2. Clipping the unsplit flat - # band to the view instead fills the whole view strip (~0.19 deg^2, 20x - # bigger) -- so a small per-cell area bound pins the seam-first fix. - for feat in fc["features"]: - assert shape(feat["geometry"]).area < 0.05 # true cell, not a band - - # ── antimeridian splitting ─────────────────────────────────────────────────── @@ -527,13 +230,11 @@ def test_shards_only(self, shardmap): out = render_shardmap(shardmap) assert _is_geojson(out["shards"]) assert out["granules"] is None - assert out["cells"] is None - def test_with_catalog_and_bbox(self, shardmap, catalog): - out = render_shardmap(shardmap, catalog, bbox=(-180, -90, 180, 90)) + def test_with_catalog(self, shardmap, catalog): + out = render_shardmap(shardmap, catalog) assert _is_geojson(out["shards"]) assert _is_geojson(out["granules"]) - assert _is_geojson(out["cells"]) def test_from_json_path(self, shardmap, tmp_path): path = tmp_path / "sm.json" @@ -548,7 +249,7 @@ def test_from_geoparquet_path(self, shardmap, catalog, tmp_path): assert len(out["granules"]["features"]) == 2 -# ── phase C: CRS selection (headless, no browser) ──────────────────────────── +# ── CRS selection (headless, no browser) ───────────────────────────────────── class TestCrsSelection: @@ -604,7 +305,7 @@ def test_bad_crs_raises(self): crs_info("EPSG:9999") -# ── phase C: CRS-aware antimeridian seam ───────────────────────────────────── +# ── CRS-aware antimeridian seam ────────────────────────────────────────────── class TestSeamAwareLayers: @@ -650,53 +351,7 @@ def test_granule_footprints_seam_flag(self): assert unsplit["type"] == "Polygon" -# ── debounce wrapper (pure stdlib; no widget stack) ────────────────────────── - - -class TestDebounce: - def test_rapid_calls_coalesce_to_one_on_loop(self): - # Importing leaflet pulls only stdlib + zagg.viz at module level (the - # ipyleaflet imports are local to the functions), so this is headless. - # The debounce now schedules on the *running* event loop (the kernel's - # main thread) rather than a background threading.Timer -- that is the - # fix for the comm-thread crash (PR #44). Drive it on a loop and assert - # the burst coalesces to one call, fired on the loop thread. - import asyncio - import threading - - from zagg.viz.leaflet import _debounce - - async def scenario(): - calls = {"n": 0, "tid": None} - - def hit(): - calls["n"] += 1 - calls["tid"] = threading.get_ident() - - deb = _debounce(0.05, hit) - for _ in range(20): - deb() # each cancels and reschedules the pending loop callback - assert calls["n"] == 0 # nothing fired yet (still within the window) - await asyncio.sleep(0.2) - return calls - - result = asyncio.run(scenario()) - assert result["n"] == 1 # the burst coalesced into a single refresh - assert result["tid"] == threading.get_ident() # ran on the loop's thread - - def test_no_loop_runs_synchronously(self): - # With no running event loop (plain script / headless), there is no comm - # to protect and no loop to schedule on, so the call runs inline. - from zagg.viz.leaflet import _debounce - - calls = {"n": 0} - deb = _debounce(0.05, lambda: calls.__setitem__("n", calls["n"] + 1)) - deb() - assert calls["n"] == 1 - deb.cancel() - - -# ── phase 2: ipyleaflet wrapper (skips when the viz extra isn't installed) ──── +# ── ipyleaflet wrapper (skips when the viz extra isn't installed) ───────────── class TestShowShardmap: @@ -716,11 +371,8 @@ def test_build_map(self, shardmap, tmp_path): shardmap.to_json(str(path)) m = show_shardmap(str(path)) assert isinstance(m, Map) - # shard layer + grid layer (+ basemap) present. + # shard layer (+ basemap) present. assert len(m.layers) >= 2 - # Debounced bounds observer is wired with a reachable cancel hook. - assert callable(getattr(m, "cancel_grid_refresh", None)) - m.cancel_grid_refresh() def test_build_map_with_catalog(self, shardmap, catalog, tmp_path): pytest.importorskip("ipyleaflet") From 33eba21cf9f29726274f0653290623586acb4172 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 23 Jun 2026 11:22:16 +0000 Subject: [PATCH 15/17] fold review: drop grid mention in notebook intro (#38) --- notebooks/shardmap_viewer.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index 29641ad..37b69d2 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -19,7 +19,7 @@ " `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines\n", " and granule footprints.\n", "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) \u2014 builds an interactive\n", - " basemap + shard layer + toggleable footprint layer + zoom-gated grid.\n", + " context basemap + shard layer + toggleable granule-footprint layer.\n", "\n", "**Polar-aware projection.** `show_shardmap` picks the display CRS from the\n", "map's extent: a polar AOI renders in NASA polar-stereographic (EPSG:3031\n", From fc50ea45292e7123792be1f7f8c8ec1d0be73942 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 27 Jun 2026 11:13:32 +0000 Subject: [PATCH 16/17] align shardmap viewer notebook with #105 binder env (#38) --- .binder/environment.yml | 9 ++-- .binder/postBuild | 12 +++-- notebooks/shardmap_viewer.ipynb | 82 ++++++++------------------------- 3 files changed, 31 insertions(+), 72 deletions(-) diff --git a/.binder/environment.yml b/.binder/environment.yml index 6c7a638..98c7a6c 100644 --- a/.binder/environment.yml +++ b/.binder/environment.yml @@ -8,10 +8,11 @@ # there is a single source of truth for the notebook dependencies (pyproject.toml). # # Only PyPI/conda-forge resolvable packages appear here: no ``git+`` URLs and no -# spherely fork. The three Binder-runnable notebooks -# (custom_aggregations, rasterized_zarr, jupyterhub_example) read the public, -# anonymous source.coop store or synthetic data, so the exact-S2 spherely -# SpatialIndex backend is never on their import path. +# spherely fork. The Binder-runnable notebooks read synthetic data, the public +# anonymous source.coop store, or anonymous CMR-STAC granule *metadata* +# (jupyterhub_example, shardmap_viewer) -- all of which use the default HEALPix +# ``mortie`` backend, so the exact-S2 spherely SpatialIndex backend is never on +# their import path. name: zagg-binder channels: - conda-forge diff --git a/.binder/postBuild b/.binder/postBuild index 9fb7827..638088d 100755 --- a/.binder/postBuild +++ b/.binder/postBuild @@ -2,11 +2,13 @@ # repo2docker runs this after building the conda environment. # # Install zagg from the checked-out repo together with its ``analysis`` extra -# (the notebook runtime: jupyter/notebook, xarray, cartopy, matplotlib, ...) and +# (the notebook runtime: jupyter/notebook, xarray, cartopy, matplotlib, ...), # the ``catalog`` extra (stac-geoparquet), which the anonymous CMR-STAC -# catalog-build cell in jupyterhub_example needs. Keeping the dependency set in +# catalog-build cells in jupyterhub_example / shardmap_viewer need, and the +# ``viz`` extra (ipyleaflet -- a pure-Python Jupyter widget) for the +# shardmap_viewer notebook's interactive map. Keeping the dependency set in # pyproject.toml's extras means the Binder image and a local -# ``pip install "zagg[analysis,catalog]"`` resolve the same way. +# ``pip install "zagg[analysis,catalog,viz]"`` resolve the same way. # # We deliberately do NOT install: # * the spherely exact-S2 fork (not on PyPI; not needed -- catalog building @@ -23,8 +25,8 @@ set -euo pipefail # irrelevant for running the example notebooks). git fetch --tags --unshallow 2>/dev/null || git fetch --tags 2>/dev/null || true if git describe --tags --abbrev=0 >/dev/null 2>&1; then - python -m pip install --no-cache-dir ".[analysis,catalog]" + python -m pip install --no-cache-dir ".[analysis,catalog,viz]" else HATCH_VCS_PRETEND_VERSION="0.0.0+binder" \ - python -m pip install --no-cache-dir ".[analysis,catalog]" + python -m pip install --no-cache-dir ".[analysis,catalog,viz]" fi diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index 37b69d2..8a99b57 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -4,58 +4,14 @@ "cell_type": "markdown", "id": "cell-0", "metadata": {}, - "source": [ - "# Shard-map viewer\n", - "\n", - "This notebook demonstrates the `zagg.viz` shard-map viewer (issue\n", - "[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\n", - "ATL06 granules** fetched anonymously from NASA CMR-STAC (no Earthdata\n", - "Login required \u2014 anonymous CMR queries are publicly accessible).\n", - "\n", - "The viewer has two layers:\n", - "\n", - "- **Headless render core** (`zagg.viz.render_shardmap` and friends) \u2014 pure\n", - " Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n", - " `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines\n", - " and granule footprints.\n", - "- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) \u2014 builds an interactive\n", - " context basemap + shard layer + toggleable granule-footprint layer.\n", - "\n", - "**Polar-aware projection.** `show_shardmap` picks the display CRS from the\n", - "map's extent: a polar AOI renders in NASA polar-stereographic (EPSG:3031\n", - "Antarctic / EPSG:3413 Arctic) with a **GIBS** basemap \u2014 undistorted at the\n", - "pole \u2014 while mid-latitude AOIs stay on Web Mercator + OpenStreetMap.\n", - "\n", - "## Install\n", - "\n", - "```bash\n", - "pip install \"zagg[viz,catalog]\" # core + ipyleaflet widget + stac-geoparquet\n", - "```\n", - "\n", - "**Requires an internet connection** for the anonymous CMR-STAC query in\n", - "section 1. No credentials are needed.\n", - "\n", - "## Areas of interest\n", - "\n", - "The notebook queries two AOIs:\n", - "\n", - "1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) \u2014 January\n", - " 2020. Dense ICESat-2 coverage near the pole; O(10\u201340) granules in two\n", - " weeks. Renders in EPSG:3031 (Antarctic Polar Stereographic).\n", - "2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) \u2014 June\n", - " 2020. Second example showing the same pipeline on an Arctic AOI\n", - " (EPSG:3413).\n", - "\n", - "Both inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\n", - "day one \u2014 see section 3 for the round-trip to disk and reload." - ] + "source": "# Shard-map viewer\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/englacial/zagg/main?urlpath=lab/tree/notebooks/shardmap_viewer.ipynb)\n\n_Runs end-to-end on [Binder](https://mybinder.org/v2/gh/englacial/zagg/main?urlpath=lab/tree/notebooks/shardmap_viewer.ipynb): it queries real ICESat-2 ATL06 granule **metadata** anonymously from NASA CMR-STAC (no Earthdata Login -- the CMR-STAC search is anonymous) and renders with `ipyleaflet`. The Binder image already provides `zagg[analysis,catalog,viz]` (incl. `stac-geoparquet` and `ipyleaflet`) via the repo's `.binder/` environment._\n\nThis notebook demonstrates the `zagg.viz` shard-map viewer (issue\n[#38](https://github.com/englacial/zagg/issues/38)) with **real ICESat-2\nATL06 granules** fetched anonymously from NASA CMR-STAC.\n\nThe viewer has two layers:\n\n- **Headless render core** (`zagg.viz.render_shardmap` and friends) — pure\n Python, no browser/widget stack. Turns a `ShardMap` (plus an optional\n `Catalog`) into WGS84 GeoJSON `FeatureCollection` dicts: shard outlines\n and granule footprints.\n- **ipyleaflet wrapper** (`zagg.viz.show_shardmap`) — builds an interactive\n context basemap + shard layer + toggleable granule-footprint layer.\n\n**Polar-aware projection.** `show_shardmap` picks the display CRS from the\nmap's extent: a polar AOI renders in NASA polar-stereographic (EPSG:3031\nAntarctic / EPSG:3413 Arctic) with a **GIBS** basemap — undistorted at the\npole — while mid-latitude AOIs stay on Web Mercator + OpenStreetMap.\n\n## Install\n\nOn Binder this is already set up by `.binder/postBuild`. Locally:\n\n```bash\npip install \"zagg[analysis,catalog,viz]\" # core + stac-geoparquet + ipyleaflet widget\n```\n\nThe `viz` extra pulls in `ipyleaflet` (the interactive map); `catalog` pulls\nin `stac-geoparquet` (the CMR-STAC catalog). **Requires an internet\nconnection** for the anonymous CMR-STAC query in section 1 — no credentials\nare needed.\n\n## Areas of interest\n\nThe notebook queries two AOIs:\n\n1. **Antarctic Peninsula** (`[-65, -70, -55, -64]` WGS84 lon/lat) — January\n 2020. Dense ICESat-2 coverage near the pole; O(10–40) granules in two\n weeks. Renders in EPSG:3031 (Antarctic Polar Stereographic).\n2. **Jakobshavn Glacier, West Greenland** (`[-52, 68, -45, 72]`) — June\n 2020. Second example showing the same pipeline on an Arctic AOI\n (EPSG:3413).\n\nBoth inputs (ShardMap JSON and STAC-geoparquet Catalog) are supported from\nday one — see section 3 for the round-trip to disk and reload." }, { "cell_type": "markdown", "id": "cell-1", "metadata": {}, "source": [ - "## 1. Antarctic Peninsula \u2014 fetch real ATL06 granules from NASA CMR-STAC\n", + "## 1. Antarctic Peninsula — fetch real ATL06 granules from NASA CMR-STAC\n", "\n", "`CMRSource` speaks directly to NASA's CMR-STAC endpoint (`requests`).\n", "No Earthdata Login or credentials are needed for anonymous granule-metadata\n", @@ -69,7 +25,7 @@ "id": "cell-2", "metadata": {}, "outputs": [], - "source": "from zagg.catalog.sources import CMRSource, Query\n\n# Antarctic Peninsula AOI \u2014 two weeks in January 2020.\nAOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\nSTART_DATE = \"2020-01-01\"\nEND_DATE = \"2020-01-15\"\n\nquery = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START_DATE,\n end_date=END_DATE,\n region=AOI_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog = CMRSource().fetch(query)\nprint(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" + "source": "from zagg.catalog.sources import CMRSource, Query\n\n# Antarctic Peninsula AOI — two weeks in January 2020.\nAOI_BBOX = (-65.0, -70.0, -55.0, -64.0) # lon_min, lat_min, lon_max, lat_max\nSTART_DATE = \"2020-01-01\"\nEND_DATE = \"2020-01-15\"\n\nquery = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START_DATE,\n end_date=END_DATE,\n region=AOI_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog = CMRSource().fetch(query)\nprint(f\"Fetched {len(catalog)} granules for {query.collection} over {AOI_BBOX}\")" }, { "cell_type": "code", @@ -78,7 +34,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Inspect the catalog schema \u2014 each row is a STAC Item with WKB geometry.\n", + "# Inspect the catalog schema — each row is a STAC Item with WKB geometry.\n", "print(\"Schema columns:\", catalog.table.schema.names[:10], \"...\")\n", "print(\"Total rows:\", catalog.table.num_rows)\n", "\n", @@ -96,14 +52,14 @@ "source": [ "## 2. Build a ShardMap on a HEALPix grid\n", "\n", - "We use a HEALPix grid (`parent_order=6, child_order=12`) \u2014 the same\n", + "We use a HEALPix grid (`parent_order=6, child_order=12`) — the same\n", "configuration as `src/zagg/configs/atl06.yaml`. The `mortie` backend is\n", "used for HEALPix intersection and requires no extra install. If `spherely`\n", "is installed, pass `backend='auto'` to prefer exact S2 intersections.\n", "\n", "`ShardMap.build` maps each grid shard (a parent-order HEALPix cell) to the\n", "set of granules whose footprint intersects it. The manifest is self-contained\n", - "\u2014 it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", + "— it stores `{\"id\", \"s3\", \"https\"}` per granule, so the aggregation runner\n", "never needs the Catalog again at run time." ] }, @@ -141,7 +97,7 @@ "# Inspect a few shard assignments.\n", "for key, shard_granules in zip(shardmap.shard_keys[:4], shardmap.granules[:4]):\n", " ids = [g[\"id\"] for g in shard_granules]\n", - " print(f\" shard {key:6d}: {len(ids)} granule(s) \u2014 {ids[:2]}{'...' if len(ids) > 2 else ''}\")" + " print(f\" shard {key:6d}: {len(ids)} granule(s) — {ids[:2]}{'...' if len(ids) > 2 else ''}\")" ] }, { @@ -194,7 +150,7 @@ "\n", "`render_shardmap` assembles every layer into one dict of GeoJSON\n", "`FeatureCollection`s. Pass a `catalog` to include granule footprints. Each\n", - "value is plain, JSON-serializable GeoJSON \u2014 no widgets required." + "value is plain, JSON-serializable GeoJSON — no widgets required." ] }, { @@ -244,19 +200,19 @@ "id": "cell-15", "metadata": {}, "source": [ - "## 5. Interactive map \u2014 Antarctic Peninsula (ipyleaflet, polar-stereographic)\n", + "## 5. Interactive map — Antarctic Peninsula (ipyleaflet, polar-stereographic)\n", "\n", "`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\n", "Catalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\n", "as file paths or in-memory objects.\n", "\n", "**Projection is auto-selected from the map's extent.** This Antarctic AOI is\n", - "entirely south of \u221260\u00b0, so the viewer renders in **EPSG:3031** (Antarctic\n", + "entirely south of −60°, so the viewer renders in **EPSG:3031** (Antarctic\n", "Polar Stereographic) with a NASA **GIBS** polar basemap instead of the\n", - "distorted Web Mercator default. An Arctic AOI (poleward of +60\u00b0) gets\n", + "distorted Web Mercator default. An Arctic AOI (poleward of +60°) gets\n", "**EPSG:3413**; mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Pass\n", "`crs=\"EPSG:3857\"` to force Mercator, or `crs=\"EPSG:3413\"`/`\"EPSG:3031\"` to\n", - "force a specific pole. Vector layers stay WGS84 GeoJSON \u2014 proj4leaflet\n", + "force a specific pole. Vector layers stay WGS84 GeoJSON — proj4leaflet\n", "reprojects them client-side.\n", "\n", "**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\n", @@ -265,11 +221,11 @@ "\n", "### Verification checklist\n", "\n", - "1. **Basemap** \u2014 the GIBS polar basemap renders, pans, zooms, and the\n", + "1. **Basemap** — the GIBS polar basemap renders, pans, zooms, and the\n", " continent is undistorted at the pole (no Web-Mercator stretching).\n", - "2. **Shard outlines** \u2014 blue polygons over the Peninsula; click one for\n", + "2. **Shard outlines** — blue polygons over the Peninsula; click one for\n", " its `shard_key` and `n_granules` properties.\n", - "3. **Granule footprints toggle** \u2014 layer-control (top-right); the ICESat-2\n", + "3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n", " track footprints appear/disappear." ] }, @@ -279,14 +235,14 @@ "id": "cell-16", "metadata": {}, "outputs": [], - "source": "from zagg.viz import show_shardmap\n\n# File-path interface \u2014 same as after `python -m zagg.catalog ...`.\n# CRS auto-selected from the AOI extent: this Antarctic map renders in EPSG:3031.\nm = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\nprint(\"display CRS:\", m.crs[\"name\"])\nprint(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\nm" + "source": "from zagg.viz import show_shardmap\n\n# File-path interface — same as after `python -m zagg.catalog ...`.\n# CRS auto-selected from the AOI extent: this Antarctic map renders in EPSG:3031.\nm = show_shardmap(str(sm_path), catalog=str(cat_path), zoom=5)\nprint(\"display CRS:\", m.crs[\"name\"])\nprint(\"layers:\", [getattr(layer, \"name\", type(layer).__name__) for layer in m.layers])\nm" }, { "cell_type": "markdown", "id": "cell-17", "metadata": {}, "source": [ - "## 6. Second AOI \u2014 Jakobshavn Glacier, West Greenland\n", + "## 6. Second AOI — Jakobshavn Glacier, West Greenland\n", "\n", "A second example using an Arctic AOI to show the same pipeline working\n", "independently on a different region. Jakobshavn Glacier is one of\n", @@ -299,7 +255,7 @@ "id": "cell-18", "metadata": {}, "outputs": [], - "source": "# Imports also available from earlier cells; repeated here for standalone re-runs.\nfrom zagg.viz import render_shardmap, show_shardmap\n# Second AOI: Jakobshavn Glacier, West Greenland \u2014 June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render \u2014 confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" + "source": "# Imports also available from earlier cells; repeated here for standalone re-runs.\nfrom zagg.viz import render_shardmap, show_shardmap\n# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render — confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" }, { "cell_type": "code", @@ -323,4 +279,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file From c3d2e1b12a671b3f2ce0920fd007970c52350cb4 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 27 Jun 2026 11:16:14 +0000 Subject: [PATCH 17/17] fold review: align section-5 install note to zagg[analysis,catalog,viz] (#38) --- notebooks/shardmap_viewer.ipynb | 32 ++------------------------------ 1 file changed, 2 insertions(+), 30 deletions(-) diff --git a/notebooks/shardmap_viewer.ipynb b/notebooks/shardmap_viewer.ipynb index 8a99b57..60f1866 100644 --- a/notebooks/shardmap_viewer.ipynb +++ b/notebooks/shardmap_viewer.ipynb @@ -199,35 +199,7 @@ "cell_type": "markdown", "id": "cell-15", "metadata": {}, - "source": [ - "## 5. Interactive map — Antarctic Peninsula (ipyleaflet, polar-stereographic)\n", - "\n", - "`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\n", - "Catalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\n", - "as file paths or in-memory objects.\n", - "\n", - "**Projection is auto-selected from the map's extent.** This Antarctic AOI is\n", - "entirely south of −60°, so the viewer renders in **EPSG:3031** (Antarctic\n", - "Polar Stereographic) with a NASA **GIBS** polar basemap instead of the\n", - "distorted Web Mercator default. An Arctic AOI (poleward of +60°) gets\n", - "**EPSG:3413**; mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Pass\n", - "`crs=\"EPSG:3857\"` to force Mercator, or `crs=\"EPSG:3413\"`/`\"EPSG:3031\"` to\n", - "force a specific pole. Vector layers stay WGS84 GeoJSON — proj4leaflet\n", - "reprojects them client-side.\n", - "\n", - "**Run in JupyterLab** with `pip install \"zagg[viz,catalog]\"` to see the live map.\n", - "Under headless `nbconvert` the Map object is constructed (no error) but\n", - "tiles won't display.\n", - "\n", - "### Verification checklist\n", - "\n", - "1. **Basemap** — the GIBS polar basemap renders, pans, zooms, and the\n", - " continent is undistorted at the pole (no Web-Mercator stretching).\n", - "2. **Shard outlines** — blue polygons over the Peninsula; click one for\n", - " its `shard_key` and `n_granules` properties.\n", - "3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n", - " track footprints appear/disappear." - ] + "source": "## 5. Interactive map — Antarctic Peninsula (ipyleaflet, polar-stereographic)\n\n`show_shardmap` builds an `ipyleaflet.Map` from the saved ShardMap and\nCatalog paths. Both ShardMap JSON and STAC-geoparquet Catalog are accepted\nas file paths or in-memory objects.\n\n**Projection is auto-selected from the map's extent.** This Antarctic AOI is\nentirely south of −60°, so the viewer renders in **EPSG:3031** (Antarctic\nPolar Stereographic) with a NASA **GIBS** polar basemap instead of the\ndistorted Web Mercator default. An Arctic AOI (poleward of +60°) gets\n**EPSG:3413**; mid-latitude AOIs stay on Web Mercator + OpenStreetMap. Pass\n`crs=\"EPSG:3857\"` to force Mercator, or `crs=\"EPSG:3413\"`/`\"EPSG:3031\"` to\nforce a specific pole. Vector layers stay WGS84 GeoJSON — proj4leaflet\nreprojects them client-side.\n\n**Run in JupyterLab** (Binder already has the deps via `.binder/postBuild`;\nlocally `pip install \"zagg[analysis,catalog,viz]\"`) to see the live map.\nUnder headless `nbconvert` the Map object is constructed (no error) but\ntiles won't display.\n\n### Verification checklist\n\n1. **Basemap** — the GIBS polar basemap renders, pans, zooms, and the\n continent is undistorted at the pole (no Web-Mercator stretching).\n2. **Shard outlines** — blue polygons over the Peninsula; click one for\n its `shard_key` and `n_granules` properties.\n3. **Granule footprints toggle** — layer-control (top-right); the ICESat-2\n track footprints appear/disappear." }, { "cell_type": "code", @@ -255,7 +227,7 @@ "id": "cell-18", "metadata": {}, "outputs": [], - "source": "# Imports also available from earlier cells; repeated here for standalone re-runs.\nfrom zagg.viz import render_shardmap, show_shardmap\n# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render — confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" + "source": "# render_shardmap / show_shardmap re-imported here for a standalone re-run of\n# this section; CMRSource / Query / ShardMap / grid still come from sections 1-2.\nfrom zagg.viz import render_shardmap, show_shardmap\n\n# Second AOI: Jakobshavn Glacier, West Greenland — June 2020.\nAOI2_BBOX = (-52.0, 68.0, -45.0, 72.0)\nSTART2, END2 = \"2020-06-01\", \"2020-06-15\"\n\nquery2 = Query(\n short_name=\"ATL06\",\n version=\"007\",\n start_date=START2,\n end_date=END2,\n region=AOI2_BBOX,\n provider=\"NSIDC_CPRD\",\n)\n\ncatalog2 = CMRSource().fetch(query2)\nshardmap2 = ShardMap.build(catalog2, grid, backend=\"mortie\")\n\nprint(f\"Greenland AOI: {len(catalog2)} granules, \"\n f\"{len(shardmap2.shard_keys)} shards\")\n\n# Headless render — confirm layers are populated.\nlayers2 = render_shardmap(shardmap2, catalog2)\nprint(\"Layer feature counts:\",\n {name: (fc and len(fc[\"features\"])) for name, fc in layers2.items()})" }, { "cell_type": "code",