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