MapLibre API

ras2cng.maplibre

Build a compact MapLibre project bundle from a ras2cng archive.

The archive remains the authoritative, queryable delivery format. This module creates a browser delivery companion: one PMTiles archive for geometry and, when requested, one PMTiles archive for raw HDF result values joined to their source model elements. It intentionally does not rasterize results; stored-map COGs are a separate, explicit publication step.

PackageSummary dataclass

Files and counts produced by :func:package_maplibre_viewer.

Source code in ras2cng/maplibre.py

@dataclass(frozen=True)
class PackageSummary:
    """Files and counts produced by :func:`package_maplibre_viewer`."""

    manifest_path: Path
    geometry_pmtiles: Path
    result_pmtiles: Path | None
    geometry_layer_count: int
    result_layer_count: int
    bounds: tuple[float, float, float, float]

TerrainPackageSummary dataclass

Browser terrain artifact produced from an archived terrain COG.

Source code in ras2cng/maplibre.py

@dataclass(frozen=True)
class TerrainPackageSummary:
    """Browser terrain artifact produced from an archived terrain COG."""

    manifest_path: Path
    pmtiles_path: Path
    source_cog: Path
    raster_stats: dict[str, float]
    max_zoom: int

RasterPackageSummary dataclass

Browser display derivative paired with an authoritative numeric COG.

Source code in ras2cng/maplibre.py

@dataclass(frozen=True)
class RasterPackageSummary:
    """Browser display derivative paired with an authoritative numeric COG."""

    manifest_path: Path
    pmtiles_path: Path
    source_cog: Path
    raster_stats: dict[str, float]
    max_zoom: int
    layer_id: str

VectorResultPackageSummary dataclass

Queryable vector Stored Map packaged for browser delivery.

Source code in ras2cng/maplibre.py

@dataclass(frozen=True)
class VectorResultPackageSummary:
    """Queryable vector Stored Map packaged for browser delivery."""

    manifest_path: Path
    pmtiles_path: Path
    source_vector: Path
    feature_count: int
    layer_id: str

package_maplibre_terrain(cog_path, viewer_dir, *, name='Terrain', layer_id=None, source_cog=None, units='ft', visible=True, max_zoom=None, scratch_dir=None, overwrite=False)

Add a RAS-styled, queryable terrain PMTiles layer to a viewer bundle.

The source COG remains the numerical source for identify queries. The PMTiles overlay is a colorized Web Mercator representation used only for display. Its highest zoom is capped at source cell resolution so no terrain detail is invented in the browser.

Source code in ras2cng/maplibre.py

def package_maplibre_terrain(
    cog_path: Path,
    viewer_dir: Path,
    *,
    name: str = "Terrain",
    layer_id: str | None = None,
    source_cog: str | None = None,
    units: str = "ft",
    visible: bool = True,
    max_zoom: int | None = None,
    scratch_dir: Path | None = None,
    overwrite: bool = False,
) -> TerrainPackageSummary:
    """Add a RAS-styled, queryable terrain PMTiles layer to a viewer bundle.

    The source COG remains the numerical source for identify queries. The
    PMTiles overlay is a colorized Web Mercator representation used only for
    display. Its highest zoom is capped at source cell resolution so no
    terrain detail is invented in the browser.
    """

    cog_path = Path(cog_path)
    viewer_dir = Path(viewer_dir)
    manifest_path = viewer_dir / "manifest.json"
    if not cog_path.is_file():
        raise FileNotFoundError(f"Terrain COG does not exist: {cog_path}")
    if not manifest_path.is_file():
        raise FileNotFoundError(f"MapLibre viewer manifest does not exist: {manifest_path}")

    manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
    tilesets = manifest.setdefault("tilesets", [])
    terrain_id = _slug(layer_id or "terrain")
    if not terrain_id:
        raise ValueError("Terrain layer ID must contain at least one letter or number.")
    existing = next((item for item in tilesets if item.get("id") == terrain_id), None)
    if existing and not overwrite:
        raise FileExistsError(
            f"Viewer already has terrain tileset {terrain_id!r}; "
            "pass overwrite=True to replace it."
        )

    for executable in (
        _gdalinfo_command(),
        _gdaldem_command(),
        _gdalwarp_command(),
        _gdal_translate_command(),
        _gdaladdo_command(),
        _pmtiles_command(),
    ):
        _require_cli(executable)

    if scratch_dir is not None:
        scratch_dir = Path(scratch_dir).resolve()
        scratch_dir.mkdir(parents=True, exist_ok=True)
        if not scratch_dir.is_dir():
            raise ValueError(f"Terrain scratch directory is not a directory: {scratch_dir}")

    stats = _raster_stats(_gdalinfo(cog_path))
    raster_metadata = _raster_source_metadata(cog_path)
    tiles_dir = viewer_dir / "tiles"
    tiles_dir.mkdir(parents=True, exist_ok=True)
    output = tiles_dir / f"{terrain_id}.pmtiles"
    if output.exists() and not overwrite:
        raise FileExistsError(f"Terrain PMTiles already exists: {output}")
    source_href = source_cog or _relative_href(cog_path, viewer_dir)

    selected_max_zoom = _render_raster_pmtiles(
        cog_path,
        output,
        ramp_writer=lambda path: _terrain_color_ramp(stats, path),
        max_zoom=max_zoom,
        scratch_dir=scratch_dir,
        prefix=terrain_id,
    )

    terrain_tileset = {
        "id": terrain_id,
        "name": name,
        "type": "raster",
        "href": f"tiles/{terrain_id}.pmtiles",
        "sourceCog": source_href,
        "bytes": output.stat().st_size,
        "tileSize": 256,
        "groupId": "ras-terrains",
        "visible": visible,
        "opacity": 1.0,
        "maxzoom": selected_max_zoom,
        "rasterStats": stats,
        "ramp": "stretched",
        "domainPolicy": "fixed",
        "sourceKind": "terrain",
        "legend": {
            "type": "continuous",
            "mode": "stretched",
            "preset": "rasmapper.terrain",
            "domainPolicy": "fixed",
            "colors": [
                f"#{red:02x}{green:02x}{blue:02x}"
                for _, red, green, blue, _ in _RAS_TERRAIN_COLORS
            ],
        },
        "queryable": True,
        "units": units,
        "storedMap": {
            "mapType": "terrain",
            "source": "HEC-RAS terrain GeoTIFF",
            "cogBytes": cog_path.stat().st_size,
        },
        **raster_metadata,
    }
    if existing:
        tilesets[tilesets.index(existing)] = terrain_tileset
    else:
        tilesets.append(terrain_tileset)
    groups = manifest.setdefault("groups", [])
    if not any(group.get("id") == "ras-terrains" for group in groups):
        groups.append({"id": "ras-terrains", "name": "Terrain", "visible": True})
    apply_manifest_v2(manifest, archive=_viewer_archive_metadata(viewer_dir))
    manifest_path.write_text(json.dumps(manifest, indent=2) + "\n", encoding="utf-8")
    return TerrainPackageSummary(
        manifest_path=manifest_path,
        pmtiles_path=output,
        source_cog=cog_path,
        raster_stats=stats,
        max_zoom=selected_max_zoom,
    )

package_maplibre_stored_map(cog_path, viewer_dir, *, plan, map_type, name=None, profile=None, geometry=None, layer_id=None, source_cog=None, units='ft', visible=False, domain_policy='fixed', max_zoom=None, scratch_dir=None, overwrite=False)

Publish one RASMapper Stored Map as display PMTiles plus numeric COG.

The numeric COG remains authoritative for Identify and future view-local styling. The PMTiles derivative is a fast default visualization and is explicitly recorded as RASMapper/RasProcess interpolation, not raw HDF computation-element data.

