[Proj] Spheroidal gnomonic projection

[Charles Karney]

The gnomonic map projection is a central projection of the sphere on a
tangent plane.  The key property of the projection is that all geodesics
map to straight lines in the projection.

This property cannot be preserved for an ellipsoid.  However, we can
obtain a projection where the projected geodesics are approximately
straight close to the center of the projection.  It is obtained as the
limit of a 2-point azimuthal projection as the two points approach one
another.  The method generalizes to any 2-dimensional surface; but I
only have the necessary formulas worked out for the ellipsoid.

This projection has the following properties:

 (1) azimuthal, all azimuths from the center are correct;
 (2) hence, all lines through the center are geodesics;
 (3) all other lines are *approximately* geodesics.

To quantify point 3, consider this projection with some specified center
on the WGS84 ellipsoid.  Take an arbitrary pair of points within r =
1000 km of the center.  Draw a straight line between these on the map;
then back project this line onto the ellipsoid.  How close is this line
to a geodesic?  I find that *at worst*,

  max error in initial/final azimuth = 1.0" = 0.39 * f*(r/a)^3
  max deviation from geodesic = 1.66 m = 0.13 * f*(r/a)^3*r

(a = major radius, f = flattening, and the formulas show the scaling of
the errors).  For comparison, for a gnomonic projection obtained by
projecting from the center of the ellipsoid (i.e., approximating the
geodesic by a great ellipse), the equivalent figures are about 100 times
worse:

  max error in initial/final azimuth = 108" = 1.0 * f*(r/a)
  max deviation from geodesic = 265 m = 0.5 * f*(r/a)*r

The method of constructing the projection entails defining two
quantities,

  * m, the reduced length.  Take two geodesics of length s which start
    from the center with azimuths differing by dalpha.  The endpoints
    are separated by m * dalpha.

  * M, the geodesic scale.  Take two geodesics of length s which are
    parallel close to the center and initially separated by dt.  The
    endpoints are separated by M * dt.

The reduced length was introduced by Gauss (1827) and Christoffel
(1868).  Helmert (1880), Sec. 6.5, gives an explicit formula for the
ellipsoid.  The geodesic scale was introduced (I think) by

  G. V. Bagratuni,
  Course in spheroidal geodesy,
  FTD-MT-64-390 (US Air Force, 1967), Sec. 17
  [Kurs sferoidicheskoi geodezii, (Moscow, 1962)].
  http://handle.dtic.mil/100.2/AD650520

(Bagratuni uses the symbol n.)  An explicit formula for the ellipsoid is
given in version 1.2 of GeographicLib by GeodesicLine::Scale.

To carry out the projection for an arbitrary point, solve the inverse
geodesic problem between the center and the point.  Let alpha be the
azimuth of the geodesic at the center and s be its length.  Compute

  rho = m/M

Then the projected point is

  x = rho * sin(alpha); y = rho * cos(alpha)

The radial scale is 1/M^2 and the transverse scale is 1/M.  For a
sphere, m = a * sin(s/a) and M = cos(s/a), thus

  rho = a * tan(s/a)

which gives the spherical gnomonic projection.

In order to reverse the projection, compute

  alpha = atan2(x, y); rho = sqrt(x^2 + y^2)

Use Newton's method to solve for s given rho; for this we need drho/ds
which is given by the radial scale = 1/M^2.  Solve the direct geodesic
problem from the center using length s and azimuth alpha to recover the
original point.

I just checked in a Gnomonic class for GeographicLib that implements the
forward and reverse projections.  It's available through svn from
SourceForge.

This projection minimizes the error in the straightness of geodesics
near the center.  If, instead, you want to minimize the error over some
region, then use a two-point azimuthal projection with two base points
judiciously chosen within the region.

Questions:

(1) Is this new?
(2) Is Bagratuni this first person to define M?
(3) Is GeographicLib the first place where the formula for M for an
    ellipsoid is given?
(4) Does anyone have a better copy of Bagratuni than the one given
    above.  (I'll even take a good copy in the Russian original.)

-- 
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662