OmniSciDB  471d68cefb
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Utm.h
Go to the documentation of this file.
1 /*
2  * Copyright 2021 OmniSci, Inc.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
23 #include "Shared/funcannotations.h"
24 #include "Shared/misc.h"
25 
26 #include <cmath>
27 #include <limits>
28 
29 #ifdef _WIN32
30 // M_PI is defined a file then isn't include in the OS release
31 // TODO look at moving 'constants.h' to be OS avaliable.
32 constexpr double M_PI{3.14159265358979323846};
33 #endif
34 
35 namespace {
36 // Naming conventions break from style guide to match equation variables.
37 constexpr double f = 1 / 298.257223563; // WGS84 Earth flattening
38 constexpr double a = 6378137; // WGS84 Equatorial radius (m)
39 constexpr double k0 = 0.9996; // Point scale on central UTM meridians
40 constexpr double E0 = 500e3; // UTM False easting (m)
41 constexpr double deg_div_rad = 180 / M_PI;
42 constexpr double rad_div_deg = M_PI / 180;
43 
44 // Formulas are from https://arxiv.org/pdf/1002.1417v3.pdf
45 // except delta coefficients are from Kawase 2011.
46 constexpr double n = f / (2 - f);
47 // Section 3 Exact Mapping calculation using elliptic integral in Mathematica:
48 // f = 10^9 / 298257223563;
49 // n = f / (2 - f);
50 // N[(2*6378137/Pi) EllipticE[4n/(1+n)^2], 20]
51 // Can be verified with online approximation https://tinyurl.com/exact-elliptic
52 // constexpr double A = 6.3674491458234153093e6;
53 constexpr double A =
54  a * shared::horner(n * n, 1, 1. / 4, 1. / 64, 1. / 256, 25. / 16384) / (1 + n);
55 constexpr double k0_A = k0 * A;
56 
57 // clang-format off
58 #ifdef __CUDACC__
59 __device__ __constant__ const
60 #else
61 constexpr
62 #endif
63 std::array<double, 9> alphas // Eqs (35)
64  { std::numeric_limits<double>::quiet_NaN()
65  , shared::horner(n, 0, 1./2, -2./3, 5./16, 41./180, -127./288, 7891./37800, 72161./387072, -18975107./50803200)
66  , shared::horner(n, 0, 0, 13./48, -3./5, 557./1440, 281./630, -1983433./1935360, 13769./28800, 148003883./174182400)
67  , shared::horner(n, 0, 0, 0, 61./240, -103./140, 15061./26880, 167603./181440, -67102379./29030400, 79682431./79833600)
68  , shared::horner(n, 0, 0, 0, 0, 49561./161280, -179./168, 6601661./7257600, 97445./49896, -40176129013./7664025600)
69  , shared::horner(n, 0, 0, 0, 0, 0, 34729./80640, -3418889./1995840, 14644087./9123840, 2605413599./622702080)
70  , shared::horner(n, 0, 0, 0, 0, 0, 0, 212378941./319334400, -30705481./10378368, 175214326799./58118860800)
71  , shared::horner(n, 0, 0, 0, 0, 0, 0, 0, 1522256789./1383782400, -16759934899./3113510400)
72  , shared::horner(n, 0, 0, 0, 0, 0, 0, 0, 0, 1424729850961./743921418240)
73  };
74 #ifdef __CUDACC__
75 __device__ __constant__ const
76 #else
77 constexpr
78 #endif
79 std::array<double, 9> betas // Eqs (36)
80  { std::numeric_limits<double>::quiet_NaN()
81  , shared::horner(n, 0, 1./2, -2./3, 37./96, -1./360, -81./512, 96199./604800, -5406467./38707200, 7944359./67737600)
82  , shared::horner(n, 0, 0, 1./48, 1./15, -437./1440, 46./105, -1118711./3870720, 51841./1209600, 24749483./348364800)
83  , shared::horner(n, 0, 0, 0, 17./480, -37./840, -209./4480, 5569./90720, 9261899./58060800, -6457463./17740800)
84  , shared::horner(n, 0, 0, 0, 0, 4397./161280, -11./504, -830251./7257600, 466511./2494800, 324154477./7664025600)
85  , shared::horner(n, 0, 0, 0, 0, 0, 4583./161280, -108847./3991680, -8005831./63866880, 22894433./124540416)
86  , shared::horner(n, 0, 0, 0, 0, 0, 0, 20648693./638668800, -16363163./518918400, -2204645983./12915302400)
87  , shared::horner(n, 0, 0, 0, 0, 0, 0, 0, 219941297./5535129600, -497323811./12454041600)
88  , shared::horner(n, 0, 0, 0, 0, 0, 0, 0, 0, 191773887257./3719607091200)
89  };
90 /* Based on A General Formula for Calculating Meridian Arc Length and its Application
91 to Coordinate Conversion in the Gauss-Kr├╝ger Projection by Kawase Dec 2011
92 https://www.gsi.go.jp/common/000062452.pdf page 8
93 Mathematica to calculate Fourier coefficients:
94 ats[x_] := 2 Sqrt[n]/(n + 1)*ArcTanh[2 Sqrt[n]/(n + 1)*Sin[x]];
95 op[j_] := If[j == 0, #, op[j - 1][Cos[w] D[#, w]]] &;
96 term[k_, x_] := (op[k - 1][ats[w]^k Cos[w]] /. {w -> x})/k!;
97 sum[j_, x_] := Sum[term[k, x], {k, 1, j}]
98 InputForm@Expand@FullSimplify@Normal@Series[sum[8, x], {n, 0, 8}]
99 */
100 #ifdef __CUDACC__
101 __device__ __constant__ const
102 #else
103 constexpr
104 #endif
105 std::array<double, 9> deltas
106  { std::numeric_limits<double>::quiet_NaN()
107  , shared::horner(n, 0, 2, -2./3, -2, 116./45, 26./45, -2854./675, 16822./4725, 189416./99225)
108  , shared::horner(n, 0, 0, 7./3, -8./5, -227./45, 2704./315, 2323./945, -31256./1575, 141514./8505)
109  , shared::horner(n, 0, 0, 0, 56./15, -136./35, -1262./105, 73814./2835, 98738./14175, -2363828./31185)
110  , shared::horner(n, 0, 0, 0, 0, 4279./630, -332./35, - 399572./14175, 11763988./155925, 14416399./935550)
111  , shared::horner(n, 0, 0, 0, 0, 0, 4174./315, -144838./6237, -2046082./31185, 258316372./1216215)
112  , shared::horner(n, 0, 0, 0, 0, 0, 0, 601676./22275, -115444544./2027025, -2155215124./14189175)
113  , shared::horner(n, 0, 0, 0, 0, 0, 0, 0, 38341552./675675, -170079376./1216215)
114  , shared::horner(n, 0, 0, 0, 0, 0, 0, 0, 0, 1383243703./11351340)
115  };
116 // clang-format on
117 
118 constexpr unsigned N = 6; // N = 6 is smallest value that passes GeoSpatial.UTMTransform
119 } // namespace
120 
121 // The small_dlambda_ functions are used when x is within 4.77 degrees of the UTM zone's
122 // central meridian. This allows for a 1.77 degree overlap with the UTM zones on its
123 // east and west sides to use the faster trig/hypertrig functions.
125  // Constructor sets xi_ and eta_ shared by both _x and _y calculations.
126  unsigned const srid_;
127  double eta_;
128  double xi_;
129  bool small_dlambda_; // if true then use faster small-angle functions on dlambda
130 
131  // double sigma_;
132  // double tau_{0};
133  // double k_; // scaling
134  // static inline double sq(double const x) { return x * x; }
135 
136  public:
137  // Required: srid in [32601,32660] U [32701,32760]
138  ALWAYS_INLINE Transform4326ToUTM(unsigned const srid, double const x, double const y)
139  : srid_(srid) {
140  unsigned const zone = srid_ % 100u;
141  double const x0 = zone * 6.0 - 183;
142  double const dlambda = (x - x0) * rad_div_deg;
143  small_dlambda_ = -1.0 / 12 <= dlambda && dlambda <= 1.0 / 12;
144  double const phi = y * rad_div_deg;
145  double const c = 2 * sqrt(n) / (1 + n); // Boost will have a constexpr sqrt
146 
147  // https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Simplified_formulae
148  double const t = sinh(atanh(sin(phi)) - c * atanh(c * sin(phi)));
149  if (small_dlambda_) {
150  eta_ = shared::fastAtanh(shared::fastSin(dlambda) / sqrt(1 + t * t));
151  xi_ = atan(t / shared::fastCos(dlambda));
152  } else {
153  eta_ = atanh(sin(dlambda) / sqrt(1 + t * t));
154  xi_ = atan(t / cos(dlambda));
155  }
156  /* Calculate scaling. May be used in the future.
157  double sigma_sum = 0;
158  for (unsigned j = N; j; --j) {
159  sigma_sum += 2 * j * alphas[j] * cos(2 * j * xi_) * cosh(2 * j * eta_);
160  }
161  sigma_ = 1 + sigma_sum;
162  for (unsigned j = N; j; --j) {
163  tau_ += 2 * j * alphas[j] * sin(2 * j * xi_) * sinh(2 * j * eta_);
164  }
165  k_ = k0_A *
166  sqrt((1 + sq(tan(phi) * (1 - n) / (1 + n))) * (sq(sigma_) + sq(tau_)) /
167  (sq(t) + sq(cos(dlambda)))) /
168  a;
169  */
170  }
171 
172  ALWAYS_INLINE double calculateX() const {
173  double sum = 0; // Sum in reverse to add smallest numbers first
174  if (small_dlambda_) {
175  for (unsigned j = N; j; --j) {
176  sum += alphas[j] * cos(2 * j * xi_) * shared::fastSinh(2 * j * eta_);
177  }
178  } else {
179  for (unsigned j = N; j; --j) {
180  sum += alphas[j] * cos(2 * j * xi_) * sinh(2 * j * eta_);
181  }
182  }
183  return E0 + k0_A * (eta_ + sum);
184  }
185 
186  // The inline specifier is required here (using ALWAYS_INLINE for safety.)
187  // Without it, GeoSpatial.UTMTransform fails with
188  // C++ exception with description "Query must run in cpu mode." thrown in the test body.
189  // and the log contains:
190  // NativeCodegen.cpp:1167 Failed to generate PTX: NVVM IR ParseError: generatePTX: use
191  // of undefined value '@_ZNK18Transform4326ToUTM10calculateYEv'
192  // %117 = call fastcc double
193  // @_ZNK18Transform4326ToUTM10calculateYEv(%class.Transform4326ToUTM* nonnull
194  // dereferenceable(24) %0) #1
195  // ^
196  //. Switching to CPU execution target.
197  ALWAYS_INLINE double calculateY() const {
198  double const N0 = (32700 < srid_) * 10e6; // UTM False northing (m)
199  double sum = 0;
200  if (small_dlambda_) {
201  for (unsigned j = N; j; --j) {
202  sum += alphas[j] * sin(2 * j * xi_) * shared::fastCosh(2 * j * eta_);
203  }
204  } else {
205  for (unsigned j = N; j; --j) {
206  sum += alphas[j] * sin(2 * j * xi_) * cosh(2 * j * eta_);
207  }
208  }
209  return N0 + k0_A * (xi_ + sum);
210  }
211 };
212 
214  // Constructor sets eta_ and xi_ shared by both _x and _y calculations.
215  unsigned const srid_;
216  double eta_;
217  double xi_;
218  bool small_eta_; // if true then use faster small-angle functions on eta_
219 
220  // static inline double sq(double const x) { return x * x; }
221 
222  public:
223  // ALWAYS_INLINE is required here.
224  // inline didn't fix similar problem as with Transform4326ToUTM::calculateY() above.
225  // Required: srid in [32601,32660] U [32701,32760]
226  ALWAYS_INLINE TransformUTMTo4326(unsigned const srid, double const x, double const y)
227  : srid_(srid) {
228  double const eta = (x - E0) / k0_A;
229  small_eta_ = -1.0 / 12 <= eta && eta <= 1.0 / 12;
230  double const N0 = (32700 < srid_) * 10e6; // UTM False northing (m)
231  double const xi = (y - N0) / k0_A;
232 
233  double eta_sum = 0;
234  double xi_sum = 0;
235  if (small_eta_) {
236  for (unsigned j = N; j; --j) {
237  eta_sum += betas[j] * cos(2 * j * xi) * shared::fastSinh(2 * j * eta);
238  }
239  for (unsigned j = N; j; --j) {
240  xi_sum += betas[j] * sin(2 * j * xi) * shared::fastCosh(2 * j * eta);
241  }
242  } else {
243  for (unsigned j = N; j; --j) {
244  eta_sum += betas[j] * cos(2 * j * xi) * sinh(2 * j * eta);
245  }
246  for (unsigned j = N; j; --j) {
247  xi_sum += betas[j] * sin(2 * j * xi) * cosh(2 * j * eta);
248  }
249  }
250  eta_ = eta - eta_sum;
251  xi_ = xi - xi_sum;
252  }
253 
254  // lambda (longitude, degrees)
255  ALWAYS_INLINE double calculateX() const {
256  unsigned const zone = srid_ % 100u;
257  double const lambda0 = zone * 6.0 - 183;
258  double const sinh_eta = small_eta_ ? shared::fastSinh(eta_) : sinh(eta_);
259  return lambda0 + atan(sinh_eta / cos(xi_)) * deg_div_rad;
260  }
261 
262  // phi (latitude, degrees)
263  ALWAYS_INLINE double calculateY() const {
264 #if 1
265  // https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Simplified_formulae
266  double const cosh_eta = small_eta_ ? shared::fastCosh(eta_) : cosh(eta_);
267  double const chi = asin(sin(xi_) / cosh_eta);
268  double sum = 0;
269  for (unsigned j = N; j; --j) {
270  sum += deltas[j] * sin(2 * j * chi);
271  }
272  return (chi + sum) * deg_div_rad;
273 #else
274  // Heavier calculation from Transverse Mercator with an accuracy of a few nanometers
275  // by Karney 3 Feb 2011. This does not make use of the constexpr delta Fourier
276  // coefficients used by Kawase 2011 above which appears to be more accurate.
277  double const e = std::sqrt(f * (2 - f));
278  double const taup =
279  std::sin(xi_) / std::sqrt(sq(std::sinh(eta_)) + sq(std::cos(xi_)));
280 
281  double tau[2]{taup, 0.0 / 0.0};
282  unsigned i = 0;
283  for (; tau[0] != tau[1]; ++i) {
284  double const tau_i = tau[i & 1];
285  double const sigma_i =
286  std::sinh(e * std::tanh(e * tau_i / std::sqrt(1 + sq(tau_i))));
287  double const taup_i =
288  tau_i * std::sqrt(1 + sq(sigma_i)) - sigma_i * std::sqrt(1 + sq(tau_i));
289  double const dtau_i = (taup - taup_i) * (1 + (1 - e * e) * sq(tau_i)) /
290  ((1 - e * e) * std::sqrt((1 + sq(taup_i)) * (1 + sq(tau_i))));
291  tau[~i & 1] = tau_i + dtau_i;
292  }
293  return std::atan(tau[0]) * deg_div_rad;
294 #endif
295  }
296 };
constexpr double A
Definition: Utm.h:53
double eta_
Definition: Utm.h:216
double fastCosh(double const x)
Definition: misc.h:223
constexpr double E0
Definition: Utm.h:40
constexpr std::array< double, 9 > alphas
Definition: Utm.h:64
constexpr std::array< double, 9 > deltas
Definition: Utm.h:106
constexpr double k0_A
Definition: Utm.h:55
double fastCos(double const x)
Definition: misc.h:214
ALWAYS_INLINE double calculateY() const
Definition: Utm.h:197
ALWAYS_INLINE TransformUTMTo4326(unsigned const srid, double const x, double const y)
Definition: Utm.h:226
constexpr double rad_div_deg
Definition: Utm.h:42
#define M_PI
Definition: constants.h:25
constexpr double a
Definition: Utm.h:38
constexpr std::array< double, 9 > betas
Definition: Utm.h:80
ALWAYS_INLINE double calculateX() const
Definition: Utm.h:172
ALWAYS_INLINE Transform4326ToUTM(unsigned const srid, double const x, double const y)
Definition: Utm.h:138
ALWAYS_INLINE double calculateY() const
Definition: Utm.h:263
constexpr double k0
Definition: Utm.h:39
unsigned const srid_
Definition: Utm.h:126
double xi_
Definition: Utm.h:217
double xi_
Definition: Utm.h:128
double eta_
Definition: Utm.h:127
ALWAYS_INLINE double calculateX() const
Definition: Utm.h:255
constexpr double horner(double const x, double const c0, COEFFICIENTS...c)
Definition: misc.h:198
constexpr unsigned N
Definition: Utm.h:118
unsigned const srid_
Definition: Utm.h:215
char * t
double fastSinh(double const x)
Definition: misc.h:241
char * f
constexpr double n
Definition: Utm.h:46
bool small_dlambda_
Definition: Utm.h:129
#define ALWAYS_INLINE
constexpr double deg_div_rad
Definition: Utm.h:41
bool small_eta_
Definition: Utm.h:218
double fastAtanh(double const x)
Definition: misc.h:208
double fastSin(double const x)
Definition: misc.h:232