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// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
8
9// This file was modified by Oracle on 2015-2022.
10// Modifications copyright (c) 2015-2022, Oracle and/or its affiliates.
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#include <boost/geometry.hpp>
22#include <geometry_test_common.hpp>
23
24#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
25#include <GeographicLib/Constants.hpp>
26#include <GeographicLib/Geodesic.hpp>
27#include <GeographicLib/PolygonArea.hpp>
28
29template <typename polygon>
30auto geographiclib_area(polygon p) {
31
32 using namespace GeographicLib;
33
34 //Geodesic geod(6371008.8, 0.0);//Constants::WGS84_f());
35 const Geodesic& geod = Geodesic::WGS84();
36 PolygonArea poly(geod);
37 for (auto it = boost::begin(boost::geometry::exterior_ring(p));
38 it != boost::end(boost::geometry::exterior_ring(p)); ++it)
39 {
40 double x = bg::get_as_radian<0>(*it) * bg::math::r2d<double>();
41 double y = bg::get_as_radian<1>(*it) * bg::math::r2d<double>();
42 poly.AddPoint( y, x);
43 }
44 double perimeter, area;
45 poly.Compute(false, true, perimeter, area);
46 // geographiclib has different orientation than BG
47 return -area;
48}
49#endif // BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
50
51namespace bg = boost::geometry;
52
53//Testing spherical and geographic strategies
54template <typename CT>
55void test_spherical_geo()
56{
57 using ct = CT;
58
59 //Geographic
60
61 using pt_geo = typename bg::model::point
62 <
63 ct, 2, bg::cs::geographic<bg::degree>
64 >;
65
66 using pt_geo_d = bg::model::point
67 <
68 double, 2, bg::cs::geographic<bg::degree>
69 >;
70 using pt_geo_ld = bg::model::point
71 <
72 long double, 2, bg::cs::geographic<bg::degree>
73 >;
74
75 bg::strategy::area::geographic<bg::strategy::vincenty, 5> area_geographic;
76 bg::strategy::area::geographic<bg::strategy::andoyer> area_a;
77 bg::strategy::area::geographic<bg::strategy::thomas> area_t;
78 bg::strategy::area::geographic<bg::strategy::vincenty> area_v;
79 bg::strategy::area::geographic<bg::strategy::karney> area_k;
80
81 bg::model::polygon<pt_geo> geometry_geo;
82
83
84 //Spherical
85
86 using pt_sph = typename bg::model::point
87 <
88 ct, 2, bg::cs::spherical_equatorial<bg::degree>
89 >;
90
91 bg::model::polygon<pt_sph> geometry;
92
93 // unit-sphere has area of 4-PI. Polygon covering 1/8 of it:
94 std::string poly = "POLYGON((0 0,0 90,90 0,0 0))";
95
96 bg::strategy::area::spherical
97 <
98 ct
99 > strategy_unary(1.0);
100
101 ct const four = 4.0;
102 ct const eight = 8.0;
103 ct expected = four * boost::geometry::math::pi<ct>() / eight;
104 bg::read_wkt(poly, geometry);
105 ct area = bg::area(geometry, strategy_unary);
106 BOOST_CHECK_CLOSE(area, expected, 0.0001);
107
108 // With strategy, radius 2 -> 4 pi r^2
109 bg::strategy::area::spherical
110 <
111 ct
112 > strategy(2.0);
113
114 area = bg::area(geometry, strategy);
115 ct const two = 2.0;
116 BOOST_CHECK_CLOSE(area, two * two * expected, 0.0001);
117
118 // Geographic total area of earth is about 510065626583900.6 (WGS84 ellipsoid)
119 // (510072000 in https://en.wikipedia.org/wiki/Earth)
120 // So the 1/8 is 6.375820332*10^13 and here we get something close to it
121 bg::read_wkt(poly, geometry_geo);
122 area = bg::area(geometry_geo, area_geographic);
123 //GeoGraphicLib gives: 63758202715511.055
124 BOOST_CHECK_CLOSE(area, 63758202715509.844, 0.0001);
125
126
127 // Wrangel Island (dateline crossing)
128 // With (spherical) Earth strategy
129 poly = "POLYGON((-178.7858 70.7852, 177.4758 71.2333, 179.7436 71.5733, -178.7858 70.7852))";
130
131 bg::strategy::area::spherical
132 <
133 ct
134 > spherical_earth(6373);
135 bg::read_wkt(poly, geometry);
136 area = bg::area(geometry, spherical_earth);
137 // SQL Server gives: 4537.9654419375
138 // PostGIS gives: 4537.9311668307
139 // Note: those are Geographic, this test is Spherical
140 BOOST_CHECK_CLOSE(area, 4506.6389, 0.001);
141
142 bg::read_wkt(poly, geometry_geo);
143 area = bg::area(geometry_geo, area_geographic);
144 BOOST_CHECK_CLOSE(area, 4537929936.5883484, 0.0001);
145#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
146 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.1);
147#endif
148 // Wrangel, more in detail
149 poly = "POLYGON((-178.568604 71.564148,-178.017548 71.449692,-177.833313 71.3461,\
150 -177.502838 71.277466 ,-177.439453 71.226929,-177.620026 71.116638,\
151 -177.9389 71.037491,-178.8186 70.979965,-179.274445 70.907761,\
152 -180 70.9972,179.678314 70.895538,179.272766 70.888596,\
153 178.791016 70.7964,178.617737 71.035538,178.872192 71.217484,\
154 179.530273 71.4383 ,-180 71.535843 ,-179.628601 71.577194,\
155 -179.305298 71.551361,-179.03421 71.597748,-178.568604 71.564148))";
156 bg::read_wkt(poly, geometry);
157 area = bg::area(geometry, spherical_earth);
158 // SQL Server gives: 7669.10402181435
159 // PostGIS gives: 7669.55565459832
160 BOOST_CHECK_CLOSE(area, 7616.523769, 0.001);
161
162 bg::read_wkt(poly, geometry_geo);
163 area = bg::area(geometry_geo, area_geographic);
164 BOOST_CHECK_CLOSE(area, 7669498457.4130802, 0.0001);
165
166 // Check more at the equator
167 /*
168 select 1,geography::STGeomFromText('POLYGON((-178.7858 10.7852 , 179.7436 11.5733 , \
169 177.4758 11.2333 , -178.7858 10.7852))',4326) .STArea()/1000000.0
170 union select 2,geography::STGeomFromText('POLYGON((-178.7858 20.7852 , 179.7436 21.5733 ,\
171 177.4758 21.2333 , -178.7858 20.7852))',4326) .STArea()/1000000.0
172 union select 3,geography::STGeomFromText('POLYGON((-178.7858 30.7852 , 179.7436 31.5733 , \
173 177.4758 31.2333 , -178.7858 30.7852))',4326) .STArea()/1000000.0
174 union select 0,geography::STGeomFromText('POLYGON((-178.7858 0.7852 , 179.7436 1.5733 , \
175 177.4758 1.2333 , -178.7858 0.7852))',4326) .STArea()/1000000.0
176 union select 4,geography::STGeomFromText('POLYGON((-178.7858 40.7852 , 179.7436 41.5733 , \
177 177.4758 41.2333 , -178.7858 40.7852))',4326) .STArea()/1000000.0
178 */
179
180 poly = "POLYGON((-178.7858 0.7852, 177.4758 1.2333, 179.7436 1.5733, -178.7858 0.7852))";
181 bg::read_wkt(poly, geometry);
182 area = bg::area(geometry, spherical_earth);
183 BOOST_CHECK_CLOSE(area, 14136.09946, 0.001); // SQL Server gives: 14064.1902284513
184 bg::read_wkt(poly, geometry_geo);
185 area = bg::area(geometry_geo, area_geographic);
186 BOOST_CHECK_CLOSE(area, 14064129044.677765, 0.0001);
187
188 poly = "POLYGON((-178.7858 10.7852, 177.4758 11.2333, 179.7436 11.5733, -178.7858 10.7852))";
189 // Geographiclib 1.51 :13696308940.35794067382812
190 bg::read_wkt(poly, geometry);
191 area = bg::area(geometry, spherical_earth);
192 BOOST_CHECK_CLOSE(area, 13760.2456, 0.001); // SQL Server gives: 13697.0941155193
193 bg::read_wkt(poly, geometry_geo);
194 area = bg::area(geometry_geo, area_geographic);
195 BOOST_CHECK_CLOSE(area, 13696308940.342197, 0.0001);
196
197 poly = "POLYGON((-178.7858 20.7852, 177.4758 21.2333, 179.7436 21.5733, -178.7858 20.7852))";
198 bg::read_wkt(poly, geometry);
199 area = bg::area(geometry, spherical_earth);
200 BOOST_CHECK_CLOSE(area, 12987.8682, 0.001); // SQL Server gives: 12944.3970990317 -> -39m^2
201 bg::read_wkt(poly, geometry_geo);
202 area = bg::area(geometry_geo, area_geographic);
203 BOOST_CHECK_CLOSE(area, 12943176284.489822, 0.0001);
204#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
205 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.00001);
206#endif
207
208 poly = "POLYGON((-178.7858 30.7852, 177.4758 31.2333, 179.7436 31.5733, -178.7858 30.7852))";
209 bg::read_wkt(poly, geometry);
210 area = bg::area(geometry, spherical_earth);
211 BOOST_CHECK_CLOSE(area, 11856.3935, 0.001); // SQL Server gives: 11838.5338423574 -> -18m^2
212 bg::read_wkt(poly, geometry_geo);
213 area = bg::area(geometry_geo, area_geographic);
214 BOOST_CHECK_CLOSE(area, 11837280445.394035, 0.0001);
215#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
216 // GEOGRAPHICLIB 1.51: 11837280466.204895
217 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.000001);
218#endif
219
220 poly = "POLYGON((-178.7858 40.7852, 177.4758 41.2333, 179.7436 41.5733, -178.7858 40.7852))";
221 bg::read_wkt(poly, geometry);
222 area = bg::area(geometry, spherical_earth);
223 BOOST_CHECK_CLOSE(area, 10404.627685523914, 0.001);
224 // SQL Server gives: 10412.0607137119, -> +8m^2
225 bg::read_wkt(poly, geometry_geo);
226 area = bg::area(geometry_geo, area_geographic);
227 BOOST_CHECK_CLOSE(area, 10411098789.387386, 0.0001);
228#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
229 // GEOGRAPHICLIB 1.51: 10411098790.577026
230 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.000001);
231#endif
232
233 // Concave
234 poly = "POLYGON((0 40,1 42,0 44,2 43,4 44,3 42,4 40,2 41,0 40))";
235 bg::read_wkt(poly, geometry);
236 area = bg::area(geometry, spherical_earth);
237 BOOST_CHECK_CLOSE(area, 73538.2958, 0.001); // SQL Server gives: 73604.2047689719
238 bg::read_wkt(poly, geometry_geo);
239 area = bg::area(geometry_geo, area_geographic);
240 BOOST_CHECK_CLOSE(area, 73604163095.545319, 0.0001);
241#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
242 // GEOGRAPHICLIB 1.51: 73604163091.274048
243 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.000001);
244#endif
245
246 // With hole POLYGON((0 40,4 40,4 44,0 44,0 40),(1 41,2 43,3 42,1 41))
247 poly = "POLYGON((0 40,0 44,4 44,4 40,0 40),(1 41,3 42,2 43,1 41))";
248 bg::read_wkt(poly, geometry);
249 area = bg::area(geometry, spherical_earth);
250 BOOST_CHECK_CLOSE(area, 133233.844876, 0.001); // SQL Server gives: 133353.335
251 bg::read_wkt(poly, geometry_geo);
252 area = bg::area(geometry_geo, area_geographic);
253 BOOST_CHECK_CLOSE(area, 133353077343.26692, 0.0001);
254#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
255 // GEOGRAPHICLIB 1.51: 147189206350.20142
256 //BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.9);
257#endif
258
259 // mean Earth's radius^2
260 double r2 = bg::math::sqr(value: bg::get_radius<0>(geometry: bg::srs::sphere<double>()));
261
262 // around 0 meridian
263 {
264 std::string poly1 = "POLYGON((-10 0,-10 10,0 10,0 0,-10 0))";
265 std::string poly2 = "POLYGON((0 0,0 10,10 10,10 0,0 0))";
266 std::string poly3 = "POLYGON((-5 0,-5 10,5 10,5 0,-5 0))";
267 bg::read_wkt(poly1, geometry);
268 ct area1 = bg::area(geometry);
269 bg::read_wkt(poly2, geometry);
270 ct area2 = bg::area(geometry);
271 bg::read_wkt(poly3, geometry);
272 ct area3 = bg::area(geometry);
273 BOOST_CHECK_CLOSE(area1, area2, 0.001);
274 BOOST_CHECK_CLOSE(area2, area3, 0.001);
275 BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1848, 0.001);
276 //geographic
277 bg::read_wkt(poly1, geometry_geo);
278 area1 = bg::area(geometry_geo, area_geographic);
279 bg::read_wkt(poly2, geometry_geo);
280 area2 = bg::area(geometry_geo, area_geographic);
281 bg::read_wkt(poly3, geometry_geo);
282 area3 = bg::area(geometry_geo, area_geographic);
283 BOOST_CHECK_CLOSE(area1, area2, 0.001);
284 BOOST_CHECK_CLOSE(area2, area3, 0.001);
285 BOOST_CHECK_CLOSE(area1, 1227877191611.3574, 0.001);
286#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
287 // GEOGRAPHICLIB 1.51: 1227877191609.6267
288 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area3, 0.000001);
289#endif
290 }
291 {
292 std::string poly1 = "POLYGON((-10 -5,-10 5,0 5,0 -5,-10 -5))";
293 std::string poly2 = "POLYGON((0 -5,0 5,10 5,10 -5,0 -5))";
294 std::string poly3 = "POLYGON((-5 -5,-5 5,5 5,5 -5,-5 -5))";
295 bg::read_wkt(poly1, geometry);
296 ct area1 = bg::area(geometry);
297 bg::read_wkt(poly2, geometry);
298 ct area2 = bg::area(geometry);
299 bg::read_wkt(poly3, geometry);
300 ct area3 = bg::area(geometry);
301 BOOST_CHECK_CLOSE(area1, area2, 0.001);
302 BOOST_CHECK_CLOSE(area2, area3, 0.001);
303 BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0261, 0.001);
304 //geographic
305 bg::read_wkt(poly1, geometry_geo);
306 area1 = bg::area(geometry_geo, area_geographic);
307 bg::read_wkt(poly2, geometry_geo);
308 area2 = bg::area(geometry_geo, area_geographic);
309 bg::read_wkt(poly3, geometry_geo);
310 area3 = bg::area(geometry_geo, area_geographic);
311 BOOST_CHECK_CLOSE(area1, area2, 0.001);
312 BOOST_CHECK_CLOSE(area2, area3, 0.001);
313 BOOST_CHECK_CLOSE(area1, 1232514639151.7168, 0.001);
314#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
315 // GEOGRAPHICLIB 1.51: 1232514639151.6357
316 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area3, 0.000001);
317#endif
318 }
319 // around 180 meridian
320 {
321 std::string poly1 = "POLYGON((-180 0,-180 10,-170 10,-170 0,-180 0))";
322 std::string poly2 = "POLYGON((175 0,175 10,-175 10,-175 0,175 0))";
323 std::string poly3 = "POLYGON((170 0,170 10,180 10,180 0,170 0))";
324 std::string poly4 = "POLYGON((170 0,170 10,-180 10,-180 0,170 0))";
325 std::string poly5 = "POLYGON((180 0,180 10,-170 10,-170 0,180 0))";
326 bg::read_wkt(poly1, geometry);
327 ct area1 = bg::area(geometry);
328 bg::read_wkt(poly2, geometry);
329 ct area2 = bg::area(geometry);
330 bg::read_wkt(poly3, geometry);
331 ct area3 = bg::area(geometry);
332 bg::read_wkt(poly4, geometry);
333 ct area4 = bg::area(geometry);
334 bg::read_wkt(poly5, geometry);
335 ct area5 = bg::area(geometry);
336 BOOST_CHECK_CLOSE(area1, area2, 0.001);
337 BOOST_CHECK_CLOSE(area2, area3, 0.001);
338 BOOST_CHECK_CLOSE(area3, area4, 0.001);
339 BOOST_CHECK_CLOSE(area4, area5, 0.001);
340 BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1833, 0.001);
341 //geographic
342 bg::read_wkt(poly1, geometry_geo);
343 area1 = bg::area(geometry_geo, area_geographic);
344 bg::read_wkt(poly2, geometry_geo);
345 area2 = bg::area(geometry_geo, area_geographic);
346 bg::read_wkt(poly3, geometry_geo);
347 area3 = bg::area(geometry_geo, area_geographic);
348 bg::read_wkt(poly4, geometry_geo);
349 area4 = bg::area(geometry_geo, area_geographic);
350 bg::read_wkt(poly5, geometry_geo);
351 area5 = bg::area(geometry_geo, area_geographic);
352 BOOST_CHECK_CLOSE(area1, area2, 0.001);
353 BOOST_CHECK_CLOSE(area2, area3, 0.001);
354 BOOST_CHECK_CLOSE(area3, area4, 0.001);
355 BOOST_CHECK_CLOSE(area4, area5, 0.001);
356 BOOST_CHECK_CLOSE(area1, 1227877191611.3574, 0.001);
357#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
358 // GEOGRAPHICLIB 1.51: 1227877191609.6267
359 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area5, 0.000001);
360#endif
361 }
362 {
363 std::string poly1 = "POLYGON((-180 -5,-180 5,-170 5,-170 -5,-180 -5))";
364 std::string poly2 = "POLYGON((175 -5,175 5,-175 5,-175 -5,175 -5))";
365 std::string poly3 = "POLYGON((170 -5,170 5,180 5,180 -5,170 -5))";
366 std::string poly4 = "POLYGON((170 -5,170 5,180 5,180 -5,170 -5))";
367 std::string poly5 = "POLYGON((180 -5,180 5,-170 5,-170 -5,180 -5))";
368 bg::read_wkt(poly1, geometry);
369 ct area1 = bg::area(geometry);
370 bg::read_wkt(poly2, geometry);
371 ct area2 = bg::area(geometry);
372 bg::read_wkt(poly3, geometry);
373 ct area3 = bg::area(geometry);
374 bg::read_wkt(poly4, geometry);
375 ct area4 = bg::area(geometry);
376 bg::read_wkt(poly5, geometry);
377 ct area5 = bg::area(geometry);
378 BOOST_CHECK_CLOSE(area1, area2, 0.001);
379 BOOST_CHECK_CLOSE(area2, area3, 0.001);
380 BOOST_CHECK_CLOSE(area3, area4, 0.001);
381 BOOST_CHECK_CLOSE(area4, area5, 0.001);
382 BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0247, 0.001);
383 //geographic
384 bg::read_wkt(poly1, geometry_geo);
385 area1 = bg::area(geometry_geo, area_geographic);
386 bg::read_wkt(poly2, geometry_geo);
387 area2 = bg::area(geometry_geo, area_geographic);
388 bg::read_wkt(poly3, geometry_geo);
389 area3 = bg::area(geometry_geo, area_geographic);
390 bg::read_wkt(poly4, geometry_geo);
391 area4 = bg::area(geometry_geo, area_geographic);
392 bg::read_wkt(poly5, geometry_geo);
393 area5 = bg::area(geometry_geo, area_geographic);
394 BOOST_CHECK_CLOSE(area1, area2, 0.001);
395 BOOST_CHECK_CLOSE(area2, area3, 0.001);
396 BOOST_CHECK_CLOSE(area3, area4, 0.001);
397 BOOST_CHECK_CLOSE(area4, area5, 0.001);
398 BOOST_CHECK_CLOSE(area1, 1232514639151.7168, 0.001);
399#ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
400 // GEOGRAPHICLIB 1.51: 1232514639151.6357
401 BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area5, 0.000001);
402#endif
403 }
404 // around poles
405 {
406 std::string poly1 = "POLYGON((0 80,-90 80,-180 80,90 80,0 80))";
407 std::string poly2 = "POLYGON((0 80,-90 80,180 80,90 80,0 80))";
408 std::string poly3 = "POLYGON((0 -80,90 -80,-180 -80,-90 -80,0 -80))";
409 std::string poly4 = "POLYGON((0 -80,90 -80,180 -80,-90 -80,0 -80))";
410 bg::read_wkt(poly1, geometry);
411 ct area1 = bg::area(geometry);
412 bg::read_wkt(poly2, geometry);
413 ct area2 = bg::area(geometry);
414 bg::read_wkt(poly3, geometry);
415 ct area3 = bg::area(geometry);
416 bg::read_wkt(poly4, geometry);
417 ct area4 = bg::area(geometry);
418 BOOST_CHECK_CLOSE(area1, area2, 0.001);
419 BOOST_CHECK_CLOSE(area2, area3, 0.001);
420 BOOST_CHECK_CLOSE(area3, area4, 0.001);
421 //geographic
422 bg::read_wkt(poly1, geometry_geo);
423 area1 = bg::area(geometry_geo, area_geographic);
424 bg::read_wkt(poly2, geometry_geo);
425 area2 = bg::area(geometry_geo, area_geographic);
426 bg::read_wkt(poly3, geometry_geo);
427 area3 = bg::area(geometry_geo, area_geographic);
428 bg::read_wkt(poly4, geometry_geo);
429 area4 = bg::area(geometry_geo, area_geographic);
430 BOOST_CHECK_CLOSE(area1, area2, 0.001);
431 BOOST_CHECK_CLOSE(area2, area3, 0.001);
432 BOOST_CHECK_CLOSE(area3, area4, 0.001);
433 }
434
435 {
436 bg::model::ring<pt_sph> aurha; // a'dam-utr-rott.-den haag-a'dam
437 std::string poly = "POLYGON((4.892 52.373,5.119 52.093,4.479 51.930,\
438 4.23 52.08,4.892 52.373))";
439 bg::read_wkt(poly, aurha);
440 /*if (polar)
441 {
442 // Create colatitudes (measured from pole)
443 for (pt_sph& p : aurha)
444 {
445 bg::set<1>(p, ct(90) - bg::get<1>(p));
446 }
447 bg::correct(aurha);
448 }*/
449 bg::strategy::area::spherical
450 <
451 ct
452 > area_spherical(6372.795);
453 area = bg::area(aurha, area_spherical);
454 BOOST_CHECK_CLOSE(area, 1476.645675, 0.0001);
455 //geographic
456 bg::read_wkt(poly, geometry_geo);
457 area = bg::area(geometry_geo, area_geographic);
458 BOOST_CHECK_CLOSE(area, 1481555970.0765088, 0.001);
459
460 // SQL Server gives: 1481.55595960659
461 // for select geography::STGeomFromText('POLYGON((4.892 52.373,4.23 52.08,
462 // 4.479 51.930,5.119 52.093,4.892 52.373))',4326).STArea()/1000000.0
463 }
464
465 {
466 bg::model::polygon<pt_sph, false> geometry_sph;
467 std::string wkt = "POLYGON((0 0, 5 0, 5 5, 0 5, 0 0))";
468 bg::read_wkt(wkt, geometry_sph);
469
470 area = bg::area(geometry_sph, bg::strategy::area::spherical<>(6371228.0));
471 BOOST_CHECK_CLOSE(area, 308932296103.83051, 0.0001);
472
473 bg::model::polygon<pt_geo, false> geometry_geo;
474 bg::read_wkt(wkt, geometry_geo);
475
476 area = bg::area(geometry_geo, bg::strategy::area::geographic<>(bg::srs::spheroid<double>(6371228.0, 6371228.0)));
477 BOOST_CHECK_CLOSE(area, 308932296103.82574, 0.001);
478 }
479
480 {
481 // small areas: ~20m^2
482 // see https://github.com/boostorg/geometry/issues/799
483 // geographiclib 1.50.1 returns: 25.5736
484
485 bg::model::polygon<pt_sph> geometry;
486 std::string wkt ="POLYGON((-0.020000 51.470027,-0.020031 51.470019,-0.020043 51.470000,\
487 -0.020031 51.469981,-0.020000 51.469973,-0.019969 51.469981,\
488 -0.019957 51.470000,-0.019969 51.470019,-0.020000 51.470027))";
489 bg::read_wkt(wkt, geometry);
490
491 bg::strategy::area::spherical
492 <
493 ct
494 > area_spherical(6371228.0);
495 area = bg::area(geometry, area_spherical);
496 BOOST_CHECK_CLOSE(area, -25.48014, 0.0001);
497
498 //geographic
499
500 bg::model::polygon<pt_geo_d> geometry_geo_d;
501
502 bg::read_wkt(wkt, geometry&: geometry_geo_d);
503 area = bg::area(geometry: geometry_geo_d, strategy: area_a);
504 BOOST_CHECK_CLOSE(area, -25.47837, 0.001);
505 area = bg::area(geometry: geometry_geo_d, strategy: area_t);
506 BOOST_CHECK_CLOSE(area, -25.57355, 0.001);
507 area = bg::area(geometry: geometry_geo_d, strategy: area_v);
508 BOOST_CHECK_CLOSE(area, -25.55881, 0.001);
509 area = bg::area(geometry: geometry_geo_d, strategy: area_k);
510 BOOST_CHECK_CLOSE(area, -25.57357, 0.001);
511
512 bg::model::polygon<pt_geo_ld> geometry_geo_ld;
513
514 bg::read_wkt(wkt, geometry&: geometry_geo_ld);
515 area = bg::area(geometry: geometry_geo_ld, strategy: area_a);
516 BOOST_CHECK_CLOSE(area, -25.57978, 0.4); // -25.478374 with vc-14.1
517 area = bg::area(geometry: geometry_geo_ld, strategy: area_t);
518 BOOST_CHECK_CLOSE(area, -25.57359, 0.001);
519 area = bg::area(geometry: geometry_geo_ld, strategy: area_v);
520 BOOST_CHECK_CLOSE(area, -25.57394, 0.06); // -25.558816 with vc-14.1
521 area = bg::area(geometry: geometry_geo_ld, strategy: area_k);
522 BOOST_CHECK_CLOSE(area, -25.57359, 0.001);
523 }
524
525 {
526 // very small areas: ~2m^2
527 // geographiclib 1.50.1 returns: 2.24581
528
529 bg::model::polygon<pt_sph> geometry;
530 std::string wkt ="POLYGON((0.020000 51.470027,0.020005 51.470027,\
531 0.020005 51.470000,0.020007 51.4699,0.020000 51.470027))";
532 bg::read_wkt(wkt, geometry);
533
534 bg::strategy::area::spherical
535 <
536 ct
537 > area_spherical(6371228.0);
538 area = bg::area(geometry, area_spherical);
539 BOOST_CHECK_CLOSE(area, 2.2375981058307093, 0.001);
540
541 //geographic
542
543 bg::model::polygon<pt_geo_d> geometry_geo_d;
544
545 bg::read_wkt(wkt, geometry&: geometry_geo_d);
546 area = bg::area(geometry: geometry_geo_d, strategy: area_a);
547 BOOST_CHECK_CLOSE(area, 2.2374430036130444, 0.001);
548 area = bg::area(geometry: geometry_geo_d, strategy: area_t);
549 BOOST_CHECK_CLOSE(area, 2.245814319435655, 0.001);
550 area = bg::area(geometry: geometry_geo_d, strategy: area_v);
551 BOOST_CHECK_CLOSE(area, 2.2374430036130444, 0.001);
552 area = bg::area(geometry: geometry_geo_d, strategy: area_k);
553 BOOST_CHECK_CLOSE(area, 2.2458140167194878, 0.001);
554
555 bg::model::polygon<pt_geo_ld> geometry_geo_ld;
556
557 bg::read_wkt(wkt, geometry&: geometry_geo_ld);
558 area = bg::area(geometry: geometry_geo_ld, strategy: area_a);
559 BOOST_CHECK_CLOSE(area, 2.2374430039839437, 0.001);
560 area = bg::area(geometry: geometry_geo_ld, strategy: area_t);
561 BOOST_CHECK_CLOSE(area, 2.2457934740573235, 0.001);
562 area = bg::area(geometry: geometry_geo_ld, strategy: area_v);
563 BOOST_CHECK_CLOSE(area, 2.2374430039839437, 0.001);
564 area = bg::area(geometry: geometry_geo_ld, strategy: area_k);
565 BOOST_CHECK_CLOSE(area, 2.2458050314935392, 0.001);
566 }
567}
568
569template <typename Poly, bool north>
570Poly generate_polygon_1()
571{
572 using PT = typename bg::point_type<Poly>::type;
573 int const sign = north ? 1 : -1;
574
575 // Rotated by 90.0 - 0.00000000000001
576 Poly sp;
577 sp.outer().push_back(PT(-9.00000000000000142e+01,sign*6.00172345808800074e+01));
578 sp.outer().push_back(PT(9.00000000000000142e+01,sign*6.00172345808800074e+01));
579 sp.outer().push_back(PT(8.87671232876712253e+01,sign*6.00000000000000071e+01));
580 sp.outer().push_back(PT(8.38356164383561691e+01,sign*6.00000000000000071e+01));
581 sp.outer().push_back(PT(7.89041095890410986e+01,sign*6.00000000000000071e+01));
582 sp.outer().push_back(PT(7.39726027397260282e+01,sign*6.00000000000000071e+01));
583 sp.outer().push_back(PT(6.90410958904109719e+01,sign*6.00000000000000071e+01));
584 sp.outer().push_back(PT(6.41095890410958873e+01,sign*6.00000000000000071e+01));
585 sp.outer().push_back(PT(5.91780821917808240e+01,sign*6.00000000000000071e+01));
586 sp.outer().push_back(PT(5.42465753424657606e+01,sign*6.00000000000000071e+01));
587 sp.outer().push_back(PT(4.93150684931506831e+01,sign*6.00000000000000071e+01));
588 sp.outer().push_back(PT(4.43835616438356126e+01,sign*6.00000000000000071e+01));
589 sp.outer().push_back(PT(3.94520547945205493e+01,sign*6.00000000000000071e+01));
590 sp.outer().push_back(PT(3.45205479452054860e+01,sign*6.00000000000000071e+01));
591 sp.outer().push_back(PT(2.95890410958904120e+01,sign*6.00000000000000071e+01));
592 sp.outer().push_back(PT(2.46575342465753415e+01,sign*6.00000000000000071e+01));
593 sp.outer().push_back(PT(1.97260273972602747e+01,sign*6.00000000000000071e+01));
594 sp.outer().push_back(PT(1.47945205479452060e+01,sign*6.00000000000000071e+01));
595 sp.outer().push_back(PT(9.86301369863013733e+00,sign*6.00000000000000071e+01));
596 sp.outer().push_back(PT(4.93150684931506866e+00,sign*6.00000000000000071e+01));
597 sp.outer().push_back(PT(0.00000000000000000e+00,sign*6.00000000000000071e+01));
598 sp.outer().push_back(PT(-9.86301369863014976e+00,sign*6.00000000000000071e+01));
599 sp.outer().push_back(PT(-1.47945205479451669e+01,sign*6.00000000000000071e+01));
600 sp.outer().push_back(PT(-1.97260273972602853e+01,sign*6.00000000000000071e+01));
601 sp.outer().push_back(PT(-2.46575342465753522e+01,sign*6.00000000000000071e+01));
602 sp.outer().push_back(PT(-2.95890410958904191e+01,sign*6.00000000000000071e+01));
603 sp.outer().push_back(PT(-3.45205479452054860e+01,sign*6.00000000000000071e+01));
604 sp.outer().push_back(PT(-3.94520547945206062e+01,sign*6.00000000000000071e+01));
605 sp.outer().push_back(PT(-4.43835616438356269e+01,sign*6.00000000000000071e+01));
606 sp.outer().push_back(PT(-4.93150684931507470e+01,sign*6.00000000000000071e+01));
607 sp.outer().push_back(PT(-5.42465753424657606e+01,sign*6.00000000000000071e+01));
608 sp.outer().push_back(PT(-5.91780821917808240e+01,sign*6.00000000000000071e+01));
609 sp.outer().push_back(PT(-6.41095890410958873e+01,sign*6.00000000000000071e+01));
610 sp.outer().push_back(PT(-6.90410958904110146e+01,sign*6.00000000000000071e+01));
611 sp.outer().push_back(PT(-7.39726027397260282e+01,sign*6.00000000000000071e+01));
612 sp.outer().push_back(PT(-7.89041095890410418e+01,sign*6.00000000000000071e+01));
613 sp.outer().push_back(PT(-8.38356164383561691e+01,sign*6.00000000000000071e+01));
614 sp.outer().push_back(PT(-8.87671232876712253e+01,sign*6.00000000000000071e+01));
615 sp.outer().push_back(PT(-9.00000000000000142e+01,sign*6.00172345808800074e+01));
616 return sp;
617}
618
619template <typename Poly, bool north>
620Poly generate_polygon_2()
621{
622 using PT = typename bg::point_type<Poly>::type;
623 int const sign = north ? 1 : -1;
624
625 // Rotated By 90
626 Poly sp;
627 sp.outer().push_back(PT(-9.00000000000000000e+01,sign*6.00172345808800074e+01));
628 sp.outer().push_back(PT(9.00000000000000000e+01,sign*6.00172345808800074e+01));
629 sp.outer().push_back(PT(8.87671232876712253e+01,sign*6.00000000000000071e+01));
630 sp.outer().push_back(PT(8.38356164383561691e+01,sign*6.00000000000000071e+01));
631 sp.outer().push_back(PT(7.89041095890410986e+01,sign*6.00000000000000071e+01));
632 sp.outer().push_back(PT(7.39726027397260282e+01,sign*6.00000000000000071e+01));
633 sp.outer().push_back(PT(6.90410958904109719e+01,sign*6.00000000000000071e+01));
634 sp.outer().push_back(PT(6.41095890410958873e+01,sign*6.00000000000000071e+01));
635 sp.outer().push_back(PT(5.91780821917808240e+01,sign*6.00000000000000071e+01));
636 sp.outer().push_back(PT(5.42465753424657606e+01,sign*6.00000000000000071e+01));
637 sp.outer().push_back(PT(4.93150684931506831e+01,sign*6.00000000000000071e+01));
638 sp.outer().push_back(PT(4.43835616438356126e+01,sign*6.00000000000000071e+01));
639 sp.outer().push_back(PT(3.94520547945205493e+01,sign*6.00000000000000071e+01));
640 sp.outer().push_back(PT(3.45205479452054860e+01,sign*6.00000000000000071e+01));
641 sp.outer().push_back(PT(2.95890410958904120e+01,sign*6.00000000000000071e+01));
642 sp.outer().push_back(PT(2.46575342465753415e+01,sign*6.00000000000000071e+01));
643 sp.outer().push_back(PT(1.97260273972602747e+01,sign*6.00000000000000071e+01));
644 sp.outer().push_back(PT(1.47945205479452060e+01,sign*6.00000000000000071e+01));
645 sp.outer().push_back(PT(9.86301369863013733e+00,sign*6.00000000000000071e+01));
646 sp.outer().push_back(PT(4.93150684931506866e+00,sign*6.00000000000000071e+01));
647 sp.outer().push_back(PT(0.00000000000000000e+00,sign*6.00000000000000071e+01));
648 sp.outer().push_back(PT(-9.86301369863014976e+00,sign*6.00000000000000071e+01));
649 sp.outer().push_back(PT(-1.47945205479451669e+01,sign*6.00000000000000071e+01));
650 sp.outer().push_back(PT(-1.97260273972602853e+01,sign*6.00000000000000071e+01));
651 sp.outer().push_back(PT(-2.46575342465753522e+01,sign*6.00000000000000071e+01));
652 sp.outer().push_back(PT(-2.95890410958904191e+01,sign*6.00000000000000071e+01));
653 sp.outer().push_back(PT(-3.45205479452054860e+01,sign*6.00000000000000071e+01));
654 sp.outer().push_back(PT(-3.94520547945206062e+01,sign*6.00000000000000071e+01));
655 sp.outer().push_back(PT(-4.43835616438356269e+01,sign*6.00000000000000071e+01));
656 sp.outer().push_back(PT(-4.93150684931507470e+01,sign*6.00000000000000071e+01));
657 sp.outer().push_back(PT(-5.42465753424657606e+01,sign*6.00000000000000071e+01));
658 sp.outer().push_back(PT(-5.91780821917808240e+01,sign*6.00000000000000071e+01));
659 sp.outer().push_back(PT(-6.41095890410958873e+01,sign*6.00000000000000071e+01));
660 sp.outer().push_back(PT(-6.90410958904110146e+01,sign*6.00000000000000071e+01));
661 sp.outer().push_back(PT(-7.39726027397260282e+01,sign*6.00000000000000071e+01));
662 sp.outer().push_back(PT(-7.89041095890410418e+01,sign*6.00000000000000071e+01));
663 sp.outer().push_back(PT(-8.38356164383561691e+01,sign*6.00000000000000071e+01));
664 sp.outer().push_back(PT(-8.87671232876712253e+01,sign*6.00000000000000071e+01));
665 sp.outer().push_back(PT(-9.00000000000000000e+01,sign*6.00172345808800074e+01));
666 return sp;
667}
668
669template <typename Poly, bool north>
670Poly generate_polygon_3()
671{
672 using PT = typename bg::point_type<Poly>::type;
673 int const sign = north ? 1 : -1;
674
675 Poly sp;
676 sp.outer().push_back(PT(-9.00000000000000142e+01,sign*6.00172345808800074e+01));
677 sp.outer().push_back(PT(9.00000000000000142e+01,sign*6.00172345808800074e+01));
678 sp.outer().push_back(PT(0.00000000000000000e+00,sign*6.00000000000000071e+01));
679 sp.outer().push_back(PT(-9.00000000000000142e+01,sign*6.00172345808800074e+01));
680 return sp;
681}
682
683template <typename Poly, bool north>
684Poly generate_polygon_4()
685{
686 using PT = typename bg::point_type<Poly>::type;
687 int const sign = north ? 1 : -1;
688
689 Poly sp;
690 sp.outer().push_back(PT(-9.00000000000000000e+01,sign*6.00172345808800074e+01));
691 sp.outer().push_back(PT(9.00000000000000000e+01,sign*6.00172345808800074e+01));
692 sp.outer().push_back(PT(0.00000000000000000e+00,sign*6.00000000000000071e+01));
693 sp.outer().push_back(PT(-9.00000000000000000e+01,sign*6.00172345808800074e+01));
694 return sp;
695}
696
697template <typename Poly, bool north>
698Poly generate_polygon_5()
699{
700 using PT = typename bg::point_type<Poly>::type;
701 int const sign = north ? 1 : -1;
702
703 Poly sp;
704 sp.outer().push_back(PT(-80,sign*6.00172345808800074e+01));
705 sp.outer().push_back(PT(80,sign*6.00172345808800074e+01));
706 sp.outer().push_back(PT(0,sign*6.00000000000000071e+01));
707 sp.outer().push_back(PT(-80,sign*6.00172345808800074e+01));
708 return sp;
709}
710
711template <typename Poly, bool north>
712Poly generate_polygon_6()
713{
714 using PT = typename bg::point_type<Poly>::type;
715 int const sign = north ? 1 : -1;
716
717 Poly sp;
718 sp.outer().push_back(PT(-8.999999999999142e+01,sign*6.00172345808800074e+01));
719 sp.outer().push_back(PT(8.999999999999142e+01,sign*6.00172345808800074e+01));
720 sp.outer().push_back(PT(0.00000000000000000e+00,sign*6.00000000000000071e+01));
721 sp.outer().push_back(PT(-8.999999999999142e+01,sign*6.00172345808800074e+01));
722 return sp;
723}
724
725template <typename Poly, bool north>
726Poly generate_polygon_7()
727{
728 using PT = typename bg::point_type<Poly>::type;
729 int const sign = north ? 1 : -1;
730
731 Poly sp;
732 sp.outer().push_back(PT(-80,sign*6.00172345808800074e+01));
733 sp.outer().push_back(PT(100,sign*6.00172345808800074e+01));
734 sp.outer().push_back(PT(10,sign*6.00000000000000071e+01));
735 sp.outer().push_back(PT(-80,sign*6.00172345808800074e+01));
736 return sp;
737}
738
739template <typename Poly, bool north>
740Poly generate_polygon_8()
741{
742 using PT = typename bg::point_type<Poly>::type;
743 int const sign = north ? 1 : -1;
744
745 Poly sp;
746 sp.outer().push_back(PT(-80,sign*6.00172345808800074e+01));
747 sp.outer().push_back(PT(100,sign*6.00172345808800074e+01));
748 sp.outer().push_back(PT(30,sign*6.00000000000000071e+01));
749 sp.outer().push_back(PT(10,sign*6.00000000000000071e+01));
750 sp.outer().push_back(PT(-80,sign*6.00172345808800074e+01));
751 return sp;
752}
753
754template <typename Poly, bool north>
755Poly generate_polygon_9()
756{
757 using PT = typename bg::point_type<Poly>::type;
758 int const sign = north ? 1 : -1;
759
760 Poly sp;
761 sp.outer().push_back(PT(-90,sign*6.00172345808800074e+01));
762 sp.outer().push_back(PT(90.0000000000000000000001,sign*6.00172345808800074e+01));
763 sp.outer().push_back(PT(0,sign*6.00000000000000071e+01));
764 sp.outer().push_back(PT(-90,sign*6.00172345808800074e+01));
765 return sp;
766}
767
768template <typename Poly, bool north>
769Poly generate_polygon_10()
770{
771 using PT = typename bg::point_type<Poly>::type;
772 int const sign = north ? 1 : -1;
773
774 Poly sp;
775 sp.outer().push_back(PT(-90,sign*6.00172345808800074e+01));
776 sp.outer().push_back(PT(89.9999999999999999,sign*6.00172345808800074e+01));
777 sp.outer().push_back(PT(0,sign*6.00000000000000071e+01));
778 sp.outer().push_back(PT(-90,sign*6.00172345808800074e+01));
779 return sp;
780}
781
782template <typename PT, bool north, typename Strategy>
783void segment_through_pole(Strategy str, std::vector<double> results)
784{
785 // special cases with segments that pass through (or nearby) a pole
786 // see issue https://github.com/boostorg/geometry/issues/1063
787
788 using polygon = boost::geometry::model::polygon
789 <
790 PT,
791 north, // clockwise
792 true, // closed
793 std::vector,
794 std::vector,
795 std::allocator,
796 std::allocator
797 >;
798
799 {
800 auto sp = generate_polygon_1<polygon, north>();
801 auto area = bg::area(sp, str);
802 BOOST_CHECK_CLOSE(area, results[0], 0.001);
803 }
804 {
805 auto sp = generate_polygon_2<polygon, north>();
806 auto area = bg::area(sp, str);
807 BOOST_CHECK_CLOSE(area, results[0], 0.001);
808 }
809 {
810 auto sp = generate_polygon_3<polygon, north>();
811 auto area = bg::area(sp, str);
812 BOOST_CHECK_CLOSE(area, results[1], 0.001);
813 }
814 {
815 auto sp = generate_polygon_4<polygon, north>();
816 auto area = bg::area(sp, str);
817 BOOST_CHECK_CLOSE(area, results[1], 0.001);
818 }
819 {
820 auto sp = generate_polygon_5<polygon, north>();
821 auto area = bg::area(sp, str);
822 BOOST_CHECK_CLOSE(area, results[2], 0.001);
823 }
824 {
825 auto sp = generate_polygon_6<polygon, north>();
826 auto area = bg::area(sp, str);
827 BOOST_CHECK_CLOSE(area, results[1], 0.001);
828 }
829 {
830 auto sp = generate_polygon_7<polygon, north>();
831 auto area = bg::area(sp, str);
832 BOOST_CHECK_CLOSE(area, results[1], 0.001);
833 }
834 {
835 auto sp = generate_polygon_8<polygon, north>();
836 auto area = bg::area(sp, str);
837 BOOST_CHECK_CLOSE(area, results[3], 0.001);
838 }
839 {
840 auto sp = generate_polygon_9<polygon, north>();
841 auto area = bg::area(sp, str);
842 BOOST_CHECK_CLOSE(area, results[1], 0.001);
843 }
844 {
845 auto sp = generate_polygon_10<polygon, north>();
846 auto area = bg::area(sp, str);
847 BOOST_CHECK_CLOSE(area, results[1], 0.001);
848 }
849}
850
851int test_main(int, char* [])
852{
853
854 test_spherical_geo<double>();
855
856 using sph_pt = typename bg::model::point
857 <
858 double, 2, bg::cs::spherical_equatorial<bg::degree>
859 >;
860 bg::strategy::area::spherical<> sph_str;
861
862 std::vector<double> results_sph {0.4204069, 0.2865233, 0.2261384, 0.3206943};
863
864 bool const north = true;
865 bool const south = false;
866 segment_through_pole<sph_pt, north>(str: sph_str, results: results_sph);
867 segment_through_pole<sph_pt, south>(str: sph_str, results: results_sph);
868
869
870 using geo_pt = typename bg::model::point
871 <
872 double, 2, bg::cs::geographic<bg::degree>
873 >;
874
875 bg::strategy::area::geographic<bg::strategy::andoyer> area_a;
876 bg::strategy::area::geographic<bg::strategy::thomas> area_t;
877 bg::strategy::area::geographic<bg::strategy::vincenty> area_v;
878 bg::strategy::area::geographic<bg::strategy::karney> area_k;
879
880 std::vector<double> results_geo {17188025916331.664, 11713812318907.75, 9245554289678.4492,
881 13111012015900.15};
882
883 segment_through_pole<geo_pt, north>(str: area_a, results: results_geo);
884 segment_through_pole<geo_pt, north>(str: area_v, results: results_geo);
885 //NOTE: thomas strategy results in wrong results in some cases
886 //segment_through_pole<geo_pt>(area_t, results_geo);
887 segment_through_pole<geo_pt, north>(str: area_k, results: results_geo);
888
889 segment_through_pole<geo_pt, south>(str: area_a, results: results_geo);
890 segment_through_pole<geo_pt, south>(str: area_v, results: results_geo);
891 segment_through_pole<geo_pt, south>(str: area_k, results: results_geo);
892
893 return 0;
894}
895

source code of boost/libs/geometry/test/algorithms/area/area_sph_geo.cpp