# Converting WGS84 to ECEF in Python?

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
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 ""
``````
• 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. Feb 28, 2017 at 22:38

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)

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()
``````

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

``````print lla_to_ecef_1( np.radians( lat),  np.radians(  lon), alt)
``````

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

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)
``````