Convert Cartesian (ECEF) Coordinates to lat, lon, alt

Convert earth-centered, earth-fixed (ECEF) coordinates to latitude, longitude, and altitude.

Description

ECEF2LLA - convert earth-centered earth-fixed (ECEF) cartesian coordinates to latitude, longitude, and altitude.

USAGE:
[lat,lon,alt] = ecef2lla(x,y,z)

lat = geodetic latitude (radians)
lon = longitude (radians)
alt = height above WGS84 ellipsoid (m)
x = ECEF X-coordinate (m)
y = ECEF Y-coordinate (m)
z = ECEF Z-coordinate (m)

Notes: This function assumes the WGS84 model. Latitude is customary geodetic (not geocentric).

Michael Kleder, April 2006

• MATLAB release: MATLAB 7.0.1 (R14SP1)

[Yunus Levent Ekinci]

very nice,thanks,good job

[Jim Gottwald]

Both ecef2lla.m and lla2ecef.m are excellent algorithms. I tested these algorithms with the Satellite Tool Kit using LLA data for 8 cities around the world. Lat and long conversions are accurate to at least 4 decimal places. Altitude conversions are accurate to within 0.01 meters. ECEF conversions are generally accurate to within 0.03 meters.

Thanks for making these available!

[Dave Shell]

VERY USEFULL Program! THANK YOU!!!

[Molly Kerr]

Thank you!

[Jacq Mihof]

hi,i was wondering if anybody had Matlab code for a cartesian to polar conversion?

[Glen Davidson]

Do you have a reference for this calculation?
I thought this method was an approximation, and recursive methods are necessary.

[The Author]

Thanks for a very good question. There are both iterative and non-iterative solutions for the conversion from ECEF Cartesian to (lat,lon,alt), so I decided to choose a non-iterative technique. If you believe that LLA2ECEF is correct, then you can verify the accuracy of ECEF2LLA by running an experiment to verify that the two functions indeed invert each other over the entire sphere. (I obtain an angular disagreement of about 0.03 picoradians and an altitude disagreement of around 0.3 micrometers. I presume these are nonzero because of roundoff errors.)
maxerr=[0 0 0];
for iter=1:1e5
    t = rand*pi-pi/2;
    n = rand*2*pi;
    a = rand * 10000-5000;
    [x,y,z]=lla2ecef(t,n,a);
    [t2,n2,a2]=ecef2lla(x,y,z);
    dif=mod(abs([t-t2 n-n2 a-a2]),2*pi);
    maxerr=max([dif;maxerr]);
end
maxerr

[Jacq Mihof]

thanks for yer help,
was wpndering if i take in a signal, x(t) = I + jQ, to my system - is it possible to split the phase and envelope of this signal in two through a cartesian to polar conversion in matlab?

[Paul Connelly]

This is a great function, and could be dropped right in where i needed it. My 'amateurs' question is if i have lat and long but no altitude but seeing as I am on the sea use alt 0. I get a value for x,y,z. If I now move in x and y but have no new value for z, Is it valid to use the same z value in the new position?
 

[imran rao]

nice

[John Lansberry]

I haven't checked out this algorithm in detail, but the best algorithm I have ever found for numerical accuracy and efficiency is by Bowring (the following code is vectorized):

  function [latitude, longitude, altitude] = ecef2geo(posvececef, degreeflag)
%-------------------------------------------------------------------------------
% function [latitude, longitude, altitude] = ecef2geo(posvececef, degreeflag)
%
% references:
%
% bowring b., 1976. transformation from spatial to geographical coordinates,
% survey review 23, pp. 323-327.
%
% bowring, b.r. (1985) the accuracy of geodetic latitude and height equations,
% survey review, vol. 28, no. 218, pp. 202-206.
%
% inputs description
%
% posvececef earth-centered earth-fixed position vector (m).
% degreeflag specifies whether output latitude and longitude are in
% degrees.
%
% outputs description
%
% latitude geodetic latitude, positive north (deg or rad).
% longitude geodetic longitude, measured east from greenwich meridian (deg
% or rad).
% altitude geodetic altitude above wgs-84 ellipsoid (m).
%
% date version author change description
% -------- ------- ------------ ----------------------------------------------
% 11-16-00 1.00 je lansberry initial version
% 01-14-02 1.01 je lansberry converted from fortran95 to matlab and
% vectorized
% 05-05-03 1.02 je lansberry modified variable names to agree with bowring
% (1976)
% 10-22-03 1.03 je lansberry modified per bowring (1985)
%-------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% wgs-84 defining parameters.
%-------------------------------------------------------------------------------
  a = 6378137.0;
  f = 1.0 / 298.257223563;

%-------------------------------------------------------------------------------
% wgs-84 derived parameters.
%-------------------------------------------------------------------------------
  one_f = 1.0 - f;

  b = a * one_f; % semi-minor axis
  e2 = f * (2.0 - f); % first eccentricity squared
  epsilon = e2 / (1.0 - e2); % second eccentricity squared
  b_a = one_f;

