[Proj] Spheroidal gnomonic projection

Tue Jun 1 11:18:50 PDT 2010

Cole did the two-point azimuthal equidistant on the Helmert ellipsoid back in the early 1900's for the Survey of Egypt.  I do not have a copy, though.  It's pretty obscure, but the intent was to provide a graphic for determining the direction to face for daily prayers.

Clifford J. Mugnier, C.P., C.M.S.
Chief of Geodesy,
Center for GeoInformatics
Department of Civil Engineering 
Patrick F. Taylor Hall 3223A
LOUISIANA STATE UNIVERSITY 
Baton Rouge, LA  70803
Voice and Facsimile:  (225) 578-8536 [Academic] 
Voice and Facsimile:  (225) 578-4578 [Research] 
Cell: (225) 238-8975 [Academic & Research]
Honorary Life Member of the 
Louisiana Society of Professional Surveyors 
Fellow Emeritus of the ASPRS 
Member of the Americas Petroleum Survey Group

________________________________

From: proj-bounces at lists.maptools.org on behalf of Charles Karney
Sent: Tue 01-Jun-10 12:42
To: proj at lists.maptools.org
Cc: deakin at rmit.edu.au; william at karney.com; kevin at karney.com; knud.poder at mail.dk; Craig.M.Rollins at nga.mil
Subject: [Proj] Spheroidal gnomonic projection

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
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20100601/90d685be/attachment.html>