178 lines
8.6 KiB
Python
178 lines
8.6 KiB
Python
"""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}
|