"""Reusable geometric coverage of verified project nodes, independent of business analysis.""" import json from math import isfinite from typing import Any from pyproj import CRS, Transformer from shapely import MultiPoint, concave_hull, get_num_coordinates from shapely.geometry import GeometryCollection, mapping, shape from shapely.ops import transform, unary_union from app.native.wndb.gis.network_views import get_node_coverage_links, get_node_geographic_positions CONCAVE_RATIO = 0.25 BUFFER_METERS = 25.0 SIMPLIFY_METERS = 2.0 MAX_COVERAGE_NODES = 10000 def exclusive_coverage_features(features: list[dict[str, Any]]) -> list[dict[str, Any]]: """Clip display envelopes, with later groups taking precedence; retain input order. This changes only derived geometry, never node membership or model values. No post-clipping buffer or simplification may reintroduce an overlap. """ occupied = GeometryCollection() result = [] summaries = {} for feature in features: properties = feature["properties"] summaries.update(properties.get("coverage_group_summary", {})) for feature in features: properties = feature["properties"] summaries[properties["coverage_group"]] = { key: properties.get(key) for key in ("node_count", "mapped_node_count", "missing_node_count") } for feature in reversed(features): original = shape(feature["geometry"]) visible = original.difference(occupied) occupied = unary_union([occupied, original]) if visible.is_empty: continue if (visible.geom_type not in {"Polygon", "MultiPolygon"} or not visible.is_valid or get_num_coordinates(visible) > 20000): raise ValueError("Invalid exclusive derived coverage") result.append({ **feature, "properties": {**feature["properties"], "overlap_policy": "higher_priority_group_first"}, "geometry": mapping(visible), }) visible_groups = {feature["properties"]["coverage_group"] for feature in result} summaries = {key: {**value, "display_area_available": key in visible_groups} for key, value in summaries.items()} for feature in result: # Retain counts even when that group's entire display envelope was # clipped away. Overlay transport deliberately sends only GeoJSON. feature["properties"]["coverage_group_summary"] = summaries return list(reversed(result)) def _components(node_ids: list[str], links: list[dict[str, Any]], closed: set[str]) -> list[list[str]]: parent = {node_id: node_id for node_id in node_ids} def root(node: str) -> str: while parent[node] != node: parent[node] = parent[parent[node]] node = parent[node] return node for link in links: start, end = link["start_node_id"], link["end_node_id"] if link["id"] not in closed and link.get("status") != "CLOSED" and start in parent and end in parent: parent[root(start)] = root(end) groups: dict[str, list[str]] = {} for node_id in node_ids: groups.setdefault(root(node_id), []).append(node_id) return list(groups.values()) def build_node_coverage( groups: list[dict[str, Any]], positions: list[dict[str, Any]], links: list[dict[str, Any]], *, split_by_connectivity: bool = False, excluded_link_ids: list[str] | None = None, concave_ratio: float = CONCAVE_RATIO, buffer_m: float = BUFFER_METERS, exclusive_groups: bool = False, ) -> dict[str, Any]: """Keep caller-defined comparison groups and derive geometry without assigning business meaning.""" coordinates = {} for row in positions: longitude, latitude = row.get("longitude"), row.get("latitude") if (isinstance(longitude, (int, float)) and not isinstance(longitude, bool) and isinstance(latitude, (int, float)) and not isinstance(latitude, bool) and isfinite(longitude) and isfinite(latitude) and -180 <= longitude <= 180 and -90 < latitude < 90): coordinates[str(row["id"])] = (longitude, latitude) features, groups_metadata = [], {} closed = set(excluded_link_ids or []) simplify_m = min(SIMPLIFY_METERS, buffer_m / 4) for group in groups: group_id, label = group["group_id"], group["title"] ids = sorted(set(group["node_ids"])) mapped = [node for node in ids if node in coordinates] missing = [node for node in ids if node not in coordinates] metadata = {"node_ids": ids, "mapped_node_count": len(mapped), "missing_node_ids": missing} groups_metadata[group_id] = {"title": label, **metadata} metadata = groups_metadata[group_id] if not mapped: metadata["component_count"] = 0 continue points = [coordinates[node] for node in mapped] longitude = sum(point[0] for point in points) / len(points) latitude = sum(point[1] for point in points) / len(points) local = CRS.from_proj4( f"+proj=aeqd +lat_0={latitude} +lon_0={longitude} +datum=WGS84 +units=m" ) project = Transformer.from_crs(4326, local, always_xy=True).transform unproject = Transformer.from_crs(local, 4326, always_xy=True).transform envelopes = [] components = _components(ids, links, closed) if split_by_connectivity else [ids] for nodes in components: component_points = [coordinates[node] for node in nodes if node in coordinates] if not component_points: continue cloud = transform(project, MultiPoint(component_points)) # A point or collinear cloud remains drawable through the same metric buffer. envelope = concave_hull(cloud, ratio=concave_ratio).buffer(buffer_m) envelopes.append(envelope) metadata["component_count"] = len(envelopes) geometry = unary_union(envelopes).simplify(simplify_m, preserve_topology=True) if geometry.is_empty or not geometry.is_valid or get_num_coordinates(geometry) > 20000: raise ValueError("Invalid derived coverage") feature = { "type": "Feature", "id": group_id, "properties": { "coverage_group": group_id, "label": label, "derived": True, "display_only": True, "node_count": len(ids), "mapped_node_count": len(mapped), "missing_node_count": len(missing), "component_count": len(envelopes), }, "geometry": mapping(transform(unproject, geometry)), } features.append(feature) if exclusive_groups: features = exclusive_coverage_features(features) visible_groups = {feature["properties"]["coverage_group"] for feature in features} for group_id, metadata in groups_metadata.items(): metadata["display_area_available"] = group_id in visible_groups return { "status": "partial" if any(group["missing_node_ids"] for group in groups_metadata.values()) else "ready" if features else "empty", "method": "concave_hull_buffer", "parameters": {"concave_ratio": concave_ratio, "buffer_m": buffer_m, "split_by_connectivity": split_by_connectivity, "excluded_link_ids": sorted(closed), "simplify_m": simplify_m, "projection": "local_azimuthal_equidistant", "exclusive_groups": exclusive_groups}, "output_crs": "EPSG:4326", "display_only": True, "position_source": "gis.node_geometries", "boundary": "节点空间覆盖示意；不推断故障、污染、停水、水力影响或用户服务边界。", "groups": groups_metadata, "map_payload": {"type": "FeatureCollection", "features": features}, } def generate_node_coverage(network: str, groups: list[dict[str, Any]], **parameters: Any) -> dict[str, Any]: node_ids = sorted({node for group in groups for node in group["node_ids"]}) if len(node_ids) > MAX_COVERAGE_NODES: raise ValueError("Node coverage supports at most 10000 distinct nodes") positions = get_node_geographic_positions(network, node_ids) links = get_node_coverage_links(network, node_ids) if parameters.get("split_by_connectivity") else [] coverage = build_node_coverage(groups, positions, links, **parameters) payload = coverage.pop("map_payload") if len(json.dumps(payload, ensure_ascii=False).encode("utf-8")) > 400 * 1024: raise ValueError("Coverage exceeds presentation budget; narrow the node selection") return {"coverage": coverage, "map_payload": payload}