The paper explores an ellipsoid-to-sphere latitude transformation with a particular combination of desirable properties when used for the internal coordinate encoding in global or wide-area 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 radius-vector of a point on the ellipsoid surface. The combination of desirable properties includes a simple, closed-form computation in both mapping directions; a straightforward derivation of the indicatrix; and a fast and convenient calculation specifically well suited for implementations with a system-wide 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 ellipsoid-to-sphere conformal mapping.

The paper includes the performance and numerical comparison of this type of mapping with its canonical alternative.

The motivation for exploring various ellipsoid-to-sphere 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 "pre-computer" 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 real-life 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 real-life activity.

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 two-dimensional 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.

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 Cauchy-Rieman 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 field-work 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 wide-area 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 radius-vector of a point on the ellipsoid surface) but, instead, a "near-conformal" ellipsoid-to-sphere latitude transformation, what kind of benefits would we obtain, and at what cost, when compared with its canonical alternative: use of the rigorous conformal ellipsoid-to-sphere 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

Computing the angle between the equatorial plane and the radius-vector of a point on the ellipsoid surface (phi_c, geocentric latitude), given the same angle of the ellipsoid normal (phi, ellipsoidal latitude) is simple:

Where a and b are semi-major and semi-minor ellipsoidal axes, and esq is the square of the ellipsoid eccentricity:

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 semi-axes (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:

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 implementation-dependent "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]).

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 "Core-2" 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:

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 near-conformal mapping would in such case be as follows:

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:

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:

Let u be the reciprocal of t:

The curvature of the prime vertical (p) in the given point will then be:

And r, curvature of the meridian:

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:

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.

Use of geocentric latitude as a near-conformal "auxiliary" sphere for topological and metric computations, and as a basis for the high-volume coordinate data storage formats results in computations that are particularly well suited for use on digital computers. Mapping transformations are significantly faster than rigorous ellipsoid-to-sphere 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 problem-local vicinity on the ellipsoid and the near-conformal sphere - is equally straightforward.

1. Bomford, G.: Geodesy London, Oxford University Press, 1975 ed.
2. Jordan, Eggert and Kneissl: Handbuch der Vermessungskunde (Band IV: Mathematische Geodaesie), 10th ed. Metzlersche Verlag., Stuttgart, 1959
3. John P. Snyder: Map Projections Used by the U.S. Geological Survey Geological Survey Bulletin 1532, Second Edition, 1982
4. Qihe Yang, John P. Snyder and Waldo R. Tobler: Map Projection Transformation Taylor & Francis, 2000