LambertConformalConic.hpp

Go to the documentation of this file.
00001 /**
00002  * \file LambertConformalConic.hpp
00003  * \brief Header for GeographicLib::LambertConformalConic class
00004  *
00005  * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP)
00011 #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP "$Id: LambertConformalConic.hpp 6816 2010-02-05 21:03:10Z karney $"
00012 
00013 #include "GeographicLib/Constants.hpp"
00014 #include <algorithm>
00015 
00016 namespace GeographicLib {
00017 
00018   /**
00019    * \brief Lambert Conformal Conic Projection
00020    *
00021    * Implementation taken from the report,
00022    * - J. P. Snyder,
00023    *   <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projections: A
00024    *   Working Manual</a>, USGS Professional Paper 1395 (1987),
00025    *   pp. 107&ndash;109.
00026    *
00027    * This is a straightforward implementation of the equations in Snyder except
00028    * that Newton's method is used to invert the projection and that several of
00029    * the formulas are modified so that the projection correctly degenerates to
00030    * the Mercator projection or the polar sterographic projection.
00031    *
00032    * The ellipsoid parameters, the standard parallels, and the scale on the
00033    * standard parallels are set in the constructor.  If the standard parallels
00034    * are both at a single pole the projection becomes the polar stereographic
00035    * projection.  If the standard parallels are symmetric around equator, the
00036    * projection becomes the Mercator projection.  The central meridian (which
00037    * is a trivial shift of the longitude) is specified as the \e lon0 argument
00038    * of the LambertConformalConic::Forward and LambertConformalConic::Reverse
00039    * functions.  The latitude of origin is taken to be the latitude of tangency
00040    * which lies between the standard parallels which is given by
00041    * LambertConformalConic::OriginLatitude.  There is no provision in this
00042    * class for specifying a false easting or false northing or a different
00043    * latitude of origin.  However these are can be simply included by the
00044    * calling function.  For example the Pennsylvania South state coordinate
00045    * system (<a href="http://www.spatialreference.org/ref/epsg/3364/">
00046    * EPSG:3364</a>) is obtained by:
00047    \verbatim
00048    const double
00049      a = GeographicLib::Constants::WGS84_a(), r = 298.257222101, // GRS80
00050      lat1 = 39 + 56/60.0, lat1 = 40 + 58/60.0, // standard parallels
00051      k1 = 1,                                   // scale
00052      lat0 = 39 + 20/60.0, lon0 = 75 + 45/60.0, // origin
00053      fe = 600000, fn = 0;                      // false easting and northing
00054    // Set up basic projection
00055    const GeographicLib::LambertConformalConic PASouth(a, r, lat1, lat2, k1);
00056    double x0, y0;
00057    {
00058      double gamma, k;
00059      // Transform origin point
00060      PASouth.Forward(lon0, lat0, lon0, x0, y0, gamma, k);
00061      x0 -= fe; y0 -= fn;         // Combine result with false origin
00062    }
00063    double lat, lon, x, y, gamma, k;
00064    // Sample conversion from geodetic to PASouth grid
00065    std::cin >> lat >> lon;
00066    PASouth.Forward(lon0, lat, lon, x, y, gamma, k);
00067    x -= x0; y -= y0;
00068    std::cout << x << " " << y << "\n";
00069    // Sample conversion from PASouth grid to geodetic
00070    std::cin >> x >> y;
00071    x += x0; y += y0;
00072    PASouth.Reverse(lon0, x, y, lat, lon, gamma, k);
00073    std::cout << lat << " " << lon << "\n";
00074    \endverbatim
00075    **********************************************************************/
00076   class LambertConformalConic {
00077   private:
00078     typedef Math::real real;
00079     const real _a, _r, _f, _e2, _e, _e2m;
00080     real _sign, _n, _nc, _lt0, _t0n, _t0nm1, _scale, _lat0, _k0;
00081     static const real eps, eps2, epsx, tol, ahypover;
00082     static const int numit = 5;
00083     static inline real sq(real x) throw() { return x * x; }
00084     // e * atanh(e * x) = log( ((1 + e*x)/(1 - e*x))^(e/2) ) if f >= 0
00085     // - sqrt(-e2) * atan( sqrt(-e2) * x)                    if f < 0
00086     inline real eatanhe(real x) const throw() {
00087       return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x);
00088     }
00089     inline real mf(real sphi, real cphi) const throw() {
00090       return cphi/std::sqrt(1 - _e2 * sq(sphi)); // Snyder's m, p 108, eq 14-15
00091     }
00092     inline real tf(real sphi, real cphi) const throw() {
00093       // Snyder's t, p 108, eq 15-9a
00094       // First factor is sqrt((1 - sphi) / (1 + sphi))
00095       // Second factor is ((1 + e * sphi)/(1 - e * sphi)) ^ (e/2)
00096       return (sphi >= 0 ? cphi / (1 + sphi) : (1 - sphi) / cphi) *
00097         std::exp(eatanhe(sphi));
00098     }
00099     inline real logtf(real sphi, real cphi) const throw() {
00100       // Return log(t).  This retains relative accuracy near the equator where
00101       // t -> 1.  This is the negative of the standard Mercator northing.  Note
00102       // that log( sqrt((1 - sin(phi)) / (1 + sin(phi))) ) = -asinh(tan(phi))
00103       return -Math::asinh(sphi/std::max(epsx, cphi)) + eatanhe(sphi);
00104     }
00105     inline real logmtf(real sphi) const throw() {
00106       // Return log(m/t).  This allows the cancellation of cphi in m and t.
00107       return std::log((1 + sphi)/std::sqrt(1 - _e2 * sq(sphi))) -
00108         eatanhe(sphi);
00109     }
00110     void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) throw();
00111   public:
00112 
00113     /**
00114      * Constructor for a ellipsoid radius \e a (meters), reciprocal flattening
00115      * \e r, standard parallel (the circle of tangency) \e stdlat (degrees),
00116      * and scale on the standard parallel \e k0.  Setting \e r = 0 implies \e r
00117      * = inf or flattening = 0 (i.e., a sphere). An exception is thrown if \e a
00118      * or \e k0 is not positive or if \e stdlat is not in the range [-90, 90].
00119      **********************************************************************/
00120     LambertConformalConic(real a, real r, real stdlat, real k0);
00121 
00122     /**
00123      * Constructor for a ellipsoid radius \e a (meters), reciprocal flattening
00124      * \e r, standard parallels \e stdlat1 (degrees) and \e stdlat2 (degrees),
00125      * and the scale on the standard parallels \e k1.  Setting \e r = 0 implies
00126      * \e r = inf or flattening = 0 (i.e., a sphere). An exception is thrown if
00127      * \e a or \e k0 is not positive or if \e stdlat1 or \e stdlat2 is not in
00128      * the range [-90, 90].  In addition, if either \e stdlat1 or \e stdlat2 is
00129      * a pole, then an exception is thrown if \e stdlat1 is not equal \e
00130      * stdlat2
00131      **********************************************************************/
00132     LambertConformalConic(real a, real r, real stdlat1, real stdlat2, real k1);
00133 
00134     /**
00135      * An alternative constructor for 2 standard parallels where the parallels
00136      * are given by their sines and cosines.  This allows parallels close to
00137      * the poles to be specified accurately.
00138      **********************************************************************/
00139     LambertConformalConic(real a, real r,
00140                           real sinlat1, real coslat1,
00141                           real sinlat2, real coslat2,
00142                           real k1);
00143 
00144     /**
00145      * Alter the scale for the projection so that on latitude \e lat, the scale
00146      * is \e k (default 1).  The allows a "latitude of true scale" to be
00147      * specified.  An exception is thrown if \e k is not positive.
00148      **********************************************************************/
00149     void SetScale(real lat, real k = real(1));
00150 
00151     /**
00152      * Convert from latitude \e lat (degrees) and longitude \e lon (degrees) to
00153      * Lambert conformal conic easting \e x (meters) and northing \e y
00154      * (meters).  The central meridian is given by \e lon0 (degrees) and the
00155      * latitude origin is given by LambertConformalConic::LatitudeOrigin().
00156      * Also return the meridian convergence \e gamma (degrees) and the scale \e
00157      * k.  No false easting or northing is added and \e lat should be in the
00158      * range [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360].
00159      * The values of \e x and \e y returned for points which project to
00160      * infinity (i.e., one or both of the poles) will be large but finite.  The
00161      * value of \e k returned for the poles in infinite (unless \e lat equals
00162      * the latitude origin).
00163      **********************************************************************/
00164     void Forward(real lon0, real lat, real lon,
00165                  real& x, real& y, real& gamma, real& k) const throw();
00166 
00167     /**
00168      * Convert from Lambert conformal conic easting \e x (meters) and northing
00169      * \e y (meters) to latitude \e lat (degrees) and longitude \e lon
00170      * (degrees) .  The central meridian is given by \e lon0 (degrees) and the
00171      * latitude origin is given by LambertConformalConic::LatitudeOrigin().
00172      * Also return the meridian convergence \e gamma (degrees) and the scale \e
00173      * k.  No false easting or northing is added.  \e lon0 should be in the
00174      * range [-180, 360].  The value of \e lon returned is in the range [-180,
00175      * 180).
00176      **********************************************************************/
00177     void Reverse(real lon0, real x, real y,
00178                  real& lat, real& lon, real& gamma, real& k) const throw();
00179 
00180     /**
00181      * The major radius of the ellipsoid (meters).  This is that value of \e a
00182      * used in the constructor.
00183      **********************************************************************/
00184     Math::real MajorRadius() const throw() { return _a; }
00185 
00186     /**
00187      * The inverse flattening of the ellipsoid.  This is that value of \e r
00188      * used in the constructor.  A value of 0 is returned for a sphere
00189      * (infinite inverse flattening).
00190      **********************************************************************/
00191     Math::real InverseFlattening() const throw() { return _r; }
00192 
00193     /**
00194      * The latitude of the origin for the projection (degrees).  This is that
00195      * latitude of minimum scale and equals the \e stdlat in the 1-parallel
00196      * constructor and lies between \e stdlat1 and \e stdlat2 in the 2-parallel
00197      * constructors.
00198      **********************************************************************/
00199     Math::real OriginLatitude() const throw() { return _lat0; }
00200 
00201     /**
00202      * The central scale for the projection.  This is that scale on the
00203      * latitude of origin..
00204      **********************************************************************/
00205     Math::real CentralScale() const throw() { return _k0; }
00206 
00207     /**
00208      * A global instantiation of LambertConformalConic with the WGS84 ellipsoid
00209      * and the UPS scale factor and \e stdlat = 0.  This degenerates to the
00210      * Mercator projection.
00211      **********************************************************************/
00212     const static LambertConformalConic Mercator;
00213   };
00214 
00215 } // namespace GeographicLib
00216 
00217 #endif

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1