osmWorkflow/blender/osmassets/osm.py

"""OSM XML parsing and the local metric projection.

Pure Python: no `bpy`, so this runs and tests outside Blender.
"""

import math
import xml.etree.ElementTree as ET


def _element_id(element):
    return element.attrib.get("id", "")


def tags(element):
    return {t.attrib.get("k", ""): t.attrib.get("v", "")
            for t in element.findall("tag")}


def _coords_from_refs(refs, nodes):
    return [nodes[r] for r in refs if r in nodes]


def _open_ring(coords):
    if len(coords) > 1 and coords[0] == coords[-1]:
        return coords[:-1]
    return coords


def _is_closed_ring(coords):
    return len(coords) >= 4 and coords[0] == coords[-1] and len(_open_ring(coords)) >= 3


def _join_member_rings(members):
    """Build closed rings from relation member coordinate runs.

    Handles the common OSM multipolygon cases: members are either already
    closed ways, or open way fragments whose endpoints can be stitched together.
    Malformed leftovers are dropped rather than aborting the whole import.
    """
    rings = []
    pending = [list(member) for member in members if len(member) >= 2]
    while pending:
        ring = pending.pop(0)
        changed = True
        while not _is_closed_ring(ring) and changed:
            changed = False
            for index, candidate in enumerate(pending):
                if ring[-1] == candidate[0]:
                    ring.extend(candidate[1:])
                elif ring[-1] == candidate[-1]:
                    ring.extend(reversed(candidate[:-1]))
                elif ring[0] == candidate[-1]:
                    ring = candidate[:-1] + ring
                elif ring[0] == candidate[0]:
                    ring = list(reversed(candidate[1:])) + ring
                else:
                    continue
                pending.pop(index)
                changed = True
                break
        if _is_closed_ring(ring):
            rings.append(ring)
    return rings


def parse_osm(path):
    root = ET.parse(path).getroot()
    bounds_node = root.find("bounds")
    if bounds_node is None:
        raise RuntimeError("OSM file does not contain a bounds element")
    bounds = {"min_lon": float(bounds_node.attrib["minlon"]),
              "min_lat": float(bounds_node.attrib["minlat"]),
              "max_lon": float(bounds_node.attrib["maxlon"]),
              "max_lat": float(bounds_node.attrib["maxlat"])}

    nodes = {}
    point_features = []
    for node in root.findall("node"):
        try:
            node_id = int(node.attrib["id"])
            coord = (float(node.attrib["lon"]), float(node.attrib["lat"]))
            node_tags = tags(node)
            nodes[node_id] = coord
            if node_tags:
                point_features.append({"id": node.attrib.get("id", ""),
                                       "coord": coord, "tags": node_tags})
        except (KeyError, ValueError):
            continue

    ways = []
    way_coords_by_id = {}
    for way in root.findall("way"):
        if way.attrib.get("action") == "delete":
            continue
        refs = []
        for ref in way.findall("nd"):
            try:
                refs.append(int(ref.attrib["ref"]))
            except (KeyError, ValueError):
                pass
        coords = _coords_from_refs(refs, nodes)
        way_id = _element_id(way)
        if len(coords) >= 2:
            way_coords_by_id[way_id] = coords
        way_tags = tags(way)
        if len(coords) >= 2:
            ways.append({"id": way_id, "coords": coords, "tags": way_tags})
    for relation in root.findall("relation"):
        if relation.attrib.get("action") == "delete":
            continue
        relation_tags = tags(relation)
        if relation_tags.get("type") != "multipolygon" or "building" not in relation_tags:
            continue
        outer_members = []
        inner_members = []
        for member in relation.findall("member"):
            if member.attrib.get("type") != "way":
                continue
            coords = way_coords_by_id.get(member.attrib.get("ref", ""))
            if not coords:
                continue
            role = member.attrib.get("role", "")
            if role == "inner":
                inner_members.append(coords)
            elif role in ("", "outer"):
                outer_members.append(coords)
        outer_rings = _join_member_rings(outer_members)
        if not outer_rings:
            continue
        inner_rings = _join_member_rings(inner_members)
        relation_id = _element_id(relation)
        for index, outer in enumerate(outer_rings):
            synthetic_id = relation_id if len(outer_rings) == 1 else f"{relation_id}:{index + 1}"
            ways.append({
                "id": synthetic_id,
                "coords": outer,
                "inner_coords": inner_rings,
                "tags": relation_tags,
                "source": "relation",
            })
    return bounds, ways, point_features


def parse_height(feature_tags, default):
    try:
        return max(0.5, float(feature_tags.get("height", default)))
    except (TypeError, ValueError):
        return default


class Projector:
    """Equirectangular projection about the centre of the OSM bounds.

    Output is metres in a local ENU frame (X east, Y north), which is what both
    the Blender scene and the Cesium GLB are authored in.
    """

    def __init__(self, bounds):
        self.bounds = bounds
        self.lon0 = (bounds["min_lon"] + bounds["max_lon"]) / 2
        self.lat0 = (bounds["min_lat"] + bounds["max_lat"]) / 2
        self.m_per_lat = 111320.0
        self.m_per_lon = 111320.0 * math.cos(math.radians(self.lat0))

    def xy(self, lon_lat):
        lon, lat = lon_lat
        return ((lon - self.lon0) * self.m_per_lon,
                (lat - self.lat0) * self.m_per_lat)

    def inside(self, lon_lat, pad=0.00035):
        lon, lat = lon_lat
        b = self.bounds
        return (b["min_lon"] - pad <= lon <= b["max_lon"] + pad and
                b["min_lat"] - pad <= lat <= b["max_lat"] + pad)

    def ring(self, coords):
        return [self.xy(c) for c in coords]