Ellipsoid RadiusVector as a NearConformal Ellipsoid to Sphere Mapping
Hrvoje Lukatela,
Geodyssey Limited
http://www.geodyssey.com/
Paper proposal, ACM International Symposium on Advances in Geographic Information Systems (ACM GIS 2007)
Abstract
The paper explores an ellipsoidtosphere latitude transformation
with a particular combination of desirable properties when used for
the internal coordinate encoding in global or widearea GIS systems.
In this transformation, the spherical latitude is obtained by the same
production as that used to compute the angle between the equatorial
plane and the radiusvector of a point on the ellipsoid surface.
The combination of desirable properties includes a simple, closedform
computation in both mapping directions; a straightforward derivation
of the indicatrix; and a fast and convenient calculation specifically
well suited for implementations with a systemwide preference for
solid vector algebra over the spherical or spheroidal trigonometry.
However, the most important property of this transformation is its
numerical closeness to the rigorous ellipsoidtosphere conformal
mapping.
The paper includes the performance and numerical comparison of this
type of mapping with its canonical alternative.
Introduction
The motivation for exploring various ellipsoidtosphere mapping
transformations is as old as is the use of an ellipsoid of rotation
as the computational geometry reference surface: spherical productions
are both simpler and faster. These two characteristics have become not
less, but more important as the computations are carried out on a
digital computer: the simplicity benefits software design, testing
and verification, and the speed becomes critical as the volume of data
encompassed in even some trivial systems outstrips the volume of
positional information handled by the national control grid projects
of the "precomputer" time. However, the numerical differences between
the ellipsoid and the spherical productions must be well defined and
understood, as the simplification must not exceed the limits imposed
by the reallife activity that is served by the computational system.
We are therefore in search of such simplifications that will result in
errors which are safely below the level of an order of magnitude more
precise than the geodetic measurements of the reallife activity.
Auxiliary latitudes
Snyder [Snyder p. 16] provides a comprehensive list, definitions
and the computational details of transformations that have been used
for the purpose outlined above. As usual in the geodetic and cartographic
practice, Snyder calls them "Auxiliary Latitudes". They are:
 Conformal
 Authalic
 Rectifying
 Geocentric
 Reduced
