<HTML dir=ltr><HEAD><TITLE>[Proj] Spheroidal gnomonic projection</TITLE>
<META content="text/html; charset=unicode" http-equiv=Content-Type>
<META name=GENERATOR content="MSHTML 8.00.7600.16535"></HEAD>
<BODY>
<DIV dir=ltr id=idOWAReplyText20813>
<DIV dir=ltr><FONT color=#000000 size=2 face="Times New Roman">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.</FONT></DIV>
<DIV dir=ltr><FONT color=#000000 size=2 face="Times New Roman"></FONT> </DIV></DIV>
<DIV dir=ltr id=idSignature49189>
<DIV><FONT color=#000000 size=2 face="Times New Roman"><SPAN style="FONT-SIZE: 10pt">Clifford J. Mugnier, C.P., C.M.S.<BR>Chief of Geodesy,</SPAN></FONT></DIV>
<DIV><FONT color=#000000 size=2 face="Times New Roman"><SPAN style="FONT-SIZE: 10pt"><SPAN style="FONT-VARIANT: small-caps; FONT-FAMILY: 'Times New Roman','serif'; COLOR: black; FONT-SIZE: 10pt; mso-ascii-theme-font: minor-latin; mso-fareast-font-family: 'Times New Roman'; mso-fareast-theme-font: minor-latin; mso-hansi-theme-font: minor-latin; mso-ansi-language: EN-US; mso-fareast-language: EN-US; mso-bidi-language: EN-US"><STRONG>Center for GeoInformatics</STRONG></SPAN></SPAN></FONT></DIV>
<DIV><FONT color=#000000 size=2 face="Times New Roman"><SPAN style="FONT-SIZE: 10pt"><SPAN style="FONT-VARIANT: small-caps; FONT-FAMILY: 'Times New Roman','serif'; COLOR: black; FONT-SIZE: 10pt; mso-ascii-theme-font: minor-latin; mso-fareast-font-family: 'Times New Roman'; mso-fareast-theme-font: minor-latin; mso-hansi-theme-font: minor-latin; mso-ansi-language: EN-US; mso-fareast-language: EN-US; mso-bidi-language: EN-US"></SPAN>Department of Civil Engineering <BR>Patrick F. Taylor Hall 3223A<BR><STRONG>LOUISIANA STATE UNIVERSITY</STRONG> <BR>Baton Rouge, LA 70803<BR>Voice and Facsimile: (225) 578-8536 [Academic] <BR>Voice and Facsimile: (225) 578-4578 [Research] </SPAN></FONT></DIV>
<DIV><FONT color=#000000 size=2 face="Times New Roman"><SPAN style="FONT-SIZE: 10pt"><SPAN style="FONT-SIZE: 10pt">Cell: (225) 238-8975 [Academic & Research]<BR></SPAN>Honorary Life Member of the <BR>Louisiana Society of Professional Surveyors <BR>Fellow Emeritus of the ASPRS <BR>Member of the Americas Petroleum Survey Group<BR></SPAN></FONT></DIV></DIV>
<DIV dir=ltr><BR>
<HR tabIndex=-1>
<FONT size=2 face=Tahoma><B>From:</B> proj-bounces@lists.maptools.org on behalf of Charles Karney<BR><B>Sent:</B> Tue 01-Jun-10 12:42<BR><B>To:</B> proj@lists.maptools.org<BR><B>Cc:</B> deakin@rmit.edu.au; william@karney.com; kevin@karney.com; knud.poder@mail.dk; Craig.M.Rollins@nga.mil<BR><B>Subject:</B> [Proj] Spheroidal gnomonic projection<BR></FONT><BR></DIV>
<DIV>
<P><FONT size=2>The gnomonic map projection is a central projection of the sphere on a<BR>tangent plane. The key property of the projection is that all geodesics<BR>map to straight lines in the projection.<BR><BR>This property cannot be preserved for an ellipsoid. However, we can<BR>obtain a projection where the projected geodesics are approximately<BR>straight close to the center of the projection. It is obtained as the<BR>limit of a 2-point azimuthal projection as the two points approach one<BR>another. The method generalizes to any 2-dimensional surface; but I<BR>only have the necessary formulas worked out for the ellipsoid.<BR><BR>This projection has the following properties:<BR><BR> (1) azimuthal, all azimuths from the center are correct;<BR> (2) hence, all lines through the center are geodesics;<BR> (3) all other lines are *approximately* geodesics.<BR><BR>To quantify point 3, consider this projection with some specified center<BR>on the WGS84 ellipsoid. Take an arbitrary pair of points within r =<BR>1000 km of the center. Draw a straight line between these on the map;<BR>then back project this line onto the ellipsoid. How close is this line<BR>to a geodesic? I find that *at worst*,<BR><BR> max error in initial/final azimuth = 1.0" = 0.39 * f*(r/a)^3<BR> max deviation from geodesic = 1.66 m = 0.13 * f*(r/a)^3*r<BR><BR>(a = major radius, f = flattening, and the formulas show the scaling of<BR>the errors). For comparison, for a gnomonic projection obtained by<BR>projecting from the center of the ellipsoid (i.e., approximating the<BR>geodesic by a great ellipse), the equivalent figures are about 100 times<BR>worse:<BR><BR> max error in initial/final azimuth = 108" = 1.0 * f*(r/a)<BR> max deviation from geodesic = 265 m = 0.5 * f*(r/a)*r<BR><BR>The method of constructing the projection entails defining two<BR>quantities,<BR><BR> * m, the reduced length. Take two geodesics of length s which start<BR> from the center with azimuths differing by dalpha. The endpoints<BR> are separated by m * dalpha.<BR><BR> * M, the geodesic scale. Take two geodesics of length s which are<BR> parallel close to the center and initially separated by dt. The<BR> endpoints are separated by M * dt.<BR><BR>The reduced length was introduced by Gauss (1827) and Christoffel<BR>(1868). Helmert (1880), Sec. 6.5, gives an explicit formula for the<BR>ellipsoid. The geodesic scale was introduced (I think) by<BR><BR> G. V. Bagratuni,<BR> Course in spheroidal geodesy,<BR> FTD-MT-64-390 (US Air Force, 1967), Sec. 17<BR> [Kurs sferoidicheskoi geodezii, (Moscow, 1962)].<BR> <A href="http://handle.dtic.mil/100.2/AD650520">http://handle.dtic.mil/100.2/AD650520</A><BR><BR>(Bagratuni uses the symbol n.) An explicit formula for the ellipsoid is<BR>given in version 1.2 of GeographicLib by GeodesicLine::Scale.<BR><BR>To carry out the projection for an arbitrary point, solve the inverse<BR>geodesic problem between the center and the point. Let alpha be the<BR>azimuth of the geodesic at the center and s be its length. Compute<BR><BR> rho = m/M<BR><BR>Then the projected point is<BR><BR> x = rho * sin(alpha); y = rho * cos(alpha)<BR><BR>The radial scale is 1/M^2 and the transverse scale is 1/M. For a<BR>sphere, m = a * sin(s/a) and M = cos(s/a), thus<BR><BR> rho = a * tan(s/a)<BR><BR>which gives the spherical gnomonic projection.<BR><BR>In order to reverse the projection, compute<BR><BR> alpha = atan2(x, y); rho = sqrt(x^2 + y^2)<BR><BR>Use Newton's method to solve for s given rho; for this we need drho/ds<BR>which is given by the radial scale = 1/M^2. Solve the direct geodesic<BR>problem from the center using length s and azimuth alpha to recover the<BR>original point.<BR><BR>I just checked in a Gnomonic class for GeographicLib that implements the<BR>forward and reverse projections. It's available through svn from<BR>SourceForge.<BR><BR>This projection minimizes the error in the straightness of geodesics<BR>near the center. If, instead, you want to minimize the error over some<BR>region, then use a two-point azimuthal projection with two base points<BR>judiciously chosen within the region.<BR><BR>Questions:<BR><BR>(1) Is this new?<BR>(2) Is Bagratuni this first person to define M?<BR>(3) Is GeographicLib the first place where the formula for M for an<BR> ellipsoid is given?<BR>(4) Does anyone have a better copy of Bagratuni than the one given<BR> above. (I'll even take a good copy in the Russian original.)<BR><BR>--<BR>Charles Karney <ckarney@sarnoff.com><BR>Sarnoff Corporation, Princeton, NJ 08543-5300<BR><BR>Tel: +1 609 734 2312<BR>Fax: +1 609 734 2662<BR>_______________________________________________<BR>Proj mailing list<BR>Proj@lists.maptools.org<BR><A href="http://lists.maptools.org/mailman/listinfo/proj">http://lists.maptools.org/mailman/listinfo/proj</A><BR></FONT></P></DIV></BODY></HTML>