229 lines
8.5 KiB
Python
229 lines
8.5 KiB
Python
#!/usr/bin/env python3
|
|
"""Normalize malformed lane-arrow meshes emitted by osm2streets."""
|
|
|
|
import argparse
|
|
import json
|
|
import math
|
|
from pathlib import Path
|
|
|
|
from osgeo import ogr
|
|
|
|
|
|
METERS_PER_DEGREE = 111320.0
|
|
ogr.UseExceptions()
|
|
|
|
|
|
def cli_args():
|
|
parser = argparse.ArgumentParser(
|
|
description="Normalize osm2streets triangulated lane-arrow polygons."
|
|
)
|
|
parser.add_argument("--input", required=True, type=Path)
|
|
parser.add_argument("--outline-simplify-meters", required=True, type=float)
|
|
return parser.parse_args()
|
|
|
|
|
|
def point_segment_distance(point, start, end, latitude):
|
|
meters_lon = METERS_PER_DEGREE * math.cos(math.radians(latitude))
|
|
px = (point[0] - start[0]) * meters_lon
|
|
py = (point[1] - start[1]) * METERS_PER_DEGREE
|
|
bx = (end[0] - start[0]) * meters_lon
|
|
by = (end[1] - start[1]) * METERS_PER_DEGREE
|
|
length_squared = bx * bx + by * by
|
|
if length_squared == 0:
|
|
return math.hypot(px, py)
|
|
projection = max(0.0, min(1.0, (px * bx + py * by) / length_squared))
|
|
return math.hypot(px - projection * bx, py - projection * by)
|
|
|
|
|
|
def simplify_ring(points, tolerance, latitude):
|
|
points = list(points)
|
|
changed = True
|
|
while changed and len(points) > 3:
|
|
changed = False
|
|
for index, point in enumerate(points):
|
|
previous = points[index - 1]
|
|
following = points[(index + 1) % len(points)]
|
|
if point_segment_distance(point, previous, following, latitude) <= tolerance:
|
|
points.pop(index)
|
|
changed = True
|
|
break
|
|
return points
|
|
|
|
|
|
def tail_edge_candidate(points):
|
|
best = None
|
|
for index in range(len(points)):
|
|
start = points[index]
|
|
end = points[(index + 1) % len(points)]
|
|
previous = points[index - 1]
|
|
following = points[(index + 2) % len(points)]
|
|
before = [start[0] - previous[0], start[1] - previous[1]]
|
|
after = [following[0] - end[0], following[1] - end[1]]
|
|
edge = [end[0] - start[0], end[1] - start[1]]
|
|
before_length = math.hypot(*before)
|
|
after_length = math.hypot(*after)
|
|
edge_length = math.hypot(*edge)
|
|
if min(before_length, after_length, edge_length) == 0:
|
|
continue
|
|
alignment = (
|
|
before[0] * after[0] + before[1] * after[1]
|
|
) / (before_length * after_length)
|
|
if before_length <= edge_length or after_length <= edge_length or alignment >= -0.9:
|
|
continue
|
|
score = -alignment * min(before_length, after_length) / edge_length
|
|
if best is None or score > best[0]:
|
|
best = (score, index, before, after, before_length, after_length)
|
|
return best
|
|
|
|
|
|
def square_arrow_tail(points, latitude):
|
|
# A normalized straight arrow has seven exterior vertices. Other arrow
|
|
# silhouettes are left untouched because their tail cannot be inferred safely.
|
|
if len(points) != 7:
|
|
return points
|
|
meters_lon = METERS_PER_DEGREE * math.cos(math.radians(latitude))
|
|
origin = points[0]
|
|
local = [
|
|
[
|
|
(point[0] - origin[0]) * meters_lon,
|
|
(point[1] - origin[1]) * METERS_PER_DEGREE,
|
|
]
|
|
for point in points
|
|
]
|
|
candidate = tail_edge_candidate(local)
|
|
if candidate is None:
|
|
return points
|
|
_, index, before, after, before_length, after_length = candidate
|
|
axis = [
|
|
before[0] / before_length - after[0] / after_length,
|
|
before[1] / before_length - after[1] / after_length,
|
|
]
|
|
axis_length = math.hypot(*axis)
|
|
if axis_length == 0:
|
|
return points
|
|
axis = [axis[0] / axis_length, axis[1] / axis_length]
|
|
end_index = (index + 1) % len(local)
|
|
midpoint = [
|
|
(local[index][0] + local[end_index][0]) / 2.0,
|
|
(local[index][1] + local[end_index][1]) / 2.0,
|
|
]
|
|
for point_index in (index, end_index):
|
|
offset = [
|
|
local[point_index][0] - midpoint[0],
|
|
local[point_index][1] - midpoint[1],
|
|
]
|
|
projection = offset[0] * axis[0] + offset[1] * axis[1]
|
|
local[point_index][0] -= projection * axis[0]
|
|
local[point_index][1] -= projection * axis[1]
|
|
points[point_index] = [
|
|
origin[0] + local[point_index][0] / meters_lon,
|
|
origin[1] + local[point_index][1] / METERS_PER_DEGREE,
|
|
]
|
|
return points
|
|
|
|
|
|
def normalize_polygon(geometry, tolerance):
|
|
if geometry.GetGeometryName() == "MULTIPOLYGON":
|
|
geometry = geometry.UnionCascaded()
|
|
if geometry is None or geometry.GetGeometryName() != "POLYGON":
|
|
raise ValueError("triangle merge did not produce a Polygon")
|
|
|
|
source_ring = geometry.GetGeometryRef(0)
|
|
points = [source_ring.GetPoint(i)[:2] for i in range(source_ring.GetPointCount() - 1)]
|
|
if len(points) <= 3:
|
|
raise ValueError("arrow exterior has too few points")
|
|
latitude = sum(point[1] for point in points) / len(points)
|
|
points = simplify_ring(points, tolerance, latitude)
|
|
points = square_arrow_tail(points, latitude)
|
|
|
|
normalized = ogr.Geometry(ogr.wkbPolygon)
|
|
outer = ogr.Geometry(ogr.wkbLinearRing)
|
|
for point in points + [points[0]]:
|
|
outer.AddPoint_2D(*point)
|
|
normalized.AddGeometry(outer)
|
|
for index in range(1, geometry.GetGeometryCount()):
|
|
normalized.AddGeometry(geometry.GetGeometryRef(index))
|
|
if normalized.IsEmpty() or not normalized.IsValid():
|
|
raise ValueError("normalized arrow geometry is invalid")
|
|
return normalized
|
|
|
|
|
|
def normalize_file(input_path, tolerance):
|
|
with input_path.open("r", encoding="utf-8") as handle:
|
|
collection = json.load(handle)
|
|
|
|
collection["features"] = merge_custom_turn_lane_arrow_components(
|
|
collection.get("features", [])
|
|
)
|
|
for index, feature in enumerate(collection["features"]):
|
|
geometry = ogr.CreateGeometryFromJson(json.dumps(feature.get("geometry", {})))
|
|
if geometry is None or geometry.IsEmpty():
|
|
raise ValueError(f"feature {index} has no usable geometry")
|
|
try:
|
|
normalized = normalize_polygon(geometry, tolerance)
|
|
except ValueError as error:
|
|
raise ValueError(f"feature {index}: {error}") from error
|
|
feature["geometry"] = json.loads(
|
|
normalized.ExportToJson(options=["COORDINATE_PRECISION=15"])
|
|
)
|
|
|
|
temp_path = input_path.with_name(f".{input_path.name}.tmp")
|
|
with temp_path.open("w", encoding="utf-8") as handle:
|
|
json.dump(collection, handle, ensure_ascii=False, separators=(",", ":"))
|
|
handle.write("\n")
|
|
temp_path.replace(input_path)
|
|
|
|
|
|
def merge_custom_turn_lane_arrow_components(features):
|
|
"""Restore each imported SVG icon before QGIS assigns its single outline."""
|
|
grouped = {}
|
|
for feature in features:
|
|
properties = feature.get("properties", {})
|
|
arrow_id = properties.get("custom_arrow_id")
|
|
if properties.get("source") == "osm_turn_lanes" and arrow_id:
|
|
grouped.setdefault(arrow_id, []).append(feature)
|
|
|
|
merged_ids = set()
|
|
output = []
|
|
for feature in features:
|
|
properties = feature.get("properties", {})
|
|
arrow_id = properties.get("custom_arrow_id")
|
|
if properties.get("source") != "osm_turn_lanes" or not arrow_id:
|
|
output.append(feature)
|
|
continue
|
|
if arrow_id in merged_ids:
|
|
continue
|
|
merged_ids.add(arrow_id)
|
|
components = grouped[arrow_id]
|
|
multi = ogr.Geometry(ogr.wkbMultiPolygon)
|
|
for component in components:
|
|
geometry = ogr.CreateGeometryFromJson(
|
|
json.dumps(component.get("geometry", {}))
|
|
)
|
|
if geometry is None or geometry.IsEmpty() or geometry.GetGeometryName() != "POLYGON":
|
|
raise ValueError(f"custom arrow {arrow_id} has an invalid component")
|
|
multi.AddGeometry(geometry)
|
|
united = multi.UnionCascaded()
|
|
if united is None or united.IsEmpty() or united.GetGeometryName() != "POLYGON":
|
|
raise ValueError(
|
|
f"custom arrow {arrow_id} components did not form one Polygon"
|
|
)
|
|
merged = dict(components[0])
|
|
merged["properties"] = dict(components[0].get("properties", {}))
|
|
merged["properties"].pop("arrow_part", None)
|
|
merged["properties"]["component_count"] = len(components)
|
|
merged["geometry"] = json.loads(united.ExportToJson())
|
|
output.append(merged)
|
|
return output
|
|
|
|
|
|
def main():
|
|
args = cli_args()
|
|
if args.outline_simplify_meters < 0:
|
|
raise ValueError("--outline-simplify-meters must be non-negative")
|
|
normalize_file(args.input.resolve(), args.outline_simplify_meters)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|