Source code in ras2cng/maplibre.py

def package_maplibre_stored_map(
    cog_path: Path,
    viewer_dir: Path,
    *,
    plan: str,
    map_type: str,
    name: str | None = None,
    profile: str | None = None,
    geometry: str | None = None,
    layer_id: str | None = None,
    source_cog: str | None = None,
    units: str = "ft",
    visible: bool = False,
    domain_policy: str = "fixed",
    max_zoom: int | None = None,
    scratch_dir: Path | None = None,
    overwrite: bool = False,
) -> RasterPackageSummary:
    """Publish one RASMapper Stored Map as display PMTiles plus numeric COG.

    The numeric COG remains authoritative for Identify and future view-local
    styling. The PMTiles derivative is a fast default visualization and is
    explicitly recorded as RASMapper/RasProcess interpolation, not raw HDF
    computation-element data.
    """

    provenance: dict[str, Any] = {
        "mapType": map_type,
        "source": "RASMapper/RasProcess Stored Map",
        "interpolationAuthority": "RASMapper/RasProcess",
    }
    return _package_maplibre_numeric_raster(
        cog_path,
        viewer_dir,
        plan=plan,
        map_type=map_type,
        name=name,
        profile=profile,
        geometry=geometry,
        layer_id=layer_id,
        source_cog=source_cog,
        units=units,
        visible=visible,
        domain_policy=domain_policy,
        max_zoom=max_zoom,
        scratch_dir=scratch_dir,
        overwrite=overwrite,
        source_kind="stored-map",
        provenance=provenance,
        result_kind="rasmapper_stored_map",
        legend_type="continuous",
        legend_mode="stretched",
    )

package_maplibre_stored_vector(vector_path, viewer_dir, *, plan, map_type, name=None, profile=None, geometry=None, layer_id=None, crs=None, visible=False, min_zoom=0, max_zoom=17, scratch_dir=None, overwrite=False)

Publish a RASMapper vector Stored Map as queryable PMTiles.

Source code in ras2cng/maplibre.py

def package_maplibre_stored_vector(
    vector_path: Path,
    viewer_dir: Path,
    *,
    plan: str,
    map_type: str,
    name: str | None = None,
    profile: str | None = None,
    geometry: str | None = None,
    layer_id: str | None = None,
    crs: str | None = None,
    visible: bool = False,
    min_zoom: int = 0,
    max_zoom: int = 17,
    scratch_dir: Path | None = None,
    overwrite: bool = False,
) -> VectorResultPackageSummary:
    """Publish a RASMapper vector Stored Map as queryable PMTiles."""

    vector_path = Path(vector_path)
    viewer_dir = Path(viewer_dir)
    manifest_path = viewer_dir / "manifest.json"
    if not vector_path.is_file():
        raise FileNotFoundError(f"Stored Map vector does not exist: {vector_path}")
    if not manifest_path.is_file():
        raise FileNotFoundError(f"MapLibre viewer manifest does not exist: {manifest_path}")

    plan_id = _slug(plan)
    if plan_id.isdigit():
        plan_id = f"p{plan_id.zfill(2)}"
    elif plan_id.startswith("p") and plan_id[1:].isdigit():
        plan_id = f"p{plan_id[1:].zfill(2)}"
    if not plan_id:
        raise ValueError("A plan identifier is required for every Stored Map vector")
    map_slug = _slug(map_type)
    profile_slug = _slug(profile or "")
    layer_id = layer_id or "-".join(
        value for value in (plan_id, map_slug, profile_slug) if value
    )
    if not layer_id:
        raise ValueError("Could not derive a Stored Map vector layer identifier")

    manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
    tilesets = manifest.setdefault("tilesets", [])
    existing = next((item for item in tilesets if item.get("id") == layer_id), None)
    if existing and not overwrite:
        raise FileExistsError(
            f"Viewer already has vector layer {layer_id}; pass overwrite=True to replace it."
        )

    _require_cli(_tippecanoe_command())
    _require_cli(_pmtiles_command())
    if vector_path.suffix.lower() in {".parquet", ".geoparquet"}:
        frame = gpd.read_parquet(vector_path)
    else:
        frame = gpd.read_file(vector_path)
    frame = _to_wgs84(frame, vector_path, crs)
    if frame.empty:
        raise ValueError(f"Stored Map vector has no publishable features: {vector_path}")

    scratch_parent = Path(scratch_dir).resolve() if scratch_dir else None
    if scratch_parent:
        scratch_parent.mkdir(parents=True, exist_ok=True)
    output = viewer_dir / "tiles" / f"{layer_id}.pmtiles"
    if output.exists() and not overwrite:
        raise FileExistsError(f"Vector PMTiles already exists: {output}")
    with tempfile.TemporaryDirectory(
        prefix=f"ras2cng-{layer_id}-",
        dir=str(scratch_parent) if scratch_parent else None,
    ) as temporary:
        work_dir = Path(temporary)
        ndgeojson = work_dir / f"{layer_id}.ndgeojson"
        feature_count, geometry_types, bounds = _write_ndgeojson(frame, ndgeojson)
        _run_tippecanoe(
            output,
            [(layer_id, ndgeojson)],
            min_zoom,
            max_zoom,
            work_dir / "tippecanoe",
        )

    provenance: dict[str, Any] = {
        "source": "RASMapper/RasProcess Stored Map",
        "interpolationAuthority": "RASMapper/RasProcess",
        "mapType": map_type,
        "plan": plan_id,
        "sourceVector": vector_path.name,
    }
    if profile:
        provenance["profile"] = profile
    if geometry:
        provenance["geometry"] = geometry
    layer = {
        "id": layer_id,
        "name": name or " ".join(value for value in (map_type, profile) if value),
        "sourceLayer": layer_id,
        "groupId": f"ras-results-{plan_id}",
        "geometryId": geometry,
        "visible": visible,
        "kind": map_slug.replace("-", "_"),
        "sourceKind": "stored-map",
        "style": {
            "fill": "#2563eb",
            "fillOpacity": 0.12,
            "line": "#1d4ed8",
            "lineWidth": 1.6,
        },
        "featureCount": feature_count,
        "geometryTypes": geometry_types,
        "bounds": list(bounds),
        "sort": 90,
        "queryable": True,
        "provenance": provenance,
    }
    tileset = {
        "id": layer_id,
        "type": "vector",
        "href": f"tiles/{layer_id}.pmtiles",
        "bytes": output.stat().st_size,
        "layers": [layer],
        "groupId": f"ras-results-{plan_id}",
        "resultKind": "stored_map",
    }
    if existing:
        tilesets[tilesets.index(existing)] = tileset
    else:
        tilesets.append(tileset)
    groups = manifest.setdefault("groups", [])
    group_id = f"ras-results-{plan_id}"
    if not any(group.get("id") == group_id for group in groups):
        groups.append(
            {
                "id": group_id,
                "name": f"Plan {plan_id}",
                "visible": False,
                "resultKind": "stored_map",
            }
        )
    apply_manifest_v2(manifest, archive=_viewer_archive_metadata(viewer_dir))
    manifest_path.write_text(json.dumps(manifest, indent=2) + "\n", encoding="utf-8")
    return VectorResultPackageSummary(
        manifest_path=manifest_path,
        pmtiles_path=output,
        source_vector=vector_path,
        feature_count=feature_count,
        layer_id=layer_id,
    )

package_maplibre_calculated_map(cog_path, viewer_dir, *, plan, recipe_id, name=None, profile=None, geometry=None, layer_id=None, source_cog=None, units=None, provenance_path=None, visible=False, domain_policy='fixed', max_zoom=None, scratch_dir=None, overwrite=False)

Publish a controlled recipe output under a plan's Calculated Layers.

