<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.&nbsp; I do not have a copy, though.&nbsp; 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>&nbsp;</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&nbsp; 70803<BR>Voice and Facsimile:&nbsp; (225) 578-8536 [Academic] <BR>Voice and Facsimile:&nbsp; (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)&nbsp;238-8975 [Academic &amp; Research]<BR></SPAN>Honorary Life Member of the <BR>Louisiana Society of Professional Surveyors <BR>Fellow&nbsp;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.&nbsp; 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.&nbsp; However, we can<BR>obtain a projection where the projected geodesics are approximately<BR>straight close to the center of the projection.&nbsp; It is obtained as the<BR>limit of a 2-point azimuthal projection as the two points approach one<BR>another.&nbsp; 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>&nbsp;(1) azimuthal, all azimuths from the center are correct;<BR>&nbsp;(2) hence, all lines through the center are geodesics;<BR>&nbsp;(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.&nbsp; Take an arbitrary pair of points within r =<BR>1000 km of the center.&nbsp; Draw a straight line between these on the map;<BR>then back project this line onto the ellipsoid.&nbsp; How close is this line<BR>to a geodesic?&nbsp; I find that *at worst*,<BR><BR>&nbsp; max error in initial/final azimuth = 1.0" = 0.39 * f*(r/a)^3<BR>&nbsp; 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).&nbsp; 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>&nbsp; max error in initial/final azimuth = 108" = 1.0 * f*(r/a)<BR>&nbsp; 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>&nbsp; * m, the reduced length.&nbsp; Take two geodesics of length s which start<BR>&nbsp;&nbsp;&nbsp; from the center with azimuths differing by dalpha.&nbsp; The endpoints<BR>&nbsp;&nbsp;&nbsp; are separated by m * dalpha.<BR><BR>&nbsp; * M, the geodesic scale.&nbsp; Take two geodesics of length s which are<BR>&nbsp;&nbsp;&nbsp; parallel close to the center and initially separated by dt.&nbsp; The<BR>&nbsp;&nbsp;&nbsp; endpoints are separated by M * dt.<BR><BR>The reduced length was introduced by Gauss (1827) and Christoffel<BR>(1868).&nbsp; Helmert (1880), Sec. 6.5, gives an explicit formula for the<BR>ellipsoid.&nbsp; The geodesic scale was introduced (I think) by<BR><BR>&nbsp; G. V. Bagratuni,<BR>&nbsp; Course in spheroidal geodesy,<BR>&nbsp; FTD-MT-64-390 (US Air Force, 1967), Sec. 17<BR>&nbsp; [Kurs sferoidicheskoi geodezii, (Moscow, 1962)].<BR>&nbsp; <A href="http://handle.dtic.mil/100.2/AD650520">http://handle.dtic.mil/100.2/AD650520</A><BR><BR>(Bagratuni uses the symbol n.)&nbsp; 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.&nbsp; Let alpha be the<BR>azimuth of the geodesic at the center and s be its length.&nbsp; Compute<BR><BR>&nbsp; rho = m/M<BR><BR>Then the projected point is<BR><BR>&nbsp; x = rho * sin(alpha); y = rho * cos(alpha)<BR><BR>The radial scale is 1/M^2 and the transverse scale is 1/M.&nbsp; For a<BR>sphere, m = a * sin(s/a) and M = cos(s/a), thus<BR><BR>&nbsp; rho = a * tan(s/a)<BR><BR>which gives the spherical gnomonic projection.<BR><BR>In order to reverse the projection, compute<BR><BR>&nbsp; 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.&nbsp; 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.&nbsp; It's available through svn from<BR>SourceForge.<BR><BR>This projection minimizes the error in the straightness of geodesics<BR>near the center.&nbsp; 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>&nbsp;&nbsp;&nbsp; ellipsoid is given?<BR>(4) Does anyone have a better copy of Bagratuni than the one given<BR>&nbsp;&nbsp;&nbsp; above.&nbsp; (I'll even take a good copy in the Russian original.)<BR><BR>--<BR>Charles Karney &lt;ckarney@sarnoff.com&gt;<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>