geometry/srs/projection.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0005
0006 // This file was modified by Oracle on 2017-2020.
0007 // Modifications copyright (c) 2017-2020, Oracle and/or its affiliates.
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0009
0010 // Use, modification and distribution is subject to the Boost Software License,
0011 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013
0014 #ifndef BOOST_GEOMETRY_SRS_PROJECTION_HPP
0015 #define BOOST_GEOMETRY_SRS_PROJECTION_HPP
0016
0017
0018 #include <memory>
0019 #include <string>
0020 #include <type_traits>
0021
0022 #include <boost/range/size.hpp>
0023 #include <boost/throw_exception.hpp>
0024
0025 #include <boost/geometry/algorithms/convert.hpp>
0026 #include <boost/geometry/algorithms/detail/convert_point_to_point.hpp>
0027
0028 #include <boost/geometry/core/coordinate_dimension.hpp>
0029 #include <boost/geometry/core/static_assert.hpp>
0030
0031 #include <boost/geometry/srs/projections/dpar.hpp>
0032 #include <boost/geometry/srs/projections/exception.hpp>
0033 #include <boost/geometry/srs/projections/factory.hpp>
0034 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0035 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0036 #include <boost/geometry/srs/projections/impl/pj_init.hpp>
0037 #include <boost/geometry/srs/projections/invalid_point.hpp>
0038 #include <boost/geometry/srs/projections/proj4.hpp>
0039 #include <boost/geometry/srs/projections/spar.hpp>
0040
0041 #include <boost/geometry/views/detail/indexed_point_view.hpp>
0042
0043
0044 namespace boost { namespace geometry
0045 {
0046
0047 namespace projections
0048 {
0049
0050 #ifndef DOXYGEN_NO_DETAIL
0051 namespace detail
0052 {
0053
0054 template <typename G1, typename G2>
0055 struct same_tags
0056     : std::is_same<geometry::tag_t<G1>, geometry::tag_t<G2>>
0057 {};
0058
0059 template <typename CT>
0060 struct promote_to_double
0061 {
0062     typedef std::conditional_t
0063         <
0064             std::is_integral<CT>::value || std::is_same<CT, float>::value,
0065             double, CT
0066         > type;
0067 };
0068
0069 // Copy coordinates of dimensions >= MinDim
0070 template <std::size_t MinDim, typename Point1, typename Point2>
0071 inline void copy_higher_dimensions(Point1 const& point1, Point2 & point2)
0072 {
0073     static const std::size_t dim1 = geometry::dimension<Point1>::value;
0074     static const std::size_t dim2 = geometry::dimension<Point2>::value;
0075     static const std::size_t lesser_dim = dim1 < dim2 ? dim1 : dim2;
0076     BOOST_GEOMETRY_STATIC_ASSERT((lesser_dim >= MinDim),
0077         "The dimension of Point1 or Point2 is too small.",
0078         Point1, Point2);
0079
0080     geometry::detail::conversion::point_to_point
0081         <
0082             Point1, Point2, MinDim, lesser_dim
0083         > ::apply(point1, point2);
0084
0085     // TODO: fill point2 with zeros if dim1 < dim2 ?
0086     // currently no need because equal dimensions are checked
0087 }
0088
0089
0090 struct forward_point_projection_policy
0091 {
0092     template <typename LL, typename XY, typename Proj>
0093     static inline bool apply(LL const& ll, XY & xy, Proj const& proj)
0094     {
0095         return proj.forward(ll, xy);
0096     }
0097 };
0098
0099 struct inverse_point_projection_policy
0100 {
0101     template <typename XY, typename LL, typename Proj>
0102     static inline bool apply(XY const& xy, LL & ll, Proj const& proj)
0103     {
0104         return proj.inverse(xy, ll);
0105     }
0106 };
0107
0108 template <typename PointPolicy>
0109 struct project_point
0110 {
0111     template <typename P1, typename P2, typename Proj>
0112     static inline bool apply(P1 const& p1, P2 & p2, Proj const& proj)
0113     {
0114         // (Geographic -> Cartesian) will be projected, rest will be copied.
0115         // So first copy third or higher dimensions
0116         projections::detail::copy_higher_dimensions<2>(p1, p2);
0117
0118         if (! PointPolicy::apply(p1, p2, proj))
0119         {
0120             // For consistency with transformation
0121             set_invalid_point(p2);
0122             return false;
0123         }
0124
0125         return true;
0126     }
0127 };
0128
0129 template <typename PointPolicy>
0130 struct project_range
0131 {
0132     template <typename Proj>
0133     struct convert_policy
0134     {
0135         explicit convert_policy(Proj const& proj)
0136             : m_proj(proj)
0137             , m_result(true)
0138         {}
0139
0140         template <typename Point1, typename Point2>
0141         inline void apply(Point1 const& point1, Point2 & point2)
0142         {
0143             if (! project_point<PointPolicy>::apply(point1, point2, m_proj) )
0144                 m_result = false;
0145         }
0146
0147         bool result() const
0148         {
0149             return m_result;
0150         }
0151
0152     private:
0153         Proj const& m_proj;
0154         bool m_result;
0155     };
0156
0157     template <typename R1, typename R2, typename Proj>
0158     static inline bool apply(R1 const& r1, R2 & r2, Proj const& proj)
0159     {
0160         return geometry::detail::conversion::range_to_range::apply(r1, r2,
0161             convert_policy<Proj>(proj)).result();
0162     }
0163 };
0164
0165 template <typename Policy>
0166 struct project_multi
0167 {
0168     template <typename G1, typename G2, typename Proj>
0169     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0170     {
0171         range::resize(g2, boost::size(g1));
0172         return apply(boost::begin(g1), boost::end(g1),
0173                      boost::begin(g2),
0174                      proj);
0175     }
0176
0177 private:
0178     template <typename It1, typename It2, typename Proj>
0179     static inline bool apply(It1 g1_first, It1 g1_last, It2 g2_first, Proj const& proj)
0180     {
0181         bool result = true;
0182         for ( ; g1_first != g1_last ; ++g1_first, ++g2_first )
0183         {
0184             if (! Policy::apply(*g1_first, *g2_first, proj))
0185             {
0186                 result = false;
0187             }
0188         }
0189         return result;
0190     }
0191 };
0192
0193 template
0194 <
0195     typename Geometry,
0196     typename PointPolicy,
0197     typename Tag = geometry::tag_t<Geometry>
0198 >
0199 struct project_geometry
0200 {};
0201
0202 template <typename Geometry, typename PointPolicy>
0203 struct project_geometry<Geometry, PointPolicy, point_tag>
0204     : project_point<PointPolicy>
0205 {};
0206
0207 template <typename Geometry, typename PointPolicy>
0208 struct project_geometry<Geometry, PointPolicy, multi_point_tag>
0209     : project_range<PointPolicy>
0210 {};
0211
0212 template <typename Geometry, typename PointPolicy>
0213 struct project_geometry<Geometry, PointPolicy, segment_tag>
0214 {
0215     template <typename G1, typename G2, typename Proj>
0216     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0217     {
0218         bool r1 = apply<0>(g1, g2, proj);
0219         bool r2 = apply<1>(g1, g2, proj);
0220         return r1 && r2;
0221     }
0222
0223 private:
0224     template <std::size_t Index, typename G1, typename G2, typename Proj>
0225     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0226     {
0227         geometry::detail::indexed_point_view<G1 const, Index> pt1(g1);
0228         geometry::detail::indexed_point_view<G2, Index> pt2(g2);
0229         return project_point<PointPolicy>::apply(pt1, pt2, proj);
0230     }
0231 };
0232
0233 template <typename Geometry, typename PointPolicy>
0234 struct project_geometry<Geometry, PointPolicy, linestring_tag>
0235     : project_range<PointPolicy>
0236 {};
0237
0238 template <typename Geometry, typename PointPolicy>
0239 struct project_geometry<Geometry, PointPolicy, multi_linestring_tag>
0240     : project_multi< project_range<PointPolicy> >
0241 {};
0242
0243 template <typename Geometry, typename PointPolicy>
0244 struct project_geometry<Geometry, PointPolicy, ring_tag>
0245     : project_range<PointPolicy>
0246 {};
0247
0248 template <typename Geometry, typename PointPolicy>
0249 struct project_geometry<Geometry, PointPolicy, polygon_tag>
0250 {
0251     template <typename G1, typename G2, typename Proj>
0252     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0253     {
0254         bool r1 = project_range
0255                     <
0256                         PointPolicy
0257                     >::apply(geometry::exterior_ring(g1),
0258                              geometry::exterior_ring(g2),
0259                              proj);
0260         bool r2 = project_multi
0261                     <
0262                         project_range<PointPolicy>
0263                     >::apply(geometry::interior_rings(g1),
0264                              geometry::interior_rings(g2),
0265                              proj);
0266         return r1 && r2;
0267     }
0268 };
0269
0270 template <typename MultiPolygon, typename PointPolicy>
0271 struct project_geometry<MultiPolygon, PointPolicy, multi_polygon_tag>
0272     : project_multi
0273         <
0274             project_geometry
0275             <
0276                 typename boost::range_value<MultiPolygon>::type,
0277                 PointPolicy,
0278                 polygon_tag
0279             >
0280         >
0281 {};
0282
0283
0284 } // namespace detail
0285 #endif // DOXYGEN_NO_DETAIL
0286
0287
0288 template <typename Params>
0289 struct dynamic_parameters
0290 {
0291     static const bool is_specialized = false;
0292 };
0293
0294 template <>
0295 struct dynamic_parameters<srs::proj4>
0296 {
0297     static const bool is_specialized = true;
0298     static inline srs::detail::proj4_parameters apply(srs::proj4 const& params)
0299     {
0300         return srs::detail::proj4_parameters(params.str());
0301     }
0302 };
0303
0304 template <typename T>
0305 struct dynamic_parameters<srs::dpar::parameters<T> >
0306 {
0307     static const bool is_specialized = true;
0308     static inline srs::dpar::parameters<T> const& apply(srs::dpar::parameters<T> const& params)
0309     {
0310         return params;
0311     }
0312 };
0313
0314
0315 // proj_wrapper class and its specializations wrapps the internal projection
0316 // representation and implements transparent creation of projection object
0317 template <typename Proj, typename CT>
0318 class proj_wrapper
0319 {
0320     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
0321         "Unknown projection definition.",
0322         Proj);
0323 };
0324
0325 template <typename CT>
0326 class proj_wrapper<srs::dynamic, CT>
0327 {
0328     // Some projections do not work with float -> wrong results
0329     // select <double> from int/float/double and else selects T
0330     typedef typename projections::detail::promote_to_double<CT>::type calc_t;
0331
0332     typedef projections::parameters<calc_t> parameters_type;
0333     typedef projections::detail::dynamic_wrapper_b<calc_t, parameters_type> vprj_t;
0334
0335 public:
0336     template
0337     <
0338         typename Params,
0339         std::enable_if_t
0340             <
0341                 dynamic_parameters<Params>::is_specialized,
0342                 int
0343             > = 0
0344     >
0345     proj_wrapper(Params const& params)
0346         : m_ptr(create(dynamic_parameters<Params>::apply(params)))
0347     {}
0348
0349     vprj_t const& proj() const { return *m_ptr; }
0350     vprj_t & mutable_proj() { return *m_ptr; }
0351
0352 private:
0353     template <typename Params>
0354     static vprj_t* create(Params const& params)
0355     {
0356         parameters_type parameters = projections::detail::pj_init<calc_t>(params);
0357
0358         vprj_t* result = projections::detail::create_new(params, parameters);
0359
0360         if (result == NULL)
0361         {
0362             if (parameters.id.is_unknown())
0363             {
0364                 BOOST_THROW_EXCEPTION(projection_not_named_exception());
0365             }
0366             else
0367             {
0368                 // TODO: handle non-string projection id
0369                 BOOST_THROW_EXCEPTION(projection_unknown_id_exception());
0370             }
0371         }
0372
0373         return result;
0374     }
0375
0376     std::shared_ptr<vprj_t> m_ptr;
0377 };
0378
0379 template <typename StaticParameters, typename CT>
0380 class static_proj_wrapper_base
0381 {
0382     typedef typename projections::detail::promote_to_double<CT>::type calc_t;
0383
0384     typedef projections::parameters<calc_t> parameters_type;
0385
0386     typedef typename srs::spar::detail::pick_proj_tag
0387         <
0388             StaticParameters
0389         >::type proj_tag;
0390
0391     typedef typename projections::detail::static_projection_type
0392         <
0393             proj_tag,
0394             typename projections::detail::static_srs_tag<StaticParameters>::type,
0395             StaticParameters,
0396             calc_t,
0397             parameters_type
0398         >::type projection_type;
0399
0400 public:
0401     projection_type const& proj() const { return m_proj; }
0402     projection_type & mutable_proj() { return m_proj; }
0403
0404 protected:
0405     explicit static_proj_wrapper_base(StaticParameters const& s_params)
0406         : m_proj(s_params,
0407                  projections::detail::pj_init<calc_t>(s_params))
0408     {}
0409
0410 private:
0411     projection_type m_proj;
0412 };
0413
0414 template <typename ...Ps, typename CT>
0415 class proj_wrapper<srs::spar::parameters<Ps...>, CT>
0416     : public static_proj_wrapper_base<srs::spar::parameters<Ps...>, CT>
0417 {
0418     typedef srs::spar::parameters<Ps...>
0419         static_parameters_type;
0420     typedef static_proj_wrapper_base
0421         <
0422             static_parameters_type,
0423             CT
0424         > base_t;
0425
0426 public:
0427     proj_wrapper()
0428         : base_t(static_parameters_type())
0429     {}
0430
0431     proj_wrapper(static_parameters_type const& s_params)
0432         : base_t(s_params)
0433     {}
0434 };
0435
0436
0437 // projection class implements transparent forward/inverse projection interface
0438 template <typename Proj, typename CT>
0439 class projection
0440     : private proj_wrapper<Proj, CT>
0441 {
0442     typedef proj_wrapper<Proj, CT> base_t;
0443
0444 public:
0445     projection()
0446     {}
0447
0448     template <typename Params>
0449     explicit projection(Params const& params)
0450         : base_t(params)
0451     {}
0452
0453     /// Forward projection, from Latitude-Longitude to Cartesian
0454     template <typename LL, typename XY>
0455     inline bool forward(LL const& ll, XY& xy) const
0456     {
0457         BOOST_GEOMETRY_STATIC_ASSERT(
0458             (projections::detail::same_tags<LL, XY>::value),
0459             "Not supported combination of Geometries.",
0460             LL, XY);
0461
0462         concepts::check_concepts_and_equal_dimensions<LL const, XY>();
0463
0464         return projections::detail::project_geometry
0465                 <
0466                     LL,
0467                     projections::detail::forward_point_projection_policy
0468                 >::apply(ll, xy, base_t::proj());
0469     }
0470
0471     /// Inverse projection, from Cartesian to Latitude-Longitude
0472     template <typename XY, typename LL>
0473     inline bool inverse(XY const& xy, LL& ll) const
0474     {
0475         BOOST_GEOMETRY_STATIC_ASSERT(
0476             (projections::detail::same_tags<XY, LL>::value),
0477             "Not supported combination of Geometries.",
0478             XY, LL);
0479
0480         concepts::check_concepts_and_equal_dimensions<XY const, LL>();
0481
0482         return projections::detail::project_geometry
0483                 <
0484                     XY,
0485                     projections::detail::inverse_point_projection_policy
0486                 >::apply(xy, ll, base_t::proj());
0487     }
0488 };
0489
0490 } // namespace projections
0491
0492
0493 namespace srs
0494 {
0495
0496
0497 /*!
0498     \brief Representation of projection
0499     \details Either dynamic or static projection representation
0500     \ingroup projection
0501     \tparam Parameters default dynamic tag or static projection parameters
0502     \tparam CT calculation type used internally
0503 */
0504 template
0505 <
0506     typename Parameters = srs::dynamic,
0507     typename CT = double
0508 >
0509 class projection
0510     : public projections::projection<Parameters, CT>
0511 {
0512     typedef projections::projection<Parameters, CT> base_t;
0513
0514 public:
0515     projection()
0516     {}
0517
0518     projection(Parameters const& parameters)
0519         : base_t(parameters)
0520     {}
0521
0522     /*!
0523     \ingroup projection
0524     \brief Initializes a projection as a string, using the format with + and =
0525     \details The projection can be initialized with a string (with the same format as the PROJ4 package) for
0526       convenient initialization from, for example, the command line
0527     \par Example
0528         <tt>srs::proj4("+proj=labrd +ellps=intl +lon_0=46d26'13.95E +lat_0=18d54S +azi=18d54 +k_0=.9995 +x_0=400000 +y_0=800000")</tt>
0529         for the Madagascar projection.
0530     */
0531     template
0532     <
0533         typename DynamicParameters,
0534         std::enable_if_t
0535             <
0536                 projections::dynamic_parameters<DynamicParameters>::is_specialized,
0537                 int
0538             > = 0
0539     >
0540     projection(DynamicParameters const& dynamic_parameters)
0541         : base_t(dynamic_parameters)
0542     {}
0543 };
0544
0545
0546 } // namespace srs
0547
0548
0549 }} // namespace boost::geometry
0550
0551
0552 #endif // BOOST_GEOMETRY_SRS_PROJECTION_HPP