Source code in ras2cng/maplibre.py

def package_maplibre_calculated_map(
    cog_path: Path,
    viewer_dir: Path,
    *,
    plan: str,
    recipe_id: str,
    name: str | None = None,
    profile: str | None = None,
    geometry: str | None = None,
    layer_id: str | None = None,
    source_cog: str | None = None,
    units: str | None = None,
    provenance_path: Path | None = None,
    visible: bool = False,
    domain_policy: str = "fixed",
    max_zoom: int | None = None,
    scratch_dir: Path | None = None,
    overwrite: bool = False,
) -> RasterPackageSummary:
    """Publish a controlled recipe output under a plan's Calculated Layers."""

    from ras2cng.raster_recipes import get_raster_recipe

    cog_path = Path(cog_path)
    recipe = get_raster_recipe(recipe_id)
    provenance_path = Path(provenance_path) if provenance_path else cog_path.with_suffix(".provenance.json")
    recipe_run: dict[str, Any] = {}
    if provenance_path.is_file():
        recipe_run = json.loads(provenance_path.read_text(encoding="utf-8"))
        recorded_id = (recipe_run.get("recipe") or {}).get("recipe_id")
        if recorded_id and recorded_id != recipe_id:
            raise ValueError(
                f"Calculated raster provenance records recipe {recorded_id!r}, not {recipe_id!r}"
            )
    profile = profile or recipe_run.get("profile")
    units = units or (recipe_run.get("output") or {}).get("units")
    if not units:
        raise ValueError("Calculated raster units are required or must exist in its provenance sidecar")
    if recipe.categorical and domain_policy != "fixed":
        raise ValueError("Categorical calculated rasters require domain_policy='fixed'")
    provenance: dict[str, Any] = {
        "mapType": recipe.recipe_id,
        "recipeId": recipe.recipe_id,
        "recipeVersion": recipe.version,
        "source": "ras2cng controlled raster recipe",
        "arithmeticAuthority": "ras2cng",
        "interpolationAuthority": "RASMapper/RasProcess source rasters",
        "parameters": recipe_run.get("parameters", dict(recipe.parameter_defaults)),
        "inputs": recipe_run.get("inputs", {}),
    }
    return _package_maplibre_numeric_raster(
        cog_path,
        viewer_dir,
        plan=plan,
        map_type=recipe.recipe_id,
        name=name or recipe.name,
        profile=profile,
        geometry=geometry,
        layer_id=layer_id,
        source_cog=source_cog,
        units=units,
        visible=visible,
        domain_policy=domain_policy,
        max_zoom=max_zoom,
        scratch_dir=scratch_dir,
        overwrite=overwrite,
        source_kind="calculated",
        provenance=provenance,
        result_kind="calculated_raster",
        legend_type="categorical" if recipe.categorical else "continuous",
        legend_mode="discrete" if recipe.categorical else "stretched",
    )

package_maplibre_viewer(archive_dir, output_dir, *, geometry_hdfs, title=None, source_project=None, crs=None, include_vector_results=False, primary_geometry=None, show_all_primary_geometry=False, min_zoom=0, max_zoom=17, scratch_dir=None)

Create a MapLibre viewer bundle from a completed ras2cng archive.

geometry_hdfs maps archive IDs such as g01 to their original HDF geometry files. Requiring that mapping ensures every model footprint in the browser bundle is produced by HdfProject.get_project_extent rather than approximated from delivery tiles.

Source code in ras2cng/maplibre.py

