GeographicLib  1.21
Gnomonic.hpp
Go to the documentation of this file.
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