1 // Boost.Geometry (aka GGL, Generic Geometry Library) 2 3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. 4 5 // Use, modification and distribution is subject to the Boost Software License, 6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 7 // http://www.boost.org/LICENSE_1_0.txt) 8 9 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP 10 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP 11 12 13 #include <boost/geometry/core/cs.hpp> 14 #include <boost/geometry/core/access.hpp> 15 #include <boost/geometry/core/radian_access.hpp> 16 17 #include <boost/geometry/util/math.hpp> 18 #include <boost/geometry/util/select_calculation_type.hpp> 19 #include <boost/geometry/util/promote_floating_point.hpp> 20 21 #include <boost/geometry/strategies/distance.hpp> 22 23 24 25 namespace boost { namespace geometry 26 { 27 28 29 namespace strategy { namespace distance 30 { 31 32 33 namespace comparable 34 { 35 36 // Comparable haversine. 37 // To compare distances, we can avoid: 38 // - multiplication with radius and 2.0 39 // - applying sqrt 40 // - applying asin (which is strictly (monotone) increasing) 41 template 42 < 43 typename RadiusType, 44 typename CalculationType = void 45 > 46 class haversine 47 { 48 public : 49 template <typename Point1, typename Point2> 50 struct calculation_type 51 : promote_floating_point 52 < 53 typename select_calculation_type 54 < 55 Point1, 56 Point2, 57 CalculationType 58 >::type 59 > 60 {}; 61 62 typedef RadiusType radius_type; 63 haversine(RadiusType const & r=1.0)64 explicit inline haversine(RadiusType const& r = 1.0) 65 : m_radius(r) 66 {} 67 68 template <typename Point1, typename Point2> 69 static inline typename calculation_type<Point1, Point2>::type apply(Point1 const & p1,Point2 const & p2)70 apply(Point1 const& p1, Point2 const& p2) 71 { 72 return calculate<typename calculation_type<Point1, Point2>::type>( 73 get_as_radian<0>(p1), get_as_radian<1>(p1), 74 get_as_radian<0>(p2), get_as_radian<1>(p2) 75 ); 76 } 77 radius() const78 inline RadiusType radius() const 79 { 80 return m_radius; 81 } 82 83 84 private : 85 template <typename R, typename T1, typename T2> calculate(T1 const & lon1,T1 const & lat1,T2 const & lon2,T2 const & lat2)86 static inline R calculate(T1 const& lon1, T1 const& lat1, 87 T2 const& lon2, T2 const& lat2) 88 { 89 return math::hav(lat2 - lat1) 90 + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1); 91 } 92 93 RadiusType m_radius; 94 }; 95 96 97 98 } // namespace comparable 99 100 /*! 101 \brief Distance calculation for spherical coordinates 102 on a perfect sphere using haversine 103 \ingroup strategies 104 \tparam RadiusType \tparam_radius 105 \tparam CalculationType \tparam_calculation 106 \author Adapted from: http://williams.best.vwh.net/avform.htm 107 \see http://en.wikipedia.org/wiki/Great-circle_distance 108 \note (from Wiki:) The great circle distance d between two 109 points with coordinates {lat1,lon1} and {lat2,lon2} is given by: 110 d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2)) 111 A mathematically equivalent formula, which is less subject 112 to rounding error for short distances is: 113 d=2*asin(sqrt((sin((lat1-lat2) / 2))^2 114 + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2)) 115 116 117 \qbk{ 118 [heading See also] 119 [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] 120 } 121 122 */ 123 template 124 < 125 typename RadiusType, 126 typename CalculationType = void 127 > 128 class haversine 129 { 130 typedef comparable::haversine<RadiusType, CalculationType> comparable_type; 131 132 public : 133 template <typename Point1, typename Point2> 134 struct calculation_type 135 : services::return_type<comparable_type, Point1, Point2> 136 {}; 137 138 typedef RadiusType radius_type; 139 140 /*! 141 \brief Constructor 142 \param radius radius of the sphere, defaults to 1.0 for the unit sphere 143 */ haversine(RadiusType const & radius=1.0)144 inline haversine(RadiusType const& radius = 1.0) 145 : m_radius(radius) 146 {} 147 148 /*! 149 \brief applies the distance calculation 150 \return the calculated distance (including multiplying with radius) 151 \param p1 first point 152 \param p2 second point 153 */ 154 template <typename Point1, typename Point2> 155 inline typename calculation_type<Point1, Point2>::type apply(Point1 const & p1,Point2 const & p2) const156 apply(Point1 const& p1, Point2 const& p2) const 157 { 158 typedef typename calculation_type<Point1, Point2>::type calculation_type; 159 calculation_type const a = comparable_type::apply(p1, p2); 160 calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a)); 161 return calculation_type(m_radius) * c; 162 } 163 164 /*! 165 \brief access to radius value 166 \return the radius 167 */ radius() const168 inline RadiusType radius() const 169 { 170 return m_radius; 171 } 172 173 private : 174 RadiusType m_radius; 175 }; 176 177 178 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 179 namespace services 180 { 181 182 template <typename RadiusType, typename CalculationType> 183 struct tag<haversine<RadiusType, CalculationType> > 184 { 185 typedef strategy_tag_distance_point_point type; 186 }; 187 188 189 template <typename RadiusType, typename CalculationType, typename P1, typename P2> 190 struct return_type<haversine<RadiusType, CalculationType>, P1, P2> 191 : haversine<RadiusType, CalculationType>::template calculation_type<P1, P2> 192 {}; 193 194 195 template <typename RadiusType, typename CalculationType> 196 struct comparable_type<haversine<RadiusType, CalculationType> > 197 { 198 typedef comparable::haversine<RadiusType, CalculationType> type; 199 }; 200 201 202 template <typename RadiusType, typename CalculationType> 203 struct get_comparable<haversine<RadiusType, CalculationType> > 204 { 205 private : 206 typedef haversine<RadiusType, CalculationType> this_type; 207 typedef comparable::haversine<RadiusType, CalculationType> comparable_type; 208 public : applyboost::geometry::strategy::distance::services::get_comparable209 static inline comparable_type apply(this_type const& input) 210 { 211 return comparable_type(input.radius()); 212 } 213 }; 214 215 template <typename RadiusType, typename CalculationType, typename P1, typename P2> 216 struct result_from_distance<haversine<RadiusType, CalculationType>, P1, P2> 217 { 218 private : 219 typedef haversine<RadiusType, CalculationType> this_type; 220 typedef typename return_type<this_type, P1, P2>::type return_type; 221 public : 222 template <typename T> applyboost::geometry::strategy::distance::services::result_from_distance223 static inline return_type apply(this_type const& , T const& value) 224 { 225 return return_type(value); 226 } 227 }; 228 229 230 // Specializations for comparable::haversine 231 template <typename RadiusType, typename CalculationType> 232 struct tag<comparable::haversine<RadiusType, CalculationType> > 233 { 234 typedef strategy_tag_distance_point_point type; 235 }; 236 237 238 template <typename RadiusType, typename CalculationType, typename P1, typename P2> 239 struct return_type<comparable::haversine<RadiusType, CalculationType>, P1, P2> 240 : comparable::haversine<RadiusType, CalculationType>::template calculation_type<P1, P2> 241 {}; 242 243 244 template <typename RadiusType, typename CalculationType> 245 struct comparable_type<comparable::haversine<RadiusType, CalculationType> > 246 { 247 typedef comparable::haversine<RadiusType, CalculationType> type; 248 }; 249 250 251 template <typename RadiusType, typename CalculationType> 252 struct get_comparable<comparable::haversine<RadiusType, CalculationType> > 253 { 254 private : 255 typedef comparable::haversine<RadiusType, CalculationType> this_type; 256 public : applyboost::geometry::strategy::distance::services::get_comparable257 static inline this_type apply(this_type const& input) 258 { 259 return input; 260 } 261 }; 262 263 264 template <typename RadiusType, typename CalculationType, typename P1, typename P2> 265 struct result_from_distance<comparable::haversine<RadiusType, CalculationType>, P1, P2> 266 { 267 private : 268 typedef comparable::haversine<RadiusType, CalculationType> strategy_type; 269 typedef typename return_type<strategy_type, P1, P2>::type return_type; 270 public : 271 template <typename T> applyboost::geometry::strategy::distance::services::result_from_distance272 static inline return_type apply(strategy_type const& strategy, T const& distance) 273 { 274 return_type const s = sin((distance / strategy.radius()) / return_type(2)); 275 return s * s; 276 } 277 }; 278 279 280 // Register it as the default for point-types 281 // in a spherical equatorial coordinate system 282 template <typename Point1, typename Point2> 283 struct default_strategy 284 < 285 point_tag, point_tag, Point1, Point2, 286 spherical_equatorial_tag, spherical_equatorial_tag 287 > 288 { 289 typedef strategy::distance::haversine<typename select_coordinate_type<Point1, Point2>::type> type; 290 }; 291 292 // Note: spherical polar coordinate system requires "get_as_radian_equatorial" 293 294 295 } // namespace services 296 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 297 298 299 }} // namespace strategy::distance 300 301 302 }} // namespace boost::geometry 303 304 305 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP 306