Download Transforming Cartesian Coordinates to Geographical Coordinates: A Comprehensive Guide and more Summaries Mathematics in PDF only on Docsity!
Transforming Cartesian coordinates X,Y,Z
to Geographical coordinates φ , λ , h
George P. Gerdan and Rodney E. Deakin Department of Land Information Royal Melbourne Institute of Technology GPO Box 2476V, MELBOURNE, VIC, 3001 AUSTRALIA
[published in The Australian Surveyor, Vol. 44, No. 1, pp. 55-63, June 1999]
ABSTRACT
In November 1994, the Intergovernmental Committee on Surveying and Mapping (ICSM) adopted a new geodetic datum for Australia and recommended its progressive implementation nationally. This datum, designated Geocentric Datum of Australia (GDA) – compatible with the Global Positioning System (GPS), a US Department of Defense satellite positioning system – will eventually replace the existing Australian Geodetic Datum (AGD). Consequently, “old” AGD coordinates need to be converted to “new” GDA coordinates, which
may require the transformation of Cartesian coordinates X,Y,Z to geographical coordinates φ , λ, h (latitude,
longitude, height) related to an ellipsoid of revolution.
Given X,Y,Z , longitude is easily derived, but not so latitude, which requires more sophisticated evaluation – usually by iterative techniques, or complicated direct methods – the height following readily once latitude is obtained. Thus Cartesian-to-geographic transformations revolve around the determination of latitude; this paper reviews published techniques, some quite recent, which may be of use to practitioners.
INTRODUCTION
Conversion from geographic to Cartesian coordinates is a simple operation since closed formulae relate X,Y,Z to
φ , λ, h. Inverse computations, on the other hand, are more complex, since no simple relationships link φ to
X,Y,Z. Various techniques for computing φ are well documented in the literature, some more efficient than
others when comparing computational time; indirect (or iterative) techniques are often preferred for their relative simplicity when compared with more complicated direct methods. This paper compares six techniques; two direct, Paul (1973) and Ozone (1985) and four indirect, Bowring (1976), Borkowski (1989), Lin & Wang (1995) and Simple Iteration. The last is described in several geodesy texts, for instance Bomford (1980, p. 679).
These are not the only techniques for computing φ given X,Y,Z; Borkowski (1989) refers to fourteen other
methods; three direct, five iterative, four using closed approximate formulae cast as series expansions and two others, as well as another direct method of his own (Borkowski 1987). Vincenty (1985) claims that the earliest solution to this problem was given by Doerrie in his book Kubische und biquadratische Gleichungen, Muenchen, 1948.
This review is not definitive, since only a small selection of the published techniques are chosen, but amongst them, Bowring’s method and Simple Iteration have often been used for comparative analysis. For instance, Borkowski (1989) compares his iterative method against nine others, including three variants of Simple Iteration, concluding that whilst simple iterative procedures are the easiest to program, his solution is more appropriate – requiring fewer iterations to reach acceptable accuracy. Laskowski (1991) compares Borkowski’s and Bowring’s methods, demonstrating that Bowring’s is faster whilst Lin & Wang (1995) find their method to be faster than Bowring’s. Lin & Wang mention a review of methods by Rapp (1984) who concluded (at the time) that Bowring’s procedure was the most efficient. Thus from the literature, it could be concluded that Lin & Wang’s method is the most efficient Cartesian-to-geographic transformation. This is in fact confirmed by computational analysis of the six techniques, which are regarded as representative of the numerous published methods.
GEODETIC COORDINATE SYSTEMS IN USE IN AUSTRALIA
In 1966, under the direction of the National Mapping Council (NMC) all geodetic surveys in Australia were recomputed and adjusted on the then new AGD, an astronomically derived topocentric datum having a physical origin near the centroid of the geodetic network and fixing an ellipsoid of revolution, the Australian National Spheroid (ANS), with respect to the Earth’s rotational axis. The national adjustment yielded an homogeneous set of geographical coordinates (latitudes and longitudes) for the geodetic network. At the same time, the NMC defined a system of rectangular grid coordinates (eastings and northings) known as the Australian Map Grid (AMG), based on a Universal Transverse Mercator (UTM) projection of AGD latitudes and longitudes.
After 1966 there were several readjustments of the national geodetic network, densified and strengthened by the inclusion of improved measurements, each readjustment referred to as a Geodetic Model of Australia (GMA). In 1984 the NMC, recognizing the eventual need for Australia to convert to a geocentric datum, adopted the latest readjustment at the time, GMA82, as an interim step in this process. This geographical coordinate set was defined as AGD84 with AMG84 grid coordinates, and to avoid confusion, earlier coordinate sets derived from the 1966 adjustment were defined as AGD66 and AMG66. Both ADG66 and AGD84 coordinates have a common datum (defined in 1966) excepting that AGD84 coordinates were derived from an adjustment which more correctly allowed for the separation between the geoid and the ANS over Australia (NMC 1986).
In 1988, the NMC was superseded by the ICSM, representing the mapping organizations of the States and Territories of the Commonwealth of Australia and New Zealand. The GDA was adopted by the ICSM in November 1994 in response to anticipated demand by major users of GPS technology such as the Australian Defence Force, the International Civil Aviation Organization, the International Hydrographic Organization and the International Association of Geodesy (Steed 1996, p. 24). The new datum is primarily based on the coordinates of eight geologically stable sites across Australia with permanent GPS tracking facilities known as the Australian Fiducial Network (AFN), supplemented by a network of seventy survey stations (covering Australia at approximately 500km intervals) which together form the Australian National Network (ANN). Geocentric Cartesian coordinates of these stations were derived from an adjustment of precise GPS observations obtained from – (i) a two week global observation period in 1992 conducted by the International GPS Geodynamics Service at approximately two hundred sites around the world (including all the AFN sites) and (ii) ICSM campaigns in 1992, ’93 and ’94 linking all AFN and ANN sites. These coordinates are related to the International Earth Rotation Service (IERS) Terrestrial Reference Frame for 1992 (ITRF92) at epoch 1994. [The epoch 1994.0 (1st^ Jan. 1994) reflects the fact that monitoring stations used by IERS are moving with respect to each other due to earth crustal motion; the epoch date indicating the datum is ITRF92 adjusted for station motion in the intervening period]. The ICSM has defined GDA94 coordinates as latitudes and longitudes related to the ellipsoid of the Geodetic Reference System 1980 (GRS80) [BG 1988, p. 348] and Map Grid Australia 1994 (MGA94) grid coordinates as a UTM projection of those latitudes and longitudes.
The GRS80 was adopted by the International Union of Geodesy and Geophysics (IUGG) at its XVIITH^ General Assembly in Canberra, December 1979, as the best representation of the size and shape of the Earth. The GRS80 is based on the theory of a geocentric equipotential ellipsoid and is defined by four physical parameters of the Earth – (i) a the equatorial radius, (ii) GM the geocentric gravitational constant, (iii) J 2 the dynamical
form factor, and (iv) ω the angular velocity – from which geometric constants, e (square of the first
eccentricity) and f (flattening) are derived. The World Geodetic System 1984 (WGS84), the datum for GPS, is also based on GRS80, except that the dynamical form factor of the Earth is expressed in a modified form, causing very small differences between derived constants of the WGS84 and GRS80 ellipsoids. These differences can be regarded as negligible for all practical purposes and GDA94 is considered compatible with WGS84.
2
Geometric parameters of relevant ellipsoids are: Ellipsoid semi-major axis a flattening f semi-minor axis b = a ( 1 − f ) ANS 6378160 m 1/298.25 6356774.7192 m GRS80 6378137 m 1/298.257222101 6356752.3141 m WGS84 6378137 m 1/298.257223563 6356752.3142 m
Note 3: Sets of transformation parameters ( ds, R (^) X , RY , RZ and X (^) T , YT , ZT ) for AGD66 and AGD84 are
available from AUSLIG.
The Cartesian-to-geographic conversion problem is embodied in (1.8) where functions of φ ( tan φ ,ν and sin φ )
appear on both sides of the equation making the evaluation of φ difficult. The following techniques, direct and
iterative, are solutions to this problem.
normal
Gr
ee
n w
ic
h
M
eri
di
an
Equ^ ato^ r
E
N
P
X
Z
O Y
H
φ
λ
M
a
P h
D
E
a
C
Q
b
ellipsoid
G
R
Figure 1 Section of ellipsoid
P
O
H
normal
N
M
D
a
auxiliary circle
ν
Z
C
P
E
PC
h Q
S
φ
θ ψ (^) R T
Figure 2 Meridian ellipse
ELLIPSOID RELATIONSHIPS
Referring to Figures 1 and 2, the following relationships will be of use in the derivations that follow.
(i) OE = OG = a and ON = b are the semi-major and semi-minor axes of the ellipsoid (ii) PH is the ellipsoidal normal and P PE = h. PE is the projection of P on the ellipsoid and PC is the projection of PE on the auxiliary circle of the meridian ellipse
(iii) H PE = ν and C PE = ρ are the radii of curvature in the prime vertical and meridian planes of P
respectively
(iv) OH =ν e^2 sin φand DH =ν e^2
(v) latitudes: PDM = φ ( geodetic ) , PC O M = ψ ( parametric ), PE O M = θ ( geocentric )
(vi) latitudes are related by: tan ψ= tanφ
b a
and tan θ= tan ψ= tan
b a
b a
2
(vii) PS = p = X^2 + Y^2 and PT = Z
(viii) the ratio PC R PE R = a b and PE Q = p E = a cos ψ , PE R = Z E = b sin ψ
(ix) perpendicular distances from C (the centre of curvature of the meridian ellipse) to the Z- and OM-axes respectively are: a b a
a e a
2 2 2
F − 3 2 3
HG^
I
KJ^
cos ψ= cosψ and −
F −
HG^
I
KJ^
a b b
b e b
2 2 2
sin^3 ψ^2 sin^3 ψ
The following relations between geometric constants of the ellipsoid are also of use.
f
a b a
flattening (2.1) b = a (1 − f ) semi-minor axis (2.2)
e a b a
(^2) f f
2 2 = 2 2
= ( − ) first eccentricity squared (2.3)
e a b b
e e
2
2 2 2
2
1 2 second eccentricity squared^ (2.4)
SIMPLE ITERATION
This technique, described in various forms in geodesy texts, owes its popularity to its programming simplicity. The basis for the method is (1.8)
tan sin
Z + e p
An approximate value φ 0 is used in the right hand side (RHS) of (1.8) to evaluate tan φ 1 , (and hence φ 1 ) on the
left hand side (LHS). This “new” value, φ 1 is then used in the RHS to give the next value, tan φ 2. This
procedure is repeated until the difference between successive LHS values, tan φ i , tan φ i + 1 , reaches an acceptable
limit; the iteration thus converges to a solution of tan φ.
The starting value φ 0 is obtained from tan θ ′ =
Z
p
, θ ′being a good approximation of the geocentric latitude θ.
Hence, since tan φ= tan θ≈ tan ′
a b
a b
2 2
2
2 θ^ ,^ using (2.4) gives
tan
Z + e ′ p
After evaluating φ the height h is calculated from (1.9).
In common with many direct solutions, the method fails when Z approaches zero, in which case q and ρ 1 both
approach unity, t 1 becomes infinitely small and ζ is indeterminate. In such cases, Paul (1973, p. 136) gives an
approximate equation for the calculation of tan φ.
OZONE’S METHOD
This method (Ozone 1985) is direct in so far as tan φ is obtained from a simple closed equation, but only after
several intermediate variables have been evaluated. An outline of the development of the necessary equations is given below.
From the ellipsoid relationships, where subscripts E refer to points on the ellipsoid
tan φ = tan
Z Z
p p
a b
E E
which, with p E = a cos ψ and Z E = b sin ψ can be rearranged to give
ap b Z a b
cos ψ sinψ
− = 2 −^2 (5.2)
Making the substitution tan
u
and using trigonometric identities leads to
cos ψ =
u u
2 2
, sin ψ =
u u
, tan ψ =
u u
and substituting expressions for cos ψ and sin ψ into (5.2) gives
ap u u
b Z u u
a b
2
2 (^12 ) 1
which can be simplified to give a quartic equation in u u^4 − 4 M u^3 − 4 N u − 1 = 0 (5.4)
where M
a p a b b Z
− ( 2 −^2
N
a p a b b Z
+ ( 2 −^2
Introducing variables I, J, K enables (5.4) to be expressed in two forms (i) the product of two quadratic equations u^2 − ( 2 M − J u ) + ( I + J ) u^2 − ( 2 M + J u ) + ( I − J ) = 0
(ii) the difference of two squares u^2 M u I J u K (5.8)
(^2 )
d −^2 +^ i −^ b +^ g =
Both (5.7) and (5.8) give
u^4 − 4 M u^3 + ( 2 I − J^2 + 4 M^2 ) u^2 − ( 2 K J + 4 M I u ) + ( I^2 − K^2 = (5.9)
Comparing coefficients of u and u^2 , and the constant terms of (5.4) and (5.9) gives J^2 = 2 I + 4 M^2 , K^2 = I^2 + 1 , J K = 2 ( N − MI ) (5.10) which can be rearranged to give a cubic equation in I
I^3 + V I − W = 0 (5.11) for which standard solutions exist for the real root of I (Spiegel 1968, p. 32)
I
V W W V W W
= F
H
I
K
+ F
H
I
K
F
H
GG
I
K
JJ −^
F
H
I
K
+ F
H
I
K
F
H
GG
I
K
JJ
3 2 3 2
(^13 ) (5.12)
where V = 4 N M + 1 (^2) )
W = 2 ( N^2 − M (5.14)
After solving for I, J and K follow from (5.10)
J = 2 I + 4 M 2 (5.15)
K
N MI
J
The parameter u can now be obtained from the real solution of (5.7), a quadratic equation in u for which the positive root is
u
M J G
where G = ( 2 M + J ) 2 − 4 ( I − K ) (5.18)
Now since tan φ = tan
a b
ψ and from (5.3), where tan ψ =
u u
then
tan ( )
a u b u
Thus, having X,Y,Z (hence p ) for a point related to an ellipsoid, tan φ is obtained by computing the following
variables in order: M and N from (5.5) and (5.6), noting that Z is taken as positive, V and W from (5.13) and
(5.14), I from (5.12), J and K from (5.15) and (5.16), G from (5.18), then u from (5.17) and finally tan φ
from (5.19).
After calculating φ , its sign is determined by inspection of Z.
BOWRING’S METHOD
This method (Bowring 1976) is iterative, but owing to the nature of the equation used can be regarded as exact, since for all practical purposes, second or third iterations are not required (Bowring 1976, p.326).
From the ellipsoid relationships, perpendicular distances from C (the centre of curvature of the meridian ellipse)
to the Z- and OM -axes are e^2 a cos^3 ψ and − e ′^2 b sin^3 ψ respectively, thus an expression for tan φ in terms of
the parametric latitude ψ is
tan
sin cos
Z b e p a e
2 3 2 3 (6.1)
This equation, having functions of φ on both sides of the equation, clearly requires an iterative solution for
tan φ. An initial value ψ 0 is obtained from tan θ ′ =
Z
p
, θ ′ being a good approximation of the geocentric
latitude θ , and since tan ψ= tan θ≈ tan ′
a b
a b
θ , then
tan ψ 0 ≈
a Z b p
Bowring (1976) shows that for all Earth-bound points ( − 5 000, m ≤ h ≤10 000, m ) the maximum error in φ ,
induced by using only a single iteration, is 0 .000 000 030 ′′. Thus for all practical purposes, the evaluation of
tan φ by (6.1), with a first approximation of ψ from (6.2), can be regarded as exact.
BORKOWSKI’S METHOD
This is an indirect method using Newton’s iterative technique to solve for the parametric latitude ψ (Borkowski
1989). An outline of the development of the necessary equations is given below.
From the ellipsoid relationships, the rectangular components of h in the meridian plane are given by
h sin φ = Z − b sin ψ (7.1)
h cos φ = p − a cos ψ (7.2)
Eliminating h by division and substituting for tan φ gives
a b
Z b p a
tan
sin cos
which can be rearranged as
2 ap sin ψ− 2 bZ cos ψ− 2 sin ψ cos ψ( a^2 − b^2 ) = 0 (7.4)
Replacing ap and bZ with q cos Ω and q sin Ω respectively, gives
2 q ( cos Ω sin ψ− sin Ωcos ψ) − 2 sin ψ cos ψ( a^2 − b^2 )= 0 (7.5)
X
a
X
a
m a
E =
2 ,^
Y
a
Y
a
m a
E =
2 ,^
Z
b
Z
b
m b
E =
which, when squared and substituted into (8.1), and simplified, gives a function of m as
f m
p
a m a
Z
b m b
F +
H
I
K
F +
H
I
K
2 2
2
2 2 2 1 0 (8.6)
and m can be obtained by Newton’s iterative technique
m m f m n n (^) f m
n n
+ 1 =^ −^ ′
where the function and its derivative are
f m p
a
m a
Z
b
m b
n n n
F +
H
I
K
F +
H
I
K
2 2
2
2 2 2 1 (8.8)
F +
H
I
K
F +
H
I
K
R
S
T
U
V
W
f m
p
a a m a
Z
b b m b
n n n
2 3
2 3 (8.9)
An initial approximation m (^) 0 can be obtained from (Lin & Wang 1995, p. 301)
m
ab a Z b p a b a Z b p a Z b p
0
2 2 2 2 2 2 2 2 2 2 4 2 4 2
(^32)
d i d i
d i
After calculating m, pE and ZE are obtained from (8.5) as
p p m a
E =
2
Z
Z
m b
E =
2
and the latitude φ and height h obtained from
tan φ =
a Z b p
E E
2 2 (8.13)
h = ± b p − p E g 2 + b Z − Z E g^2 (8.14)
Note: h is negative if ( p + Z ) is less than ( p (^) E + ZE )
Thus, having X,Y,Z (hence p ) for a point related to an ellipsoid, m is obtained by iteration after first calculating a starting value m from (8.10) then using (8.7), with (8.8) and (8.9) calculated using successively improved
values of m. Having calculated m the coordinates of
0
PE on the ellipsoid are obtained and φ follows from (8.13).
The attraction of this method, compared to the other iterative techniques ( Simple Iteration, Bowring’s and
Borkowski’s ) is that no trigonometric functions are used in the calculation of tan φ or h.
COMPUTER TESTING OF THE TRANSFORMATION METHODS
The preference for one transformation method over another depends on the answer to three questions relevant to computer evaluation – (i) is it easy to code as a computer algorithm? (ii) is it accurate? and (iii) is it fast? In answer to the first two questions; all the methods are simple to code as computer algorithms, requiring similar numbers of lines of source code and each yielding accurate results, providing that Paul’s and Ozone’s direct methods are not used when Z = 0 or very small, in which cases they either fail or give unreliable results. The answer to the last question, often the determining factor in the adoption of a particular method, depends on two factors, (a) the number square roots, cube roots and trigonometric functions requiring evaluation and (b) the
number of iterations required for the indirect methods. With regard to these factors, Lin & Wang’s method
requires only a single evaluation of a square root to compute tan φ and, if points are restricted to an Earth-bound
region ( , , the analysis shows that Simple Iteration is the only indirect method actually
requiring multiple iterations to reach an acceptable accuracy for
− 5 000 m ≤ h ≤10 000 , m )
φ , the other indirect methods ( Bowring’s,
Borkowski’s and Lin & Wang’s ) giving, for all practical purposes, “exact” solutions for tan φ without iteration.
In order to evaluate accuracy and relative speeds, each method was coded as a computer algorithm (using Borland’s C++ version 5.0A compiler on a PC with a 133MHz Pentium processor) and tested over a grid of points at 0.1° intervals between latitudes 5°S and 50°S and longitudes 110°E and 160°E on the GRS80 ellipsoid, each assigned an ellipsoidal height h (225,450 points in total, covering the Australian continent, Tasmania and offshore islands). For each algorithm, the following procedure was adopted
= 10 000, m
- compute X,Y,Z coordinates for each point ( φ , λ, h ) using (1.1), (1.2) and (1.3) [these Cartesian and
geographic values regarded as “exact”] then, using the X,Y,Z coordinates
- compute the longitude λ (^) C using (1.7)
- compute the latitude φ (^) C of each point using the relevant equations for each method
- compute the height hC using (1.9) for each method, excepting Lin and Wang’s , where (8.14) was used
- compare φ (^) C and hC with the exact values to determine the magnitude of the “geographical errors”
- compute X,Y,Z again, using (1.1), (1.2) and (1.3) with computed values φ (^) C , λ C , hC and compare with exact
values to determine the magnitude of the “Cartesian errors”.
Note: Simple Iteration was the only indirect method where iteration was used, the others, Bowring’s, Borkowski’s and Lin & Wang’s, were used without iteration.
For each algorithm, the time required to perform the Cartesian-to-geographic transformations for all points was recorded (in processor “clock-ticks”) and converted to relative speed units, with Bowring’s method given a value of 50. Results of the comparative tests are given in Table 1.
Conversion Maximum Errors (magnitudes) Relative
Method^ φ^ (sec)^ h m (^ )^ X^ ( m^ )^ Y^ (^ m )^ Z^ ( m^ )^ Speed
Lin & Wang 6.87e-11 2.53e-09 9.31e-10 9.31e-10 1.86e-09 42 Bowring 2.88e-08 1.01e-06 9.31e-10 9.31e-10 1.32e-06 50 Paul 3.90e-07 1.14e-06 9.31e-10 9.31e-10 1.21e-05 52 Ozone 6.87e-11 2.42e-09 9.31e-10 9.31e-10 2.33e-09 56 Simple Iteration 1.35e-05 6.12e-05 9.31e-10 9.31e-10 4.19e-04 62 Borkowski 2.88e-08 1.01e-06 9.31e-10 9.31e-10 1.32e-06 63
Table 1 Cartesian-to-geographic conversion algorithms ranked from fastest ( Lin & Wang ) to slowest ( Borkowski ).
CONCLUSION AND RECOMMENDATION
This paper has presented the mathematical background to six Cartesian-to-geographic transformation techniques, two direct and four indirect. In theory, the indirect methods require iteration for an acceptable
solution of φ , but in practice, by restricting the use of the techniques to Earth-bound regions
, only the Simple Iteration technique actually requires iteration, the other indirect
methods giving, for all practical purposes, “exact” solutions for tan
( − 5 000, m ≤ h ≤10 000, m )
φ without iteration. All the methods were
coded as computer algorithms, tested (on a region covering Australia) and found to give acceptable results for
the computation of φ and h. Lin & Wang’s method is appreciably faster than the other five methods and is
recommended for use in Australia.
