将[经度, 纬度, 高度]格式的数据转换为[x, y, z]格式

Question

我有一个关于航班的数据集。飞机和地面站之间通过通信，利用GPS我可以得到飞机当前位置和地面站位置之间的距离，还有它们的经度、纬度和高度。

现在我需要把这些经度、纬度和高度的数据转换成一个x,y,z的坐标系统。

因为地面站的位置是固定的，我想把它放在坐标系统的原点，也就是(0,0,0)。我在寻找一个可以处理这种转换的Python模块。

作为例子，我提供了地面站的实际数据和飞机的一个位置，这些数据需要被转换：

地面站 [应该是(0, 0, 0)]
纬度: 41.947.694
经度: 3.209.083
高度(米): 379.41

飞机位置:
纬度: 419.763.015.750
经度: 34.890.851.772
高度(米): 971.32

如果我理解正确的话，把地面站放到原点，我需要从飞机的位置中减去地面站的位置，这样就变成：

地面站:
纬度: 0
经度: 0
高度(米): 0

飞机位置:
纬度: 419.721.068.056
经度: 34.887.642.689
高度(米): 591.91

因为我对处理这些数据还很陌生，尝试了一些简单的解决方案但没有找到能转换这些数据的方法，所以希望你们能帮我找到解决方案。

我需要转换后的数据来模拟3D中的椭球体，以找到可能的交点。其实也可以用其他方法来实现，比如用Matlab也可以，虽然我之前没有用过。

Answer 1

将经纬度和高度的GPS坐标转换成“xyz”坐标，这些坐标是基于WGS84参考椭球体的地理坐标（这里的垂直方向是指向椭球体的法线，而不是指向地球中心的方向），这个过程在很多编程库中都有实现。既然提到了Matlab，它的实现方式是geodetic2ecef。

这个过程在很多书籍和网站上都有描述，甚至在维基百科上也能找到，所以我们可以放心地从头开始实现它。

如果你把站点和飞机的坐标都转换成这个ECEF的x,y,z坐标，你就可以计算它们之间的差值，这样你就得到了一个以站点为(0, 0, 0)的XYZ坐标系统。

需要注意的是，正如评论中提到的，你不能直接对纬度或经度进行有意义的差值计算。要找出两点之间的距离，我们只能像处理向量或笛卡尔坐标那样进行大小的相减；纬度和经度并不是线性相关的。

现在，之前计算出的xyz差值本身并没有太大意义，因为它们的参考方向是全局的——z轴是与地球自转轴平行的，而xz平面是与本初子午线平行的。

为此，我建议对这些差值进行“本地化”，通过添加一次（双重）旋转，将这些全局方向调整到与参考站相关的方向：新的z轴指向地理垂直方向（向天顶），新的x轴指向站点的北方，y轴则指向东方。

同样，这些计算在文献和网上都有详细记录，所以在这里从头实现它们是安全的。

下面是一个例子，可能的Python实现，虽然效率不高，因为它重复计算了三角函数，但我的目的是尽量让它清晰易懂，而不是追求效率。

import math

a = 6378137  # semi major axis, WGS84
# f = 1 / 298.257223563 # WGS84
e2 = 6.69437999014e-3  # eccentricity, WGS84 e2 = 1 - (1-f)*(1-f)


def gps_to_geocentric_ecef(lat, lon, alt):
    # converts geodetic latitude, longitude, altitude based on WGS84 (GPS coordinates)
    # to geocentric x, y, z in meters (ECEF), sqrt(x^2+y^2+z^2) = geocentric radius
    # |z| distance to equatorial plane, |x| distance to the plane of the prime meridian
    # https://en.wikipedia.org/wiki/Geodetic_coordinates#Conversion
    # https://github.com/proj4js/proj4js/blob/master/lib/datumUtils.js#L68C5-L70C44
    # verify: https://tool-online.com/en/coordinate-converter.php (right select WORLD and XYZ (geocentric))
    # or: https://www.oc.nps.edu/oc2902w/coord/llhxyz.htm - paste lat, lon, height, then press "LLH to ECEF"

    lat = math.radians(lat)
    lon = math.radians(lon)
    sin_lat = math.sin(lat)
    cos_lat = math.cos(lat)
    cos_lon = math.cos(lon)
    sin_lon = math.sin(lon)

    rn = a / (math.sqrt(1 - e2 * sin_lat * sin_lat))
    x = (rn + alt) * cos_lat * cos_lon
    y = (rn + alt) * cos_lat * sin_lon
    z = ((rn * (1 - e2)) + alt) * sin_lat

    return [x, y, z]


