src/map/proj/polarstereographic.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 "polarstereographic.h"
  46
  47
  48 #define POLAR_POW(EsSin) pow((1.0 - EsSin) / (1.0 + EsSin), _es_OVER_2)
  49
  50 PolarStereographic::PolarStereographic(const Ellipsoid &ellipsoid,
  51   double latitudeOrigin, double longitudeOrigin,
  52   double falseEasting, double falseNorthing)
  53 {
  54         _two_a = 2.0 * ellipsoid.radius();
  55
  56         if (longitudeOrigin > M_PI)
  57                 longitudeOrigin -= 2*M_PI;
  58
  59         if (latitudeOrigin < 0) {
  60                 _southernHemisphere = 1;
  61                 _originLatitude = -deg2rad(latitudeOrigin);
  62                 _originLongitude = -deg2rad(longitudeOrigin);
  63         } else {
  64                 _southernHemisphere = 0;
  65                 _originLatitude = deg2rad(latitudeOrigin);
  66                 _originLongitude = deg2rad(longitudeOrigin);
  67         }
  68         _falseEasting = falseEasting;
  69         _falseNorthing = falseNorthing;
  70
  71         double es2 = ellipsoid.es();
  72         _es = sqrt(es2);
  73         _es_OVER_2 = _es / 2.0;
  74
  75         if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10) {
  76                 double slat = sin(_originLatitude);
  77                 double essin = _es * slat;
  78                 double pow_es = POLAR_POW(essin);
  79                 double clat = cos(_originLatitude);
  80                 _mc = clat / sqrt(1.0 - essin * essin);
  81                 _a_mc = ellipsoid.radius() * _mc;
  82                 _tc = tan(M_PI_4 - _originLatitude / 2.0) / pow_es;
  83         } else {
  84                 double one_PLUS_es = 1.0 + _es;
  85                 double one_MINUS_es = 1.0 - _es;
  86                 _e4 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es,
  87                   one_MINUS_es));
  88         }
  89 }
  90
  91 PointD PolarStereographic::ll2xy(const Coordinates &c) const
  92 {
  93         double Easting, Northing;
  94         double Longitude = deg2rad(c.lon());
  95         double Latitude = deg2rad(c.lat());
  96
  97
  98         if (fabs(fabs(Latitude) - M_PI_2) < 1.0e-10) {
  99                 Easting = 0.0;
 100                 Northing = 0.0;
 101         } else {
 102                 if (_southernHemisphere != 0) {
 103                         Longitude *= -1.0;
 104                         Latitude *= -1.0;
 105                 }
 106
 107                 double dlam = Longitude - _originLongitude;
 108                 if (dlam > M_PI)
 109                         dlam -= 2*M_PI;
 110                 if (dlam < -M_PI)
 111                         dlam += 2*M_PI;
 112
 113                 double slat = sin(Latitude);
 114                 double essin = _es * slat;
 115                 double pow_es = POLAR_POW(essin);
 116                 double t = tan(M_PI_4 - Latitude / 2.0) / pow_es;
 117
 118                 double rho;
 119                 if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10)
 120                         rho = _a_mc * t / _tc;
 121                 else
 122                         rho = _two_a * t / _e4;
 123
 124                 Easting = rho * sin(dlam) + _falseEasting;
 125
 126                 if (_southernHemisphere != 0) {
 127                         Easting *= -1.0;
 128                         Northing = rho * cos(dlam) + _falseNorthing;
 129                 } else
 130                         Northing = -rho * cos(dlam) + _falseNorthing;
 131         }
 132
 133         return PointD(Easting, Northing);
 134 }
 135
 136 Coordinates PolarStereographic::xy2ll(const PointD &p) const
 137 {
 138         double Latitude, Longitude;
 139         double dy = p.y() - _falseNorthing;
 140         double dx = p.x() - _falseEasting;
 141
 142         if ((dy == 0.0) && (dx == 0.0)) {
 143                 Latitude = M_PI_2;
 144                 Longitude = _originLongitude;
 145         } else {
 146                 if (_southernHemisphere != 0) {
 147                         dy *= -1.0;
 148                         dx *= -1.0;
 149                 }
 150
 151                 double rho = sqrt(dx * dx + dy * dy);
 152                 double t;
 153                 if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10)
 154                         t = rho * _tc / (_a_mc);
 155                 else
 156                         t = rho * _e4 / (_two_a);
 157                 double PHI = M_PI_2 - 2.0 * atan(t);
 158
 159                 double tempPHI = 0.0;
 160                 while (fabs(PHI - tempPHI) > 1.0e-10) {
 161                         tempPHI = PHI;
 162                         double sin_PHI = sin(PHI);
 163                         double essin =  _es * sin_PHI;
 164                         double pow_es = POLAR_POW(essin);
 165                         PHI = M_PI_2 - 2.0 * atan(t * pow_es);
 166                 }
 167                 Latitude = PHI;
 168                 Longitude = _originLongitude + atan2(dx, -dy);
 169
 170                 if (Longitude > M_PI)
 171                         Longitude -= 2*M_PI;
 172                 else if (Longitude < -M_PI)
 173                         Longitude += 2*M_PI;
 174
 175                 if (Latitude > M_PI_2)
 176                         Latitude = M_PI_2;
 177                 else if (Latitude < -M_PI_2)
 178                         Latitude = -M_PI_2;
 179
 180                 if (Longitude > M_PI)
 181                         Longitude = M_PI;
 182                 else if (Longitude < -M_PI)
 183                         Longitude = -M_PI;
 184         }
 185
 186         if (_southernHemisphere != 0) {
 187                 Latitude *= -1.0;
 188                 Longitude *= -1.0;
 189         }
 190
 191         return Coordinates(rad2deg(Longitude), rad2deg(Latitude));
 192 }
 193
 194 bool PolarStereographic::operator==(const CT &ct) const
 195 {
 196         const PolarStereographic *other
 197           = dynamic_cast<const PolarStereographic*>(&ct);
 198         return (other != 0 && _originLatitude == other->_originLatitude
 199           && _originLongitude == other->_originLongitude
 200           && _falseEasting == other->_falseEasting
 201           && _falseNorthing == other->_falseNorthing && _two_a == other->_two_a
 202           && _es == other->_es && _southernHemisphere == other->_southernHemisphere);
 203 }