#!/usr/bin/env python3 from __future__ import annotations import argparse from collections import defaultdict import json from pathlib import Path from typing import Any import mapbox_vector_tile import yaml from shapely import wkb from shapely.geometry import Point, box from shapely.ops import transform, unary_union from cm93db.config import load_config from cm93db.database import as_wkb_sql, connect, mysql_compat_mode_for_connection from cm95pbf.identity import ( FidManifest, ROLE_PRIMARY, feature_source_key, ) from cm95pbf.policy import ( expanded_mercator_bounds, lonlat_to_mercator, mercator_to_lonlat, xyz_tile_bounds, ) from cm95pbf.presentation_audit import ( LIGHT_PRESENTATION_SIGNATURE_FIELDS, canonical_light_presentation, light_presentation_signature, parse_light_presentations, ) from cm95pbf.s52 import S52LookupRegistry SCHEMA = "navsea.cm95.s52-presentation-conservation-audit/v1" NAV_COMPONENT_CLASSES = ( "LIGHTS", "BCNLAT", "BCNSPP", "BOYSPP", "RTPBCN", "TOWERS", "_extgn", ) BASE_PRECEDENCE = { "BCNLAT": 0, "BCNSPP": 1, "BOYSPP": 2, "RTPBCN": 3, "TOWERS": 4, "_extgn": 5, } def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser( description=( "Audit multi-LIGHTS presentation conservation directly from " "CM95 source rows to decoded delivery PBF." ) ) parser.add_argument("--dataset-id", type=int, required=True) parser.add_argument("--plan", required=True) parser.add_argument("--aoi-bbox", required=True, help="west,south,east,north") parser.add_argument("--manifest", required=True) parser.add_argument("--lookups", required=True) parser.add_argument("--pbf-dir", required=True) parser.add_argument("--output", required=True) parser.add_argument("--point-table", default="Paper") parser.add_argument("--boundary-table", default="Plain") parser.add_argument("--mariner-mode", default="runtime") return parser.parse_args() def parse_bbox(raw: str) -> tuple[float, float, float, float]: values = tuple(float(item.strip()) for item in raw.split(",")) if len(values) != 4: raise SystemExit("--aoi-bbox must be west,south,east,north") west, south, east, north = values if west >= east or south >= north: raise SystemExit("invalid --aoi-bbox") return west, south, east, north def json_attrs(value: Any) -> dict[str, Any]: if isinstance(value, dict): return dict(value) if isinstance(value, (bytes, bytearray, memoryview)): value = bytes(value).decode("utf-8") try: obj = json.loads(value or "{}") except Exception: return {} return dict(obj) if isinstance(obj, dict) else {} def _source_expected_legacy( dataset_id: int, bbox: tuple[float, float, float, float], allowed_scales: set[int], manifest_path: Path, registry: S52LookupRegistry, ) -> tuple[dict[int, dict[str, Any]], list[dict[str, Any]]]: west, south, east, north = bbox cfg = load_config() with connect(cfg) as conn: with conn.cursor() as cur: placeholders = ",".join(["%s"] * len(NAV_COMPONENT_CLASSES)) cur.execute( f""" SELECT f.id, f.cell_id, c.source_relative_path, f.source_record_index, f.object_class, f.geometry_type, f.native_scale, f.attributes_json, f.bbox_min_lon, f.bbox_min_lat FROM chart_feature f JOIN chart_cell c ON c.id=f.cell_id WHERE f.dataset_id=%s AND f.geometry_type='POINT' AND f.object_class IN ({placeholders}) AND f.native_scale IN ({",".join(["%s"] * len(allowed_scales))}) AND f.bbox_max_lon >= %s AND f.bbox_min_lon <= %s AND f.bbox_max_lat >= %s AND f.bbox_min_lat <= %s ORDER BY f.cell_id, f.source_record_index """, ( dataset_id, *NAV_COMPONENT_CLASSES, *sorted(allowed_scales), west, east, south, north, ), ) rows = [dict(row) for row in cur.fetchall()] groups: dict[tuple[int, int, str], list[dict[str, Any]]] = defaultdict(list) for row in rows: row["attrs"] = json_attrs(row.pop("attributes_json", {})) relation = row["attrs"].get("_cm93_relation20_raw") if relation is None: continue groups[(int(row["native_scale"]), int(row["cell_id"]), str(relation))].append(row) expected: dict[int, dict[str, Any]] = {} source_failures: list[dict[str, Any]] = [] with FidManifest(manifest_path) as manifest: for (native_scale, cell_id, relation), members in groups.items(): lights = [row for row in members if row["object_class"] == "LIGHTS"] bases = [ row for row in members if row["object_class"] in BASE_PRECEDENCE ] if len(lights) < 2 or not bases: continue base = sorted( bases, key=lambda row: ( BASE_PRECEDENCE.get(row["object_class"], 99), int(row["source_record_index"]), ), )[0] source_key = feature_source_key( dataset_id, str(base["source_relative_path"]), int(base["source_record_index"]), str(base["geometry_type"]), ) try: fid = manifest.mvt_id(source_key, ROLE_PRIMARY) except Exception as exc: source_failures.append( { "cell_id": cell_id, "relation20": relation, "base_source_record_index": base["source_record_index"], "failure": f"manifest:{type(exc).__name__}:{exc}", } ) continue presentations: list[dict[str, Any]] = [] for light in sorted( lights, key=lambda row: int(row["source_record_index"]) ): props = registry.properties( str(light["object_class"]), str(light["geometry_type"]), dict(light["attrs"]), ) presentation = { key: value for key, value in props.items() if str(key).startswith("s52_") } presentation["source_record_index"] = int( light["source_record_index"] ) if presentation.get("s52_csp_resolved") != "LIGHTS05": source_failures.append( { "feature_id": fid, "cell_id": cell_id, "relation20": relation, "source_record_index": light["source_record_index"], "failure": "source LIGHTS did not resolve to LIGHTS05", } ) presentations.append(presentation) canonical = [ canonical_light_presentation(item) for item in presentations ] signatures = [ light_presentation_signature(item) for item in presentations ] expected[fid] = { "feature_id": fid, "cell_id": cell_id, "relation20": relation, "base_source_record_index": int(base["source_record_index"]), "base_object_class": str(base["object_class"]), "native_scale": native_scale, "presentation_count": len(presentations), "signatures": sorted(signatures), "presentations": sorted( canonical, key=lambda item: str(item.get("source_record_index")), ), } return expected, source_failures def _decode_wkb(value: Any): if value is None: return None if isinstance(value, memoryview): value = value.tobytes() elif isinstance(value, bytearray): value = bytes(value) return wkb.loads(value) def _to_mercator_geom(geom): return transform( lambda x, y, z=None: lonlat_to_mercator(x, y), geom, ) def _tile_xy_for_lonlat(lon: float, lat: float, z: int) -> tuple[int, int]: import math n = 1 << z x = int((lon + 180.0) / 360.0 * n) lat = max(min(lat, 85.05112878), -85.05112878) lat_rad = math.radians(lat) y = int( ( 1.0 - math.asinh(math.tan(lat_rad)) / math.pi ) / 2.0 * n ) return x, y def _plan_tile_scope( plan: dict[str, Any], policy: dict[str, Any], ) -> tuple[ dict[int, dict[tuple[int, int], tuple[float, float, float, float]]], tuple[float, float, float, float], ]: buffer_px = float(policy.get("buffer_tile_pixels", 8)) by_zoom: dict[ int, dict[tuple[int, int], tuple[float, float, float, float]], ] = defaultdict(dict) broad: list[float] | None = None for job in plan.get("jobs") or []: z = int(job["z"]) x = int(job["x"]) y = int(job["y"]) tile = xyz_tile_bounds(z, x, y) bounds = expanded_mercator_bounds(tile, buffer_px) by_zoom[z][(x, y)] = bounds if broad is None: broad = list(bounds) else: broad[0] = min(broad[0], bounds[0]) broad[1] = min(broad[1], bounds[1]) broad[2] = max(broad[2], bounds[2]) broad[3] = max(broad[3], bounds[3]) if broad is None: raise RuntimeError("production plan contains no jobs") return dict(by_zoom), tuple(broad) def _load_coverage_unions( conn, dataset_id: int, scales: set[int], broad_bbox_wgs: tuple[float, float, float, float], broad_box_m, ) -> dict[int, Any]: west, south, east, north = broad_bbox_wgs db_mode = mysql_compat_mode_for_connection(conn) coverage_wkb_sql = as_wkb_sql(db_mode, "geom") result: dict[int, Any] = {} for scale in sorted(scales): with conn.cursor() as cur: cur.execute( """ SELECT id, bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat FROM chart_cell WHERE dataset_id=%s AND native_scale=%s AND bbox_max_lon >= %s AND bbox_min_lon <= %s AND bbox_max_lat >= %s AND bbox_min_lat <= %s """, ( dataset_id, scale, west, east, south, north, ), ) cells = [dict(row) for row in cur.fetchall()] if not cells: continue ids = [int(row["id"]) for row in cells] placeholders = ",".join(["%s"] * len(ids)) with conn.cursor() as cur: cur.execute( f""" SELECT cell_id, {coverage_wkb_sql} AS geom_wkb FROM chart_coverage WHERE dataset_id=%s AND cell_id IN ({placeholders}) """, (dataset_id, *ids), ) coverage_rows = [dict(row) for row in cur.fetchall()] by_cell: dict[int, list[Any]] = defaultdict(list) for row in coverage_rows: geom = _decode_wkb(row.get("geom_wkb")) if geom is None or geom.is_empty: continue by_cell[int(row["cell_id"])].append( _to_mercator_geom(geom) ) # Keep the independent audit scope identical to OpenCPN/Builder # coverage semantics: only explicit M_COVR rows define valid coverage. # Never substitute the chart_cell index bbox for a missing M_COVR. pieces = [ unary_union(by_cell[cid]) for cid in sorted(by_cell) if by_cell[cid] ] if pieces: result[scale] = unary_union(pieces).intersection( broad_box_m ) return result def _visible_regions( plan: dict[str, Any], coverage_by_scale: dict[int, Any], ) -> dict[tuple[int, int], Any]: visible: dict[tuple[int, int], Any] = {} for z in plan.get("zooms") or []: z = int(z) ladder = [ int(value) for value in ( (plan.get("source_stats") or {}) .get(str(z), {}) .get("scale_ladder") or [] ) ] occupied = None for scale in ladder: coverage = coverage_by_scale.get(scale) if coverage is None or coverage.is_empty: visible[(z, scale)] = None continue region = ( coverage if occupied is None else coverage.difference(occupied) ) visible[(z, scale)] = region occupied = ( coverage if occupied is None else unary_union([occupied, coverage]) ) return visible def source_expected( dataset_id: int, plan: dict[str, Any], manifest_path: Path, registry: S52LookupRegistry, ) -> tuple[dict[int, dict[str, Any]], list[dict[str, Any]]]: policy_path = Path(str(plan.get("policy_path") or "")) if not policy_path.is_file(): raise RuntimeError(f"production policy missing: {policy_path}") policy = yaml.safe_load(policy_path.read_text(encoding="utf-8")) or {} jobs_by_zoom, broad_m = _plan_tile_scope(plan, policy) broad_box_m = box(*broad_m) west, south = mercator_to_lonlat(broad_m[0], broad_m[1]) east, north = mercator_to_lonlat(broad_m[2], broad_m[3]) broad_bbox_wgs = (west, south, east, north) allowed_scales = { int(scale) for stats in (plan.get("source_stats") or {}).values() for scale in (stats.get("scale_ladder") or []) } ladder_by_zoom = { int(z): [ int(scale) for scale in ( (plan.get("source_stats") or {}) .get(str(z), {}) .get("scale_ladder") or [] ) ] for z in (plan.get("zooms") or []) } buffer_px = float(policy.get("buffer_tile_pixels", 8)) neighbor_radius = max(1, int(buffer_px // 256) + 1) cfg = load_config() with connect(cfg) as conn: coverage_by_scale = _load_coverage_unions( conn, dataset_id, allowed_scales, broad_bbox_wgs, broad_box_m, ) visible_by_zoom_scale = _visible_regions( plan, coverage_by_scale, ) db_mode = mysql_compat_mode_for_connection(conn) point_wkb_sql = as_wkb_sql(db_mode, "g.geom") with conn.cursor() as cur: class_placeholders = ",".join( ["%s"] * len(NAV_COMPONENT_CLASSES) ) scale_placeholders = ",".join( ["%s"] * len(allowed_scales) ) cur.execute( f""" SELECT f.id, f.cell_id, c.source_relative_path, f.source_record_index, f.object_class, f.geometry_type, f.native_scale, f.attributes_json, {point_wkb_sql} AS geom_wkb FROM chart_feature f JOIN chart_cell c ON c.id=f.cell_id JOIN chart_point_geometry g ON g.feature_id=f.id WHERE f.dataset_id=%s AND f.geometry_type='POINT' AND f.object_class IN ({class_placeholders}) AND f.native_scale IN ({scale_placeholders}) AND f.bbox_max_lon >= %s AND f.bbox_min_lon <= %s AND f.bbox_max_lat >= %s AND f.bbox_min_lat <= %s ORDER BY f.native_scale, f.cell_id, f.source_record_index """, ( dataset_id, *NAV_COMPONENT_CLASSES, *sorted(allowed_scales), west, east, south, north, ), ) rows = [dict(row) for row in cur.fetchall()] source_failures: list[dict[str, Any]] = [] groups: dict[ tuple[int, int, str], list[dict[str, Any]], ] = defaultdict(list) for row in rows: row["attrs"] = json_attrs(row.pop("attributes_json", {})) relation = row["attrs"].get("_cm93_relation20_raw") if relation is None: continue geom = _decode_wkb(row.pop("geom_wkb", None)) if geom is None or geom.is_empty or geom.geom_type != "Point": source_failures.append( { "feature_id": int(row["id"]), "failure": "source_navigation_component_not_point", } ) continue row["geom_wgs84"] = geom row["geom_m"] = _to_mercator_geom(geom) row["visible_tiles"] = set() scale = int(row["native_scale"]) lon = float(geom.x) lat = float(geom.y) px = float(row["geom_m"].x) py = float(row["geom_m"].y) for z, ladder in ladder_by_zoom.items(): if scale not in ladder: continue visible_region = visible_by_zoom_scale.get((z, scale)) if ( visible_region is None or visible_region.is_empty or not visible_region.intersects(row["geom_m"]) ): continue x0, y0 = _tile_xy_for_lonlat(lon, lat, z) planned = jobs_by_zoom.get(z) or {} for dx in range(-neighbor_radius, neighbor_radius + 1): for dy in range(-neighbor_radius, neighbor_radius + 1): x = x0 + dx y = y0 + dy bounds = planned.get((x, y)) if bounds is None: continue if ( bounds[0] <= px <= bounds[2] and bounds[1] <= py <= bounds[3] ): row["visible_tiles"].add(f"{z}/{x}/{y}.pbf") groups[ ( scale, int(row["cell_id"]), str(relation), ) ].append(row) expected: dict[int, dict[str, Any]] = {} with FidManifest(manifest_path) as manifest: for (native_scale, cell_id, relation), members in groups.items(): selected_by_tile: dict[str, list[dict[str, Any]]] = defaultdict(list) for member in members: for tile_rel in member["visible_tiles"]: selected_by_tile[tile_rel].append(member) for tile_rel, selected in sorted(selected_by_tile.items()): lights = [ row for row in selected if row["object_class"] == "LIGHTS" ] bases = [ row for row in selected if row["object_class"] in BASE_PRECEDENCE ] if len(lights) < 2 or not bases: continue base = sorted( bases, key=lambda row: ( BASE_PRECEDENCE.get(row["object_class"], 99), int(row["source_record_index"]), ), )[0] source_key = feature_source_key( dataset_id, str(base["source_relative_path"]), int(base["source_record_index"]), str(base["geometry_type"]), ) try: fid = manifest.mvt_id(source_key, ROLE_PRIMARY) except Exception as exc: source_failures.append( { "tile": tile_rel, "cell_id": cell_id, "relation20": relation, "base_source_record_index": base[ "source_record_index" ], "failure": ( f"manifest:{type(exc).__name__}:{exc}" ), } ) continue presentations: list[dict[str, Any]] = [] for light in sorted( lights, key=lambda row: int(row["source_record_index"]), ): props = registry.properties( str(light["object_class"]), str(light["geometry_type"]), dict(light["attrs"]), ) presentation = { key: value for key, value in props.items() if str(key).startswith("s52_") } presentation["source_record_index"] = int( light["source_record_index"] ) if ( presentation.get("s52_csp_resolved") != "LIGHTS05" ): source_failures.append( { "feature_id": fid, "tile": tile_rel, "cell_id": cell_id, "relation20": relation, "source_record_index": light[ "source_record_index" ], "failure": ( "source LIGHTS did not resolve " "to LIGHTS05" ), } ) presentations.append(presentation) canonical = [ canonical_light_presentation(item) for item in presentations ] signatures = sorted( light_presentation_signature(item) for item in presentations ) canonical = sorted( canonical, key=lambda item: str( item.get("source_record_index") ), ) state = expected.setdefault( fid, { "feature_id": fid, "cell_id": cell_id, "relation20": relation, "base_source_record_index": int( base["source_record_index"] ), "base_object_class": str(base["object_class"]), "native_scale": native_scale, "presentation_count": len(presentations), "signatures": signatures, "presentations": canonical, "occurrences": [], }, ) if state["signatures"] != signatures: source_failures.append( { "feature_id": fid, "tile": tile_rel, "failure": ( "source_expected_occurrence_drift" ), "canonical_signatures": state[ "signatures" ], "tile_signatures": signatures, } ) state["occurrences"].append(tile_rel) for state in expected.values(): state["occurrences"] = sorted( set(state.get("occurrences") or []) ) state["occurrence_count"] = len(state["occurrences"]) state["occurrence_examples"] = state["occurrences"][:8] return expected, source_failures def pbf_actual( pbf_dir: Path, ) -> tuple[dict[int, dict[str, Any]], list[dict[str, Any]]]: actual: dict[int, dict[str, Any]] = {} failures: list[dict[str, Any]] = [] for path in sorted(pbf_dir.rglob("*.pbf")): rel = path.relative_to(pbf_dir).as_posix() decoded = mapbox_vector_tile.decode(path.read_bytes()) layer = decoded.get("navigation_marks") or {} for feature in layer.get("features") or []: fid = feature.get("id") props = feature.get("properties") or {} raw = props.get("s52_light_presentations") if fid is None or raw in (None, ""): continue try: presentations = parse_light_presentations(raw) except Exception as exc: failures.append( { "pbf": rel, "feature_id": fid, "failure": ( "invalid s52_light_presentations:" f"{type(exc).__name__}:{exc}" ), } ) continue canonical = [ canonical_light_presentation(item) for item in presentations ] signatures = sorted( light_presentation_signature(item) for item in presentations ) occurrence = { "pbf": rel, "presentation_count": len(presentations), "signatures": signatures, "presentations": sorted( canonical, key=lambda item: str(item.get("source_record_index")), ), } fid_int = int(fid) state = actual.setdefault( fid_int, { "feature_id": fid_int, "presentation_count": len(presentations), "signatures": signatures, "presentations": occurrence["presentations"], "occurrence_count": 0, "occurrences": [], "occurrence_examples": [], "occurrence_drift": [], }, ) state["occurrence_count"] += 1 state["occurrences"].append(rel) if len(state["occurrence_examples"]) < 8: state["occurrence_examples"].append(rel) if state["signatures"] != signatures: state["occurrence_drift"].append(occurrence) for state in actual.values(): state["occurrences"] = sorted(set(state.get("occurrences") or [])) state["occurrence_count"] = len(state["occurrences"]) state["occurrence_examples"] = state["occurrences"][:8] return actual, failures def compare_contract( expected: dict[int, dict[str, Any]], actual: dict[int, dict[str, Any]], source_failures: list[dict[str, Any]], pbf_failures: list[dict[str, Any]], ) -> dict[str, Any]: missing = sorted(set(expected) - set(actual)) unexpected = sorted(set(actual) - set(expected)) signature_mismatches: list[dict[str, Any]] = [] occurrence_drift: list[dict[str, Any]] = [] occurrence_scope_mismatches: list[dict[str, Any]] = [] for fid in sorted(set(expected) & set(actual)): exp = expected[fid] got = actual[fid] if exp["signatures"] != got["signatures"]: signature_mismatches.append( { "feature_id": fid, "cell_id": exp["cell_id"], "relation20": exp["relation20"], "expected_count": exp["presentation_count"], "actual_count": got["presentation_count"], "expected_signatures": exp["signatures"], "actual_signatures": got["signatures"], "expected_presentations": exp["presentations"], "actual_presentations": got["presentations"], } ) expected_occurrences = set(exp.get("occurrences") or []) actual_occurrences = set(got.get("occurrences") or []) if expected_occurrences != actual_occurrences: occurrence_scope_mismatches.append( { "feature_id": fid, "expected_occurrence_count": len(expected_occurrences), "actual_occurrence_count": len(actual_occurrences), "missing_occurrences": sorted( expected_occurrences - actual_occurrences )[:100], "unexpected_occurrences": sorted( actual_occurrences - expected_occurrences )[:100], } ) if got.get("occurrence_drift"): occurrence_drift.append( { "feature_id": fid, "canonical_signatures": got["signatures"], "drift_examples": got["occurrence_drift"][:8], } ) failures: list[str] = [] if source_failures: failures.append("source_resolution_failures") if pbf_failures: failures.append("pbf_parse_failures") if missing: failures.append("missing_multi_light_groups") if unexpected: failures.append("unexpected_multi_light_groups") if signature_mismatches: failures.append("presentation_signature_mismatch") if occurrence_drift: failures.append("tile_occurrence_presentation_drift") if occurrence_scope_mismatches: failures.append("tile_occurrence_scope_mismatch") return { "status": "PASS" if not failures else "FAIL", "expected_multi_light_groups": len(expected), "actual_multi_light_groups": len(actual), "expected_presentation_records": sum( int(item["presentation_count"]) for item in expected.values() ), "actual_presentation_records": sum( int(item["presentation_count"]) for item in actual.values() ), "matched_group_count": ( len(set(expected) & set(actual)) - len(signature_mismatches) ), "missing_group_count": len(missing), "unexpected_group_count": len(unexpected), "signature_mismatch_count": len(signature_mismatches), "occurrence_drift_count": len(occurrence_drift), "occurrence_scope_mismatch_count": len(occurrence_scope_mismatches), "source_resolution_failure_count": len(source_failures), "pbf_parse_failure_count": len(pbf_failures), "failures": failures, "missing_groups": [ expected[fid] for fid in missing[:100] ], "unexpected_groups": [ actual[fid] for fid in unexpected[:100] ], "signature_mismatches": signature_mismatches[:100], "occurrence_drift": occurrence_drift[:100], "occurrence_scope_mismatches": occurrence_scope_mismatches[:100], "source_resolution_failures": source_failures[:100], "pbf_parse_failures": pbf_failures[:100], } def main() -> None: args = parse_args() bbox = parse_bbox(args.aoi_bbox) plan = json.loads(Path(args.plan).read_text(encoding="utf-8")) allowed_scales: set[int] = set() for item in (plan.get("source_stats") or {}).values(): for scale in item.get("scale_ladder") or []: allowed_scales.add(int(scale)) if not allowed_scales: raise SystemExit("plan contains no source scale ladder") registry = S52LookupRegistry( Path(args.lookups), point_table=args.point_table, boundary_table=args.boundary_table, mariner_mode=args.mariner_mode, ) expected, source_failures = source_expected( args.dataset_id, plan, Path(args.manifest), registry, ) actual, pbf_failures = pbf_actual(Path(args.pbf_dir)) contract = compare_contract( expected, actual, source_failures, pbf_failures, ) result = { "schema": SCHEMA, "dataset_id": args.dataset_id, "plan": str(Path(args.plan)), "plan_id": plan.get("plan_id"), "scope_mode": "production-plan-visible-footprint-v1", "policy": str(plan.get("policy_path") or ""), "allowed_native_scales": sorted(allowed_scales), "aoi_bbox_wgs84": list(bbox), "manifest": str(Path(args.manifest)), "lookups": str(Path(args.lookups)), "pbf_dir": str(Path(args.pbf_dir)), "signature_fields": list(LIGHT_PRESENTATION_SIGNATURE_FIELDS), **contract, } output = Path(args.output) output.parent.mkdir(parents=True, exist_ok=True) output.write_text( json.dumps(result, ensure_ascii=False, indent=2), encoding="utf-8", ) print(json.dumps(result, ensure_ascii=False, indent=2)) if result["status"] != "PASS": raise SystemExit(2) if __name__ == "__main__": main()