"""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 mapping 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 _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, ) -> 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) 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"}, "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}