#!/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()