6

I am trying to convert some coordinates geodetic coordinates (WGS84) to ECEF. I've tried my own implementation in Python and using pyproj and am getting difference results. Could someone please help me out as to where the issue is:

Here's my code:

import numpy as np

def lla_to_ecef_1(lat, lon, alt):
    # see http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html
    rad = np.float64(6378137.0)        # Radius of the Earth (in meters)
    f = np.float64(1.0/298.257223563)  # Flattening factor WGS84 Model
    cosLat = np.cos(lat)
    sinLat = np.sin(lat)
    FF     = (1.0-f)**2
    C      = 1/np.sqrt(cosLat**2 + FF * sinLat**2)
    S      = C * FF

    x = (rad * C + alt)*cosLat * np.cos(lon)
    y = (rad * C + alt)*cosLat * np.sin(lon)
    z = (rad * S + alt)*sinLat
    return x, y, z

def lla_to_ecef_2(lat, lon, alt):
    import pyproj
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    x, y, z = pyproj.transform(lla, ecef, lon, lat, alt, radians=False)
    return x, y, z

coords = [
  (37.4001100556,  -79.1539111111,  208.38),
  (37.3996955278,  -79.153841,  208.48),
  (37.3992233889,  -79.15425175,  208.18),
  (37.3989114167,  -79.1532775833,  208.48),
  (37.3993285556,  -79.1533773333,  208.28),
  (37.3992801667,  -79.1537883611,  208.38),
  (37.3992441111,  -79.1540981944,  208.48),
  (37.3992616389,  -79.1539428889,  208.58),
  (37.3993530278,  -79.1531711944,  208.28),
  (37.4001223889,  -79.1538085556,  208.38),
  (37.3992922222,  -79.15368575,  208.28),
  (37.3998074167,  -79.1529132222,  208.18),
  (37.400068,  -79.1542711389,  208.48),
  (37.3997516389,  -79.1533794444,  208.38),
  (37.3988933333,  -79.1534320556,  208.38),
  (37.3996279444,  -79.154401,  208.58),
]

for lat, lon, alt in coords:
  print lla_to_ecef_1(lat, lon, alt)
  print lla_to_ecef_2(lat, lon, alt)
  print ""
1
  • Maybe double-check your C and S equations. What you have as FF (1.0 - f)**2, I've usually seen as e2 where e2 is f*(2.0-f). C is 1.0/np.sqrt(1.0 - e2*sinlat*sinlat). Including the results of one of your test points using the 2 functions might be useful too.
    – mkennedy
    Feb 28, 2017 at 22:38

3 Answers 3

8

Here is a test(seems both methods work with same precision):

import math

import pyproj

coords = [
  (37.4001100556,  -79.1539111111,  208.38),
  (37.3996955278,  -79.153841,  208.48),
  (37.3992233889,  -79.15425175,  208.18),
  (37.3989114167,  -79.1532775833,  208.48),
  (37.3993285556,  -79.1533773333,  208.28),
  (37.3992801667,  -79.1537883611,  208.38),
  (37.3992441111,  -79.1540981944,  208.48),
  (37.3992616389,  -79.1539428889,  208.58),
  (37.3993530278,  -79.1531711944,  208.28),
  (37.4001223889,  -79.1538085556,  208.38),
  (37.3992922222,  -79.15368575,  208.28),
  (37.3998074167,  -79.1529132222,  208.18),
  (37.400068,  -79.1542711389,  208.48),
  (37.3997516389,  -79.1533794444,  208.38),
  (37.3988933333,  -79.1534320556,  208.38),
  (37.3996279444,  -79.154401,  208.58),
]

def gps_to_ecef_pyproj(lat, lon, alt):
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    x, y, z = pyproj.transform(lla, ecef, lon, lat, alt, radians=False)

    return x, y, z

def gps_to_ecef_custom(lat, lon, alt):
    rad_lat = lat * (math.pi / 180.0)
    rad_lon = lon * (math.pi / 180.0)

    a = 6378137.0
    finv = 298.257223563
    f = 1 / finv
    e2 = 1 - (1 - f) * (1 - f)
    v = a / math.sqrt(1 - e2 * math.sin(rad_lat) * math.sin(rad_lat))

    x = (v + alt) * math.cos(rad_lat) * math.cos(rad_lon)
    y = (v + alt) * math.cos(rad_lat) * math.sin(rad_lon)
    z = (v * (1 - e2) + alt) * math.sin(rad_lat)

    return x, y, z

def run_test():

    for pt in coords:
        print('pyproj', gps_to_ecef_pyproj(pt[0], pt[1], pt[2]))
        print('custom', gps_to_ecef_custom(pt[0], pt[1], pt[2]))


run_test()
4

np.cos() and np.sin() need 'radians' arguments. That's all.

print lla_to_ecef_1( np.radians( lat),  np.radians(  lon), alt)
0
3

The answer from @mrgloom is unfortunately now deprecated (part of the pyproj 2 upgrade)... here is another way that will work (not complete for the example from OP)

lat = 37.4001100556
lon = -79.1539111111
hae = 208.38
transformer = pyproj.Transformer.from_crs(
    {"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
    {"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
    )
x ,y, z = transformer.transform(lon,lat,hae,radians = False)
print(x, y, z)

954652.1879304583 -4982682.9811522495 3852889.356975911

more details: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1%20x%20,y,%20z%20=%20pyproj.transform(lla,%20ecef,%20posn_llh[1],posn_llh[0],posn_llh[2],radians%20=%20False)

Upgrading to pyproj 2 from pyproj 1 We recommended using the pyproj.transformer.Transformer and pyproj.crs.CRS in place of the pyproj.Proj and pyproj.transformer.transform().

Also see: Getting Started

Optimize Transformations

Warning

pyproj.transformer.transform() and pyproj.transformer.itransform() are deprecated.

pyproj 1 style:

> from functools import partial 
> from pyproj import Proj, transform
> proj_4326 = Proj(init="epsg:4326") 
> proj_3857 = Proj(init="epsg:3857")
> transformer = partial(transform, proj_4326, proj_3857) 
> transformer(12,
  1. pyproj 2 style:
> from pyproj import Transformer transformer =
> Transformer.from_crs("epsg:4326", "epsg:3857")
> transformer.transform(12, 12)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.