"""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: """WGS84 ECEF to local ENU projection about the OSM bounds centre. 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. Cesium places the GLB with eastNorthUpToFixedFrame, so using the same ellipsoid transform is required to keep route coordinates aligned across the whole scene. """ WGS84_A = 6378137.0 WGS84_E2 = 6.6943799901413165e-3 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._lon0_rad = math.radians(self.lon0) self._lat0_rad = math.radians(self.lat0) self._sin_lon0 = math.sin(self._lon0_rad) self._cos_lon0 = math.cos(self._lon0_rad) self._sin_lat0 = math.sin(self._lat0_rad) self._cos_lat0 = math.cos(self._lat0_rad) self._origin_ecef = self._ecef(self._lon0_rad, self._lat0_rad) denominator = math.sqrt(1.0 - self.WGS84_E2 * self._sin_lat0 ** 2) prime_vertical_radius = self.WGS84_A / denominator meridional_radius = self.WGS84_A * (1.0 - self.WGS84_E2) / denominator ** 3 radians_per_degree = math.pi / 180.0 self.m_per_lon = prime_vertical_radius * self._cos_lat0 * radians_per_degree self.m_per_lat = meridional_radius * radians_per_degree def xy(self, lon_lat): lon, lat = lon_lat x, y, z = self._ecef(math.radians(lon), math.radians(lat)) dx = x - self._origin_ecef[0] dy = y - self._origin_ecef[1] dz = z - self._origin_ecef[2] east = -self._sin_lon0 * dx + self._cos_lon0 * dy north = (-self._sin_lat0 * self._cos_lon0 * dx - self._sin_lat0 * self._sin_lon0 * dy + self._cos_lat0 * dz) return east, north def _ecef(self, lon_rad, lat_rad): sin_lat = math.sin(lat_rad) cos_lat = math.cos(lat_rad) radius = self.WGS84_A / math.sqrt(1.0 - self.WGS84_E2 * sin_lat ** 2) return (radius * cos_lat * math.cos(lon_rad), radius * cos_lat * math.sin(lon_rad), radius * (1.0 - self.WGS84_E2) * sin_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]