src/map/proj/mercator.cpp

   1 /*
   2  * Based on libgeotrans with the following Source Code Disclaimer:
   3
   4 1. The GEOTRANS source code ("the software") is provided free of charge by
   5 the National Imagery and Mapping Agency (NIMA) of the United States
   6 Department of Defense. Although NIMA makes no copyright claim under Title 17
   7 U.S.C., NIMA claims copyrights in the source code under other legal regimes.
   8 NIMA hereby grants to each user of the software a license to use and
   9 distribute the software, and develop derivative works.
  10
  11 2. Warranty Disclaimer: The software was developed to meet only the internal
  12 requirements of the U.S. National Imagery and Mapping Agency. The software
  13 is provided "as is," and no warranty, express or implied, including but not
  14 limited to the implied warranties of merchantability and fitness for
  15 particular purpose or arising by statute or otherwise in law or from a
  16 course of dealing or usage in trade, is made by NIMA as to the accuracy and
  17 functioning of the software.
  18
  19 3. NIMA and its personnel are not required to provide technical support or
  20 general assistance with respect to the software.
  21
  22 4. Neither NIMA nor its personnel will be liable for any claims, losses, or
  23 damages arising from or connected with the use of the software. The user
  24 agrees to hold harmless the United States National Imagery and Mapping
  25 Agency. The user's sole and exclusive remedy is to stop using the software.
  26
  27 5. NIMA requests that products developed using the software credit the
  28 source of the software with the following statement, "The product was
  29 developed using GEOTRANS, a product of the National Imagery and Mapping
  30 Agency and U.S. Army Engineering Research and Development Center."
  31
  32 6. For any products developed using the software, NIMA requires a disclaimer
  33 that use of the software does not indicate endorsement or approval of the
  34 product by the Secretary of Defense or the National Imagery and Mapping
  35 Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
  36 the National Imagery and Mapping Agency, the initials "NIMA", the seal of
  37 the National Imagery and Mapping Agency, or any colorable imitation thereof
  38 shall not be used to imply approval, endorsement, or authorization of a
  39 product without prior written permission from United States Secretary of
  40 Defense.
  41
  42 */
  43
  44 #include "map/ellipsoid.h"
  45 #include "mercator.h"
  46
  47 Mercator::Mercator(const Ellipsoid &ellipsoid, double latitudeOrigin,
  48   double longitudeOrigin, double falseEasting, double falseNorthing)
  49 {
  50         double es = ellipsoid.es();
  51         double es2;
  52         double es3;
  53         double es4;
  54         double sin_olat;
  55
  56         _latitudeOrigin = deg2rad(latitudeOrigin);
  57         _longitudeOrigin = deg2rad(longitudeOrigin);
  58         if (_longitudeOrigin > M_PI)
  59                 _longitudeOrigin -= 2 * M_PI;
  60         _falseNorthing = falseNorthing;
  61         _falseEasting = falseEasting;
  62
  63         _a = ellipsoid.radius();
  64         _e = sqrt(es);
  65
  66         sin_olat = sin(_latitudeOrigin);
  67         _scaleFactor = 1.0 / (sqrt(1.e0 - es * sin_olat * sin_olat)
  68           / cos(_latitudeOrigin));
  69         es2 = es * es;
  70         es3 = es2 * es;
  71         es4 = es3 * es;
  72         _ab = es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0 + 13.e0 * es4 / 360.e0;
  73         _bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0 + 811.e0 * es4 / 11520.e0;
  74         _cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
  75         _db = 4279.e0 * es4 / 161280.e0;
  76 }
  77
  78 PointD Mercator::ll2xy(const Coordinates &c) const
  79 {
  80         double lon = deg2rad(c.lon());
  81         double lat = deg2rad(c.lat());
  82         double ctanz2;
  83         double e_x_sinlat;
  84         double delta_lon;
  85         double tan_temp;
  86         double pow_temp;
  87
  88         if (lon > M_PI)
  89                 lon -= 2 * M_PI;
  90         e_x_sinlat = _e * sin(lat);
  91         tan_temp = tan(M_PI_4 + lat / 2.e0);
  92         pow_temp = pow((1.e0 - e_x_sinlat) / (1.e0 + e_x_sinlat), _e / 2.e0);
  93         ctanz2 = tan_temp * pow_temp;
  94         delta_lon = lon - _longitudeOrigin;
  95         if (delta_lon > M_PI)
  96           delta_lon -= 2 * M_PI;
  97         if (delta_lon < -M_PI)
  98           delta_lon += 2 * M_PI;
  99
 100         return PointD(_scaleFactor * _a * delta_lon + _falseEasting,
 101           _scaleFactor * _a * log(ctanz2) + _falseNorthing);
 102 }
 103
 104 Coordinates Mercator::xy2ll(const PointD &p) const
 105 {
 106         double dx;
 107         double dy;
 108         double xphi;
 109         double lat, lon;
 110
 111         dy = p.y() - _falseNorthing;
 112         dx = p.x() - _falseEasting;
 113         lon = _longitudeOrigin + dx / (_scaleFactor * _a);
 114         xphi = M_PI_2 - 2.e0 * atan(1.e0 / exp(dy / (_scaleFactor * _a)));
 115         lat = xphi + _ab * sin(2.e0 * xphi) + _bb * sin(4.e0 * xphi)
 116           + _cb * sin(6.e0 * xphi) + _db * sin(8.e0 * xphi);
 117         if (lon > M_PI)
 118                 lon -= 2 * M_PI;
 119         if (lon < -M_PI)
 120                 lon += 2 * M_PI;
 121
 122         return Coordinates(rad2deg(lon), rad2deg(lat));
 123 }
 124
 125 bool Mercator::operator==(const CT &ct) const
 126 {
 127         const Mercator *other = dynamic_cast<const Mercator*>(&ct);
 128         return (other != 0 && _a == other->_a && _e == other->_e
 129           && _latitudeOrigin == other->_latitudeOrigin
 130           && _longitudeOrigin == other->_longitudeOrigin
 131           && _falseNorthing == other->_falseNorthing
 132           && _falseEasting == other->_falseEasting);
 133 }