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, 2016, 2017.
10 // Modifications copyright (c) 2015-2017, Oracle and/or its affiliates.
11
12 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
13 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
14
15 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
16 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
17
18 // Use, modification and distribution is subject to the Boost Software License,
19 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
20 // http://www.boost.org/LICENSE_1_0.txt)
21
22 #include <boost/geometry.hpp>
23 #include <geometry_test_common.hpp>
24
25 namespace bg = boost::geometry;
26
27 //Testing spherical and geographic strategies
28 template <typename CT>
test_spherical_geo()29 void test_spherical_geo()
30 {
31 typedef CT ct;
32
33 //Geographic
34
35 typedef typename bg::model::point
36 <
37 ct, 2, bg::cs::geographic<bg::degree>
38 > pt_geo;
39
40 bg::strategy::area::geographic
41 <
42 bg::strategy::vincenty,
43 5
44 > area_geographic;
45
46 bg::model::polygon<pt_geo> geometry_geo;
47
48 //Spherical
49
50 typedef typename bg::model::point
51 <
52 ct, 2, bg::cs::spherical_equatorial<bg::degree>
53 > pt;
54 bg::model::polygon<pt> geometry;
55
56 // unit-sphere has area of 4-PI. Polygon covering 1/8 of it:
57 // calculations splitted for ttmath
58 std::string poly = "POLYGON((0 0,0 90,90 0,0 0))";
59
60 bg::strategy::area::spherical
61 <
62 ct
63 > strategy_unary(1.0);
64
65 ct const four = 4.0;
66 ct const eight = 8.0;
67 ct expected = four * boost::geometry::math::pi<ct>() / eight;
68 bg::read_wkt(poly, geometry);
69 ct area = bg::area(geometry, strategy_unary);
70 BOOST_CHECK_CLOSE(area, expected, 0.0001);
71
72 // With strategy, radius 2 -> 4 pi r^2
73 bg::strategy::area::spherical
74 <
75 ct
76 > strategy(2.0);
77
78 area = bg::area(geometry, strategy);
79 ct const two = 2.0;
80 BOOST_CHECK_CLOSE(area, two * two * expected, 0.0001);
81
82 // Geographic total area of earth is about 510065626583900.6 (WGS84 ellipsoid)
83 // (510072000 in https://en.wikipedia.org/wiki/Earth)
84 // So the 1/8 is 6.375820332*10^13 and here we get something close to it
85 bg::read_wkt(poly, geometry_geo);
86 area = bg::area(geometry_geo, area_geographic);
87 //GeoGraphicLib gives: 63758202715511.055
88 BOOST_CHECK_CLOSE(area, 63758202715509.844, 0.0001);
89
90
91 // Wrangel Island (dateline crossing)
92 // With (spherical) Earth strategy
93 poly = "POLYGON((-178.7858 70.7852, 177.4758 71.2333, 179.7436 71.5733, -178.7858 70.7852))";
94
95 bg::strategy::area::spherical
96 <
97 ct
98 > spherical_earth(6373);
99 bg::read_wkt(poly, geometry);
100 area = bg::area(geometry, spherical_earth);
101 // SQL Server gives: 4537.9654419375
102 // PostGIS gives: 4537.9311668307
103 // Note: those are Geographic, this test is Spherical
104 BOOST_CHECK_CLOSE(area, 4506.6389, 0.001);
105
106 bg::read_wkt(poly, geometry_geo);
107 area = bg::area(geometry_geo, area_geographic);
108 BOOST_CHECK_CLOSE(area, 4537929936.5349159, 0.0001);
109
110 // Wrangel, more in detail
111 poly = "POLYGON((-178.568604 71.564148,-178.017548 71.449692,-177.833313 71.3461,\
112 -177.502838 71.277466 ,-177.439453 71.226929,-177.620026 71.116638,\
113 -177.9389 71.037491,-178.8186 70.979965,-179.274445 70.907761,\
114 -180 70.9972,179.678314 70.895538,179.272766 70.888596,\
115 178.791016 70.7964,178.617737 71.035538,178.872192 71.217484,\
116 179.530273 71.4383 ,-180 71.535843 ,-179.628601 71.577194,\
117 -179.305298 71.551361,-179.03421 71.597748,-178.568604 71.564148))";
118 bg::read_wkt(poly, geometry);
119 area = bg::area(geometry, spherical_earth);
120 // SQL Server gives: 7669.10402181435
121 // PostGIS gives: 7669.55565459832
122 BOOST_CHECK_CLOSE(area, 7616.523769, 0.001);
123
124 bg::read_wkt(poly, geometry_geo);
125 area = bg::area(geometry_geo, area_geographic);
126 BOOST_CHECK_CLOSE(area, 7669498457.4130802, 0.0001);
127
128 // Check more at the equator
129 /*
130 select 1,geography::STGeomFromText('POLYGON((-178.7858 10.7852 , 179.7436 11.5733 , \
131 177.4758 11.2333 , -178.7858 10.7852))',4326) .STArea()/1000000.0
132 union select 2,geography::STGeomFromText('POLYGON((-178.7858 20.7852 , 179.7436 21.5733 ,\
133 177.4758 21.2333 , -178.7858 20.7852))',4326) .STArea()/1000000.0
134 union select 3,geography::STGeomFromText('POLYGON((-178.7858 30.7852 , 179.7436 31.5733 , \
135 177.4758 31.2333 , -178.7858 30.7852))',4326) .STArea()/1000000.0
136 union select 0,geography::STGeomFromText('POLYGON((-178.7858 0.7852 , 179.7436 1.5733 , \
137 177.4758 1.2333 , -178.7858 0.7852))',4326) .STArea()/1000000.0
138 union select 4,geography::STGeomFromText('POLYGON((-178.7858 40.7852 , 179.7436 41.5733 , \
139 177.4758 41.2333 , -178.7858 40.7852))',4326) .STArea()/1000000.0
140 */
141
142 poly = "POLYGON((-178.7858 0.7852, 177.4758 1.2333, 179.7436 1.5733, -178.7858 0.7852))";
143 bg::read_wkt(poly, geometry);
144 area = bg::area(geometry, spherical_earth);
145 BOOST_CHECK_CLOSE(area, 14136.09946, 0.001); // SQL Server gives: 14064.1902284513
146 bg::read_wkt(poly, geometry_geo);
147 area = bg::area(geometry_geo, area_geographic);
148 BOOST_CHECK_CLOSE(area, 14064129044.674297, 0.0001);
149
150 poly = "POLYGON((-178.7858 10.7852, 177.4758 11.2333, 179.7436 11.5733, -178.7858 10.7852))";
151 bg::read_wkt(poly, geometry);
152 area = bg::area(geometry, spherical_earth);
153 BOOST_CHECK_CLOSE(area, 13760.2456, 0.001); // SQL Server gives: 13697.0941155193
154 bg::read_wkt(poly, geometry_geo);
155 area = bg::area(geometry_geo, area_geographic);
156 BOOST_CHECK_CLOSE(area, 13696308940.315653, 0.0001);
157
158 poly = "POLYGON((-178.7858 20.7852, 177.4758 21.2333, 179.7436 21.5733, -178.7858 20.7852))";
159 bg::read_wkt(poly, geometry);
160 area = bg::area(geometry, spherical_earth);
161 BOOST_CHECK_CLOSE(area, 12987.8682, 0.001); // SQL Server gives: 12944.3970990317 -> -39m^2
162 bg::read_wkt(poly, geometry_geo);
163 area = bg::area(geometry_geo, area_geographic);
164 BOOST_CHECK_CLOSE(area, 12943176284.560806, 0.0001);
165
166 poly = "POLYGON((-178.7858 30.7852, 177.4758 31.2333, 179.7436 31.5733, -178.7858 30.7852))";
167 bg::read_wkt(poly, geometry);
168 area = bg::area(geometry, spherical_earth);
169 BOOST_CHECK_CLOSE(area, 11856.3935, 0.001); // SQL Server gives: 11838.5338423574 -> -18m^2
170 bg::read_wkt(poly, geometry_geo);
171 area = bg::area(geometry_geo, area_geographic);
172 BOOST_CHECK_CLOSE(area, 11837280445.349375, 0.0001);
173
174 poly = "POLYGON((-178.7858 40.7852, 177.4758 41.2333, 179.7436 41.5733, -178.7858 40.7852))";
175 bg::read_wkt(poly, geometry);
176 area = bg::area(geometry, spherical_earth);
177 BOOST_CHECK_CLOSE(area, 10404.627685523914, 0.001);
178 // SQL Server gives: 10412.0607137119, -> +8m^2
179 bg::read_wkt(poly, geometry_geo);
180 area = bg::area(geometry_geo, area_geographic);
181 BOOST_CHECK_CLOSE(area, 10411098789.39222, 0.0001);
182
183 // Concave
184 poly = "POLYGON((0 40,1 42,0 44,2 43,4 44,3 42,4 40,2 41,0 40))";
185 bg::read_wkt(poly, geometry);
186 area = bg::area(geometry, spherical_earth);
187 BOOST_CHECK_CLOSE(area, 73538.2958, 0.001); // SQL Server gives: 73604.2047689719
188 bg::read_wkt(poly, geometry_geo);
189 area = bg::area(geometry_geo, area_geographic);
190 BOOST_CHECK_CLOSE(area, 73604208172.719223, 0.0001);
191
192 // With hole POLYGON((0 40,4 40,4 44,0 44,0 40),(1 41,2 43,3 42,1 41))
193 poly = "POLYGON((0 40,0 44,4 44,4 40,0 40),(1 41,3 42,2 43,1 41))";
194 bg::read_wkt(poly, geometry);
195 area = bg::area(geometry, spherical_earth);
196 BOOST_CHECK_CLOSE(area, 133233.844876, 0.001); // SQL Server gives: 133353.335
197 bg::read_wkt(poly, geometry_geo);
198 area = bg::area(geometry_geo, area_geographic);
199 BOOST_CHECK_CLOSE(area, 133353077343.10347, 0.0001);
200
201 // mean Earth's radius^2
202 double r2 = bg::math::sqr(bg::get_radius<0>(bg::srs::sphere<double>()));
203
204 // around 0 meridian
205 {
206 std::string poly1 = "POLYGON((-10 0,-10 10,0 10,0 0,-10 0))";
207 std::string poly2 = "POLYGON((0 0,0 10,10 10,10 0,0 0))";
208 std::string poly3 = "POLYGON((-5 0,-5 10,5 10,5 0,-5 0))";
209 bg::read_wkt(poly1, geometry);
210 ct area1 = bg::area(geometry);
211 bg::read_wkt(poly2, geometry);
212 ct area2 = bg::area(geometry);
213 bg::read_wkt(poly3, geometry);
214 ct area3 = bg::area(geometry);
215 BOOST_CHECK_CLOSE(area1, area2, 0.001);
216 BOOST_CHECK_CLOSE(area2, area3, 0.001);
217 BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1848, 0.001);
218 //geographic
219 bg::read_wkt(poly1, geometry_geo);
220 area1 = bg::area(geometry_geo, area_geographic);
221 bg::read_wkt(poly2, geometry_geo);
222 area2 = bg::area(geometry_geo, area_geographic);
223 bg::read_wkt(poly3, geometry_geo);
224 area3 = bg::area(geometry_geo, area_geographic);
225 BOOST_CHECK_CLOSE(area1, area2, 0.001);
226 BOOST_CHECK_CLOSE(area2, area3, 0.001);
227 BOOST_CHECK_CLOSE(area1, 1227877191611.2805, 0.001);
228 }
229 {
230 std::string poly1 = "POLYGON((-10 -5,-10 5,0 5,0 -5,-10 -5))";
231 std::string poly2 = "POLYGON((0 -5,0 5,10 5,10 -5,0 -5))";
232 std::string poly3 = "POLYGON((-5 -5,-5 5,5 5,5 -5,-5 -5))";
233 bg::read_wkt(poly1, geometry);
234 ct area1 = bg::area(geometry);
235 bg::read_wkt(poly2, geometry);
236 ct area2 = bg::area(geometry);
237 bg::read_wkt(poly3, geometry);
238 ct area3 = bg::area(geometry);
239 BOOST_CHECK_CLOSE(area1, area2, 0.001);
240 BOOST_CHECK_CLOSE(area2, area3, 0.001);
241 BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0261, 0.001);
242 //geographic
243 bg::read_wkt(poly1, geometry_geo);
244 area1 = bg::area(geometry_geo, area_geographic);
245 bg::read_wkt(poly2, geometry_geo);
246 area2 = bg::area(geometry_geo, area_geographic);
247 bg::read_wkt(poly3, geometry_geo);
248 area3 = bg::area(geometry_geo, area_geographic);
249 BOOST_CHECK_CLOSE(area1, area2, 0.001);
250 BOOST_CHECK_CLOSE(area2, area3, 0.001);
251 BOOST_CHECK_CLOSE(area1, 1232514639151.6477, 0.001);
252 }
253 // around 180 meridian
254 {
255 std::string poly1 = "POLYGON((-180 0,-180 10,-170 10,-170 0,-180 0))";
256 std::string poly2 = "POLYGON((175 0,175 10,-175 10,-175 0,175 0))";
257 std::string poly3 = "POLYGON((170 0,170 10,180 10,180 0,170 0))";
258 std::string poly4 = "POLYGON((170 0,170 10,-180 10,-180 0,170 0))";
259 std::string poly5 = "POLYGON((180 0,180 10,-170 10,-170 0,180 0))";
260 bg::read_wkt(poly1, geometry);
261 ct area1 = bg::area(geometry);
262 bg::read_wkt(poly2, geometry);
263 ct area2 = bg::area(geometry);
264 bg::read_wkt(poly3, geometry);
265 ct area3 = bg::area(geometry);
266 bg::read_wkt(poly4, geometry);
267 ct area4 = bg::area(geometry);
268 bg::read_wkt(poly5, geometry);
269 ct area5 = bg::area(geometry);
270 BOOST_CHECK_CLOSE(area1, area2, 0.001);
271 BOOST_CHECK_CLOSE(area2, area3, 0.001);
272 BOOST_CHECK_CLOSE(area3, area4, 0.001);
273 BOOST_CHECK_CLOSE(area4, area5, 0.001);
274 BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1833, 0.001);
275 //geographic
276 bg::read_wkt(poly1, geometry_geo);
277 area1 = bg::area(geometry_geo, area_geographic);
278 bg::read_wkt(poly2, geometry_geo);
279 area2 = bg::area(geometry_geo, area_geographic);
280 bg::read_wkt(poly3, geometry_geo);
281 area3 = bg::area(geometry_geo, area_geographic);
282 bg::read_wkt(poly4, geometry_geo);
283 area4 = bg::area(geometry_geo, area_geographic);
284 bg::read_wkt(poly5, geometry_geo);
285 area5 = bg::area(geometry_geo, area_geographic);
286 BOOST_CHECK_CLOSE(area1, area2, 0.001);
287 BOOST_CHECK_CLOSE(area2, area3, 0.001);
288 BOOST_CHECK_CLOSE(area3, area4, 0.001);
289 BOOST_CHECK_CLOSE(area4, area5, 0.001);
290 BOOST_CHECK_CLOSE(area1, 1227877191611.2805, 0.001);
291 }
292 {
293 std::string poly1 = "POLYGON((-180 -5,-180 5,-170 5,-170 -5,-180 -5))";
294 std::string poly2 = "POLYGON((175 -5,175 5,-175 5,-175 -5,175 -5))";
295 std::string poly3 = "POLYGON((170 -5,170 5,180 5,180 -5,170 -5))";
296 std::string poly4 = "POLYGON((170 -5,170 5,180 5,180 -5,170 -5))";
297 std::string poly5 = "POLYGON((180 -5,180 5,-170 5,-170 -5,180 -5))";
298 bg::read_wkt(poly1, geometry);
299 ct area1 = bg::area(geometry);
300 bg::read_wkt(poly2, geometry);
301 ct area2 = bg::area(geometry);
302 bg::read_wkt(poly3, geometry);
303 ct area3 = bg::area(geometry);
304 bg::read_wkt(poly4, geometry);
305 ct area4 = bg::area(geometry);
306 bg::read_wkt(poly5, geometry);
307 ct area5 = bg::area(geometry);
308 BOOST_CHECK_CLOSE(area1, area2, 0.001);
309 BOOST_CHECK_CLOSE(area2, area3, 0.001);
310 BOOST_CHECK_CLOSE(area3, area4, 0.001);
311 BOOST_CHECK_CLOSE(area4, area5, 0.001);
312 BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0247, 0.001);
313 //geographic
314 bg::read_wkt(poly1, geometry_geo);
315 area1 = bg::area(geometry_geo, area_geographic);
316 bg::read_wkt(poly2, geometry_geo);
317 area2 = bg::area(geometry_geo, area_geographic);
318 bg::read_wkt(poly3, geometry_geo);
319 area3 = bg::area(geometry_geo, area_geographic);
320 bg::read_wkt(poly4, geometry_geo);
321 area4 = bg::area(geometry_geo, area_geographic);
322 bg::read_wkt(poly5, geometry_geo);
323 area5 = bg::area(geometry_geo, area_geographic);
324 BOOST_CHECK_CLOSE(area1, area2, 0.001);
325 BOOST_CHECK_CLOSE(area2, area3, 0.001);
326 BOOST_CHECK_CLOSE(area3, area4, 0.001);
327 BOOST_CHECK_CLOSE(area4, area5, 0.001);
328 BOOST_CHECK_CLOSE(area1, 1232514639151.6477, 0.001);
329 }
330 // around poles
331 {
332 std::string poly1 = "POLYGON((0 80,-90 80,-180 80,90 80,0 80))";
333 std::string poly2 = "POLYGON((0 80,-90 80,180 80,90 80,0 80))";
334 std::string poly3 = "POLYGON((0 -80,90 -80,-180 -80,-90 -80,0 -80))";
335 std::string poly4 = "POLYGON((0 -80,90 -80,180 -80,-90 -80,0 -80))";
336 bg::read_wkt(poly1, geometry);
337 ct area1 = bg::area(geometry);
338 bg::read_wkt(poly2, geometry);
339 ct area2 = bg::area(geometry);
340 bg::read_wkt(poly3, geometry);
341 ct area3 = bg::area(geometry);
342 bg::read_wkt(poly4, geometry);
343 ct area4 = bg::area(geometry);
344 BOOST_CHECK_CLOSE(area1, area2, 0.001);
345 BOOST_CHECK_CLOSE(area2, area3, 0.001);
346 BOOST_CHECK_CLOSE(area3, area4, 0.001);
347 //geographic
348 bg::read_wkt(poly1, geometry_geo);
349 area1 = bg::area(geometry_geo, area_geographic);
350 bg::read_wkt(poly2, geometry_geo);
351 area2 = bg::area(geometry_geo, area_geographic);
352 bg::read_wkt(poly3, geometry_geo);
353 area3 = bg::area(geometry_geo, area_geographic);
354 bg::read_wkt(poly4, geometry_geo);
355 area4 = bg::area(geometry_geo, area_geographic);
356 BOOST_CHECK_CLOSE(area1, area2, 0.001);
357 BOOST_CHECK_CLOSE(area2, area3, 0.001);
358 BOOST_CHECK_CLOSE(area3, area4, 0.001);
359 }
360
361 {
362 bg::model::ring<pt> aurha; // a'dam-utr-rott.-den haag-a'dam
363 std::string poly = "POLYGON((4.892 52.373,5.119 52.093,4.479 51.930,\
364 4.23 52.08,4.892 52.373))";
365 bg::read_wkt(poly, aurha);
366 /*if (polar)
367 {
368 // Create colatitudes (measured from pole)
369 BOOST_FOREACH(pt& p, aurha)
370 {
371 bg::set<1>(p, ct(90) - bg::get<1>(p));
372 }
373 bg::correct(aurha);
374 }*/
375 bg::strategy::area::spherical
376 <
377 ct
378 > area_spherical(6372.795);
379 area = bg::area(aurha, area_spherical);
380 BOOST_CHECK_CLOSE(area, 1476.645675, 0.0001);
381 //geographic
382 bg::read_wkt(poly, geometry_geo);
383 area = bg::area(geometry_geo, area_geographic);
384 BOOST_CHECK_CLOSE(area, 1481555970.0765088, 0.001);
385
386 // SQL Server gives: 1481.55595960659
387 // for select geography::STGeomFromText('POLYGON((4.892 52.373,4.23 52.08,
388 // 4.479 51.930,5.119 52.093,4.892 52.373))',4326).STArea()/1000000.0
389 }
390
391 {
392 bg::model::polygon<pt, false> geometry_sph;
393 std::string wkt = "POLYGON((0 0, 5 0, 5 5, 0 5, 0 0))";
394 bg::read_wkt(wkt, geometry_sph);
395
396 area = bg::area(geometry_sph, bg::strategy::area::spherical<>(6371228.0));
397 BOOST_CHECK_CLOSE(area, 308932296103.83051, 0.0001);
398
399 bg::model::polygon<pt_geo, false> geometry_geo;
400 bg::read_wkt(wkt, geometry_geo);
401
402 area = bg::area(geometry_geo, bg::strategy::area::geographic<>(bg::srs::spheroid<double>(6371228.0, 6371228.0)));
403 BOOST_CHECK_CLOSE(area, 308932296103.82574, 0.001);
404 }
405 }
406
test_main(int,char * [])407 int test_main(int, char* [])
408 {
409
410 test_spherical_geo<double>();
411
412 return 0;
413 }
414