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