Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:36:41

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002 
0003 // Copyright (c) 2017-2020, Oracle and/or its affiliates.
0004 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0005 
0006 // Use, modification and distribution is subject to the Boost Software License,
0007 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0008 // http://www.boost.org/LICENSE_1_0.txt)
0009 
0010 #ifndef BOOST_GEOMETRY_SRS_TRANSFORMATION_HPP
0011 #define BOOST_GEOMETRY_SRS_TRANSFORMATION_HPP
0012 
0013 
0014 #include <string>
0015 #include <type_traits>
0016 
0017 #include <boost/range/size.hpp>
0018 #include <boost/throw_exception.hpp>
0019 
0020 #include <boost/geometry/algorithms/convert.hpp>
0021 
0022 #include <boost/geometry/core/coordinate_dimension.hpp>
0023 #include <boost/geometry/core/static_assert.hpp>
0024 
0025 #include <boost/geometry/geometries/point.hpp>
0026 #include <boost/geometry/geometries/multi_point.hpp>
0027 #include <boost/geometry/geometries/linestring.hpp>
0028 #include <boost/geometry/geometries/ring.hpp>
0029 #include <boost/geometry/geometries/segment.hpp>
0030 
0031 #include <boost/geometry/srs/projection.hpp>
0032 #include <boost/geometry/srs/projections/grids.hpp>
0033 #include <boost/geometry/srs/projections/impl/pj_transform.hpp>
0034 
0035 #include <boost/geometry/views/detail/indexed_point_view.hpp>
0036 
0037 
0038 namespace boost { namespace geometry
0039 {
0040 
0041 namespace projections { namespace detail
0042 {
0043 
0044 template <typename T1, typename T2>
0045 inline bool same_object(T1 const& , T2 const& )
0046 {
0047     return false;
0048 }
0049 
0050 template <typename T>
0051 inline bool same_object(T const& o1, T const& o2)
0052 {
0053     return boost::addressof(o1) == boost::addressof(o2);
0054 }
0055 
0056 template
0057 <
0058     typename PtIn,
0059     typename PtOut,
0060     bool SameUnits = std::is_same
0061                         <
0062                             typename geometry::detail::cs_angular_units<PtIn>::type,
0063                             typename geometry::detail::cs_angular_units<PtOut>::type
0064                         >::value
0065 >
0066 struct transform_geometry_point_coordinates
0067 {
0068     static inline void apply(PtIn const& in, PtOut & out, bool /*enable_angles*/)
0069     {
0070         geometry::set<0>(out, geometry::get<0>(in));
0071         geometry::set<1>(out, geometry::get<1>(in));
0072     }
0073 };
0074 
0075 template <typename PtIn, typename PtOut>
0076 struct transform_geometry_point_coordinates<PtIn, PtOut, false>
0077 {
0078     static inline void apply(PtIn const& in, PtOut & out, bool enable_angles)
0079     {
0080         if (enable_angles)
0081         {
0082             geometry::set_from_radian<0>(out, geometry::get_as_radian<0>(in));
0083             geometry::set_from_radian<1>(out, geometry::get_as_radian<1>(in));
0084         }
0085         else
0086         {
0087             geometry::set<0>(out, geometry::get<0>(in));
0088             geometry::set<1>(out, geometry::get<1>(in));
0089         }
0090     }
0091 };
0092 
0093 template <typename Geometry, typename CT>
0094 struct transform_geometry_point
0095 {
0096     typedef typename geometry::point_type<Geometry>::type point_type;
0097 
0098     typedef geometry::model::point
0099         <
0100             typename select_most_precise
0101                 <
0102                     typename geometry::coordinate_type<point_type>::type,
0103                     CT
0104                 >::type,
0105             geometry::dimension<point_type>::type::value,
0106             typename geometry::coordinate_system<point_type>::type
0107         > type;
0108 
0109     template <typename PtIn, typename PtOut>
0110     static inline void apply(PtIn const& in, PtOut & out, bool enable_angles)
0111     {
0112         transform_geometry_point_coordinates<PtIn, PtOut>::apply(in, out, enable_angles);
0113         projections::detail::copy_higher_dimensions<2>(in, out);
0114     }
0115 };
0116 
0117 template <typename Geometry, typename CT>
0118 struct transform_geometry_range_base
0119 {
0120     struct convert_strategy
0121     {
0122         convert_strategy(bool enable_angles)
0123             : m_enable_angles(enable_angles)
0124         {}
0125 
0126         template <typename PtIn, typename PtOut>
0127         void apply(PtIn const& in, PtOut & out)
0128         {
0129             transform_geometry_point<Geometry, CT>::apply(in, out, m_enable_angles);
0130         }
0131 
0132         bool m_enable_angles;
0133     };
0134 
0135     template <typename In, typename Out>
0136     static inline void apply(In const& in, Out & out, bool enable_angles)
0137     {
0138         // Change the order and/or closure if needed
0139         // In - arbitrary geometry
0140         // Out - either Geometry or std::vector
0141         // So the order and closure of In and Geometry shoudl be compared
0142         // std::vector's order is assumed to be the same as of Geometry
0143         geometry::detail::conversion::range_to_range
0144             <
0145                 In,
0146                 Out,
0147                 geometry::point_order<In>::value != geometry::point_order<Out>::value
0148             >::apply(in, out, convert_strategy(enable_angles));
0149     }
0150 };
0151 
0152 template
0153 <
0154     typename Geometry,
0155     typename CT,
0156     typename Tag = typename geometry::tag<Geometry>::type
0157 >
0158 struct transform_geometry
0159 {};
0160 
0161 template <typename Point, typename CT>
0162 struct transform_geometry<Point, CT, point_tag>
0163     : transform_geometry_point<Point, CT>
0164 {};
0165 
0166 template <typename Segment, typename CT>
0167 struct transform_geometry<Segment, CT, segment_tag>
0168 {
0169     typedef geometry::model::segment
0170         <
0171             typename transform_geometry_point<Segment, CT>::type
0172         > type;
0173 
0174     template <typename In, typename Out>
0175     static inline void apply(In const& in, Out & out, bool enable_angles)
0176     {
0177         apply<0>(in, out, enable_angles);
0178         apply<1>(in, out, enable_angles);
0179     }
0180 
0181 private:
0182     template <std::size_t Index, typename In, typename Out>
0183     static inline void apply(In const& in, Out & out, bool enable_angles)
0184     {
0185         geometry::detail::indexed_point_view<In const, Index> in_pt(in);
0186         geometry::detail::indexed_point_view<Out, Index> out_pt(out);
0187         transform_geometry_point<Segment, CT>::apply(in_pt, out_pt, enable_angles);
0188     }
0189 };
0190 
0191 template <typename MultiPoint, typename CT>
0192 struct transform_geometry<MultiPoint, CT, multi_point_tag>
0193     : transform_geometry_range_base<MultiPoint, CT>
0194 {
0195     typedef model::multi_point
0196         <
0197             typename transform_geometry_point<MultiPoint, CT>::type
0198         > type;
0199 };
0200 
0201 template <typename LineString, typename CT>
0202 struct transform_geometry<LineString, CT, linestring_tag>
0203     : transform_geometry_range_base<LineString, CT>
0204 {
0205     typedef model::linestring
0206         <
0207             typename transform_geometry_point<LineString, CT>::type
0208         > type;
0209 };
0210 
0211 template <typename Ring, typename CT>
0212 struct transform_geometry<Ring, CT, ring_tag>
0213     : transform_geometry_range_base<Ring, CT>
0214 {
0215     typedef model::ring
0216         <
0217             typename transform_geometry_point<Ring, CT>::type,
0218             geometry::point_order<Ring>::value == clockwise,
0219             geometry::closure<Ring>::value == closed
0220         > type;
0221 };
0222 
0223 
0224 template
0225 <
0226     typename OutGeometry,
0227     typename CT,
0228     bool EnableTemporary = ! std::is_same
0229                                 <
0230                                     typename select_most_precise
0231                                         <
0232                                             typename geometry::coordinate_type<OutGeometry>::type,
0233                                             CT
0234                                         >::type,
0235                                     typename geometry::coordinate_type<OutGeometry>::type
0236                                 >::value
0237 >
0238 struct transform_geometry_wrapper
0239 {
0240     typedef transform_geometry<OutGeometry, CT> transform;
0241     typedef typename transform::type type;
0242 
0243     template <typename InGeometry>
0244     transform_geometry_wrapper(InGeometry const& in, OutGeometry & out, bool input_angles)
0245         : m_out(out)
0246     {
0247         transform::apply(in, m_temp, input_angles);
0248     }
0249 
0250     type & get() { return m_temp; }
0251     void finish() { geometry::convert(m_temp, m_out); } // this is always copy 1:1 without changing the order or closure
0252 
0253 private:
0254     type m_temp;
0255     OutGeometry & m_out;
0256 };
0257 
0258 template
0259 <
0260     typename OutGeometry,
0261     typename CT
0262 >
0263 struct transform_geometry_wrapper<OutGeometry, CT, false>
0264 {
0265     typedef transform_geometry<OutGeometry, CT> transform;
0266     typedef OutGeometry type;
0267 
0268     transform_geometry_wrapper(OutGeometry const& in, OutGeometry & out, bool input_angles)
0269         : m_out(out)
0270     {
0271         if (! same_object(in, out))
0272             transform::apply(in, out, input_angles);
0273     }
0274 
0275     template <typename InGeometry>
0276     transform_geometry_wrapper(InGeometry const& in, OutGeometry & out, bool input_angles)
0277         : m_out(out)
0278     {
0279         transform::apply(in, out, input_angles);
0280     }
0281 
0282     OutGeometry & get() { return m_out; }
0283     void finish() {}
0284 
0285 private:
0286     OutGeometry & m_out;
0287 };
0288 
0289 template <typename CT>
0290 struct transform_range
0291 {
0292     template
0293     <
0294         typename Proj1, typename Par1,
0295         typename Proj2, typename Par2,
0296         typename RangeIn, typename RangeOut,
0297         typename Grids
0298     >
0299     static inline bool apply(Proj1 const& proj1, Par1 const& par1,
0300                              Proj2 const& proj2, Par2 const& par2,
0301                              RangeIn const& in, RangeOut & out,
0302                              Grids const& grids1, Grids const& grids2)
0303     {
0304         // NOTE: this has to be consistent with pj_transform()
0305         bool const input_angles = !par1.is_geocent && par1.is_latlong;
0306 
0307         transform_geometry_wrapper<RangeOut, CT> wrapper(in, out, input_angles);
0308 
0309         bool res = true;
0310         try
0311         {
0312             res = pj_transform(proj1, par1, proj2, par2, wrapper.get(), grids1, grids2);
0313         }
0314         catch (projection_exception const&)
0315         {
0316             res = false;
0317         }
0318         catch(...)
0319         {
0320             BOOST_RETHROW
0321         }
0322 
0323         wrapper.finish();
0324 
0325         return res;
0326     }
0327 };
0328 
0329 template <typename Policy>
0330 struct transform_multi
0331 {
0332     template
0333     <
0334         typename Proj1, typename Par1,
0335         typename Proj2, typename Par2,
0336         typename MultiIn, typename MultiOut,
0337         typename Grids
0338     >
0339     static inline bool apply(Proj1 const& proj1, Par1 const& par1,
0340                              Proj2 const& proj2, Par2 const& par2,
0341                              MultiIn const& in, MultiOut & out,
0342                              Grids const& grids1, Grids const& grids2)
0343     {
0344         if (! same_object(in, out))
0345             range::resize(out, boost::size(in));
0346 
0347         return apply(proj1, par1, proj2, par2,
0348                      boost::begin(in), boost::end(in),
0349                      boost::begin(out),
0350                      grids1, grids2);
0351     }
0352 
0353 private:
0354     template
0355     <
0356         typename Proj1, typename Par1,
0357         typename Proj2, typename Par2,
0358         typename InIt, typename OutIt,
0359         typename Grids
0360     >
0361     static inline bool apply(Proj1 const& proj1, Par1 const& par1,
0362                              Proj2 const& proj2, Par2 const& par2,
0363                              InIt in_first, InIt in_last, OutIt out_first,
0364                              Grids const& grids1, Grids const& grids2)
0365     {
0366         bool res = true;
0367         for ( ; in_first != in_last ; ++in_first, ++out_first )
0368         {
0369             if ( ! Policy::apply(proj1, par1, proj2, par2, *in_first, *out_first, grids1, grids2) )
0370             {
0371                 res = false;
0372             }
0373         }
0374         return res;
0375     }
0376 };
0377 
0378 template
0379 <
0380     typename Geometry,
0381     typename CT,
0382     typename Tag = typename geometry::tag<Geometry>::type
0383 >
0384 struct transform
0385     : not_implemented<Tag>
0386 {};
0387 
0388 template <typename Point, typename CT>
0389 struct transform<Point, CT, point_tag>
0390 {
0391     template
0392     <
0393         typename Proj1, typename Par1,
0394         typename Proj2, typename Par2,
0395         typename PointIn, typename PointOut,
0396         typename Grids
0397     >
0398     static inline bool apply(Proj1 const& proj1, Par1 const& par1,
0399                              Proj2 const& proj2, Par2 const& par2,
0400                              PointIn const& in, PointOut & out,
0401                              Grids const& grids1, Grids const& grids2)
0402     {
0403         // NOTE: this has to be consistent with pj_transform()
0404         bool const input_angles = !par1.is_geocent && par1.is_latlong;
0405 
0406         transform_geometry_wrapper<PointOut, CT> wrapper(in, out, input_angles);
0407 
0408         typedef typename transform_geometry_wrapper<PointOut, CT>::type point_type;
0409         point_type * ptr = boost::addressof(wrapper.get());
0410 
0411         std::pair<point_type *, point_type *> range = std::make_pair(ptr, ptr + 1);
0412 
0413         bool res = true;
0414         try
0415         {
0416             res = pj_transform(proj1, par1, proj2, par2, range, grids1, grids2);
0417         }
0418         catch (projection_exception const&)
0419         {
0420             res = false;
0421         }
0422         catch(...)
0423         {
0424             BOOST_RETHROW
0425         }
0426 
0427         wrapper.finish();
0428 
0429         return res;
0430     }
0431 };
0432 
0433 template <typename MultiPoint, typename CT>
0434 struct transform<MultiPoint, CT, multi_point_tag>
0435     : transform_range<CT>
0436 {};
0437 
0438 template <typename Segment, typename CT>
0439 struct transform<Segment, CT, segment_tag>
0440 {
0441     template
0442     <
0443         typename Proj1, typename Par1,
0444         typename Proj2, typename Par2,
0445         typename SegmentIn, typename SegmentOut,
0446         typename Grids
0447     >
0448     static inline bool apply(Proj1 const& proj1, Par1 const& par1,
0449                              Proj2 const& proj2, Par2 const& par2,
0450                              SegmentIn const& in, SegmentOut & out,
0451                              Grids const& grids1, Grids const& grids2)
0452     {
0453         // NOTE: this has to be consistent with pj_transform()
0454         bool const input_angles = !par1.is_geocent && par1.is_latlong;
0455 
0456         transform_geometry_wrapper<SegmentOut, CT> wrapper(in, out, input_angles);
0457 
0458         typedef typename geometry::point_type
0459             <
0460                 typename transform_geometry_wrapper<SegmentOut, CT>::type
0461             >::type point_type;
0462 
0463         point_type points[2];
0464 
0465         geometry::detail::assign_point_from_index<0>(wrapper.get(), points[0]);
0466         geometry::detail::assign_point_from_index<1>(wrapper.get(), points[1]);
0467 
0468         std::pair<point_type*, point_type*> range = std::make_pair(points, points + 2);
0469 
0470         bool res = true;
0471         try
0472         {
0473             res = pj_transform(proj1, par1, proj2, par2, range, grids1, grids2);
0474         }
0475         catch (projection_exception const&)
0476         {
0477             res = false;
0478         }
0479         catch(...)
0480         {
0481             BOOST_RETHROW
0482         }
0483 
0484         geometry::detail::assign_point_to_index<0>(points[0], wrapper.get());
0485         geometry::detail::assign_point_to_index<1>(points[1], wrapper.get());
0486 
0487         wrapper.finish();
0488 
0489         return res;
0490     }
0491 };
0492 
0493 template <typename Linestring, typename CT>
0494 struct transform<Linestring, CT, linestring_tag>
0495     : transform_range<CT>
0496 {};
0497 
0498 template <typename MultiLinestring, typename CT>
0499 struct transform<MultiLinestring, CT, multi_linestring_tag>
0500     : transform_multi<transform_range<CT> >
0501 {};
0502 
0503 template <typename Ring, typename CT>
0504 struct transform<Ring, CT, ring_tag>
0505     : transform_range<CT>
0506 {};
0507 
0508 template <typename Polygon, typename CT>
0509 struct transform<Polygon, CT, polygon_tag>
0510 {
0511     template
0512     <
0513         typename Proj1, typename Par1,
0514         typename Proj2, typename Par2,
0515         typename PolygonIn, typename PolygonOut,
0516         typename Grids
0517     >
0518     static inline bool apply(Proj1 const& proj1, Par1 const& par1,
0519                              Proj2 const& proj2, Par2 const& par2,
0520                              PolygonIn const& in, PolygonOut & out,
0521                              Grids const& grids1, Grids const& grids2)
0522     {
0523         bool r1 = transform_range
0524                     <
0525                         CT
0526                     >::apply(proj1, par1, proj2, par2,
0527                              geometry::exterior_ring(in),
0528                              geometry::exterior_ring(out),
0529                              grids1, grids2);
0530         bool r2 = transform_multi
0531                     <
0532                         transform_range<CT>
0533                      >::apply(proj1, par1, proj2, par2,
0534                               geometry::interior_rings(in),
0535                               geometry::interior_rings(out),
0536                               grids1, grids2);
0537         return r1 && r2;
0538     }
0539 };
0540 
0541 template <typename MultiPolygon, typename CT>
0542 struct transform<MultiPolygon, CT, multi_polygon_tag>
0543     : transform_multi
0544         <
0545             transform
0546                 <
0547                     typename boost::range_value<MultiPolygon>::type,
0548                     CT,
0549                     polygon_tag
0550                 >
0551         >
0552 {};
0553 
0554 
0555 }} // namespace projections::detail
0556 
0557 namespace srs
0558 {
0559 
0560 
0561 /*!
0562     \brief Representation of projection
0563     \details Either dynamic or static projection representation
0564     \ingroup projection
0565     \tparam Proj1 default_dynamic or static projection parameters
0566     \tparam Proj2 default_dynamic or static projection parameters
0567     \tparam CT calculation type used internally
0568 */
0569 template
0570 <
0571     typename Proj1 = srs::dynamic,
0572     typename Proj2 = srs::dynamic,
0573     typename CT = double
0574 >
0575 class transformation
0576 {
0577     typedef typename projections::detail::promote_to_double<CT>::type calc_t;
0578 
0579 public:
0580     // Both static and default constructed
0581     transformation()
0582     {}
0583 
0584     // First dynamic, second static and default constructed
0585     template
0586     <
0587         typename Parameters1,
0588         std::enable_if_t
0589             <
0590                 std::is_same<Proj1, srs::dynamic>::value
0591              && projections::dynamic_parameters<Parameters1>::is_specialized,
0592                 int
0593             > = 0
0594     >
0595     explicit transformation(Parameters1 const& parameters1)
0596         : m_proj1(parameters1)
0597     {}
0598 
0599     // First static, second static and default constructed
0600     explicit transformation(Proj1 const& parameters1)
0601         : m_proj1(parameters1)
0602     {}
0603 
0604     // Both dynamic
0605     template
0606     <
0607         typename Parameters1, typename Parameters2,
0608         std::enable_if_t
0609             <
0610                 std::is_same<Proj1, srs::dynamic>::value
0611              && std::is_same<Proj2, srs::dynamic>::value
0612              && projections::dynamic_parameters<Parameters1>::is_specialized
0613              && projections::dynamic_parameters<Parameters2>::is_specialized,
0614                 int
0615             > = 0
0616     >
0617     transformation(Parameters1 const& parameters1,
0618                    Parameters2 const& parameters2)
0619         : m_proj1(parameters1)
0620         , m_proj2(parameters2)
0621     {}
0622 
0623     // First dynamic, second static
0624     template
0625     <
0626         typename Parameters1,
0627         std::enable_if_t
0628             <
0629                 std::is_same<Proj1, srs::dynamic>::value
0630              && projections::dynamic_parameters<Parameters1>::is_specialized,
0631                 int
0632             > = 0
0633     >
0634     transformation(Parameters1 const& parameters1,
0635                    Proj2 const& parameters2)
0636         : m_proj1(parameters1)
0637         , m_proj2(parameters2)
0638     {}
0639 
0640     // First static, second dynamic
0641     template
0642     <
0643         typename Parameters2,
0644         std::enable_if_t
0645             <
0646                 std::is_same<Proj2, srs::dynamic>::value
0647              && projections::dynamic_parameters<Parameters2>::is_specialized,
0648                 int
0649             > = 0
0650     >
0651     transformation(Proj1 const& parameters1,
0652                    Parameters2 const& parameters2)
0653         : m_proj1(parameters1)
0654         , m_proj2(parameters2)
0655     {}
0656 
0657     // Both static
0658     transformation(Proj1 const& parameters1,
0659                    Proj2 const& parameters2)
0660         : m_proj1(parameters1)
0661         , m_proj2(parameters2)
0662     {}
0663 
0664     template <typename GeometryIn, typename GeometryOut>
0665     bool forward(GeometryIn const& in, GeometryOut & out) const
0666     {
0667         return forward(in, out, transformation_grids<detail::empty_grids_storage>());
0668     }
0669 
0670     template <typename GeometryIn, typename GeometryOut>
0671     bool inverse(GeometryIn const& in, GeometryOut & out) const
0672     {
0673         return inverse(in, out, transformation_grids<detail::empty_grids_storage>());
0674     }
0675 
0676     template <typename GeometryIn, typename GeometryOut, typename GridsStorage>
0677     bool forward(GeometryIn const& in, GeometryOut & out,
0678                  transformation_grids<GridsStorage> const& grids) const
0679     {
0680         BOOST_GEOMETRY_STATIC_ASSERT(
0681             (projections::detail::same_tags<GeometryIn, GeometryOut>::value),
0682             "Not supported combination of Geometries.",
0683             GeometryIn, GeometryOut);
0684 
0685         return projections::detail::transform
0686                 <
0687                     GeometryOut,
0688                     calc_t
0689                 >::apply(m_proj1.proj(), m_proj1.proj().params(),
0690                          m_proj2.proj(), m_proj2.proj().params(),
0691                          in, out,
0692                          grids.src_grids,
0693                          grids.dst_grids);
0694     }
0695 
0696     template <typename GeometryIn, typename GeometryOut, typename GridsStorage>
0697     bool inverse(GeometryIn const& in, GeometryOut & out,
0698                  transformation_grids<GridsStorage> const& grids) const
0699     {
0700         BOOST_GEOMETRY_STATIC_ASSERT(
0701             (projections::detail::same_tags<GeometryIn, GeometryOut>::value),
0702             "Not supported combination of Geometries.",
0703             GeometryIn, GeometryOut);
0704 
0705         return projections::detail::transform
0706                 <
0707                     GeometryOut,
0708                     calc_t
0709                 >::apply(m_proj2.proj(), m_proj2.proj().params(),
0710                          m_proj1.proj(), m_proj1.proj().params(),
0711                          in, out,
0712                          grids.dst_grids,
0713                          grids.src_grids);
0714     }
0715 
0716     template <typename GridsStorage>
0717     inline transformation_grids<GridsStorage> initialize_grids(GridsStorage & grids_storage) const
0718     {
0719         transformation_grids<GridsStorage> result(grids_storage);
0720 
0721         using namespace projections::detail;
0722         pj_gridlist_from_nadgrids(m_proj1.proj().params(),
0723                                   result.src_grids);
0724         pj_gridlist_from_nadgrids(m_proj2.proj().params(),
0725                                   result.dst_grids);
0726 
0727         return result;
0728     }
0729 
0730 private:
0731     projections::proj_wrapper<Proj1, CT> m_proj1;
0732     projections::proj_wrapper<Proj2, CT> m_proj2;
0733 };
0734 
0735 
0736 } // namespace srs
0737 
0738 
0739 }} // namespace boost::geometry
0740 
0741 
0742 #endif // BOOST_GEOMETRY_SRS_TRANSFORMATION_HPP