1 // Boost.Geometry
2 // This file is manually converted from PROJ4
3
4 // This file was modified by Oracle on 2017, 2018.
5 // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
6 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
7
8 // Use, modification and distribution is subject to the Boost Software License,
9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
10 // http://www.boost.org/LICENSE_1_0.txt)
11
12 // This file is converted from PROJ4, http://trac.osgeo.org/proj
13 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
14 // PROJ4 is maintained by Frank Warmerdam
15 // This file was converted to Geometry Library by Adam Wulkiewicz
16
17 // Original copyright notice:
18
19 // Copyright (c) 2000, Frank Warmerdam
20
21 // Permission is hereby granted, free of charge, to any person obtaining a
22 // copy of this software and associated documentation files (the "Software"),
23 // to deal in the Software without restriction, including without limitation
24 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
25 // and/or sell copies of the Software, and to permit persons to whom the
26 // Software is furnished to do so, subject to the following conditions:
27
28 // The above copyright notice and this permission notice shall be included
29 // in all copies or substantial portions of the Software.
30
31 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
32 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
33 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
34 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
35 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
36 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
37 // DEALINGS IN THE SOFTWARE.
38
39 #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
40 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
41
42
43 #include <boost/geometry/core/access.hpp>
44 #include <boost/geometry/core/coordinate_dimension.hpp>
45 #include <boost/geometry/core/radian_access.hpp>
46
47 #include <boost/geometry/srs/projections/impl/geocent.hpp>
48 #include <boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp>
49 #include <boost/geometry/srs/projections/impl/projects.hpp>
50 #include <boost/geometry/srs/projections/invalid_point.hpp>
51
52 #include <boost/geometry/util/range.hpp>
53
54 #include <cstring>
55 #include <cmath>
56
57
58 namespace boost { namespace geometry { namespace projections
59 {
60
61 namespace detail
62 {
63
64 // -----------------------------------------------------------
65 // Boost.Geometry helpers begin
66 // -----------------------------------------------------------
67
68 template
69 <
70 typename Point,
71 bool HasCoord2 = (dimension<Point>::value > 2)
72 >
73 struct z_access
74 {
75 typedef typename coordinate_type<Point>::type type;
getboost::geometry::projections::detail::z_access76 static inline type get(Point const& point)
77 {
78 return geometry::get<2>(point);
79 }
setboost::geometry::projections::detail::z_access80 static inline void set(Point & point, type const& h)
81 {
82 return geometry::set<2>(point, h);
83 }
84 };
85
86 template <typename Point>
87 struct z_access<Point, false>
88 {
89 typedef typename coordinate_type<Point>::type type;
getboost::geometry::projections::detail::z_access90 static inline type get(Point const& )
91 {
92 return type(0);
93 }
setboost::geometry::projections::detail::z_access94 static inline void set(Point & , type const& )
95 {}
96 };
97
98 template <typename XYorXYZ>
99 inline typename z_access<XYorXYZ>::type
get_z(XYorXYZ const & xy_or_xyz)100 get_z(XYorXYZ const& xy_or_xyz)
101 {
102 return z_access<XYorXYZ>::get(xy_or_xyz);
103 }
104
105 template <typename XYorXYZ>
set_z(XYorXYZ & xy_or_xyz,typename z_access<XYorXYZ>::type const & z)106 inline void set_z(XYorXYZ & xy_or_xyz,
107 typename z_access<XYorXYZ>::type const& z)
108 {
109 return z_access<XYorXYZ>::set(xy_or_xyz, z);
110 }
111
112 template
113 <
114 typename Range,
115 bool AddZ = (dimension<typename boost::range_value<Range>::type>::value < 3)
116 >
117 struct range_wrapper
118 {
119 typedef Range range_type;
120 typedef typename boost::range_value<Range>::type point_type;
121 typedef typename coordinate_type<point_type>::type coord_t;
122
range_wrapperboost::geometry::projections::detail::range_wrapper123 range_wrapper(Range & range)
124 : m_range(range)
125 {}
126
get_rangeboost::geometry::projections::detail::range_wrapper127 range_type & get_range() { return m_range; }
128
get_zboost::geometry::projections::detail::range_wrapper129 coord_t get_z(std::size_t i) { return detail::get_z(range::at(m_range, i)); }
set_zboost::geometry::projections::detail::range_wrapper130 void set_z(std::size_t i, coord_t const& v) { return detail::set_z(range::at(m_range, i), v); }
131
132 private:
133 Range & m_range;
134 };
135
136 template <typename Range>
137 struct range_wrapper<Range, true>
138 {
139 typedef Range range_type;
140 typedef typename boost::range_value<Range>::type point_type;
141 typedef typename coordinate_type<point_type>::type coord_t;
142
range_wrapperboost::geometry::projections::detail::range_wrapper143 range_wrapper(Range & range)
144 : m_range(range)
145 , m_zs(boost::size(range), coord_t(0))
146 {}
147
get_rangeboost::geometry::projections::detail::range_wrapper148 range_type & get_range() { return m_range; }
149
get_zboost::geometry::projections::detail::range_wrapper150 coord_t get_z(std::size_t i) { return m_zs[i]; }
set_zboost::geometry::projections::detail::range_wrapper151 void set_z(std::size_t i, coord_t const& v) { m_zs[i] = v; }
152
153 private:
154 Range & m_range;
155 std::vector<coord_t> m_zs;
156 };
157
158 // -----------------------------------------------------------
159 // Boost.Geometry helpers end
160 // -----------------------------------------------------------
161
162 template <typename Par>
Dx_BF(Par const & defn)163 inline typename Par::type Dx_BF(Par const& defn) { return defn.datum_params[0]; }
164 template <typename Par>
Dy_BF(Par const & defn)165 inline typename Par::type Dy_BF(Par const& defn) { return defn.datum_params[1]; }
166 template <typename Par>
Dz_BF(Par const & defn)167 inline typename Par::type Dz_BF(Par const& defn) { return defn.datum_params[2]; }
168 template <typename Par>
Rx_BF(Par const & defn)169 inline typename Par::type Rx_BF(Par const& defn) { return defn.datum_params[3]; }
170 template <typename Par>
Ry_BF(Par const & defn)171 inline typename Par::type Ry_BF(Par const& defn) { return defn.datum_params[4]; }
172 template <typename Par>
Rz_BF(Par const & defn)173 inline typename Par::type Rz_BF(Par const& defn) { return defn.datum_params[5]; }
174 template <typename Par>
M_BF(Par const & defn)175 inline typename Par::type M_BF(Par const& defn) { return defn.datum_params[6]; }
176
177 /*
178 ** This table is intended to indicate for any given error code in
179 ** the range 0 to -56, whether that error will occur for all locations (ie.
180 ** it is a problem with the coordinate system as a whole) in which case the
181 ** value would be 0, or if the problem is with the point being transformed
182 ** in which case the value is 1.
183 **
184 ** At some point we might want to move this array in with the error message
185 ** list or something, but while experimenting with it this should be fine.
186 **
187 **
188 ** NOTE (2017-10-01): Non-transient errors really should have resulted in a
189 ** PJ==0 during initialization, and hence should be handled at the level
190 ** before calling pj_transform. The only obvious example of the contrary
191 ** appears to be the PJD_ERR_GRID_AREA case, which may also be taken to
192 ** mean "no grids available"
193 **
194 **
195 */
196
197 static const int transient_error[60] = {
198 /* 0 1 2 3 4 5 6 7 8 9 */
199 /* 0 to 9 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
200 /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,
201 /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
202 /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
203 /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
204 /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0 };
205
206
207 template <typename T, typename Range>
208 inline int pj_geocentric_to_geodetic( T const& a, T const& es,
209 Range & range );
210 template <typename T, typename Range>
211 inline int pj_geodetic_to_geocentric( T const& a, T const& es,
212 Range & range );
213
214 /************************************************************************/
215 /* pj_transform() */
216 /* */
217 /* Currently this function doesn't recognise if two projections */
218 /* are identical (to short circuit reprojection) because it is */
219 /* difficult to compare PJ structures (since there are some */
220 /* projection specific components). */
221 /************************************************************************/
222
223 template <
224 typename SrcPrj,
225 typename DstPrj2,
226 typename Par,
227 typename Range,
228 typename Grids
229 >
pj_transform(SrcPrj const & srcprj,Par const & srcdefn,DstPrj2 const & dstprj,Par const & dstdefn,Range & range,Grids const & srcgrids,Grids const & dstgrids)230 inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn,
231 DstPrj2 const& dstprj, Par const& dstdefn,
232 Range & range,
233 Grids const& srcgrids,
234 Grids const& dstgrids)
235
236 {
237 typedef typename boost::range_value<Range>::type point_type;
238 typedef typename coordinate_type<point_type>::type coord_t;
239 static const std::size_t dimension = geometry::dimension<point_type>::value;
240 std::size_t point_count = boost::size(range);
241 bool result = true;
242
243 /* -------------------------------------------------------------------- */
244 /* Transform unusual input coordinate axis orientation to */
245 /* standard form if needed. */
246 /* -------------------------------------------------------------------- */
247 // Ignored
248
249 /* -------------------------------------------------------------------- */
250 /* Transform Z to meters if it isn't already. */
251 /* -------------------------------------------------------------------- */
252 if( srcdefn.vto_meter != 1.0 && dimension > 2 )
253 {
254 for( std::size_t i = 0; i < point_count; i++ )
255 {
256 point_type & point = geometry::range::at(range, i);
257 set_z(point, get_z(point) * srcdefn.vto_meter);
258 }
259 }
260
261 /* -------------------------------------------------------------------- */
262 /* Transform geocentric source coordinates to lat/long. */
263 /* -------------------------------------------------------------------- */
264 if( srcdefn.is_geocent )
265 {
266 // Point should be cartesian 3D (ECEF)
267 if (dimension < 3)
268 BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) );
269 //return PJD_ERR_GEOCENTRIC;
270
271 if( srcdefn.to_meter != 1.0 )
272 {
273 for(std::size_t i = 0; i < point_count; i++ )
274 {
275 point_type & point = range::at(range, i);
276 if( ! is_invalid_point(point) )
277 {
278 set<0>(point, get<0>(point) * srcdefn.to_meter);
279 set<1>(point, get<1>(point) * srcdefn.to_meter);
280 }
281 }
282 }
283
284 range_wrapper<Range, false> rng(range);
285 int err = pj_geocentric_to_geodetic( srcdefn.a_orig, srcdefn.es_orig,
286 rng );
287 if( err != 0 )
288 BOOST_THROW_EXCEPTION( projection_exception(err) );
289 //return err;
290
291 // NOTE: here 3D cartesian ECEF is converted into 3D geodetic LLH
292 }
293
294 /* -------------------------------------------------------------------- */
295 /* Transform source points to lat/long, if they aren't */
296 /* already. */
297 /* -------------------------------------------------------------------- */
298 else if( !srcdefn.is_latlong )
299 {
300 // Point should be cartesian 2D or 3D (map projection)
301
302 /* Check first if projection is invertible. */
303 /*if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL))
304 {
305 pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
306 pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR,
307 "pj_transform(): source projection not invertible" );
308 return -17;
309 }*/
310
311 /* If invertible - First try inv3d if defined */
312 //if (srcdefn->inv3d != NULL)
313 //{
314 // /* Three dimensions must be defined */
315 // if ( z == NULL)
316 // {
317 // pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
318 // return PJD_ERR_GEOCENTRIC;
319 // }
320
321 // for (i=0; i < point_count; i++)
322 // {
323 // XYZ projected_loc;
324 // XYZ geodetic_loc;
325
326 // projected_loc.u = x[point_offset*i];
327 // projected_loc.v = y[point_offset*i];
328 // projected_loc.w = z[point_offset*i];
329
330 // if (projected_loc.u == HUGE_VAL)
331 // continue;
332
333 // geodetic_loc = pj_inv3d(projected_loc, srcdefn);
334 // if( srcdefn->ctx->last_errno != 0 )
335 // {
336 // if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
337 // && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
338 // && (srcdefn->ctx->last_errno > 0
339 // || srcdefn->ctx->last_errno < -44 || point_count == 1
340 // || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
341 // return srcdefn->ctx->last_errno;
342 // else
343 // {
344 // geodetic_loc.u = HUGE_VAL;
345 // geodetic_loc.v = HUGE_VAL;
346 // geodetic_loc.w = HUGE_VAL;
347 // }
348 // }
349
350 // x[point_offset*i] = geodetic_loc.u;
351 // y[point_offset*i] = geodetic_loc.v;
352 // z[point_offset*i] = geodetic_loc.w;
353
354 // }
355
356 //}
357 //else
358 {
359 /* Fallback to the original PROJ.4 API 2d inversion - inv */
360 for( std::size_t i = 0; i < point_count; i++ )
361 {
362 point_type & point = range::at(range, i);
363
364 if( is_invalid_point(point) )
365 continue;
366
367 try
368 {
369 pj_inv(srcprj, srcdefn, point, point);
370 }
371 catch(projection_exception const& e)
372 {
373 if( (e.code() != 33 /*EDOM*/
374 && e.code() != 34 /*ERANGE*/ )
375 && (e.code() > 0
376 || e.code() < -44 /*|| point_count == 1*/
377 || transient_error[-e.code()] == 0) ) {
378 BOOST_RETHROW
379 } else {
380 set_invalid_point(point);
381 result = false;
382 if (point_count == 1)
383 return result;
384 }
385 }
386 }
387 }
388 }
389
390 /* -------------------------------------------------------------------- */
391 /* But if they are already lat long, adjust for the prime */
392 /* meridian if there is one in effect. */
393 /* -------------------------------------------------------------------- */
394 if( srcdefn.from_greenwich != 0.0 )
395 {
396 for( std::size_t i = 0; i < point_count; i++ )
397 {
398 point_type & point = range::at(range, i);
399
400 if( ! is_invalid_point(point) )
401 set<0>(point, get<0>(point) + srcdefn.from_greenwich);
402 }
403 }
404
405 /* -------------------------------------------------------------------- */
406 /* Do we need to translate from geoid to ellipsoidal vertical */
407 /* datum? */
408 /* -------------------------------------------------------------------- */
409 /*if( srcdefn->has_geoid_vgrids && z != NULL )
410 {
411 if( pj_apply_vgridshift( srcdefn, "sgeoidgrids",
412 &(srcdefn->vgridlist_geoid),
413 &(srcdefn->vgridlist_geoid_count),
414 0, point_count, point_offset, x, y, z ) != 0 )
415 return pj_ctx_get_errno(srcdefn->ctx);
416 }*/
417
418 /* -------------------------------------------------------------------- */
419 /* Convert datums if needed, and possible. */
420 /* -------------------------------------------------------------------- */
421 if ( ! pj_datum_transform( srcdefn, dstdefn, range, srcgrids, dstgrids ) )
422 {
423 result = false;
424 }
425
426 /* -------------------------------------------------------------------- */
427 /* Do we need to translate from ellipsoidal to geoid vertical */
428 /* datum? */
429 /* -------------------------------------------------------------------- */
430 /*if( dstdefn->has_geoid_vgrids && z != NULL )
431 {
432 if( pj_apply_vgridshift( dstdefn, "sgeoidgrids",
433 &(dstdefn->vgridlist_geoid),
434 &(dstdefn->vgridlist_geoid_count),
435 1, point_count, point_offset, x, y, z ) != 0 )
436 return dstdefn->ctx->last_errno;
437 }*/
438
439 /* -------------------------------------------------------------------- */
440 /* But if they are staying lat long, adjust for the prime */
441 /* meridian if there is one in effect. */
442 /* -------------------------------------------------------------------- */
443 if( dstdefn.from_greenwich != 0.0 )
444 {
445 for( std::size_t i = 0; i < point_count; i++ )
446 {
447 point_type & point = range::at(range, i);
448
449 if( ! is_invalid_point(point) )
450 set<0>(point, get<0>(point) - dstdefn.from_greenwich);
451 }
452 }
453
454 /* -------------------------------------------------------------------- */
455 /* Transform destination latlong to geocentric if required. */
456 /* -------------------------------------------------------------------- */
457 if( dstdefn.is_geocent )
458 {
459 // Point should be cartesian 3D (ECEF)
460 if (dimension < 3)
461 BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) );
462 //return PJD_ERR_GEOCENTRIC;
463
464 // NOTE: In the original code the return value of the following
465 // function is not checked
466 range_wrapper<Range, false> rng(range);
467 int err = pj_geodetic_to_geocentric( dstdefn.a_orig, dstdefn.es_orig,
468 rng );
469 if( err == -14 )
470 result = false;
471 else
472 BOOST_THROW_EXCEPTION( projection_exception(err) );
473
474 if( dstdefn.fr_meter != 1.0 )
475 {
476 for( std::size_t i = 0; i < point_count; i++ )
477 {
478 point_type & point = range::at(range, i);
479 if( ! is_invalid_point(point) )
480 {
481 set<0>(point, get<0>(point) * dstdefn.fr_meter);
482 set<1>(point, get<1>(point) * dstdefn.fr_meter);
483 }
484 }
485 }
486 }
487
488 /* -------------------------------------------------------------------- */
489 /* Transform destination points to projection coordinates, if */
490 /* desired. */
491 /* -------------------------------------------------------------------- */
492 else if( !dstdefn.is_latlong )
493 {
494
495 //if( dstdefn->fwd3d != NULL)
496 //{
497 // for( i = 0; i < point_count; i++ )
498 // {
499 // XYZ projected_loc;
500 // LPZ geodetic_loc;
501
502 // geodetic_loc.u = x[point_offset*i];
503 // geodetic_loc.v = y[point_offset*i];
504 // geodetic_loc.w = z[point_offset*i];
505
506 // if (geodetic_loc.u == HUGE_VAL)
507 // continue;
508
509 // projected_loc = pj_fwd3d( geodetic_loc, dstdefn);
510 // if( dstdefn->ctx->last_errno != 0 )
511 // {
512 // if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
513 // && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
514 // && (dstdefn->ctx->last_errno > 0
515 // || dstdefn->ctx->last_errno < -44 || point_count == 1
516 // || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
517 // return dstdefn->ctx->last_errno;
518 // else
519 // {
520 // projected_loc.u = HUGE_VAL;
521 // projected_loc.v = HUGE_VAL;
522 // projected_loc.w = HUGE_VAL;
523 // }
524 // }
525
526 // x[point_offset*i] = projected_loc.u;
527 // y[point_offset*i] = projected_loc.v;
528 // z[point_offset*i] = projected_loc.w;
529 // }
530
531 //}
532 //else
533 {
534 for(std::size_t i = 0; i < point_count; i++ )
535 {
536 point_type & point = range::at(range, i);
537
538 if( is_invalid_point(point) )
539 continue;
540
541 try {
542 pj_fwd(dstprj, dstdefn, point, point);
543 } catch (projection_exception const& e) {
544
545 if( (e.code() != 33 /*EDOM*/
546 && e.code() != 34 /*ERANGE*/ )
547 && (e.code() > 0
548 || e.code() < -44 /*|| point_count == 1*/
549 || transient_error[-e.code()] == 0) ) {
550 BOOST_RETHROW
551 } else {
552 set_invalid_point(point);
553 result = false;
554 if (point_count == 1)
555 return result;
556 }
557 }
558 }
559 }
560 }
561
562 /* -------------------------------------------------------------------- */
563 /* If a wrapping center other than 0 is provided, rewrap around */
564 /* the suggested center (for latlong coordinate systems only). */
565 /* -------------------------------------------------------------------- */
566 else if( dstdefn.is_latlong && dstdefn.is_long_wrap_set )
567 {
568 for( std::size_t i = 0; i < point_count; i++ )
569 {
570 point_type & point = range::at(range, i);
571 coord_t x = get_as_radian<0>(point);
572
573 if( is_invalid_point(point) )
574 continue;
575
576 // TODO - units-dependant constants could be used instead
577 while( x < dstdefn.long_wrap_center - math::pi<coord_t>() )
578 x += math::two_pi<coord_t>();
579 while( x > dstdefn.long_wrap_center + math::pi<coord_t>() )
580 x -= math::two_pi<coord_t>();
581
582 set_from_radian<0>(point, x);
583 }
584 }
585
586 /* -------------------------------------------------------------------- */
587 /* Transform Z from meters if needed. */
588 /* -------------------------------------------------------------------- */
589 if( dstdefn.vto_meter != 1.0 && dimension > 2 )
590 {
591 for( std::size_t i = 0; i < point_count; i++ )
592 {
593 point_type & point = geometry::range::at(range, i);
594 set_z(point, get_z(point) * dstdefn.vfr_meter);
595 }
596 }
597
598 /* -------------------------------------------------------------------- */
599 /* Transform normalized axes into unusual output coordinate axis */
600 /* orientation if needed. */
601 /* -------------------------------------------------------------------- */
602 // Ignored
603
604 return result;
605 }
606
607 /************************************************************************/
608 /* pj_geodetic_to_geocentric() */
609 /************************************************************************/
610
611 template <typename T, typename Range, bool AddZ>
pj_geodetic_to_geocentric(T const & a,T const & es,range_wrapper<Range,AddZ> & range_wrapper)612 inline int pj_geodetic_to_geocentric( T const& a, T const& es,
613 range_wrapper<Range, AddZ> & range_wrapper )
614
615 {
616 //typedef typename boost::range_iterator<Range>::type iterator;
617 typedef typename boost::range_value<Range>::type point_type;
618 //typedef typename coordinate_type<point_type>::type coord_t;
619
620 Range & rng = range_wrapper.get_range();
621 std::size_t point_count = boost::size(rng);
622
623 int ret_errno = 0;
624
625 T const b = (es == 0.0) ? a : a * sqrt(1-es);
626
627 GeocentricInfo<T> gi;
628 if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 )
629 {
630 return error_geocentric;
631 }
632
633 for( std::size_t i = 0 ; i < point_count ; ++i )
634 {
635 point_type & point = range::at(rng, i);
636
637 if( is_invalid_point(point) )
638 continue;
639
640 T X = 0, Y = 0, Z = 0;
641 if( pj_Convert_Geodetic_To_Geocentric( gi,
642 get_as_radian<0>(point),
643 get_as_radian<1>(point),
644 range_wrapper.get_z(i), // Height
645 X, Y, Z ) != 0 )
646 {
647 ret_errno = error_lat_or_lon_exceed_limit;
648 set_invalid_point(point);
649 /* but keep processing points! */
650 }
651 else
652 {
653 set<0>(point, X);
654 set<1>(point, Y);
655 range_wrapper.set_z(i, Z);
656 }
657 }
658
659 return ret_errno;
660 }
661
662 /************************************************************************/
663 /* pj_geodetic_to_geocentric() */
664 /************************************************************************/
665
666 template <typename T, typename Range, bool AddZ>
pj_geocentric_to_geodetic(T const & a,T const & es,range_wrapper<Range,AddZ> & range_wrapper)667 inline int pj_geocentric_to_geodetic( T const& a, T const& es,
668 range_wrapper<Range, AddZ> & range_wrapper )
669
670 {
671 //typedef typename boost::range_iterator<Range>::type iterator;
672 typedef typename boost::range_value<Range>::type point_type;
673 //typedef typename coordinate_type<point_type>::type coord_t;
674
675 Range & rng = range_wrapper.get_range();
676 std::size_t point_count = boost::size(rng);
677
678 T const b = (es == 0.0) ? a : a * sqrt(1-es);
679
680 GeocentricInfo<T> gi;
681 if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 )
682 {
683 return error_geocentric;
684 }
685
686 for( std::size_t i = 0 ; i < point_count ; ++i )
687 {
688 point_type & point = range::at(rng, i);
689
690 if( is_invalid_point(point) )
691 continue;
692
693 T Longitude = 0, Latitude = 0, Height = 0;
694 pj_Convert_Geocentric_To_Geodetic( gi,
695 get<0>(point),
696 get<1>(point),
697 range_wrapper.get_z(i), // z
698 Longitude, Latitude, Height );
699
700 set_from_radian<0>(point, Longitude);
701 set_from_radian<1>(point, Latitude);
702 range_wrapper.set_z(i, Height); // Height
703 }
704
705 return 0;
706 }
707
708 /************************************************************************/
709 /* pj_compare_datums() */
710 /* */
711 /* Returns TRUE if the two datums are identical, otherwise */
712 /* FALSE. */
713 /************************************************************************/
714
715 template <typename Par>
pj_compare_datums(Par & srcdefn,Par & dstdefn)716 inline bool pj_compare_datums( Par & srcdefn, Par & dstdefn )
717 {
718 if( srcdefn.datum_type != dstdefn.datum_type )
719 {
720 return false;
721 }
722 else if( srcdefn.a_orig != dstdefn.a_orig
723 || math::abs(srcdefn.es_orig - dstdefn.es_orig) > 0.000000000050 )
724 {
725 /* the tolerance for es is to ensure that GRS80 and WGS84 are
726 considered identical */
727 return false;
728 }
729 else if( srcdefn.datum_type == datum_3param )
730 {
731 return (srcdefn.datum_params[0] == dstdefn.datum_params[0]
732 && srcdefn.datum_params[1] == dstdefn.datum_params[1]
733 && srcdefn.datum_params[2] == dstdefn.datum_params[2]);
734 }
735 else if( srcdefn.datum_type == datum_7param )
736 {
737 return (srcdefn.datum_params[0] == dstdefn.datum_params[0]
738 && srcdefn.datum_params[1] == dstdefn.datum_params[1]
739 && srcdefn.datum_params[2] == dstdefn.datum_params[2]
740 && srcdefn.datum_params[3] == dstdefn.datum_params[3]
741 && srcdefn.datum_params[4] == dstdefn.datum_params[4]
742 && srcdefn.datum_params[5] == dstdefn.datum_params[5]
743 && srcdefn.datum_params[6] == dstdefn.datum_params[6]);
744 }
745 else if( srcdefn.datum_type == datum_gridshift )
746 {
747 return srcdefn.nadgrids == dstdefn.nadgrids;
748 }
749 else
750 return true;
751 }
752
753 /************************************************************************/
754 /* pj_geocentic_to_wgs84() */
755 /************************************************************************/
756
757 template <typename Par, typename Range, bool AddZ>
pj_geocentric_to_wgs84(Par const & defn,range_wrapper<Range,AddZ> & range_wrapper)758 inline int pj_geocentric_to_wgs84( Par const& defn,
759 range_wrapper<Range, AddZ> & range_wrapper )
760
761 {
762 typedef typename boost::range_value<Range>::type point_type;
763 typedef typename coordinate_type<point_type>::type coord_t;
764
765 Range & rng = range_wrapper.get_range();
766 std::size_t point_count = boost::size(rng);
767
768 if( defn.datum_type == datum_3param )
769 {
770 for(std::size_t i = 0; i < point_count; i++ )
771 {
772 point_type & point = range::at(rng, i);
773
774 if( is_invalid_point(point) )
775 continue;
776
777 set<0>(point, get<0>(point) + Dx_BF(defn));
778 set<1>(point, get<1>(point) + Dy_BF(defn));
779 range_wrapper.set_z(i, range_wrapper.get_z(i) + Dz_BF(defn));
780 }
781 }
782 else if( defn.datum_type == datum_7param )
783 {
784 for(std::size_t i = 0; i < point_count; i++ )
785 {
786 point_type & point = range::at(rng, i);
787
788 if( is_invalid_point(point) )
789 continue;
790
791 coord_t x = get<0>(point);
792 coord_t y = get<1>(point);
793 coord_t z = range_wrapper.get_z(i);
794
795 coord_t x_out, y_out, z_out;
796
797 x_out = M_BF(defn)*( x - Rz_BF(defn)*y + Ry_BF(defn)*z) + Dx_BF(defn);
798 y_out = M_BF(defn)*( Rz_BF(defn)*x + y - Rx_BF(defn)*z) + Dy_BF(defn);
799 z_out = M_BF(defn)*(-Ry_BF(defn)*x + Rx_BF(defn)*y + z) + Dz_BF(defn);
800
801 set<0>(point, x_out);
802 set<1>(point, y_out);
803 range_wrapper.set_z(i, z_out);
804 }
805 }
806
807 return 0;
808 }
809
810 /************************************************************************/
811 /* pj_geocentic_from_wgs84() */
812 /************************************************************************/
813
814 template <typename Par, typename Range, bool AddZ>
pj_geocentric_from_wgs84(Par const & defn,range_wrapper<Range,AddZ> & range_wrapper)815 inline int pj_geocentric_from_wgs84( Par const& defn,
816 range_wrapper<Range, AddZ> & range_wrapper )
817
818 {
819 typedef typename boost::range_value<Range>::type point_type;
820 typedef typename coordinate_type<point_type>::type coord_t;
821
822 Range & rng = range_wrapper.get_range();
823 std::size_t point_count = boost::size(rng);
824
825 if( defn.datum_type == datum_3param )
826 {
827 for(std::size_t i = 0; i < point_count; i++ )
828 {
829 point_type & point = range::at(rng, i);
830
831 if( is_invalid_point(point) )
832 continue;
833
834 set<0>(point, get<0>(point) - Dx_BF(defn));
835 set<1>(point, get<1>(point) - Dy_BF(defn));
836 range_wrapper.set_z(i, range_wrapper.get_z(i) - Dz_BF(defn));
837 }
838 }
839 else if( defn.datum_type == datum_7param )
840 {
841 for(std::size_t i = 0; i < point_count; i++ )
842 {
843 point_type & point = range::at(rng, i);
844
845 if( is_invalid_point(point) )
846 continue;
847
848 coord_t x = get<0>(point);
849 coord_t y = get<1>(point);
850 coord_t z = range_wrapper.get_z(i);
851
852 coord_t x_tmp = (x - Dx_BF(defn)) / M_BF(defn);
853 coord_t y_tmp = (y - Dy_BF(defn)) / M_BF(defn);
854 coord_t z_tmp = (z - Dz_BF(defn)) / M_BF(defn);
855
856 x = x_tmp + Rz_BF(defn)*y_tmp - Ry_BF(defn)*z_tmp;
857 y = -Rz_BF(defn)*x_tmp + y_tmp + Rx_BF(defn)*z_tmp;
858 z = Ry_BF(defn)*x_tmp - Rx_BF(defn)*y_tmp + z_tmp;
859
860 set<0>(point, x);
861 set<1>(point, y);
862 range_wrapper.set_z(i, z);
863 }
864 }
865
866 return 0;
867 }
868
869
pj_datum_check_error(int err)870 inline bool pj_datum_check_error(int err)
871 {
872 return err != 0 && (err > 0 || transient_error[-err] == 0);
873 }
874
875 /************************************************************************/
876 /* pj_datum_transform() */
877 /* */
878 /* The input should be long/lat/z coordinates in radians in the */
879 /* source datum, and the output should be long/lat/z */
880 /* coordinates in radians in the destination datum. */
881 /************************************************************************/
882
883 template <typename Par, typename Range, typename Grids>
pj_datum_transform(Par const & srcdefn,Par const & dstdefn,Range & range,Grids const & srcgrids,Grids const & dstgrids)884 inline bool pj_datum_transform(Par const& srcdefn,
885 Par const& dstdefn,
886 Range & range,
887 Grids const& srcgrids,
888 Grids const& dstgrids)
889
890 {
891 typedef typename Par::type calc_t;
892
893 // This has to be consistent with default spheroid and pj_ellps
894 // TODO: Define in one place
895 static const calc_t wgs84_a = 6378137.0;
896 static const calc_t wgs84_b = 6356752.3142451793;
897 static const calc_t wgs84_es = 1. - (wgs84_b * wgs84_b) / (wgs84_a * wgs84_a);
898
899 bool result = true;
900
901 calc_t src_a, src_es, dst_a, dst_es;
902
903 /* -------------------------------------------------------------------- */
904 /* We cannot do any meaningful datum transformation if either */
905 /* the source or destination are of an unknown datum type */
906 /* (ie. only a +ellps declaration, no +datum). This is new */
907 /* behavior for PROJ 4.6.0. */
908 /* -------------------------------------------------------------------- */
909 if( srcdefn.datum_type == datum_unknown
910 || dstdefn.datum_type == datum_unknown )
911 return result;
912
913 /* -------------------------------------------------------------------- */
914 /* Short cut if the datums are identical. */
915 /* -------------------------------------------------------------------- */
916 if( pj_compare_datums( srcdefn, dstdefn ) )
917 return result;
918
919 src_a = srcdefn.a_orig;
920 src_es = srcdefn.es_orig;
921
922 dst_a = dstdefn.a_orig;
923 dst_es = dstdefn.es_orig;
924
925 /* -------------------------------------------------------------------- */
926 /* Create a temporary Z array if one is not provided. */
927 /* -------------------------------------------------------------------- */
928
929 range_wrapper<Range> z_range(range);
930
931 /* -------------------------------------------------------------------- */
932 /* If this datum requires grid shifts, then apply it to geodetic */
933 /* coordinates. */
934 /* -------------------------------------------------------------------- */
935 if( srcdefn.datum_type == datum_gridshift )
936 {
937 try {
938 pj_apply_gridshift_2<false>( srcdefn, range, srcgrids );
939 } catch (projection_exception const& e) {
940 if (pj_datum_check_error(e.code())) {
941 BOOST_RETHROW
942 }
943 }
944
945 src_a = wgs84_a;
946 src_es = wgs84_es;
947 }
948
949 if( dstdefn.datum_type == datum_gridshift )
950 {
951 dst_a = wgs84_a;
952 dst_es = wgs84_es;
953 }
954
955 /* ==================================================================== */
956 /* Do we need to go through geocentric coordinates? */
957 /* ==================================================================== */
958 if( src_es != dst_es || src_a != dst_a
959 || srcdefn.datum_type == datum_3param
960 || srcdefn.datum_type == datum_7param
961 || dstdefn.datum_type == datum_3param
962 || dstdefn.datum_type == datum_7param)
963 {
964 /* -------------------------------------------------------------------- */
965 /* Convert to geocentric coordinates. */
966 /* -------------------------------------------------------------------- */
967 int err = pj_geodetic_to_geocentric( src_a, src_es, z_range );
968 if (pj_datum_check_error(err))
969 BOOST_THROW_EXCEPTION( projection_exception(err) );
970 else if (err != 0)
971 result = false;
972
973 /* -------------------------------------------------------------------- */
974 /* Convert between datums. */
975 /* -------------------------------------------------------------------- */
976 if( srcdefn.datum_type == datum_3param
977 || srcdefn.datum_type == datum_7param )
978 {
979 try {
980 pj_geocentric_to_wgs84( srcdefn, z_range );
981 } catch (projection_exception const& e) {
982 if (pj_datum_check_error(e.code())) {
983 BOOST_RETHROW
984 }
985 }
986 }
987
988 if( dstdefn.datum_type == datum_3param
989 || dstdefn.datum_type == datum_7param )
990 {
991 try {
992 pj_geocentric_from_wgs84( dstdefn, z_range );
993 } catch (projection_exception const& e) {
994 if (pj_datum_check_error(e.code())) {
995 BOOST_RETHROW
996 }
997 }
998 }
999
1000 /* -------------------------------------------------------------------- */
1001 /* Convert back to geodetic coordinates. */
1002 /* -------------------------------------------------------------------- */
1003 err = pj_geocentric_to_geodetic( dst_a, dst_es, z_range );
1004 if (pj_datum_check_error(err))
1005 BOOST_THROW_EXCEPTION( projection_exception(err) );
1006 else if (err != 0)
1007 result = false;
1008 }
1009
1010 /* -------------------------------------------------------------------- */
1011 /* Apply grid shift to destination if required. */
1012 /* -------------------------------------------------------------------- */
1013 if( dstdefn.datum_type == datum_gridshift )
1014 {
1015 try {
1016 pj_apply_gridshift_2<true>( dstdefn, range, dstgrids );
1017 } catch (projection_exception const& e) {
1018 if (pj_datum_check_error(e.code()))
1019 BOOST_RETHROW
1020 }
1021 }
1022
1023 return result;
1024 }
1025
1026 } // namespace detail
1027
1028 }}} // namespace boost::geometry::projections
1029
1030 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
1031