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