strategies/spherical/distance_haversine.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018.
0006 // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
0007
0008 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Use, modification and distribution is subject to the Boost Software License,
0012 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0013 // http://www.boost.org/LICENSE_1_0.txt)
0014
0015 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
0016 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
0017
0018
0019 #include <boost/geometry/core/access.hpp>
0020 #include <boost/geometry/core/coordinate_promotion.hpp>
0021 #include <boost/geometry/core/cs.hpp>
0022 #include <boost/geometry/core/radian_access.hpp>
0023
0024 #include <boost/geometry/srs/sphere.hpp>
0025
0026 #include <boost/geometry/strategies/distance.hpp>
0027 #include <boost/geometry/strategies/spherical/get_radius.hpp>
0028
0029 #include <boost/geometry/util/math.hpp>
0030 #include <boost/geometry/util/select_calculation_type.hpp>
0031
0032
0033 namespace boost { namespace geometry
0034 {
0035
0036
0037 namespace strategy { namespace distance
0038 {
0039
0040
0041 namespace comparable
0042 {
0043
0044 // Comparable haversine.
0045 // To compare distances, we can avoid:
0046 // - multiplication with radius and 2.0
0047 // - applying sqrt
0048 // - applying asin (which is strictly (monotone) increasing)
0049 template
0050 <
0051     typename RadiusTypeOrSphere = double,
0052     typename CalculationType = void
0053 >
0054 class haversine
0055 {
0056 public :
0057     template <typename Point1, typename Point2>
0058     struct calculation_type
0059         : promote_floating_point
0060           <
0061               typename select_calculation_type
0062                   <
0063                       Point1,
0064                       Point2,
0065                       CalculationType
0066                   >::type
0067           >
0068     {};
0069
0070     typedef typename strategy_detail::get_radius
0071         <
0072             RadiusTypeOrSphere
0073         >::type radius_type;
0074
0075     inline haversine()
0076         : m_radius(1.0)
0077     {}
0078
0079     template <typename RadiusOrSphere>
0080     explicit inline haversine(RadiusOrSphere const& radius_or_sphere)
0081         : m_radius(strategy_detail::get_radius
0082                     <
0083                         RadiusOrSphere
0084                     >::apply(radius_or_sphere))
0085     {}
0086
0087     template <typename Point1, typename Point2>
0088     static inline typename calculation_type<Point1, Point2>::type
0089     apply(Point1 const& p1, Point2 const& p2)
0090     {
0091         return calculate<typename calculation_type<Point1, Point2>::type>(
0092                    get_as_radian<0>(p1), get_as_radian<1>(p1),
0093                    get_as_radian<0>(p2), get_as_radian<1>(p2)
0094                );
0095     }
0096
0097     inline radius_type radius() const
0098     {
0099         return m_radius;
0100     }
0101
0102
0103 private :
0104     template <typename R, typename T1, typename T2>
0105     static inline R calculate(T1 const& lon1, T1 const& lat1,
0106                               T2 const& lon2, T2 const& lat2)
0107     {
0108         return math::hav(lat2 - lat1)
0109                 + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1);
0110     }
0111
0112     radius_type m_radius;
0113 };
0114
0115
0116
0117 } // namespace comparable
0118
0119 /*!
0120 \brief Distance calculation for spherical coordinates
0121 on a perfect sphere using haversine
0122 \ingroup strategies
0123 \tparam RadiusTypeOrSphere \tparam_radius_or_sphere
0124 \tparam CalculationType \tparam_calculation
0125 \author Adapted from: http://williams.best.vwh.net/avform.htm
0126 \see http://en.wikipedia.org/wiki/Great-circle_distance
0127 \note (from Wiki:) The great circle distance d between two
0128 points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
0129     d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))
0130 A mathematically equivalent formula, which is less subject
0131     to rounding error for short distances is:
0132     d=2*asin(sqrt((sin((lat1-lat2) / 2))^2
0133     + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2))
0134 \qbk{
0135 [heading See also]
0136 [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
0137 }
0138 */
0139 template
0140 <
0141     typename RadiusTypeOrSphere = double,
0142     typename CalculationType = void
0143 >
0144 class haversine
0145 {
0146     typedef comparable::haversine<RadiusTypeOrSphere, CalculationType> comparable_type;
0147
0148 public :
0149     template <typename Point1, typename Point2>
0150     struct calculation_type
0151         : services::return_type<comparable_type, Point1, Point2>
0152     {};
0153
0154     typedef typename strategy_detail::get_radius
0155         <
0156             RadiusTypeOrSphere
0157         >::type radius_type;
0158
0159     /*!
0160     \brief Default constructor, radius set to 1.0 for the unit sphere
0161     */
0162     inline haversine()
0163         : m_radius(1.0)
0164     {}
0165
0166     /*!
0167     \brief Constructor
0168     \param radius_or_sphere radius of the sphere or sphere model
0169     */
0170     template <typename RadiusOrSphere>
0171     explicit inline haversine(RadiusOrSphere const& radius_or_sphere)
0172         : m_radius(strategy_detail::get_radius
0173                     <
0174                         RadiusOrSphere
0175                     >::apply(radius_or_sphere))
0176     {}
0177
0178     /*!
0179     \brief applies the distance calculation
0180     \return the calculated distance (including multiplying with radius)
0181     \param p1 first point
0182     \param p2 second point
0183     */
0184     template <typename Point1, typename Point2>
0185     inline typename calculation_type<Point1, Point2>::type
0186     apply(Point1 const& p1, Point2 const& p2) const
0187     {
0188         typedef typename calculation_type<Point1, Point2>::type calculation_type;
0189         calculation_type const a = comparable_type::apply(p1, p2);
0190         calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a));
0191         return calculation_type(m_radius) * c;
0192     }
0193
0194     /*!
0195     \brief access to radius value
0196     \return the radius
0197     */
0198     inline radius_type radius() const
0199     {
0200         return m_radius;
0201     }
0202
0203 private :
0204     radius_type m_radius;
0205 };
0206
0207
0208 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0209 namespace services
0210 {
0211
0212 template <typename RadiusType, typename CalculationType>
0213 struct tag<haversine<RadiusType, CalculationType> >
0214 {
0215     typedef strategy_tag_distance_point_point type;
0216 };
0217
0218
0219 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
0220 struct return_type<haversine<RadiusType, CalculationType>, P1, P2>
0221     : haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
0222 {};
0223
0224
0225 template <typename RadiusType, typename CalculationType>
0226 struct comparable_type<haversine<RadiusType, CalculationType> >
0227 {
0228     typedef comparable::haversine<RadiusType, CalculationType> type;
0229 };
0230
0231
0232 template <typename RadiusType, typename CalculationType>
0233 struct get_comparable<haversine<RadiusType, CalculationType> >
0234 {
0235 private :
0236     typedef haversine<RadiusType, CalculationType> this_type;
0237     typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
0238 public :
0239     static inline comparable_type apply(this_type const& input)
0240     {
0241         return comparable_type(input.radius());
0242     }
0243 };
0244
0245 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
0246 struct result_from_distance<haversine<RadiusType, CalculationType>, P1, P2>
0247 {
0248 private :
0249     typedef haversine<RadiusType, CalculationType> this_type;
0250     typedef typename return_type<this_type, P1, P2>::type return_type;
0251 public :
0252     template <typename T>
0253     static inline return_type apply(this_type const& , T const& value)
0254     {
0255         return return_type(value);
0256     }
0257 };
0258
0259
0260 // Specializations for comparable::haversine
0261 template <typename RadiusType, typename CalculationType>
0262 struct tag<comparable::haversine<RadiusType, CalculationType> >
0263 {
0264     typedef strategy_tag_distance_point_point type;
0265 };
0266
0267
0268 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
0269 struct return_type<comparable::haversine<RadiusType, CalculationType>, P1, P2>
0270     : comparable::haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
0271 {};
0272
0273
0274 template <typename RadiusType, typename CalculationType>
0275 struct comparable_type<comparable::haversine<RadiusType, CalculationType> >
0276 {
0277     typedef comparable::haversine<RadiusType, CalculationType> type;
0278 };
0279
0280
0281 template <typename RadiusType, typename CalculationType>
0282 struct get_comparable<comparable::haversine<RadiusType, CalculationType> >
0283 {
0284 private :
0285     typedef comparable::haversine<RadiusType, CalculationType> this_type;
0286 public :
0287     static inline this_type apply(this_type const& input)
0288     {
0289         return input;
0290     }
0291 };
0292
0293
0294 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
0295 struct result_from_distance<comparable::haversine<RadiusType, CalculationType>, P1, P2>
0296 {
0297 private :
0298     typedef comparable::haversine<RadiusType, CalculationType> strategy_type;
0299     typedef typename return_type<strategy_type, P1, P2>::type return_type;
0300 public :
0301     template <typename T>
0302     static inline return_type apply(strategy_type const& strategy, T const& distance)
0303     {
0304         return_type const s = sin((distance / strategy.radius()) / return_type(2));
0305         return s * s;
0306     }
0307 };
0308
0309
0310 // Register it as the default for point-types
0311 // in a spherical equatorial coordinate system
0312 template <typename Point1, typename Point2>
0313 struct default_strategy
0314     <
0315         point_tag, point_tag, Point1, Point2,
0316         spherical_equatorial_tag, spherical_equatorial_tag
0317     >
0318 {
0319     typedef strategy::distance::haversine<typename select_coordinate_type<Point1, Point2>::type> type;
0320 };
0321
0322 // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
0323
0324
0325 } // namespace services
0326 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0327
0328
0329 }} // namespace strategy::distance
0330
0331
0332 }} // namespace boost::geometry
0333
0334
0335 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP