enhancement(roadmap): GIS gaps surfaced by Digital-Earth mapLibre-parity audit (conversions / analysis / catalogs)

[MAfarrag]

Context

Who's filing & why. Digital-Earth is a geospatial
visualization package built on top of pyramids: pyramids is the GIS engine (I/O, formats, reprojection,
analysis, tiling, cloud); Digital-Earth only turns pyramids objects into maps/figures. By design Digital-Earth
must not reimplement GIS or import competitor libraries (rasterio / xarray / fiona / GDAL-directly /
earthengine). So when Digital-Earth needs a GIS capability pyramids does not expose, it is filed here.

How this list was produced. While building Digital-Earth's map/web visualization tiers we audited pyramids
0.34.0 against the modern web-GIS feature surface. ~90% was already present (formats incl.
KML/GPX/FlatGeobuf/PMTiles-read, COG/XYZ tiling + overviews, reproject/resample/crop, the full vector op set,
slope/aspect/hillshade/contour/zonal, S3/GCS/Azure/STAC/WMS/WFS). The items below are the remaining gaps. The
PY-* tags are Digital-Earth's internal IDs — ignore them.

How to work this issue (please read before implementing)

• One item = one PR. This is an umbrella; do not implement all of it at once. Start with item 1
  (terrain-RGB) — it is the only item that blocks already-shipped Digital-Earth code; the rest are future-facing.
• Follow pyramids' existing conventions (verified against the 0.34.0 source):
  • Raster behaviour is implemented in an engine mixin under src/pyramids/dataset/engines/ and exposed via a
    thin delegator on the Dataset facade (src/pyramids/dataset/dataset.py). Tiling/I/O lives in
    engines/io.py (to_xyz, get_tile, create_overviews) and engines/cog.py (to_cog); cell/grid in
    engines/cell.py (get_cell_polygons); reprojection in engines/spatial.py (to_crs); vectorize in
    engines/vectorize.py. Vector behaviour lives on FeatureCollection
    (src/pyramids/feature/collection.py: to_parquet, to_file, …).
  • Writers take path first, then keyword-only options after *, and return Path (see to_cog). Reads
    return the pyramids object. Reuse the existing kwarg vocabulary (resampling, out_dtype, nodata, bands).
  • Add tests (the repo's pytest style) and a Google-style docstring with a runnable doctest matching the
    existing methods (e.g. to_xyz/to_cog use a Dataset.create_from_array(...) example).
• Do not assume an external dependency is acceptable — where one is needed (routing, H3, EE) it is called out
  per item; confirm the dependency policy in a comment before adding it.

---

Item 1 — Dataset.to_terrain_rgb(...) · PY-IO.9 · PRIORITY (blocks shipped DE code)

What. Encode a single-band elevation raster (a DEM, metres) into terrain-RGB raster tiles — elevation
packed into the R/G/B channels — so browser/GPU engines (MapLibre, deck.gl, Cesium) can render 3-D terrain.

Where it lives. New method in a raster I/O engine (alongside to_xyz/create_overviews in
src/pyramids/dataset/engines/io.py), delegated from Dataset (dataset.py) like the other to_* writers.

Signature (mirror to_cog's keyword-only style):

def to_terrain_rgb(
    self,
    path: str | Path,
    *,
    encoding: str = "mapbox",      # "mapbox" (Terrain-RGB) or "terrarium" (Mapzen)
    tiles: bool = True,            # True -> write an XYZ {z}/{x}/{y}.png pyramid; False -> one RGB(A) raster
    min_zoom: int = 0,
    max_zoom: int | None = None,   # default: derived from the source resolution
    tile_size: int = 256,
    base_val: float = -10000.0,    # mapbox only: elevation that maps to RGB (0,0,0)
    interval: float = 0.1,         # mapbox only: metres per encoded unit
    resampling: str = "bilinear",  # reproject/downsample resampling
    band: int = 0,
) -> Path:

Behaviour (be exact — do not infer):

1. Take band band of the source as float elevation in metres. Reproject to EPSG:3857 if not already
   (reuse engines/spatial.py::to_crs).
2. Encode each pixel to 8-bit R, G, B:
   • encoding="mapbox" (Mapbox Terrain-RGB): with v = round((height - base_val) / interval):
     R = (v >> 16) & 255, G = (v >> 8) & 255, B = v & 255.
     Decoder (for the round-trip test): height = base_val + (R*65536 + G*256 + B) * interval.
   • encoding="terrarium" (Mapzen): with v = height + 32768:
     R = floor(v / 256), G = floor(v) % 256, B = floor((v - floor(v)) * 256).
     Decoder: height = (R*256 + G + B/256) - 32768.
3. NoData pixels: emit a fully-transparent pixel — write RGBA with A=0 there (A=255 elsewhere). If the
   source has no nodata, write RGB.
4. Output:
   • tiles=True -> an XYZ pyramid of tile_sizextile_size PNGs under path/{z}/{x}/{y}.png, zooms
     min_zoom..max_zoom (reuse the tiling path used by to_xyz/get_tile); default max_zoom from the native
     pixel size.
   • tiles=False -> a single 3- or 4-band PNG/GeoTIFF at path.
5. Return the output Path (the tile-root dir, or the file).

Edge cases to handle (named, not assumed): values outside the encodable range must be clamped (mapbox:
v clamped to [0, 2**24 - 1]); max_zoom < min_zoom -> ValueError; a multi-band source without an explicit
band -> use band 0 and document it; a non-projected source is reprojected (do not error).

Acceptance criteria:

• Round-trip: decoding a written tile/raster with the formula above recovers the source elevation within
  interval (mapbox) / 1/256 m (terrarium).
• tiles=True writes a valid {z}/{x}/{y}.png XYZ layout; tiles load in a MapLibre raster-dem source.
• NoData becomes transparent (alpha 0); a fully-valid DEM yields opaque RGB.
• to_terrain_rgb is exposed on Dataset, has a Google-style docstring + runnable doctest, and tests.

References: Mapbox Terrain-RGB v1 spec; Mapzen Terrarium spec; rio-rgbify is a reference encoder (do not
depend on it — implement with numpy + the existing tiler).

---

Item 2 — FeatureCollection.interpolate_to_raster(...) · PY-AN.15

What. Interpolate scattered point measurements into a continuous raster surface (IDW + ordinary Kriging).
FeatureCollection.interpolate today is line/1-D interpolation — this is point->grid.

Where. New method on FeatureCollection (src/pyramids/feature/collection.py); returns a Dataset.

Signature:

def interpolate_to_raster(
    self, column: str, *, method: str = "idw", cell_size: float | None = None,
    bounds: tuple[float, float, float, float] | None = None, power: float = 2.0,
    n_neighbors: int | None = None, variogram: str = "spherical", nodata: float = -9999.0,
) -> "Dataset":

Behaviour: read column as the z-value at each point geometry; build a grid (from bounds+cell_size, else
the layer extent and a sensible default cell size); method="idw" -> inverse-distance weighting with exponent
power (optionally limited to n_neighbors); method="kriging" -> ordinary kriging with the named variogram.
Output a single-band Dataset in the layer's CRS with nodata outside the convex hull.

Edge cases: <3 points -> ValueError; non-numeric / all-NaN column -> ValueError; duplicate coordinates
-> average them. Acceptance: IDW evaluated at a sample point ~= that point's value; output CRS == input CRS;
the grid honours bounds/cell_size. Reference: classic IDW/kriging formulations; scipy / pykrige as
references — confirm dependency policy before adding pykrige.

---

Item 3 — Network analysis (isochrones / service areas / OD matrices) · PY-AN.17

What. Routing-based analysis. Needs a routing engine (OSRM / Valhalla / pgRouting / OSMnx) — so the first
deliverable is a design decision, not code: does pyramids wrap a routing backend, or is this out of scope?
Please comment with a decision before implementing. If in scope, suggested surface:
isochrones(points, *, minutes, profile="driving", backend=...) -> FeatureCollection (polygons) and
od_matrix(origins, destinations, *, profile) -> DataFrame. Acceptance to be defined once the backend is
chosen. Reference: OSRM /table + isochrone, Valhalla isochrone API, OSMnx.

---

Item 4 — H3 hexagonal grids · PY-AN.18

What. Uber H3 indexing/binning. Where. FeatureCollection. Signatures:

def to_h3(self, resolution: int) -> "FeatureCollection"            # add an `h3` index column per point
def h3_bin(self, resolution: int, *, agg: str = "count", column: str | None = None) -> "FeatureCollection"

Behaviour: to_h3 adds the H3 cell index (resolution 0-15) for each point geometry. h3_bin groups points by
cell and returns one hexagon polygon per occupied cell with the aggregate (count, or agg over column).
Output CRS EPSG:4326 (H3 is lat/lng). Edge cases: non-point geometry -> ValueError; resolution outside
0-15 -> ValueError. Acceptance: cell count matches the h3 library for a known fixture; hex polygons are
valid and closed. Reference: h3-py (this needs the h3 dependency — confirm policy first).

---

Item 5 — Spatial statistics (autocorrelation / hotspots) · PY-AN.19

What. Global/local Moran's I, Getis-Ord Gi* hotspots, LISA clusters over a FeatureCollection column. Needs a
spatial-weights concept (contiguity / k-nearest / distance-band) — define that first. Suggested surface:
spatial_autocorrelation(column, *, weights="queen") -> dict (global I + p-value) and
hotspots(column, *, weights="queen") -> FeatureCollection (per-feature Gi* z-score + significance class).
Acceptance: results match esda / libpysal on a known fixture within tolerance. Reference: PySAL esda,
libpysal (confirm dependency policy). Scope precisely before building — this is the largest item (L).

---

Item 6 — Google Earth Engine reader · PY-CL.7

What. Read Earth Engine assets/collections into pyramids Dataset / DatasetCollection. Must live in
pyramids (Digital-Earth forbids xee / earthengine-api). Suggested surface: a pyramids.ee module —
pyramids.ee.read(asset_id, *, region, scale, bands=None) -> Dataset and a collection reader for image
collections. Requires earthengine-api (+ auth) — confirm dependency/auth policy before implementing.
Acceptance: a small public EE image reads into a Dataset with correct CRS/transform/bands. Reference:
earthengine-api, xee (as behaviour reference only).

---

Convenience wrappers (the capability already exists in GDAL/pyramids — just surface it)

Each is small (~half a day). Add method + test + docstring, following the conventions above.

• PY-IO.8 — FeatureCollection.to_pmtiles(path) -> Path. GDAL's PMTiles driver already writes; wrap it like
  the existing to_file / to_parquet. Acceptance: output reopens via read_file and round-trips geometry +
  attributes.
• PY-IO.13 — read_gpx_layers(path) -> dict[str, FeatureCollection]. The GPX driver already exposes
  waypoints / tracks / routes sub-layers; return all present ones keyed by name. Acceptance: a
  multi-layer GPX yields the present sub-layers; missing ones are omitted.
• PY-IO.14 — vector-tile (MVT) writer. Overlaps PY-IO.8; consider folding into that PR. Needed only for
  very large served vector layers.
• PY-AN.8 — FeatureCollection.fishnet(bounds, cell_size, *, crs=None) -> FeatureCollection. Arbitrary-extent
  vector grid of square polygons. Dataset.get_cell_polygons (engines/cell.py) already does the raster-aligned
  case — this is the vector / arbitrary-bbox analogue. Acceptance: cell count == ceil(w/cs) * ceil(h/cs); all
  polygons valid and non-overlapping.
• PY-CL.4 — ArcGIS FeatureServer helper. The ESRIJSON driver already reads FeatureServer output; add a URL
  convenience + result pagination (FeatureServer caps records per request). Acceptance: a paged layer reads all
  features, not just the first page.
• PY-CL.6 — Planetary Computer signing helper. STAC access works (pyramids.stac); confirm or add a
  planetary-computer-style token / URL-signing helper so signed asset hrefs read via the existing readers.

---

Suggested order

1. Item 1 (terrain-RGB) — unblocks shipped Digital-Earth code. Do this first.
2. The convenience wrappers (PY-IO.8, .13, .14, PY-AN.8, PY-CL.4, .6) — quick wins.
3. Items 2, 4, 5 (interpolation, H3, spatial stats) — each its own PR; confirm any new dependency first.
4. Items 3 and 6 (routing, EE) — decision-first: comment with the dependency/scope decision before coding.

[MAfarrag]

FeatureCollection

Scoping the feature/ subpackage items from this umbrella against docs/SCOPE.md (the module is the
fiona/geopandas niche — fully in scope). Verified against the current source: none of these exist on
FeatureCollection yet, and the generic point→grid interpolation engine (Dataset.from_points /
ops/interpolate.py) plus the PMTiles/MVT/GPX/ESRIJSON GDAL drivers are already present, so most of this is thin
wrapping rather than new machinery.

Will implement — in scope, no new dependency (one PR each)

• fishnet(bounds, cell_size, *, crs=None) -> FeatureCollection (PY-AN.8) — arbitrary-extent vector grid of
  square polygons; the vector/bbox analogue of Dataset.get_cell_polygons. Same family as the new
  voronoi/quadtree.
• interpolate_to_raster(column, *, method="idw", cell_size=None, bounds=None, power=2.0, n_neighbors=None, nodata=-9999.0) -> Dataset (Item 2, IDW only) — thin facade over the existing Dataset.from_points
  (gdal.Grid invdist/invdistnn). Point-only; <3 points / non-numeric column → ValueError; coincident
  points averaged.
• to_pmtiles(path, *, min_zoom=0, max_zoom=None, layer_name=None) -> Path (PY-IO.8) — OGR PMTiles driver,
  mirrors to_file.
• to_mvt(path, *, …) -> Path (PY-IO.14) — OGR MVT driver ({z}/{x}/{y}.pbf); shares a private tiled-vector
  helper with to_pmtiles (folded into the same PR).
• read_gpx_layers(path) -> dict[str, FeatureCollection] (PY-IO.13) — returns the present
  waypoints/routes/tracks/… sub-layers via list_layers + read_file(layer=…).
• from_featureserver(url, *, where="1=1", max_records=None, page_size=1000) -> FeatureCollection (PY-CL.4)
  — generic ArcGIS REST read over the ESRIJSON driver with resultOffset pagination (provider-agnostic, not
  EO-agency auth).

All follow the repo conventions: writers take path first, keyword-only opts after *, return Path; readers
return the pyramids object; Google-style docstring + runnable doctest + pytest tests.

Will implement only after a dependency decision (borderline)

• to_h3(resolution) / h3_bin(resolution, *, agg="count", column=None) (Item 4) — needs h3 (optional
  extra [h3]). The binning itself is the same generic point-aggregation as quadtree; only the H3 standard +
  dependency is specialized. Confirm the h3 dependency before I start.
• kriging branch of interpolate_to_raster (method="kriging", Item 2 cont.) — needs pykrige (optional
  extra [kriging]). Confirm pykrige before I start. IDW ships without it.

Both would go through the optional-extras helper pattern (base/_utils.py), as extras, never core.

Out of scope — not taking these here

• Network analysis / routing (Item 3) — distinct domain needing an external routing backend; downstream, not
  the fiona/geopandas niche.
• Spatial statistics (Item 5, Moran's I / Getis-Ord / LISA) — the PySAL/esda domain; belongs downstream
  (pyramids for I/O, PySAL for stats).
• Planetary Computer signing (PY-CL.6) — provider-specific EO signing; docs/SCOPE.md S3 already
  relocated PlanetaryComputerSigner to earthlens.stac. The generic Signer protocol + BearerTokenSigner
  stay.

Suggested order

fishnet → format wrappers (to_pmtiles+to_mvt, then read_gpx_layers) → IDW interpolate_to_raster →
from_featureserver → (after dependency approval) H3 → kriging. One item = one PR.

[MAfarrag]

Item 1 — Dataset.to_terrain_rgb (PY-IO.9) is implemented.

Tracked in its own issue #593 and addressed by PR #590.

• New IO.to_terrain_rgb in src/pyramids/dataset/engines/io.py, delegated from the Dataset facade.
• Mapbox Terrain-RGB and Mapzen Terrarium encodings (numpy bit-packing, out-of-range values clamped); reprojects
  to EPSG:3857; no-data → transparent RGBA, otherwise plain RGB. No new dependency.
• tiles=False writes one RGB(A) raster (PNG by suffix, else GeoTIFF); tiles=True writes a {z}/{x}/{y}.png
  XYZ pyramid, with max_zoom derived from the source resolution when unset.
• Covered by tests/dataset/test_terrain_rgb.py (round-trips, no-data transparency, reprojection, tile layout,
  validation errors) with a Google-style docstring + doctest.

The remaining items (2–6 and the convenience wrappers) are unaffected and still open here.