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 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:
phi_c = arctan((1 - esq) * tan(phi))
Where a and b are semi-major and semi-minor 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 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:
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 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:
Near-conformal mapping: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:
ellipsoid to sphere: 4.578The 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.437The 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.
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.