]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Boost.Geometry (aka GGL, Generic Geometry Library) |
2 | // Unit Test | |
3 | ||
4 | // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. | |
5 | // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. | |
6 | // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. | |
7 | ||
8 | // This file was modified by Oracle on 2014, 2015. | |
9 | // Modifications copyright (c) 2014-2015 Oracle and/or its affiliates. | |
10 | ||
11 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
12 | ||
13 | // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library | |
14 | // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. | |
15 | ||
16 | // Use, modification and distribution is subject to the Boost Software License, | |
17 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
18 | // http://www.boost.org/LICENSE_1_0.txt) | |
19 | ||
20 | ||
21 | #include <geometry_test_common.hpp> | |
22 | ||
23 | #include <boost/concept_check.hpp> | |
24 | ||
25 | #include <boost/geometry/strategies/geographic/distance_vincenty.hpp> | |
26 | #include <boost/geometry/strategies/geographic/side_vincenty.hpp> | |
27 | #include <boost/geometry/algorithms/detail/vincenty_inverse.hpp> | |
28 | #include <boost/geometry/algorithms/detail/vincenty_direct.hpp> | |
29 | ||
30 | #include <boost/geometry/core/srs.hpp> | |
31 | #include <boost/geometry/strategies/strategies.hpp> | |
32 | #include <boost/geometry/algorithms/assign.hpp> | |
33 | #include <boost/geometry/geometries/point.hpp> | |
34 | #include <test_common/test_point.hpp> | |
35 | ||
36 | #ifdef HAVE_TTMATH | |
37 | # include <boost/geometry/extensions/contrib/ttmath_stub.hpp> | |
38 | #endif | |
39 | ||
40 | template <typename T> | |
41 | void normalize_deg(T & deg) | |
42 | { | |
43 | while ( deg > T(180) ) | |
44 | deg -= T(360); | |
45 | while ( deg <= T(-180) ) | |
46 | deg += T(360); | |
47 | } | |
48 | ||
49 | template <typename T> | |
50 | T difference_deg(T const& a1, T const& a2) | |
51 | { | |
52 | T d = a1 - a2; | |
53 | normalize_deg(d); | |
54 | return d; | |
55 | } | |
56 | ||
57 | template <typename T> | |
58 | void check_deg(std::string const& name, T const& a1, T const& a2, T const& percent, T const& error) | |
59 | { | |
60 | T diff = bg::math::abs(difference_deg(a1, a2)); | |
61 | ||
62 | if ( bg::math::equals(a1, T(0)) || bg::math::equals(a2, T(0)) ) | |
63 | { | |
64 | if ( diff > error ) | |
65 | { | |
66 | BOOST_ERROR(name << " - the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << error << "}"); | |
67 | } | |
68 | } | |
69 | else | |
70 | { | |
71 | T greater = (std::max)(bg::math::abs(a1), bg::math::abs(a2)); | |
72 | ||
73 | if ( diff > greater * percent / T(100) ) | |
74 | { | |
75 | BOOST_ERROR(name << " the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << percent << "}%"); | |
76 | } | |
77 | } | |
78 | } | |
79 | ||
80 | double azimuth(double deg, double min, double sec) | |
81 | { | |
82 | min = fabs(min); | |
83 | sec = fabs(sec); | |
84 | ||
85 | if ( deg < 0 ) | |
86 | { | |
87 | min = -min; | |
88 | sec = -sec; | |
89 | } | |
90 | ||
91 | return deg + min/60.0 + sec/3600.0; | |
92 | } | |
93 | ||
94 | double azimuth(double deg, double min) | |
95 | { | |
96 | return azimuth(deg, min, 0.0); | |
97 | } | |
98 | ||
99 | template <typename P> | |
100 | bool non_precise_ct() | |
101 | { | |
102 | typedef typename bg::coordinate_type<P>::type ct; | |
103 | return boost::is_integral<ct>::value || boost::is_float<ct>::value; | |
104 | } | |
105 | ||
106 | template <typename P1, typename P2, typename Spheroid> | |
107 | void test_vincenty(double lon1, double lat1, double lon2, double lat2, | |
108 | double expected_distance, | |
109 | double expected_azimuth_12, | |
110 | double /*expected_azimuth_21*/, | |
111 | Spheroid const& spheroid) | |
112 | { | |
113 | typedef typename bg::promote_floating_point | |
114 | < | |
115 | typename bg::select_calculation_type<P1, P2, void>::type | |
116 | >::type calc_t; | |
117 | ||
118 | calc_t tolerance = non_precise_ct<P1>() || non_precise_ct<P2>() ? | |
119 | 5.0 : 0.001; | |
120 | calc_t error = non_precise_ct<P1>() || non_precise_ct<P2>() ? | |
121 | 1e-5 : 1e-12; | |
122 | ||
123 | // formula | |
124 | { | |
125 | double const d2r = bg::math::d2r<double>(); | |
126 | double const r2d = bg::math::r2d<double>(); | |
127 | ||
128 | typedef bg::detail::vincenty_inverse<calc_t, true, true> inverse_formula; | |
129 | typename inverse_formula::result_type | |
130 | result_i = inverse_formula::apply(lon1 * d2r, | |
131 | lat1 * d2r, | |
132 | lon2 * d2r, | |
133 | lat2 * d2r, | |
134 | spheroid); | |
135 | calc_t dist = result_i.distance; | |
136 | calc_t az12 = result_i.azimuth; | |
137 | //calc_t az21 = vi.azimuth21(); | |
138 | ||
139 | calc_t az12_deg = az12 * r2d; | |
140 | //calc_t az21_deg = az21 * r2d; | |
141 | ||
142 | BOOST_CHECK_CLOSE(dist, calc_t(expected_distance), tolerance); | |
143 | check_deg("az12_deg", az12_deg, calc_t(expected_azimuth_12), tolerance, error); | |
144 | //check_deg("az21_deg", az21_deg, calc_t(expected_azimuth_21), tolerance, error); | |
145 | ||
146 | typedef bg::detail::vincenty_direct<calc_t> direct_formula; | |
147 | typename direct_formula::result_type | |
148 | result_d = direct_formula::apply(lon1 * d2r, | |
149 | lat1 * d2r, | |
150 | dist, | |
151 | az12, | |
152 | spheroid); | |
153 | calc_t direct_lon2 = result_d.lon2; | |
154 | calc_t direct_lat2 = result_d.lat2; | |
155 | //calc_t direct_az21 = vd.azimuth21(); | |
156 | ||
157 | calc_t direct_lon2_deg = direct_lon2 * r2d; | |
158 | calc_t direct_lat2_deg = direct_lat2 * r2d; | |
159 | //calc_t direct_az21_deg = direct_az21 * r2d; | |
160 | ||
161 | check_deg("direct_lon2_deg", direct_lon2_deg, calc_t(lon2), tolerance, error); | |
162 | check_deg("direct_lat2_deg", direct_lat2_deg, calc_t(lat2), tolerance, error); | |
163 | //check_deg("direct_az21_deg", direct_az21_deg, az21_deg, tolerance, error); | |
164 | } | |
165 | ||
166 | // distance strategy | |
167 | { | |
168 | typedef bg::strategy::distance::vincenty<Spheroid> vincenty_type; | |
169 | ||
170 | BOOST_CONCEPT_ASSERT( | |
171 | ( | |
172 | bg::concepts::PointDistanceStrategy<vincenty_type, P1, P2>) | |
173 | ); | |
174 | ||
175 | vincenty_type vincenty(spheroid); | |
176 | typedef typename bg::strategy::distance::services::return_type<vincenty_type, P1, P2>::type return_type; | |
177 | ||
178 | P1 p1; | |
179 | P2 p2; | |
180 | ||
181 | bg::assign_values(p1, lon1, lat1); | |
182 | bg::assign_values(p2, lon2, lat2); | |
183 | ||
184 | BOOST_CHECK_CLOSE(vincenty.apply(p1, p2), return_type(expected_distance), tolerance); | |
185 | BOOST_CHECK_CLOSE(bg::distance(p1, p2, vincenty), return_type(expected_distance), tolerance); | |
186 | } | |
187 | } | |
188 | ||
189 | template <typename P1, typename P2> | |
190 | void test_vincenty(double lon1, double lat1, double lon2, double lat2, | |
191 | double expected_distance, | |
192 | double expected_azimuth_12, | |
193 | double expected_azimuth_21) | |
194 | { | |
195 | test_vincenty<P1, P2>(lon1, lat1, lon2, lat2, | |
196 | expected_distance, expected_azimuth_12, expected_azimuth_21, | |
197 | bg::srs::spheroid<double>()); | |
198 | } | |
199 | ||
200 | template <typename PS, typename P> | |
201 | void test_side(double lon1, double lat1, | |
202 | double lon2, double lat2, | |
203 | double lon, double lat, | |
204 | int expected_side) | |
205 | { | |
206 | // Set radius type, but for integer coordinates we want to have floating point radius type | |
207 | typedef typename bg::promote_floating_point | |
208 | < | |
209 | typename bg::coordinate_type<PS>::type | |
210 | >::type rtype; | |
211 | ||
212 | typedef bg::srs::spheroid<rtype> stype; | |
213 | ||
214 | typedef bg::strategy::side::vincenty<stype> strategy_type; | |
215 | ||
216 | strategy_type strategy; | |
217 | ||
218 | PS p1, p2; | |
219 | P p; | |
220 | ||
221 | bg::assign_values(p1, lon1, lat1); | |
222 | bg::assign_values(p2, lon2, lat2); | |
223 | bg::assign_values(p, lon, lat); | |
224 | ||
225 | int side = strategy.apply(p1, p2, p); | |
226 | ||
227 | BOOST_CHECK_EQUAL(side, expected_side); | |
228 | } | |
229 | ||
230 | template <typename P1, typename P2> | |
231 | void test_all() | |
232 | { | |
233 | // See: | |
234 | // - http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp | |
235 | // - http://www.ga.gov.au/geodesy/datums/vincenty_direct.jsp | |
236 | // Values in the comments below was calculated using the above pages | |
237 | // in some cases distances may be different, previously used values was left | |
238 | ||
239 | // use km | |
240 | double gda_a = 6378.1370; | |
241 | double gda_f = 1.0 / 298.25722210; | |
242 | double gda_b = gda_a * ( 1.0 - gda_f ); | |
243 | bg::srs::spheroid<double> gda_spheroid(gda_a, gda_b); | |
244 | ||
245 | // Test fractional coordinates only for non-integral types | |
246 | if ( BOOST_GEOMETRY_CONDITION( | |
247 | ! boost::is_integral<typename bg::coordinate_type<P1>::type>::value | |
248 | && ! boost::is_integral<typename bg::coordinate_type<P2>::type>::value ) ) | |
249 | { | |
250 | // Flinders Peak -> Buninyong | |
251 | test_vincenty<P1, P2>(azimuth(144,25,29.52440), azimuth(-37,57,3.72030), | |
252 | azimuth(143,55,35.38390), azimuth(-37,39,10.15610), | |
253 | 54.972271, azimuth(306,52,5.37), azimuth(127,10,25.07), | |
254 | gda_spheroid); | |
255 | } | |
256 | ||
257 | // Lodz -> Trondheim | |
258 | test_vincenty<P1, P2>(azimuth(19,28), azimuth(51,47), | |
259 | azimuth(10,21), azimuth(63,23), | |
260 | 1399.032724, azimuth(340,54,25.14), azimuth(153,10,0.19), | |
261 | gda_spheroid); | |
262 | // London -> New York | |
263 | test_vincenty<P1, P2>(azimuth(0,7,39), azimuth(51,30,26), | |
264 | azimuth(-74,0,21), azimuth(40,42,46), | |
265 | 5602.044851, azimuth(288,31,36.82), azimuth(51,10,33.43), | |
266 | gda_spheroid); | |
267 | ||
268 | // Shanghai -> San Francisco | |
269 | test_vincenty<P1, P2>(azimuth(121,30), azimuth(31,12), | |
270 | azimuth(-122,25), azimuth(37,47), | |
271 | 9899.698550, azimuth(45,12,44.76), azimuth(309,50,20.88), | |
272 | gda_spheroid); | |
273 | ||
274 | test_vincenty<P1, P2>(0, 0, 0, 50, 5540.847042, 0, 180, gda_spheroid); // N | |
275 | test_vincenty<P1, P2>(0, 0, 0, -50, 5540.847042, 180, 0, gda_spheroid); // S | |
276 | test_vincenty<P1, P2>(0, 0, 50, 0, 5565.974540, 90, -90, gda_spheroid); // E | |
277 | test_vincenty<P1, P2>(0, 0, -50, 0, 5565.974540, -90, 90, gda_spheroid); // W | |
278 | ||
279 | test_vincenty<P1, P2>(0, 0, 50, 50, 7284.879297, azimuth(32,51,55.87), azimuth(237,24,50.12), gda_spheroid); // NE | |
280 | ||
281 | // The original distance values, azimuths calculated using the web form mentioned above | |
282 | // Using default spheroid units (meters) | |
283 | test_vincenty<P1, P2>(0, 89, 1, 80, 1005153.5769, azimuth(178,53,23.85), azimuth(359,53,18.35)); // sub-polar | |
284 | test_vincenty<P1, P2>(4, 52, 4, 52, 0.0, 0, 0); // no point difference | |
285 | test_vincenty<P1, P2>(4, 52, 3, 40, 1336039.890, azimuth(183,41,29.08), azimuth(2,58,5.13)); // normal case | |
286 | ||
287 | test_side<P1, P2>(0, 0, 0, 1, 0, 2, 0); | |
288 | test_side<P1, P2>(0, 0, 0, 1, 0, -2, 0); | |
289 | test_side<P1, P2>(10, 0, 10, 1, 10, 2, 0); | |
290 | test_side<P1, P2>(10, 0, 10, -1, 10, 2, 0); | |
291 | ||
292 | test_side<P1, P2>(10, 0, 10, 1, 0, 2, 1); // left | |
293 | test_side<P1, P2>(10, 0, 10, -1, 0, 2, -1); // right | |
294 | ||
295 | test_side<P1, P2>(-10, -10, 10, 10, 10, 0, -1); // right | |
296 | test_side<P1, P2>(-10, -10, 10, 10, -10, 0, 1); // left | |
297 | test_side<P1, P2>(170, -10, -170, 10, -170, 0, -1); // right | |
298 | test_side<P1, P2>(170, -10, -170, 10, 170, 0, 1); // left | |
299 | } | |
300 | ||
301 | template <typename P> | |
302 | void test_all() | |
303 | { | |
304 | test_all<P, P>(); | |
305 | } | |
306 | ||
307 | int test_main(int, char* []) | |
308 | { | |
309 | //test_all<float[2]>(); | |
310 | //test_all<double[2]>(); | |
311 | test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >(); | |
312 | test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >(); | |
313 | test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >(); | |
314 | ||
315 | #if defined(HAVE_TTMATH) | |
316 | test_all<bg::model::point<ttmath::Big<1,4>, 2, bg::cs::geographic<bg::degree> > >(); | |
317 | test_all<bg::model::point<ttmath_big, 2, bg::cs::geographic<bg::degree> > >(); | |
318 | #endif | |
319 | ||
320 | ||
321 | return 0; | |
322 | } |