Conformal refers to a generic mathematical property: it requires
that an infinitesimally small twodimensional object, mapped from
one surface to another (in our case from an ellipsoid to the sphere)
may change the scale but must preserve the shape. This is the same
as stating that the mapping scales in two principal directions
(the meridian and the prime vertical) must be the same. [footnote:
[Jordan] provides a comprehensive treatment of the geodetic fundamentals
of the subject of this paper; [Qihe] of its cartographic equivalent)].
Authalic mapping preserves the area of small shapes, and
rectifying preserves the scale along the meridian. Geocentric
latitude forms the basis of the mapping discussed further in this text.
Reduced (or, as also termed in [Snyder], "parametric")
latitude is numerically quite close to the true ellipsoidal latitude,
and is defined, similarly to the geocentric, by a simple
geometrical construct.
Conformal mapping
The first in Snyder's list of Auxiliary Latitudes is also the most
important one: not only in the context of the historical development
of Geodesy and Cartography, but also because it provides a fundamental
building block of a large number of numerical recipes in the current
use. The derivation [see Qihe] is based on the CauchyRieman conditions
of general conformal mapping from one surface to another, and its
standard implementation has been brought into the geodetic practice
by C. F. Gauss.
Important practical reasons led to the wide adoption of the conformal
ellipsoid mapping (in both variants, planar and spherical): classical
geodetic fieldwork was performed overwhelmingly by the measurement of
angles. If the sides of the triangles (typically in the creation of
control networks) or the legs of the traverses (in the subsequent
engineering surveying) are relatively short (relatively, that is, to
the radius of the planet), then the angles can be used in computations
on the conformal sphere as they are measured on the ellipsoid. For
the computations in the plane (for instance, in a widearea surveying
and mapping planar coordinate system based on the UTM projection) the
triangles may require nothing but a trivial, location independent
adjustment for the spherical excess.
Snyder's discussion of Auxiliary Latitudes concludes [Snyder, p.22]
with a numerical table of "corrections" (actually, a table of
differences) between the true ellipsoidal latitude and each of the
auxiliary ones.
The table is particularly interesting in one detail: a large degree
of numerical closeness between the conformal latitude and the
geocentric one. This closeness prompts the question: if we
consider the geocentric latitude not to be the geometric entity that
it actually is (a direction of the radiusvector of a point on the
ellipsoid surface) but, instead, a "nearconformal"
ellipsoidtosphere latitude transformation, what kind of benefits
would we obtain, and at what cost, when compared with its canonical
alternative: use of the rigorous conformal ellipsoidtosphere mapping.
When modified, to address primarily the difference between the
conformal and geocentric latitudes, computed for WGS84 ellipsoid
(Snyder's examples are for Clarke 1866 ellipsoid), and including
the second differences, the table looks as follows:
Ellipsoid Conformal Geocentric d dd

n0:00:00.000 n0:00:00.000 n0:00:00.000 0.000 0.000
n5:00:00.000 n4:58:00.107 n4:58:00.106 0.001 0.001
n10:00:00.000 n9:56:03.827 n9:56:03.819 0.008 0.007
n15:00:00.000 n14:54:14.667 n14:54:14.642 0.026 0.018
n20:00:00.000 n19:52:35.925 n19:52:35.868 0.058 0.032
n25:00:00.000 n24:51:10.590 n24:51:10.485 0.105 0.047
n30:00:00.000 n29:50:01.255 n29:50:01.089 0.166 0.061
n35:00:00.000 n34:49:10.037 n34:49:09.799 0.238 0.072
n40:00:00.000 n39:48:38.512 n39:48:38.198 0.314 0.076
n45:00:00.000 n44:48:27.663 n44:48:27.276 0.386 0.072
n50:00:00.000 n49:48:37.849 n49:48:37.402 0.447 0.061
n55:00:00.000 n54:49:08.792 n54:49:08.304 0.489 0.041
n60:00:00.000 n59:49:59.578 n59:49:59.074 0.504 0.015
n65:00:00.000 n64:51:08.683 n64:51:08.194 0.489 0.015
n70:00:00.000 n69:52:34.018 n69:52:33.576 0.442 0.047
n75:00:00.000 n74:54:12.990 n74:54:12.627 0.363 0.078
n80:00:00.000 n79:56:02.582 n79:56:02.324 0.258 0.105
n85:00:00.000 n84:57:59.445 n84:57:59.310 0.134 0.124
n90:00:00.000 n90:00:00.000 n90:00:00.000 0.000 0.134
The transformation
Computing the angle between the equatorial plane and the radiusvector
of a point on the ellipsoid surface (phi_c, geocentric latitude),
given the same angle of the ellipsoid normal (phi, ellipsoidal
latitude) is simple:
phi_c = arctan((1  esq) * tan(phi))
Where a and b are semimajor and semiminor ellipsoidal
axes, and esq is the square of the ellipsoid eccentricity:
esq = ((a * a)  (b * b)) / (a * a)
It is generally prudent to avoid the use of angular forms of ellipsoid
coordinates in the computational systems that use modern digital
computers, and instead use the vector notation (i.e., the direction
cosines). Indeed, [Bomford] states as early as 1975 "...such form
is symmetrical and much preferred for computations...". If the system
records and operates on the sine (s) and cosine (c) of
latitude, and if instead of eccentricity, the system uses the squares
of the major and minor ellipsoidal semiaxes (asq and bsq,
respectively); and if the system will use the results of the mapping in
the vector form (u and v, respectively, as the sine and
cosine of the spherical latitude), the symmetrical form, (perhaps as
suggested by Bomford?), will have the following form:
u = bsq * s
v = asq * c
The above discussion considers only the meridian ellipse, and ignores
the fact that in practical implementations the cosine of the latitude
will not be known or recorded as an independent value, but will
instead be a combination of two direction cosines of the longitude of
some point on the ellipsoid. This, however, means only that another
square root operation will be required to "fold" and "unfold" the
longitude into and from the meridian ellipse: spherical longitude in
all such mappings is equal to the ellipsoidal longitude.
To implement such latitude mapping and its inverse in computer code
becomes quite a bit simpler than doing the same for the rigorous
conformal transformation (cf. [Quihe], p. 79, production 3.5.1).
Indeed, most implementations will, instead of the transcendental and
natural logarithm based production mentioned, use an expansion based on
the series of terms of the ascending powers of ellipsoid eccentricity
(Quihe, p81, production 3.5.0, which includes the terms with
e**8 ). Unlike rigorous conformal mapping, the
transformation productions are closed in both directions: this means
that there will be no implementationdependent "drift" if one
point is, for some reason, repeatedly transferred between the two
surfaces. In contrast, in the case of the rigorous conformal sphere
inverse mapping, there is no option but to use either the iteration
(as in [Quihe], p.86) or the expansion [Snyder, p.18]).
Execution speed comparisons
The following values have been obtained by creating an array of ten
million points, randomly positioned on an ellipsoid and then each
subjected to the transformation in a tight loop. The coding was done
in C language, with gcc compiler, and the execution time was measured
on a PC with Intel Pentium "Core2" CPU, running at 1.8 GHz clock rate.
The times are in seconds.
The first set of values assumes that both given and computed
coordinates are in the vector form:
Nearconformal mapping:
ellipsoid to sphere: 1.297
sphere to ellipsoid: 1.234
Rigorous conformal mapping:
ellipsoid to sphere: 6.359
sphere to ellipsoid: 30.329
Conformal sphere to ellipsoid mapping in the above was implemented using
iteration, with the closeness criterion of about 1/10th of a millimeter
on the Earth's surface.
It is, however, reasonable to assume that a GIS system will be designed
so that ellipsoidal values will be imported and exported from the system
in their conventional, angular form, and the spherical points used for
internal computations and data storage will be in the vector (direction
cosine) form. The above values for the nearconformal mapping would in
such case be as follows:
ellipsoid to sphere: 4.578
sphere to ellipsoid: 3.640
The simple transformation of coordinates from angular into vector
form and its inverse requires evaluation of sines and cosines. Such
transformation will benefit greatly if the evaluation of the sine and
cosine of an angle in a single hardware instruction (a feature available
on many floating point instruction sets) is used by the code. For the
same set of points as above, the values are:
Angle to vector: 3.437
Vector to angle: 3.422
The indicatrix
The computation of the two principal radii of curvature on ellipsoid
is simplified by the introduction of t, the free term of the
tangential plane equation. If the ellipsoid point coordinates (i.e.,
the direction of the normal to the ellipsoid) are given in their
normalized vector form (i,j and k), t
can be computed in a "symmetrical form" similar to the ones used
above:
t = asq * i * i + asq * j * j + bsq * k * k
Let u be the reciprocal of t:
u = 1 / t
The curvature of the prime vertical (p) in the given point will
then be:
p = asq * u
And r, curvature of the meridian:
r = p * bsq * u * u
General curvature q on the ellipsoid (i.e., one in the azimuth
defined by two direction cosines in the tangential plane, f
and g can then be computed as follows:
q = (f * f) / p + (g * g) / r
The above productions can be used to determine the elements of the
indicatrix for any point on the ellipsoid: the radius of curvature of
the sphere is 1, and for the purpose of the computation of indicatrix
axes, such sphere can be safely considered to be conformal.
Conclusion
Use of geocentric latitude as a nearconformal "auxiliary" sphere
for topological and metric computations, and as a basis for the
highvolume coordinate data storage formats results in computations
that are particularly well suited for use on digital computers.
Mapping transformations are significantly faster than rigorous
ellipsoidtosphere conformal mapping. If the ellipsoidal
computations are performed on the coordinates in their vector form,
the parameters of the indicatrix  the basis for the scaling of metric
results between the problemlocal vicinity on the ellipsoid and the
nearconformal sphere  is equally straightforward.
References

Bomford, G.: Geodesy
London, Oxford University Press, 1975 ed.

Jordan, Eggert and Kneissl: Handbuch der Vermessungskunde
(Band IV: Mathematische Geodaesie), 10th ed.
Metzlersche Verlag., Stuttgart, 1959

John P. Snyder: Map Projections Used by the U.S. Geological Survey
Geological Survey Bulletin 1532, Second Edition, 1982

Qihe Yang, John P. Snyder and Waldo R. Tobler:
Map Projection Transformation Taylor & Francis, 2000
...
