210 lines
7.5 KiB
Python
210 lines
7.5 KiB
Python
"""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]
|