Direct conversion of latitude and height from one ellipsoid to another

I suggest a method for converting geodetic height and latitude from one oblate ellipsoid of revolution to another having the same center and symmetry axis. Unlike other approaches, the method here does not obtain height and latitude from Earth-centered, Earth-fixed (ECEF) Cartesian coordinates; this feature allows height conversion with high accuracy even in cases where the data format limits the precision of the latitude data. Height and latitude conversions can be expressed as a Fourier series in even multiples of latitude, with height changes having only cosines and latitude changes having only sines. The absolute difference of the flattenings of the two ellipsoids furnishes a simple upper bound on the maximum absolute latitude change. Conversions between the TOPEX and WGS84 ellipsoids, a practical necessity for the inter-mission comparison of satellite laser and radar altimeter data, illustrate the findings. Because the differences in the flattenings and semi-major axes of these ellipsoids are small, truncating the Fourier series after the term in twice the latitude gives an approximate conversion with an error less than 9 ×\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 10–12 radians of latitude and about 6 ×\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 10–6 m of height.


Introduction and motivation
Satellite laser and radar altimeter mission data report Earth's ocean and ice surface heights referred to various ellipsoids. The TOPEX/POSEIDON and Jason series of radar altimeters, and the ICESat laser altimeter, refer height to the TOPEX ellipsoid, while the Sentinel-3 radar and ICESat-2 laser altimeters refer height to the WGS84 ellipsoid. Height and latitude on one ellipsoid differ from height and latitude on the other, and the differences are functions of latitude because these two ellipsoids have different shapes as well as different sizes (Table 1). In order to use the current and historic constellation of altimeters to observe temporal changes in sea level and/or ice height, one must first put all the observations into the same geodetic ellipsoidal height system. Progress in climate science requires measurement of sea level rise rates with errors less than 0. . The conversion of height from one ellipsoid to another should therefore be made with an accuracy of at least 0.1 mm, and ideally much better than that.
The usual approach to height conversion (Kotsakis 2008), like more general conversions (Badekas 1969;Deakin 2004;Reit 1998, Soler 1976Watson 2006), moves from one geodetic system through Earth-centered, Earth-fixed (ECEF) Cartesian coordinates and then to the second geodetic system. The calculation of ECEF from geodetic coordinates is elementary, but calculating geodetic from ECEF requires either the solution of a quartic equation or iteration of a sequence of successive approximations for latitude, φ, and height, h [e.g., Bowring (1976), Fok andIz (2003), Fukushima (1999), Vermeille (2002), Zhu (1994)]. The process is computationally expensive, and unless an algorithm is particulary robust and intelligently coded, it may have errors that vary systematically with latitude (Fukushima 1999;. The accuracy of height and/or latitude conversions computed by moving from (h, φ) E1 to ECEF and then from ECEF to (h, φ) E2 (the subscripts indicate positions referred to ellipsoid E1 and ellipsoid E2) is limited. The magnitude of an ECEF vector is more than 6 × 10 10 larger than 0.1 mm, and computing small differences between large quantities is numerically inaccurate (Acton 1990). The accu- racy of the ECEF computed from (h, φ) E1 values given in the data may also be limited by the data product format. The TOPEX altimeter data product format (AVISO 1996), for example, rounds the geo-location information to the nearest integer 1 mm, which limits the accuracy of the height conversion that can be calculated via ECEF (Fig. 1).
The above considerations motivate the present paper. The investigation here applies when E1 and E2 are oblate ellipsoids of revolution having the same center and symmetry axis, so that ECEF x, y, z coordinates are invariant during the conversion from E1 to E2. In this case, longitude need not be considered; the useful ECEF coordinates are ( p, z), in which p x 2 + y 2 1/2 .

Definitions: ellipsoid specification, ECEF Cartesian and geodetic coordinates, orthogonality
An oblate ellipsoid of revolution is defined by two real positive constants, at least one of which must have length units. For purposes of illustration, this paper will take the WGS84 ellipsoid as E1 and the TOPEX ellipsoid as E2, as then the change in height, h h E2 − h E1 , is positive. The defining parameters of these ellipsoids are their semimajor axis, a, and reciprocal flattening, 1/ f (Table 1, AVISO 1996;NIMA 2000). The expressions here are written in terms of a and f (a − b)/a. The semi-minor axis, b, and the eccentricity, e a 2 − b 2 1/2 /a, may also be used where they simplify an expression. Define Treating p and z as constants one finds ∂h/∂φ 0, identically, because on any one ellipsoid h and φ are orthogonal coordinates. Orthogonality guarantees that on any one ellipsoid, infinitesimal changes in h and φ are independent, but it does not guarantee that the finite h and φ which occur while changing ellipsoids may be handled separately. Therefore, in the next section, I find a strict bound on| φ|.

Latitude conversion
The geocentric latitude of a point, θ arctan(z/ p), does not change as the specification of the ellipsoid changes. The ECEF coordinate equations on any given ellipsoid yield and if the point is on the ellipsoid this reduces to tan (θ ) 1 − e 2 tan(φ), in which case one has also tan(ρ) (1 − f ) tan(φ), and tan (θ ) (1 − f ) tan(ρ), where ρ is the parametric, or "reduced", latitude. All of these equations have the form tan(u K ) α J K tan(u J ), where u J and u K are different kinds of latitudes; as such, their difference u K − u J is zero at the equator and the poles, and therefore, the difference can be developed in a Fourier sine series in sin(2nu J ) with terms for positive integer n. This Fourier sine series is obtained as follows.
The right-most equality is obtained from the one on its left by multipling numerator and denominator by cos 2 (u J ), using double angle formulae, and setting β J K (α J K − 1)/ (α J K + 1). Then and so one has The order of the subscripts on α J K and β J K in the above indicates that they express latitude u K as a sine series in latitude u J . Reversing the roles of the "input" and "output" one has α K J 1/α J K , and β K J −β J K . Consequently, one may express one or the other of the differences ±(u K − u J ) as a series which is majorized by an alternating power series, and therefore truncating the series at the n th term gives the latitude difference with an error smaller than |β J K | n /n.
If φ E1 and φ E2 are the geodetic latitudes, referred to ellipsoids E1 and E2, of a point with geocentric latitude θ , then their difference can be expressed as using beta values β E1 and β E2 that express geodetic latitude on the subscripted ellipsoids in terms of a sine series in geocentric latitude. β E1 , and likewise for E2. At the altitude of the Sentinel-3 satellites κ ≈ 8/9, while if h is the height of sea level above or below the ellipsoid then κ ≈ 1, with an error less than 2 × 10 -5 . β ≈ e 2 /2, to first order in e 2 , and the alternating series argument yields a simple bound on |φ E2 − φ E1 |: As an example, f TOPEX − f WGS84 ≈ 2.5 × 10 -9 , and positive; therefore, the difference in latitude on the two ellipsoids is less than a few nanoradians, with TOPEX latitude more polar than WGS84 latitude in middle latitudes. A latitude change of this size represents a north-south displacement of less than 2 cm on the Earth's surface, which may be negligible for many applications. The displacement is an apparent one, caused by the change in flattening; the ( p, z) position of the data remains invariant as the ellipsoid changes.
In theory, the kappa values on the right-hand side depend, through N and h, on the φ values on the left-hand side, but for practical purposes, one may treat both κ E1 and κ E2 as constants and use only the first term in the series. The upper bound estimate, shown as the dashed curve in Fig. 2, is obtained with κ E1 κ E2 1. The solid curve in Fig. 2 is obtained with κ WGS84 1 and κ TOPEX a/ a + h , in which a is N TOPEX evaluated at 45˚latitude, a a/ 1 − e 2 /2, and h δ 0 is the mean value of the h from WGS84 to TOPEX. The solid line approximation differs from the total latitude change by less than 9 × 10 -12 radians, an invisible difference in Fig. 2. If this level of latitude error is tolerable, then one need not compute θ tan −1 (z/ p), as only sin(2θ ) ( pz)/ p 2 + z 2 is needed.

Height conversion
To second order the Taylor series expressing h h E2 −h E1 in terms of a a E2 − a E1 and f f E2 − f E1 has the form because ∂ n h/∂a n 0 for n > 1. The partial derivatives, to be evaluated at E1, are The S k factors are of order 1, so the first-order terms in h have magnitudes of order a and b f /2; the second-order terms are of order a f /4 and a( f ) 2 /8. For the conversion from WGS84 to TOPEX, the maximum magnitudes of these terms are about 0.7 m, 16 mm, 1.8 nm, and 5 pm, respectively. These are shown in Fig. 3, with one or two subscripts The first four terms in the Taylor series expansion of the height change in moving from the WGS84 to the TOPEX ellipsoid. The subscripts on h indicate the derivatives of the corresponding Taylor series term. The y axis labels mm, nm, pm indicate millimeters, nanometers, and picometers, respectively on h indicating the partial derivatives corresponding to that term, e.g. h a f a f ∂ 2 h/∂a∂ f . The Fourier series for h δ 0 + δ 1 cos(2φ) + δ 2 cos (4φ)+· · ·. The Fourier series for S 1 and S −1 each have an infinite number of terms with coefficients formed by recursions on complete elliptic integrals of the first and second kinds (Cvijovic, 2010). The Fourier series coefficients for h cannot be expressed analytically, but one may evalute h at a uniformly spaced sequence of latitudes and then compute the coefficients by least-squares. The first three coefficients obtained this way from E1 WGS84 and E2 TOPEX are given in Table 2. Figure 3 shows that three of the first four terms in the Taylor series expansion for h have a component in cos(2φ), and the Fourier series terms in higher frequencies are likely to be very small compared to the cos(2φ) term; therefore, one may expect that there is a good approximation in the form h ≈ δ 0 + δ 1 cos(2φ). It can be found most accurately by the least-squares approach above, but the following approach is nearly as accurate and yields an alternative expression.

Simple approximation of 1h
In the equatorial plane, z 0, S 1, and h(φ 0) (a E1 − a E2 ). On the polar axis, p 0, S (1 − f ), and h (|φ| π/2) (b E1 − b E2 ). Using these "boundary conditions" to set the two constants in h ≈ δ 0 + δ 1 cos(2φ) one obtains a result which can be rearranged into the form h (a E1 − a E2 )cos 2 φ +(b E1 − b E2 )sin 2 φ. An IDL script evaluating this equation is included in a set of IDL scripts furnished by the U.S. National Snow and Ice Data Center and known as the ICESat/GLAS Tools; according to the "README.txt" file included with the tools, this form was found empirically by Terry Horan. For the transformation from WGS84 to TOPEX, this form has a maximum error of 12.48 µm.

Notes on computation
Height conversions in this paper use S k ( f , φ) for k ∈ {−1, 1}, and if picometer accuracy is desired, then also for k −3. If the computing environment has "log1p" to evaluate log(1 + ξ ) when |ξ | 1, then S k exp klog1p −e 2 sin 2 φ /2 can be used for all k. If the environment has "hypot" to evaluate ξ 2 + η 2 1/2 , then S hypot [cosφ, (1 − f )sinφ] may be more accurate than S sqrt 1 − e 2 sin 2 φ . One should not compute S k ( f , φ) by expanding the square root in a binomial series in e 2 sin 2 φ as the series truncation error will increase with latitude. All three methods of computing S are available in the author's computing environment and give S −1 with relative error magnitude equal to 2 −52 , or about 2.22 × 10 -16 .
Above, approximations were compared to a "truth" obtained by conversions through ECEF as (φ, h) E1 → ( p, z) → (φ, h) E2 . The many algorithms for estimating (φ, h) from ( p, z) distribute error across latitude differently (Fukushima 1999), and error distribution in iterative methods depends on whether the stopping criterion examines h or p or z or φ (Fok and Iz 2003). I used the Matlab mapping toolbox routine "geodetic2ecef", and as I didn't know what algorithm it used, I also coded my own routines for the simple iteration and Bowring's (1976) methods, and checked that all three algorithms agreed to within 2 nm in h and 3 femtoradians in φ.

Concluding remarks
Transforming geodetic height and latitude coordinates from one ellipsoid to another by passing through ECEF coordinates can introduce error (Fig. 1), is computationally slow, and may distribute error unevenly across latitude (Fukushima 1999(Fukushima , 2006Fok and Iz 2003). I have presented a method for coordinate conversion that does not obtain (φ, h) from ( p, z), and allows the conversion of either φ or h to be done independently of the other. Each coordinate's conversion is expressible as a Fourier series in even multiples of latitude, with height changes having only cosines and latitude changes having only sines. The derivation of the sine series for latitude allows an upper bound on the maximum latitude change to be stated simply in terms of the flattenings of the two ellipsoids. In converting between modern Earth ellipsoids, it will be sufficient for most applications to use the Fourier series only as far as the term in twice the latitude. For the conversion between WGS84 and TOPEX, this gives the latitude conversion with an error less than 9 × 10 -12 radians and the height conversion within about 6 microns.