Simple WGS 84 - ECEF conversion functions.

One of the most unexpectedly complicated things I have found during my research into determining the orbit of a meteor using Doppler measurements is the conversion between geographical coordinates and rectangular, or Cartesian coordinates, also called ECEF (Earth Centered Earth Fixed) coordinates:

At first sight it seems easy: The earth is a sphere, so it only need a few sines and cosines. But later, when you need more precision, you discover earth is not a sphere, it is an ellipsoid, so things start to become problematic.
Later, you discover the ellipsoid you are using is not the ellipsoid used by the rest of the world. Most of the world uses WGS 84 (including GPS) and your coordinates are GPS based, so you need to use the right formulas to do the conversion.

There is a lot of literature on Internet about this, but as soon you start to read you discover a big problem. Conversion from longitude, latitude and altitude is easy and straightforward, for example in Octave code:

% WGS84 Lon-Lat-Alt coordinates to earth centered X-Y-Z
function xyz = lla2xyz(lla)
  a=6378137;            % radius
  e=8.1819190842622e-2; % eccentricity
  asq=a*a;
  esq=e*e;
  lon=lla(1);
  lat=lla(2);
  alt=lla(3);
  N=a/sqrt(1-esq*sin(lat)^2);
  x=(N+alt)*cos(lat)*cos(lon);
  y=(N+alt)*cos(lat)*sin(lon);
  z=((1-esq)*N+alt)*sin(lat);
  xyz=[x y z];
endfunction;

Note: I wrote the functions in Octave/Matlab code without advanced functions to make the task of porting to other languages easy. For example, many sqrt instructions can be written in a more clear way using Octave's norm function.

Surprisingly, WGS 84 has no inverse model: you can calculate Cartesian coordinates easily from WGS 84 latitude, longitude and altitude. But going from Cartesian coordinates to WGS 84 latitude, longitude and altitude is very tricky.

There are two methods to make the conversion. The first one involves making iterations, so conversion precision is function of the number of iterations and convergence criteria. I want a simple and fast conversion, so I discarded this option.

The other method is a simple algorithm:

% X-Y-Z earth centered coordinates to WGS 84 Lon-Lat-Alt

function lla=xyz2lla(xyz)
  a=6378137;            % radius
  e=8.1819190842622e-2; % eccentricity
  asq=a*a;
  esq=e*e;
  x=xyz(1);
  y=xyz(2);
  z=xyz(3);
  b=sqrt(asq*(1-esq));
  bsq=b*b;
  ep=sqrt((asq-bsq)/bsq);
  p=sqrt(x*x+y*y);
  th=atan2(a*z,b*p);
  lon=atan2(y,x);
  lat=atan2((z+ep*ep*b*sin(th)^3),(p-esq*a*cos(th)^3));
  N=a/(sqrt(1-esq*sin(lat)^2));
  alt=p/cos(lat)-N;
  lla=[lon lat alt];
endfunction;

Well, I just said there is no inverse model, and now I show you a simple algorithm to make the conversion. What is the catch? The catch is this algorithm only works for points situated at altitude = 0, or really close to 0. If you play with it, you will find out it works fine finding latitude and longitude, but fails with altitude.

The solution I found to solve this problem is:

• Get longitude, latitude and altitude with xyz2lla using desired x, y, z values.
• Discard altitude and get ground coordinates of calculated latitude and longitude using lla2xyz(lon,lat,0)
• Call this new ground Cartesian coordinates gx, gy, gz
• Get the module of vectors [x y z] and [gx gy gz].
• Altitude is the distance between points x,y,z and gx,gy,gz. Subtract length of vector [gx gy gz] to length of vector [x y z] to get altitude.

This is an obvious solution, but I didn't found it on the Internet, so excuse me if I have just reinvent the wheel.

In Octave code can be implemented in this way:

% Earth centered X-Y-Z coordinates to WGS84 Lon-Lat-Asl

function lla=xyz2lla(xyz)
  a=6378137;            % radius
  e=8.1819190842622e-2; % eccentricity
  asq=a*a;
  esq=e*e;
  x=xyz(1);
  y=xyz(2);
  z=xyz(3);
  b=sqrt(asq*(1-esq));
  bsq=b*b;
  ep=sqrt((asq-bsq)/bsq);
  p=sqrt(x*x+y*y);
  th=atan2(a*z,b*p);
  lon=atan2(y,x);
  lat=atan2((z+ep*ep*b*sin(th)^3),(p-esq*a*cos(th)^3));
  N=a/(sqrt(1-esq*sin(lat)^2));
  %alt=p/cos(lat)-N; % Not needed here
  g=lla2xyz([lon lat 0]);
  gm=sqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3));
  am=sqrt(x*x+y*y+z*z);
  alt=am-gm;
  lla=[lon lat alt];
