将CH1903+/LV95（EPSG:2056）DEM转换为WGS84并做垂直校正

我正在尝试将瑞士Alti3D DEM的高程数据与WGS84（EPSG:4326）格式的数据对比，最初用GDAL的gdal.Warp函数将DEM重投影到WGS84并匹配参考栅格的分辨率，代码如下：

from osgeo import gdal

def match_resolution(reference_raster: str, target_raster: str, resampling_alg: str = "cubic"):
    """
    将目标栅格重投影以匹配参考栅格的分辨率，使用GDAL Warp函数。
    
    参数:
    - reference_raster: 参考栅格文件路径，目标栅格将匹配其分辨率。
    - target_raster: 需要重投影的目标栅格文件路径。
    - resampling_alg: 重采样算法，默认'cubic'，可选'nearest'、'bilinear'等。
    """
    
    # 打开参考栅格获取分辨率信息
    ref_ds = gdal.Open(reference_raster)
    if ref_ds is None:
        raise ValueError(f"无法打开参考栅格: {reference_raster}")
    
    ref_geotransform = ref_ds.GetGeoTransform()
    ref_proj = ref_ds.GetProjection()
    ref_x_res = ref_geotransform[1]  # X方向像素分辨率
    ref_y_res = abs(ref_geotransform[5])  # Y方向像素分辨率
    
    ref_ds = None
    
    # 设置Warp参数
    warp_options = gdal.WarpOptions(
        format="GTiff",
        xRes=ref_x_res,
        yRes=ref_y_res,
        resampleAlg=resampling_alg,
        dstSRS=ref_proj
    )
    output_raster = target_raster.replace(".tif", "_rm-12m-4326.tif")
    gdal.Warp(destNameOrDestDS=output_raster, srcDSOrSrcDSTab=target_raster, options=warp_options)
    
    print(f"重投影后的栅格已保存至: {output_raster}")

运行后结果出现约50米的偏差，发现问题在于瑞士DEM数据需要垂直校正，但找不到对应的校正网格，也不清楚如何用Python结合GDAL/PROJ解决。

GDAL识别的瑞士Alti3D DEM完整投影信息如下：

'PROJCS["CH1903+ / LV95",GEOGCS["CH1903+",DATUM["CH1903+",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6150"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4150"]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.9524055555556],PARAMETER["longitude_of_center",7.43958333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2056"]]'

解决方案 #

1. 获取官方垂直校正网格

瑞士联邦地形局提供了用于LV95/LV03与WGS84之间垂直转换的CHENyx06校正网格（.gtx格式），获取后需放置到PROJ的网格搜索路径中（默认路径为proj/share/proj/，或通过PROJ_LIB环境变量指定自定义路径）。

2. 修改重投影代码，启用垂直基准转换

默认的gdal.Warp仅处理水平坐标转换，需明确指定垂直基准才能触发校正。瑞士Alti3D DEM的垂直基准是LN02（EPSG:5728），目标WGS84的垂直基准是WGS84大地高（EPSG:4979），修改后的代码如下：

from osgeo import gdal, osr

def warp_with_vertical_correction(reference_raster: str, target_raster: str, resampling_alg: str = "cubic"):
    # 读取参考栅格参数
    ref_ds = gdal.Open(reference_raster)
    if ref_ds is None:
        raise ValueError(f"无法打开参考栅格: {reference_raster}")
    
    ref_geotransform = ref_ds.GetGeoTransform()
    ref_x_res = ref_geotransform[1]
    ref_y_res = abs(ref_geotransform[5])
    ref_proj = ref_ds.GetProjection()
    ref_ds = None

    # 定义源坐标系（LV95水平基准 + LN02垂直基准）
    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(2056)  # LV95水平坐标
    src_srs.SetVertCS("LN02", osr.SRS_VC_GEOIDAL, osr.SRS_UOM_METRE, "EPSG:5728")

    # 定义目标坐标系（WGS84水平基准 + WGS84大地高）
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromWkt(ref_proj)  # 继承参考栅格的WGS84水平投影
    dst_srs.SetVertCS("WGS84", osr.SRS_VC_ELLIPSOIDAL, osr.SRS_UOM_METRE, "EPSG:4979")

    # 配置Warp选项，启用垂直转换
    warp_options = gdal.WarpOptions(
        format="GTiff",
        xRes=ref_x_res,
        yRes=ref_y_res,
        resampleAlg=resampling_alg,
        srcSRS=src_srs.ExportToWkt(),
        dstSRS=dst_srs.ExportToWkt(),
        dstNodata=-9999,
        warpOptions=["VECTORIZE_ALL"]  # 确保垂直转换应用到所有像素
    )

    output_raster = target_raster.replace(".tif", "_rm-12m-4326-corrected.tif")
    gdal.Warp(destNameOrDestDS=output_raster, srcDSOrSrcDSTab=target_raster, options=warp_options)
    print(f"校正并重投影后的栅格已保存至: {output_raster}")

3. 验证转换有效性

可以通过PROJ命令验证垂直转换是否配置成功：

projinfo -s EPSG:2056+5728 -t EPSG:4326+4979

若输出中包含CHENyx06相关的转换步骤，说明校正网格已正确加载。

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