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