299 lines
9.8 KiB
Python
299 lines
9.8 KiB
Python
import platform
|
|
import math
|
|
from typing import Any
|
|
import pyclipper
|
|
from ..core.database import read, try_read, read_all
|
|
from .network_views import get_boundary_link_ids, get_topology_rows
|
|
|
|
|
|
def from_postgis_polygon(polygon: str) -> list[tuple[float, float]]:
|
|
boundary = polygon.lower().removeprefix('polygon((').removesuffix('))').split(',')
|
|
xys = []
|
|
for pt in boundary:
|
|
xy = pt.split(' ')
|
|
xys.append((float(xy[0]), float(xy[1])))
|
|
return xys
|
|
|
|
|
|
def to_postgis_polygon(boundary: list[tuple[float, float]]) -> str:
|
|
polygon = ''
|
|
for pt in boundary:
|
|
polygon += f'{pt[0]} {pt[1]},'
|
|
return str(f'polygon(({polygon[:-1]}))')
|
|
|
|
|
|
def to_postgis_linestring(boundary: list[tuple[float, float]]) -> str:
|
|
line = ''
|
|
for pt in boundary:
|
|
line += f'{pt[0]} {pt[1]},'
|
|
return str(f'linestring({line[:-1]})')
|
|
|
|
|
|
def get_nodes_in_boundary(name: str, boundary: list[tuple[float, float]]) -> list[str]:
|
|
rows = read_all(
|
|
name,
|
|
"select node_id from gis.node_geometries "
|
|
"where st_intersects(geom, st_geomfromtext(%s, 900914)) order by node_id",
|
|
(to_postgis_polygon(boundary),),
|
|
)
|
|
return [str(row["node_id"]) for row in rows]
|
|
|
|
|
|
def _get_links_on_boundary(name: str, nodes: list[str]) -> list[str]:
|
|
return get_boundary_link_ids(name, nodes)
|
|
|
|
|
|
def get_nodes_in_region(name: str, region_id: str) -> list[str]:
|
|
stored = read_all(
|
|
name,
|
|
"select node_id from gis.region_nodes where region_id = %s order by node_id",
|
|
(region_id,),
|
|
)
|
|
if stored:
|
|
return [str(row["node_id"]) for row in stored]
|
|
rows = read_all(
|
|
name,
|
|
"select n.node_id from gis.node_geometries n join gis.regions r "
|
|
"on st_intersects(n.geom, r.boundary) where r.id = %s order by n.node_id",
|
|
(region_id,),
|
|
)
|
|
return [str(row["node_id"]) for row in rows]
|
|
|
|
|
|
def get_links_on_region_boundary(name: str, region_id: str) -> list[str]:
|
|
nodes = get_nodes_in_region(name, region_id)
|
|
return _get_links_on_boundary(name, nodes)
|
|
|
|
|
|
def calculate_convex_hull(name: str, nodes: list[str]) -> list[tuple[float, float]]:
|
|
row = read(
|
|
name,
|
|
"select st_astext(st_convexhull(st_collect(geom))) as boundary "
|
|
"from gis.node_geometries where node_id = any(%s)",
|
|
(nodes,),
|
|
)
|
|
return from_postgis_polygon(str(row["boundary"]))
|
|
|
|
|
|
def _verify_platform():
|
|
_platform = platform.system()
|
|
if _platform != "Windows":
|
|
raise Exception(f'Platform {_platform} unsupported (not yet)')
|
|
|
|
|
|
def _normal(v: tuple[float, float]) -> tuple[float, float]:
|
|
l = math.sqrt(v[0] * v[0] + v[1] * v[1])
|
|
return (v[0] / l, v[1] / l)
|
|
|
|
|
|
def _angle(v: tuple[float, float]) -> float:
|
|
if v[0] >= 0 and v[1] >= 0:
|
|
return math.asin(v[1])
|
|
elif v[0] <= 0 and v[1] >= 0:
|
|
return math.pi - math.asin(v[1])
|
|
elif v[0] <= 0 and v[1] <= 0:
|
|
return math.asin(-v[1]) + math.pi
|
|
elif v[0] >= 0 and v[1] <= 0:
|
|
return math.pi * 2 - math.asin(-v[1])
|
|
return 0
|
|
|
|
|
|
def _angle_of_node_link(node: str, link: str, nodes, links) -> float:
|
|
n1 = node
|
|
n2 = links[link]['node1'] if n1 == links[link]['node2'] else links[link]['node2']
|
|
x1, y1 = nodes[n1]['x'], nodes[n1]['y']
|
|
x2, y2 = nodes[n2]['x'], nodes[n2]['y']
|
|
if y1 == y2:
|
|
v = ((x2 - x1) / abs(x2 - x1), 0.0)
|
|
else:
|
|
v = _normal((x2 - x1, y2 - y1))
|
|
return _angle(v)
|
|
|
|
|
|
class Topology:
|
|
def __init__(self, db: str, nodes: list[str]) -> None:
|
|
node_rows, link_rows = get_topology_rows(db, nodes)
|
|
|
|
self._nodes: dict[str, Any] = {
|
|
str(row["id"]): {
|
|
"x": float(row["x"]),
|
|
"y": float(row["y"]),
|
|
"type": str(row["node_type"]),
|
|
"links": [],
|
|
}
|
|
for row in node_rows
|
|
}
|
|
self._node_list = list(self._nodes)
|
|
self._max_x_node = max(
|
|
self._nodes,
|
|
key=lambda node_id: self._nodes[node_id]["x"],
|
|
default="",
|
|
)
|
|
|
|
self._links: dict[str, Any] = {}
|
|
for row in link_rows:
|
|
link_id = str(row["id"])
|
|
node1 = str(row["start_node_id"])
|
|
node2 = str(row["end_node_id"])
|
|
if node1 not in self._nodes or node2 not in self._nodes:
|
|
continue
|
|
self._links[link_id] = {
|
|
"node1": node1,
|
|
"node2": node2,
|
|
"length": float(row["length"]),
|
|
}
|
|
self._nodes[node1]["links"].append(link_id)
|
|
self._nodes[node2]["links"].append(link_id)
|
|
self._link_list = list(self._links)
|
|
|
|
def nodes(self):
|
|
return self._nodes
|
|
|
|
def node_list(self):
|
|
return self._node_list
|
|
|
|
def max_x_node(self):
|
|
return self._max_x_node
|
|
|
|
def links(self):
|
|
return self._links
|
|
|
|
def link_list(self):
|
|
return self._link_list
|
|
|
|
|
|
def _calculate_boundary(cursor: str, t_nodes: dict[str, Any], t_links: dict[str, Any]) -> tuple[list[str], dict[str, list[str]], list[tuple[float, float]]]:
|
|
in_angle = 0
|
|
|
|
vertices: list[str] = []
|
|
path: dict[str, list[str]] = {}
|
|
while True:
|
|
# prevent duplicated node
|
|
if len(vertices) > 0 and cursor == vertices[-1]:
|
|
break
|
|
|
|
# prevent duplicated path
|
|
if len(vertices) >= 3 and vertices[0] == vertices[-1] and vertices[1] == cursor:
|
|
break
|
|
|
|
vertices.append(cursor)
|
|
|
|
sorted_links = []
|
|
overlapped_link = ''
|
|
for link in t_nodes[cursor]['links']:
|
|
angle = _angle_of_node_link(cursor, link, t_nodes, t_links)
|
|
if angle == in_angle:
|
|
overlapped_link = link
|
|
continue
|
|
sorted_links.append((angle, link))
|
|
|
|
# work into a branch, return
|
|
if len(sorted_links) == 0:
|
|
path[overlapped_link] = []
|
|
cursor = vertices[-2]
|
|
in_angle = _angle_of_node_link(cursor, overlapped_link, t_nodes, t_links)
|
|
continue
|
|
|
|
sorted_links = sorted(sorted_links, key=lambda s:s[0])
|
|
out_link = sorted_links[0][1]
|
|
for angle, link in sorted_links:
|
|
if angle > in_angle:
|
|
out_link = link
|
|
break
|
|
|
|
path[out_link] = []
|
|
cursor = t_links[out_link]['node1'] if cursor == t_links[out_link]['node2'] else t_links[out_link]['node2']
|
|
in_angle = _angle_of_node_link(cursor, out_link, t_nodes, t_links)
|
|
|
|
boundary: list[tuple[float, float]] = []
|
|
for node in vertices:
|
|
boundary.append((t_nodes[node]['x'], t_nodes[node]['y']))
|
|
|
|
return (vertices, path, boundary)
|
|
|
|
|
|
def _collect_new_links(in_links: dict[str, list[str]], t_nodes: dict[str, Any], t_links: dict[str, Any], new_nodes: dict[str, Any], new_links: dict[str, Any]) -> tuple[dict[str, Any], dict[str, Any]]:
|
|
for link, pts in in_links.items():
|
|
node1 = t_links[link]['node1']
|
|
node2 = t_links[link]['node2']
|
|
x1, x2 = t_nodes[node1]['x'], t_nodes[node2]['x']
|
|
y1, y2 = t_nodes[node1]['y'], t_nodes[node2]['y']
|
|
|
|
if node1 not in new_nodes:
|
|
new_nodes[node1] = { 'x': x1, 'y': y1, 'links': [] }
|
|
if node2 not in new_nodes:
|
|
new_nodes[node2] = { 'x': x2, 'y': y2, 'links': [] }
|
|
|
|
x_delta = x2 - x1
|
|
y_delta = y2 - y1
|
|
use_x = abs(x_delta) > abs(y_delta)
|
|
|
|
if len(pts) == 0:
|
|
new_links[link] = t_links[link]
|
|
else:
|
|
sorted_nodes: list[tuple[float, str]] = []
|
|
sorted_nodes.append((0.0, node1))
|
|
sorted_nodes.append((1.0, node2))
|
|
i = 0
|
|
for pt in pts:
|
|
x, y = new_nodes[pt]['x'], new_nodes[pt]['y']
|
|
percent = ((x - x1) / x_delta) if use_x else ((y - y1) / y_delta)
|
|
sorted_nodes.append((percent, pt))
|
|
i += 1
|
|
sorted_nodes = sorted(sorted_nodes, key=lambda s:s[0])
|
|
|
|
for i in range(1, len(sorted_nodes)):
|
|
l = sorted_nodes[i - 1][1]
|
|
r = sorted_nodes[i][1]
|
|
new_link = f'LINK_[{l}]_[{r}]'
|
|
new_links[new_link] = { 'node1': l, 'node2': r }
|
|
|
|
return (new_nodes, new_links)
|
|
|
|
|
|
def calculate_boundary(name: str, nodes: list[str], accurate = False) -> list[tuple[float, float]]:
|
|
topology = Topology(name, nodes)
|
|
t_nodes = topology.nodes()
|
|
t_links = topology.links()
|
|
_, _, boundary = _calculate_boundary(topology.max_x_node(), t_nodes, t_links)
|
|
return boundary
|
|
|
|
|
|
'''
|
|
# CClipper2.dll
|
|
# int inflate_paths(double* path, size_t size, double delta, int jt, int et, double miter_limit, int precision, double arc_tolerance, double** out_path, size_t* out_size);
|
|
# int simplify_paths(double* path, size_t size, double epsilon, int is_closed_path, double** out_path, size_t* out_size);
|
|
# void free_paths(double** paths);
|
|
'''
|
|
def inflate_boundary(name: str, boundary: list[tuple[float, float]], delta: float = 0.5) -> list[tuple[float, float]]:
|
|
if boundary[0] == boundary[-1]:
|
|
del(boundary[-1])
|
|
|
|
precision = 2
|
|
scale = 10 ** precision
|
|
path = [(round(x * scale), round(y * scale)) for x, y in boundary]
|
|
|
|
offset = pyclipper.PyclipperOffset(miter_limit=2.0)
|
|
offset.AddPath(path, pyclipper.JT_SQUARE, pyclipper.ET_CLOSEDPOLYGON)
|
|
solutions = offset.Execute(round(delta * scale))
|
|
if len(solutions) == 0:
|
|
return []
|
|
|
|
result: list[tuple[float, float]] = []
|
|
for x, y in solutions[0]:
|
|
result.append((x / scale, y / scale))
|
|
result.append(result[0])
|
|
return result
|
|
|
|
|
|
def inflate_region(name: str, region_id: str, delta: float = 0.5) -> list[tuple[float, float]]:
|
|
r = try_read(name, "select id, st_astext(boundary) as boundary_geom from gis.regions where id = %s", (region_id,))
|
|
if r == None:
|
|
return []
|
|
boundary = from_postgis_polygon(str(r['boundary_geom']))
|
|
return inflate_boundary(name, boundary, delta)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
_verify_platform()
|