fix: scope S-52 conservation audit to production footprint

docs/CM95_S52_语义守恒与运行时能力审计标准_v1.md

@@ -138,6 +138,18 @@ Source→PBF 审计必须绑定具体 Production Plan：
 
 不得把 AOI 内但不属于该 plan scale ladder 的 Source record 计入 expected set。
 
+仅按 AOI bbox + scale ladder 过滤仍然不够。Source→PBF expected scope 必须重放 Production 的实际可见范围语义：
+
+- 按每个计划 tile 计算与生产一致的 `buffer_tile_pixels` expanded bounds；
+- 按 zoom 的 scale ladder 顺序计算 coverage；
+- finer scale coverage 先占用区域，coarser scale 只允许进入 `coverage - occupied` 的剩余可见区；
+- AOI 边缘允许因 tile buffer 合法带入外侧 Source component；
+- 同一 relation group 必须在“该 tile、该 native_scale、该 visible region”实际 selected 的成员集合上重新选择 base 与 LIGHTS；
+- 不得先在全局 relation group 上固定一个 base，再套到所有 tile；
+- expected tile occurrence set 必须与 PBF actual occurrence set 完全一致。
+
+因此 audit 的 scope truth 是 **Production Plan + policy + coverage + tile buffer**，而不是裸 AOI bbox。AOI bbox 只用于说明工程范围，不能单独定义 expected set。
+
 审计 artifact 必须记录 plan_id 和 PBF directory。Finalize 必须验证 artifact 的 plan_id / PBF directory 与当前 release 一致，禁止拿其他 release 的 PASS 报告复用。
 ## 5. Multi-LIGHTS 专项规则
 

src/cm95pbf/release.py

@@ -97,6 +97,7 @@ def _external_audit_contract(
         "actual_presentation_records",
         "signature_mismatch_count",
         "occurrence_drift_count",
+        "occurrence_scope_mismatch_count",
         "required_runtime_kinds",
         "missing_runtime_kinds",
         "noncompliant_runtime_kinds",

tools/cm95_s52_presentation_conservation_audit.py

@@ -8,14 +8,24 @@ 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 connect
+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,
@@ -87,7 +97,7 @@ def json_attrs(value: Any) -> dict[str, Any]:
     return dict(obj) if isinstance(obj, dict) else {}
 
 
-def source_expected(
+def _source_expected_legacy(
     dataset_id: int,
     bbox: tuple[float, float, float, float],
     allowed_scales: set[int],
@@ -136,18 +146,18 @@ def source_expected(
             )
             rows = [dict(row) for row in cur.fetchall()]
 
-    groups: dict[tuple[int, str], list[dict[str, Any]]] = defaultdict(list)
+    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["cell_id"]), str(relation))].append(row)
+        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 (cell_id, relation), members in groups.items():
+        for (native_scale, cell_id, relation), members in groups.items():
             lights = [row for row in members if row["object_class"] == "LIGHTS"]
             bases = [
                 row
@@ -225,7 +235,7 @@ def source_expected(
                 "relation20": relation,
                 "base_source_record_index": int(base["source_record_index"]),
                 "base_object_class": str(base["object_class"]),
-                "native_scale": int(base["native_scale"]),
+                "native_scale": native_scale,
                 "presentation_count": len(presentations),
                 "signatures": sorted(signatures),
                 "presentations": sorted(
@@ -236,6 +246,515 @@ def source_expected(
     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)
+            )
+
+        pieces = []
+        for cell in cells:
+            cid = int(cell["id"])
+            if by_cell.get(cid):
+                pieces.append(unary_union(by_cell[cid]))
+                continue
+            minx, miny = lonlat_to_mercator(
+                float(cell["bbox_min_lon"]),
+                float(cell["bbox_min_lat"]),
+            )
+            maxx, maxy = lonlat_to_mercator(
+                float(cell["bbox_max_lon"]),
+                float(cell["bbox_max_lat"]),
+            )
+            pieces.append(box(minx, miny, maxx, maxy))
+
+        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]]]:
@@ -293,15 +812,21 @@ def pbf_actual(
                     "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
 
 
@@ -315,6 +840,7 @@ def compare_contract(
     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]
@@ -333,6 +859,22 @@ def compare_contract(
                     "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(
                 {
@@ -355,6 +897,8 @@ def compare_contract(
         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",
@@ -373,6 +917,7 @@ def compare_contract(
         "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,
@@ -384,6 +929,7 @@ def compare_contract(
         ],
         "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],
     }
@@ -407,8 +953,7 @@ def main() -> None:
     )
     expected, source_failures = source_expected(
         args.dataset_id,
-        bbox,
+        plan,
-        allowed_scales,
         Path(args.manifest),
         registry,
     )
@@ -424,6 +969,8 @@ def main() -> None:
         "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)),