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