feat: osm2streets上游箭头修复

scripts/normalize-lane-arrows.py

@@ -0,0 +1,182 @@
+#!/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)
+
+    for index, feature in enumerate(collection.get("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 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()