00001 /** 00002 * \file Gnomonic.hpp 00003 * \brief Header for GeographicLib::Gnomonic class 00004 * 00005 * Copyright (c) Charles Karney (2010-2011) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_GNOMONIC_HPP) 00011 #define GEOGRAPHICLIB_GNOMONIC_HPP 1 00012 00013 #include <GeographicLib/Geodesic.hpp> 00014 #include <GeographicLib/GeodesicLine.hpp> 00015 #include <GeographicLib/Constants.hpp> 00016 00017 namespace GeographicLib { 00018 00019 /** 00020 * \brief %Gnomonic projection 00021 * 00022 * %Gnomonic projection centered at an arbitrary position \e C on the 00023 * ellipsoid. This projection is derived in Section 8 of 00024 * - C. F. F. Karney, 00025 * <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> 00026 * Algorithms for geodesics</a>, 00027 * J. Geodesy <b>87</b>, 43--55 (2013); 00028 * DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> 00029 * 10.1007/s00190-012-0578-z</a>; 00030 * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html"> 00031 * geod-addenda.html</a>. 00032 * . 00033 * The projection of \e P is defined as follows: compute the geodesic line 00034 * from \e C to \e P; compute the reduced length \e m12, geodesic scale \e 00035 * M12, and ρ = <i>m12</i>/\e M12; finally \e x = ρ sin \e azi1; \e 00036 * y = ρ cos \e azi1, where \e azi1 is the azimuth of the geodesic at \e 00037 * C. The Gnomonic::Forward and Gnomonic::Reverse methods also return the 00038 * azimuth \e azi of the geodesic at \e P and reciprocal scale \e rk in the 00039 * azimuthal direction. The scale in the radial direction if 00040 * 1/<i>rk</i><sup>2</sup>. 00041 * 00042 * For a sphere, ρ is reduces to \e a tan(<i>s12</i>/<i>a</i>), where \e 00043 * s12 is the length of the geodesic from \e C to \e P, and the gnomonic 00044 * projection has the property that all geodesics appear as straight lines. 00045 * For an ellipsoid, this property holds only for geodesics interesting the 00046 * centers. However geodesic segments close to the center are approximately 00047 * straight. 00048 * 00049 * Consider a geodesic segment of length \e l. Let \e T be the point on the 00050 * geodesic (extended if necessary) closest to \e C the center of the 00051 * projection and \e t be the distance \e CT. To lowest order, the maximum 00052 * deviation (as a true distance) of the corresponding gnomonic line segment 00053 * (i.e., with the same end points) from the geodesic is<br> 00054 * <br> 00055 * (<i>K</i>(<i>T</i>) - <i>K</i>(<i>C</i>)) 00056 * <i>l</i><sup>2</sup> \e t / 32.<br> 00057 * <br> 00058 * where \e K is the Gaussian curvature. 00059 * 00060 * This result applies for any surface. For an ellipsoid of revolution, 00061 * consider all geodesics whose end points are within a distance \e r of \e 00062 * C. For a given \e r, the deviation is maximum when the latitude of \e C 00063 * is 45°, when endpoints are a distance \e r away, and when their 00064 * azimuths from the center are ± 45° or ± 135°. 00065 * To lowest order in \e r and the flattening \e f, the deviation is \e f 00066 * (<i>r</i>/2<i>a</i>)<sup>3</sup> \e r. 00067 * 00068 * The conversions all take place using a Geodesic object (by default 00069 * Geodesic::WGS84()). For more information on geodesics see \ref geodesic. 00070 * 00071 * <b>CAUTION:</b> The definition of this projection for a sphere is 00072 * standard. However, there is no standard for how it should be extended to 00073 * an ellipsoid. The choices are: 00074 * - Declare that the projection is undefined for an ellipsoid. 00075 * - Project to a tangent plane from the center of the ellipsoid. This 00076 * causes great ellipses to appear as straight lines in the projection; 00077 * i.e., it generalizes the spherical great circle to a great ellipse. 00078 * This was proposed by independently by Bowring and Williams in 1997. 00079 * - Project to the conformal sphere with the constant of integration chosen 00080 * so that the values of the latitude match for the center point and 00081 * perform a central projection onto the plane tangent to the conformal 00082 * sphere at the center point. This causes normal sections through the 00083 * center point to appear as straight lines in the projection; i.e., it 00084 * generalizes the spherical great circle to a normal section. This was 00085 * proposed by I. G. Letoval'tsev, Generalization of the %Gnomonic 00086 * Projection for a Spheroid and the Principal Geodetic Problems Involved 00087 * in the Alignment of Surface Routes, Geodesy and Aerophotography (5), 00088 * 271--274 (1963). 00089 * - The projection given here. This causes geodesics close to the center 00090 * point to appear as straight lines in the projection; i.e., it 00091 * generalizes the spherical great circle to a geodesic. 00092 * 00093 * Example of use: 00094 * \include example-Gnomonic.cpp 00095 * 00096 * <a href="GeodesicProj.1.html">GeodesicProj</a> is a command-line utility 00097 * providing access to the functionality of AzimuthalEquidistant, Gnomonic, 00098 * and CassiniSoldner. 00099 **********************************************************************/ 00100 00101 class GEOGRAPHICLIB_EXPORT Gnomonic { 00102 private: 00103 typedef Math::real real; 00104 real eps0_, eps_; 00105 Geodesic _earth; 00106 real _a, _f; 00107 static const int numit_ = 10; 00108 public: 00109 00110 /** 00111 * Constructor for Gnomonic. 00112 * 00113 * @param[in] earth the Geodesic object to use for geodesic calculations. 00114 * By default this uses the WGS84 ellipsoid. 00115 **********************************************************************/ 00116 explicit Gnomonic(const Geodesic& earth = Geodesic::WGS84()); 00117 00118 /** 00119 * Forward projection, from geographic to gnomonic. 00120 * 00121 * @param[in] lat0 latitude of center point of projection (degrees). 00122 * @param[in] lon0 longitude of center point of projection (degrees). 00123 * @param[in] lat latitude of point (degrees). 00124 * @param[in] lon longitude of point (degrees). 00125 * @param[out] x easting of point (meters). 00126 * @param[out] y northing of point (meters). 00127 * @param[out] azi azimuth of geodesic at point (degrees). 00128 * @param[out] rk reciprocal of azimuthal scale at point. 00129 * 00130 * \e lat0 and \e lat should be in the range [−90°, 90°] and 00131 * \e lon0 and \e lon should be in the range [−540°, 540°). 00132 * The scale of the projection is 1/<i>rk</i><sup>2</sup> in the "radial" 00133 * direction, \e azi clockwise from true north, and is 1/\e rk in the 00134 * direction perpendicular to this. If the point lies "over the horizon", 00135 * i.e., if \e rk ≤ 0, then NaNs are returned for \e x and \e y (the 00136 * correct values are returned for \e azi and \e rk). A call to Forward 00137 * followed by a call to Reverse will return the original (\e lat, \e lon) 00138 * (to within roundoff) provided the point in not over the horizon. 00139 **********************************************************************/ 00140 void Forward(real lat0, real lon0, real lat, real lon, 00141 real& x, real& y, real& azi, real& rk) const; 00142 00143 /** 00144 * Reverse projection, from gnomonic to geographic. 00145 * 00146 * @param[in] lat0 latitude of center point of projection (degrees). 00147 * @param[in] lon0 longitude of center point of projection (degrees). 00148 * @param[in] x easting of point (meters). 00149 * @param[in] y northing of point (meters). 00150 * @param[out] lat latitude of point (degrees). 00151 * @param[out] lon longitude of point (degrees). 00152 * @param[out] azi azimuth of geodesic at point (degrees). 00153 * @param[out] rk reciprocal of azimuthal scale at point. 00154 * 00155 * \e lat0 should be in the range [−90°, 90°] and \e 00156 * lon0 should be in the range [−540°, 540°). \e lat 00157 * will be in the range [−90°, 90°] and \e lon will 00158 * be in the range [−180°, 180°). The scale of the 00159 * projection is 1/<i>rk</i><sup>2</sup> in the "radial" direction, \e azi 00160 * clockwise from true north, and is 1/\e rk in the direction perpendicular 00161 * to this. Even though all inputs should return a valid \e lat and \e 00162 * lon, it's possible that the procedure fails to converge for very large 00163 * \e x or \e y; in this case NaNs are returned for all the output 00164 * arguments. A call to Reverse followed by a call to Forward will return 00165 * the original (\e x, \e y) (to roundoff). 00166 **********************************************************************/ 00167 void Reverse(real lat0, real lon0, real x, real y, 00168 real& lat, real& lon, real& azi, real& rk) const; 00169 00170 /** 00171 * Gnomonic::Forward without returning the azimuth and scale. 00172 **********************************************************************/ 00173 void Forward(real lat0, real lon0, real lat, real lon, 00174 real& x, real& y) const { 00175 real azi, rk; 00176 Forward(lat0, lon0, lat, lon, x, y, azi, rk); 00177 } 00178 00179 /** 00180 * Gnomonic::Reverse without returning the azimuth and scale. 00181 **********************************************************************/ 00182 void Reverse(real lat0, real lon0, real x, real y, 00183 real& lat, real& lon) const { 00184 real azi, rk; 00185 Reverse(lat0, lon0, x, y, lat, lon, azi, rk); 00186 } 00187 00188 /** \name Inspector functions 00189 **********************************************************************/ 00190 ///@{ 00191 /** 00192 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00193 * the value inherited from the Geodesic object used in the constructor. 00194 **********************************************************************/ 00195 Math::real MajorRadius() const { return _earth.MajorRadius(); } 00196 00197 /** 00198 * @return \e f the flattening of the ellipsoid. This is the value 00199 * inherited from the Geodesic object used in the constructor. 00200 **********************************************************************/ 00201 Math::real Flattening() const { return _earth.Flattening(); } 00202 ///@} 00203 00204 /// \cond SKIP 00205 /** 00206 * <b>DEPRECATED</b> 00207 * @return \e r the inverse flattening of the ellipsoid. 00208 **********************************************************************/ 00209 Math::real InverseFlattening() const 00210 { return _earth.InverseFlattening(); } 00211 /// \endcond 00212 }; 00213 00214 } // namespace GeographicLib 00215 00216 #endif // GEOGRAPHICLIB_GNOMONIC_HPP