def package_maplibre_viewer(
    archive_dir: Path,
    output_dir: Path,
    *,
    geometry_hdfs: Mapping[str, Path],
    title: str | None = None,
    source_project: str | None = None,
    crs: str | None = None,
    include_vector_results: bool = False,
    primary_geometry: str | None = None,
    show_all_primary_geometry: bool = False,
    min_zoom: int = 0,
    max_zoom: int = 17,
    scratch_dir: Path | None = None,
) -> PackageSummary:
    """Create a MapLibre viewer bundle from a completed ras2cng archive.

    ``geometry_hdfs`` maps archive IDs such as ``g01`` to their original HDF
    geometry files. Requiring that mapping ensures every model footprint in the
    browser bundle is produced by ``HdfProject.get_project_extent`` rather than
    approximated from delivery tiles.
    """

    archive_dir = Path(archive_dir)
    output_dir = Path(output_dir)
    archive_manifest_path = archive_dir / "manifest.json"
    if not archive_manifest_path.is_file():
        raise FileNotFoundError(f"ras2cng archive manifest not found: {archive_manifest_path}")
    if output_dir.exists() and any(output_dir.iterdir()):
        raise FileExistsError(f"MapLibre output directory is not empty: {output_dir}")

    archive = json.loads(archive_manifest_path.read_text(encoding="utf-8"))
    geometry_entries = archive.get("geometry", [])
    if not geometry_entries:
        raise ValueError(f"Archive has no geometry entries: {archive_manifest_path}")

    missing_hdfs = [entry["geom_id"] for entry in geometry_entries if entry["geom_id"] not in geometry_hdfs]
    if missing_hdfs:
        raise ValueError(
            "Missing geometry HDF mapping(s) required for API-derived footprints: "
            + ", ".join(missing_hdfs)
        )
    for geom_id, hdf_path in geometry_hdfs.items():
        if not Path(hdf_path).is_file():
            raise FileNotFoundError(f"Geometry HDF for {geom_id} does not exist: {hdf_path}")

    metadata = _project_metadata(archive_dir)
    project_crs = crs or metadata.get("crs") or archive.get("project", {}).get("crs")
    _require_cli(_tippecanoe_command())
    _require_cli(_pmtiles_command())
    output_dir.mkdir(parents=True, exist_ok=True)
    tiles_dir = output_dir / "tiles"
    viewer_title = title or metadata.get("title") or archive.get("project", {}).get("name") or archive_dir.name
    source_project = source_project or metadata.get("href") or "../project.json"
    if scratch_dir is not None:
        scratch_dir = Path(scratch_dir).resolve()
        scratch_dir.mkdir(parents=True, exist_ok=True)
        if not scratch_dir.is_dir():
            raise ValueError(f"MapLibre scratch directory is not a directory: {scratch_dir}")

    geometry_cache: dict[tuple[str, str], gpd.GeoDataFrame] = {}
    result_geometry_keys = {
        (str(plan.get("geom_id", "")).lower(), str(variable["geometry_filter"]))
        for plan in archive.get("results", [])
        for variable in plan.get("variables", [])
        if include_vector_results and variable.get("geometry_filter")
    }
    geometry_overview_sources: list[tuple[str, Path]] = []
    geometry_detail_sources: list[tuple[str, Path]] = []
    geometry_overview_layers: list[dict[str, Any]] = []
    geometry_detail_layers: list[dict[str, Any]] = []
    geometry_layers: list[dict[str, Any]] = []
    groups: list[dict[str, Any]] = []
    extent_features: list[dict[str, Any]] = []
    all_bounds: list[Sequence[float]] = []

    with tempfile.TemporaryDirectory(
        prefix="ras2cng-maplibre-",
        dir=str(scratch_dir) if scratch_dir is not None else None,
    ) as temporary:
        work_dir = Path(temporary)
        for geometry_index, entry in enumerate(geometry_entries):
            geom_id = entry["geom_id"].lower()
            group_id = f"ras-geometry-{geom_id}"
            group_layers: list[dict[str, Any]] = []
            groups.append(
                {
                    "id": group_id,
                    "name": _geometry_display_label(
                        geom_id,
                        entry.get("geom_title"),
                    ),
                    "visible": geometry_index == 0,
                }
            )

            extent = _extent_from_hdf(Path(geometry_hdfs[entry["geom_id"]]), project_crs)
            extent["geometry_id"] = geom_id
            extent_source = f"{group_id}-model-extents"
            extent_path = work_dir / "geometry" / f"{extent_source}.ndgeojson"
            count, geometry_types, bounds = _write_ndgeojson(extent, extent_path)
            geometry_overview_sources.append((extent_source, extent_path))
            all_bounds.append(bounds)
            group_layers.append(
                {
                    "id": extent_source,
                    "name": _display_name("model_extents"),
                    "sourceLayer": extent_source,
                    "groupId": group_id,
                    "visible": False,
                    "kind": "model_extents",
                    "style": _GEOMETRY_STYLES["model_extents"].copy(),
                    "featureCount": count,
                    "geometryTypes": geometry_types,
                    "bounds": bounds,
                    "sort": 0,
                    "queryable": True,
                    "extentSource": "HdfProject.get_project_extent(geometry_type='footprint')",
                }
            )
            extent_features.extend(extent.iterfeatures(drop_id=True, na="null", show_bbox=False))

            archive_geometry_path = archive_dir / entry["parquet"]
            for layer in entry.get("layers", []):
                kind = layer.get("layer") or layer.get("filter_value")
                filter_value = layer.get("filter_value") or kind
                if not kind or not filter_value:
                    continue
                source_layer = f"{group_id}-{_slug(kind)}"
                source_path = work_dir / "geometry" / f"{source_layer}.ndgeojson"
                cache_geometry = (geom_id, kind) in result_geometry_keys
                gdf: gpd.GeoDataFrame | None = None
                if _is_detail_geometry(kind) and not cache_geometry:
                    count, geometry_types, bounds = _stream_dense_layer_ndgeojson(
                        archive_geometry_path,
                        filter_value,
                        kind,
                        source_path,
                        project_crs,
                    )
                else:
                    gdf = _to_wgs84(
                        _read_layer(archive_geometry_path, filter_value),
                        archive_geometry_path,
                        project_crs,
                    )
                    if gdf.empty:
                        continue
                    if cache_geometry:
                        geometry_cache[(geom_id, kind)] = gdf
                    count, geometry_types, bounds = _write_ndgeojson(gdf, source_path)
                source_layers = (
                    geometry_detail_sources if _is_detail_geometry(kind) else geometry_overview_sources
                )
                source_layers.append((source_layer, source_path))
                all_bounds.append(bounds)
                group_layers.append(
                    {
                        "id": source_layer,
                        "name": _display_name(kind),
                        "sourceLayer": source_layer,
                        "groupId": group_id,
                        "visible": False,
                        "kind": kind,
                        "style": _GEOMETRY_STYLES.get(
                            kind,
                            {"fill": "#94a3b8", "fillOpacity": 0.12, "line": "#475569", "lineWidth": 1.0},
                        ).copy(),
                        "featureCount": count,
                        "geometryTypes": geometry_types,
                        "bounds": bounds,
                        "sort": _geometry_sort(kind),
                        "queryable": True,
                    }
                )
                if gdf is not None and not cache_geometry:
                    del gdf

            geometry_layers.extend(group_layers)
            for layer in group_layers:
                target_layers = (
                    geometry_detail_layers
                    if _is_detail_geometry(layer["kind"])
                    else geometry_overview_layers
                )
                target_layers.append(layer)

        for terrain_entry in archive.get("terrain_sources", []):
            terrain_name = str(terrain_entry.get("terrain_name") or "Terrain")
            terrain_slug = _slug(terrain_name) or "terrain"
            group_id = f"ras-terrain-sources-{terrain_slug}"
            groups.append(
                {
                    "id": group_id,
                    "name": f"{terrain_name} Sources",
                    "visible": False,
                }
            )
            for layer in terrain_entry.get("layers", []):
                kind = str(layer.get("layer") or "")
                parquet_href = layer.get("parquet")
                if not kind or not parquet_href:
                    continue
                source_path_archive = archive_dir / parquet_href
                gdf = _to_wgs84(
                    gpd.read_parquet(source_path_archive),
                    source_path_archive,
                    project_crs,
                )
                if gdf.empty:
                    continue
                source_layer = f"{group_id}-{_slug(kind)}"
                source_path = work_dir / "terrain-sources" / f"{source_layer}.ndgeojson"
                count, geometry_types, bounds = _write_ndgeojson(gdf, source_path)
                geometry_overview_sources.append((source_layer, source_path))
                all_bounds.append(bounds)
                viewer_layer = {
                    "id": source_layer,
                    "name": _display_name(kind),
                    "sourceLayer": source_layer,
                    "groupId": group_id,
                    "visible": False,
                    "kind": kind,
                    "sourceKind": "terrain-source",
                    "style": _GEOMETRY_STYLES[kind].copy(),
                    "featureCount": count,
                    "geometryTypes": geometry_types,
                    "bounds": bounds,
                    "sort": _geometry_sort(kind),
                    "queryable": True,
                    "provenance": {
                        "source": "Native HEC-RAS terrain TIFF members",
                        "terrain": terrain_name,
                    },
                }
                geometry_overview_layers.append(viewer_layer)
                geometry_layers.append(viewer_layer)

        for terrain_entry in archive.get("terrain_modifications", []):
            terrain_name = str(terrain_entry.get("terrain_name") or "Terrain")
            terrain_slug = _slug(terrain_name) or "terrain"
            group_id = f"ras-terrain-modifications-{terrain_slug}"
            groups.append(
                {
                    "id": group_id,
                    "name": f"{terrain_name} Modifications",
                    "visible": False,
                }
            )
            for layer in terrain_entry.get("layers", []):
                kind = str(layer.get("layer") or "")
                parquet_href = layer.get("parquet")
                if not kind or not parquet_href:
                    continue
                source_path_archive = archive_dir / parquet_href
                gdf = _to_wgs84(
                    gpd.read_parquet(source_path_archive),
                    source_path_archive,
                    project_crs,
                )
                if gdf.empty:
                    continue
                source_layer = f"{group_id}-{_slug(kind)}"
                source_path = work_dir / "terrain-modifications" / f"{source_layer}.ndgeojson"
                count, geometry_types, bounds = _write_ndgeojson(gdf, source_path)
                geometry_overview_sources.append((source_layer, source_path))
                all_bounds.append(bounds)
                viewer_layer = {
                    "id": source_layer,
                    "name": _display_name(kind),
                    "sourceLayer": source_layer,
                    "groupId": group_id,
                    "visible": False,
                    "kind": kind,
                    "sourceKind": "terrain-modification",
                    "style": _GEOMETRY_STYLES[kind].copy(),
                    "featureCount": count,
                    "geometryTypes": geometry_types,
                    "bounds": bounds,
                    "sort": _geometry_sort(kind),
                    "queryable": True,
                    "provenance": {
                        "source": "HEC-RAS terrain modification HDF",
                        "terrain": terrain_name,
                        "sourceHdf": Path(str(terrain_entry.get("source_hdf") or "")).name,
                    },
                }
                geometry_overview_layers.append(viewer_layer)
                geometry_layers.append(viewer_layer)

        geometry_pmtiles = tiles_dir / "geometry.pmtiles"
        _run_tippecanoe(
            geometry_pmtiles,
            geometry_overview_sources,
            min_zoom,
            max_zoom,
            work_dir / "tippecanoe-overview",
        )
        geometry_detail_pmtiles: Path | None = None
        if geometry_detail_sources:
            geometry_detail_pmtiles = tiles_dir / "geometry-detail.pmtiles"
            _run_tippecanoe(
                geometry_detail_pmtiles,
                geometry_detail_sources,
                max(min_zoom, 13),
                max_zoom,
                work_dir / "tippecanoe-detail",
            )

        result_sources: list[tuple[str, Path]] = []
        result_layers: list[dict[str, Any]] = []
        if include_vector_results:
            for plan in archive.get("results", []):
                plan_id = str(plan.get("plan_id", "plan")).lower()
                result_group_id = f"ras-results-{plan_id}"
                plan_layers: list[dict[str, Any]] = []
                for variable in plan.get("variables", []):
                    variable_path = variable.get("parquet") or plan.get("parquet")
                    variable_filter = (
                        variable.get("filter_value")
                        if not variable.get("parquet") and plan.get("parquet")
                        else None
                    )
                    geometry_kind = variable.get("geometry_filter")
                    index_column = str(variable.get("index_column") or "")
                    join_columns = variable.get("join_columns") or {}
                    geom_id = str(plan.get("geom_id", "")).lower()
                    if not variable_path or not geometry_kind or not (index_column or join_columns):
                        continue
                    geometry = geometry_cache.get((geom_id, geometry_kind))
                    if geometry is None:
                        continue
                    raw_path = archive_dir / variable_path
                    variable_name = variable.get("variable") or variable.get("filter_value") or raw_path.stem
                    profile_column = str(variable.get("profile_column") or "")
                    profiles: list[Any] = [None]
                    if profile_column:
                        profile_columns = [profile_column]
                        if variable_filter:
                            profile_columns.append("layer")
                        profile_frame = pd.read_parquet(raw_path, columns=profile_columns)
                        if variable_filter:
                            profile_frame = profile_frame.loc[
                                profile_frame["layer"] == variable_filter
                            ]
                        profile_values = profile_frame[profile_column]
                        profiles = list(pd.unique(profile_values.dropna()))
                    for profile_index, profile in enumerate(profiles):
                        filters: dict[str, Any] = {}
                        if variable_filter:
                            filters["layer"] = variable_filter
                        if profile_column:
                            filters[profile_column] = profile
                        joined = _join_raw_result(
                            raw_path,
                            geometry,
                            index_column,
                            join_columns=join_columns,
                            filters=filters or None,
                        )
                        if joined.empty:
                            continue
                        profile_suffix = f"-{_slug(str(profile))}" if profile is not None else ""
                        source_layer = f"{result_group_id}-{_slug(variable_name)}{profile_suffix}"
                        source_path = work_dir / "results" / f"{source_layer}.ndgeojson"
                        count, geometry_types, bounds = _write_ndgeojson(joined, source_path)
                        result_sources.append((source_layer, source_path))
                        all_bounds.append(bounds)
                        layer_name = _display_name(variable_name)
                        if profile is not None:
                            layer_name = f"{layer_name} - {profile}"
                        raw_result = {
                            "source": variable.get("source") or "Raw HEC-RAS HDF summary result values",
                            "plan": plan_id,
                            "variable": variable_name,
                            "geometryJoin": geometry_kind,
                            "archiveParquet": variable_path,
                        }
                        if index_column:
                            raw_result["indexColumn"] = index_column
                        if join_columns:
                            raw_result["joinColumns"] = join_columns
                        if variable_filter:
                            raw_result["archiveFilter"] = {
                                "column": "layer",
                                "value": variable_filter,
                            }
                        if profile is not None:
                            raw_result["profile"] = profile
                        plan_layers.append(
                            {
                                "id": source_layer,
                                "name": layer_name,
                                "sourceLayer": source_layer,
                                "groupId": result_group_id,
                                "visible": False,
                                "kind": f"{plan_id}_{variable_name}{profile_suffix}",
                                "style": _result_style(variable_name),
                                "featureCount": count,
                                "geometryTypes": geometry_types,
                                "bounds": bounds,
                                "sort": 100 + profile_index,
                                "queryable": True,
                                "rawResult": raw_result,
                            }
                        )
                if plan_layers:
                    groups.append(
                        {
                            "id": result_group_id,
                            "name": f"Vector Results {plan_id}",
                            "visible": False,
                            "resultKind": "raw_hdf",
                        }
                    )
                    result_layers.extend(plan_layers)

        result_pmtiles: Path | None = None
        if result_sources:
            result_pmtiles = tiles_dir / "results.pmtiles"
            _run_tippecanoe(
                result_pmtiles,
                result_sources,
                min_zoom,
                max_zoom,
                work_dir / "tippecanoe-results",
            )

    final_bounds = _merge_bounds(all_bounds)
    center = [
        (final_bounds[0] + final_bounds[2]) / 2.0,
        (final_bounds[1] + final_bounds[3]) / 2.0,
    ]
    geometry_tilesets: list[dict[str, Any]] = [
        {
            "id": "geometry",
            "type": "vector",
            "href": "tiles/geometry.pmtiles",
            "bytes": geometry_pmtiles.stat().st_size,
            "layers": geometry_overview_layers,
        }
    ]
    if geometry_detail_pmtiles:
        geometry_tilesets.append(
            {
                "id": "geometry-detail",
                "type": "vector",
                "href": "tiles/geometry-detail.pmtiles",
                "bytes": geometry_detail_pmtiles.stat().st_size,
                "layers": geometry_detail_layers,
                "minzoom": max(min_zoom, 13),
            }
        )

    manifest: dict[str, Any] = {
        "schema": "rascommander.maplibre.project/1",
        "generatedBy": "ras2cng maplibre",
        "sourceProject": source_project,
        "title": viewer_title,
        "bounds": list(final_bounds),
        "center": center,
        "zoom": _default_zoom(final_bounds),
        "sourceCrs": project_crs,
        "tilesets": geometry_tilesets,
        "groups": groups,
        "notes": (
            "Geometry is delivered as PMTiles. Vector Results are raw HEC-RAS "
            "HDF summary values joined to their source geometry for display; "
            "they are not RASMapper-interpolated raster results."
        ),
    }
    primary_geometry_group_id = None
    if primary_geometry:
        normalized_primary = _slug(primary_geometry)
        primary_geometry_group_id = (
            normalized_primary
            if normalized_primary.startswith("ras-geometry-")
            else f"ras-geometry-{normalized_primary}"
        )
    apply_maplibre_default_visibility(
        manifest,
        primary_geometry_group_id=(
            primary_geometry_group_id or _preferred_result_geometry_group_id(archive)
        ),
        show_all_primary_geometry=show_all_primary_geometry,
    )
    if result_pmtiles:
        manifest["tilesets"].append(
            {
                "id": "results",
                "type": "vector",
                "href": "tiles/results.pmtiles",
                "bytes": result_pmtiles.stat().st_size,
                "layers": result_layers,
                "resultKind": "raw_hdf",
            }
        )
    apply_manifest_v2(manifest, archive=archive)

    (output_dir / "model_extent.geojson").write_text(
        json.dumps(
            {
                "type": "FeatureCollection",
                "name": f"{_slug(viewer_title)}-model-extents",
                "features": extent_features,
            },
            indent=2,
            default=str,
        )
        + "\n",
        encoding="utf-8",
    )
    manifest_path = output_dir / "manifest.json"
    manifest_path.write_text(json.dumps(manifest, indent=2) + "\n", encoding="utf-8")
    return PackageSummary(
        manifest_path=manifest_path,
        geometry_pmtiles=geometry_pmtiles,
        result_pmtiles=result_pmtiles,
        geometry_layer_count=len(geometry_layers),
        result_layer_count=len(result_layers),
        bounds=final_bounds,
    )

