diff options
author | DongHun Kwak <dh0128.kwak@samsung.com> | 2019-12-05 15:12:59 +0900 |
---|---|---|
committer | DongHun Kwak <dh0128.kwak@samsung.com> | 2019-12-05 15:12:59 +0900 |
commit | b8cf34c691623e4ec329053cbbf68522a855882d (patch) | |
tree | 34da08632a99677f6b79ecb65e5b655a5b69a67f /boost/geometry/srs/projections/impl/pj_transform.hpp | |
parent | 3fdc3e5ee96dca5b11d1694975a65200787eab86 (diff) | |
download | boost-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.hpp | 1009 |
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 |