File Exchange

image thumbnail

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

version 1.0 (1.41 KB) by

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

46 Downloads

Updated

No License

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

Comments and Ratings (20)

as ss

as ss (view profile)

What is lat is 0?

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

How do you handle the exception with div 0 ?

R.

peter pansen

works as it is supposed to! thanks!

is fa

zalmay amiri

b d

P. McNamara

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

mafalda sapo

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

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.

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

imran rao

nice

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?
 

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?

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

Glen Davidson

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

Jacq Mihof

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

Molly Kerr

Thank you!

Dave Shell

VERY USEFULL Program! THANK YOU!!!

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!

Yunus Levent Ekinci

very nice,thanks,good job

Updates

Corrected a numerical instability near exact poles; now accurate to about 2 millimeters.

MATLAB Release
MATLAB 7.0.1 (R14SP1)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video