apply_maplibre_default_visibility(manifest, *, primary_geometry_group_id=None, show_all_primary_geometry=False)

Apply the standard initial geometry view to a MapLibre manifest.

The viewer should begin with one geometry configuration, its authoritative model-limit footprint, and enough model context to orient a reviewer. A 1D geometry uses centerlines; a 2D geometry uses mesh context plus mesh refinement controls when they are present. Dense faces, cross sections, boundary conditions, and structures remain opt-in.

Source code in ras2cng/maplibre.py

def apply_maplibre_default_visibility(
    manifest: dict[str, Any],
    *,
    primary_geometry_group_id: str | None = None,
    show_all_primary_geometry: bool = False,
) -> None:
    """Apply the standard initial geometry view to a MapLibre manifest.

    The viewer should begin with one geometry configuration, its authoritative
    model-limit footprint, and enough model context to orient a reviewer. A
    1D geometry uses centerlines; a 2D geometry uses mesh context plus mesh
    refinement controls when they are present. Dense faces, cross sections,
    boundary conditions, and structures remain opt-in.
    """

    geometry_groups: dict[str, list[dict[str, Any]]] = {}
    for tileset in manifest.get("tilesets", []):
        if tileset.get("type") != "vector":
            continue
        for layer in tileset.get("layers", []):
            group_id = str(layer.get("groupId") or "")
            if group_id.startswith("ras-geometry-"):
                geometry_groups.setdefault(group_id, []).append(layer)

    if not geometry_groups:
        return

    if primary_geometry_group_id not in geometry_groups:
        configured_groups = [
            str(group.get("id"))
            for group in manifest.get("groups", [])
            if group.get("visible") and str(group.get("id")) in geometry_groups
        ]
        primary_geometry_group_id = configured_groups[0] if configured_groups else next(iter(geometry_groups))

    for group_id, layers in geometry_groups.items():
        for layer in layers:
            layer["visible"] = False
        for group in manifest.get("groups", []):
            if group.get("id") == group_id:
                group["visible"] = group_id == primary_geometry_group_id

    primary_layers = geometry_groups[primary_geometry_group_id]
    if show_all_primary_geometry:
        for layer in primary_layers:
            layer["visible"] = True
        return

    kinds = {str(layer.get("kind") or "") for layer in primary_layers}
    is_2d = bool({"mesh_areas", "mesh_cells", "mesh_faces", "breaklines", "refinement_regions"} & kinds)
    default_kinds = {"model_extents", "pipe_conduits", "pipe_nodes"}
    if is_2d:
        default_kinds.update({"mesh_areas", "mesh_cells", "breaklines", "refinement_regions"})
    else:
        default_kinds.update({"centerlines", "river_centerlines", "river_reaches"})

    for layer in primary_layers:
        if layer.get("kind") in default_kinds:
            layer["visible"] = True

