python 坐标系转化

python ECEF转GPS

• 代码如下：

def ecef2gps(point_3d):
    # pyproj库函数
    # ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    # lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    # lon, lat, alt = pyproj.transform(ecef, lla, point_3d[0], point_3d[1], point_3d[2], radians=False)

    # 公式计算
    WGS84_A = 6378137.0
    WGS84_E = 0.0818191908

    b = (WGS84_A * WGS84_A * (1 - WGS84_E * WGS84_E))  ** 0.5
    ep = ((WGS84_A * WGS84_A - b * b) / (b * b)) ** 0.5
    p = math.hypot(point_3d[0], point_3d[1])
    th = math.atan2(WGS84_A * point_3d[2], b * p)
    lon = math.atan2(point_3d[1], point_3d[0])
    lat = math.atan2((point_3d[2] + ep * ep * b * (math.sin(th) **3 )), (p - WGS84_E * WGS84_E * WGS84_A * (math.cos(th) ** 3)))
    N = WGS84_A / ((1 - WGS84_E * WGS84_E * math.sin(lat) * math.sin(lat)) ** 0.5)
    alt = p / math.cos(lat) - N
    lat = RAD2DEG(lat)
    lon = RAD2DEG(lon)

    print('--------------------gps---------------------')
    print(lat,lon,alt)
    print('--------------------------------------------')

    return lat,lon,alt

python GPS转ECEF

• 代码如下：

def gps2ecef(gps):
    WGS84_A = 6378137.0
    WGS84_B = 6356752.314245

    clat = math.cos(deg2rad(gps[0]))
    slat = math.sin(deg2rad(gps[0]))
    clon = math.cos(deg2rad(gps[1]))
    slon = math.sin(deg2rad(gps[1]))
    a2 = WGS84_A * WGS84_A
    b2 = WGS84_B * WGS84_B

    L = 1.0 / ((a2 * clat * clat + b2 * slat * slat) ** 0.5)
    x = (a2 * L + gps[2]) * clat * clon
    y = (a2 * L + gps[2]) * clat * slon
    z = (b2 * L + gps[2]) * slat

    ecef = [x, y, z]

    print('------------------ecef--------------------------')
    print(ecef)
    return ecef

• 弧度角度转化

def rad2deg(x):
    return x * 180.0 / math.pi


def deg2rad(x):
    return x * math.pi / 180.0