def rotate_ecef_to_local_xyz(dx, dy, dz, lat, lon):
    # rotates the coordinate differences dx, dy, dz fromm the global frame (ECEF) to
    # a local frame that has the geodetic normal to (lat, lon) as z-axis, while
    # the new x-axis is towards North and the y-axis towards East

    lat = math.radians(lat)
    lon = math.radians(lon)

    cos_lon = math.cos(lon)
    sin_lon = math.sin(lon)
    cos_lat = math.cos(lat)
    sin_lat = math.sin(lat)
    dx, dy, dz = (- dx * cos_lon * sin_lat - dy * sin_lon * sin_lat + dz * cos_lat,
                  - dx * sin_lon + dy * cos_lon,
                  dx * cos_lon * cos_lat + dy * sin_lon * cos_lat + dz * sin_lat)

    return [dx, dy, dz]


def main():
    bgr_lat, bgr_lon, bgr_alt = [41.947694, 3.209083, 379.41]
    airplane_gps = [
        [41.9763015750,  3.4890851772, 971.32]
    ]

    # ecef coordinates for the station:
    bgr_xyz = gps_to_geocentric_ecef(bgr_lat, bgr_lon, bgr_alt)

    # ecef coordinates for the (first) airplane:
    plane_xyz = gps_to_geocentric_ecef(*airplane_gps[0])

    # xyz position of the airplane relative to the station, ecef orientation:
    diff_xyz = [plane_xyz[i] - bgr_xyz[i] for i in range(3)]

    # xyz position of the airplane relative to the station, station-local orientation
    local_xyz = rotate_ecef_to_local_xyz(*diff_xyz, bgr_lat, bgr_lon)

    print(local_xyz)  # final result, in meters (x-towards North, y-towards East, z-towards Zenith)


if __name__ == '__main__':
    main()

由于高度是以米为单位给出的，我也使用了WGS84的半长轴（也是以米为单位），结果也是以米为单位：这个平面距离北方3215米，距离东方23210米，高度为549米（这是从通过参考站BGR的切平面测量的高度）。

虽然我通过多种方式仔细验证了这段代码，但我还是希望没有人会基于它来飞飞机 :)

Answer 2

你可以使用pyproj这个库，并定义一下

import pyproj
import math

P = pyproj.Proj(proj='utm', zone=32, ellps='WGS84', preserve_units=True)
G = pyproj.Geod(ellps='WGS84')


base_x_offset = 19930.039590830565
base_y_offset = 4660215.067408236

def lonlat_to_xy(Lon, Lat):
    p = P(Lon, Lat)    
    return (p[0] - base_x_offset, p[1] - base_y_offset)    

def xy_to_lonlat(x,y):
    return P(x + base_x_offset, y + base_y_offset, inverse=True)    

def distance(Lat1, Lon1, Lat2, Lon2):
    return G.inv(Lon1, Lat1, Lon2, Lat2)[2]

然后用下面的代码来检查：

base_lon = 3.209083
base_lat = 41.947694

plane_lon = 3.4890851772
plane_lat = 41.9763015750

p_base = lonlat_to_xy(base_lon, base_lat)
print(p_base)

p_plane = lonlat_to_xy(plane_lon, plane_lat)
print(p_plane)

还有

xy_to_lonlat( lonlat_to_xy(plane_lon, plane_lat)[0] , lonlat_to_xy(plane_lon, plane_lat)[1] )

xy_to_lonlat( lonlat_to_xy(base_lon, base_lat)[0] , lonlat_to_xy(base_lon, base_lat)[1] )

不过，这个方法没有处理z坐标。而且在把一个球体转换成平面地图的时候，可能会出现很多问题。