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
0161             <
0162                 R1, R2,
0163                 geometry::point_order<R1>::value != geometry::point_order<R2>::value
0164             >::apply(r1, r2, convert_policy<Proj>(proj)).result();
0165     }
0166 };
0167
0168 template <typename Policy>
0169 struct project_multi
0170 {
0171     template <typename G1, typename G2, typename Proj>
0172     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0173     {
0174         range::resize(g2, boost::size(g1));
0175         return apply(boost::begin(g1), boost::end(g1),
0176                      boost::begin(g2),
0177                      proj);
0178     }
0179
0180 private:
0181     template <typename It1, typename It2, typename Proj>
0182     static inline bool apply(It1 g1_first, It1 g1_last, It2 g2_first, Proj const& proj)
0183     {
0184         bool result = true;
0185         for ( ; g1_first != g1_last ; ++g1_first, ++g2_first )
0186         {
0187             if (! Policy::apply(*g1_first, *g2_first, proj))
0188             {
0189                 result = false;
0190             }
0191         }
0192         return result;
0193     }
0194 };
0195
0196 template
0197 <
0198     typename Geometry,
0199     typename PointPolicy,
0200     typename Tag = geometry::tag_t<Geometry>
0201 >
0202 struct project_geometry
0203 {};
0204
0205 template <typename Geometry, typename PointPolicy>
0206 struct project_geometry<Geometry, PointPolicy, point_tag>
0207     : project_point<PointPolicy>
0208 {};
0209
0210 template <typename Geometry, typename PointPolicy>
0211 struct project_geometry<Geometry, PointPolicy, multi_point_tag>
0212     : project_range<PointPolicy>
0213 {};
0214
0215 template <typename Geometry, typename PointPolicy>
0216 struct project_geometry<Geometry, PointPolicy, segment_tag>
0217 {
0218     template <typename G1, typename G2, typename Proj>
0219     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0220     {
0221         bool r1 = apply<0>(g1, g2, proj);
0222         bool r2 = apply<1>(g1, g2, proj);
0223         return r1 && r2;
0224     }
0225
0226 private:
0227     template <std::size_t Index, typename G1, typename G2, typename Proj>
0228     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0229     {
0230         geometry::detail::indexed_point_view<G1 const, Index> pt1(g1);
0231         geometry::detail::indexed_point_view<G2, Index> pt2(g2);
0232         return project_point<PointPolicy>::apply(pt1, pt2, proj);
0233     }
0234 };
0235
0236 template <typename Geometry, typename PointPolicy>
0237 struct project_geometry<Geometry, PointPolicy, linestring_tag>
0238     : project_range<PointPolicy>
0239 {};
0240
0241 template <typename Geometry, typename PointPolicy>
0242 struct project_geometry<Geometry, PointPolicy, multi_linestring_tag>
0243     : project_multi< project_range<PointPolicy> >
0244 {};
0245
0246 template <typename Geometry, typename PointPolicy>
0247 struct project_geometry<Geometry, PointPolicy, ring_tag>
0248     : project_range<PointPolicy>
0249 {};
0250
0251 template <typename Geometry, typename PointPolicy>
0252 struct project_geometry<Geometry, PointPolicy, polygon_tag>
0253 {
0254     template <typename G1, typename G2, typename Proj>
0255     static inline bool apply(G1 const& g1, G2 & g2, Proj const& proj)
0256     {
0257         bool r1 = project_range
0258                     <
0259                         PointPolicy
0260                     >::apply(geometry::exterior_ring(g1),
0261                              geometry::exterior_ring(g2),
0262                              proj);
0263         bool r2 = project_multi
0264                     <
0265                         project_range<PointPolicy>
0266                     >::apply(geometry::interior_rings(g1),
0267                              geometry::interior_rings(g2),
0268                              proj);
0269         return r1 && r2;
0270     }
0271 };
0272
0273 template <typename MultiPolygon, typename PointPolicy>
0274 struct project_geometry<MultiPolygon, PointPolicy, multi_polygon_tag>
0275     : project_multi
0276         <
0277             project_geometry
0278             <
0279                 typename boost::range_value<MultiPolygon>::type,
0280                 PointPolicy,
0281                 polygon_tag
0282             >
0283         >
0284 {};
0285
0286
0287 } // namespace detail
0288 #endif // DOXYGEN_NO_DETAIL
0289
0290
0291 template <typename Params>
0292 struct dynamic_parameters
0293 {
0294     static const bool is_specialized = false;
0295 };
0296
0297 template <>
0298 struct dynamic_parameters<srs::proj4>
0299 {
0300     static const bool is_specialized = true;
0301     static inline srs::detail::proj4_parameters apply(srs::proj4 const& params)
0302     {
0303         return srs::detail::proj4_parameters(params.str());
0304     }
0305 };
0306
0307 template <typename T>
0308 struct dynamic_parameters<srs::dpar::parameters<T> >
0309 {
0310     static const bool is_specialized = true;
0311     static inline srs::dpar::parameters<T> const& apply(srs::dpar::parameters<T> const& params)
0312     {
0313         return params;
0314     }
0315 };
0316
0317
0318 // proj_wrapper class and its specializations wrapps the internal projection
0319 // representation and implements transparent creation of projection object
0320 template <typename Proj, typename CT>
0321 class proj_wrapper
0322 {
0323     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
0324         "Unknown projection definition.",
0325         Proj);
0326 };
0327
0328 template <typename CT>
0329 class proj_wrapper<srs::dynamic, CT>
0330 {
0331     // Some projections do not work with float -> wrong results
0332     // select <double> from int/float/double and else selects T
0333     typedef typename projections::detail::promote_to_double<CT>::type calc_t;
0334
0335     typedef projections::parameters<calc_t> parameters_type;
0336     typedef projections::detail::dynamic_wrapper_b<calc_t, parameters_type> vprj_t;
0337
0338 public:
0339     template
0340     <
0341         typename Params,
0342         std::enable_if_t
0343             <
0344                 dynamic_parameters<Params>::is_specialized,
0345                 int
0346             > = 0
0347     >
0348     proj_wrapper(Params const& params)
0349         : m_ptr(create(dynamic_parameters<Params>::apply(params)))
0350     {}
0351
0352     vprj_t const& proj() const { return *m_ptr; }
0353     vprj_t & mutable_proj() { return *m_ptr; }
0354
0355 private:
0356     template <typename Params>
0357     static vprj_t* create(Params const& params)
0358     {
0359         parameters_type parameters = projections::detail::pj_init<calc_t>(params);
0360
0361         vprj_t* result = projections::detail::create_new(params, parameters);
0362
0363         if (result == NULL)
0364         {
0365             if (parameters.id.is_unknown())
0366             {
0367                 BOOST_THROW_EXCEPTION(projection_not_named_exception());
0368             }
0369             else
0370             {
0371                 // TODO: handle non-string projection id
0372                 BOOST_THROW_EXCEPTION(projection_unknown_id_exception());
0373             }
0374         }
0375
0376         return result;
0377     }
0378
0379     std::shared_ptr<vprj_t> m_ptr;
0380 };
0381
0382 template <typename StaticParameters, typename CT>
0383 class static_proj_wrapper_base
0384 {
0385     typedef typename projections::detail::promote_to_double<CT>::type calc_t;
0386
0387     typedef projections::parameters<calc_t> parameters_type;
0388
0389     typedef typename srs::spar::detail::pick_proj_tag
0390         <
0391             StaticParameters
0392         >::type proj_tag;
0393
0394     typedef typename projections::detail::static_projection_type
0395         <
0396             proj_tag,
0397             typename projections::detail::static_srs_tag<StaticParameters>::type,
0398             StaticParameters,
0399             calc_t,
0400             parameters_type
0401         >::type projection_type;
0402
0403 public:
0404     projection_type const& proj() const { return m_proj; }
0405     projection_type & mutable_proj() { return m_proj; }
0406
0407 protected:
0408     explicit static_proj_wrapper_base(StaticParameters const& s_params)
0409         : m_proj(s_params,
0410                  projections::detail::pj_init<calc_t>(s_params))
0411     {}
0412
0413 private:
0414     projection_type m_proj;
0415 };
0416
0417 template <typename ...Ps, typename CT>
0418 class proj_wrapper<srs::spar::parameters<Ps...>, CT>
0419     : public static_proj_wrapper_base<srs::spar::parameters<Ps...>, CT>
0420 {
0421     typedef srs::spar::parameters<Ps...>
0422         static_parameters_type;
0423     typedef static_proj_wrapper_base
0424         <
0425             static_parameters_type,
0426             CT
0427         > base_t;
0428
0429 public:
0430     proj_wrapper()
0431         : base_t(static_parameters_type())
0432     {}
0433
0434     proj_wrapper(static_parameters_type const& s_params)
0435         : base_t(s_params)
0436     {}
0437 };
0438
0439
0440 // projection class implements transparent forward/inverse projection interface
0441 template <typename Proj, typename CT>
0442 class projection
0443     : private proj_wrapper<Proj, CT>
0444 {
0445     typedef proj_wrapper<Proj, CT> base_t;
0446
0447 public:
0448     projection()
0449     {}
0450
0451     template <typename Params>
0452     explicit projection(Params const& params)
0453         : base_t(params)
0454     {}
0455
0456     /// Forward projection, from Latitude-Longitude to Cartesian
0457     template <typename LL, typename XY>
0458     inline bool forward(LL const& ll, XY& xy) const
0459     {
0460         BOOST_GEOMETRY_STATIC_ASSERT(
0461             (projections::detail::same_tags<LL, XY>::value),
0462             "Not supported combination of Geometries.",
0463             LL, XY);
0464
0465         concepts::check_concepts_and_equal_dimensions<LL const, XY>();
0466
0467         return projections::detail::project_geometry
0468                 <
0469                     LL,
0470                     projections::detail::forward_point_projection_policy
0471                 >::apply(ll, xy, base_t::proj());
0472     }
0473
0474     /// Inverse projection, from Cartesian to Latitude-Longitude
0475     template <typename XY, typename LL>
0476     inline bool inverse(XY const& xy, LL& ll) const
0477     {
0478         BOOST_GEOMETRY_STATIC_ASSERT(
0479             (projections::detail::same_tags<XY, LL>::value),
0480             "Not supported combination of Geometries.",
0481             XY, LL);
0482
0483         concepts::check_concepts_and_equal_dimensions<XY const, LL>();
0484
0485         return projections::detail::project_geometry
0486                 <
0487                     XY,
0488                     projections::detail::inverse_point_projection_policy
0489                 >::apply(xy, ll, base_t::proj());
0490     }
0491 };
0492
0493 } // namespace projections
0494
0495
0496 namespace srs
0497 {
0498
0499
0500 /*!
0501     \brief Representation of projection
0502     \details Either dynamic or static projection representation
0503     \ingroup projection
0504     \tparam Parameters default dynamic tag or static projection parameters
0505     \tparam CT calculation type used internally
0506 */
0507 template
0508 <
0509     typename Parameters = srs::dynamic,
0510     typename CT = double
0511 >
0512 class projection
0513     : public projections::projection<Parameters, CT>
0514 {
0515     typedef projections::projection<Parameters, CT> base_t;
0516
0517 public:
0518     projection()
0519     {}
0520
0521     projection(Parameters const& parameters)
0522         : base_t(parameters)
0523     {}
0524
0525     /*!
0526     \ingroup projection
0527     \brief Initializes a projection as a string, using the format with + and =
0528     \details The projection can be initialized with a string (with the same format as the PROJ4 package) for
0529       convenient initialization from, for example, the command line
0530     \par Example
0531         <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>
0532         for the Madagascar projection.
0533     */
0534     template
0535     <
0536         typename DynamicParameters,
0537         std::enable_if_t
0538             <
0539                 projections::dynamic_parameters<DynamicParameters>::is_specialized,
0540                 int
0541             > = 0
0542     >
0543     projection(DynamicParameters const& dynamic_parameters)
0544         : base_t(dynamic_parameters)
0545     {}
0546 };
0547
0548
0549 } // namespace srs
0550
0551
0552 }} // namespace boost::geometry
0553
0554
0555 #endif // BOOST_GEOMETRY_SRS_PROJECTION_HPP