No BSD License  

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

4.0

4.0 | 15 ratings Rate this file 53 Downloads (last 30 days) File Size: 1.41 KB File ID: #7941

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

by

 

28 Jun 2005 (Updated )

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

| Watch this File

File Information
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)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (20)
07 Oct 2009 as ss  
25 Aug 2009 Ramon Nogueron

What is lat is 0?

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

How do you handle the exception with div 0 ?

R.

06 Aug 2008 peter pansen

works as it is supposed to! thanks!

11 May 2008 is fa  
30 Apr 2008 zalmay amiri  
03 Apr 2008 b d  
06 Nov 2007 P. McNamara

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

04 Oct 2007 mafalda sapo

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

08 Apr 2007 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.

19 Mar 2007 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

23 Feb 2007 imran rao

nice

24 May 2006 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?

07 Feb 2006 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?

27 Dec 2005 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

18 Dec 2005 Glen Davidson

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

05 Dec 2005 Jacq Mihof

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

01 Dec 2005 Molly Kerr

Thank you!

23 Sep 2005 Dave Shell

VERY USEFULL Program! THANK YOU!!!

13 Sep 2005 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!

26 Aug 2005 Yunus Levent Ekinci

very nice,thanks,good job

Updates
26 Apr 2006

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

Contact us