src/transversemercator.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 "rd.h"
  46 #include "ellipsoid.h"
  47 #include "transversemercator.h"
  48
  49
  50 #define SPHSN(lat) \
  51         ((double)(_e.radius() / sqrt(1.e0 - _es * pow(sin(lat), 2))))
  52 #define SPHTMD(lat) \
  53         ((double)(_ap * lat - _bp * sin(2.e0 * lat) + _cp * sin(4.e0 * lat) \
  54           - _dp * sin(6.e0 * lat) + _ep * sin(8.e0 * lat)))
  55 #define DENOM(lat) \
  56         ((double)(sqrt(1.e0 - _es * pow(sin(lat),2))))
  57 #define SPHSR(lat) \
  58         ((double)(_e.radius() * (1.e0 - _es) / pow(DENOM(lat), 3)))
  59
  60
  61 TransverseMercator::TransverseMercator(const Ellipsoid &ellipsoid,
  62   double latitudeOrigin, double longitudeOrigin, double scale,
  63   double falseEasting, double falseNorthing)
  64 {
  65         double tn, tn2, tn3, tn4, tn5;
  66         double b;
  67
  68
  69         _e = ellipsoid;
  70         _longitudeOrigin = deg2rad(longitudeOrigin);
  71         _latitudeOrigin = deg2rad(latitudeOrigin);
  72         _scale = scale;
  73         _falseEasting = falseEasting;
  74         _falseNorthing = falseNorthing;
  75
  76         _es = 2 * _e.flattening() - _e.flattening() * _e.flattening();
  77         _ebs = (1 / (1 - _es)) - 1;
  78
  79         b = _e.radius() * (1 - _e.flattening());
  80
  81         tn = (_e.radius() - b) / (_e.radius() + b);
  82         tn2 = tn * tn;
  83         tn3 = tn2 * tn;
  84         tn4 = tn3 * tn;
  85         tn5 = tn4 * tn;
  86
  87         _ap = _e.radius() * (1.e0 - tn + 5.e0 * (tn2 - tn3) / 4.e0 + 81.e0
  88           * (tn4 - tn5) / 64.e0);
  89         _bp = 3.e0 * _e.radius() * (tn - tn2 + 7.e0 * (tn3 - tn4) / 8.e0 + 55.e0
  90           * tn5 / 64.e0 ) / 2.e0;
  91         _cp = 15.e0 * _e.radius() * (tn2 - tn3 + 3.e0 * (tn4 - tn5 ) / 4.e0) / 16.0;
  92         _dp = 35.e0 * _e.radius() * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
  93         _ep = 315.e0 * _e.radius() * (tn4 - tn5) / 512.e0;
  94 }
  95
  96 QPointF TransverseMercator::ll2xy(const Coordinates &c) const
  97 {
  98         double rl;
  99         double cl, c2, c3, c5, c7;
 100         double dlam;
 101         double eta, eta2, eta3, eta4;
 102         double sl, sn;
 103         double t, tan2, tan3, tan4, tan5, tan6;
 104         double t1, t2, t3, t4, t5, t6, t7, t8, t9;
 105         double tmd, tmdo;
 106         double x, y;
 107
 108
 109         dlam = deg2rad(c.lon()) - _longitudeOrigin;
 110
 111         if (dlam > M_PI)
 112                 dlam -= (2 * M_PI);
 113         if (dlam < -M_PI)
 114                 dlam += (2 * M_PI);
 115         if (fabs(dlam) < 2.e-10)
 116                 dlam = 0.0;
 117
 118         rl = deg2rad(c.lat());
 119         sl = sin(rl);
 120         cl = cos(rl);
 121         c2 = cl * cl;
 122         c3 = c2 * cl;
 123         c5 = c3 * c2;
 124         c7 = c5 * c2;
 125         t = sl / cl;
 126         tan2 = t * t;
 127         tan3 = tan2 * t;
 128         tan4 = tan3 * t;
 129         tan5 = tan4 * t;
 130         tan6 = tan5 * t;
 131         eta = _ebs * c2;
 132         eta2 = eta * eta;
 133         eta3 = eta2 * eta;
 134         eta4 = eta3 * eta;
 135
 136         sn = SPHSN(rl);
 137         tmd = SPHTMD(rl);
 138         tmdo = SPHTMD (_latitudeOrigin);
 139
 140
 141         t1 = (tmd - tmdo) * _scale;
 142         t2 = sn * sl * cl * _scale / 2.e0;
 143         t3 = sn * sl * c3 * _scale * (5.e0 - tan2 + 9.e0 * eta + 4.e0 * eta2)
 144           / 24.e0;
 145         t4 = sn * sl * c5 * _scale * (61.e0 - 58.e0 * tan2 + tan4 + 270.e0 * eta
 146           - 330.e0 * tan2 * eta + 445.e0 * eta2 + 324.e0 * eta3 - 680.e0 * tan2
 147           * eta2 + 88.e0 * eta4 - 600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4)
 148           / 720.e0;
 149         t5 = sn * sl * c7 * _scale * (1385.e0 - 3111.e0 * tan2 + 543.e0 * tan4
 150           - tan6) / 40320.e0;
 151
 152         y = _falseNorthing + t1 + pow(dlam, 2.e0) * t2 + pow(dlam, 4.e0) * t3
 153           + pow(dlam, 6.e0) * t4 + pow(dlam, 8.e0) * t5;
 154
 155
 156         t6 = sn * cl * _scale;
 157         t7 = sn * c3 * _scale * (1.e0 - tan2 + eta) /6.e0;
 158         t8 = sn * c5 * _scale * (5.e0 - 18.e0 * tan2 + tan4 + 14.e0 * eta - 58.e0
 159           * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3 - 64.e0 * tan2 * eta2 - 24.e0
 160           * tan2 * eta3) / 120.e0;
 161         t9 = sn * c7 * _scale * (61.e0 - 479.e0 * tan2 + 179.e0 * tan4 - tan6)
 162           / 5040.e0;
 163
 164         x = _falseEasting + dlam * t6 + pow(dlam, 3.e0) * t7 + pow(dlam, 5.e0)
 165           * t8 + pow(dlam, 7.e0) * t9;
 166
 167         return QPointF(x, y);
 168 }
 169
 170 Coordinates TransverseMercator::xy2ll(const QPointF &p) const
 171 {
 172         double cl;
 173         double de;
 174         double dlam;
 175         double eta, eta2, eta3, eta4;
 176         double ftphi;
 177         double sn;
 178         double sr;
 179         double t, tan2, tan4;
 180         double t10, t11, t12, t13, t14, t15, t16, t17;
 181         double tmd, tmdo;
 182         double lat, lon;
 183
 184
 185         tmdo = SPHTMD(_latitudeOrigin);
 186         tmd = tmdo + (p.y() - _falseNorthing) / _scale;
 187
 188         sr = SPHSR(0.e0);
 189         ftphi = tmd / sr;
 190
 191         for (int i = 0; i < 5 ; i++) {
 192                 t10 = SPHTMD(ftphi);
 193                 sr = SPHSR(ftphi);
 194                 ftphi = ftphi + (tmd - t10) / sr;
 195         }
 196
 197         sr = SPHSR(ftphi);
 198         sn = SPHSN(ftphi);
 199
 200     cl = cos(ftphi);
 201
 202     t = tan(ftphi);
 203     tan2 = t * t;
 204     tan4 = tan2 * tan2;
 205     eta = _ebs * pow(cl, 2);
 206     eta2 = eta * eta;
 207     eta3 = eta2 * eta;
 208     eta4 = eta3 * eta;
 209     de = p.x() - _falseEasting;
 210         if (fabs(de) < 0.0001)
 211                 de = 0.0;
 212
 213         t10 = t / (2.e0 * sr * sn * pow(_scale, 2));
 214         t11 = t * (5.e0  + 3.e0 * tan2 + eta - 4.e0 * pow(eta, 2) - 9.e0 * tan2
 215           * eta) / (24.e0 * sr * pow(sn, 3) * pow(_scale, 4));
 216         t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4 - 252.e0 * tan2
 217           * eta - 3.e0 * eta2 + 100.e0 * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
 218           * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2 + 84.e0 * tan2 * eta3 - 192.e0
 219           * tan2 * eta4) / (720.e0 * sr * pow(sn, 5) * pow(_scale, 6));
 220         t13 = t * (1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0 * pow(t,6))
 221           / (40320.e0 * sr * pow(sn, 7) * pow(_scale, 8));
 222         lat = ftphi - pow(de, 2) * t10 + pow(de, 4) * t11 - pow(de, 6) * t12
 223           + pow(de, 8) * t13;
 224
 225         t14 = 1.e0 / (sn * cl * _scale);
 226         t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn, 3) * cl * pow(_scale, 3));
 227         t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2 + 8.e0 * tan2 * eta
 228           + 24.e0 * tan4 - 4.e0 * eta3 + 4.e0 * tan2 * eta2 + 24.e0 * tan2 * eta3)
 229           / (120.e0 * pow(sn, 5) * cl * pow(_scale, 5));
 230         t17 = (61.e0 +  662.e0 * tan2 + 1320.e0 * tan4 + 720.e0 * pow(t,6))
 231           / (5040.e0 * pow(sn, 7) * cl * pow(_scale, 7));
 232
 233         dlam = de * t14 - pow(de, 3) * t15 + pow(de, 5) * t16 - pow(de, 7) * t17;
 234
 235         lon = _longitudeOrigin + dlam;
 236         while (lat > deg2rad(90.0)) {
 237                 lat = M_PI - lat;
 238                 lon += M_PI;
 239                 if (lon > M_PI)
 240                         lon -= (2 * M_PI);
 241         }
 242
 243         while (lat < deg2rad(-90.0)) {
 244                 lat = - (lat + M_PI);
 245                 lon += M_PI;
 246                 if (lon > M_PI)
 247                         lon -= (2 * M_PI);
 248         }
 249
 250         if (lon > (2 * M_PI))
 251                 lon -= (2 * M_PI);
 252         if (lon < -M_PI)
 253                 lon += (2 * M_PI);
 254
 255         return Coordinates(rad2deg(lon), rad2deg(lat));
 256 }