Viewer Manifest Contract

ras2cng.viewer_manifest

Build and validate the RAS Commander MapLibre viewer manifest contract.

apply_manifest_v2(manifest, *, archive=None)

Add the v2 semantic contract while retaining v1 compatibility fields.

Source code in ras2cng/viewer_manifest.py

def apply_manifest_v2(
    manifest: dict[str, Any],
    *,
    archive: Mapping[str, Any] | None = None,
) -> dict[str, Any]:
    """Add the v2 semantic contract while retaining v1 compatibility fields."""

    compatibility = manifest.get("compatibility") or {}
    legacy_schema = str(
        compatibility.get("legacySchema")
        or manifest.get("schema")
        or LEGACY_MAPLIBRE_SCHEMA
    )
    existing_layers = manifest.get("layers") if isinstance(manifest.get("layers"), dict) else {}
    resources: dict[str, dict[str, Any]] = {}
    layers: dict[str, dict[str, Any]] = {}
    legends: dict[str, dict[str, Any]] = deepcopy(manifest.get("legends") or {})

    for tileset in manifest.get("tilesets", []):
        display_resource, numeric_resource = _add_tileset_resources(resources, tileset)
        if tileset.get("type") == "vector":
            for legacy_layer in tileset.get("layers", []):
                layer_id = str(legacy_layer.get("id") or "")
                if not layer_id:
                    raise ValueError("Every vector viewer layer requires an id.")
                if layer_id in layers:
                    raise ValueError(f"Duplicate viewer layer id: {layer_id}")
                layers[layer_id] = _vector_layer_record(
                    legacy_layer,
                    tileset,
                    display_resource,
                )
        elif tileset.get("type") == "raster":
            layer_id = str(tileset.get("id") or "")
            if not layer_id:
                raise ValueError("Every raster viewer layer requires an id.")
            if layer_id in layers:
                raise ValueError(f"Duplicate viewer layer id: {layer_id}")
            legend_id, legend = _raster_legend(tileset)
            legends[legend_id] = legend
            layers[layer_id] = _raster_layer_record(
                tileset,
                display_resource,
                numeric_resource,
                legend_id,
            )

    raster_query = manifest.get("rasterQuery") or {}
    for resource in resources.values():
        if resource.get("type") != "cog":
            continue
        if raster_query.get("sourceCrs"):
            resource.setdefault("crs", raster_query["sourceCrs"])
        if raster_query.get("sourceProj4"):
            resource.setdefault("proj4", raster_query["sourceProj4"])

    # Incremental terrain/result packaging rebuilds the semantic contract from
    # compatibility tilesets. Retain context already enriched by the archive
    # so a later layer import cannot erase HEC-RAS plan and geometry titles.
    for layer_id, layer in layers.items():
        previous = existing_layers.get(layer_id)
        if not isinstance(previous, dict):
            continue
        _copy_present(
            previous,
            layer,
            "plan",
            "planTitle",
            "geometry",
            "geometryTitle",
            "terrain",
        )

    _add_hybrid_basemap(resources, layers)
    _enrich_layer_context(layers, archive)
    associations = _build_associations(manifest, layers, archive)
    tree = _build_tree(manifest, layers, archive, associations)
    interaction = _build_interaction(manifest, layers)
    capabilities = _build_capabilities(layers, archive)

    provenance = deepcopy(manifest.get("provenance") or {})
    provenance.setdefault("generatedBy", manifest.get("generatedBy", "ras2cng maplibre"))
    provenance.setdefault("sourceProject", manifest.get("sourceProject"))
    provenance.setdefault("sourceCrs", manifest.get("sourceCrs"))
    if archive:
        provenance.setdefault("archiveSchemaVersion", archive.get("schema_version"))
    provenance.setdefault(
        "resultSemantics",
        {
            "rawHdf": "Values at HEC-RAS computation elements; no surface interpolation.",
            "storedMap": "Raster surface generated by RASMapper/RasProcess.",
        },
    )

    manifest["schema"] = MAPLIBRE_SCHEMA
    manifest["resources"] = resources
    manifest["layers"] = layers
    manifest["tree"] = tree
    manifest["associations"] = associations
    manifest["legends"] = legends
    manifest["interaction"] = interaction
    manifest["capabilities"] = capabilities
    manifest.setdefault("timeAxes", {})
    manifest["provenance"] = provenance
    manifest["compatibility"] = {
        "legacySchema": legacy_schema,
        "legacyFields": ["tilesets", "groups"],
        "legacyViewerSupported": True,
    }
    validate_manifest_v2(manifest)
    return manifest

validate_manifest_v2(manifest)

Raise ValueError when required v2 references are inconsistent.

Source code in ras2cng/viewer_manifest.py

