Converting ECEF to WGS84
Mike Mason
I am trying to convert ECEF (x,y,z earth centered earth fixed)
coordinates to WGS84 coordinates. Does anyone know if MatLab has a
function that does this. If not how about some good recomendations on
books and websites or maybe someone will give me an example to look at.
Everthing that I have seen so far looks contradictory.
Thanks a lot.
Mike
Tom L. Davis
<Michael...@West.Raytheon.com> wrote:
Does this do what you want?
Tom
function p2 = xform( op, p1, c1, c2 )
%XFORM Transform between earth coordinate systems.
% xform( 'init', N ) allocate N NED coord systems
% xform( 'new', P1, C1, Cn ) new Cn NED at P1(C1)
% P2 = xform( 'pos', P1, C1, C2 ) P2(C2) == P1(C1)
% V2 = xform( 'vel', V1, C1, C2 ) V2(C2) == V1(C1)
% K2 = xform( 'cov', K1, C1, C2 ) K2(C2) == K1(C1)
%
% where Pi are point vectors (1x3), Vi are velocity vectors (1x3)
% and Ki are covariance matrices (3x3). Ci are coordinate system
% references with the following numbering scheme:
% Cn = -1 LLA p = [ lat lon alt ] in radians and meters
% = 0 ECEF p = [ x y z ] in meters
% = 1 NED # 1 p = [ x y z ] in meters
% = 2 NED # 2 p = [ x y z ] in meters
% = N NED # N p = [ x y z ] in meters
%
% For example, place the third NED coordinate system origin at
% N30 00 00.0, W90 00 00.0 and 100m elevation and then find
% the ECEF coordinates of a point 200 meters to the North:
% xform( 'new', [pi/6 -pi/2 100], -1, 3 );
% p = xform( 'pos', [200 0 0], 3, 0 );
global XFORM
switch op
% P2 = xform( 'pos', P1, C1, C2 ) P2(C2) <== P1(C1)
% Transform position P1 in coord system C1
% to position P2 in coord system C2
%
case 'pos'
% transform P1(C1) to ECEF
%
if c1 < 0, % p1 is LLA
p0 = lla2ecef(p1)';
elseif c1 == 0, % p1 is ECEF
p0 = p1';
else % p1 is NED
p0 = XFORM(c1).Cen * p1' + ... % XFORM(c1).org;
repmat( XFORM(c1).org, 1, size(p1,1) );
end
% transform P0 to C2
%
if c2 < 0, % p2 is LLA
p2 = ecef2lla(p0');
elseif c2 == 0, % p2 is ECEF
p2 = p0';
else % p2 is NED
p2 = XFORM(c2).Cen' * ...
( p0 - repmat( XFORM(c2).org, 1, size(p0,2) ) );
p2 = p2';
end
% V2 = xform( 'pos', V1, C1, C2 ) V2(C2) <== V1(C1)
% Transform velocity V1 in coord system C1
% to velocity V2 in coord system C2
% Velocities can be in any units
%
case 'vel'
% transform V1(C1) to ECEF
%
if c1 < 0,
error('LLA not valid for ''vel''');
elseif c1 == 0, % v1 is ECEF
p0 = p1';
else % v1 is NED
p0 = XFORM(c1).Cen * p1';
end
% transform P0 to C2
%
if c2 < 0,
error('LLA not valid for ''vel''');
elseif c2 == 0, % v2 is ECEF
p2 = p0';
else % v2 is NED
p2 = XFORM(c2).Cen' * p0;
p2 = p2';
end
% P2 = xform( 'cov', P1, C1, C2 ) P2(C2) <== P1(C1)
% Transform covariance P1 in coord system C1
% to covariance P2 in coord system C2
%
case 'cov'
% transform P1(C1) to ECEF
%
if c1 < 0,
error('LLA not valid for ''cov''');
elseif c1 == 0, % p1 is ECEF
p0 = p1;
else % p1 is NED
p0 = XFORM(c1).Cen * p1 * XFORM(c1).Cen';
end
% transform P0 to C2
%
if c2 < 0,
error('LLA not valid for ''cov''');
elseif c2 == 0, % p2 is ECEF
p2 = p0;
else % p2 is NED
p2 = XFORM(c2).Cen' * p0 * XFORM(c2).Cen;
end
% xform( 'new', P1, C1, Ci ) -- define new NED # Ci at P1(C1)
%
case 'new'
if isempty(XFORM),
xform init;
end
% transform P1(C1) to ECEF
%
if c1 < 0, % p1 is LLA
lat = p1(1);
lon = p1(2);
XFORM(c2).org = lla2ecef(p1)';
elseif c1 == 0, % p1 is ECEF
XFORM(c2).org = p1';
x = ecef2lla(p1);
lat = x(1);
lon = x(2);
else % p1 is NED
p0 = XFORM(c1).Cen * p1' + XFORM(c1).org;
XFORM(c2).org = p0;
lat = p0(1);
lon = p0(2);
end
% Transformation matrix:
% For column vectors: p.ecef = Cen * p.ned
% For row vectors: p.ecef = p.ned * Cen'
%
XFORM(c2).Cen = [ -sin(lat)*cos(lon), -sin(lon), -cos(lat)*cos(lon);
-sin(lat)*sin(lon), cos(lon), -cos(lat)*sin(lon);
cos(lat), 0, -sin(lat) ];
% xform( 'init', N ) -- allocate N NED coord systems
%
case 'init'
if nargin < 2,
p1 = 10; % default is 10 NED coord systems
end
XFORM = struct( 'org', cell(p1,1) ,'Cen', eye(3) );
[XFORM.org] = deal( zeros(3,1) );
xform( 'new', [0 0 0], -1, 1 );
otherwise
error( 'Bad op code' );
end
% LLA to ECEF using WGS-84
% from "Understanding GPS Principles and Applications"
% section 2.2.3
%
function ecef = lla2ecef( lat_lon_alt )
a = 6378.137e3; % earth radius @ equator
f = 1 / 298.257223563; % flattening
b = a * (1-f); % earth radius @ poles
e2 = 1 - (1-f)^2; % eccentricity squared
[n,m] = size( lat_lon_alt );
if m ~= 3,
error( 'Input must have three columns' );
end
ecef = zeros( n, m );
lat = lat_lon_alt(:,1);
lon = lat_lon_alt(:,2);
alt = lat_lon_alt(:,3);
ecef = [ a*cos(lon) ./ sqrt( 1 + (1-e2)*tan(lat).^2 ) +
alt.*cos(lon).*cos(lat), ...
a*sin(lon) ./ sqrt( 1 + (1-e2)*tan(lat).^2 ) +
alt.*sin(lon).*cos(lat), ...
a*(1-e2)*sin(lat) ./ sqrt( 1 - e2*sin(lat).^2 ) + alt.*sin(lat) ];
% ECEF to LLA using WGS-84
% from "Understanding GPS Principles and Applications"
% section 2.2.3
%
function lla = ecef2lla( ecef )
a = 6378.137e3; % earth radius @ equator
f = 1 / 298.257223563; % flattening
b = a * (1-f); % earth radius @ poles
e2 = 1 - (1-f)^2; % eccentricity squared
ep = sqrt((a/b)^2-1); % second eccentricity
[n,m] = size( ecef );
if m ~= 3,
error( 'Input must have three columns' );
end
lat_lon_alt = zeros( n, m );
x = ecef(:,1);
y = ecef(:,2);
z = ecef(:,3);
r = sqrt( x.^2 + y.^2 );
E2 = a^2 - b^2;
F = 54 * ( b*z).^2;
G = r.^2 + (1-e2) * z.^2 - e2 * E2;
c = ( e2^2 * F .* r.^2 ) ./ G.^3;
s = ( 1 + c + sqrt(c.^2 + 2*c) ).^(1/3);
P = F ./ ( 3 * (s + 1./s + 1).^2 .* G.^2 );
Q = sqrt( 1 + 2*e2^2 * P );
temp = a^2/2*(1+1./Q ) - (1-e2)*P.*z.^2 ./ (Q.*(1+Q)) - P.*r.^2/2;
r0 = -e2*P.*r ./ (1+Q) + sqrt( max(temp,0) );
U = sqrt( ( r - e2*r0).^2 + z.^2 );
V = sqrt( ( r - e2*r0).^2 + (1-e2)*z.^2 );
z0 = b^2 * z ./ (a*V);
lla = [ atan2( ( z + ep^2*z0 ), r ), ... % lat
atan2( y, x ), ... % long
U.*( 1 - b^2./(a*V) ) ]; % alt
return
% Ellipsoid name Semi-major axis Inverse flattenning
% (a) (1/f)
%
% Airy 6377563.3960 299.324964600
% Airy_modified 6377340.1890 299.324964600
% Australian_National 6378160.0000 298.250000000
% Bessel 1841 6377397.1550 299.152812800
% Bessel 1841 in Namibia 6377483.8650 299.152812800
% Clarke 1866 6378206.4000 294.978698200
% Clarke 1880 6378249.1450 293.465000000
% Everest (Sabah & Sarawak) 6377298.5560 300.801700000
% Everest 1830 6377276.3450 300.801700000
% Everest 1948 6377304.0630 300.801700000
% Everest 1956 6377301.2430 300.801700000
% Everest_Modified 6377304.0630 300.801700000
% GRS-80 6378137.0000 298.257222101
% Helmert 1906 6378200.0000 298.300000000
% Hough 6378270.0000 297.000000000
% International 6378388.0000 297.000000000
% Krassovsky 6378245.0000 298.300000000
% Modified Fisher 1960 6378155.0000 298.300000000
% SGS 85 6378136.0000 298.257000000
% South America 1969 6378160.0000 298.250000000
% WGS-72 6378135.0000 298.260000000
% WGS-84 6378137.0000 298.257223563
% e2 = 1 - (b/a)^2; % eccentricity squared
% f = 1 - b/a; % flattening
% ep = sqrt((a/b)^2-1); % second eccentricity
% ep = 0.0820944379496; % second eccentricity
% e2 = ep^2 / (ep^2+1); % eccentricity squared
% b = a / sqrt( ep^2 + 1 ); % earth radius @ poles
% invf = a / (a-b); % inverse flattening
Brian.Farrelly
Mike
I have functions that will convert between lat/lon and geocentric
coordinates and also a function for datum shift. The latter is
useful if you have positions in a datum other than WGS84. For
example in the offshore oil industry in Europe we use still
use ED50.
If you use these in conjunction with M_Map (a free and very good
mapping package) which can convert lat/lon to UTM for a range
of ellipsoids including the one used in WGS84 then you should
be OK.
See www2.ocgy.ubc.ca/~rich/map.html for M_Map.
What is the contradiction to which you are referring?
Brian
--
mailto:Brian.F...@nho.hydro.com Norsk Hydro Research Centre
phone +47 55 99 68 74 ((( Postboks 7190
fax +47 55 99 69 70 2oooS N-5020 Bergen
home +47 55 13 78 49 HYDRO Norway
rpb
You stated that you had functions to do coordinate transformations. Would
you mind posting them here, or @ least a URL to them? Many others (myself
included) would be interested in these transformations. Oh, and also thanks
for the link to M_Map - will check that out when I get time.
Thanking you in advance,
Robert Brookover
Mike Mason
I have been given some contradictory material to work with (basically
someone photocopied a bit of a chapter and never bothered to say from
which book). Anyway the chapter has some circular reasoning. It uses
height to find phi and then uses phi to find the height. I found some
better sources to work and M_Map should be very helpful.
Once again thanks,
Mike
Brian.Farrelly
Mike,
This may not be as contradictory as it looks at first glance!
What you are looking at is probably an iterative solution.
Using a starting value for h make a first estimate of phi.
Then using the estimate of phi make a refined estimate of h.
And so on. In some cases only one iteration is necessary and
then the logic may not be apparent.
Be careful about which height is being considered. What you
are being asked for is probably the height above the reference
ellipsoid. What we usually think of as height is height above mean
sea level. For most purposes height above the ellipsoid is the
sum of height above mean sea level and the height of the geoid above
the reference ellipsoid.
When you change datums the reference ellipsoid changes so the height
over the ellipsoid changes even though height above mean sea level
is unchanged.
Photocopies without the source are extremely frustrating. I have
learned from bitter experience to copy the title page and ensure that
it is attached to the relevant pages.
Brian.Farrelly
>
> Brian,
> You stated that you had functions to do coordinate transformations. Would
> you mind posting them here, or @ least a URL to them? Many others (myself
> included) would be interested in these transformations. Oh, and also thanks
> for the link to M_Map - will check that out when I get time.
>
Robert,
I will send you the routines plus an example by e-mail asap.
I don't want to post them on the group as the result will be
a bit too long.
If you don't get them within a week remind me. I am notoriously
forgetful.
Brian.Farrelly
>
> Brian,
> You stated that you had functions to do coordinate transformations. Would
> you mind posting them here, or @ least a URL to them? Many others (myself
> included) would be interested in these transformations. Oh, and also thanks
> for the link to M_Map - will check that out when I get time.
>
Robert,
I am having problems sending you these functions by e-mail.
I get error messages back from the address in your posting on
the group and suspect you have some anti-spam mechanism
in action!
Can you send me an e-mail to the address below so that I can
try again?
rpb
I sent you an email from one of my other accts. I will check w/ my IT ppl to see
what the problem is with my "main" work acct.
Anyway, I sincerely appreciate your efforts - which go above & beyond the MATLAB
newsgroup!!
Cheers,
Rob