// 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/coordinate_dimension.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/srs/projections/impl/geocent.hpp>
#include <boost/geometry/srs/projections/impl/pj_apply_gridshift.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
// -----------------------------------------------------------
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 -56, 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.
**
**
** NOTE (2017-10-01): Non-transient errors really should have resulted in a
** PJ==0 during initialization, and hence should be handled at the level
** before calling pj_transform. The only obvious example of the contrary
** appears to be the PJD_ERR_GRID_AREA case, which may also be taken to
** mean "no grids available"
**
**
*/
static const int transient_error[60] = {
/* 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,
/* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 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,
typename Grids
>
inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn,
DstPrj2 const& dstprj, Par const& dstdefn,
Range & range,
Grids const& srcgrids,
Grids const& dstgrids)
{
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(error_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, srcgrids, dstgrids ) )
{
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(error_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 error_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 = error_lat_or_lon_exceed_limit;
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 error_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 == datum_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 == datum_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 == datum_gridshift )
{
return srcdefn.nadgrids == dstdefn.nadgrids;
}
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 == datum_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 == datum_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 == datum_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 == datum_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, typename Grids>
inline bool pj_datum_transform(Par const& srcdefn,
Par const& dstdefn,
Range & range,
Grids const& srcgrids,
Grids const& dstgrids)
{
typedef typename Par::type calc_t;
// This has to be consistent with default spheroid and pj_ellps
// TODO: Define in one place
static const calc_t wgs84_a = 6378137.0;
static const calc_t wgs84_b = 6356752.3142451793;
static const calc_t wgs84_es = 1. - (wgs84_b * wgs84_b) / (wgs84_a * wgs84_a);
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 == datum_unknown
|| dstdefn.datum_type == datum_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 == datum_gridshift )
{
try {
pj_apply_gridshift_2<false>( srcdefn, range, srcgrids );
} catch (projection_exception const& e) {
if (pj_datum_check_error(e.code())) {
BOOST_RETHROW
}
}
src_a = wgs84_a;
src_es = wgs84_es;
}
if( dstdefn.datum_type == datum_gridshift )
{
dst_a = wgs84_a;
dst_es = wgs84_es;
}
/* ==================================================================== */
/* Do we need to go through geocentric coordinates? */
/* ==================================================================== */
if( src_es != dst_es || src_a != dst_a
|| srcdefn.datum_type == datum_3param
|| srcdefn.datum_type == datum_7param
|| dstdefn.datum_type == datum_3param
|| dstdefn.datum_type == datum_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 == datum_3param
|| srcdefn.datum_type == datum_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 == datum_3param
|| dstdefn.datum_type == datum_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 == datum_gridshift )
{
try {
pj_apply_gridshift_2<true>( dstdefn, range, dstgrids );
} 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