%-------------------------------------------------------------------------------
% extract ecef components from input position vector.
%-------------------------------------------------------------------------------
  x = posvececef(:, 1);
  y = posvececef(:, 2);
  z = posvececef(:, 3);

%-------------------------------------------------------------------------------
% initialize outputs.
%-------------------------------------------------------------------------------
  latitude = zeros(size(x));
  longitude = latitude;
  altitude = latitude;

%-------------------------------------------------------------------------------
% quick check for all components zero.
%-------------------------------------------------------------------------------
  ii0 = (x == 0 & y == 0 & z == 0);

  if any(ii0),
    latitude(ii0) = 0;
    longitude(ii0) = 0;
    altitude(ii0) = 0;
  end

%-------------------------------------------------------------------------------
% quick calculations at poles.
%-------------------------------------------------------------------------------
  ii1 = (x == 0 & y == 0 & z ~= 0);

  if any(ii1),
    latitude(ii1) = sign(z(ii1)) * pi / 2;
    longitude(ii1) = 0;
    altitude(ii1) = abs(z(ii1)) - b;
  end

%-------------------------------------------------------------------------------
% quick calculations at equator.
%-------------------------------------------------------------------------------
  ii2 = (~ii0 & ~ii1 & z == 0.0);

  if any(ii2),
    longitude(ii2) = atan2(y(ii2), x(ii2));
    latitude(ii2) = 0;
    p = sqrt(x(ii2).^2 + y(ii2).^2);
    altitude(ii2) = p - a;
  end

%-------------------------------------------------------------------------------
% main algorithm. in bowring (1985), u is the parametric latitude. it is crucial
% to maintain the appropriate signs for the sin(u) and sin(lat) in the equations
% below.
%-------------------------------------------------------------------------------
  ii = ~ii0 & ~ii1 & ~ii2;

  if any(ii),

    p2 = x(ii).^2 + y(ii).^2;
    r2 = p2 + z(ii).^2;
    p = sqrt(p2);
    r = sqrt(r2);

%-------------------------------------------------------------------------------
% equation (17) from bowring (1985), shown to improve numerical accuracy in lat
%-------------------------------------------------------------------------------
    tanu = b_a * (z(ii) ./ p) .* (1 + epsilon * b ./ r);
    tan2u = tanu .* tanu;

%-------------------------------------------------------------------------------
% avoid trigonometric functions for determining cos3u and sin3u
%-------------------------------------------------------------------------------
    cos2u = 1.0 ./ (1.0 + tan2u);
    cosu = sqrt(cos2u);
    cos3u = cos2u .* cosu;

    sinu = tanu .* cosu;
    sin2u = 1.0 - cos2u;
    sin3u = sin2u .* sinu;

%-------------------------------------------------------------------------------
% equation (18) from bowring (1985)
%-------------------------------------------------------------------------------
    tanlat = (z(ii) + epsilon * b * sin3u) ./ (p - e2 * a * cos3u);

    tan2lat = tanlat .* tanlat;
    cos2lat = 1.0 ./ (1.0 + tan2lat);
    sin2lat = 1.0 - cos2lat;

    coslat = sqrt(cos2lat);
    sinlat = tanlat .* coslat;

    longitude(ii) = atan2(y(ii), x(ii));
    latitude(ii) = atan(tanlat);

%-------------------------------------------------------------------------------
% equation (7) from bowring (1985), shown to be numerically superior to other
% height equations. note that equation (7) from bowring (1985) writes the last
% term as a^2 / nu, but this reduces to a * sqrt(1 - e^2 * sin(lat)^2), because
% nu = a / sqrt(1 - e^2 * sin(lat)^2).
%-------------------------------------------------------------------------------
    altitude(ii) = p .* coslat + z(ii) .* sinlat - a * sqrt(1.0 - e2 * sin2lat);

  end

% longitude = unwrap(longitude);

%-------------------------------------------------------------------------------
% convert outputs if necessary.
%-------------------------------------------------------------------------------
  if nargin == 2 & degreeflag == 1,
    radtodeg = 180 / pi;
    latitude = latitude * radtodeg;
    longitude = longitude * radtodeg;
  end

[Sajan Mushini]

Its not giving precise values. for minute changes in X,Y and Z, its not able to get the precise values. Hope you can look into this. I think, we should change some values to double float.

[mafalda sapo]

this program is a piece of geat!!!!!!!!!!!!

[P. McNamara]

Be careful. Also we are looking at your downloads records.

[b d]

 

[zalmay amiri]

 

[is fa]

 

[peter pansen]

works as it is supposed to! thanks!

[Ramon Nogueron]

What is lat is 0?

alt = p./cos(lat)-N;

How do you handle the exception with div 0 ?

R.

[as ss]