def validate_manifest_v2(manifest: Mapping[str, Any]) -> None:
    """Raise ``ValueError`` when required v2 references are inconsistent."""

    if manifest.get("schema") != MAPLIBRE_SCHEMA:
        raise ValueError(f"Viewer manifest schema must be {MAPLIBRE_SCHEMA!r}.")
    resources = manifest.get("resources")
    layers = manifest.get("layers")
    tree = manifest.get("tree")
    if not isinstance(resources, Mapping):
        raise ValueError("Viewer manifest resources must be an object.")
    if not isinstance(layers, Mapping):
        raise ValueError("Viewer manifest layers must be an object.")
    if not isinstance(tree, list):
        raise ValueError("Viewer manifest tree must be an array.")

    expected_roots = [root_id for root_id, _ in ROOT_DEFINITIONS]
    observed_roots = [node.get("id") for node in tree]
    if observed_roots != expected_roots:
        raise ValueError(
            "Viewer manifest roots must be ordered as " + ", ".join(expected_roots)
        )

    for layer_id, layer in layers.items():
        resource_id = layer.get("resource")
        if resource_id not in resources:
            raise ValueError(f"Layer {layer_id!r} references missing resource {resource_id!r}.")
        query = layer.get("query") or {}
        numeric_resource = query.get("numericResource")
        if numeric_resource and numeric_resource not in resources:
            raise ValueError(
                f"Layer {layer_id!r} references missing numeric resource {numeric_resource!r}."
            )
        domain_policy = (layer.get("style") or {}).get("domainPolicy", "fixed")
        if domain_policy == "current-view":
            legend_id = (layer.get("style") or {}).get("legendRef")
            legend = (manifest.get("legends") or {}).get(legend_id, {})
            if legend.get("type") == "categorical":
                raise ValueError(
                    f"Categorical layer {layer_id!r} cannot use current-view styling."
                )
            numeric = resources.get(numeric_resource, {})
            if numeric.get("type") != "cog":
                raise ValueError(
                    f"Current-view layer {layer_id!r} requires an authoritative numeric COG."
                )
            if not numeric.get("serviceAsset") or not numeric.get("serviceRevision"):
                raise ValueError(
                    f"Current-view layer {layer_id!r} requires serviceAsset and serviceRevision metadata."
                )
            service = (manifest.get("services") or {}).get("numericRaster") or {}
            if not all(
                service.get(key)
                for key in ("baseUrl", "statisticsPath", "samplePath", "tilePath")
            ):
                raise ValueError(
                    f"Current-view layer {layer_id!r} requires the numericRaster service contract."
                )

    tree_layer_ids: set[str] = set()

    def visit(node: Mapping[str, Any]) -> None:
        layer_id = node.get("layerId")
        if layer_id:
            if layer_id not in layers:
                raise ValueError(f"Tree references missing layer {layer_id!r}.")
            if layer_id in tree_layer_ids:
                raise ValueError(f"Tree references layer {layer_id!r} more than once.")
            tree_layer_ids.add(str(layer_id))
        for child in node.get("children", []):
            visit(child)

    for root in tree:
        visit(root)

    if tree_layer_ids != set(layers):
        missing = sorted(set(layers) - tree_layer_ids)
        raise ValueError("Viewer layers missing from semantic tree: " + ", ".join(missing))

    interaction = manifest.get("interaction") or {}
    active_layer = interaction.get("activeLayerId")
    if active_layer and active_layer not in layers:
        raise ValueError(f"Active layer {active_layer!r} does not exist.")
    for layer_id in interaction.get("pinnedLayerIds", []):
        if layer_id not in layers:
            raise ValueError(f"Pinned layer {layer_id!r} does not exist.")

Publication Validation

ras2cng.publication

Validate an Example Library viewer bundle before public catalog admission.

validate_example_publication(viewer_manifest, archive_manifest=None, *, check_files=True, check_http_ranges=False)

Run the strict RAS Commander Example Library publication gate.

Source code in ras2cng/publication.py

