python代码实现xyz点文件高程基准转换（从wgs84到85国家高程）

前言

从nasa上下载的12.5m的dem在和国内的dem融合过程中发现他们之间存在高程差值，nasa上12.5m的dem的高程基准是WGS84，而国内的数据大多采用的是1985国家高程基准；即前一个测量的是椭球高，后一个测量的是大地高。一般来讲在我国东北地区，椭球高>大地高。

目前由于我的数据需要，只用处理xyz点文件即可，即使用python代码批量将所有点的z值从基于WGS84的高程基准转为基于85国家高程基准。

一、基本原理

大地高 (Orthometric height)= 椭球高 (Ellipsodial height) – 大地水准面差距(Geoid height)

大地高：1985国家基准

椭球高：WGS84

大地水准面差距：可从全球大地水准面模型中获得

全球大地水准面模型主要有：EGM2008、EGM84、EGM96，我在编写python代码的时候采用了EGM2008模型，使用的模型精度是2.5'。EGM2008的模型精度有1'、2.5'、5'，数值越大模型精度越小。我根据我的工作需要选择了2.5'精度的。

EGM2008模型下载地址：https://www.agisoft.com/zh-cn/downloads/geoids/

二、python代码

import rasterio
import pandas as pd
import numpy as np

# 1. 配置文件路径
xyz_file_path = "path/to/your_file_xyz4326.xyz"  # 输入文件路径（EPSG:4326）
egm2008_tif_path = "path/to/egm2008.tif"  # EGM2008 文件路径
output_4326_path = "path/to/output_file_xyz4326.xyz"  # 输出文件路径（EPSG:4326）

# 2. 加载 EGM2008 数据
with rasterio.open(egm2008_tif_path) as geoid:
    geoid_data = geoid.read(1)  # 将 EGM2008 数据载入内存
    transform = geoid.transform  # 获取栅格的变换参数

# 获取大地水准面偏差 N 值的函数
def get_geoid_undulation(lon, lat):
    """获取大地水准面偏差 N 值"""
    col, row = ~transform * (lon, lat)  # 经纬度转栅格索引
    col, row = int(col), int(row)      # 索引取整
    if 0 <= col < geoid_data.shape[1] and 0 <= row < geoid_data.shape[0]:
        return geoid_data[row, col]
    else:
        raise ValueError(f"坐标超出范围：lon={lon}, lat={lat}")

# 3. 主处理逻辑
print("正在读取 XYZ 文件...")
xyz_data = pd.read_csv(xyz_file_path, sep='\s+', header=None, names=["Lon", "Lat", "Z"])

# 处理高程基准转换
print("正在进行高程基准转换...")
results_4326 = []

for index, row in xyz_data.iterrows():
    try:
        lon, lat, z = row["Lon"], row["Lat"], row["Z"]

        # 获取 N 值
        N = get_geoid_undulation(lon, lat)

        # 转换为 1985 国家高程基准
        H = z - N

        # 保存 EPSG:4326 的结果
        results_4326.append([lon, lat, H])

    except Exception as e:
        print(f"处理点 ({lon}, {lat}, {z}) 时发生错误: {e}")
        continue

# 保存结果到文件
print("正在保存结果到文件...")
# 保存 EPSG:4326
pd.DataFrame(results_4326, columns=["Lon", "Lat", "Z"]).to_csv(output_4326_path, sep=" ", header=False, index=False)

print(f"处理完成！结果已保存到：\nEPSG:4326 文件：{output_4326_path}")

运行该代码只需要修改 #1.配置文件路径 。

 输入的xyz文件的投影坐标系需要是EPSG:4326,不是的话需要自己转一下。输出的xyz文件也是epsg4326坐标系下的。