- Generated by
1.9.4
|
GeographicLib 2.1.2
|
In addition to solving the geodesic problem for the triaxial ellipsoid, Jacobi (1839) briefly mentions the problem of the conformal projection of ellipsoid. He covers this in greater detail in Vorlesungen über Dynamik, §28, which is now available in an English translation: Lectures on Dynamics (errata).
It is convenient to rotate Jacobi's ellipsoidal coordinate system so that \((X,Y,Z)\) are respectively the middle, large, and small axes of the ellipsoid. The longitude \(\omega\) is shifted so that \((\beta=0,\omega=0)\) is the point \((b,0,0)\). The coordinates are thus defined by
\[ \begin{align} X &= b \cos\beta \cos\omega, \\ Y &= a \sin\omega \frac{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}} {\sqrt{a^2 - c^2}}, \\ Z &= c \sin\beta \frac{\sqrt{a^2\cos^2\omega + b^2\sin^2\omega - c^2}} {\sqrt{a^2 - c^2}}. \end{align} \]
After this change, the large principal ellipse is the equator, \(\beta=0\), while the small principal ellipse is the prime meridian, \(\omega=0\). The four umbilic points, \(\left|\omega\right| = \left|\beta\right| = \frac12\pi\), lie on middle principal ellipse in the plane \(X=0\).
Jacobi gives the following conformal mapping of the triaxial ellipsoid onto a plane
\[ \begin{align} x &= \frac{\sqrt{a^2-c^2}}b \int \frac {\sqrt{a^2 \cos^2\omega + b^2 \sin^2\omega}} {\sqrt{a^2 \cos^2\omega + b^2 \sin^2\omega - c^2}}\, d\omega, \\ y &= \frac{\sqrt{a^2-c^2}}b \int \frac {\sqrt{b^2 \sin^2\beta + c^2 \cos^2\beta}} {\sqrt{a^2 - b^2 \sin^2\beta - c^2 \cos^2\beta}}\, d\beta. \end{align} \]
The scale of the projection is
\[ k = \frac{\sqrt{a^2-c^2}} {b\sqrt{a^2 \cos^2\omega + b^2 (\sin^2\omega-\sin^2\beta) - c^2 \cos^2\beta}}. \]
I have scaled the Jacobi's projection by a constant factor,
\[ \frac{\sqrt{a^2-c^2}}{2b}, \]
so that it reduces to the familiar formulas in the case of an oblate ellipsoid.
The projection may be expressed in terms of elliptic integrals,
\[ \begin{align} x&=(1+e_a^2)\,\Pi(\omega',-e_a^2, \cos\nu),\\ y&=(1-e_c^2)\,\Pi(\beta' , e_c^2, \sin\nu),\\ k&=\frac{\sqrt{e_a^2+e_c^2}}{b\sqrt{e_a^2\cos^2\omega+e_c^2\cos^2\beta}}, \end{align} \]
where
\[ \begin{align} e_a &= \frac{\sqrt{a^2-b^2}}b, \qquad e_c = \frac{\sqrt{b^2-c^2}}b, \\ \tan\omega' &= \frac ba \tan\omega = \frac 1{\sqrt{1+e_a^2}} \tan\omega, \\ \tan\beta' &= \frac bc \tan\beta = \frac 1{\sqrt{1-e_c^2}} \tan\beta, \\ \tan\nu &= \frac ac \frac{\sqrt{b^2-c^2}}{\sqrt{a^2-b^2}} =\frac{e_c}{e_a}\frac{\sqrt{1+e_a^2}}{\sqrt{1-e_c^2}}, \end{align} \]
\(\nu\) is the geographic latitude of the umbilic point at \(\beta = \omega = \frac12\pi\) (the angle a normal at the umbilic point makes with the equatorial plane), and \(\Pi(\phi,\alpha^2,k)\) is the elliptic integral of the third kind, https://dlmf.nist.gov/19.2.E7. This allows the projection to be numerically computed using EllipticFunction::Pi(real phi) const.
Nyrtsov, et al.,
also expressed the projection in terms of elliptic integrals. However, their expressions don't allow the limits of ellipsoids of revolution to be readily recovered. The relations https://dlmf.nist.gov/19.7.E5 can be used to put their results in the form given here.
\(x\) (resp. \(y\)) depends on \(\omega\) (resp. \(\beta\)) alone, so that latitude-longitude grid maps to straight lines in the projection. In this sense, the Jacobi projection is the natural generalization of the Mercator projection for the triaxial ellipsoid. (See below for the limit \(a\rightarrow b\), which makes this connection explicit.)
In the general case (all the axes are different), the scale diverges only at the umbilic points. The behavior of these singularities is illustrated by the complex function
\[ f(z;e) = \cosh^{-1}(z/e) - \log(2/e). \]
For \(e > 0\), this function has two square root singularities at \(\pm e\), corresponding to the two northern umbilic points. Plotting contours of its real (resp. imaginary) part gives the behavior of the lines of constant latitude (resp. longitude) near the north pole in Fig. 1. If we pass to the limit \(e\rightarrow 0\), then \( f(z;e)\rightarrow\log z\), and the two branch points merge yielding a stronger (logarithmic) singularity at \(z = 0\), concentric circles of latitude, and radial lines of longitude.
Again in the general case, the extents of \(x\) and \(y\) are finite,
\[ \begin{align} x\bigl(\tfrac12\pi\bigr) &=(1+e_a^2)\,\Pi(-e_a^2, \cos\nu),\\ y\bigl(\tfrac12\pi\bigr) &=(1-e_c^2)\,\Pi(e_c^2, \sin\nu), \end{align} \]
where \(\Pi(\alpha^2,k)\) is the complete elliptic integral of the third kind, https://dlmf.nist.gov/19.2.E8.
In particular, if we substitute values appropriate for the earth,
\[ \begin{align} a&=(6378137+35)\,\mathrm m,\\ b&=(6378137-35)\,\mathrm m,\\ c&=6356752\,\mathrm m,\\ \end{align} \]
we have
\[ \begin{align} x\bigl({\textstyle\frac12}\pi\bigr) &= 1.5720928 = \hphantom{0}90.07428^\circ,\\ y\bigl({\textstyle\frac12}\pi\bigr) &= 4.2465810 = 243.31117^\circ.\\ \end{align} \]
The projection may be inverted (to give \(\omega\) in terms of \(x\) and \(\beta\) in terms of \(y\)) by using Newton's method to find the root of, for example, \(x(\omega) - x_0 = 0\). The derivative of the elliptic integral is, of course, just given by its defining relation.
If rhumb lines are defined as curves with a constant bearing relative to the ellipsoid coordinates, then these are straight lines in the Jacobi projection. A rhumb line which passes over an umbilic point immediately retraces its path. A rhumb line which crosses the line joining the two northerly umbilic points starts traveling south with a reversed heading (e.g., a NE heading becomes a SW heading). This behavior is preserved in the limit \(a\rightarrow b\) (although the longitude becomes indeterminate in this limit).
Oblate ellipsoid, \(a\rightarrow b\). The coordinate system is
\[ \begin{align} X &= b \cos\beta \cos\omega, \\ Y &= b \cos\beta \sin\omega, \\ Z &= c \sin\beta. \end{align} \]
Thus \(\beta\) (resp. \(\beta'\)) is the parametric (resp. geographic) latitude and \(\omega=\omega'\) is the longitude; the quantity \(e_c\) is the eccentricity of the ellipsoid. Using https://dlmf.nist.gov/19.6.E12 and https://dlmf.nist.gov/19.2.E19 the projection reduces to the normal Mercator projection for an oblate ellipsoid,
\[ \begin{align} x &= \omega,\\ y &= \sinh^{-1}\tan\beta' - e_c \tanh^{-1}(e_c\sin\beta'),\\ k &= \frac1{b\cos\beta}. \end{align} \]
Prolate ellipsoid, \(c\rightarrow b\). The coordinate system is
\[ \begin{align} X &= b \cos\omega \cos\beta, \\ Y &= a \sin\omega, \\ Z &= b \cos\omega \sin\beta. \end{align} \]
Thus \(\omega\) (resp. \(\omega'\)) now plays the role of the parametric (resp. geographic) latitude and while \(\beta=\beta'\) is the longitude. Using https://dlmf.nist.gov/19.6.E12 https://dlmf.nist.gov/19.2.E19 and https://dlmf.nist.gov/19.2.E18 the projection reduces to similar expressions with the roles of \(\beta\) and \(\omega\) switched,
\[ \begin{align} x &= \sinh^{-1}\tan\omega' + e_a \tan^{-1}(e_a\sin\omega'),\\ y &= \beta,\\ k &= \frac1{b\cos\omega}. \end{align} \]
Sphere, \(a\rightarrow b\) and \(c\rightarrow b\). This is a non-uniform limit depending on the parameter \(\nu\),
\[ \begin{align} X &= b \cos\omega \cos\beta, \\ Y &= b \sin\omega \sqrt{1 - \sin^2\nu\sin^2\beta}, \\ Z &= b \sin\beta \sqrt{1 - \cos^2\nu\sin^2\omega}. \end{align} \]
Using https://dlmf.nist.gov/19.6.E13 the projection can be put in terms of the elliptic integral of the first kind, https://dlmf.nist.gov/19.2.E4
\[ \begin{align} x &= F(\omega, \cos\nu), \\ y &= F(\beta, \sin\nu), \\ k &= \frac1{b \sqrt{\cos^2\nu\cos^2\omega + \sin^2\nu\cos^2\beta}}. \end{align} \]
Obtaining the limit of a sphere via an oblate (resp. prolate) ellipsoid corresponds to setting \(\nu = \frac12\pi\) (resp. \(\nu =0\)). In these limits, the elliptic integral reduces to elementary functions https://dlmf.nist.gov/19.6.E7 and https://dlmf.nist.gov/19.6.E8
\[ \begin{align} F(\phi, 0) &= \phi, \\ F(\phi, 1) &= \mathop{\mathrm{gd}}\nolimits^{-1}\phi = \sinh^{-1}\tan\phi. \end{align} \]
The spherical limit gives the projection found by É. Guyou in
who apparently derived it without realizing that it is just a special case of the projection Jacobi had given some 40 years earlier. Guyou's name is usually associated with the particular choice, \(\nu=\frac14\pi\), in which case the hemisphere \(\left|\omega\right|\le\frac12\pi\) is mapped into a square. However, by varying \(\nu\in[0,\frac12\pi]\) the hemisphere can be mapped into a rectangle with any aspect ratio, \(K(\cos\nu) : K(\sin\nu)\), where \(K(k)\) is the complete elliptic integral of the first find, https://dlmf.nist.gov/19.2.E8.
An essential tool in deriving conformal projections of an ellipsoid of revolution is the conformal mapping of the ellipsoid onto a sphere. This allows conformal projections of the sphere to be generalized to the case of an ellipsoid of revolution. This conformal mapping is obtained by using the ellipsoidal Mercator projection to map the ellipsoid to the plane and then using the spherical Mercator projection to map the plane onto the sphere.
A similar construction is possible for a triaxial ellipsoid. Map each octant of the ellipsoid onto a rectangle using the Jacobi projection. The aspect ratio of this rectangle is
\[ (1+e_a^2)\,\Pi(-e_a^2, \cos\nu) : (1-e_c^2)\,\Pi( e_c^2, \sin\nu). \]
Find the value of \(\nu'\) such that this ratio equals
\[ K(\cos\nu') : K(\sin\nu'). \]
Map the rectangle onto the equivalent octant of the sphere using Guyou's projection with parameter \(\nu = \nu'\). This reduces to the standard construction in the limit of an ellipsoid of revolution.
The JacobiConformal class provides an implementation of the Jacobi conformal projection is given here. NOTE: This is just sample code. It is not part of GeographicLib itself.
1.9.4