def validate_example_publication(
    viewer_manifest: Path | Mapping[str, Any],
    archive_manifest: Path | Mapping[str, Any] | None = None,
    *,
    check_files: bool = True,
    check_http_ranges: bool = False,
) -> PublicationReport:
    """Run the strict RAS Commander Example Library publication gate."""

    manifest, manifest_path = _load_document(viewer_manifest)
    archive, archive_path = (
        _load_document(archive_manifest)
        if archive_manifest is not None
        else (None, None)
    )
    base_dir = manifest_path.parent if manifest_path else None
    report = PublicationReport(str(manifest_path or "<mapping>"))

    try:
        validate_manifest_v2(manifest)
    except Exception as error:
        report.add("error", "manifest.v2", str(error))

    resources = manifest.get("resources") if isinstance(manifest.get("resources"), Mapping) else {}
    layers = manifest.get("layers") if isinstance(manifest.get("layers"), Mapping) else {}
    completed_plan_ids = {
        str(plan.get("plan_id"))
        for plan in (archive or {}).get("results", [])
        if isinstance(plan, Mapping) and plan.get("completed") is True and plan.get("plan_id")
    }
    source_crs = (manifest.get("provenance") or {}).get("sourceCrs") or manifest.get("sourceCrs")
    if not source_crs:
        report.add("error", "project.crs", "The viewer manifest has no validated project CRS.")

    _validate_no_local_paths(manifest, report)
    _validate_resources(
        resources,
        report,
        base_dir=base_dir,
        check_files=check_files,
        check_http_ranges=check_http_ranges,
    )
    _validate_extent_color_service(manifest, resources, layers, report)

    geometry_layers = {
        layer_id: layer for layer_id, layer in layers.items()
        if layer.get("sourceKind") == "geometry"
    }
    extent_layers = {
        layer_id: layer for layer_id, layer in geometry_layers.items()
        if layer.get("role") == _MODEL_EXTENT_ROLE
    }
    geometry_ids = sorted(
        {
            str(layer.get("geometry"))
            for layer in geometry_layers.values()
            if layer.get("geometry")
        }
    )
    visible_geometry_ids = {
        str(layer.get("geometry"))
        for layer in geometry_layers.values()
        if layer.get("geometry") and layer.get("visible") is True
    }
    if geometry_ids and len(visible_geometry_ids) != 1:
        report.add(
            "error",
            "defaults.geometry",
            "Exactly one geometry must be enabled by default; enabled geometries: "
            + (", ".join(sorted(visible_geometry_ids)) or "none"),
        )
    if not extent_layers:
        report.add("error", "geometry.extent", "No API-derived Model Extents layer is published.")
    for geometry_id in geometry_ids:
        matches = [layer for layer in extent_layers.values() if layer.get("geometry") == geometry_id]
        if not matches:
            report.add(
                "error",
                "geometry.extent",
                f"Geometry {geometry_id} has no Model Extents layer.",
                geometry_id,
            )
    if extent_layers and not any(layer.get("visible") is True for layer in extent_layers.values()):
        report.add("error", "defaults.extent", "Model Extents must be enabled by default.")
    for layer_id, layer in extent_layers.items():
        if not _valid_wgs84_bounds(layer.get("bounds")):
            report.add("error", "geometry.extent-bounds", "Model Extents has invalid WGS84 bounds.", layer_id)

    plan_ids = sorted(
        {
            str(layer.get("plan"))
            for layer in layers.values()
            if layer.get("plan")
        }
    )
    raw_layers = {
        layer_id: layer for layer_id, layer in layers.items()
        if layer.get("sourceKind") == "raw-hdf"
    }
    stored_layers = {
        layer_id: layer for layer_id, layer in layers.items()
        if layer.get("sourceKind") == "stored-map"
    }
    calculated_layers = {
        layer_id: layer for layer_id, layer in layers.items()
        if layer.get("sourceKind") == "calculated"
    }
    terrain_layers = {
        layer_id: layer for layer_id, layer in layers.items()
        if layer.get("sourceKind") == "terrain"
    }
    for layer_id, layer in {**raw_layers, **stored_layers}.items():
        missing_titles = [
            key
            for key in ("planTitle", "geometryTitle")
            if not str(layer.get(key) or "").strip()
        ]
        if missing_titles:
            report.add(
                "error",
                "results.metadata",
                "Result layer is missing " + ", ".join(missing_titles) + ".",
                layer_id,
            )
        if layer.get("visible") is True:
            report.add(
                "error",
                "defaults.results",
                "Result layers must be disabled by default.",
                layer_id,
            )
    if not plan_ids:
        report.add("error", "results.plan", "No result plan is published.")
    admission_plan_ids = sorted(set(plan_ids) | completed_plan_ids)
    geometry_2d_ids = {
        str(layer.get("geometry"))
        for layer in geometry_layers.values()
        if layer.get("role") in _TWO_DIMENSIONAL_ROLES and layer.get("geometry")
    }
    archive_plans = {
        str(plan.get("plan_id")): plan
        for plan in (archive or {}).get("results", [])
        if isinstance(plan, Mapping) and plan.get("plan_id")
    }
    stored_map_exempt_plans: set[str] = set()
    for plan_id in admission_plan_ids:
        plan_raw = [layer for layer in raw_layers.values() if layer.get("plan") == plan_id]
        plan_stored = [layer for layer in stored_layers.values() if layer.get("plan") == plan_id]
        if not plan_raw:
            report.add(
                "error",
                "results.raw-hdf",
                "No raw HDF vector result layers are published.",
                plan_id,
            )
        expected_variables = _mappable_archive_variables(archive_plans.get(plan_id, {}))
        published_variables = {
            str((layer.get("provenance") or {}).get("variable") or "").casefold()
            for layer in plan_raw
            if (layer.get("provenance") or {}).get("variable")
        }
        for normalized_name, variable_name in expected_variables.items():
            if normalized_name not in published_variables:
                report.add(
                    "error",
                    "results.raw-variable",
                    f"Joinable raw HDF result variable {variable_name} is not published.",
                    f"{plan_id}:{variable_name}",
                )
        plan_geometry = str(archive_plans.get(plan_id, {}).get("geom_id") or "")
        stored_maps_applicable = bool(terrain_layers) or (
            plan_geometry in geometry_2d_ids
            if plan_geometry
            else bool(geometry_2d_ids)
        )
        if not plan_stored and stored_maps_applicable:
            report.add(
                "error",
                "results.stored-map",
                "No RASMapper Stored Map rasters are published.",
                plan_id,
            )
        elif not plan_stored:
            stored_map_exempt_plans.add(plan_id)
            report.add(
                "warning",
                "results.stored-map-not-applicable",
                "Pure 1D plan has no project terrain; continuous RASMapper Stored Map rasters are not applicable.",
                plan_id,
            )
        elif stored_maps_applicable:
            published_map_types = {
                map_type
                for layer in plan_stored
                if (
                    map_type := stored_map_type_key(
                        str((layer.get("provenance") or {}).get("mapType") or layer.get("role") or "")
                    )
                )
            }
            if any(
                _is_authoritative_depth_derived_boundary(layer)
                for layer in calculated_layers.values()
                if layer.get("plan") == plan_id
            ):
                published_map_types.add("inundation_boundary")
            missing_map_types = sorted(
                REQUIRED_STORED_MAP_TYPE_KEYS - published_map_types
            )
            if missing_map_types:
                report.add(
                    "error",
                    "results.stored-map-type",
                    "Complete Stored Map set is missing: "
                    + ", ".join(missing_map_types),
                    plan_id,
                )

    for layer_id, layer in raw_layers.items():
        query = layer.get("query") or {}
        provenance = layer.get("provenance") or {}
        if query.get("enabled") is not True:
            report.add("error", "results.raw-query", "Raw HDF layer is not queryable.", layer_id)
        if provenance.get("interpolationAuthority") != "none":
            report.add(
                "error",
                "results.raw-provenance",
                "Raw HDF values must explicitly declare no interpolation authority.",
                layer_id,
            )

    for layer_id, layer in stored_layers.items():
        query = layer.get("query") or {}
        provenance = layer.get("provenance") or {}
        numeric_id = query.get("numericResource")
        display_resource = resources.get(layer.get("resource")) or {}
        is_vector_stored_map = display_resource.get("type") == "vector-pmtiles"
        if (
            not is_vector_stored_map
            and (not numeric_id or (resources.get(numeric_id) or {}).get("type") != "cog")
        ):
            report.add(
                "error",
                "results.numeric-cog",
                "Stored Map has no authoritative numeric COG resource.",
                layer_id,
            )
        if is_vector_stored_map and query.get("enabled") is not True:
            report.add(
                "error",
                "results.vector-query",
                "Vector Stored Map is not queryable.",
                layer_id,
            )
        if provenance.get("interpolationAuthority") != "RASMapper/RasProcess":
            report.add(
                "error",
                "results.stored-provenance",
                "Stored Map must identify RASMapper/RasProcess as interpolation authority.",
                layer_id,
            )

    has_2d = any(layer.get("role") in _TWO_DIMENSIONAL_ROLES for layer in geometry_layers.values())
    if has_2d and not terrain_layers:
        report.add("error", "terrain.required", "A 2D model has no published terrain layer.")
    if has_2d and terrain_layers and not any(layer.get("visible") is True for layer in terrain_layers.values()):
        report.add("error", "defaults.terrain", "A 2D project terrain must be enabled by default.")
    for layer_id, layer in terrain_layers.items():
        query = layer.get("query") or {}
        numeric_id = query.get("numericResource")
        if not numeric_id or (resources.get(numeric_id) or {}).get("type") != "cog":
            report.add(
                "error",
                "terrain.numeric-cog",
                "Terrain has no associated queryable numeric COG.",
                layer_id,
            )

    basemap = next(
        (layer for layer in layers.values() if layer.get("role") == "basemap"),
        None,
    )
    if not basemap or basemap.get("visible") is not True:
        report.add("error", "defaults.basemap", "Hybrid satellite imagery must be enabled by default.")

    if archive is None:
        report.add(
            "error",
            "archive.required",
            "The archive manifest is required to verify successful plan completion.",
        )
    else:
        if not completed_plan_ids:
            report.add("error", "results.completed", "No successfully computed plan is recorded.")
        for plan_id in plan_ids:
            if plan_id not in completed_plan_ids:
                report.add(
                    "error",
                    "results.completed",
                    "Published result plan is not recorded as successfully computed.",
                    plan_id,
                )
        if terrain_layers:
            _validate_archive_terrain_policy(
                archive,
                report,
                base_dir=archive_path.parent if archive_path else None,
                check_files=check_files,
            )

    model_bounds = [layer.get("bounds") for layer in extent_layers.values() if _valid_wgs84_bounds(layer.get("bounds"))]
    numeric_stored_layers = {
        layer_id: layer
        for layer_id, layer in stored_layers.items()
        if (layer.get("query") or {}).get("numericResource")
    }
    for layer_id, layer in {**terrain_layers, **numeric_stored_layers}.items():
        numeric_id = (layer.get("query") or {}).get("numericResource")
        numeric = resources.get(numeric_id) or {}
        if not numeric.get("crs"):
            report.add("error", "raster.crs", "Numeric COG has no source CRS metadata.", layer_id)
        bounds = numeric.get("bounds")
        if not _valid_wgs84_bounds(bounds):
            report.add("error", "raster.bounds", "Numeric COG has invalid WGS84 bounds.", layer_id)
        elif model_bounds and not any(_bounds_intersect(bounds, extent) for extent in model_bounds):
            report.add("error", "raster.location", "Numeric COG does not intersect Model Extents.", layer_id)

    report.counts = {
        "resources": len(resources),
        "layers": len(layers),
        "geometries": len(geometry_ids),
        "plans": len(plan_ids),
        "completed_plans": len(completed_plan_ids),
        "raw_results": len(raw_layers),
        "stored_maps": len(stored_layers),
        "stored_map_exempt_plans": len(stored_map_exempt_plans),
        "terrains": len(terrain_layers),
    }
    return report