diff options
Diffstat (limited to 'boost/geometry/srs/projections/impl')
27 files changed, 2634 insertions, 594 deletions
diff --git a/boost/geometry/srs/projections/impl/aasincos.hpp b/boost/geometry/srs/projections/impl/aasincos.hpp index 4678029bc0..de064f6a01 100644 --- a/boost/geometry/srs/projections/impl/aasincos.hpp +++ b/boost/geometry/srs/projections/impl/aasincos.hpp @@ -3,6 +3,10 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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) @@ -66,7 +70,7 @@ inline T aasin(T const& v) { if (av > aasincos::ONE_TOL<T>()) { - BOOST_THROW_EXCEPTION( projection_exception(-19) ); + BOOST_THROW_EXCEPTION( projection_exception(error_acos_asin_arg_too_large) ); } return (v < 0.0 ? -geometry::math::half_pi<T>() : geometry::math::half_pi<T>()); } @@ -83,7 +87,7 @@ inline T aacos(T const& v) { if (av > aasincos::ONE_TOL<T>()) { - BOOST_THROW_EXCEPTION( projection_exception(-19) ); + BOOST_THROW_EXCEPTION( projection_exception(error_acos_asin_arg_too_large) ); } return (v < 0.0 ? geometry::math::pi<T>() : 0.0); } diff --git a/boost/geometry/srs/projections/impl/base_dynamic.hpp b/boost/geometry/srs/projections/impl/base_dynamic.hpp index 3979c34300..5eafbe5bea 100644 --- a/boost/geometry/srs/projections/impl/base_dynamic.hpp +++ b/boost/geometry/srs/projections/impl/base_dynamic.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -15,6 +15,7 @@ #include <string> +#include <boost/geometry/srs/projections/exception.hpp> #include <boost/geometry/srs/projections/impl/projects.hpp> namespace boost { namespace geometry { namespace projections diff --git a/boost/geometry/srs/projections/impl/dms_parser.hpp b/boost/geometry/srs/projections/impl/dms_parser.hpp index bbecc9b1a2..21e87d1f67 100644 --- a/boost/geometry/srs/projections/impl/dms_parser.hpp +++ b/boost/geometry/srs/projections/impl/dms_parser.hpp @@ -41,18 +41,12 @@ #include <string> -#include <boost/static_assert.hpp> - -#if !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) -#include <boost/lexical_cast.hpp> -#endif // !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) - #include <boost/algorithm/string.hpp> - #include <boost/config.hpp> +#include <boost/static_assert.hpp> #include <boost/geometry/core/cs.hpp> - +#include <boost/geometry/srs/projections/str_cast.hpp> #include <boost/geometry/util/math.hpp> namespace boost { namespace geometry { namespace projections @@ -129,27 +123,23 @@ struct dms_parser bool has_dms[3]; dms_value() -#ifndef BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX +#if !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && (!defined(_MSC_VER) || (_MSC_VER >= 1900)) // workaround for VC++ 12 (aka 2013) : dms{0, 0, 0} , has_dms{false, false, false} -#endif + {} +#else { -#ifdef BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX std::fill(dms, dms + 3, T(0)); std::fill(has_dms, has_dms + 3, false); -#endif } +#endif }; template <size_t I> static inline void assign_dms(dms_value& dms, std::string& value, bool& has_value) { -#if !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) - dms.dms[I] = boost::lexical_cast<T>(value.c_str()); -#else // !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) - dms.dms[I] = std::atof(value.c_str()); -#endif // !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) + dms.dms[I] = geometry::str_cast<T>(value); dms.has_dms[I] = true; has_value = false; value.clear(); diff --git a/boost/geometry/srs/projections/impl/function_overloads.hpp b/boost/geometry/srs/projections/impl/function_overloads.hpp index 43d33c6bea..5056ecb4db 100644 --- a/boost/geometry/srs/projections/impl/function_overloads.hpp +++ b/boost/geometry/srs/projections/impl/function_overloads.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -26,12 +26,6 @@ inline T atan2(T const& a, T const& b) using std::atan2; return atan2(a, b); } -template <typename T> -inline T pow(T const& a, T const& b) -{ - using std::pow; - return pow(a, b); -} */ template <typename T> diff --git a/boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp b/boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp new file mode 100644 index 0000000000..ee3aa6ed23 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp @@ -0,0 +1,429 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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: +// Author: Frank Warmerdam, warmerdam@pobox.com + +// 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_APPLY_GRIDSHIFT_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP + + +#include <boost/geometry/core/assert.hpp> +#include <boost/geometry/core/radian_access.hpp> + +#include <boost/geometry/srs/projections/impl/pj_gridlist.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +// Originally implemented in nad_intr.c +template <typename CalcT> +inline void nad_intr(CalcT in_lon, CalcT in_lat, + CalcT & out_lon, CalcT & out_lat, + pj_ctable const& ct) +{ + pj_ctable::lp_t frct; + pj_ctable::ilp_t indx; + boost::int32_t in; + + indx.lam = int_floor(in_lon /= ct.del.lam); + indx.phi = int_floor(in_lat /= ct.del.phi); + frct.lam = in_lon - indx.lam; + frct.phi = in_lat - indx.phi; + // TODO: implement differently + out_lon = out_lat = HUGE_VAL; + if (indx.lam < 0) { + if (indx.lam == -1 && frct.lam > 0.99999999999) { + ++indx.lam; + frct.lam = 0.; + } else + return; + } else if ((in = indx.lam + 1) >= ct.lim.lam) { + if (in == ct.lim.lam && frct.lam < 1e-11) { + --indx.lam; + frct.lam = 1.; + } else + return; + } + if (indx.phi < 0) { + if (indx.phi == -1 && frct.phi > 0.99999999999) { + ++indx.phi; + frct.phi = 0.; + } else + return; + } else if ((in = indx.phi + 1) >= ct.lim.phi) { + if (in == ct.lim.phi && frct.phi < 1e-11) { + --indx.phi; + frct.phi = 1.; + } else + return; + } + boost::int32_t index = indx.phi * ct.lim.lam + indx.lam; + pj_ctable::flp_t const& f00 = ct.cvs[index++]; + pj_ctable::flp_t const& f10 = ct.cvs[index]; + index += ct.lim.lam; + pj_ctable::flp_t const& f11 = ct.cvs[index--]; + pj_ctable::flp_t const& f01 = ct.cvs[index]; + CalcT m00, m10, m01, m11; + m11 = m10 = frct.lam; + m00 = m01 = 1. - frct.lam; + m11 *= frct.phi; + m01 *= frct.phi; + frct.phi = 1. - frct.phi; + m00 *= frct.phi; + m10 *= frct.phi; + out_lon = m00 * f00.lam + m10 * f10.lam + + m01 * f01.lam + m11 * f11.lam; + out_lat = m00 * f00.phi + m10 * f10.phi + + m01 * f01.phi + m11 * f11.phi; +} + +// Originally implemented in nad_cvt.c +template <bool Inverse, typename CalcT> +inline void nad_cvt(CalcT const& in_lon, CalcT const& in_lat, + CalcT & out_lon, CalcT & out_lat, + pj_gi const& gi) +{ + static const int max_iterations = 10; + static const CalcT tol = 1e-12; + static const CalcT toltol = tol * tol; + static const CalcT pi = math::pi<CalcT>(); + + // horizontal grid expected + BOOST_GEOMETRY_ASSERT_MSG(gi.format != pj_gi::gtx, + "Vertical grid cannot be used in horizontal shift."); + + pj_ctable const& ct = gi.ct; + + // TODO: implement differently + if (in_lon == HUGE_VAL) + { + out_lon = HUGE_VAL; + out_lat = HUGE_VAL; + return; + } + + // normalize input to ll origin + pj_ctable::lp_t tb; + tb.lam = in_lon - ct.ll.lam; + tb.phi = in_lat - ct.ll.phi; + tb.lam = adjlon (tb.lam - pi) + pi; + + pj_ctable::lp_t t; + nad_intr(tb.lam, tb.phi, t.lam, t.phi, ct); + if (t.lam == HUGE_VAL) + { + out_lon = HUGE_VAL; + out_lat = HUGE_VAL; + return; + } + + if (! Inverse) + { + out_lon = in_lon - t.lam; + out_lat = in_lat - t.phi; + return; + } + + t.lam = tb.lam + t.lam; + t.phi = tb.phi - t.phi; + + int i = max_iterations; + pj_ctable::lp_t del, dif; + do + { + nad_intr(t.lam, t.phi, del.lam, del.phi, ct); + + // This case used to return failure, but I have + // changed it to return the first order approximation + // of the inverse shift. This avoids cases where the + // grid shift *into* this grid came from another grid. + // While we aren't returning optimally correct results + // I feel a close result in this case is better than + // no result. NFW + // To demonstrate use -112.5839956 49.4914451 against + // the NTv2 grid shift file from Canada. + if (del.lam == HUGE_VAL) + { + // Inverse grid shift iteration failed, presumably at grid edge. Using first approximation. + break; + } + + dif.lam = t.lam - del.lam - tb.lam; + dif.phi = t.phi + del.phi - tb.phi; + t.lam -= dif.lam; + t.phi -= dif.phi; + + } + while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); // prob. slightly faster than hypot() + + if (i==0) + { + // Inverse grid shift iterator failed to converge. + out_lon = HUGE_VAL; + out_lat = HUGE_VAL; + return; + } + + out_lon = adjlon (t.lam + ct.ll.lam); + out_lat = t.phi + ct.ll.phi; +} + + +/************************************************************************/ +/* find_grid() */ +/* */ +/* Determine which grid is the correct given an input coordinate. */ +/************************************************************************/ + +// Originally find_ctable() +// here divided into grid_disjoint(), find_grid() and load_grid() + +template <typename T> +inline bool grid_disjoint(T const& lam, T const& phi, + pj_ctable const& ct) +{ + double epsilon = (fabs(ct.del.phi)+fabs(ct.del.lam))/10000.0; + return ct.ll.phi - epsilon > phi + || ct.ll.lam - epsilon > lam + || (ct.ll.phi + (ct.lim.phi-1) * ct.del.phi + epsilon < phi) + || (ct.ll.lam + (ct.lim.lam-1) * ct.del.lam + epsilon < lam); +} + +template <typename T> +inline pj_gi * find_grid(T const& lam, + T const& phi, + std::vector<pj_gi>::iterator first, + std::vector<pj_gi>::iterator last) +{ + pj_gi * gip = NULL; + + for( ; first != last ; ++first ) + { + // skip tables that don't match our point at all. + if (! grid_disjoint(lam, phi, first->ct)) + { + // skip vertical grids + if (first->format != pj_gi::gtx) + { + gip = boost::addressof(*first); + break; + } + } + } + + // If we didn't find a child then nothing more to do + if( gip == NULL ) + return gip; + + // Otherwise use the child, first checking it's children + pj_gi * child = find_grid(lam, phi, first->children.begin(), first->children.end()); + if (child != NULL) + gip = child; + + return gip; +} + +template <typename T> +inline pj_gi * find_grid(T const& lam, + T const& phi, + pj_gridinfo & grids, + std::vector<std::size_t> const& gridindexes) +{ + pj_gi * gip = NULL; + + // keep trying till we find a table that works + for (std::size_t i = 0 ; i < gridindexes.size() ; ++i) + { + pj_gi & gi = grids[gridindexes[i]]; + + // skip tables that don't match our point at all. + if (! grid_disjoint(lam, phi, gi.ct)) + { + // skip vertical grids + if (gi.format != pj_gi::gtx) + { + gip = boost::addressof(gi); + break; + } + } + } + + if (gip == NULL) + return gip; + + // If we have child nodes, check to see if any of them apply. + pj_gi * child = find_grid(lam, phi, gip->children.begin(), gip->children.end()); + if (child != NULL) + gip = child; + + // if we get this far we have found a suitable grid + return gip; +} + + +template <typename StreamPolicy> +inline bool load_grid(StreamPolicy const& stream_policy, pj_gi_load & gi) +{ + // load the grid shift info if we don't have it. + if (gi.ct.cvs.empty()) + { + typename StreamPolicy::stream_type is; + stream_policy.open(is, gi.gridname); + + if (! pj_gridinfo_load(is, gi)) + { + //pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return false; + } + } + + return true; +} + + +/************************************************************************/ +/* pj_apply_gridshift_3() */ +/* */ +/* This is the real workhorse, given a gridlist. */ +/************************************************************************/ + +template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range> +inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy, + Range & range, + srs::grids & grids, + std::vector<std::size_t> const& gridindexes) +{ + typedef typename boost::range_size<Range>::type size_type; + + // If the grids are empty the indexes are as well + if (gridindexes.empty()) + { + //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + //return PJD_ERR_FAILED_TO_LOAD_GRID; + return false; + } + + size_type point_count = boost::size(range); + + for (size_type i = 0 ; i < point_count ; ++i) + { + typename boost::range_reference<Range>::type + point = range::at(range, i); + + CalcT in_lon = geometry::get_as_radian<0>(point); + CalcT in_lat = geometry::get_as_radian<1>(point); + + pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes); + + if ( gip != NULL ) + { + // load the grid shift info if we don't have it. + if (! gip->ct.cvs.empty() || load_grid(stream_policy, *gip)) + { + // TODO: use set_invalid_point() or similar mechanism + CalcT out_lon = HUGE_VAL; + CalcT out_lat = HUGE_VAL; + + nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip); + + // TODO: check differently + if ( out_lon != HUGE_VAL ) + { + geometry::set_from_radian<0>(point, out_lon); + geometry::set_from_radian<1>(point, out_lat); + } + } + } + } + + return true; +} + + +/************************************************************************/ +/* pj_apply_gridshift_2() */ +/* */ +/* This implementation uses the gridlist from a coordinate */ +/* system definition. If the gridlist has not yet been */ +/* populated in the coordinate system definition we set it up */ +/* now. */ +/************************************************************************/ + +template <bool Inverse, typename Par, typename Range, typename ProjGrids> +inline bool pj_apply_gridshift_2(Par const& defn, Range & range, ProjGrids const& grids) +{ + /*if( defn->catalog_name != NULL ) + return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset, + x, y, z );*/ + + /*std::vector<std::size_t> gridindexes; + pj_gridlist_from_nadgrids(pj_get_param_s(defn.params, "nadgrids"), + grids.storage_ptr->stream_policy, + grids.storage_ptr->grids, + gridindexes);*/ + + BOOST_GEOMETRY_ASSERT(grids.storage_ptr != NULL); + + // At this point the grids should be initialized + if (grids.hindexes.empty()) + return false; + + return pj_apply_gridshift_3 + < + Inverse, typename Par::type + >(grids.storage_ptr->stream_policy, + range, + grids.storage_ptr->hgrids, + grids.hindexes); +} + +template <bool Inverse, typename Par, typename Range> +inline bool pj_apply_gridshift_2(Par const& , Range & , srs::detail::empty_projection_grids const& ) +{ + return false; +} + + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP diff --git a/boost/geometry/srs/projections/impl/pj_apply_gridshift_shared.hpp b/boost/geometry/srs/projections/impl/pj_apply_gridshift_shared.hpp new file mode 100644 index 0000000000..2f0d63169e --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_apply_gridshift_shared.hpp @@ -0,0 +1,157 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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: +// Author: Frank Warmerdam, warmerdam@pobox.com + +// 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_APPLY_GRIDSHIFT_SHARED_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_SHARED_HPP + + +#include <boost/geometry/core/assert.hpp> +#include <boost/geometry/core/radian_access.hpp> + +#include <boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp> +#include <boost/geometry/srs/projections/impl/pj_gridlist_shared.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range> +inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy, + Range & range, + srs::shared_grids & grids, + std::vector<std::size_t> const& gridindexes) +{ + typedef typename boost::range_size<Range>::type size_type; + + // If the grids are empty the indexes are as well + if (gridindexes.empty()) + { + //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + //return PJD_ERR_FAILED_TO_LOAD_GRID; + return false; + } + + size_type point_count = boost::size(range); + + // local storage + pj_gi_load local_gi; + + for (size_type i = 0 ; i < point_count ; ) + { + bool load_needed = false; + + CalcT in_lon = 0; + CalcT in_lat = 0; + + { + boost::shared_lock<boost::shared_mutex> lock(grids.mutex); + + for ( ; i < point_count ; ++i ) + { + typename boost::range_reference<Range>::type + point = range::at(range, i); + + in_lon = geometry::get_as_radian<0>(point); + in_lat = geometry::get_as_radian<1>(point); + + pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes); + + if (gip == NULL) + { + // do nothing + } + else if (! gip->ct.cvs.empty()) + { + // TODO: use set_invalid_point() or similar mechanism + CalcT out_lon = HUGE_VAL; + CalcT out_lat = HUGE_VAL; + + nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip); + + // TODO: check differently + if (out_lon != HUGE_VAL) + { + geometry::set_from_radian<0>(point, out_lon); + geometry::set_from_radian<1>(point, out_lat); + } + } + else + { + // loading is needed + local_gi = *gip; + load_needed = true; + break; + } + } + } + + if (load_needed) + { + if (load_grid(stream_policy, local_gi)) + { + boost::unique_lock<boost::shared_mutex> lock(grids.mutex); + + // check again in case other thread already loaded the grid. + pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes); + + if (gip != NULL && gip->ct.cvs.empty()) + { + // swap loaded local storage with empty grid + local_gi.swap(*gip); + } + } + else + { + ++i; + } + } + } + + return true; +} + + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_SHARED_HPP diff --git a/boost/geometry/srs/projections/impl/pj_auth.hpp b/boost/geometry/srs/projections/impl/pj_auth.hpp index 899d0b25b1..8a9e34a09c 100644 --- a/boost/geometry/srs/projections/impl/pj_auth.hpp +++ b/boost/geometry/srs/projections/impl/pj_auth.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -48,14 +48,22 @@ namespace boost { namespace geometry { namespace projections { namespace detail { -static const int APA_SIZE = 3; +template <typename T> +struct apa +{ + static const std::size_t size = 3; + + T const& operator[](size_t i) const { return data[i]; } + T & operator[](size_t i) { return data[i]; } + +private: + T data[3]; +}; /* determine latitude from authalic latitude */ template <typename T> -inline bool pj_authset(T const& es, T* APA) +inline detail::apa<T> pj_authset(T const& es) { - BOOST_GEOMETRY_ASSERT(0 != APA); - static const T P00 = .33333333333333333333; static const T P01 = .17222222222222222222; static const T P02 = .10257936507936507936; @@ -64,29 +72,28 @@ inline bool pj_authset(T const& es, T* APA) static const T P20 = .01641501294219154443; T t = 0; + detail::apa<T> apa; - // if (APA = (double *)pj_malloc(APA_SIZE * sizeof(double))) { - APA[0] = es * P00; + apa[0] = es * P00; t = es * es; - APA[0] += t * P01; - APA[1] = t * P10; + apa[0] += t * P01; + apa[1] = t * P10; t *= es; - APA[0] += t * P02; - APA[1] += t * P11; - APA[2] = t * P20; + apa[0] += t * P02; + apa[1] += t * P11; + apa[2] = t * P20; } - return true; + + return apa; } template <typename T> -inline T pj_authlat(T const& beta, const T* APA) +inline T pj_authlat(T const& beta, detail::apa<T> const& apa) { - BOOST_GEOMETRY_ASSERT(0 != APA); - T const t = beta + beta; - return(beta + APA[0] * sin(t) + APA[1] * sin(t + t) + APA[2] * sin(t + t + t)); + return(beta + apa[0] * sin(t) + apa[1] * sin(t + t) + apa[2] * sin(t + t + t)); } } // namespace detail diff --git a/boost/geometry/srs/projections/impl/pj_datum_set.hpp b/boost/geometry/srs/projections/impl/pj_datum_set.hpp index 5301dc7cfe..622efe3403 100644 --- a/boost/geometry/srs/projections/impl/pj_datum_set.hpp +++ b/boost/geometry/srs/projections/impl/pj_datum_set.hpp @@ -75,7 +75,7 @@ inline void pj_datum_add_defn(BGParams const& , std::vector<pvalue<T> >& pvalues /* definition will last into the pj_ell_set() function called */ /* after this one. */ /* -------------------------------------------------------------------- */ - std::string name = pj_param(pvalues, "sdatum").s; + std::string name = pj_get_param_s(pvalues, "datum"); if(! name.empty()) { /* find the datum definition */ @@ -91,19 +91,19 @@ inline void pj_datum_add_defn(BGParams const& , std::vector<pvalue<T> >& pvalues if (index == -1) { - BOOST_THROW_EXCEPTION( projection_exception(-9) ); + BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) ); } - if(! pj_datums[index].ellipse_id.empty()) + pj_datums_type const& datum = pj_datums[index]; + + if(! datum.ellipse_id.empty()) { - std::string entry("ellps="); - entry +=pj_datums[index].ellipse_id; - pvalues.push_back(pj_mkparam<T>(entry)); + pvalues.push_back(pj_mkparam<T>("ellps", datum.ellipse_id)); } - if(! pj_datums[index].defn.empty()) + if(! datum.defn_n.empty() && ! datum.defn_v.empty()) { - pvalues.push_back(pj_mkparam<T>(pj_datums[index].defn)); + pvalues.push_back(pj_mkparam<T>(datum.defn_n, datum.defn_v)); } } } @@ -129,11 +129,12 @@ inline void pj_datum_add_defn(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAI || ! boost::is_same<typename datum_traits::ellps_type, void>::value; BOOST_MPL_ASSERT_MSG((not_set_or_known), UNKNOWN_DATUM, (bg_parameters_type)); - std::string defn = datum_traits::definition(); + std::string def_n = datum_traits::def_n(); + std::string def_v = datum_traits::def_v(); - if (! defn.empty()) + if (! def_n.empty() && ! def_v.empty()) { - pvalues.push_back(pj_mkparam<T>(defn)); + pvalues.push_back(pj_mkparam<T>(def_n, def_v)); } } @@ -146,21 +147,21 @@ inline void pj_datum_set(BGParams const& bg_params, std::vector<pvalue<T> >& pva { static const T SEC_TO_RAD = detail::SEC_TO_RAD<T>(); - projdef.datum_type = PJD_UNKNOWN; + projdef.datum_type = datum_unknown; pj_datum_add_defn(bg_params, pvalues); /* -------------------------------------------------------------------- */ /* Check for nadgrids parameter. */ /* -------------------------------------------------------------------- */ - std::string nadgrids = pj_param(pvalues, "snadgrids").s; - std::string towgs84 = pj_param(pvalues, "stowgs84").s; + std::string nadgrids = pj_get_param_s(pvalues, "nadgrids"); + std::string towgs84 = pj_get_param_s(pvalues, "towgs84"); if(! nadgrids.empty()) { /* We don't actually save the value separately. It will continue to exist int he param list for use in pj_apply_gridshift.c */ - projdef.datum_type = PJD_GRIDSHIFT; + projdef.datum_type = datum_gridshift; } /* -------------------------------------------------------------------- */ @@ -187,7 +188,7 @@ inline void pj_datum_set(BGParams const& bg_params, std::vector<pvalue<T> >& pva || projdef.datum_params[5] != 0.0 || projdef.datum_params[6] != 0.0 ) { - projdef.datum_type = PJD_7PARAM; + projdef.datum_type = datum_7param; /* transform from arc seconds to radians */ projdef.datum_params[3] *= SEC_TO_RAD; @@ -199,7 +200,7 @@ inline void pj_datum_set(BGParams const& bg_params, std::vector<pvalue<T> >& pva } else { - projdef.datum_type = PJD_3PARAM; + projdef.datum_type = datum_3param; } /* Note that pj_init() will later switch datum_type to diff --git a/boost/geometry/srs/projections/impl/pj_datums.hpp b/boost/geometry/srs/projections/impl/pj_datums.hpp index 55da24a2ca..ab2fcb2357 100644 --- a/boost/geometry/srs/projections/impl/pj_datums.hpp +++ b/boost/geometry/srs/projections/impl/pj_datums.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -41,53 +41,70 @@ #include <boost/geometry/srs/projections/impl/projects.hpp> +#include <string> + namespace boost { namespace geometry { namespace projections { namespace detail { +// Originally defined in projects.h +struct pj_datums_type +{ + std::string id; /* datum keyword */ + std::string defn_n; /* e.g. "to_wgs84" */ + std::string defn_v; /* e.g. "0,0,0" */ + std::string ellipse_id; /* ie from ellipse table */ + std::string comments; /* EPSG code, etc */ +}; + +// Originally defined in projects.h +struct pj_prime_meridians_type +{ + std::string id; /* prime meridian keyword */ + std::string defn; /* offset from greenwich in DMS format. */ +}; + /* * The ellipse code must match one from pj_ellps.c. The datum id should * be kept to 12 characters or less if possible. Use the official OGC * datum name for the comments if available. */ -static const PJ_DATUMS pj_datums[] = +static const pj_datums_type pj_datums[] = { - /* id definition ellipse comments */ - /* -- ---------- ------- -------- */ - {"WGS84", "towgs84=0,0,0", + {"WGS84", "towgs84", "0,0,0", "WGS84", ""}, - {"GGRS87", "towgs84=-199.87,74.79,246.62", + {"GGRS87", "towgs84", "-199.87,74.79,246.62", "GRS80", "Greek_Geodetic_Reference_System_1987"}, - {"NAD83", "towgs84=0,0,0", + {"NAD83", "towgs84", "0,0,0", "GRS80", "North_American_Datum_1983"}, - {"NAD27", "nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat", + {"NAD27", "nadgrids", "@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat", "clrk66", "North_American_Datum_1927"}, - {"potsdam", "towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7", + {"potsdam", "towgs84", "598.1,73.7,418.2,0.202,0.045,-2.455,6.7", "bessel", "Potsdam Rauenberg 1950 DHDN"}, - {"carthage", "towgs84=-263.0,6.0,431.0", + {"carthage", "towgs84", "-263.0,6.0,431.0", "clrk80ign", "Carthage 1934 Tunisia"}, - {"hermannskogel", "towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232", + {"hermannskogel", "towgs84", "577.326,90.129,463.919,5.137,1.474,5.297,2.4232", "bessel", "Hermannskogel"}, - {"ire65", "towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", + {"ire65", "towgs84", "482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", "mod_airy", "Ireland 1965"}, - {"nzgd49", "towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", + {"nzgd49", "towgs84", "59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", "intl", "New Zealand Geodetic Datum 1949"}, - {"OSGB36", "towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", + {"OSGB36", "towgs84", "446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", "airy", "Airy 1830"} }; -static const PJ_PRIME_MERIDIANS pj_prime_meridians[] = +static const pj_prime_meridians_type pj_prime_meridians[] = { /* id definition */ /* -- ---------- */ diff --git a/boost/geometry/srs/projections/impl/pj_ell_set.hpp b/boost/geometry/srs/projections/impl/pj_ell_set.hpp index 6f5f14b780..6e0b817adb 100644 --- a/boost/geometry/srs/projections/impl/pj_ell_set.hpp +++ b/boost/geometry/srs/projections/impl/pj_ell_set.hpp @@ -14,7 +14,9 @@ // 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 -// PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam) +// PROJ4 is converted to Geometry Library by +// Barend Gehrels (Geodan, Amsterdam) +// Adam Wulkiewicz // Original copyright notice: @@ -78,12 +80,12 @@ inline void pj_ell_set(BGParams const& /*bg_params*/, std::vector<pvalue<T> >& p a = es = 0.; /* R takes precedence */ - if (pj_param(parameters, "tR").i) - a = pj_param(parameters, "dR").f; - else { /* probable elliptical figure */ + if (pj_param_f(parameters, "R", a)) { + /* empty */ + } else { /* probable elliptical figure */ /* check if ellps present and temporarily append its values to pl */ - name = pj_param(parameters, "sellps").s; + name = pj_get_param_s(parameters, "ellps"); if (! name.empty()) { const int n = sizeof(pj_ellps) / sizeof(pj_ellps[0]); @@ -97,59 +99,56 @@ inline void pj_ell_set(BGParams const& /*bg_params*/, std::vector<pvalue<T> >& p } if (index == -1) { - BOOST_THROW_EXCEPTION( projection_exception(-9) ); + BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) ); } - parameters.push_back(pj_mkparam<T>(pj_ellps[index].major)); - parameters.push_back(pj_mkparam<T>(pj_ellps[index].ell)); + pj_ellps_type const& pj_ellp = pj_ellps[index]; + parameters.push_back(pj_mkparam<T>("a", pj_ellp.major_v)); + parameters.push_back(pj_mkparam<T>(pj_ellp.ell_n, pj_ellp.ell_v)); } - a = pj_param(parameters, "da").f; - if (pj_param(parameters, "tes").i) /* eccentricity squared */ - es = pj_param(parameters, "des").f; - else if (pj_param(parameters, "te").i) { /* eccentricity */ - e = pj_param(parameters, "de").f; + a = pj_get_param_f(parameters, "a"); + if (pj_param_f(parameters, "es", es)) {/* eccentricity squared */ + /* empty */ + } else if (pj_param_f(parameters, "e", e)) { /* eccentricity */ es = e * e; - } else if (pj_param(parameters, "trf").i) { /* recip flattening */ - es = pj_param(parameters, "drf").f; + } else if (pj_param_f(parameters, "rf", es)) { /* recip flattening */ if (!es) { - BOOST_THROW_EXCEPTION( projection_exception(-10) ); + BOOST_THROW_EXCEPTION( projection_exception(error_rev_flattening_is_zero) ); } es = 1./ es; es = es * (2. - es); - } else if (pj_param(parameters, "tf").i) { /* flattening */ - es = pj_param(parameters, "df").f; + } else if (pj_param_f(parameters, "f", es)) { /* flattening */ es = es * (2. - es); - } else if (pj_param(parameters, "tb").i) { /* minor axis */ - b = pj_param(parameters, "db").f; + } else if (pj_param_f(parameters, "b", b)) { /* minor axis */ es = 1. - (b * b) / (a * a); } /* else es == 0. and sphere of radius a */ if (!b) b = a * sqrt(1. - es); /* following options turn ellipsoid into equivalent sphere */ - if (pj_param(parameters, "bR_A").i) { /* sphere--area of ellipsoid */ + if (pj_get_param_b(parameters, "R_A")) { /* sphere--area of ellipsoid */ a *= 1. - es * (SIXTH<T>() + es * (RA4<T>() + es * RA6<T>())); es = 0.; - } else if (pj_param(parameters, "bR_V").i) { /* sphere--vol. of ellipsoid */ + } else if (pj_get_param_b(parameters, "R_V")) { /* sphere--vol. of ellipsoid */ a *= 1. - es * (SIXTH<T>() + es * (RV4<T>() + es * RV6<T>())); es = 0.; - } else if (pj_param(parameters, "bR_a").i) { /* sphere--arithmetic mean */ + } else if (pj_get_param_b(parameters, "R_a")) { /* sphere--arithmetic mean */ a = .5 * (a + b); es = 0.; - } else if (pj_param(parameters, "bR_g").i) { /* sphere--geometric mean */ + } else if (pj_get_param_b(parameters, "R_g")) { /* sphere--geometric mean */ a = sqrt(a * b); es = 0.; - } else if (pj_param(parameters, "bR_h").i) { /* sphere--harmonic mean */ + } else if (pj_get_param_b(parameters, "R_h")) { /* sphere--harmonic mean */ a = 2. * a * b / (a + b); es = 0.; } else { - int i = pj_param(parameters, "tR_lat_a").i; + T tmp; + int i = pj_param_r(parameters, "R_lat_a", tmp); if (i || /* sphere--arith. */ - pj_param(parameters, "tR_lat_g").i) { /* or geom. mean at latitude */ - T tmp; + pj_param_r(parameters, "R_lat_g", tmp)) { /* or geom. mean at latitude */ - tmp = sin(pj_param(parameters, i ? "rR_lat_a" : "rR_lat_g").f); + tmp = sin(tmp); if (geometry::math::abs(tmp) > geometry::math::half_pi<T>()) { - BOOST_THROW_EXCEPTION( projection_exception(-11) ); + BOOST_THROW_EXCEPTION( projection_exception(error_ref_rad_larger_than_90) ); } tmp = 1. - es * tmp * tmp; a *= i ? .5 * (1. - es + tmp) / ( tmp * sqrt(tmp)) : @@ -161,10 +160,10 @@ inline void pj_ell_set(BGParams const& /*bg_params*/, std::vector<pvalue<T> >& p /* some remaining checks */ if (es < 0.) { - BOOST_THROW_EXCEPTION( projection_exception(-12) ); + BOOST_THROW_EXCEPTION( projection_exception(error_es_less_than_zero) ); } if (a <= 0.) { - BOOST_THROW_EXCEPTION( projection_exception(-13) ); + BOOST_THROW_EXCEPTION( projection_exception(error_major_axis_not_given) ); } } @@ -192,13 +191,85 @@ inline void pj_ell_set(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> c /* some remaining checks */ if (es < 0.) { - BOOST_THROW_EXCEPTION( projection_exception(-12) ); + BOOST_THROW_EXCEPTION( projection_exception(error_es_less_than_zero) ); } if (a <= 0.) { - BOOST_THROW_EXCEPTION( projection_exception(-13) ); + BOOST_THROW_EXCEPTION( projection_exception(error_major_axis_not_given) ); } } +template <typename T> +inline void pj_calc_ellipsoid_params(parameters<T> & p, T const& a, T const& es) { +/**************************************************************************************** + Calculate a large number of ancillary ellipsoidal parameters, in addition to + the two traditional PROJ defining parameters: Semimajor axis, a, and the + eccentricity squared, es. + + Most of these parameters are fairly cheap to compute in comparison to the overall + effort involved in initializing a PJ object. They may, however, take a substantial + part of the time taken in computing an individual point transformation. + + So by providing them up front, we can amortize the (already modest) cost over all + transformations carried out over the entire lifetime of a PJ object, rather than + incur that cost for every single transformation. + + Most of the parameter calculations here are based on the "angular eccentricity", + i.e. the angle, measured from the semiminor axis, of a line going from the north + pole to one of the foci of the ellipsoid - or in other words: The arc sine of the + eccentricity. + + The formulae used are mostly taken from: + + Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991). + Columbus, Ohio: Dept. of Geodetic Science + and Surveying, Ohio State University. + +****************************************************************************************/ + + p.a = a; + p.es = es; + + /* Compute some ancillary ellipsoidal parameters */ + if (p.e==0) + p.e = sqrt(p.es); /* eccentricity */ + //p.alpha = asin (p.e); /* angular eccentricity */ + + /* second eccentricity */ + //p.e2 = tan (p.alpha); + //p.e2s = p.e2 * p.e2; + + /* third eccentricity */ + //p.e3 = (0!=p.alpha)? sin (p.alpha) / sqrt(2 - sin (p.alpha)*sin (p.alpha)): 0; + //p.e3s = p.e3 * p.e3; + + /* flattening */ + //if (0==p.f) + // p.f = 1 - cos (p.alpha); /* = 1 - sqrt (1 - PIN->es); */ + //p.rf = p.f != 0.0 ? 1.0/p.f: HUGE_VAL; + + /* second flattening */ + //p.f2 = (cos(p.alpha)!=0)? 1/cos (p.alpha) - 1: 0; + //p.rf2 = p.f2 != 0.0 ? 1/p.f2: HUGE_VAL; + + /* third flattening */ + //p.n = pow (tan (p.alpha/2), 2); + //p.rn = p.n != 0.0 ? 1/p.n: HUGE_VAL; + + /* ...and a few more */ + //if (0==p.b) + // p.b = (1 - p.f)*p.a; + //p.rb = 1. / p.b; + p.ra = 1. / p.a; + + p.one_es = 1. - p.es; + if (p.one_es == 0.) { + BOOST_THROW_EXCEPTION( projection_exception(error_eccentricity_is_one) ); + } + + p.rone_es = 1./p.one_es; +} + + } // namespace detail }}} // namespace boost::geometry::projections diff --git a/boost/geometry/srs/projections/impl/pj_ellps.hpp b/boost/geometry/srs/projections/impl/pj_ellps.hpp index 586802778c..7a10b7de67 100644 --- a/boost/geometry/srs/projections/impl/pj_ellps.hpp +++ b/boost/geometry/srs/projections/impl/pj_ellps.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -41,55 +41,67 @@ #include <boost/geometry/srs/projections/impl/projects.hpp> +#include <string> + namespace boost { namespace geometry { namespace projections { namespace detail { -static const PJ_ELLPS pj_ellps[] = +// Originally defined in projects.h +struct pj_ellps_type +{ + std::string id; /* ellipse keyword name */ + std::string major_v; /* a's value */ + std::string ell_n; /* elliptical parameter name */ + std::string ell_v; /* elliptical parameter value */ + std::string name; /* comments */ +}; + +static const pj_ellps_type pj_ellps[] = { - {"MERIT", "a=6378137.0", "rf=298.257", "MERIT 1983"}, - {"SGS85", "a=6378136.0", "rf=298.257", "Soviet Geodetic System 85"}, - {"GRS80", "a=6378137.0", "rf=298.257222101", "GRS 1980(IUGG, 1980)"}, - {"IAU76", "a=6378140.0", "rf=298.257", "IAU 1976"}, - {"airy", "a=6377563.396", "b=6356256.910", "Airy 1830"}, - {"APL4.9", "a=6378137.0.", "rf=298.25", "Appl. Physics. 1965"}, - {"NWL9D", "a=6378145.0.", "rf=298.25", "Naval Weapons Lab., 1965"}, - {"mod_airy", "a=6377340.189", "b=6356034.446", "Modified Airy"}, - {"andrae", "a=6377104.43", "rf=300.0", "Andrae 1876 (Den., Iclnd.)"}, - {"aust_SA", "a=6378160.0", "rf=298.25", "Australian Natl & S. Amer. 1969"}, - {"GRS67", "a=6378160.0", "rf=298.2471674270", "GRS 67(IUGG 1967)"}, - {"bessel", "a=6377397.155", "rf=299.1528128", "Bessel 1841"}, - {"bess_nam", "a=6377483.865", "rf=299.1528128", "Bessel 1841 (Namibia)"}, - {"clrk66", "a=6378206.4", "b=6356583.8", "Clarke 1866"}, - {"clrk80", "a=6378249.145", "rf=293.4663", "Clarke 1880 mod."}, - {"clrk80ign", "a=6378249.2", "rf=293.4660212936269", "Clarke 1880 (IGN)."}, - {"CPM", "a=6375738.7", "rf=334.29", "Comm. des Poids et Mesures 1799"}, - {"delmbr", "a=6376428.", "rf=311.5", "Delambre 1810 (Belgium)"}, - {"engelis", "a=6378136.05", "rf=298.2566", "Engelis 1985"}, - {"evrst30", "a=6377276.345", "rf=300.8017", "Everest 1830"}, - {"evrst48", "a=6377304.063", "rf=300.8017", "Everest 1948"}, - {"evrst56", "a=6377301.243", "rf=300.8017", "Everest 1956"}, - {"evrst69", "a=6377295.664", "rf=300.8017", "Everest 1969"}, - {"evrstSS", "a=6377298.556", "rf=300.8017", "Everest (Sabah & Sarawak)"}, - {"fschr60", "a=6378166.", "rf=298.3", "Fischer (Mercury Datum) 1960"}, - {"fschr60m", "a=6378155.", "rf=298.3", "Modified Fischer 1960"}, - {"fschr68", "a=6378150.", "rf=298.3", "Fischer 1968"}, - {"helmert", "a=6378200.", "rf=298.3", "Helmert 1906"}, - {"hough", "a=6378270.0", "rf=297.", "Hough"}, - {"intl", "a=6378388.0", "rf=297.", "International 1909 (Hayford)"}, - {"krass", "a=6378245.0", "rf=298.3", "Krassovsky, 1942"}, - {"kaula", "a=6378163.", "rf=298.24", "Kaula 1961"}, - {"lerch", "a=6378139.", "rf=298.257", "Lerch 1979"}, - {"mprts", "a=6397300.", "rf=191.", "Maupertius 1738"}, - {"new_intl", "a=6378157.5", "b=6356772.2", "New International 1967"}, - {"plessis", "a=6376523.", "b=6355863.", "Plessis 1817 (France)"}, - {"SEasia", "a=6378155.0", "b=6356773.3205", "Southeast Asia"}, - {"walbeck", "a=6376896.0", "b=6355834.8467", "Walbeck"}, - {"WGS60", "a=6378165.0", "rf=298.3", "WGS 60"}, - {"WGS66", "a=6378145.0", "rf=298.25", "WGS 66"}, - {"WGS72", "a=6378135.0", "rf=298.26", "WGS 72"}, - {"WGS84", "a=6378137.0", "rf=298.257223563", "WGS 84"}, - {"sphere", "a=6370997.0", "b=6370997.0", "Normal Sphere (r=6370997)"} + {"MERIT", "6378137.0", "rf", "298.257", "MERIT 1983"}, + {"SGS85", "6378136.0", "rf", "298.257", "Soviet Geodetic System 85"}, + {"GRS80", "6378137.0", "rf", "298.257222101", "GRS 1980(IUGG, 1980)"}, + {"IAU76", "6378140.0", "rf", "298.257", "IAU 1976"}, + {"airy", "6377563.396", "b", "6356256.910", "Airy 1830"}, + {"APL4.9", "6378137.0.", "rf", "298.25", "Appl. Physics. 1965"}, + {"NWL9D", "6378145.0.", "rf", "298.25", "Naval Weapons Lab., 1965"}, + {"mod_airy", "6377340.189", "b", "6356034.446", "Modified Airy"}, + {"andrae", "6377104.43", "rf", "300.0", "Andrae 1876 (Den., Iclnd.)"}, + {"aust_SA", "6378160.0", "rf", "298.25", "Australian Natl & S. Amer. 1969"}, + {"GRS67", "6378160.0", "rf", "298.2471674270", "GRS 67(IUGG 1967)"}, + {"bessel", "6377397.155", "rf", "299.1528128", "Bessel 1841"}, + {"bess_nam", "6377483.865", "rf", "299.1528128", "Bessel 1841 (Namibia)"}, + {"clrk66", "6378206.4", "b", "6356583.8", "Clarke 1866"}, + {"clrk80", "6378249.145", "rf", "293.4663", "Clarke 1880 mod."}, + {"clrk80ign", "6378249.2", "rf", "293.4660212936269", "Clarke 1880 (IGN)."}, + {"CPM", "6375738.7", "rf", "334.29", "Comm. des Poids et Mesures 1799"}, + {"delmbr", "6376428.", "rf", "311.5", "Delambre 1810 (Belgium)"}, + {"engelis", "6378136.05", "rf", "298.2566", "Engelis 1985"}, + {"evrst30", "6377276.345", "rf", "300.8017", "Everest 1830"}, + {"evrst48", "6377304.063", "rf", "300.8017", "Everest 1948"}, + {"evrst56", "6377301.243", "rf", "300.8017", "Everest 1956"}, + {"evrst69", "6377295.664", "rf", "300.8017", "Everest 1969"}, + {"evrstSS", "6377298.556", "rf", "300.8017", "Everest (Sabah & Sarawak)"}, + {"fschr60", "6378166.", "rf", "298.3", "Fischer (Mercury Datum) 1960"}, + {"fschr60m", "6378155.", "rf", "298.3", "Modified Fischer 1960"}, + {"fschr68", "6378150.", "rf", "298.3", "Fischer 1968"}, + {"helmert", "6378200.", "rf", "298.3", "Helmert 1906"}, + {"hough", "6378270.0", "rf", "297.", "Hough"}, + {"intl", "6378388.0", "rf", "297.", "International 1909 (Hayford)"}, + {"krass", "6378245.0", "rf", "298.3", "Krassovsky, 1942"}, + {"kaula", "6378163.", "rf", "298.24", "Kaula 1961"}, + {"lerch", "6378139.", "rf", "298.257", "Lerch 1979"}, + {"mprts", "6397300.", "rf", "191.", "Maupertius 1738"}, + {"new_intl", "6378157.5", "b", "6356772.2", "New International 1967"}, + {"plessis", "6376523.", "b", "6355863.", "Plessis 1817 (France)"}, + {"SEasia", "6378155.0", "b", "6356773.3205", "Southeast Asia"}, + {"walbeck", "6376896.0", "b", "6355834.8467", "Walbeck"}, + {"WGS60", "6378165.0", "rf", "298.3", "WGS 60"}, + {"WGS66", "6378145.0", "rf", "298.25", "WGS 66"}, + {"WGS72", "6378135.0", "rf", "298.26", "WGS 72"}, + {"WGS84", "6378137.0", "rf", "298.257223563", "WGS 84"}, + {"sphere", "6370997.0", "b", "6370997.0", "Normal Sphere (r=6370997)"} }; } // namespace detail diff --git a/boost/geometry/srs/projections/impl/pj_fwd.hpp b/boost/geometry/srs/projections/impl/pj_fwd.hpp index 73101b7f40..c1725b0a81 100644 --- a/boost/geometry/srs/projections/impl/pj_fwd.hpp +++ b/boost/geometry/srs/projections/impl/pj_fwd.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -69,7 +69,7 @@ inline void pj_fwd(Prj const& prj, P const& par, LL const& ll, XY& xy) /* check for forward and latitude or longitude overange */ if (t > EPS || geometry::math::abs(lp_lon) > 10.) { - BOOST_THROW_EXCEPTION( projection_exception(-14) ); + BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) ); } if (geometry::math::abs(t) <= EPS) diff --git a/boost/geometry/srs/projections/impl/pj_gauss.hpp b/boost/geometry/srs/projections/impl/pj_gauss.hpp index 2c90870434..94b8f89862 100644 --- a/boost/geometry/srs/projections/impl/pj_gauss.hpp +++ b/boost/geometry/srs/projections/impl/pj_gauss.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -45,13 +45,11 @@ namespace boost { namespace geometry { namespace projections { -namespace detail { namespace gauss { +namespace detail { -static const int MAX_ITER = 20; - template <typename T> -struct GAUSS +struct gauss { T C; T K; @@ -62,13 +60,13 @@ struct GAUSS template <typename T> inline T srat(T const& esinp, T const& exp) { - return (pow((1.0 - esinp) / (1.0 + esinp), exp)); + return (math::pow((T(1) - esinp) / (T(1) + esinp), exp)); } template <typename T> -inline GAUSS<T> gauss_ini(T const& e, T const& phi0, T& chi, T& rc) +inline gauss<T> gauss_ini(T const& e, T const& phi0, T& chi, T& rc) { - static const T FORTPI = detail::FORTPI<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); using std::asin; using std::cos; @@ -80,7 +78,7 @@ inline GAUSS<T> gauss_ini(T const& e, T const& phi0, T& chi, T& rc) T cphi = 0; T es = 0; - GAUSS<T> en; + gauss<T> en; es = e * e; en.e = e; sphi = sin(phi0); @@ -91,38 +89,41 @@ inline GAUSS<T> gauss_ini(T const& e, T const& phi0, T& chi, T& rc) en.C = sqrt(1.0 + es * cphi * cphi / (1.0 - es)); chi = asin(sphi / en.C); en.ratexp = 0.5 * en.C * e; - en.K = tan(0.5 * chi + FORTPI) - / (pow(tan(0.5 * phi0 + FORTPI), en.C) * srat(en.e * sphi, en.ratexp)); + en.K = tan(0.5 * chi + fourth_pi) + / (math::pow(tan(T(0.5) * phi0 + fourth_pi), en.C) * srat(en.e * sphi, en.ratexp)); return en; } template <typename T> -inline void gauss(GAUSS<T> const& en, T& lam, T& phi) +inline void gauss_fwd(gauss<T> const& en, T& lam, T& phi) { - static const T FORTPI = detail::FORTPI<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); + static const T half_pi = detail::half_pi<T>(); - phi = 2.0 * atan(en.K * pow(tan(0.5 * phi + FORTPI), en.C) - * srat(en.e * sin(phi), en.ratexp) ) - geometry::math::half_pi<T>(); + phi = T(2) * atan(en.K * math::pow(tan(T(0.5) * phi + fourth_pi), en.C) + * srat(en.e * sin(phi), en.ratexp) ) - half_pi; lam *= en.C; } template <typename T> -inline void inv_gauss(GAUSS<T> const& en, T& lam, T& phi) +inline void gauss_inv(gauss<T> const& en, T& lam, T& phi) { - static const T FORTPI = detail::FORTPI<T>(); - static const T DEL_TOL = 1e-14; + static const int max_iter = 20; + static const T fourth_pi = detail::fourth_pi<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T del_tol = 1e-14; lam /= en.C; - const T num = pow(tan(0.5 * phi + FORTPI) / en.K, 1.0 / en.C); + const T num = math::pow(tan(T(0.5) * phi + fourth_pi) / en.K, T(1) / en.C); int i = 0; - for (i = MAX_ITER; i; --i) + for (i = max_iter; i; --i) { - const T elp_phi = 2.0 * atan(num * srat(en.e * sin(phi), - 0.5 * en.e)) - geometry::math::half_pi<T>(); + const T elp_phi = 2.0 * atan(num * srat(en.e * sin(phi), - 0.5 * en.e)) - half_pi; - if (geometry::math::abs(elp_phi - phi) < DEL_TOL) + if (geometry::math::abs(elp_phi - phi) < del_tol) { break; } @@ -132,11 +133,12 @@ inline void inv_gauss(GAUSS<T> const& en, T& lam, T& phi) /* convergence failed */ if (!i) { - BOOST_THROW_EXCEPTION( projection_exception(-17) ); + BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) ); } } -}} // namespace detail::gauss +} // namespace detail + }}} // namespace boost::geometry::projections #endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GAUSS_HPP diff --git a/boost/geometry/srs/projections/impl/pj_gridinfo.hpp b/boost/geometry/srs/projections/impl/pj_gridinfo.hpp new file mode 100644 index 0000000000..af1fe471d8 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_gridinfo.hpp @@ -0,0 +1,960 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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: +// Author: Frank Warmerdam, warmerdam@pobox.com + +// 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_GRIDINFO_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP + + +#include <boost/algorithm/string.hpp> + +#include <boost/geometry/core/assert.hpp> + +#include <boost/cstdint.hpp> + +#include <algorithm> +#include <string> +#include <vector> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +/************************************************************************/ +/* swap_words() */ +/* */ +/* Convert the byte order of the given word(s) in place. */ +/************************************************************************/ + +inline bool is_lsb() +{ + static const int byte_order_test = 1; + static bool result = (1 == ((const unsigned char *) (&byte_order_test))[0]); + return result; +} + +inline void swap_words( char *data, int word_size, int word_count ) +{ + for (int word = 0; word < word_count; word++) + { + for (int i = 0; i < word_size/2; i++) + { + std::swap(data[i], data[word_size-i-1]); + } + + data += word_size; + } +} + +inline bool cstr_equal(const char * s1, const char * s2, std::size_t n) +{ + return std::equal(s1, s1 + n, s2); +} + +struct is_trimmable_char +{ + inline bool operator()(char ch) + { + return ch == '\n' || ch == ' '; + } +}; + +// structs originally defined in projects.h + +struct pj_ctable +{ + struct lp_t { double lam, phi; }; + struct flp_t { float lam, phi; }; + struct ilp_t { boost::int32_t lam, phi; }; + + std::string id; // ascii info + lp_t ll; // lower left corner coordinates + lp_t del; // size of cells + ilp_t lim; // limits of conversion matrix + std::vector<flp_t> cvs; // conversion matrix + + inline void swap(pj_ctable & r) + { + id.swap(r.id); + std::swap(ll, r.ll); + std::swap(del, r.del); + std::swap(lim, r.lim); + cvs.swap(r.cvs); + } +}; + +struct pj_gi_load +{ + enum format_t { missing = 0, ntv1, ntv2, gtx, ctable, ctable2 }; + typedef boost::long_long_type offset_t; + + explicit pj_gi_load(std::string const& gname = "", + format_t f = missing, + offset_t off = 0, + bool swap = false) + : gridname(gname) + , format(f) + , grid_offset(off) + , must_swap(swap) + {} + + std::string gridname; // identifying name of grid, eg "conus" or ntv2_0.gsb + + format_t format; // format of this grid, ie "ctable", "ntv1", "ntv2" or "missing". + + offset_t grid_offset; // offset in file, for delayed loading + bool must_swap; // only for NTv2 + + pj_ctable ct; + + inline void swap(pj_gi_load & r) + { + gridname.swap(r.gridname); + std::swap(format, r.format); + std::swap(grid_offset, r.grid_offset); + std::swap(must_swap, r.must_swap); + ct.swap(r.ct); + } + +}; + +struct pj_gi + : pj_gi_load +{ + explicit pj_gi(std::string const& gname = "", + pj_gi_load::format_t f = missing, + pj_gi_load::offset_t off = 0, + bool swap = false) + : pj_gi_load(gname, f, off, swap) + {} + + std::vector<pj_gi> children; + + inline void swap(pj_gi & r) + { + pj_gi_load::swap(r); + children.swap(r.children); + } +}; + +typedef std::vector<pj_gi> pj_gridinfo; + + +/************************************************************************/ +/* pj_gridinfo_load_ctable() */ +/* */ +/* Load the data portion of a ctable formatted grid. */ +/************************************************************************/ + +// Originally nad_ctable_load() defined in nad_init.c +template <typename IStream> +bool pj_gridinfo_load_ctable(IStream & is, pj_gi_load & gi) +{ + pj_ctable & ct = gi.ct; + + // Move the input stream by the size of the proj4 original CTABLE + std::size_t header_size = 80 + + 2 * sizeof(pj_ctable::lp_t) + + sizeof(pj_ctable::ilp_t) + + sizeof(pj_ctable::flp_t*); + is.seekg(header_size); + + // read all the actual shift values + std::size_t a_size = ct.lim.lam * ct.lim.phi; + ct.cvs.resize(a_size); + + std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size; + is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size); + + if (is.fail() || is.gcount() != ch_size) + { + ct.cvs.clear(); + //ctable loading failed on fread() - binary incompatible? + return false; + } + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_load_ctable2() */ +/* */ +/* Load the data portion of a ctable2 formatted grid. */ +/************************************************************************/ + +// Originally nad_ctable2_load() defined in nad_init.c +template <typename IStream> +bool pj_gridinfo_load_ctable2(IStream & is, pj_gi_load & gi) +{ + pj_ctable & ct = gi.ct; + + is.seekg(160); + + // read all the actual shift values + std::size_t a_size = ct.lim.lam * ct.lim.phi; + ct.cvs.resize(a_size); + + std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size; + is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size); + + if (is.fail() || is.gcount() != ch_size) + { + //ctable2 loading failed on fread() - binary incompatible? + ct.cvs.clear(); + return false; + } + + if (! is_lsb()) + { + swap_words(reinterpret_cast<char*>(&ct.cvs[0]), 4, (int)a_size * 2); + } + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_load_ntv1() */ +/* */ +/* NTv1 format. */ +/* We process one line at a time. Note that the array storage */ +/* direction (e-w) is different in the NTv1 file and what */ +/* the CTABLE is supposed to have. The phi/lam are also */ +/* reversed, and we have to be aware of byte swapping. */ +/************************************************************************/ + +// originally in pj_gridinfo_load() function +template <typename IStream> +inline bool pj_gridinfo_load_ntv1(IStream & is, pj_gi_load & gi) +{ + static const double s2r = math::d2r<double>() / 3600.0; + + std::size_t const r_size = gi.ct.lim.lam * 2; + std::size_t const ch_size = sizeof(double) * r_size; + + is.seekg(gi.grid_offset); + + std::vector<double> row_buf(r_size); + gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi); + + for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ ) + { + is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size); + + if (is.fail() || is.gcount() != ch_size) + { + gi.ct.cvs.clear(); + return false; + } + + if (is_lsb()) + swap_words(reinterpret_cast<char*>(&row_buf[0]), 8, (int)r_size); + + // convert seconds to radians + for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ ) + { + pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)]; + + cvs.phi = (float) (row_buf[i*2] * s2r); + cvs.lam = (float) (row_buf[i*2+1] * s2r); + } + } + + return true; +} + +/* -------------------------------------------------------------------- */ +/* pj_gridinfo_load_ntv2() */ +/* */ +/* NTv2 format. */ +/* We process one line at a time. Note that the array storage */ +/* direction (e-w) is different in the NTv2 file and what */ +/* the CTABLE is supposed to have. The phi/lam are also */ +/* reversed, and we have to be aware of byte swapping. */ +/* -------------------------------------------------------------------- */ + +// originally in pj_gridinfo_load() function +template <typename IStream> +inline bool pj_gridinfo_load_ntv2(IStream & is, pj_gi_load & gi) +{ + static const double s2r = math::d2r<double>() / 3600.0; + + std::size_t const r_size = gi.ct.lim.lam * 4; + std::size_t const ch_size = sizeof(float) * r_size; + + is.seekg(gi.grid_offset); + + std::vector<float> row_buf(r_size); + gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi); + + for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ ) + { + is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size); + + if (is.fail() || is.gcount() != ch_size) + { + gi.ct.cvs.clear(); + return false; + } + + if (gi.must_swap) + { + swap_words(reinterpret_cast<char*>(&row_buf[0]), 4, (int)r_size); + } + + // convert seconds to radians + for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ ) + { + pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)]; + + // skip accuracy values + cvs.phi = (float) (row_buf[i*4] * s2r); + cvs.lam = (float) (row_buf[i*4+1] * s2r); + } + } + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_load_gtx() */ +/* */ +/* GTX format. */ +/************************************************************************/ + +// originally in pj_gridinfo_load() function +template <typename IStream> +inline bool pj_gridinfo_load_gtx(IStream & is, pj_gi_load & gi) +{ + boost::int32_t words = gi.ct.lim.lam * gi.ct.lim.phi; + std::size_t const ch_size = sizeof(float) * words; + + is.seekg(gi.grid_offset); + + // TODO: Consider changing this unintuitive code + // NOTE: Vertical shift data (one float per point) is stored in a container + // holding horizontal shift data (two floats per point). + gi.ct.cvs.resize((words + 1) / 2); + + is.read(reinterpret_cast<char*>(&gi.ct.cvs[0]), ch_size); + + if (is.fail() || is.gcount() != ch_size) + { + gi.ct.cvs.clear(); + return false; + } + + if (is_lsb()) + { + swap_words(reinterpret_cast<char*>(&gi.ct.cvs[0]), 4, words); + } + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_load() */ +/* */ +/* This function is intended to implement delayed loading of */ +/* the data contents of a grid file. The header and related */ +/* stuff are loaded by pj_gridinfo_init(). */ +/************************************************************************/ + +template <typename IStream> +inline bool pj_gridinfo_load(IStream & is, pj_gi_load & gi) +{ + if (! gi.ct.cvs.empty()) + { + return true; + } + + if (! is.is_open()) + { + return false; + } + + // Original platform specific CTable format. + if (gi.format == pj_gi::ctable) + { + return pj_gridinfo_load_ctable(is, gi); + } + // CTable2 format. + else if (gi.format == pj_gi::ctable2) + { + return pj_gridinfo_load_ctable2(is, gi); + } + // NTv1 format. + else if (gi.format == pj_gi::ntv1) + { + return pj_gridinfo_load_ntv1(is, gi); + } + // NTv2 format. + else if (gi.format == pj_gi::ntv2) + { + return pj_gridinfo_load_ntv2(is, gi); + } + // GTX format. + else if (gi.format == pj_gi::gtx) + { + return pj_gridinfo_load_gtx(is, gi); + } + else + { + return false; + } +} + +/************************************************************************/ +/* pj_gridinfo_parent() */ +/* */ +/* Seek a parent grid file by name from a grid list */ +/************************************************************************/ + +template <typename It> +inline It pj_gridinfo_parent(It first, It last, std::string const& name) +{ + for ( ; first != last ; ++first) + { + if (first->ct.id == name) + return first; + + It parent = pj_gridinfo_parent(first->children.begin(), first->children.end(), name); + if( parent != first->children.end() ) + return parent; + } + + return last; +} + +/************************************************************************/ +/* pj_gridinfo_init_ntv2() */ +/* */ +/* Load a ntv2 (.gsb) file. */ +/************************************************************************/ + +template <typename IStream> +inline bool pj_gridinfo_init_ntv2(std::string const& gridname, + IStream & is, + pj_gridinfo & gridinfo) +{ + BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 ); + BOOST_STATIC_ASSERT( sizeof(double) == 8 ); + + static const double s2r = math::d2r<double>() / 3600.0; + + std::size_t gridinfo_orig_size = gridinfo.size(); + + // Read the overview header. + char header[11*16]; + + is.read(header, sizeof(header)); + if( is.fail() ) + { + return false; + } + + bool must_swap = (header[8] == 11) + ? !is_lsb() + : is_lsb(); + + // NOTE: This check is not implemented in proj4 + if (! cstr_equal(header + 56, "SECONDS", 7)) + { + return false; + } + + // Byte swap interesting fields if needed. + if( must_swap ) + { + swap_words( header+8, 4, 1 ); + swap_words( header+8+16, 4, 1 ); + swap_words( header+8+32, 4, 1 ); + swap_words( header+8+7*16, 8, 1 ); + swap_words( header+8+8*16, 8, 1 ); + swap_words( header+8+9*16, 8, 1 ); + swap_words( header+8+10*16, 8, 1 ); + } + + // Get the subfile count out ... all we really use for now. + boost::int32_t num_subfiles; + memcpy( &num_subfiles, header+8+32, 4 ); + + // Step through the subfiles, creating a PJ_GRIDINFO for each. + for( boost::int32_t subfile = 0; subfile < num_subfiles; subfile++ ) + { + // Read header. + is.read(header, sizeof(header)); + if( is.fail() ) + { + return false; + } + + if(! cstr_equal(header, "SUB_NAME", 8)) + { + return false; + } + + // Byte swap interesting fields if needed. + if( must_swap ) + { + swap_words( header+8+16*4, 8, 1 ); + swap_words( header+8+16*5, 8, 1 ); + swap_words( header+8+16*6, 8, 1 ); + swap_words( header+8+16*7, 8, 1 ); + swap_words( header+8+16*8, 8, 1 ); + swap_words( header+8+16*9, 8, 1 ); + swap_words( header+8+16*10, 4, 1 ); + } + + // Initialize a corresponding "ct" structure. + pj_ctable ct; + pj_ctable::lp_t ur; + + ct.id = std::string(header + 8, 8); + + ct.ll.lam = - *((double *) (header+7*16+8)); /* W_LONG */ + ct.ll.phi = *((double *) (header+4*16+8)); /* S_LAT */ + + ur.lam = - *((double *) (header+6*16+8)); /* E_LONG */ + ur.phi = *((double *) (header+5*16+8)); /* N_LAT */ + + ct.del.lam = *((double *) (header+9*16+8)); + ct.del.phi = *((double *) (header+8*16+8)); + + ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1; + ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1; + + ct.ll.lam *= s2r; + ct.ll.phi *= s2r; + ct.del.lam *= s2r; + ct.del.phi *= s2r; + + boost::int32_t gs_count; + memcpy( &gs_count, header + 8 + 16*10, 4 ); + if( gs_count != ct.lim.lam * ct.lim.phi ) + { + return false; + } + + //ct.cvs.clear(); + + // Create a new gridinfo for this if we aren't processing the + // 1st subfile, and initialize our grid info. + + // Attach to the correct list or sublist. + + // TODO is offset needed? + pj_gi gi(gridname, pj_gi::ntv2, is.tellg(), must_swap); + gi.ct = ct; + + if( subfile == 0 ) + { + gridinfo.push_back(gi); + } + else if( cstr_equal(header+24, "NONE", 4) ) + { + gridinfo.push_back(gi); + } + else + { + pj_gridinfo::iterator git = pj_gridinfo_parent(gridinfo.begin() + gridinfo_orig_size, + gridinfo.end(), + std::string((const char*)header+24, 8)); + + if( git == gridinfo.end() ) + { + gridinfo.push_back(gi); + } + else + { + git->children.push_back(gi); + } + } + + // Seek past the data. + is.seekg(gs_count * 16, std::ios::cur); + } + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_init_ntv1() */ +/* */ +/* Load an NTv1 style Canadian grid shift file. */ +/************************************************************************/ + +template <typename IStream> +inline bool pj_gridinfo_init_ntv1(std::string const& gridname, + IStream & is, + pj_gridinfo & gridinfo) +{ + BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 ); + BOOST_STATIC_ASSERT( sizeof(double) == 8 ); + + static const double d2r = math::d2r<double>(); + + // Read the header. + char header[176]; + + is.read(header, sizeof(header)); + if( is.fail() ) + { + return false; + } + + // Regularize fields of interest. + if( is_lsb() ) + { + swap_words( header+8, 4, 1 ); + swap_words( header+24, 8, 1 ); + swap_words( header+40, 8, 1 ); + swap_words( header+56, 8, 1 ); + swap_words( header+72, 8, 1 ); + swap_words( header+88, 8, 1 ); + swap_words( header+104, 8, 1 ); + } + + // NTv1 grid shift file has wrong record count, corrupt? + if( *((boost::int32_t *) (header+8)) != 12 ) + { + return false; + } + + // NOTE: This check is not implemented in proj4 + if (! cstr_equal(header + 120, "SECONDS", 7)) + { + return false; + } + + // Fill in CTABLE structure. + pj_ctable ct; + pj_ctable::lp_t ur; + + ct.id = "NTv1 Grid Shift File"; + + ct.ll.lam = - *((double *) (header+72)); + ct.ll.phi = *((double *) (header+24)); + ur.lam = - *((double *) (header+56)); + ur.phi = *((double *) (header+40)); + ct.del.lam = *((double *) (header+104)); + ct.del.phi = *((double *) (header+88)); + ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1; + ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1; + + ct.ll.lam *= d2r; + ct.ll.phi *= d2r; + ct.del.lam *= d2r; + ct.del.phi *= d2r; + //ct.cvs.clear(); + + // is offset needed? + gridinfo.push_back(pj_gi(gridname, pj_gi::ntv1, is.tellg())); + gridinfo.back().ct = ct; + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_init_gtx() */ +/* */ +/* Load a NOAA .gtx vertical datum shift file. */ +/************************************************************************/ + +template <typename IStream> +inline bool pj_gridinfo_init_gtx(std::string const& gridname, + IStream & is, + pj_gridinfo & gridinfo) +{ + BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 ); + BOOST_STATIC_ASSERT( sizeof(double) == 8 ); + + static const double d2r = math::d2r<double>(); + + // Read the header. + char header[40]; + + is.read(header, sizeof(header)); + if( is.fail() ) + { + return false; + } + + // Regularize fields of interest and extract. + double xorigin, yorigin, xstep, ystep; + boost::int32_t rows, columns; + + if( is_lsb() ) + { + swap_words( header+0, 8, 4 ); + swap_words( header+32, 4, 2 ); + } + + memcpy( &yorigin, header+0, 8 ); + memcpy( &xorigin, header+8, 8 ); + memcpy( &ystep, header+16, 8 ); + memcpy( &xstep, header+24, 8 ); + + memcpy( &rows, header+32, 4 ); + memcpy( &columns, header+36, 4 ); + + // gtx file header has invalid extents, corrupt? + if( xorigin < -360 || xorigin > 360 + || yorigin < -90 || yorigin > 90 ) + { + return false; + } + + // Fill in CTABLE structure. + pj_ctable ct; + + ct.id = "GTX Vertical Grid Shift File"; + + ct.ll.lam = xorigin; + ct.ll.phi = yorigin; + ct.del.lam = xstep; + ct.del.phi = ystep; + ct.lim.lam = columns; + ct.lim.phi = rows; + + // some GTX files come in 0-360 and we shift them back into the + // expected -180 to 180 range if possible. This does not solve + // problems with grids spanning the dateline. + if( ct.ll.lam >= 180.0 ) + ct.ll.lam -= 360.0; + + if( ct.ll.lam >= 0.0 && ct.ll.lam + ct.del.lam * ct.lim.lam > 180.0 ) + { + //"This GTX spans the dateline! This will cause problems." ); + } + + ct.ll.lam *= d2r; + ct.ll.phi *= d2r; + ct.del.lam *= d2r; + ct.del.phi *= d2r; + //ct.cvs.clear(); + + // is offset needed? + gridinfo.push_back(pj_gi(gridname, pj_gi::gtx, 40)); + gridinfo.back().ct = ct; + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_init_ctable2() */ +/* */ +/* Read the header portion of a "ctable2" format grid. */ +/************************************************************************/ + +// Originally nad_ctable2_init() defined in nad_init.c +template <typename IStream> +inline bool pj_gridinfo_init_ctable2(std::string const& gridname, + IStream & is, + pj_gridinfo & gridinfo) +{ + BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 ); + BOOST_STATIC_ASSERT( sizeof(double) == 8 ); + + char header[160]; + + is.read(header, sizeof(header)); + if( is.fail() ) + { + return false; + } + + if( !is_lsb() ) + { + swap_words( header + 96, 8, 4 ); + swap_words( header + 128, 4, 2 ); + } + + // ctable2 - wrong header! + if (! cstr_equal(header, "CTABLE V2", 9)) + { + return false; + } + + // read the table header + pj_ctable ct; + + ct.id = std::string(header + 16, std::find(header + 16, header + 16 + 80, '\0')); + //memcpy( &ct.ll.lam, header + 96, 8 ); + //memcpy( &ct.ll.phi, header + 104, 8 ); + //memcpy( &ct.del.lam, header + 112, 8 ); + //memcpy( &ct.del.phi, header + 120, 8 ); + //memcpy( &ct.lim.lam, header + 128, 4 ); + //memcpy( &ct.lim.phi, header + 132, 4 ); + memcpy( &ct.ll, header + 96, 40 ); + + // do some minimal validation to ensure the structure isn't corrupt + if ( ct.lim.lam < 1 || ct.lim.lam > 100000 + || ct.lim.phi < 1 || ct.lim.phi > 100000 ) + { + return false; + } + + // trim white space and newlines off id + boost::trim_right_if(ct.id, is_trimmable_char()); + + //ct.cvs.clear(); + + gridinfo.push_back(pj_gi(gridname, pj_gi::ctable2)); + gridinfo.back().ct = ct; + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_init_ctable() */ +/* */ +/* Read the header portion of a "ctable" format grid. */ +/************************************************************************/ + +// Originally nad_ctable_init() defined in nad_init.c +template <typename IStream> +inline bool pj_gridinfo_init_ctable(std::string const& gridname, + IStream & is, + pj_gridinfo & gridinfo) +{ + BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 ); + BOOST_STATIC_ASSERT( sizeof(double) == 8 ); + + // 80 + 2*8 + 2*8 + 2*4 + char header[120]; + + // NOTE: in proj4 data is loaded directly into CTABLE + + is.read(header, sizeof(header)); + if( is.fail() ) + { + return false; + } + + // NOTE: in proj4 LSB is not checked here + + // read the table header + pj_ctable ct; + + ct.id = std::string(header, std::find(header, header + 80, '\0')); + memcpy( &ct.ll, header + 80, 40 ); + + // do some minimal validation to ensure the structure isn't corrupt + if ( ct.lim.lam < 1 || ct.lim.lam > 100000 + || ct.lim.phi < 1 || ct.lim.phi > 100000 ) + { + return false; + } + + // trim white space and newlines off id + boost::trim_right_if(ct.id, is_trimmable_char()); + + //ct.cvs.clear(); + + gridinfo.push_back(pj_gi(gridname, pj_gi::ctable)); + gridinfo.back().ct = ct; + + return true; +} + +/************************************************************************/ +/* pj_gridinfo_init() */ +/* */ +/* Open and parse header details from a datum gridshift file */ +/* returning a list of PJ_GRIDINFOs for the grids in that */ +/* file. This superceeds use of nad_init() for modern */ +/* applications. */ +/************************************************************************/ + +template <typename IStream> +inline bool pj_gridinfo_init(std::string const& gridname, + IStream & is, + pj_gridinfo & gridinfo) +{ + char header[160]; + + // Check if the stream is opened. + if (! is.is_open()) { + return false; + } + + // Load a header, to determine the file type. + is.read(header, sizeof(header)); + + if ( is.fail() ) { + return false; + } + + is.seekg(0); + + // Determine file type. + if ( cstr_equal(header + 0, "HEADER", 6) + && cstr_equal(header + 96, "W GRID", 6) + && cstr_equal(header + 144, "TO NAD83 ", 16) ) + { + return pj_gridinfo_init_ntv1(gridname, is, gridinfo); + } + else if( cstr_equal(header + 0, "NUM_OREC", 8) + && cstr_equal(header + 48, "GS_TYPE", 7) ) + { + return pj_gridinfo_init_ntv2(gridname, is, gridinfo); + } + else if( boost::algorithm::ends_with(gridname, "gtx") + || boost::algorithm::ends_with(gridname, "GTX") ) + { + return pj_gridinfo_init_gtx(gridname, is, gridinfo); + } + else if( cstr_equal(header + 0, "CTABLE V2", 9) ) + { + return pj_gridinfo_init_ctable2(gridname, is, gridinfo); + } + else + { + return pj_gridinfo_init_ctable(gridname, is, gridinfo); + } +} + + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP diff --git a/boost/geometry/srs/projections/impl/pj_gridlist.hpp b/boost/geometry/srs/projections/impl/pj_gridlist.hpp new file mode 100644 index 0000000000..43fc3b70ab --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_gridlist.hpp @@ -0,0 +1,181 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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: +// Author: Frank Warmerdam, warmerdam@pobox.com + +// 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_GRIDLIST_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDLIST_HPP + + +#include <boost/geometry/srs/projections/grids.hpp> +#include <boost/geometry/srs/projections/impl/pj_gridinfo.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +/************************************************************************/ +/* pj_gridlist_merge_grid() */ +/* */ +/* Find/load the named gridfile and merge it into the */ +/* last_nadgrids_list. */ +/************************************************************************/ + +// Originally one function, here divided into several functions +// with overloads for various types of grids and stream policies + +inline bool pj_gridlist_find_all(std::string const& gridname, + pj_gridinfo const& grids, + std::vector<std::size_t> & gridindexes) +{ + bool result = false; + for (std::size_t i = 0 ; i < grids.size() ; ++i) + { + if (grids[i].gridname == gridname) + { + result = true; + gridindexes.push_back(i); + } + } + return result; +} + +// Fill container with sequentially increasing numbers +inline void pj_gridlist_add_seq_inc(std::vector<std::size_t> & gridindexes, + std::size_t first, std::size_t last) +{ + gridindexes.reserve(gridindexes.size() + (last - first)); + for ( ; first < last ; ++first) + { + gridindexes.push_back(first); + } +} + +// Generic stream policy and standard grids +template <typename StreamPolicy> +inline bool pj_gridlist_merge_gridfile(std::string const& gridname, + StreamPolicy const& stream_policy, + srs::grids & grids, + std::vector<std::size_t> & gridindexes) +{ + // Try to find in the existing list of loaded grids. Add all + // matching grids as with NTv2 we can get many grids from one + // file (one shared gridname). + if (pj_gridlist_find_all(gridname, grids.gridinfo, gridindexes)) + return true; + + std::size_t orig_size = grids.gridinfo.size(); + + // Try to load the named grid. + typename StreamPolicy::stream_type is; + stream_policy.open(is, gridname); + + if (! pj_gridinfo_init(gridname, is, grids.gridinfo)) + { + return false; + } + + // Add the grid now that it is loaded. + pj_gridlist_add_seq_inc(gridindexes, orig_size, grids.gridinfo.size()); + + return true; +} + + +/************************************************************************/ +/* pj_gridlist_from_nadgrids() */ +/* */ +/* This functions loads the list of grids corresponding to a */ +/* particular nadgrids string into a list, and returns it. The */ +/* list is kept around till a request is made with a different */ +/* string in order to cut down on the string parsing cost, and */ +/* the cost of building the list of tables each time. */ +/************************************************************************/ + +template <typename StreamPolicy, typename Grids> +inline void pj_gridlist_from_nadgrids(std::string const& nadgrids, + StreamPolicy const& stream_policy, + Grids & grids, + std::vector<std::size_t> & gridindexes) + +{ + // Loop processing names out of nadgrids one at a time. + for (std::string::size_type i = 0 ; i < nadgrids.size() ; ) + { + bool required = true; + + if( nadgrids[i] == '@' ) + { + required = false; + ++i; + } + + std::string::size_type end = nadgrids.find(',', i); + std::string name = nadgrids.substr(i, end - i); + + i = end; + if (end != std::string::npos) + ++i; + + if ( ! pj_gridlist_merge_gridfile(name, stream_policy, grids, gridindexes) + && required ) + { + BOOST_THROW_EXCEPTION( projection_exception(error_failed_to_load_grid) ); + } + } +} + +template <typename Par, typename GridsStorage> +inline void pj_gridlist_from_nadgrids(Par const& defn, srs::projection_grids<GridsStorage> & grids) +{ + BOOST_GEOMETRY_ASSERT(grids.storage_ptr != NULL); + + pj_gridlist_from_nadgrids(pj_get_param_s(defn.params, "nadgrids"), + grids.storage_ptr->stream_policy, + grids.storage_ptr->hgrids, + grids.hindexes); +} + + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDLIST_HPP diff --git a/boost/geometry/srs/projections/impl/pj_gridlist_shared.hpp b/boost/geometry/srs/projections/impl/pj_gridlist_shared.hpp new file mode 100644 index 0000000000..88753d034d --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_gridlist_shared.hpp @@ -0,0 +1,121 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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: +// Author: Frank Warmerdam, warmerdam@pobox.com + +// 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_GRIDLIST_SHARED_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDLIST_SHARED_HPP + + +#include <boost/geometry/srs/projections/shared_grids.hpp> +#include <boost/geometry/srs/projections/impl/pj_gridinfo.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + + +/************************************************************************/ +/* pj_gridlist_merge_grid() */ +/* */ +/* Find/load the named gridfile and merge it into the */ +/* last_nadgrids_list. */ +/************************************************************************/ + +// Generic stream policy and shared grids +template <typename StreamPolicy> +inline bool pj_gridlist_merge_gridfile(std::string const& gridname, + StreamPolicy const& stream_policy, + srs::shared_grids & grids, + std::vector<std::size_t> & gridindexes) +{ + // Try to find in the existing list of loaded grids. Add all + // matching grids as with NTv2 we can get many grids from one + // file (one shared gridname). + { + boost::shared_lock<boost::shared_mutex> lock(grids.mutex); + + if (pj_gridlist_find_all(gridname, grids.gridinfo, gridindexes)) + return true; + } + + // Try to load the named grid. + typename StreamPolicy::stream_type is; + stream_policy.open(is, gridname); + + pj_gridinfo new_grids; + + if (! pj_gridinfo_init(gridname, is, new_grids)) + { + return false; + } + + // Add the grid now that it is loaded. + + std::size_t orig_size = 0; + std::size_t new_size = 0; + + { + boost::unique_lock<boost::shared_mutex> lock(grids.mutex); + + // Try to find in the existing list of loaded grids again + // in case other thread already added it. + if (pj_gridlist_find_all(gridname, grids.gridinfo, gridindexes)) + return true; + + orig_size = grids.gridinfo.size(); + new_size = orig_size + new_grids.size(); + + grids.gridinfo.resize(new_size); + for (std::size_t i = 0 ; i < new_grids.size() ; ++ i) + new_grids[i].swap(grids.gridinfo[i + orig_size]); + } + + pj_gridlist_add_seq_inc(gridindexes, orig_size, new_size); + + return true; +} + + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDLIST_SHARED_HPP diff --git a/boost/geometry/srs/projections/impl/pj_init.hpp b/boost/geometry/srs/projections/impl/pj_init.hpp index 8b93343bcd..1d2db1af24 100644 --- a/boost/geometry/srs/projections/impl/pj_init.hpp +++ b/boost/geometry/srs/projections/impl/pj_init.hpp @@ -44,7 +44,6 @@ #include <vector> #include <boost/algorithm/string.hpp> -#include <boost/lexical_cast.hpp> #include <boost/range.hpp> #include <boost/type_traits/is_same.hpp> @@ -71,21 +70,21 @@ namespace detail template <typename BGParams, typename T> inline void pj_push_defaults(BGParams const& /*bg_params*/, parameters<T>& pin) { - pin.params.push_back(pj_mkparam<T>("ellps=WGS84")); + pin.params.push_back(pj_mkparam<T>("ellps", "WGS84")); if (pin.name == "aea") { - pin.params.push_back(pj_mkparam<T>("lat_1=29.5")); - pin.params.push_back(pj_mkparam<T>("lat_2=45.5 ")); + pin.params.push_back(pj_mkparam<T>("lat_1", "29.5")); + pin.params.push_back(pj_mkparam<T>("lat_2", "45.5 ")); } else if (pin.name == "lcc") { - pin.params.push_back(pj_mkparam<T>("lat_1=33")); - pin.params.push_back(pj_mkparam<T>("lat_2=45")); + pin.params.push_back(pj_mkparam<T>("lat_1", "33")); + pin.params.push_back(pj_mkparam<T>("lat_2", "45")); } else if (pin.name == "lagrng") { - pin.params.push_back(pj_mkparam<T>("W=2")); + pin.params.push_back(pj_mkparam<T>("W", "2")); } } @@ -100,21 +99,21 @@ inline void pj_push_defaults(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL >::type proj_tag; // statically defaulting to WGS84 - //pin.params.push_back(pj_mkparam("ellps=WGS84")); + //pin.params.push_back(pj_mkparam("ellps", "WGS84")); if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::aea>::value))) { - pin.params.push_back(pj_mkparam<T>("lat_1=29.5")); - pin.params.push_back(pj_mkparam<T>("lat_2=45.5 ")); + pin.params.push_back(pj_mkparam<T>("lat_1", "29.5")); + pin.params.push_back(pj_mkparam<T>("lat_2", "45.5 ")); } else if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::lcc>::value))) { - pin.params.push_back(pj_mkparam<T>("lat_1=33")); - pin.params.push_back(pj_mkparam<T>("lat_2=45")); + pin.params.push_back(pj_mkparam<T>("lat_1", "33")); + pin.params.push_back(pj_mkparam<T>("lat_2", "45")); } else if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::lagrng>::value))) { - pin.params.push_back(pj_mkparam<T>("W=2")); + pin.params.push_back(pj_mkparam<T>("W", "2")); } } @@ -128,7 +127,7 @@ inline void pj_init_units(std::vector<pvalue<T> > const& params, T const& default_fr_meter) { std::string s; - std::string units = pj_param(params, sunits).s; + std::string units = pj_get_param_s(params, sunits); if (! units.empty()) { const int n = sizeof(pj_units) / sizeof(pj_units[0]); @@ -142,14 +141,14 @@ inline void pj_init_units(std::vector<pvalue<T> > const& params, } if (index == -1) { - BOOST_THROW_EXCEPTION( projection_exception(-7) ); + BOOST_THROW_EXCEPTION( projection_exception(error_unknow_unit_id) ); } s = pj_units[index].to_meter; } if (s.empty()) { - s = pj_param(params, sto_meter).s; + s = pj_get_param_s(params, sto_meter); } if (! s.empty()) @@ -157,21 +156,21 @@ inline void pj_init_units(std::vector<pvalue<T> > const& params, std::size_t const pos = s.find('/'); if (pos == std::string::npos) { - to_meter = lexical_cast<T>(s); + to_meter = geometry::str_cast<T>(s); } else { - T const numerator = lexical_cast<T>(s.substr(0, pos)); - T const denominator = lexical_cast<T>(s.substr(pos + 1)); + T const numerator = geometry::str_cast<T>(s.substr(0, pos)); + T const denominator = geometry::str_cast<T>(s.substr(pos + 1)); if (numerator == 0.0 || denominator == 0.0) { - BOOST_THROW_EXCEPTION( projection_exception(-99) ); + BOOST_THROW_EXCEPTION( projection_exception(error_unit_factor_less_than_0) ); } to_meter = numerator / denominator; } if (to_meter == 0.0) { - BOOST_THROW_EXCEPTION( projection_exception(-99) ); + BOOST_THROW_EXCEPTION( projection_exception(error_unit_factor_less_than_0) ); } fr_meter = 1. / to_meter; } @@ -200,15 +199,16 @@ inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool pin.params.push_back(pj_mkparam<T>(*it)); } + // maybe TODO: handle "init" parameter /* check if +init present */ - if (pj_param(pin.params, "tinit").i) - { - // maybe TODO: handle "init" parameter - //if (!(curr = get_init(&arguments, curr, pj_param(pin.params, "sinit").s))) - } + //std::string sinit; + //if (pj_param_s(pin.params, "init", sinit)) + //{ + // //if (!(curr = get_init(&arguments, curr, sinit))) + //} // find projection -> implemented in projection factory - pin.name = pj_param(pin.params, "sproj").s; + pin.name = pj_get_param_s(pin.params, "proj"); // exception thrown in projection<> // TODO: consider throwing here both projection_unknown_id_exception and // projection_not_named_exception in order to throw before other exceptions @@ -217,7 +217,7 @@ inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool // set defaults, unless inhibited // GL-Addition, if use_defaults is false then defaults are ignored - if (use_defaults && ! pj_param(pin.params, "bno_defs").i) + if (use_defaults && ! pj_get_param_b(pin.params, "no_defs")) { // proj4 gets defaults from "proj_def.dat", file of 94/02/23 with a few defaults. // Here manually @@ -246,61 +246,59 @@ inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool pin.ra = 1. / pin.a; pin.one_es = 1. - pin.es; if (pin.one_es == 0.) { - BOOST_THROW_EXCEPTION( projection_exception(-6) ); + BOOST_THROW_EXCEPTION( projection_exception(error_eccentricity_is_one) ); } pin.rone_es = 1./pin.one_es; /* Now that we have ellipse information check for WGS84 datum */ - if( pin.datum_type == PJD_3PARAM + if( pin.datum_type == datum_3param && pin.datum_params[0] == 0.0 && pin.datum_params[1] == 0.0 && pin.datum_params[2] == 0.0 && pin.a == 6378137.0 && geometry::math::abs(pin.es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/ { - pin.datum_type = PJD_WGS84; + pin.datum_type = datum_wgs84; } /* set pin.geoc coordinate system */ - pin.geoc = (pin.es && pj_param(pin.params, "bgeoc").i); + pin.geoc = (pin.es && pj_get_param_b(pin.params, "geoc")); /* over-ranging flag */ - pin.over = pj_param(pin.params, "bover").i; + pin.over = pj_get_param_b(pin.params, "over"); /* longitude center for wrapping */ - pin.is_long_wrap_set = pj_param(pin.params, "tlon_wrap").i != 0; - if (pin.is_long_wrap_set) - pin.long_wrap_center = pj_param(pin.params, "rlon_wrap").f; + pin.is_long_wrap_set = pj_param_r(pin.params, "lon_wrap", pin.long_wrap_center); /* central meridian */ - pin.lam0 = pj_param(pin.params, "rlon_0").f; + pin.lam0 = pj_get_param_r(pin.params, "lon_0"); /* central latitude */ - pin.phi0 = pj_param(pin.params, "rlat_0").f; + pin.phi0 = pj_get_param_r(pin.params, "lat_0"); /* false easting and northing */ - pin.x0 = pj_param(pin.params, "dx_0").f; - pin.y0 = pj_param(pin.params, "dy_0").f; + pin.x0 = pj_get_param_f(pin.params, "x_0"); + pin.y0 = pj_get_param_f(pin.params, "y_0"); /* general scaling factor */ - if (pj_param(pin.params, "tk_0").i) - pin.k0 = pj_param(pin.params, "dk_0").f; - else if (pj_param(pin.params, "tk").i) - pin.k0 = pj_param(pin.params, "dk").f; - else + if (pj_param_f(pin.params, "k_0", pin.k0)) { + /* empty */ + } else if (pj_param_f(pin.params, "k", pin.k0)) { + /* empty */ + } else pin.k0 = 1.; if (pin.k0 <= 0.) { - BOOST_THROW_EXCEPTION( projection_exception(-31) ); + BOOST_THROW_EXCEPTION( projection_exception(error_k_less_than_zero) ); } /* set units */ - pj_init_units(pin.params, "sunits", "sto_meter", + pj_init_units(pin.params, "units", "to_meter", pin.to_meter, pin.fr_meter, 1., 1.); - pj_init_units(pin.params, "svunits", "svto_meter", + pj_init_units(pin.params, "vunits", "vto_meter", pin.vto_meter, pin.vfr_meter, pin.to_meter, pin.fr_meter); /* prime meridian */ - std::string pm = pj_param(pin.params, "spm").s; + std::string pm = pj_get_param_s(pin.params, "pm"); if (! pm.empty()) { std::string value; @@ -317,8 +315,8 @@ inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool dms_parser<T, true> parser; - // TODO: Handle case when lexical_cast is not used consistently. - // This should probably be done in dms_parser. + // TODO: Is this try-catch needed? + // In other cases the bad_str_cast exception is simply thrown BOOST_TRY { if (value.empty()) { @@ -327,9 +325,9 @@ inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool pin.from_greenwich = parser.apply(value).angle(); } } - BOOST_CATCH(boost::bad_lexical_cast const&) + BOOST_CATCH(geometry::bad_str_cast const&) { - BOOST_THROW_EXCEPTION( projection_exception(-46) ); + BOOST_THROW_EXCEPTION( projection_exception(error_unknown_prime_meridian) ); } BOOST_CATCH_END } diff --git a/boost/geometry/srs/projections/impl/pj_mlfn.hpp b/boost/geometry/srs/projections/impl/pj_mlfn.hpp index 1317dd2e2f..04f0d19442 100644 --- a/boost/geometry/srs/projections/impl/pj_mlfn.hpp +++ b/boost/geometry/srs/projections/impl/pj_mlfn.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -36,10 +36,17 @@ // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER // DEALINGS IN THE SOFTWARE. +/* meridional distance for ellipsoid and inverse +** 8th degree - accurate to < 1e-5 meters when used in conjunction +** with typical major axis values. +** Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds. +*/ + #ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_MLFN_HPP #define BOOST_GEOMETRY_PROJECTIONS_PJ_MLFN_HPP +#include <cstdlib> #include <boost/geometry/util/math.hpp> @@ -48,10 +55,20 @@ namespace boost { namespace geometry { namespace projections { namespace detail { -static const int EN_SIZE = 5; +template <typename T> +struct en +{ + static const std::size_t size = 5; + + T const& operator[](size_t i) const { return data[i]; } + T & operator[](size_t i) { return data[i]; } + +private: + T data[5]; +}; template <typename T> -inline bool pj_enfn(T const& es, T* en) +inline en<T> pj_enfn(T const& es) { static const T C00 = 1.; static const T C02 = .25; @@ -66,9 +83,9 @@ inline bool pj_enfn(T const& es, T* en) static const T C68 = .00569661458333333333; static const T C88 = .3076171875; - T t; //, *en; + T t; + detail::en<T> en; - //if (en = (double *)pj_malloc(EN_SIZE * sizeof(double))) { en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08))); en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08))); @@ -76,12 +93,12 @@ inline bool pj_enfn(T const& es, T* en) en[3] = (t *= es) * (C66 - es * C68); en[4] = t * es * C88; } - // return en; - return true; + + return en; } template <typename T> -inline T pj_mlfn(T const& phi, T sphi, T cphi, const T *en) +inline T pj_mlfn(T const& phi, T sphi, T cphi, detail::en<T> const& en) { cphi *= sphi; sphi *= sphi; @@ -90,13 +107,8 @@ inline T pj_mlfn(T const& phi, T sphi, T cphi, const T *en) } template <typename T> -inline T pj_inv_mlfn(T const& arg, T const& es, const T *en) +inline T pj_inv_mlfn(T const& arg, T const& es, detail::en<T> const& en) { - /* meridinal distance for ellipsoid and inverse - ** 8th degree - accurate to < 1e-5 meters when used in conjuction - ** with typical major axis values. - ** Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds. - */ static const T EPS = 1e-11; static const int MAX_ITER = 10; @@ -111,7 +123,7 @@ inline T pj_inv_mlfn(T const& arg, T const& es, const T *en) if (geometry::math::abs(t) < EPS) return phi; } - BOOST_THROW_EXCEPTION( projection_exception(-17) ); + BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) ); return phi; } diff --git a/boost/geometry/srs/projections/impl/pj_param.hpp b/boost/geometry/srs/projections/impl/pj_param.hpp index 4f33ad837f..7648055414 100644 --- a/boost/geometry/srs/projections/impl/pj_param.hpp +++ b/boost/geometry/srs/projections/impl/pj_param.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -43,15 +43,30 @@ #include <string> #include <vector> +#include <boost/geometry/srs/projections/exception.hpp> + #include <boost/geometry/srs/projections/impl/dms_parser.hpp> #include <boost/geometry/srs/projections/impl/projects.hpp> +#include <boost/mpl/assert.hpp> +#include <boost/type_traits/is_integral.hpp> + namespace boost { namespace geometry { namespace projections { namespace detail { +/* create pvalue list entry */ +template <typename T> +inline pvalue<T> pj_mkparam(std::string const& name, std::string const& value) +{ + pvalue<T> newitem; + newitem.param = name; + newitem.s = value; + //newitem.used = false; + return newitem; +} /* create pvalue list entry */ template <typename T> @@ -67,92 +82,146 @@ inline pvalue<T> pj_mkparam(std::string const& str) name.erase(loc); } + return pj_mkparam<T>(name, value); +} - pvalue<T> newitem; - newitem.param = name; - newitem.s = value; - newitem.used = 0; - newitem.i = atoi(value.c_str()); - newitem.f = atof(value.c_str()); - return newitem; +/* input exists */ +template <typename T> +inline typename std::vector<pvalue<T> >::const_iterator + pj_param_find(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + typedef typename std::vector<pvalue<T> >::const_iterator iterator; + for (iterator it = pl.begin(); it != pl.end(); it++) + { + if (it->param == name) + { + //it->used = true; + return it; + } + // TODO: needed for pipeline + /*else if (it->param == "step") + { + return pl.end(); + }*/ + } + + return pl.end(); } -/************************************************************************/ -/* pj_param() */ -/* */ -/* Test for presence or get pvalue value. The first */ -/* character in `opt' is a pvalue type which can take the */ -/* values: */ -/* */ -/* `t' - test for presence, return TRUE/FALSE in pvalue.i */ -/* `i' - integer value returned in pvalue.i */ -/* `d' - simple valued real input returned in pvalue.f */ -/* `r' - degrees (DMS translation applied), returned as */ -/* radians in pvalue.f */ -/* `s' - string returned in pvalue.s */ -/* `b' - test for t/T/f/F, return in pvalue.i */ -/* */ -/************************************************************************/ +/* input exists */ +template <typename T> +inline bool pj_param_exists(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + return pj_param_find(pl, name) != pl.end(); +} +/* integer input */ template <typename T> -inline pvalue<T> pj_param(std::vector<pvalue<T> > const& pl, std::string opt) +inline bool pj_param_i(std::vector<pvalue<T> > const& pl, std::string const& name, int & par) { - char type = opt[0]; - opt.erase(opt.begin()); + typename std::vector<pvalue<T> >::const_iterator it = pj_param_find(pl, name); + if (it != pl.end()) + { + par = geometry::str_cast<int>(it->s); + return true; + } + return false; +} - pvalue<T> value; +/* floating point input */ +template <typename T> +inline bool pj_param_f(std::vector<pvalue<T> > const& pl, std::string const& name, T & par) +{ + typename std::vector<pvalue<T> >::const_iterator it = pj_param_find(pl, name); + if (it != pl.end()) + { + par = geometry::str_cast<T>(it->s); + return true; + } + return false; +} - /* simple linear lookup */ - typedef typename std::vector<pvalue<T> >::const_iterator iterator; - for (iterator it = pl.begin(); it != pl.end(); it++) +/* radians input */ +template <typename T> +inline bool pj_param_r(std::vector<pvalue<T> > const& pl, std::string const& name, T & par) +{ + typename std::vector<pvalue<T> >::const_iterator it = pj_param_find(pl, name); + if (it != pl.end()) + { + dms_parser<T, true> parser; + par = parser.apply(it->s.c_str()).angle(); + return true; + } + return false; +} + +/* string input */ +template <typename T> +inline bool pj_param_s(std::vector<pvalue<T> > const& pl, std::string const& name, std::string & par) +{ + typename std::vector<pvalue<T> >::const_iterator it = pj_param_find(pl, name); + if (it != pl.end()) + { + par = it->s; + return true; + } + return false; +} + +/* bool input */ +template <typename T> +inline bool pj_get_param_b(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + typename std::vector<pvalue<T> >::const_iterator it = pj_param_find(pl, name); + if (it != pl.end()) { - if (it->param == opt) + switch (it->s[0]) { - //it->used = 1; - switch (type) - { - case 't': - value.i = 1; - break; - case 'i': /* integer input */ - value.i = atoi(it->s.c_str()); - break; - case 'd': /* simple real input */ - value.f = atof(it->s.c_str()); - break; - case 'r': /* degrees input */ - { - dms_parser<T, true> parser; - value.f = parser.apply(it->s.c_str()).angle(); - } - break; - case 's': /* char string */ - value.s = it->s; - break; - case 'b': /* boolean */ - switch (it->s[0]) - { - case 'F': case 'f': - value.i = 0; - break; - case '\0': case 'T': case 't': - value.i = 1; - break; - default: - value.i = 0; - break; - } - break; - } - return value; + case '\0': case 'T': case 't': + return true; + case 'F': case 'f': + return false; + default: + BOOST_THROW_EXCEPTION( projection_exception(error_invalid_boolean_param) ); + return false; } + } + return false; +} - } +// NOTE: In the original code, in pl_ell_set.c there is a function pj_get_param +// which behavior is similar to pj_param but it doesn't set `user` member to TRUE +// while pj_param does in the original code. In Boost.Geometry this member is not used. +template <typename T> +inline int pj_get_param_i(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + int res = 0; + pj_param_i(pl, name, res); + return res; +} - value.i = 0; - value.f = 0.0; - value.s = ""; - return value; +template <typename T> +inline T pj_get_param_f(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + T res = 0; + pj_param_f(pl, name, res); + return res; +} + +template <typename T> +inline T pj_get_param_r(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + T res = 0; + pj_param_r(pl, name, res); + return res; +} + +template <typename T> +inline std::string pj_get_param_s(std::vector<pvalue<T> > const& pl, std::string const& name) +{ + std::string res; + pj_param_s(pl, name, res); + return res; } } // namespace detail diff --git a/boost/geometry/srs/projections/impl/pj_phi2.hpp b/boost/geometry/srs/projections/impl/pj_phi2.hpp index 71f0cf1249..868a8c659b 100644 --- a/boost/geometry/srs/projections/impl/pj_phi2.hpp +++ b/boost/geometry/srs/projections/impl/pj_phi2.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -58,12 +58,12 @@ inline T pj_phi2(T const& ts, T const& e) i = N_ITER; do { con = e * sin (Phi); - dphi = geometry::math::half_pi<T>() - 2. * atan (ts * pow((1. - con) / - (1. + con), eccnth)) - Phi; + dphi = geometry::math::half_pi<T>() - 2. * atan (ts * math::pow((T(1) - con) / + (T(1) + con), eccnth)) - Phi; Phi += dphi; } while ( geometry::math::abs(dphi) > TOL && --i); if (i <= 0) - BOOST_THROW_EXCEPTION( projection_exception(-18) ); + BOOST_THROW_EXCEPTION( projection_exception(error_non_con_inv_phi2) ); return Phi; } diff --git a/boost/geometry/srs/projections/impl/pj_strerrno.hpp b/boost/geometry/srs/projections/impl/pj_strerrno.hpp index 22e0c48af4..418bea9266 100644 --- a/boost/geometry/srs/projections/impl/pj_strerrno.hpp +++ b/boost/geometry/srs/projections/impl/pj_strerrno.hpp @@ -86,6 +86,17 @@ pj_err_list[] = { "point not within available datum shift grids", /* -48 */ "invalid sweep axis, choose x or y", /* -49 */ "malformed pipeline", /* -50 */ + "unit conversion factor must be > 0", /* -51 */ + "invalid scale", /* -52 */ + "non-convergent computation", /* -53 */ + "missing required arguments", /* -54 */ + "lat_0 = 0", /* -55 */ + "ellipsoidal usage unsupported", /* -56 */ + "only one +init allowed for non-pipeline operations", /* -57 */ + "argument not numerical or out of range", /* -58 */ + + /* When adding error messages, remember to update ID defines in + projects.h, and transient_error array in pj_transform */ }; inline std::string pj_generic_strerrno(std::string const& msg, int err) diff --git a/boost/geometry/srs/projections/impl/pj_transform.hpp b/boost/geometry/srs/projections/impl/pj_transform.hpp index 8c2095642f..bccf150db0 100644 --- a/boost/geometry/srs/projections/impl/pj_transform.hpp +++ b/boost/geometry/srs/projections/impl/pj_transform.hpp @@ -44,6 +44,7 @@ #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> @@ -157,14 +158,6 @@ private: // 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> @@ -182,22 +175,32 @@ 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. +** 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[50] = { +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 }; + /* 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> @@ -216,10 +219,18 @@ inline int pj_geodetic_to_geocentric( T const& a, T const& es, /* projection specific components). */ /************************************************************************/ -template <typename SrcPrj, typename DstPrj2, typename Par, typename Range> +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) + Range & range, + Grids const& srcgrids, + Grids const& dstgrids) { typedef typename boost::range_value<Range>::type point_type; @@ -253,7 +264,7 @@ inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn, { // Point should be cartesian 3D (ECEF) if (dimension < 3) - BOOST_THROW_EXCEPTION( projection_exception(PJD_ERR_GEOCENTRIC) ); + BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) ); //return PJD_ERR_GEOCENTRIC; if( srcdefn.to_meter != 1.0 ) @@ -406,7 +417,7 @@ inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn, /* -------------------------------------------------------------------- */ /* Convert datums if needed, and possible. */ /* -------------------------------------------------------------------- */ - if ( ! pj_datum_transform( srcdefn, dstdefn, range ) ) + if ( ! pj_datum_transform( srcdefn, dstdefn, range, srcgrids, dstgrids ) ) { result = false; } @@ -446,7 +457,7 @@ inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn, { // Point should be cartesian 3D (ECEF) if (dimension < 3) - BOOST_THROW_EXCEPTION( projection_exception(PJD_ERR_GEOCENTRIC) ); + BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) ); //return PJD_ERR_GEOCENTRIC; // NOTE: In the original code the return value of the following @@ -615,7 +626,7 @@ inline int pj_geodetic_to_geocentric( T const& a, T const& es, GeocentricInfo<T> gi; if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 ) { - return PJD_ERR_GEOCENTRIC; + return error_geocentric; } for( std::size_t i = 0 ; i < point_count ; ++i ) @@ -632,7 +643,7 @@ inline int pj_geodetic_to_geocentric( T const& a, T const& es, range_wrapper.get_z(i), // Height X, Y, Z ) != 0 ) { - ret_errno = -14; + ret_errno = error_lat_or_lon_exceed_limit; set_invalid_point(point); /* but keep processing points! */ } @@ -668,7 +679,7 @@ inline int pj_geocentric_to_geodetic( T const& a, T const& es, GeocentricInfo<T> gi; if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 ) { - return PJD_ERR_GEOCENTRIC; + return error_geocentric; } for( std::size_t i = 0 ; i < point_count ; ++i ) @@ -714,13 +725,13 @@ inline bool pj_compare_datums( Par & srcdefn, Par & dstdefn ) considered identical */ return false; } - else if( srcdefn.datum_type == PJD_3PARAM ) + 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 == PJD_7PARAM ) + else if( srcdefn.datum_type == datum_7param ) { return (srcdefn.datum_params[0] == dstdefn.datum_params[0] && srcdefn.datum_params[1] == dstdefn.datum_params[1] @@ -730,10 +741,10 @@ inline bool pj_compare_datums( Par & srcdefn, Par & dstdefn ) && srcdefn.datum_params[5] == dstdefn.datum_params[5] && srcdefn.datum_params[6] == dstdefn.datum_params[6]); } - else if( srcdefn.datum_type == PJD_GRIDSHIFT ) + else if( srcdefn.datum_type == datum_gridshift ) { - return pj_param(srcdefn.params,"snadgrids").s - == pj_param(dstdefn.params,"snadgrids").s; + return pj_get_param_s(srcdefn.params,"nadgrids") + == pj_get_param_s(dstdefn.params,"nadgrids"); } else return true; @@ -754,7 +765,7 @@ inline int pj_geocentric_to_wgs84( Par const& defn, Range & rng = range_wrapper.get_range(); std::size_t point_count = boost::size(rng); - if( defn.datum_type == PJD_3PARAM ) + if( defn.datum_type == datum_3param ) { for(std::size_t i = 0; i < point_count; i++ ) { @@ -768,7 +779,7 @@ inline int pj_geocentric_to_wgs84( Par const& defn, range_wrapper.set_z(i, range_wrapper.get_z(i) + Dz_BF(defn)); } } - else if( defn.datum_type == PJD_7PARAM ) + else if( defn.datum_type == datum_7param ) { for(std::size_t i = 0; i < point_count; i++ ) { @@ -811,7 +822,7 @@ inline int pj_geocentric_from_wgs84( Par const& defn, Range & rng = range_wrapper.get_range(); std::size_t point_count = boost::size(rng); - if( defn.datum_type == PJD_3PARAM ) + if( defn.datum_type == datum_3param ) { for(std::size_t i = 0; i < point_count; i++ ) { @@ -825,7 +836,7 @@ inline int pj_geocentric_from_wgs84( Par const& defn, range_wrapper.set_z(i, range_wrapper.get_z(i) - Dz_BF(defn)); } } - else if( defn.datum_type == PJD_7PARAM ) + else if( defn.datum_type == datum_7param ) { for(std::size_t i = 0; i < point_count; i++ ) { @@ -869,11 +880,17 @@ inline bool pj_datum_check_error(int err) /* 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 ) +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) { + static const double wgs84_a = 6378137.0; + static const double wgs84_es = 0.0066943799901413165; + typedef typename Par::type calc_t; bool result = true; @@ -885,8 +902,8 @@ inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn, /* (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 ) + if( srcdefn.datum_type == datum_unknown + || dstdefn.datum_type == datum_unknown ) return result; /* -------------------------------------------------------------------- */ @@ -911,34 +928,34 @@ inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn, /* If this datum requires grid shifts, then apply it to geodetic */ /* coordinates. */ /* -------------------------------------------------------------------- */ - /*if( srcdefn.datum_type == PJD_GRIDSHIFT ) + if( srcdefn.datum_type == datum_gridshift ) { try { - pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z ); + pj_apply_gridshift_2<false>( srcdefn, range, srcgrids ); } catch (projection_exception const& e) { if (pj_datum_check_error(e.code())) { BOOST_RETHROW } } - src_a = SRS_WGS84_SEMIMAJOR; - src_es = SRS_WGS84_ESQUARED; + src_a = wgs84_a; + src_es = wgs84_es; } - if( dstdefn.datum_type == PJD_GRIDSHIFT ) + if( dstdefn.datum_type == datum_gridshift ) { - dst_a = SRS_WGS84_SEMIMAJOR; - dst_es = SRS_WGS84_ESQUARED; - }*/ + 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 == PJD_3PARAM - || srcdefn.datum_type == PJD_7PARAM - || dstdefn.datum_type == PJD_3PARAM - || dstdefn.datum_type == PJD_7PARAM) + || srcdefn.datum_type == datum_3param + || srcdefn.datum_type == datum_7param + || dstdefn.datum_type == datum_3param + || dstdefn.datum_type == datum_7param) { /* -------------------------------------------------------------------- */ /* Convert to geocentric coordinates. */ @@ -952,8 +969,8 @@ inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn, /* -------------------------------------------------------------------- */ /* Convert between datums. */ /* -------------------------------------------------------------------- */ - if( srcdefn.datum_type == PJD_3PARAM - || srcdefn.datum_type == PJD_7PARAM ) + if( srcdefn.datum_type == datum_3param + || srcdefn.datum_type == datum_7param ) { try { pj_geocentric_to_wgs84( srcdefn, z_range ); @@ -964,8 +981,8 @@ inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn, } } - if( dstdefn.datum_type == PJD_3PARAM - || dstdefn.datum_type == PJD_7PARAM ) + if( dstdefn.datum_type == datum_3param + || dstdefn.datum_type == datum_7param ) { try { pj_geocentric_from_wgs84( dstdefn, z_range ); @@ -989,15 +1006,15 @@ inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn, /* -------------------------------------------------------------------- */ /* Apply grid shift to destination if required. */ /* -------------------------------------------------------------------- */ - /*if( dstdefn.datum_type == PJD_GRIDSHIFT ) + if( dstdefn.datum_type == datum_gridshift ) { try { - pj_apply_gridshift_2( dstdefn, 1, point_count, point_offset, x, y, z ); + pj_apply_gridshift_2<true>( dstdefn, range, dstgrids ); } catch (projection_exception const& e) { if (pj_datum_check_error(e.code())) BOOST_RETHROW } - }*/ + } return result; } diff --git a/boost/geometry/srs/projections/impl/pj_tsfn.hpp b/boost/geometry/srs/projections/impl/pj_tsfn.hpp index d39a0fb2fc..46d3edcce0 100644 --- a/boost/geometry/srs/projections/impl/pj_tsfn.hpp +++ b/boost/geometry/srs/projections/impl/pj_tsfn.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -50,7 +50,7 @@ namespace detail { { sinphi *= e; return (tan (.5 * (geometry::math::half_pi<T>() - phi)) / - pow((1. - sinphi) / (1. + sinphi), .5 * e)); + math::pow((T(1) - sinphi) / (T(1) + sinphi), T(0.5) * e)); } } // namespace detail diff --git a/boost/geometry/srs/projections/impl/pj_units.hpp b/boost/geometry/srs/projections/impl/pj_units.hpp index 269b8ff92e..98e1fa4c9f 100644 --- a/boost/geometry/srs/projections/impl/pj_units.hpp +++ b/boost/geometry/srs/projections/impl/pj_units.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -44,11 +44,19 @@ namespace boost { namespace geometry { namespace projections { namespace detail { +// Originally defined in projects.h +struct pj_units_type +{ + std::string id; /* units keyword */ + std::string to_meter; /* multiply by value to get meters */ + std::string name; /* comments */ +}; + /* Field 2 that contains the multiplier to convert named units to meters ** may be expressed by either a simple floating point constant or a ** numerator/denomenator values (e.g. 1/1000) */ -static const PJ_UNITS pj_units[] = +static const pj_units_type pj_units[] = { { "km", "1000.", "Kilometer" }, { "m", "1.", "Meter" }, diff --git a/boost/geometry/srs/projections/impl/pj_zpoly1.hpp b/boost/geometry/srs/projections/impl/pj_zpoly1.hpp index 613b6d2b35..4cb1cb1f6c 100644 --- a/boost/geometry/srs/projections/impl/pj_zpoly1.hpp +++ b/boost/geometry/srs/projections/impl/pj_zpoly1.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -52,10 +52,10 @@ namespace boost { namespace geometry { namespace projections { namespace detail ** n should always be >= 1 though no checks are made */ template <typename T> - inline COMPLEX<T> - pj_zpoly1(COMPLEX<T> z, const COMPLEX<T> *C, int n) + inline pj_complex<T> + pj_zpoly1(pj_complex<T> z, const pj_complex<T> *C, int n) { - COMPLEX<T> a; + pj_complex<T> a; T t; a = *(C += n); @@ -71,14 +71,14 @@ namespace boost { namespace geometry { namespace projections { namespace detail /* evaluate complex polynomial and derivative */ template <typename T> - inline COMPLEX<T> - pj_zpolyd1(COMPLEX<T> z, const COMPLEX<T> *C, int n, COMPLEX<T> *der) + inline pj_complex<T> + pj_zpolyd1(pj_complex<T> z, const pj_complex<T> *C, int n, pj_complex<T> *der) { T t; bool first = true; - COMPLEX<T> a = *(C += n); - COMPLEX<T> b = a; + pj_complex<T> a = *(C += n); + pj_complex<T> b = a; while (n-- > 0) { if (first) diff --git a/boost/geometry/srs/projections/impl/proj_mdist.hpp b/boost/geometry/srs/projections/impl/proj_mdist.hpp index 1c325df1c6..5f67546f54 100644 --- a/boost/geometry/srs/projections/impl/proj_mdist.hpp +++ b/boost/geometry/srs/projections/impl/proj_mdist.hpp @@ -3,6 +3,10 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2018. +// Modifications copyright (c) 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) @@ -43,23 +47,23 @@ namespace boost { namespace geometry { namespace projections { namespace detail { - static const int MDIST_MAX_ITER = 20; - template <typename T> - struct MDIST + struct mdist { + static const int static_size = 20; + int nb; T es; T E; - T b[MDIST_MAX_ITER]; + T b[static_size]; }; - template <typename CT> - inline bool proj_mdist_ini(CT const& es, MDIST<CT>& b) + template <typename T> + inline bool proj_mdist_ini(T const& es, mdist<T>& b) { - CT numf, numfi, twon1, denf, denfi, ens, T, twon; - CT den, El, Es; - CT E[MDIST_MAX_ITER]; + T numf, numfi, twon1, denf, denfi, ens, t, twon; + T den, El, Es; + T E[mdist<T>::static_size]; int i, j; /* generate E(e^2) and its terms E[] */ @@ -68,12 +72,13 @@ namespace detail denf = 1.; twon = 4.; Es = El = E[0] = 1.; - for (i = 1; i < MDIST_MAX_ITER ; ++i) + for (i = 1; i < mdist<T>::static_size ; ++i) { numf *= (twon1 * twon1); den = twon * denf * denf * twon1; - T = numf/den; - Es -= (E[i] = T * ens); + t = numf/den; + E[i] = t * ens; + Es -= E[i]; ens *= es; twon *= 4.; denf *= ++denfi; @@ -103,7 +108,7 @@ namespace detail } template <typename T> - inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, MDIST<T> const& b) + inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, mdist<T> const& b) { T sc, sum, sphi2, D; int i; @@ -117,14 +122,14 @@ namespace detail } template <typename T> - inline T proj_inv_mdist(T const& dist, MDIST<T> const& b) + inline T proj_inv_mdist(T const& dist, mdist<T> const& b) { static const T TOL = 1e-14; T s, t, phi, k; int i; k = 1./(1.- b.es); - i = MDIST_MAX_ITER; + i = mdist<T>::static_size; phi = dist; while ( i-- ) { s = sin(phi); @@ -135,7 +140,7 @@ namespace detail return phi; } /* convergence failed */ - BOOST_THROW_EXCEPTION( projection_exception(-17) ); + BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) ); } } // namespace detail diff --git a/boost/geometry/srs/projections/impl/projects.hpp b/boost/geometry/srs/projections/impl/projects.hpp index 232ae67ae9..2f94ffd021 100644 --- a/boost/geometry/srs/projections/impl/projects.hpp +++ b/boost/geometry/srs/projections/impl/projects.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// 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, @@ -44,8 +44,8 @@ #include <string> #include <vector> -#include <boost/geometry/srs/projections/exception.hpp> -#include <boost/math/constants/constants.hpp> +#include <boost/config.hpp> +#include <boost/geometry/srs/projections/constants.hpp> #include <boost/mpl/if.hpp> #include <boost/type_traits/is_pod.hpp> @@ -57,199 +57,172 @@ namespace boost { namespace geometry { namespace projections namespace detail { -/* some useful constants */ -template <typename T> -inline T ONEPI() { return boost::math::constants::pi<T>(); } -template <typename T> -inline T HALFPI() { return boost::math::constants::half_pi<T>(); } -template <typename T> -inline T FORTPI() { return boost::math::constants::pi<T>() / T(4); } -template <typename T> -inline T TWOPI() { return boost::math::constants::two_pi<T>(); } -template <typename T> -inline T TWO_D_PI() { return boost::math::constants::two_div_pi<T>(); } -template <typename T> -inline T HALFPI_SQR() { return 2.4674011002723396547086227499689; } -template <typename T> -inline T PI_SQR() { return boost::math::constants::pi_sqr<T>(); } -template <typename T> -inline T THIRD() { return 0.3333333333333333333333333333333; } -template <typename T> -inline T TWOTHIRD() { return 0.6666666666666666666666666666666; } -template <typename T> -inline T PI_HALFPI() { return 4.7123889803846898576939650749193; } -template <typename T> -inline T TWOPI_HALFPI() { return 7.8539816339744830961566084581988; } -template <typename T> -inline T PI_DIV_3() { return 1.0471975511965977461542144610932; } - /* datum_type values */ -static const int PJD_UNKNOWN = 0; -static const int PJD_3PARAM = 1; -static const int PJD_7PARAM = 2; -static const int PJD_GRIDSHIFT = 3; -static const int PJD_WGS84 = 4; /* WGS84 (or anything considered equivelent) */ +enum datum_type +{ + datum_unknown = 0, + datum_3param = 1, + datum_7param = 2, + datum_gridshift = 3, + datum_wgs84 = 4 /* WGS84 (or anything considered equivelent) */ +}; /* library errors */ -static const int PJD_ERR_GEOCENTRIC = -45; -static const int PJD_ERR_AXIS = -47; -static const int PJD_ERR_GRID_AREA = -48; -static const int PJD_ERR_CATALOG = -49; +enum error_type +{ + error_no_args = -1, + error_no_option_in_init_file = -2, + error_no_colon_in_init_string = -3, + error_proj_not_named = -4, + error_unknown_projection_id = -5, + error_eccentricity_is_one = -6, + error_unknow_unit_id = -7, + error_invalid_boolean_param = -8, + error_unknown_ellp_param = -9, + error_rev_flattening_is_zero = -10, + error_ref_rad_larger_than_90 = -11, + error_es_less_than_zero = -12, + error_major_axis_not_given = -13, + error_lat_or_lon_exceed_limit = -14, + error_invalid_x_or_y = -15, + error_wrong_format_dms_value = -16, + error_non_conv_inv_meri_dist = -17, + error_non_con_inv_phi2 = -18, + error_acos_asin_arg_too_large = -19, + error_tolerance_condition = -20, + error_conic_lat_equal = -21, + error_lat_larger_than_90 = -22, + error_lat1_is_zero = -23, + error_lat_ts_larger_than_90 = -24, + error_control_point_no_dist = -25, + error_no_rotation_proj = -26, + error_w_or_m_zero_or_less = -27, + error_lsat_not_in_range = -28, + error_path_not_in_range = -29, + error_h_less_than_zero = -30, + error_k_less_than_zero = -31, + error_lat_1_or_2_zero_or_90 = -32, + error_lat_0_or_alpha_eq_90 = -33, + error_ellipsoid_use_required = -34, + error_invalid_utm_zone = -35, + error_tcheby_val_out_of_range = -36, + error_failed_to_find_proj = -37, + error_failed_to_load_grid = -38, + error_invalid_m_or_n = -39, + error_n_out_of_range = -40, + error_lat_1_2_unspecified = -41, + error_abs_lat1_eq_abs_lat2 = -42, + error_lat_0_half_pi_from_mean = -43, + error_unparseable_cs_def = -44, + error_geocentric = -45, + error_unknown_prime_meridian = -46, + error_axis = -47, + error_grid_area = -48, + error_invalid_sweep_axis = -49, + error_malformed_pipeline = -50, + error_unit_factor_less_than_0 = -51, + error_invalid_scale = -52, + error_non_convergent = -53, + error_missing_args = -54, + error_lat_0_is_zero = -55, + error_ellipsoidal_unsupported = -56, + error_too_many_inits = -57, + error_invalid_arg = -58 +}; template <typename T> struct pvalue { std::string param; - int used; - - int i; - T f; std::string s; + //int used; }; -template <typename T> -struct pj_const_pod -{ - int over; /* over-range flag */ - int geoc; /* geocentric latitude flag */ - int is_latlong; /* proj=latlong ... not really a projection at all */ - int is_geocent; /* proj=geocent ... not really a projection at all */ - T - a, /* major axis or radius if es==0 */ - a_orig, /* major axis before any +proj related adjustment */ - es, /* e ^ 2 */ - es_orig, /* es before any +proj related adjustment */ - e, /* eccentricity */ - ra, /* 1/A */ - one_es, /* 1 - e^2 */ - rone_es, /* 1/one_es */ - lam0, phi0, /* central longitude, latitude */ - x0, y0, /* easting and northing */ - k0, /* general scaling factor */ - to_meter, fr_meter, /* cartesian scaling */ - vto_meter, vfr_meter; /* Vertical scaling. Internal unit [m] */ - - int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ - T datum_params[7]; - T from_greenwich; /* prime meridian offset (in radians) */ - T long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ - bool is_long_wrap_set; - - // Initialize all variables to zero - pj_const_pod() - { - std::memset(this, 0, sizeof(pj_const_pod)); - } -}; +// Originally defined in proj_internal.h +//enum pj_io_units { +// pj_io_units_whatever = 0, /* Doesn't matter (or depends on pipeline neighbours) */ +// pj_io_units_classic = 1, /* Scaled meters (right), projected system */ +// pj_io_units_projected = 2, /* Meters, projected system */ +// pj_io_units_cartesian = 3, /* Meters, 3D cartesian system */ +// pj_io_units_angular = 4 /* Radians */ +//}; + +// Originally defined in proj_internal.h +/* Maximum latitudinal overshoot accepted */ +//static const double pj_epsilon_lat = 1e-12; template <typename T> -struct pj_const_non_pod +struct pj_consts { - int over; /* over-range flag */ - int geoc; /* geocentric latitude flag */ - int is_latlong; /* proj=latlong ... not really a projection at all */ - int is_geocent; /* proj=geocent ... not really a projection at all */ - T - a, /* major axis or radius if es==0 */ - a_orig, /* major axis before any +proj related adjustment */ - es, /* e ^ 2 */ - es_orig, /* es before any +proj related adjustment */ - e, /* eccentricity */ - ra, /* 1/A */ - one_es, /* 1 - e^2 */ - rone_es, /* 1/one_es */ - lam0, phi0, /* central longitude, latitude */ - x0, y0, /* easting and northing */ - k0, /* general scaling factor */ - to_meter, fr_meter, /* cartesian scaling */ - vto_meter, vfr_meter; /* Vertical scaling. Internal unit [m] */ - - int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ - T datum_params[7]; - T from_greenwich; /* prime meridian offset (in radians) */ - T long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ - bool is_long_wrap_set; - - // Initialize all variables to zero - pj_const_non_pod() - : over(0), geoc(0), is_latlong(0), is_geocent(0) - , a(0), a_orig(0), es(0), es_orig(0), e(0), ra(0) - , one_es(0), rone_es(0), lam0(0), phi0(0), x0(0), y0(0), k0(0) - , to_meter(0), fr_meter(0), vto_meter(0), vfr_meter(0) - , datum_type(PJD_UNKNOWN) - , from_greenwich(0), long_wrap_center(0), is_long_wrap_set(false) - { - datum_params[0] = 0; - datum_params[1] = 0; - datum_params[2] = 0; - datum_params[3] = 0; - datum_params[4] = 0; - datum_params[5] = 0; - datum_params[6] = 0; - } -}; + // E L L I P S O I D P A R A M E T E R S -template <typename T> -struct pj_const - : boost::mpl::if_c - < - boost::is_pod<T>::value, - pj_const_pod<T>, - pj_const_non_pod<T> - >::type -{}; + T a; /* semimajor axis (radius if eccentricity==0) */ + T ra; /* 1/a */ -// PROJ4 complex. Might be replaced with std::complex -template <typename T> -struct COMPLEX { T r, i; }; + T e; /* first eccentricity */ + T es; /* first eccentricity squared */ + T one_es; /* 1 - e^2 */ + T rone_es; /* 1/one_es */ -struct PJ_ELLPS -{ - std::string id; /* ellipse keyword name */ - std::string major; /* a= value */ - std::string ell; /* elliptical parameter */ - std::string name; /* comments */ -}; + T es_orig, a_orig; /* es and a before any +proj related adjustment */ -struct PJ_DATUMS -{ - std::string id; /* datum keyword */ - std::string defn; /* ie. "to_wgs84=..." */ - std::string ellipse_id; /* ie from ellipse table */ - std::string comments; /* EPSG code, etc */ -}; + // C O O R D I N A T E H A N D L I N G -struct PJ_PRIME_MERIDIANS -{ - std::string id; /* prime meridian keyword */ - std::string defn; /* offset from greenwich in DMS format. */ -}; + int over; /* over-range flag */ + int geoc; /* geocentric latitude flag */ + int is_latlong; /* proj=latlong ... not really a projection at all */ + int is_geocent; /* proj=geocent ... not really a projection at all */ + //int need_ellps; /* 0 for operations that are purely cartesian */ -struct PJ_UNITS -{ - std::string id; /* units keyword */ - std::string to_meter; /* multiply by value to get meters */ - std::string name; /* comments */ -}; + //enum pj_io_units left; /* Flags for input/output coordinate types */ + //enum pj_io_units right; -template <typename T> -struct DERIVS -{ - T x_l, x_p; /* derivatives of x for lambda-phi */ - T y_l, y_p; /* derivatives of y for lambda-phi */ + // C A R T O G R A P H I C O F F S E T S + + T lam0, phi0; /* central longitude, latitude */ + T x0, y0/*, z0, t0*/; /* false easting and northing (and height and time) */ + + // S C A L I N G + + T k0; /* general scaling factor */ + T to_meter, fr_meter; /* cartesian scaling */ + T vto_meter, vfr_meter; /* Vertical scaling. Internal unit [m] */ + + // D A T U M S A N D H E I G H T S Y S T E M S + + detail::datum_type datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ + T datum_params[7]; /* Parameters for 3PARAM and 7PARAM */ + + T from_greenwich; /* prime meridian offset (in radians) */ + T long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ + bool is_long_wrap_set; + + // Initialize all variables + pj_consts() + : a(0), ra(0) + , e(0), es(0), one_es(0), rone_es(0) + , es_orig(0), a_orig(0) + , over(0), geoc(0), is_latlong(0), is_geocent(0) + //, need_ellps(1) + //, left(PJ_IO_UNITS_ANGULAR), right(PJ_IO_UNITS_CLASSIC) + , lam0(0), phi0(0) + , x0(0), y0(0)/*, z0(0), t0(0)*/ + , k0(0) , to_meter(0), fr_meter(0), vto_meter(0), vfr_meter(0) + , datum_type(datum_unknown) +#if !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && (!defined(_MSC_VER) || (_MSC_VER >= 1900)) // workaround for VC++ 12 (aka 2013) + , datum_params{0, 0, 0, 0, 0, 0, 0} +#endif + , from_greenwich(0), long_wrap_center(0), is_long_wrap_set(false) + { +#if defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) || (defined(_MSC_VER) && (_MSC_VER < 1900)) // workaround for VC++ 12 (aka 2013) + std::fill(datum_params, datum_params + 7, T(0)); +#endif + } }; +// PROJ4 complex. Might be replaced with std::complex template <typename T> -struct FACTORS -{ - DERIVS<T> der; - T h, k; /* meridinal, parallel scales */ - T omega, thetap; /* angular distortion, theta prime */ - T conv; /* convergence */ - T s; /* areal scale factor */ - T a, b; /* max-min scale error */ - int code; /* info as to analytics, see following */ -}; +struct pj_complex { T r, i; }; } // namespace detail #endif // DOXYGEN_NO_DETAIL @@ -260,7 +233,7 @@ struct FACTORS \ingroup projection */ template <typename T> -struct parameters : public detail::pj_const<T> +struct parameters : public detail::pj_consts<T> { typedef T type; |