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