1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 // Unit Test
3
4 // Copyright (c) 2007-2016 Barend Gehrels, Amsterdam, the Netherlands.
5 // Copyright (c) 2008-2016 Bruno Lalande, Paris, France.
6 // Copyright (c) 2009-2016 Mateusz Loskot, London, UK.
7
8 // This file was modified by Oracle on 2014-2017.
9 // Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
10
11 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
12 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
13
14 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
15 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
16
17 // Use, modification and distribution is subject to the Boost Software License,
18 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
19 // http://www.boost.org/LICENSE_1_0.txt)
20
21
22 #include <geometry_test_common.hpp>
23
24 #include <boost/concept_check.hpp>
25
26 #include <boost/geometry/algorithms/assign.hpp>
27 #include <boost/geometry/algorithms/distance.hpp>
28 #include <boost/geometry/geometries/point.hpp>
29 #include <boost/geometry/srs/spheroid.hpp>
30 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
31 #include <boost/geometry/strategies/geographic/distance_andoyer.hpp>
32 #include <boost/geometry/strategies/geographic/side_andoyer.hpp>
33
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
41
make_deg(double deg,double min,double sec)42 double make_deg(double deg, double min, double sec)
43 {
44 return deg + min / 60.0 + sec / 3600.0;
45 }
46
to_rad(double deg)47 double to_rad(double deg)
48 {
49 return bg::math::pi<double>() * deg / 180.0;
50 }
51
to_deg(double rad)52 double to_deg(double rad)
53 {
54 return 180.0 * rad / bg::math::pi<double>();
55 }
56
normlized_deg(double deg)57 double normlized_deg(double deg)
58 {
59 if (deg > 180)
60 return deg - 360;
61 else if (deg < -180)
62 return deg + 360;
63 else
64 return deg;
65 }
66
67
68 template <typename P1, typename P2>
test_distance(double lon1,double lat1,double lon2,double lat2,double expected_km)69 void test_distance(double lon1, double lat1, double lon2, double lat2, double expected_km)
70 {
71 // Set radius type, but for integer coordinates we want to have floating point radius type
72 typedef typename bg::promote_floating_point
73 <
74 typename bg::coordinate_type<P1>::type
75 >::type rtype;
76
77 typedef bg::srs::spheroid<rtype> stype;
78
79 typedef bg::strategy::distance::andoyer<stype> andoyer_type;
80 typedef bg::strategy::distance::geographic<bg::strategy::andoyer, stype> geographic_type;
81 typedef bg::formula::andoyer_inverse<rtype, true, false> andoyer_inverse_type;
82
83 BOOST_CONCEPT_ASSERT
84 (
85 (bg::concepts::PointDistanceStrategy<andoyer_type, P1, P2>)
86 );
87
88 andoyer_type andoyer;
89 geographic_type geographic;
90 typedef typename bg::strategy::distance
91 ::services::return_type<andoyer_type, P1, P2>::type return_type;
92
93 P1 p1;
94 P2 p2;
95
96 bg::assign_values(p1, lon1, lat1);
97 bg::assign_values(p2, lon2, lat2);
98
99 return_type d_strategy = andoyer.apply(p1, p2);
100 return_type d_strategy2 = geographic.apply(p1, p2);
101 return_type d_function = bg::distance(p1, p2, andoyer);
102
103 double diff = bg::math::longitude_distance_signed<bg::degree>(lon1, lon2);
104 return_type d_formula;
105
106 // if the points lay on a meridian, distance strategy calls the special formula
107 // for meridian distance that returns different result than andoyer formula
108 // for nearly antipodal points
109 if (bg::math::equals(diff, 0.0)
110 || bg::math::equals(bg::math::abs(diff), 180.0))
111 {
112 d_formula = d_strategy;
113 }
114 else
115 {
116 d_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1),
117 to_rad(lon2), to_rad(lat2),
118 stype()).distance;
119 }
120
121 BOOST_CHECK_CLOSE(d_strategy / 1000.0, expected_km, 0.001);
122 BOOST_CHECK_CLOSE(d_strategy2 / 1000.0, expected_km, 0.001);
123 BOOST_CHECK_CLOSE(d_function / 1000.0, expected_km, 0.001);
124 BOOST_CHECK_CLOSE(d_formula / 1000.0, expected_km, 0.001);
125 }
126
127 template <typename PS, typename P>
test_azimuth(double lon1,double lat1,double lon2,double lat2,double expected_azimuth_deg)128 void test_azimuth(double lon1, double lat1,
129 double lon2, double lat2,
130 double expected_azimuth_deg)
131 {
132 // Set radius type, but for integer coordinates we want to have floating point radius type
133 typedef typename bg::promote_floating_point
134 <
135 typename bg::coordinate_type<PS>::type
136 >::type rtype;
137
138 typedef bg::srs::spheroid<rtype> stype;
139 typedef bg::formula::andoyer_inverse<rtype, false, true> andoyer_inverse_type;
140
141 rtype a_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1), to_rad(lon2), to_rad(lat2), stype()).azimuth;
142
143 rtype azimuth_deg = to_deg(a_formula);
144
145 if (bg::math::equals(azimuth_deg, -180.0))
146 azimuth_deg = 180.0;
147 if (bg::math::equals(expected_azimuth_deg, -180.0))
148 expected_azimuth_deg = 180.0;
149
150 if (bg::math::equals(expected_azimuth_deg, 0.0))
151 {
152 BOOST_CHECK(bg::math::equals(azimuth_deg, expected_azimuth_deg));
153 }
154 else
155 {
156 BOOST_CHECK_CLOSE(azimuth_deg, expected_azimuth_deg, 0.001);
157 }
158 }
159
160 template <typename P1, typename P2>
test_distazi(double lon1,double lat1,double lon2,double lat2,double expected_km,double expected_azimuth_deg)161 void test_distazi(double lon1, double lat1, double lon2, double lat2,
162 double expected_km, double expected_azimuth_deg)
163 {
164 test_distance<P1, P2>(lon1, lat1, lon2, lat2, expected_km);
165 test_azimuth<P1, P2>(lon1, lat1, lon2, lat2, expected_azimuth_deg);
166 }
167
168 // requires SW->NE
169 template <typename P1, typename P2>
test_distazi_symm(double lon1,double lat1,double lon2,double lat2,double expected_km,double expected_azimuth_deg,bool is_antipodal=false)170 void test_distazi_symm(double lon1, double lat1, double lon2, double lat2,
171 double expected_km, double expected_azimuth_deg,
172 bool is_antipodal = false)
173 {
174 double d180 = is_antipodal ? 0 : 180;
175 test_distazi<P1, P2>(lon1, lat1, lon2, lat2, expected_km, expected_azimuth_deg);
176 test_distazi<P1, P2>(-lon1, lat1, -lon2, lat2, expected_km, -expected_azimuth_deg);
177 test_distazi<P1, P2>(lon1, -lat1, lon2, -lat2, expected_km, d180 - expected_azimuth_deg);
178 test_distazi<P1, P2>(-lon1, -lat1, -lon2, -lat2, expected_km, -d180 + expected_azimuth_deg);
179 }
180
181 template <typename P1, typename P2>
test_distazi_symmNS(double lon1,double lat1,double lon2,double lat2,double expected_km,double expected_azimuth_deg)182 void test_distazi_symmNS(double lon1, double lat1, double lon2, double lat2,
183 double expected_km, double expected_azimuth_deg)
184 {
185 test_distazi<P1, P2>(lon1, lat1, lon2, lat2, expected_km, expected_azimuth_deg);
186 test_distazi<P1, P2>(lon1, -lat1, lon2, -lat2, expected_km, 180 - expected_azimuth_deg);
187 }
188
189
190 template <typename PS, typename P>
test_side(double lon1,double lat1,double lon2,double lat2,double lon,double lat,int expected_side)191 void test_side(double lon1, double lat1,
192 double lon2, double lat2,
193 double lon, double lat,
194 int expected_side)
195 {
196 // Set radius type, but for integer coordinates we want to have floating point radius type
197 typedef typename bg::promote_floating_point
198 <
199 typename bg::coordinate_type<PS>::type
200 >::type rtype;
201
202 typedef bg::srs::spheroid<rtype> stype;
203
204 typedef bg::strategy::side::andoyer<stype> strategy_type;
205 typedef bg::strategy::side::geographic<bg::strategy::andoyer, stype> strategy2_type;
206
207 strategy_type strategy;
208 strategy2_type strategy2;
209
210 PS p1, p2;
211 P p;
212
213 bg::assign_values(p1, lon1, lat1);
214 bg::assign_values(p2, lon2, lat2);
215 bg::assign_values(p, lon, lat);
216
217 int side = strategy.apply(p1, p2, p);
218 int side2 = strategy2.apply(p1, p2, p);
219
220 BOOST_CHECK_EQUAL(side, expected_side);
221 BOOST_CHECK_EQUAL(side2, expected_side);
222 }
223
224 template <typename P1, typename P2>
test_all()225 void test_all()
226 {
227 // polar
228 test_distazi<P1, P2>(0, 90, 1, 80,
229 1116.814237, 179);
230
231 // no point difference
232 test_distazi<P1, P2>(4, 52, 4, 52,
233 0.0, 0.0);
234
235 // normal cases
236 test_distazi<P1, P2>(4, 52, 3, 40,
237 1336.039890, -176.3086);
238 test_distazi<P1, P2>(3, 52, 4, 40,
239 1336.039890, 176.3086);
240 test_distazi<P1, P2>(make_deg(17, 19, 43.28),
241 make_deg(40, 30, 31.151),
242 18, 40,
243 80.323245,
244 make_deg(134, 27, 50.05));
245
246 // antipodal
247 // ok? in those cases shorter path would pass through a pole
248 // but 90 or -90 would be consistent with distance?
249 test_distazi<P1, P2>(0, 0, 180, 0, 20003.9, 0.0);
250 test_distazi<P1, P2>(0, 0, -180, 0, 20003.9, 0.0);
251 test_distazi<P1, P2>(-90, 0, 90, 0, 20003.9, 0.0);
252 test_distazi<P1, P2>(90, 0, -90, 0, 20003.9, 0.0);
253
254 // 0, 45, 90 ...
255 for (int i = 0 ; i < 360 ; i += 45)
256 {
257 // 0 45 90 ...
258 double l = normlized_deg(i);
259 // -1 44 89 ...
260 double l1 = normlized_deg(i - 1);
261 // 1 46 91 ...
262 double l2 = normlized_deg(i + 1);
263
264 // near equator
265 test_distazi_symm<P1, P2>(l1, -1, l2, 1, 313.7956, 45.1964);
266
267 // near poles
268 test_distazi_symmNS<P1, P2>(l, -89.5, l, 89.5, 19892.2, 0.0);
269 test_distazi_symmNS<P1, P2>(l, -89.6, l, 89.6, 19914.6, 0.0);
270 test_distazi_symmNS<P1, P2>(l, -89.7, l, 89.7, 19936.9, 0.0);
271 test_distazi_symmNS<P1, P2>(l, -89.8, l, 89.8, 19959.2, 0.0);
272 test_distazi_symmNS<P1, P2>(l, -89.9, l, 89.9, 19981.6, 0.0);
273 test_distazi_symmNS<P1, P2>(l, -89.99, l, 89.99, 20001.7, 0.0);
274 test_distazi_symmNS<P1, P2>(l, -89.999, l, 89.999, 20003.7, 0.0);
275 // antipodal
276 test_distazi_symmNS<P1, P2>(l, -90, l, 90, 20003.9, 0.0);
277
278 test_distazi_symm<P1, P2>(normlized_deg(l-10.0), -10.0, normlized_deg(l+135), 45, 14892.1, 34.1802);
279 test_distazi_symm<P1, P2>(normlized_deg(l-30.0), -30.0, normlized_deg(l+135), 45, 17890.7, 33.7002);
280 test_distazi_symm<P1, P2>(normlized_deg(l-40.0), -40.0, normlized_deg(l+135), 45, 19319.7, 33.4801);
281 test_distazi_symm<P1, P2>(normlized_deg(l-41.0), -41.0, normlized_deg(l+135), 45, 19459.1, 33.2408);
282 test_distazi_symm<P1, P2>(normlized_deg(l-42.0), -42.0, normlized_deg(l+135), 45, 19597.8, 32.7844);
283 test_distazi_symm<P1, P2>(normlized_deg(l-43.0), -43.0, normlized_deg(l+135), 45, 19735.8, 31.7784);
284 test_distazi_symm<P1, P2>(normlized_deg(l-44.0), -44.0, normlized_deg(l+135), 45, 19873.0, 28.5588);
285 test_distazi_symm<P1, P2>(normlized_deg(l-44.1), -44.1, normlized_deg(l+135), 45, 19886.7, 27.8304);
286 test_distazi_symm<P1, P2>(normlized_deg(l-44.2), -44.2, normlized_deg(l+135), 45, 19900.4, 26.9173);
287 test_distazi_symm<P1, P2>(normlized_deg(l-44.3), -44.3, normlized_deg(l+135), 45, 19914.1, 25.7401);
288 test_distazi_symm<P1, P2>(normlized_deg(l-44.4), -44.4, normlized_deg(l+135), 45, 19927.7, 24.1668);
289 test_distazi_symm<P1, P2>(normlized_deg(l-44.5), -44.5, normlized_deg(l+135), 45, 19941.4, 21.9599);
290 test_distazi_symm<P1, P2>(normlized_deg(l-44.6), -44.6, normlized_deg(l+135), 45, 19955.0, 18.6438);
291 test_distazi_symm<P1, P2>(normlized_deg(l-44.7), -44.7, normlized_deg(l+135), 45, 19968.6, 13.1096);
292 test_distazi_symm<P1, P2>(normlized_deg(l-44.8), -44.8, normlized_deg(l+135), 45, 19982.3, 2.0300);
293 // nearly antipodal
294 test_distazi_symm<P1, P2>(normlized_deg(l-44.9), -44.9, normlized_deg(l+135), 45, 19995.9, 0.0);
295 test_distazi_symm<P1, P2>(normlized_deg(l-44.95), -44.95, normlized_deg(l+135), 45, 20002.7, 0.0);
296 test_distazi_symm<P1, P2>(normlized_deg(l-44.99), -44.99, normlized_deg(l+135), 45, 20008.1, 0.0);
297 test_distazi_symm<P1, P2>(normlized_deg(l-44.999), -44.999, normlized_deg(l+135), 45, 20009.4, 0.0);
298 // antipodal
299 test_distazi_symm<P1, P2>(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20003.92, 0.0, true);
300 }
301
302 /* SQL Server gives:
303 1116.82586908528, 0, 1336.02721932545
304
305 with:
306 SELECT 0.001 * geography::STGeomFromText('POINT(0 90)', 4326).STDistance(geography::STGeomFromText('POINT(1 80)', 4326))
307 union SELECT 0.001 * geography::STGeomFromText('POINT(4 52)', 4326).STDistance(geography::STGeomFromText('POINT(4 52)', 4326))
308 union SELECT 0.001 * geography::STGeomFromText('POINT(4 52)', 4326).STDistance(geography::STGeomFromText('POINT(3 40)', 4326))
309 */
310
311
312 test_side<P1, P2>(0, 0, 0, 1, 0, 2, 0);
313 test_side<P1, P2>(0, 0, 0, 1, 0, -2, 0);
314 test_side<P1, P2>(10, 0, 10, 1, 10, 2, 0);
315 test_side<P1, P2>(10, 0, 10, -1, 10, 2, 0);
316
317 test_side<P1, P2>(10, 0, 10, 1, 0, 2, 1); // left
318 test_side<P1, P2>(10, 0, 10, -1, 0, 2, -1); // right
319
320 test_side<P1, P2>(-10, -10, 10, 10, 10, 0, -1); // right
321 test_side<P1, P2>(-10, -10, 10, 10, -10, 0, 1); // left
322 test_side<P1, P2>(170, -10, -170, 10, -170, 0, -1); // right
323 test_side<P1, P2>(170, -10, -170, 10, 170, 0, 1); // left
324 }
325
326 template <typename P>
test_all()327 void test_all()
328 {
329 test_all<P, P>();
330 }
331
test_main(int,char * [])332 int test_main(int, char* [])
333 {
334 //test_all<float[2]>();
335 //test_all<double[2]>();
336 //test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >();
337 //test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >();
338 test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
339
340 #if defined(HAVE_TTMATH)
341 test_all<bg::model::point<ttmath::Big<1,4>, 2, bg::cs::geographic<bg::degree> > >();
342 test_all<bg::model::point<ttmath_big, 2, bg::cs::geographic<bg::degree> > >();
343 #endif
344
345 return 0;
346 }
347