1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 // This file is manually converted from PROJ4
3
4 // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
5
6 // This file was modified by Oracle on 2017-2020.
7 // Modifications copyright (c) 2017-2020, Oracle and/or its affiliates.
8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
9
10 // Use, modification and distribution is subject to the Boost Software License,
11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
12 // http://www.boost.org/LICENSE_1_0.txt)
13
14 // This file is converted from PROJ4, http://trac.osgeo.org/proj
15 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
16 // PROJ4 is maintained by Frank Warmerdam
17 // PROJ4 is converted to Geometry Library by
18 // Barend Gehrels (Geodan, Amsterdam)
19 // Adam Wulkiewicz
20
21 // Original copyright notice:
22
23 // Permission is hereby granted, free of charge, to any person obtaining a
24 // copy of this software and associated documentation files (the "Software"),
25 // to deal in the Software without restriction, including without limitation
26 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
27 // and/or sell copies of the Software, and to permit persons to whom the
28 // Software is furnished to do so, subject to the following conditions:
29
30 // The above copyright notice and this permission notice shall be included
31 // in all copies or substantial portions of the Software.
32
33 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
34 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
35 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
36 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
37 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
38 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
39 // DEALINGS IN THE SOFTWARE.
40
41 #ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP
42 #define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP
43
44 #include <string>
45 #include <type_traits>
46 #include <vector>
47
48 #include <boost/geometry/core/static_assert.hpp>
49 #include <boost/geometry/formulas/eccentricity_sqr.hpp>
50 #include <boost/geometry/util/math.hpp>
51
52 #include <boost/geometry/srs/projections/dpar.hpp>
53 #include <boost/geometry/srs/projections/impl/pj_datum_set.hpp>
54 #include <boost/geometry/srs/projections/impl/pj_ellps.hpp>
55 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
56 #include <boost/geometry/srs/projections/proj4.hpp>
57 #include <boost/geometry/srs/projections/spar.hpp>
58
59
60 namespace boost { namespace geometry { namespace projections {
61
62 namespace detail {
63
64 /* set ellipsoid parameters a and es */
65 template <typename T>
SIXTH()66 inline T SIXTH() { return .1666666666666666667; } /* 1/6 */
67 template <typename T>
RA4()68 inline T RA4() { return .04722222222222222222; } /* 17/360 */
69 template <typename T>
RA6()70 inline T RA6() { return .02215608465608465608; } /* 67/3024 */
71 template <typename T>
RV4()72 inline T RV4() { return .06944444444444444444; } /* 5/72 */
73 template <typename T>
RV6()74 inline T RV6() { return .04243827160493827160; } /* 55/1296 */
75
76 template <typename T>
pj_ell_b_to_es(T const & a,T const & b)77 inline T pj_ell_b_to_es(T const& a, T const& b)
78 {
79 return 1. - (b * b) / (a * a);
80 }
81
82 /************************************************************************/
83 /* pj_ell_init_ellps() */
84 /************************************************************************/
85
86 // Originally a part of pj_ell_set()
87 template <typename T>
pj_ell_init_ellps(srs::detail::proj4_parameters const & params,T & a,T & b)88 inline bool pj_ell_init_ellps(srs::detail::proj4_parameters const& params, T &a, T &b)
89 {
90 /* check if ellps present and temporarily append its values to pl */
91 std::string name = pj_get_param_s(params, "ellps");
92 if (! name.empty())
93 {
94 const pj_ellps_type<T>* pj_ellps = pj_get_ellps<T>().first;
95 const int n = pj_get_ellps<T>().second;
96 int index = -1;
97 for (int i = 0; i < n && index == -1; i++)
98 {
99 if(pj_ellps[i].id == name)
100 {
101 index = i;
102 }
103 }
104
105 if (index == -1) {
106 BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) );
107 }
108
109 pj_ellps_type<T> const& pj_ellp = pj_ellps[index];
110 a = pj_ellp.a;
111 b = pj_ellp.b;
112
113 return true;
114 }
115
116 return false;
117 }
118
119 template <typename T>
pj_ell_init_ellps(srs::dpar::parameters<T> const & params,T & a,T & b)120 inline bool pj_ell_init_ellps(srs::dpar::parameters<T> const& params, T &a, T &b)
121 {
122 /* check if ellps present and temporarily append its values to pl */
123 typename srs::dpar::parameters<T>::const_iterator
124 it = pj_param_find(params, srs::dpar::ellps);
125 if (it != params.end())
126 {
127 if (it->template is_value_set<int>())
128 {
129 const pj_ellps_type<T>* pj_ellps = pj_get_ellps<T>().first;
130 const int n = pj_get_ellps<T>().second;
131 int i = it->template get_value<int>();
132
133 if (i < 0 || i >= n) {
134 BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) );
135 }
136
137 pj_ellps_type<T> const& pj_ellp = pj_ellps[i];
138 a = pj_ellp.a;
139 b = pj_ellp.b;
140 }
141 else if (it->template is_value_set<T>())
142 {
143 a = it->template get_value<T>();
144 b = a;
145 }
146 else if (it->template is_value_set<srs::spheroid<T> >())
147 {
148 srs::spheroid<T> const& s = it->template get_value<srs::spheroid<T> >();
149 a = geometry::get_radius<0>(s);
150 b = geometry::get_radius<2>(s);
151 }
152 else
153 {
154 BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) );
155 }
156
157 return true;
158 }
159
160 return false;
161 }
162
163 template
164 <
165 typename Params,
166 int I = geometry::tuples::find_index_if
167 <
168 Params,
169 srs::spar::detail::is_param_tr<srs::spar::detail::ellps_traits>::pred
170 >::value,
171 int N = geometry::tuples::size<Params>::value
172 >
173 struct pj_ell_init_ellps_static
174 {
175 template <typename T>
applyboost::geometry::projections::detail::pj_ell_init_ellps_static176 static bool apply(Params const& params, T &a, T &b)
177 {
178 typedef typename geometry::tuples::element<I, Params>::type param_type;
179 typedef srs::spar::detail::ellps_traits<param_type> traits_type;
180 typedef typename traits_type::template model_type<T>::type model_type;
181
182 param_type const& param = geometry::tuples::get<I>(params);
183 model_type const& model = traits_type::template model<T>(param);
184
185 a = geometry::get_radius<0>(model);
186 b = geometry::get_radius<2>(model);
187
188 return true;
189 }
190 };
191 template <typename Params, int N>
192 struct pj_ell_init_ellps_static<Params, N, N>
193 {
194 template <typename T>
applyboost::geometry::projections::detail::pj_ell_init_ellps_static195 static bool apply(Params const& , T & , T & )
196 {
197 return false;
198 }
199 };
200
201 template <typename T, typename ...Ps>
pj_ell_init_ellps(srs::spar::parameters<Ps...> const & params,T & a,T & b)202 inline bool pj_ell_init_ellps(srs::spar::parameters<Ps...> const& params,
203 T &a, T &b)
204 {
205 return pj_ell_init_ellps_static
206 <
207 srs::spar::parameters<Ps...>
208 >::apply(params, a, b);
209 }
210
211 /************************************************************************/
212 /* pj_ell_init() */
213 /************************************************************************/
214
215 /* initialize geographic shape parameters */
216 // This function works differently than the original pj_ell_set().
217 // It doesn't push parameters defined in ellps into params list.
218 // Instead it tries to use size (a, R) and shape (es, e, rf, f, b) parameters
219 // and then if needed falls back to ellps, then to datum and then to the default WGS84
220 template <typename Params, typename T>
pj_ell_init(Params const & params,T & a,T & es)221 inline void pj_ell_init(Params const& params, T &a, T &es)
222 {
223 /* check for varying forms of ellipsoid input */
224 a = es = 0.;
225
226 /* R takes precedence */
227 if (pj_param_f<srs::spar::r>(params, "R", srs::dpar::r, a)) {
228 /* empty */
229 } else { /* probable elliptical figure */
230
231 // Set ellipsoid's size parameter
232 a = pj_get_param_f<T, srs::spar::a>(params, "a", srs::dpar::a);
233 bool is_a_set = a != 0.0;
234
235 // Set ellipsoid's shape parameter
236 T b = 0.0;
237 bool is_ell_set = false;
238 if (pj_param_f<srs::spar::es>(params, "es", srs::dpar::es, es)) {/* eccentricity squared */
239 /* empty */
240 is_ell_set = true;
241 } else if (pj_param_f<srs::spar::e>(params, "e", srs::dpar::e, es)) { /* eccentricity */
242 es = es * es;
243 is_ell_set = true;
244 } else if (pj_param_f<srs::spar::rf>(params, "rf", srs::dpar::rf, es)) { /* recip flattening */
245 if (es == 0.0) {
246 BOOST_THROW_EXCEPTION( projection_exception(error_rev_flattening_is_zero) );
247 }
248 es = 1./ es;
249 es = es * (2. - es);
250 is_ell_set = true;
251 } else if (pj_param_f<srs::spar::f>(params, "f", srs::dpar::f, es)) { /* flattening */
252 es = es * (2. - es);
253 is_ell_set = true;
254 } else if (pj_param_f<srs::spar::b>(params, "b", srs::dpar::b, b)) { /* minor axis */
255 es = pj_ell_b_to_es(a, b);
256 is_ell_set = true;
257 } /* else es == 0. and sphere of radius a */
258
259 // NOTE: Below when ellps is used to initialize a and es
260 // b is not set because it only has sense together with a
261 // but a could have been set separately before, e.g. consider passing:
262 // a=1 ellps=airy (a=6377563.396 b=6356256.910)
263 // after setting size parameter a and shape parameter from ellps
264 // b has to be recalculated
265
266 // If ellipsoid's parameters are not set directly
267 // use ellps parameter
268 if (! is_a_set || ! is_ell_set) {
269 T ellps_a = 0, ellps_b = 0;
270 if (pj_ell_init_ellps(params, ellps_a, ellps_b)) {
271 if (! is_a_set) {
272 a = ellps_a;
273 is_a_set = true;
274 }
275 if (! is_ell_set) {
276 es = pj_ell_b_to_es(ellps_a, ellps_b);
277 is_ell_set = true;
278 }
279 }
280 }
281
282 // If ellipsoid's parameters are not set
283 // use ellps defined by datum parameter
284 if (! is_a_set || ! is_ell_set)
285 {
286 const pj_datums_type<T>* datum = pj_datum_find_datum<T>(params);
287 if (datum != NULL)
288 {
289 pj_ellps_type<T> const& pj_ellp = pj_get_ellps<T>().first[datum->ellps];
290 if (! is_a_set) {
291 a = pj_ellp.a;
292 is_a_set = true;
293 }
294 if (! is_ell_set) {
295 es = pj_ell_b_to_es(pj_ellp.a, pj_ellp.b);
296 is_ell_set = true;
297 }
298 }
299 }
300
301 // If ellipsoid's parameters are still not set
302 // use default WGS84
303 if ((! is_a_set || ! is_ell_set)
304 && ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs))
305 {
306 pj_ellps_type<T> const& pj_ellp = pj_get_ellps<T>().first[srs::dpar::ellps_wgs84];
307 if (! is_a_set) {
308 a = pj_ellp.a;
309 is_a_set = true;
310 }
311 if (! is_ell_set) {
312 es = pj_ell_b_to_es(pj_ellp.a, pj_ellp.b);
313 is_ell_set = true;
314 }
315 }
316
317 if (b == 0.0)
318 b = a * sqrt(1. - es);
319
320 /* following options turn ellipsoid into equivalent sphere */
321 if (pj_get_param_b<srs::spar::r_au>(params, "R_A", srs::dpar::r_au)) { /* sphere--area of ellipsoid */
322 a *= 1. - es * (SIXTH<T>() + es * (RA4<T>() + es * RA6<T>()));
323 es = 0.;
324 } else if (pj_get_param_b<srs::spar::r_v>(params, "R_V", srs::dpar::r_v)) { /* sphere--vol. of ellipsoid */
325 a *= 1. - es * (SIXTH<T>() + es * (RV4<T>() + es * RV6<T>()));
326 es = 0.;
327 } else if (pj_get_param_b<srs::spar::r_a>(params, "R_a", srs::dpar::r_a)) { /* sphere--arithmetic mean */
328 a = .5 * (a + b);
329 es = 0.;
330 } else if (pj_get_param_b<srs::spar::r_g>(params, "R_g", srs::dpar::r_g)) { /* sphere--geometric mean */
331 a = sqrt(a * b);
332 es = 0.;
333 } else if (pj_get_param_b<srs::spar::r_h>(params, "R_h", srs::dpar::r_h)) { /* sphere--harmonic mean */
334 a = 2. * a * b / (a + b);
335 es = 0.;
336 } else {
337 T tmp;
338 bool i = pj_param_r<srs::spar::r_lat_a>(params, "R_lat_a", srs::dpar::r_lat_a, tmp);
339 if (i || /* sphere--arith. */
340 pj_param_r<srs::spar::r_lat_g>(params, "R_lat_g", srs::dpar::r_lat_g, tmp)) { /* or geom. mean at latitude */
341
342 tmp = sin(tmp);
343 if (geometry::math::abs(tmp) > geometry::math::half_pi<T>()) {
344 BOOST_THROW_EXCEPTION( projection_exception(error_ref_rad_larger_than_90) );
345 }
346 tmp = 1. - es * tmp * tmp;
347 a *= i ? .5 * (1. - es + tmp) / ( tmp * sqrt(tmp)) :
348 sqrt(1. - es) / tmp;
349 es = 0.;
350 }
351 }
352 }
353
354 /* some remaining checks */
355 if (es < 0.) {
356 BOOST_THROW_EXCEPTION( projection_exception(error_es_less_than_zero) );
357 }
358 if (a <= 0.) {
359 BOOST_THROW_EXCEPTION( projection_exception(error_major_axis_not_given) );
360 }
361 }
362
363 template <typename Params>
364 struct static_srs_tag_check_nonexpanded
365 {
366 typedef std::conditional_t
367 <
368 geometry::tuples::exists_if
369 <
370 Params, srs::spar::detail::is_param_t<srs::spar::r>::pred
371 >::value
372 || geometry::tuples::exists_if
373 <
374 Params, srs::spar::detail::is_param<srs::spar::r_au>::pred
375 >::value
376 || geometry::tuples::exists_if
377 <
378 Params, srs::spar::detail::is_param<srs::spar::r_v>::pred
379 >::value
380 || geometry::tuples::exists_if
381 <
382 Params, srs::spar::detail::is_param<srs::spar::r_a>::pred
383 >::value
384 || geometry::tuples::exists_if
385 <
386 Params, srs::spar::detail::is_param<srs::spar::r_g>::pred
387 >::value
388 || geometry::tuples::exists_if
389 <
390 Params, srs::spar::detail::is_param<srs::spar::r_h>::pred
391 >::value
392 || geometry::tuples::exists_if
393 <
394 Params, srs::spar::detail::is_param_t<srs::spar::r_lat_a>::pred
395 >::value
396 || geometry::tuples::exists_if
397 <
398 Params, srs::spar::detail::is_param_t<srs::spar::r_lat_g>::pred
399 >::value,
400 srs_sphere_tag,
401 // NOTE: The assumption here is that if the user defines either one of:
402 // b, es, e, f, rf parameters then he wants to define spheroid, not sphere
403 std::conditional_t
404 <
405 geometry::tuples::exists_if
406 <
407 Params, srs::spar::detail::is_param_t<srs::spar::b>::pred
408 >::value
409 || geometry::tuples::exists_if
410 <
411 Params, srs::spar::detail::is_param_t<srs::spar::es>::pred
412 >::value
413 || geometry::tuples::exists_if
414 <
415 Params, srs::spar::detail::is_param_t<srs::spar::e>::pred
416 >::value
417 || geometry::tuples::exists_if
418 <
419 Params, srs::spar::detail::is_param_t<srs::spar::rf>::pred
420 >::value
421 || geometry::tuples::exists_if
422 <
423 Params, srs::spar::detail::is_param_t<srs::spar::f>::pred
424 >::value,
425 srs_spheroid_tag,
426 void
427 >
428 > type;
429 };
430
431 template <typename Params>
432 struct static_srs_tag_check_ellps
433 {
434 typedef typename geometry::tag
435 <
436 typename srs::spar::detail::ellps_traits
437 <
438 typename geometry::tuples::find_if
439 <
440 Params,
441 srs::spar::detail::is_param_tr<srs::spar::detail::ellps_traits>::pred
442 >::type
443 >::template model_type<double>::type // dummy type
444 >::type type;
445 };
446
447 template <typename Params>
448 struct static_srs_tag_check_datum
449 {
450 typedef typename geometry::tag
451 <
452 typename srs::spar::detail::ellps_traits
453 <
454 typename srs::spar::detail::datum_traits
455 <
456 typename geometry::tuples::find_if
457 <
458 Params,
459 srs::spar::detail::is_param_tr<srs::spar::detail::datum_traits>::pred
460 >::type
461 >::ellps_type
462 >::template model_type<double>::type // dummy type
463 >::type type;
464 };
465
466 template
467 <
468 typename Params,
469 typename NonExpandedTag = typename static_srs_tag_check_nonexpanded
470 <
471 Params
472 >::type,
473 typename EllpsTag = typename static_srs_tag_check_ellps
474 <
475 Params
476 >::type,
477 typename DatumTag = typename static_srs_tag_check_datum
478 <
479 Params
480 >::type
481 >
482 struct static_srs_tag
483 {
484 // User passed one of the non-ellps, non-datum parameters
485 typedef NonExpandedTag type;
486 };
487
488 template <typename Params, typename EllpsTag, typename DatumTag>
489 struct static_srs_tag<Params, void, EllpsTag, DatumTag>
490 {
491 // User didn't pass neither one of the non-ellps, non-datum parameters
492 // but passed ellps
493 typedef EllpsTag type;
494 };
495
496 template <typename Params, typename DatumTag>
497 struct static_srs_tag<Params, void, void, DatumTag>
498 {
499 // User didn't pass neither one of the non-ellps, non-datum parameters
500 // nor ellps parameter but passed datum parameter
501 typedef DatumTag type;
502 };
503
504 template <typename Params>
505 struct static_srs_tag<Params, void, void, void>
506 {
507 // User didn't pass any parameter defining model
508 // so use default or generate error
509 typedef std::conditional_t
510 <
511 geometry::tuples::exists_if
512 <
513 Params, srs::spar::detail::is_param<srs::spar::no_defs>::pred
514 >::value,
515 void,
516 srs_spheroid_tag // WGS84
517 > type;
518
519 static const bool is_found = ! std::is_void<type>::value;
520 BOOST_GEOMETRY_STATIC_ASSERT((is_found),
521 "Expected ellipsoid or sphere definition.",
522 Params);
523 };
524
525
526 template <typename T>
pj_calc_ellipsoid_params(parameters<T> & p,T const & a,T const & es)527 inline void pj_calc_ellipsoid_params(parameters<T> & p, T const& a, T const& es) {
528 /****************************************************************************************
529 Calculate a large number of ancillary ellipsoidal parameters, in addition to
530 the two traditional PROJ defining parameters: Semimajor axis, a, and the
531 eccentricity squared, es.
532
533 Most of these parameters are fairly cheap to compute in comparison to the overall
534 effort involved in initializing a PJ object. They may, however, take a substantial
535 part of the time taken in computing an individual point transformation.
536
537 So by providing them up front, we can amortize the (already modest) cost over all
538 transformations carried out over the entire lifetime of a PJ object, rather than
539 incur that cost for every single transformation.
540
541 Most of the parameter calculations here are based on the "angular eccentricity",
542 i.e. the angle, measured from the semiminor axis, of a line going from the north
543 pole to one of the foci of the ellipsoid - or in other words: The arc sine of the
544 eccentricity.
545
546 The formulae used are mostly taken from:
547
548 Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
549 Columbus, Ohio: Dept. of Geodetic Science
550 and Surveying, Ohio State University.
551
552 ****************************************************************************************/
553
554 p.a = a;
555 p.es = es;
556
557 /* Compute some ancillary ellipsoidal parameters */
558 if (p.e==0)
559 p.e = sqrt(p.es); /* eccentricity */
560 //p.alpha = asin (p.e); /* angular eccentricity */
561
562 /* second eccentricity */
563 //p.e2 = tan (p.alpha);
564 //p.e2s = p.e2 * p.e2;
565
566 /* third eccentricity */
567 //p.e3 = (0!=p.alpha)? sin (p.alpha) / sqrt(2 - sin (p.alpha)*sin (p.alpha)): 0;
568 //p.e3s = p.e3 * p.e3;
569
570 /* flattening */
571 //if (0==p.f)
572 // p.f = 1 - cos (p.alpha); /* = 1 - sqrt (1 - PIN->es); */
573 //p.rf = p.f != 0.0 ? 1.0/p.f: HUGE_VAL;
574
575 /* second flattening */
576 //p.f2 = (cos(p.alpha)!=0)? 1/cos (p.alpha) - 1: 0;
577 //p.rf2 = p.f2 != 0.0 ? 1/p.f2: HUGE_VAL;
578
579 /* third flattening */
580 //p.n = pow (tan (p.alpha/2), 2);
581 //p.rn = p.n != 0.0 ? 1/p.n: HUGE_VAL;
582
583 /* ...and a few more */
584 //if (0==p.b)
585 // p.b = (1 - p.f)*p.a;
586 //p.rb = 1. / p.b;
587 p.ra = 1. / p.a;
588
589 p.one_es = 1. - p.es;
590 if (p.one_es == 0.) {
591 BOOST_THROW_EXCEPTION( projection_exception(error_eccentricity_is_one) );
592 }
593
594 p.rone_es = 1./p.one_es;
595 }
596
597
598 } // namespace detail
599 }}} // namespace boost::geometry::projections
600
601 #endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP
602