src/map/proj/lambertconic.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 <cmath>
  45 #include "map/ellipsoid.h"
  46 #include "lambertconic.h"
  47
  48
  49 #define LAMBERT_m(clat, essin) (clat / sqrt(1.0 - essin * essin))
  50 #define LAMBERT2_t(lat, essin, es_over_2) \
  51         (tan(M_PI_4 - lat / 2) * pow((1.0 + essin) / (1.0 - essin), es_over_2))
  52 #define LAMBERT1_t(lat, essin, es_over_2) \
  53         (tan(M_PI_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_over_2))
  54
  55
  56 LambertConic1::LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
  57   double longitudeOrigin, double scale, double falseEasting,
  58   double falseNorthing)
  59 {
  60         double e_sin;
  61         double m0;
  62         double lat_orig;
  63
  64
  65         lat_orig = deg2rad(latitudeOrigin);
  66         _longitudeOrigin = deg2rad(longitudeOrigin);
  67         if (_longitudeOrigin > M_PI)
  68                 _longitudeOrigin -= 2 * M_PI;
  69
  70         _falseEasting = falseEasting;
  71         _falseNorthing = falseNorthing;
  72
  73         _e = sqrt(ellipsoid.es());
  74         _e_over_2 = _e / 2.0;
  75
  76         _n = sin(lat_orig);
  77
  78         e_sin = _e * sin(lat_orig);
  79         m0 = LAMBERT_m(cos(lat_orig), e_sin);
  80         _t0 = LAMBERT1_t(lat_orig, e_sin, _e_over_2);
  81
  82         _rho0 = ellipsoid.radius() * scale * m0 / _n;
  83
  84         _rho_olat = _rho0;
  85 }
  86
  87 PointD LambertConic1::ll2xy(const Coordinates &c) const
  88 {
  89         double t;
  90         double rho;
  91         double dlam;
  92         double theta;
  93         double lat = deg2rad(c.lat());
  94
  95
  96         if (fabs(fabs(lat) - M_PI_2) > 1.0e-10) {
  97                 t = LAMBERT1_t(lat, _e * sin(lat), _e_over_2);
  98                 rho = _rho0 * pow(t / _t0, _n);
  99         } else
 100                 rho = 0.0;
 101
 102         dlam = deg2rad(c.lon()) - _longitudeOrigin;
 103
 104         if (dlam > M_PI)
 105                 dlam -= 2 * M_PI;
 106         if (dlam < -M_PI)
 107                 dlam += 2 * M_PI;
 108
 109         theta = _n * dlam;
 110
 111         return PointD(rho * sin(theta) + _falseEasting, _rho_olat - rho
 112           * cos(theta) + _falseNorthing);
 113 }
 114
 115 Coordinates LambertConic1::xy2ll(const PointD &p) const
 116 {
 117         double dx;
 118         double dy;
 119         double rho;
 120         double rho_olat_minus_dy;
 121         double t;
 122         double PHI;
 123         double es_sin;
 124         double tempPHI = 0.0;
 125         double theta = 0.0;
 126         double tolerance = 4.85e-10;
 127         int count = 30;
 128         double lat, lon;
 129
 130
 131         dy = p.y() - _falseNorthing;
 132         dx = p.x() - _falseEasting;
 133         rho_olat_minus_dy = _rho_olat - dy;
 134         rho = sqrt(dx * dx + (rho_olat_minus_dy) * (rho_olat_minus_dy));
 135
 136         if (_n < 0.0) {
 137                 rho *= -1.0;
 138                 dx *= -1.0;
 139                 rho_olat_minus_dy *= -1.0;
 140         }
 141
 142         if (rho != 0.0) {
 143                 theta = atan2(dx, rho_olat_minus_dy) / _n;
 144                 t = _t0 * pow(rho / _rho0, 1 / _n);
 145                 PHI = M_PI_2 - 2.0 * atan(t);
 146                 while (fabs(PHI - tempPHI) > tolerance && count) {
 147                         tempPHI = PHI;
 148                         es_sin = _e * sin(PHI);
 149                         PHI = M_PI_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin),
 150                           _e_over_2));
 151                         count--;
 152                 }
 153
 154                 if (!count)
 155                         return Coordinates();
 156
 157                 lat = PHI;
 158                 lon = theta + _longitudeOrigin;
 159
 160                 if (fabs(lat) < 2.0e-7)
 161                         lat = 0.0;
 162                 if (lat > M_PI_2)
 163                         lat = M_PI_2;
 164                 else if (lat < -M_PI_2)
 165                         lat = -M_PI_2;
 166
 167                 if (lon > M_PI) {
 168                         if (lon - M_PI < 3.5e-6)
 169                                 lon = M_PI;
 170                         else
 171                                 lon -= 2 * M_PI;
 172                 }
 173                 if (lon < -M_PI) {
 174                         if (fabs(lon + M_PI) < 3.5e-6)
 175                                 lon = -M_PI;
 176                         else
 177                                 lon += 2 * M_PI;
 178                 }
 179
 180                 if (fabs(lon) < 2.0e-7)
 181                         lon = 0.0;
 182                 if (lon > M_PI)
 183                         lon = M_PI;
 184                 else if (lon < -M_PI)
 185                         lon = -M_PI;
 186         } else {
 187                 if (_n > 0.0)
 188                         lat = M_PI_2;
 189                 else
 190                         lat = -M_PI_2;
 191                 lon = _longitudeOrigin;
 192         }
 193
 194         return Coordinates(rad2deg(lon), rad2deg(lat));
 195 }
 196
 197 bool LambertConic1::operator==(const CT &ct) const
 198 {
 199         const LambertConic1 *other = dynamic_cast<const LambertConic1*>(&ct);
 200         return (other != 0 && _longitudeOrigin == other->_longitudeOrigin
 201           && _falseEasting == other->_falseEasting
 202           && _falseNorthing == other->_falseNorthing && _e == other->_e
 203           && _n == other->_n && _rho0 == other->_rho0);
 204 }
 205
 206
 207 LambertConic2::LambertConic2(const Ellipsoid &ellipsoid,
 208   double standardParallel1, double standardParallel2, double latitudeOrigin,
 209   double longitudeOrigin, double falseEasting, double falseNorthing)
 210 {
 211         double e, e_over_2, e_sin;
 212         double lat0;
 213         double k0;
 214         double t0;
 215         double t1, t2;
 216         double t_olat;
 217         double m0;
 218         double m1;
 219         double m2;
 220         double n;
 221         double const_value;
 222         double sp1, sp2;
 223         double lat_orig;
 224
 225
 226         lat_orig = deg2rad(latitudeOrigin);
 227         sp1 = deg2rad(standardParallel1);
 228         sp2 = deg2rad(standardParallel2);
 229
 230         if (fabs(sp1 - sp2) > 1.0e-10) {
 231                 e = sqrt(ellipsoid.es());
 232                 e_over_2 = e / 2.0;
 233
 234                 e_sin = e * sin(lat_orig);
 235                 t_olat = LAMBERT2_t(lat_orig, e_sin, e_over_2);
 236
 237                 e_sin = e * sin(sp1);
 238                 m1 = LAMBERT_m(cos(sp1), e_sin);
 239                 t1 = LAMBERT2_t(sp1, e_sin, e_over_2);
 240
 241                 e_sin = e * sin(sp2);
 242                 m2 = LAMBERT_m(cos(sp2), e_sin);
 243                 t2 = LAMBERT2_t(sp2, e_sin, e_over_2);
 244
 245                 n = log(m1 / m2) / log(t1 / t2);
 246
 247                 lat0 = asin(n);
 248
 249                 e_sin = e * sin(lat0);
 250                 m0 = LAMBERT_m(cos(lat0), e_sin);
 251                 t0 = LAMBERT2_t(lat0, e_sin, e_over_2);
 252
 253                 k0 = (m1 / m0) * (pow(t0 / t1, n));
 254
 255                 const_value = ((ellipsoid.radius() * m2) / (n * pow(t2, n)));
 256
 257                 falseNorthing += (const_value * pow(t_olat, n)) - (const_value
 258                   * pow(t0, n));
 259         } else {
 260                 lat0 = sp1;
 261                 k0 = 1.0;
 262         }
 263
 264         _lc1 = LambertConic1(ellipsoid, rad2deg(lat0), longitudeOrigin, k0,
 265           falseEasting, falseNorthing);
 266 }
 267
 268 PointD LambertConic2::ll2xy(const Coordinates &c) const
 269 {
 270         return _lc1.ll2xy(c);
 271 }
 272
 273 Coordinates LambertConic2::xy2ll(const PointD &p) const
 274 {
 275         return _lc1.xy2ll(p);
 276 }
 277
 278 bool LambertConic2::operator==(const CT &ct) const
 279 {
 280         const LambertConic2 *other = dynamic_cast<const LambertConic2*>(&ct);
 281         return (other != 0 && _lc1 == other->_lc1);
 282 }