[convert_to_local算法]与[经纬高转东北天算法]精度对比

1 介绍

convert_to_local算法如下，

//C++版
void CHAD2BDUtil::convert_to_local(
    const double start_long, const double start_lat,
    const double in_long, const double in_lat,
    double &out_x, double &out_y)
{
    double start_lon_rad = start_long * PI / 180.0;
    double start_lat_rad = start_lat * PI / 180.0;
    double lon_rad = in_long * PI / 180.0;
    double lat_rad = in_lat * PI / 180.0;
    double r = 6371.004 * 1.0e3;

    out_x = r * cos(lat_rad) * sin(lon_rad - start_lon_rad);
    out_y = r * sin(lat_rad - start_lat_rad);
}
//python版
def convert_to_local(start_long: float, start_lat: float, in_long: float, in_lat: float) -> list():
    out_x = 0.0
    out_y = 0.0
    start_lon_rad = start_long * math.pi / 180.0
    start_lat_rad = start_lat * math.pi / 180.0
    lon_rad = in_long * math.pi / 180.0
    lat_rad = in_lat * math.pi / 180.0
    r = 6371.004 * 1000.0

    out_x = r * math.cos(lat_rad) * math.sin(lon_rad - start_lon_rad)
    out_y = r * math.sin(lat_rad - start_lat_rad)


    return [out_x, out_y]

经纬高转东北天算法如下，

import math

#经度纬度高度转ECEF的X,Y,Z
#注意输入的经度和纬度的单位是度，高度单位是米；输出的X,Y,Z的单位是米
def lla2xyz(longitude: float, latitude: float, h: float) -> list[float]:
    longitude = longitude / 180.0 * math.pi
    latitude = latitude / 180.0 * math.pi
    
    a = 6378137.0
    b = 6356752.3142
    e2 = (a**2 - b**2) / a**2
    N = a / math.sqrt(1 - e2 * math.sin(latitude) * math.sin(latitude))
    x = (N + h) * math.cos(latitude) * math.cos(longitude)
    y = (N + h) * math.cos(latitude) * math.sin(longitude)
    z = (N * (1 - e2) + h) * math.sin(latitude)
    return [x, y, z]

#输入2个ECEF坐标，输出一个向量(以第1个点为原点)
#longitude1和latitude1分别表示第1个点的经度和纬度，x1,y1和z1表示第1个点的ECEF坐标
#x2,y2,z2表示第2个点的ECEF坐标
#返回一个列表，里面有三个元素，表示以第1个点为原点，东北天为xyz轴方向，在此坐标系下，第2个点的坐标
def xyz2enu(longitude1: float, latitude1: float, x1: float, y1: float, z1: float, x2: float, y2: float, z2: float) -> list[float]:
    longitude1 = longitude1 / 180.0 * math.pi
    latitude1 = latitude1 / 180.0 * math.pi

    de = -1 * math.sin(longitude1) * (x2 - x1) + math.cos(longitude1) * (y2 - y1)
    dn = -1 * math.sin(latitude1) * math.cos(longitude1) * (x2 - x1) - math.sin(latitude1) * math.sin(longitude1) * (y2 - y1) + math.cos(latitude1) * (z2 - z1)
    du = math.cos(latitude1) * math.cos(longitude1) * (x2 - x1) + math.cos(latitude1) * math.sin(longitude1) * (y2 - y1) + math.sin(latitude1) * (z2 - z1)
    return [de, dn, du]

def lla2enu(start_longi: float, start_lati: float, start_h: float, longi: float, lati: float, h: float) -> list[float]:
    x1, y1, z1 = lla2xyz(start_longi, start_lati, start_h)
    x2, y2, z2 = lla2xyz(longi, lati, h)
    de, dn, du = xyz2enu(start_longi, start_lati, x1, y1, z1, x2, y2, z2)
    return [de, dn, du]

2 对比

输入同一组数据，共383个几何点，QGIS显示的形状为，

使用[convert_to_local算法]的结果为，

使用[经纬高转东北天算法]的结果为，

3 结论

[convert_to_local算法]与[经纬高转东北天算法]在轨迹显示的总体效果上相差不大。