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/docs/quickstart.md b/docs/quickstart.md index 7b06385..e9190ff 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..60f1866 --- /dev/null +++ b/notebooks/shardmap_viewer.ipynb @@ -0,0 +1,254 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "cell-0", + "metadata": {}, + "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 — 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": "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", + "execution_count": null, + "id": "cell-3", + "metadata": {}, + "outputs": [], + "source": [ + "# 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: {(rec['https'] or 'None')[: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.grids import HealpixGrid\n", + "\n", + "# 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": [ + "# 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 ''}\")" + ] + }, + { + "cell_type": "markdown", + "id": "cell-7", + "metadata": {}, + "source": [ + "## 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. Each\n", + "value is plain, JSON-serializable GeoJSON — no widgets required." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cell-10", + "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()}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cell-11", + "metadata": {}, + "outputs": [], + "source": [ + "# 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", + "\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 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-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** (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", + "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": "cell-17", + "metadata": {}, + "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": "# 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", + "execution_count": null, + "id": "cell-19", + "metadata": {}, + "outputs": [], + "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": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.12.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 41e2990..4b4cc9f 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/__init__.py b/src/zagg/viz/__init__.py new file mode 100644 index 0000000..c1f67b7 --- /dev/null +++ b/src/zagg/viz/__init__.py @@ -0,0 +1,64 @@ +"""Shard-map viewer (issue #38). + +Two layers: + +- :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**. 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 ( + granule_footprints, + grid_from_signature, + render_shardmap, + shard_outlines, +) + + +def show_shardmap(shardmap_path, catalog=None, **kwargs): + """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. + 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", + "grid_from_signature", + "show_shardmap", +] diff --git a/src/zagg/viz/crs.py b/src/zagg/viz/crs.py new file mode 100644 index 0000000..f86ca24 --- /dev/null +++ b/src/zagg/viz/crs.py @@ -0,0 +1,186 @@ +"""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. 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 " + "+ellps=WGS84 +datum=WGS84 +units=m +no_defs" + ), + "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 " + "+datum=WGS84 +units=m +no_defs" + ), + "origin": _GIBS_ORIGIN, + "bounds": _GIBS_BOUNDS, + "resolutions": _GIBS_RESOLUTIONS, +} + +# 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. + + 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 + 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 new file mode 100644 index 0000000..3e2c021 --- /dev/null +++ b/src/zagg/viz/leaflet.py @@ -0,0 +1,151 @@ +"""ipyleaflet wrapper for the shard-map viewer (issue #38). + +Builds an interactive map from a saved ShardMap: a basemap, the shard-outline +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** +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]``. +""" + +from __future__ import annotations + +from zagg.viz.crs import crs_info, is_polar, pick_crs +from zagg.viz.shardmap import ( + _load_catalog, + _load_shardmap, + granule_footprints, + shard_outlines, +) + +# 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} + + +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 _leaflet_crs(projection: dict): + """A proj4leaflet ``ipyleaflet.projections`` CRS dict from a projection def. + + 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"], + "custom": True, + "proj4def": projection["proj4def"], + "origin": projection["origin"], + "bounds": projection["bounds"], + "resolutions": projection["resolutions"], + } + + +def show_shardmap( + shardmap_path, + catalog=None, + *, + 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 + 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. + zoom : int + Initial map zoom. + basemap : ipyleaflet basemap, optional + 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 + ------- + ipyleaflet.Map + 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 + + shardmap = _load_shardmap(shardmap_path) + + 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), split_seam=split_seam) + footprint_layer = GeoJSON( + data=footprint_fc, style=_FOOTPRINT_STYLE, name="granule footprints" + ) + m.add(footprint_layer) + + m.add(LayersControl(position="topright")) + return m + + +__all__ = ["show_shardmap"] diff --git a/src/zagg/viz/shardmap.py b/src/zagg/viz/shardmap.py new file mode 100644 index 0000000..7419454 --- /dev/null +++ b/src/zagg/viz/shardmap.py @@ -0,0 +1,351 @@ +"""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). + +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 + ``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``). + +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`` off the grid, which is 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_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.""" + 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`` (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 + + 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). 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) + + 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 + parts.append(_polygon_rings(sub, shift=shift)) + + if not 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 _plain_geometry(geom) -> dict: + """GeoJSON geometry for a Polygon/MultiPolygon with no antimeridian split. + + 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: + 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 _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, *, 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 + 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. + split_seam : bool + Split polygons crossing +-180 into a ``MultiPolygon``. Set False under a + polar-stereographic CRS, where there is no antimeridian seam. + + 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, split_seam=split_seam), + {"shard_key": _jsonable(key), "n_granules": len(granules)}, + ) + ) + return _collection(features) + + +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 + exterior ring), with the granule id under ``properties``. + + Parameters + ---------- + 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 + ------- + 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), split_seam=split_seam), + {"id": rec["id"]}, + ) + ) + return _collection(features) + + +# ── top-level assembly ─────────────────────────────────────────────────────── + + +def render_shardmap(shardmap, catalog=None) -> dict: + """Assemble the 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. + + Returns + ------- + dict + ``{"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} + if catalog is not None: + out["granules"] = granule_footprints(_load_catalog(catalog)) + 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", + "render_shardmap", +] diff --git a/tests/test_viz.py b/tests/test_viz.py new file mode 100644 index 0000000..a35c7c2 --- /dev/null +++ b/tests/test_viz.py @@ -0,0 +1,423 @@ +"""Tests for the shard-map viewer (issue #38). + +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 + +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, +) +from zagg.viz.crs import crs_info, is_polar, pick_crs, shardmap_bbox +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"}, + ) + + +@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: + 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" + + +# ── 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" + + 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 ───────────────────────────────────────────────── + + +class TestRenderShardmap: + def test_shards_only(self, shardmap): + out = render_shardmap(shardmap) + assert _is_geojson(out["shards"]) + assert out["granules"] is None + + def test_with_catalog(self, shardmap, catalog): + out = render_shardmap(shardmap, catalog) + assert _is_geojson(out["shards"]) + assert _is_geojson(out["granules"]) + + 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 + + +# ── 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") + + +# ── 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): + # 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" + + +# ── 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 (+ 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 + + 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"