解决FlightRadar高度转WGS84椭球后在Cesium中低于地面的问题

我从FlightRadar导出了包含飞行高度的航班CSV文件，想在Cesium中使用。但文件里地面高度始终为0，导致Cesium中显示的飞行高度低于地面，所以需要考虑WGS84椭球和地形高程两个因素。我用下面的代码把英尺单位的平均海平面（MSL）高度转换成米单位的椭球高度：

from pyproj import Transformer
from pygeodesy import ellipsoidalVincenty as ev

def feet_msl_to_meters_ellipsoid(lat, lon, height_ft_msl):

    height_m_msl = height_ft_msl * 0.3048
    point = ev.LatLon(lat, lon)
    geoid_height = point.height  # Default EGM96
    height_ellipsoid = height_m_msl + geoid_height

    return height_ellipsoid

但仍然有部分点显示在地面下方，我觉得这和地形高程有关。请问怎么通过坐标考虑地形高程因素，或者还有其他问题吗？

更新：附上问题可视化图片：

1. 核心差异：椭球高度 vs 地形上方高度 #

你当前的转换仅完成了平均海平面（MSL）到WGS84椭球高度的转换，但Cesium显示的地形是实际地表的椭球高程。航班的MSL高度是相对于海平面的高度，要让航班始终显示在地形上方，需要明确：
最终椭球高度 = 地形椭球高程 + 航班相对于地面的安全高度
而FlightRadar导出的高度多为MSL高度，因此需要先将其转换为「相对于地面的高度」，再叠加地形的椭球高程。

2. 获取地形椭球高程的两种方法 #

方法一：离线预计算（适合批量处理CSV）

使用SRTM地形数据集（如30米分辨率GeoTIFF），提前批量查询每个经纬度对应的地形高程（注意SRTM数据是MSL高度，需转换为椭球高度）：

from osgeo import gdal
from pyproj import Transformer

# 加载本地SRTM地形文件
ds = gdal.Open("srtm_terrain.tif")
band = ds.GetRasterBand(1)
gt = ds.GetGeoTransform()

# 定义EGM96转WGS84椭球的转换器
transformer = Transformer.from_crs("EPSG:4326+5773", "EPSG:4326+4978", always_xy=True)

def get_terrain_ellipsoid(lat, lon):
    # 经纬度转图像像素坐标
    x = int((lon - gt[0]) / gt[1])
    y = int((lat - gt[3]) / gt[5])
    # 获取地形MSL高度
    terrain_msl = band.ReadAsArray(x, y, 1, 1)[0][0]
    # 转换为椭球高度
    _, _, terrain_ellipsoid = transformer.transform(lon, lat, terrain_msl)
    return terrain_ellipsoid

之后对每个航班点：

# 假设csv中每行数据为(lat, lon, height_ft_msl)
terrain_ellip = get_terrain_ellipsoid(lat, lon)
flight_ellip = feet_msl_to_meters_ellipsoid(lat, lon, height_ft_msl)
# 确保航班高度高于地形（加10米安全冗余）
final_height = max(flight_ellip, terrain_ellip + 10)

方法二：Cesium在线实时修正

如果不想离线处理，可在Cesium中通过API实时获取地形高程并调整：

const viewer = new Cesium.Viewer('cesiumContainer');
const terrainProvider = viewer.terrainProvider;

// 假设flightPoints是包含经纬度、椭球高度的数组
Cesium.sampleTerrainMostDetailed(terrainProvider, 
  flightPoints.map(p => Cesium.Cartographic.fromDegrees(p.lon, p.lat))
).then(sampledPoints => {
  sampledPoints.forEach((sampled, index) => {
    const flightHeight = flightPoints[index].ellipsoidHeight;
    // 若航班高度低于地形，调整到地形上方
    if (flightHeight < sampled.height) {
      flightPoints[index].ellipsoidHeight = sampled.height + 10;
    }
  });
  // 更新Cesium中的航班可视化
  updateFlightLayer(flightPoints);
});

3. 检查FlightRadar高度数据定义 #

• 滑行阶段的航班高度可能标记为0（MSL），但当地地形高程可能高于海平面，这类点需单独处理（直接设为地形高程+5米）。
• 确认导出的高度是气压高度（MSL）还是无线电高度（相对地面）：如果是无线电高度，直接叠加地形椭球高程即可。

4. 升级大地水准面模型 #

你当前使用的EGM96模型部分区域误差较大，可尝试用更精准的EGM2008模型转换：

from pyproj import Transformer

def feet_msl_to_meters_ellipsoid_egm2008(lat, lon, height_ft_msl):
    height_m_msl = height_ft_msl * 0.3048
    # EGM2008转WGS84椭球
    transformer = Transformer.from_crs("EPSG:4326+3855", "EPSG:4326+4978", always_xy=True)
    _, _, height_ellipsoid = transformer.transform(lon, lat, height_m_msl)
    return height_ellipsoid

内容的提问来源于stack exchange，提问作者Fish1996