]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/boost/geometry/srs/projections/impl/pj_transform.hpp
projects / ceph.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
import new upstream nautilus stable release 14.2.8
[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / impl / pj_transform.hpp
CommitLineData
11fdf7f2
TL		 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>
92f5a8d4 44#include <boost/geometry/core/coordinate_dimension.hpp>
11fdf7f2
TL		 45#include <boost/geometry/core/radian_access.hpp>
46
47#include <boost/geometry/srs/projections/impl/geocent.hpp>
92f5a8d4 48#include <boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp>
11fdf7f2
TL		 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
58namespace boost { namespace geometry { namespace projections
59{
60
61namespace detail
62{
63
64// -----------------------------------------------------------
65// Boost.Geometry helpers begin
66// -----------------------------------------------------------
67
68template
69<
70 typename Point,
71 bool HasCoord2 = (dimension<Point>::value > 2)
72>
73struct z_access
74{
75 typedef typename coordinate_type<Point>::type type;
76 static inline type get(Point const& point)
77 {
78 return geometry::get<2>(point);
79 }
80 static inline void set(Point & point, type const& h)
81 {
82 return geometry::set<2>(point, h);
83 }
84};
85
86template <typename Point>
87struct z_access<Point, false>
88{
89 typedef typename coordinate_type<Point>::type type;
90 static inline type get(Point const& )
91 {
92 return type(0);
93 }
94 static inline void set(Point & , type const& )
95 {}
96};
97
98template <typename XYorXYZ>
99inline typename z_access<XYorXYZ>::type
100get_z(XYorXYZ const& xy_or_xyz)
101{
102 return z_access<XYorXYZ>::get(xy_or_xyz);
103}
104
105template <typename XYorXYZ>
106inline 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
112template
113<
114 typename Range,
115 bool AddZ = (dimension<typename boost::range_value<Range>::type>::value < 3)
116>
117struct 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
123 range_wrapper(Range & range)
124 : m_range(range)
125 {}
126
127 range_type & get_range() { return m_range; }
128
129 coord_t get_z(std::size_t i) { return detail::get_z(range::at(m_range, i)); }
130 void set_z(std::size_t i, coord_t const& v) { return detail::set_z(range::at(m_range, i), v); }
131
132private:
133 Range & m_range;
134};
135
136template <typename Range>
137struct 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
143 range_wrapper(Range & range)
144 : m_range(range)
145 , m_zs(boost::size(range), coord_t(0))
146 {}
147
148 range_type & get_range() { return m_range; }
149
150 coord_t get_z(std::size_t i) { return m_zs[i]; }
151 void set_z(std::size_t i, coord_t const& v) { m_zs[i] = v; }
152
153private:
154 Range & m_range;
155 std::vector<coord_t> m_zs;
156};
157
158// -----------------------------------------------------------
159// Boost.Geometry helpers end
160// -----------------------------------------------------------
161
11fdf7f2
TL		 162template <typename Par>
163inline typename Par::type Dx_BF(Par const& defn) { return defn.datum_params[0]; }
164template <typename Par>
165inline typename Par::type Dy_BF(Par const& defn) { return defn.datum_params[1]; }
166template <typename Par>
167inline typename Par::type Dz_BF(Par const& defn) { return defn.datum_params[2]; }
168template <typename Par>
169inline typename Par::type Rx_BF(Par const& defn) { return defn.datum_params[3]; }
170template <typename Par>
171inline typename Par::type Ry_BF(Par const& defn) { return defn.datum_params[4]; }
172template <typename Par>
173inline typename Par::type Rz_BF(Par const& defn) { return defn.datum_params[5]; }
174template <typename Par>
175inline 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
92f5a8d4 179** the range 0 to -56, whether that error will occur for all locations (ie.
11fdf7f2
TL		 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.
92f5a8d4
TL		 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**
11fdf7f2
TL		 195*/
196
92f5a8d4 197static const int transient_error[60] = {
11fdf7f2
TL		 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,
92f5a8d4
TL		 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 };
11fdf7f2
TL		 205
206
207template <typename T, typename Range>
208inline int pj_geocentric_to_geodetic( T const& a, T const& es,
209 Range & range );
210template <typename T, typename Range>
211inline 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
92f5a8d4
TL		 223template <
224 typename SrcPrj,
225 typename DstPrj2,
226 typename Par,
227 typename Range,
228 typename Grids
229>
11fdf7f2
TL		 230inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn,
231 DstPrj2 const& dstprj, Par const& dstdefn,
92f5a8d4
TL		 232 Range & range,
233 Grids const& srcgrids,
234 Grids const& dstgrids)
11fdf7f2
TL		 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)
92f5a8d4 268 BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) );
11fdf7f2
TL		 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/* -------------------------------------------------------------------- */
92f5a8d4 421 if ( ! pj_datum_transform( srcdefn, dstdefn, range, srcgrids, dstgrids ) )
11fdf7f2
TL		 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)
92f5a8d4 461 BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) );
11fdf7f2
TL		 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
611template <typename T, typename Range, bool AddZ>
612inline 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 {
92f5a8d4 630 return error_geocentric;
11fdf7f2
TL		 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 {
92f5a8d4 647 ret_errno = error_lat_or_lon_exceed_limit;
11fdf7f2
TL		 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
666template <typename T, typename Range, bool AddZ>
667inline 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 {
92f5a8d4 683 return error_geocentric;
11fdf7f2
TL		 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
715template <typename Par>
716inline 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 }
92f5a8d4 729 else if( srcdefn.datum_type == datum_3param )
11fdf7f2
TL		 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 }
92f5a8d4 735 else if( srcdefn.datum_type == datum_7param )
11fdf7f2
TL		 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 }
92f5a8d4 745 else if( srcdefn.datum_type == datum_gridshift )
11fdf7f2 746 {
92f5a8d4 747 return srcdefn.nadgrids == dstdefn.nadgrids;
11fdf7f2
TL		 748 }
749 else
750 return true;
751}
752
753/************************************************************************/
754/* pj_geocentic_to_wgs84() */
755/************************************************************************/
756
757template <typename Par, typename Range, bool AddZ>
758inline 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
92f5a8d4 768 if( defn.datum_type == datum_3param )
11fdf7f2
TL		 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 }
92f5a8d4 782 else if( defn.datum_type == datum_7param )
11fdf7f2
TL		 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
814template <typename Par, typename Range, bool AddZ>
815inline 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
92f5a8d4 825 if( defn.datum_type == datum_3param )
11fdf7f2
TL		 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 }
92f5a8d4 839 else if( defn.datum_type == datum_7param )
11fdf7f2
TL		 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
870inline 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
92f5a8d4
TL		 883template <typename Par, typename Range, typename Grids>
884inline bool pj_datum_transform(Par const& srcdefn,
885 Par const& dstdefn,
886 Range & range,
887 Grids const& srcgrids,
888 Grids const& dstgrids)
11fdf7f2
TL		 889
890{
891 typedef typename Par::type calc_t;
92f5a8d4
TL		 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
11fdf7f2
TL		 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/* -------------------------------------------------------------------- */
92f5a8d4
TL		 909 if( srcdefn.datum_type == datum_unknown
910 || dstdefn.datum_type == datum_unknown )
11fdf7f2
TL		 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/* -------------------------------------------------------------------- */
92f5a8d4 935 if( srcdefn.datum_type == datum_gridshift )
11fdf7f2
TL		 936 {
937 try {
92f5a8d4 938 pj_apply_gridshift_2<false>( srcdefn, range, srcgrids );
11fdf7f2
TL		 939 } catch (projection_exception const& e) {
940 if (pj_datum_check_error(e.code())) {
941 BOOST_RETHROW
942 }
943 }
944
92f5a8d4
TL		 945 src_a = wgs84_a;
946 src_es = wgs84_es;
11fdf7f2
TL		 947 }
948
92f5a8d4 949 if( dstdefn.datum_type == datum_gridshift )
11fdf7f2 950 {
92f5a8d4
TL		 951 dst_a = wgs84_a;
952 dst_es = wgs84_es;
953 }
11fdf7f2
TL		 954
955/* ==================================================================== */
956/* Do we need to go through geocentric coordinates? */
957/* ==================================================================== */
958 if( src_es != dst_es || src_a != dst_a
92f5a8d4
TL		 959 || srcdefn.datum_type == datum_3param
960 || srcdefn.datum_type == datum_7param
961 || dstdefn.datum_type == datum_3param
962 || dstdefn.datum_type == datum_7param)
11fdf7f2
TL		 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/* -------------------------------------------------------------------- */
92f5a8d4
TL		 976 if( srcdefn.datum_type == datum_3param
977 || srcdefn.datum_type == datum_7param )
11fdf7f2
TL		 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
92f5a8d4
TL		 988 if( dstdefn.datum_type == datum_3param
989 || dstdefn.datum_type == datum_7param )
11fdf7f2
TL		 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/* -------------------------------------------------------------------- */
92f5a8d4 1013 if( dstdefn.datum_type == datum_gridshift )
11fdf7f2
TL		 1014 {
1015 try {
92f5a8d4 1016 pj_apply_gridshift_2<true>( dstdefn, range, dstgrids );
11fdf7f2
TL		 1017 } catch (projection_exception const& e) {
1018 if (pj_datum_check_error(e.code()))
1019 BOOST_RETHROW
1020 }
92f5a8d4 1021 }
11fdf7f2
TL		 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
Ceph