endfunction;

To test this function, I wrote this small Octave script to create 100000 conversions. It generates points at random latitude, longitude and altitude, converts them to X-Y-Z values, and then, it uses the X-Y-Z values to calculate latitude, longitude and altitude using my modified xyz2lla function.. At the end, a small report is presented:

#!/usr/bin/octave -q
function r=deg2rad(d)
  r=d*pi/180;
endfunction;

function d=rad2deg(r)
  d=r*180/pi;
endfunction;

% WGS84 Lon-Lat-Alt coordinates to earth centered X-Y-Z
function xyz = lla2xyz(lla)
  a=6378137;            % radius
  e=8.1819190842622e-2; % eccentricity
  asq=a*a;
  esq=e*e;
  lon=lla(1);
  lat=lla(2);
  alt=lla(3);
  N=a/sqrt(1-esq*sin(lat)^2);
  x=(N+alt)*cos(lat)*cos(lon);
  y=(N+alt)*cos(lat)*sin(lon);
  z=((1-esq)*N+alt)*sin(lat);
  xyz=[x y z];
endfunction;

% Earth centered X-Y-Z coordinates to WGS84 Lon-Lat-Asl
function lla=xyz2lla(xyz)
  a=6378137;            % radius
  e=8.1819190842622e-2; % eccentricity
  asq=a*a;
  esq=e*e;
  x=xyz(1);
  y=xyz(2);
  z=xyz(3);
  b=sqrt(asq*(1-esq));
  bsq=b*b;
  ep=sqrt((asq-bsq)/bsq);
  p=sqrt(x*x+y*y);
  th=atan2(a*z,b*p);
  lon=atan2(y,x);
  lat=atan2((z+ep*ep*b*sin(th)^3),(p-esq*a*cos(th)^3));
  N=a/(sqrt(1-esq*sin(lat)^2));
  %alt=p/cos(lat)-N; % Not needed here
  g=lla2xyz([lon lat 0]);
  gm=sqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3));
  am=sqrt(x*x+y*y+z*z);
  alt=am-gm;
  lla=[lon lat alt];
endfunction;

iterations=100000;
sumlon=0;
sumlat=0;
sumalt=0;
for i=1:iterations
  lon=-180+rand*360; % -180...+180
  lat=-90+rand*180;  % -90...+90
  alt=-500000+rand*1500000; % -500...+1000 km
  la=deg2rad(lat);
  lo=deg2rad(lon);
  c=lla2xyz([lo la alt]);
  lla=xyz2lla(c);
  printf('Input :  lon=%f lat=%f alt=%f\n',lon,lat,alt);
  printf('Output:    x=%f   y=%f   z=%f\n',c(1),c(2),c(3));
  printf('Reverse: lon=%f lat=%f alt=%f\n',rad2deg(lla(1)),rad2deg(lla(2)),lla(3));
  sumlon=sumlon+(lla(1)-lo)^2;
  sumlat=sumlat+(lla(2)-la)^2;
  sumalt=sumalt+(lla(3)-alt)^2;
  printf('Altitude error in xyz2lla: %f meters\n',lla(3)-alt);
  printf('\n');
endfor;
rmslon=sqrt(sumlon/iterations);
rmslat=sqrt(sumlat/iterations);
rmsalt=sqrt(sumalt/iterations);
printf('Sumary:\n');
printf('Number of iterations: %d\n',iterations);
printf('RMS error in longitude: %f radians, %f degrees.\n',rmslon,rad2deg(rmslon));
printf('RMS error in latitude:  %f radians, %f degrees.\n',rmslat,rad2deg(rmslat))
printf('RMS error in altitude:  %f meters\n',rmsalt);

And the output of this little script in my system is:

As you can see, the RMS error in altitude is a little bit under 1.6 meters, not bad for such simple function.

Some interesting readings about this:

https://geodesy.curtin.edu.au/local/docs/1Featherstone295.pdf

https://microem.ru/files/2012/08/GPS.G1-X-00006.pdf

https://gist.github.com/klucar/1536194 the functions I ported to Octave.

I hope this helps to someone who need it.

Miguel A. Vallejo, EA4EOZ