fix: 修正 Cesium 巡航车道中心对齐

blender/osmassets/osm.py

@@ -149,23 +149,55 @@ def parse_height(feature_tags, default):
 
 
 class Projector:
-    """Equirectangular projection about the centre of the OSM bounds.
+    """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.
+    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.m_per_lat = 111320.0
-        self.m_per_lon = 111320.0 * math.cos(math.radians(self.lat0))
+        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
-        return ((lon - self.lon0) * self.m_per_lon,
-                (lat - self.lat0) * self.m_per_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