summaryrefslogtreecommitdiff
path: root/boost/geometry/srs/projections/impl/pj_transform.hpp
diff options
context:
space:
mode:
authorDongHun Kwak <dh0128.kwak@samsung.com>2019-12-05 15:12:59 +0900
committerDongHun Kwak <dh0128.kwak@samsung.com>2019-12-05 15:12:59 +0900
commitb8cf34c691623e4ec329053cbbf68522a855882d (patch)
tree34da08632a99677f6b79ecb65e5b655a5b69a67f /boost/geometry/srs/projections/impl/pj_transform.hpp
parent3fdc3e5ee96dca5b11d1694975a65200787eab86 (diff)
downloadboost-b8cf34c691623e4ec329053cbbf68522a855882d.tar.gz
boost-b8cf34c691623e4ec329053cbbf68522a855882d.tar.bz2
boost-b8cf34c691623e4ec329053cbbf68522a855882d.zip
Imported Upstream version 1.67.0upstream/1.67.0
Diffstat (limited to 'boost/geometry/srs/projections/impl/pj_transform.hpp')
-rw-r--r--boost/geometry/srs/projections/impl/pj_transform.hpp1009
1 files changed, 1009 insertions, 0 deletions
diff --git a/boost/geometry/srs/projections/impl/pj_transform.hpp b/boost/geometry/srs/projections/impl/pj_transform.hpp
new file mode 100644
index 0000000000..8c2095642f
--- /dev/null
+++ b/boost/geometry/srs/projections/impl/pj_transform.hpp
@@ -0,0 +1,1009 @@
+// Boost.Geometry
+// This file is manually converted from PROJ4
+
+// This file was modified by Oracle on 2017, 2018.
+// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+// This file is converted from PROJ4, http://trac.osgeo.org/proj
+// PROJ4 is originally written by Gerald Evenden (then of the USGS)
+// PROJ4 is maintained by Frank Warmerdam
+// This file was converted to Geometry Library by Adam Wulkiewicz
+
+// Original copyright notice:
+
+// Copyright (c) 2000, Frank Warmerdam
+
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+
+#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
+#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
+
+
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+
+#include <boost/geometry/srs/projections/impl/geocent.hpp>
+#include <boost/geometry/srs/projections/impl/projects.hpp>
+#include <boost/geometry/srs/projections/invalid_point.hpp>
+
+#include <boost/geometry/util/range.hpp>
+
+#include <cstring>
+#include <cmath>
+
+
+namespace boost { namespace geometry { namespace projections
+{
+
+namespace detail
+{
+
+// -----------------------------------------------------------
+// Boost.Geometry helpers begin
+// -----------------------------------------------------------
+
+template
+<
+ typename Point,
+ bool HasCoord2 = (dimension<Point>::value > 2)
+>
+struct z_access
+{
+ typedef typename coordinate_type<Point>::type type;
+ static inline type get(Point const& point)
+ {
+ return geometry::get<2>(point);
+ }
+ static inline void set(Point & point, type const& h)
+ {
+ return geometry::set<2>(point, h);
+ }
+};
+
+template <typename Point>
+struct z_access<Point, false>
+{
+ typedef typename coordinate_type<Point>::type type;
+ static inline type get(Point const& )
+ {
+ return type(0);
+ }
+ static inline void set(Point & , type const& )
+ {}
+};
+
+template <typename XYorXYZ>
+inline typename z_access<XYorXYZ>::type
+get_z(XYorXYZ const& xy_or_xyz)
+{
+ return z_access<XYorXYZ>::get(xy_or_xyz);
+}
+
+template <typename XYorXYZ>
+inline void set_z(XYorXYZ & xy_or_xyz,
+ typename z_access<XYorXYZ>::type const& z)
+{
+ return z_access<XYorXYZ>::set(xy_or_xyz, z);
+}
+
+template
+<
+ typename Range,
+ bool AddZ = (dimension<typename boost::range_value<Range>::type>::value < 3)
+>
+struct range_wrapper
+{
+ typedef Range range_type;
+ typedef typename boost::range_value<Range>::type point_type;
+ typedef typename coordinate_type<point_type>::type coord_t;
+
+ range_wrapper(Range & range)
+ : m_range(range)
+ {}
+
+ range_type & get_range() { return m_range; }
+
+ coord_t get_z(std::size_t i) { return detail::get_z(range::at(m_range, i)); }
+ void set_z(std::size_t i, coord_t const& v) { return detail::set_z(range::at(m_range, i), v); }
+
+private:
+ Range & m_range;
+};
+
+template <typename Range>
+struct range_wrapper<Range, true>
+{
+ typedef Range range_type;
+ typedef typename boost::range_value<Range>::type point_type;
+ typedef typename coordinate_type<point_type>::type coord_t;
+
+ range_wrapper(Range & range)
+ : m_range(range)
+ , m_zs(boost::size(range), coord_t(0))
+ {}
+
+ range_type & get_range() { return m_range; }
+
+ coord_t get_z(std::size_t i) { return m_zs[i]; }
+ void set_z(std::size_t i, coord_t const& v) { m_zs[i] = v; }
+
+private:
+ Range & m_range;
+ std::vector<coord_t> m_zs;
+};
+
+// -----------------------------------------------------------
+// Boost.Geometry helpers end
+// -----------------------------------------------------------
+
+/*#ifndef SRS_WGS84_SEMIMAJOR
+#define SRS_WGS84_SEMIMAJOR 6378137.0
+#endif
+
+#ifndef SRS_WGS84_ESQUARED
+#define SRS_WGS84_ESQUARED 0.0066943799901413165
+#endif*/
+
+template <typename Par>
+inline typename Par::type Dx_BF(Par const& defn) { return defn.datum_params[0]; }
+template <typename Par>
+inline typename Par::type Dy_BF(Par const& defn) { return defn.datum_params[1]; }
+template <typename Par>
+inline typename Par::type Dz_BF(Par const& defn) { return defn.datum_params[2]; }
+template <typename Par>
+inline typename Par::type Rx_BF(Par const& defn) { return defn.datum_params[3]; }
+template <typename Par>
+inline typename Par::type Ry_BF(Par const& defn) { return defn.datum_params[4]; }
+template <typename Par>
+inline typename Par::type Rz_BF(Par const& defn) { return defn.datum_params[5]; }
+template <typename Par>
+inline typename Par::type M_BF(Par const& defn) { return defn.datum_params[6]; }
+
+/*
+** This table is intended to indicate for any given error code in
+** the range 0 to -44, whether that error will occur for all locations (ie.
+** it is a problem with the coordinate system as a whole) in which case the
+** value would be 0, or if the problem is with the point being transformed
+** in which case the value is 1.
+**
+** At some point we might want to move this array in with the error message
+** list or something, but while experimenting with it this should be fine.
+*/
+
+static const int transient_error[50] = {
+ /* 0 1 2 3 4 5 6 7 8 9 */
+ /* 0 to 9 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,
+ /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
+ /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 };
+
+
+template <typename T, typename Range>
+inline int pj_geocentric_to_geodetic( T const& a, T const& es,
+ Range & range );
+template <typename T, typename Range>
+inline int pj_geodetic_to_geocentric( T const& a, T const& es,
+ Range & range );
+
+/************************************************************************/
+/* pj_transform() */
+/* */
+/* Currently this function doesn't recognise if two projections */
+/* are identical (to short circuit reprojection) because it is */
+/* difficult to compare PJ structures (since there are some */
+/* projection specific components). */
+/************************************************************************/
+
+template <typename SrcPrj, typename DstPrj2, typename Par, typename Range>
+inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn,
+ DstPrj2 const& dstprj, Par const& dstdefn,
+ Range & range)
+
+{
+ typedef typename boost::range_value<Range>::type point_type;
+ typedef typename coordinate_type<point_type>::type coord_t;
+ static const std::size_t dimension = geometry::dimension<point_type>::value;
+ std::size_t point_count = boost::size(range);
+ bool result = true;
+
+/* -------------------------------------------------------------------- */
+/* Transform unusual input coordinate axis orientation to */
+/* standard form if needed. */
+/* -------------------------------------------------------------------- */
+ // Ignored
+
+/* -------------------------------------------------------------------- */
+/* Transform Z to meters if it isn't already. */
+/* -------------------------------------------------------------------- */
+ if( srcdefn.vto_meter != 1.0 && dimension > 2 )
+ {
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = geometry::range::at(range, i);
+ set_z(point, get_z(point) * srcdefn.vto_meter);
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Transform geocentric source coordinates to lat/long. */
+/* -------------------------------------------------------------------- */
+ if( srcdefn.is_geocent )
+ {
+ // Point should be cartesian 3D (ECEF)
+ if (dimension < 3)
+ BOOST_THROW_EXCEPTION( projection_exception(PJD_ERR_GEOCENTRIC) );
+ //return PJD_ERR_GEOCENTRIC;
+
+ if( srcdefn.to_meter != 1.0 )
+ {
+ for(std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+ if( ! is_invalid_point(point) )
+ {
+ set<0>(point, get<0>(point) * srcdefn.to_meter);
+ set<1>(point, get<1>(point) * srcdefn.to_meter);
+ }
+ }
+ }
+
+ range_wrapper<Range, false> rng(range);
+ int err = pj_geocentric_to_geodetic( srcdefn.a_orig, srcdefn.es_orig,
+ rng );
+ if( err != 0 )
+ BOOST_THROW_EXCEPTION( projection_exception(err) );
+ //return err;
+
+ // NOTE: here 3D cartesian ECEF is converted into 3D geodetic LLH
+ }
+
+/* -------------------------------------------------------------------- */
+/* Transform source points to lat/long, if they aren't */
+/* already. */
+/* -------------------------------------------------------------------- */
+ else if( !srcdefn.is_latlong )
+ {
+ // Point should be cartesian 2D or 3D (map projection)
+
+ /* Check first if projection is invertible. */
+ /*if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL))
+ {
+ pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
+ pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR,
+ "pj_transform(): source projection not invertible" );
+ return -17;
+ }*/
+
+ /* If invertible - First try inv3d if defined */
+ //if (srcdefn->inv3d != NULL)
+ //{
+ // /* Three dimensions must be defined */
+ // if ( z == NULL)
+ // {
+ // pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
+ // return PJD_ERR_GEOCENTRIC;
+ // }
+
+ // for (i=0; i < point_count; i++)
+ // {
+ // XYZ projected_loc;
+ // XYZ geodetic_loc;
+
+ // projected_loc.u = x[point_offset*i];
+ // projected_loc.v = y[point_offset*i];
+ // projected_loc.w = z[point_offset*i];
+
+ // if (projected_loc.u == HUGE_VAL)
+ // continue;
+
+ // geodetic_loc = pj_inv3d(projected_loc, srcdefn);
+ // if( srcdefn->ctx->last_errno != 0 )
+ // {
+ // if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
+ // && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
+ // && (srcdefn->ctx->last_errno > 0
+ // || srcdefn->ctx->last_errno < -44 || point_count == 1
+ // || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
+ // return srcdefn->ctx->last_errno;
+ // else
+ // {
+ // geodetic_loc.u = HUGE_VAL;
+ // geodetic_loc.v = HUGE_VAL;
+ // geodetic_loc.w = HUGE_VAL;
+ // }
+ // }
+
+ // x[point_offset*i] = geodetic_loc.u;
+ // y[point_offset*i] = geodetic_loc.v;
+ // z[point_offset*i] = geodetic_loc.w;
+
+ // }
+
+ //}
+ //else
+ {
+ /* Fallback to the original PROJ.4 API 2d inversion - inv */
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ try
+ {
+ pj_inv(srcprj, srcdefn, point, point);
+ }
+ catch(projection_exception const& e)
+ {
+ if( (e.code() != 33 /*EDOM*/
+ && e.code() != 34 /*ERANGE*/ )
+ && (e.code() > 0
+ || e.code() < -44 /*|| point_count == 1*/
+ || transient_error[-e.code()] == 0) ) {
+ BOOST_RETHROW
+ } else {
+ set_invalid_point(point);
+ result = false;
+ if (point_count == 1)
+ return result;
+ }
+ }
+ }
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* But if they are already lat long, adjust for the prime */
+/* meridian if there is one in effect. */
+/* -------------------------------------------------------------------- */
+ if( srcdefn.from_greenwich != 0.0 )
+ {
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+
+ if( ! is_invalid_point(point) )
+ set<0>(point, get<0>(point) + srcdefn.from_greenwich);
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Do we need to translate from geoid to ellipsoidal vertical */
+/* datum? */
+/* -------------------------------------------------------------------- */
+ /*if( srcdefn->has_geoid_vgrids && z != NULL )
+ {
+ if( pj_apply_vgridshift( srcdefn, "sgeoidgrids",
+ &(srcdefn->vgridlist_geoid),
+ &(srcdefn->vgridlist_geoid_count),
+ 0, point_count, point_offset, x, y, z ) != 0 )
+ return pj_ctx_get_errno(srcdefn->ctx);
+ }*/
+
+/* -------------------------------------------------------------------- */
+/* Convert datums if needed, and possible. */
+/* -------------------------------------------------------------------- */
+ if ( ! pj_datum_transform( srcdefn, dstdefn, range ) )
+ {
+ result = false;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Do we need to translate from ellipsoidal to geoid vertical */
+/* datum? */
+/* -------------------------------------------------------------------- */
+ /*if( dstdefn->has_geoid_vgrids && z != NULL )
+ {
+ if( pj_apply_vgridshift( dstdefn, "sgeoidgrids",
+ &(dstdefn->vgridlist_geoid),
+ &(dstdefn->vgridlist_geoid_count),
+ 1, point_count, point_offset, x, y, z ) != 0 )
+ return dstdefn->ctx->last_errno;
+ }*/
+
+/* -------------------------------------------------------------------- */
+/* But if they are staying lat long, adjust for the prime */
+/* meridian if there is one in effect. */
+/* -------------------------------------------------------------------- */
+ if( dstdefn.from_greenwich != 0.0 )
+ {
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+
+ if( ! is_invalid_point(point) )
+ set<0>(point, get<0>(point) - dstdefn.from_greenwich);
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Transform destination latlong to geocentric if required. */
+/* -------------------------------------------------------------------- */
+ if( dstdefn.is_geocent )
+ {
+ // Point should be cartesian 3D (ECEF)
+ if (dimension < 3)
+ BOOST_THROW_EXCEPTION( projection_exception(PJD_ERR_GEOCENTRIC) );
+ //return PJD_ERR_GEOCENTRIC;
+
+ // NOTE: In the original code the return value of the following
+ // function is not checked
+ range_wrapper<Range, false> rng(range);
+ int err = pj_geodetic_to_geocentric( dstdefn.a_orig, dstdefn.es_orig,
+ rng );
+ if( err == -14 )
+ result = false;
+ else
+ BOOST_THROW_EXCEPTION( projection_exception(err) );
+
+ if( dstdefn.fr_meter != 1.0 )
+ {
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+ if( ! is_invalid_point(point) )
+ {
+ set<0>(point, get<0>(point) * dstdefn.fr_meter);
+ set<1>(point, get<1>(point) * dstdefn.fr_meter);
+ }
+ }
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Transform destination points to projection coordinates, if */
+/* desired. */
+/* -------------------------------------------------------------------- */
+ else if( !dstdefn.is_latlong )
+ {
+
+ //if( dstdefn->fwd3d != NULL)
+ //{
+ // for( i = 0; i < point_count; i++ )
+ // {
+ // XYZ projected_loc;
+ // LPZ geodetic_loc;
+
+ // geodetic_loc.u = x[point_offset*i];
+ // geodetic_loc.v = y[point_offset*i];
+ // geodetic_loc.w = z[point_offset*i];
+
+ // if (geodetic_loc.u == HUGE_VAL)
+ // continue;
+
+ // projected_loc = pj_fwd3d( geodetic_loc, dstdefn);
+ // if( dstdefn->ctx->last_errno != 0 )
+ // {
+ // if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
+ // && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
+ // && (dstdefn->ctx->last_errno > 0
+ // || dstdefn->ctx->last_errno < -44 || point_count == 1
+ // || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
+ // return dstdefn->ctx->last_errno;
+ // else
+ // {
+ // projected_loc.u = HUGE_VAL;
+ // projected_loc.v = HUGE_VAL;
+ // projected_loc.w = HUGE_VAL;
+ // }
+ // }
+
+ // x[point_offset*i] = projected_loc.u;
+ // y[point_offset*i] = projected_loc.v;
+ // z[point_offset*i] = projected_loc.w;
+ // }
+
+ //}
+ //else
+ {
+ for(std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ try {
+ pj_fwd(dstprj, dstdefn, point, point);
+ } catch (projection_exception const& e) {
+
+ if( (e.code() != 33 /*EDOM*/
+ && e.code() != 34 /*ERANGE*/ )
+ && (e.code() > 0
+ || e.code() < -44 /*|| point_count == 1*/
+ || transient_error[-e.code()] == 0) ) {
+ BOOST_RETHROW
+ } else {
+ set_invalid_point(point);
+ result = false;
+ if (point_count == 1)
+ return result;
+ }
+ }
+ }
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* If a wrapping center other than 0 is provided, rewrap around */
+/* the suggested center (for latlong coordinate systems only). */
+/* -------------------------------------------------------------------- */
+ else if( dstdefn.is_latlong && dstdefn.is_long_wrap_set )
+ {
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(range, i);
+ coord_t x = get_as_radian<0>(point);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ // TODO - units-dependant constants could be used instead
+ while( x < dstdefn.long_wrap_center - math::pi<coord_t>() )
+ x += math::two_pi<coord_t>();
+ while( x > dstdefn.long_wrap_center + math::pi<coord_t>() )
+ x -= math::two_pi<coord_t>();
+
+ set_from_radian<0>(point, x);
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Transform Z from meters if needed. */
+/* -------------------------------------------------------------------- */
+ if( dstdefn.vto_meter != 1.0 && dimension > 2 )
+ {
+ for( std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = geometry::range::at(range, i);
+ set_z(point, get_z(point) * dstdefn.vfr_meter);
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Transform normalized axes into unusual output coordinate axis */
+/* orientation if needed. */
+/* -------------------------------------------------------------------- */
+ // Ignored
+
+ return result;
+}
+
+/************************************************************************/
+/* pj_geodetic_to_geocentric() */
+/************************************************************************/
+
+template <typename T, typename Range, bool AddZ>
+inline int pj_geodetic_to_geocentric( T const& a, T const& es,
+ range_wrapper<Range, AddZ> & range_wrapper )
+
+{
+ //typedef typename boost::range_iterator<Range>::type iterator;
+ typedef typename boost::range_value<Range>::type point_type;
+ //typedef typename coordinate_type<point_type>::type coord_t;
+
+ Range & rng = range_wrapper.get_range();
+ std::size_t point_count = boost::size(rng);
+
+ int ret_errno = 0;
+
+ T const b = (es == 0.0) ? a : a * sqrt(1-es);
+
+ GeocentricInfo<T> gi;
+ if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 )
+ {
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+ for( std::size_t i = 0 ; i < point_count ; ++i )
+ {
+ point_type & point = range::at(rng, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ T X = 0, Y = 0, Z = 0;
+ if( pj_Convert_Geodetic_To_Geocentric( gi,
+ get_as_radian<0>(point),
+ get_as_radian<1>(point),
+ range_wrapper.get_z(i), // Height
+ X, Y, Z ) != 0 )
+ {
+ ret_errno = -14;
+ set_invalid_point(point);
+ /* but keep processing points! */
+ }
+ else
+ {
+ set<0>(point, X);
+ set<1>(point, Y);
+ range_wrapper.set_z(i, Z);
+ }
+ }
+
+ return ret_errno;
+}
+
+/************************************************************************/
+/* pj_geodetic_to_geocentric() */
+/************************************************************************/
+
+template <typename T, typename Range, bool AddZ>
+inline int pj_geocentric_to_geodetic( T const& a, T const& es,
+ range_wrapper<Range, AddZ> & range_wrapper )
+
+{
+ //typedef typename boost::range_iterator<Range>::type iterator;
+ typedef typename boost::range_value<Range>::type point_type;
+ //typedef typename coordinate_type<point_type>::type coord_t;
+
+ Range & rng = range_wrapper.get_range();
+ std::size_t point_count = boost::size(rng);
+
+ T const b = (es == 0.0) ? a : a * sqrt(1-es);
+
+ GeocentricInfo<T> gi;
+ if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 )
+ {
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+ for( std::size_t i = 0 ; i < point_count ; ++i )
+ {
+ point_type & point = range::at(rng, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ T Longitude = 0, Latitude = 0, Height = 0;
+ pj_Convert_Geocentric_To_Geodetic( gi,
+ get<0>(point),
+ get<1>(point),
+ range_wrapper.get_z(i), // z
+ Longitude, Latitude, Height );
+
+ set_from_radian<0>(point, Longitude);
+ set_from_radian<1>(point, Latitude);
+ range_wrapper.set_z(i, Height); // Height
+ }
+
+ return 0;
+}
+
+/************************************************************************/
+/* pj_compare_datums() */
+/* */
+/* Returns TRUE if the two datums are identical, otherwise */
+/* FALSE. */
+/************************************************************************/
+
+template <typename Par>
+inline bool pj_compare_datums( Par & srcdefn, Par & dstdefn )
+{
+ if( srcdefn.datum_type != dstdefn.datum_type )
+ {
+ return false;
+ }
+ else if( srcdefn.a_orig != dstdefn.a_orig
+ || math::abs(srcdefn.es_orig - dstdefn.es_orig) > 0.000000000050 )
+ {
+ /* the tolerance for es is to ensure that GRS80 and WGS84 are
+ considered identical */
+ return false;
+ }
+ else if( srcdefn.datum_type == PJD_3PARAM )
+ {
+ return (srcdefn.datum_params[0] == dstdefn.datum_params[0]
+ && srcdefn.datum_params[1] == dstdefn.datum_params[1]
+ && srcdefn.datum_params[2] == dstdefn.datum_params[2]);
+ }
+ else if( srcdefn.datum_type == PJD_7PARAM )
+ {
+ return (srcdefn.datum_params[0] == dstdefn.datum_params[0]
+ && srcdefn.datum_params[1] == dstdefn.datum_params[1]
+ && srcdefn.datum_params[2] == dstdefn.datum_params[2]
+ && srcdefn.datum_params[3] == dstdefn.datum_params[3]
+ && srcdefn.datum_params[4] == dstdefn.datum_params[4]
+ && srcdefn.datum_params[5] == dstdefn.datum_params[5]
+ && srcdefn.datum_params[6] == dstdefn.datum_params[6]);
+ }
+ else if( srcdefn.datum_type == PJD_GRIDSHIFT )
+ {
+ return pj_param(srcdefn.params,"snadgrids").s
+ == pj_param(dstdefn.params,"snadgrids").s;
+ }
+ else
+ return true;
+}
+
+/************************************************************************/
+/* pj_geocentic_to_wgs84() */
+/************************************************************************/
+
+template <typename Par, typename Range, bool AddZ>
+inline int pj_geocentric_to_wgs84( Par const& defn,
+ range_wrapper<Range, AddZ> & range_wrapper )
+
+{
+ typedef typename boost::range_value<Range>::type point_type;
+ typedef typename coordinate_type<point_type>::type coord_t;
+
+ Range & rng = range_wrapper.get_range();
+ std::size_t point_count = boost::size(rng);
+
+ if( defn.datum_type == PJD_3PARAM )
+ {
+ for(std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(rng, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ set<0>(point, get<0>(point) + Dx_BF(defn));
+ set<1>(point, get<1>(point) + Dy_BF(defn));
+ range_wrapper.set_z(i, range_wrapper.get_z(i) + Dz_BF(defn));
+ }
+ }
+ else if( defn.datum_type == PJD_7PARAM )
+ {
+ for(std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(rng, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ coord_t x = get<0>(point);
+ coord_t y = get<1>(point);
+ coord_t z = range_wrapper.get_z(i);
+
+ coord_t x_out, y_out, z_out;
+
+ x_out = M_BF(defn)*( x - Rz_BF(defn)*y + Ry_BF(defn)*z) + Dx_BF(defn);
+ y_out = M_BF(defn)*( Rz_BF(defn)*x + y - Rx_BF(defn)*z) + Dy_BF(defn);
+ z_out = M_BF(defn)*(-Ry_BF(defn)*x + Rx_BF(defn)*y + z) + Dz_BF(defn);
+
+ set<0>(point, x_out);
+ set<1>(point, y_out);
+ range_wrapper.set_z(i, z_out);
+ }
+ }
+
+ return 0;
+}
+
+/************************************************************************/
+/* pj_geocentic_from_wgs84() */
+/************************************************************************/
+
+template <typename Par, typename Range, bool AddZ>
+inline int pj_geocentric_from_wgs84( Par const& defn,
+ range_wrapper<Range, AddZ> & range_wrapper )
+
+{
+ typedef typename boost::range_value<Range>::type point_type;
+ typedef typename coordinate_type<point_type>::type coord_t;
+
+ Range & rng = range_wrapper.get_range();
+ std::size_t point_count = boost::size(rng);
+
+ if( defn.datum_type == PJD_3PARAM )
+ {
+ for(std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(rng, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ set<0>(point, get<0>(point) - Dx_BF(defn));
+ set<1>(point, get<1>(point) - Dy_BF(defn));
+ range_wrapper.set_z(i, range_wrapper.get_z(i) - Dz_BF(defn));
+ }
+ }
+ else if( defn.datum_type == PJD_7PARAM )
+ {
+ for(std::size_t i = 0; i < point_count; i++ )
+ {
+ point_type & point = range::at(rng, i);
+
+ if( is_invalid_point(point) )
+ continue;
+
+ coord_t x = get<0>(point);
+ coord_t y = get<1>(point);
+ coord_t z = range_wrapper.get_z(i);
+
+ coord_t x_tmp = (x - Dx_BF(defn)) / M_BF(defn);
+ coord_t y_tmp = (y - Dy_BF(defn)) / M_BF(defn);
+ coord_t z_tmp = (z - Dz_BF(defn)) / M_BF(defn);
+
+ x = x_tmp + Rz_BF(defn)*y_tmp - Ry_BF(defn)*z_tmp;
+ y = -Rz_BF(defn)*x_tmp + y_tmp + Rx_BF(defn)*z_tmp;
+ z = Ry_BF(defn)*x_tmp - Rx_BF(defn)*y_tmp + z_tmp;
+
+ set<0>(point, x);
+ set<1>(point, y);
+ range_wrapper.set_z(i, z);
+ }
+ }
+
+ return 0;
+}
+
+
+inline bool pj_datum_check_error(int err)
+{
+ return err != 0 && (err > 0 || transient_error[-err] == 0);
+}
+
+/************************************************************************/
+/* pj_datum_transform() */
+/* */
+/* The input should be long/lat/z coordinates in radians in the */
+/* source datum, and the output should be long/lat/z */
+/* coordinates in radians in the destination datum. */
+/************************************************************************/
+
+template <typename Par, typename Range>
+inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn,
+ Range & range )
+
+{
+ typedef typename Par::type calc_t;
+ bool result = true;
+
+ calc_t src_a, src_es, dst_a, dst_es;
+
+/* -------------------------------------------------------------------- */
+/* We cannot do any meaningful datum transformation if either */
+/* the source or destination are of an unknown datum type */
+/* (ie. only a +ellps declaration, no +datum). This is new */
+/* behavior for PROJ 4.6.0. */
+/* -------------------------------------------------------------------- */
+ if( srcdefn.datum_type == PJD_UNKNOWN
+ || dstdefn.datum_type == PJD_UNKNOWN )
+ return result;
+
+/* -------------------------------------------------------------------- */
+/* Short cut if the datums are identical. */
+/* -------------------------------------------------------------------- */
+ if( pj_compare_datums( srcdefn, dstdefn ) )
+ return result;
+
+ src_a = srcdefn.a_orig;
+ src_es = srcdefn.es_orig;
+
+ dst_a = dstdefn.a_orig;
+ dst_es = dstdefn.es_orig;
+
+/* -------------------------------------------------------------------- */
+/* Create a temporary Z array if one is not provided. */
+/* -------------------------------------------------------------------- */
+
+ range_wrapper<Range> z_range(range);
+
+/* -------------------------------------------------------------------- */
+/* If this datum requires grid shifts, then apply it to geodetic */
+/* coordinates. */
+/* -------------------------------------------------------------------- */
+ /*if( srcdefn.datum_type == PJD_GRIDSHIFT )
+ {
+ try {
+ pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z );
+ } catch (projection_exception const& e) {
+ if (pj_datum_check_error(e.code())) {
+ BOOST_RETHROW
+ }
+ }
+
+ src_a = SRS_WGS84_SEMIMAJOR;
+ src_es = SRS_WGS84_ESQUARED;
+ }
+
+ if( dstdefn.datum_type == PJD_GRIDSHIFT )
+ {
+ dst_a = SRS_WGS84_SEMIMAJOR;
+ dst_es = SRS_WGS84_ESQUARED;
+ }*/
+
+/* ==================================================================== */
+/* Do we need to go through geocentric coordinates? */
+/* ==================================================================== */
+ if( src_es != dst_es || src_a != dst_a
+ || srcdefn.datum_type == PJD_3PARAM
+ || srcdefn.datum_type == PJD_7PARAM
+ || dstdefn.datum_type == PJD_3PARAM
+ || dstdefn.datum_type == PJD_7PARAM)
+ {
+/* -------------------------------------------------------------------- */
+/* Convert to geocentric coordinates. */
+/* -------------------------------------------------------------------- */
+ int err = pj_geodetic_to_geocentric( src_a, src_es, z_range );
+ if (pj_datum_check_error(err))
+ BOOST_THROW_EXCEPTION( projection_exception(err) );
+ else if (err != 0)
+ result = false;
+
+/* -------------------------------------------------------------------- */
+/* Convert between datums. */
+/* -------------------------------------------------------------------- */
+ if( srcdefn.datum_type == PJD_3PARAM
+ || srcdefn.datum_type == PJD_7PARAM )
+ {
+ try {
+ pj_geocentric_to_wgs84( srcdefn, z_range );
+ } catch (projection_exception const& e) {
+ if (pj_datum_check_error(e.code())) {
+ BOOST_RETHROW
+ }
+ }
+ }
+
+ if( dstdefn.datum_type == PJD_3PARAM
+ || dstdefn.datum_type == PJD_7PARAM )
+ {
+ try {
+ pj_geocentric_from_wgs84( dstdefn, z_range );
+ } catch (projection_exception const& e) {
+ if (pj_datum_check_error(e.code())) {
+ BOOST_RETHROW
+ }
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Convert back to geodetic coordinates. */
+/* -------------------------------------------------------------------- */
+ err = pj_geocentric_to_geodetic( dst_a, dst_es, z_range );
+ if (pj_datum_check_error(err))
+ BOOST_THROW_EXCEPTION( projection_exception(err) );
+ else if (err != 0)
+ result = false;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Apply grid shift to destination if required. */
+/* -------------------------------------------------------------------- */
+ /*if( dstdefn.datum_type == PJD_GRIDSHIFT )
+ {
+ try {
+ pj_apply_gridshift_2( dstdefn, 1, point_count, point_offset, x, y, z );
+ } catch (projection_exception const& e) {
+ if (pj_datum_check_error(e.code()))
+ BOOST_RETHROW
+ }
+ }*/
+
+ return result;
+}
+
+} // namespace detail
+
+}}} // namespace boost::geometry::projections
+
+#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP