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–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