ceph/src/boost/libs/geometry/test/algorithms/area/area_sph_geo.cpp

   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-2021.
  10 // Modifications copyright (c) 2015-2021, 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
  29 template <typename polygon>
  30 auto 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
  51 namespace bg = boost::geometry;
  52
  53 //Testing spherical and geographic strategies
  54 template <typename CT>
  55 void test_spherical_geo()
  56 {
  57     typedef CT ct;
  58
  59     //Geographic
  60
  61     typedef typename bg::model::point
  62         <
  63             ct, 2, bg::cs::geographic<bg::degree>
  64         > pt_geo;
  65
  66     bg::strategy::area::geographic
  67         <
  68             bg::strategy::vincenty,
  69             5
  70         > area_geographic;
  71
  72     bg::model::polygon<pt_geo> geometry_geo;
  73
  74     //Spherical
  75
  76     typedef typename bg::model::point
  77         <
  78             ct, 2, bg::cs::spherical_equatorial<bg::degree>
  79         > pt;
  80     bg::model::polygon<pt> geometry;
  81
  82     // unit-sphere has area of 4-PI. Polygon covering 1/8 of it:
  83     std::string poly = "POLYGON((0 0,0 90,90 0,0 0))";
  84
  85     bg::strategy::area::spherical
  86         <
  87             ct
  88         > strategy_unary(1.0);
  89
  90     ct const four = 4.0;
  91     ct const eight = 8.0;
  92     ct expected = four * boost::geometry::math::pi<ct>() / eight;
  93     bg::read_wkt(poly, geometry);
  94     ct area = bg::area(geometry, strategy_unary);
  95     BOOST_CHECK_CLOSE(area, expected, 0.0001);
  96
  97     // With strategy, radius 2 -> 4 pi r^2
  98     bg::strategy::area::spherical
  99         <
 100             ct
 101         > strategy(2.0);
 102
 103     area = bg::area(geometry, strategy);
 104     ct const two = 2.0;
 105     BOOST_CHECK_CLOSE(area, two * two * expected, 0.0001);
 106
 107     // Geographic total area of earth is about 510065626583900.6 (WGS84 ellipsoid)
 108     // (510072000 in https://en.wikipedia.org/wiki/Earth)
 109     // So the 1/8 is 6.375820332*10^13 and here we get something close to it
 110     bg::read_wkt(poly, geometry_geo);
 111     area = bg::area(geometry_geo, area_geographic);
 112     //GeoGraphicLib gives: 63758202715511.055
 113     BOOST_CHECK_CLOSE(area, 63758202715509.844, 0.0001);
 114
 115
 116     // Wrangel Island (dateline crossing)
 117     // With (spherical) Earth strategy
 118     poly = "POLYGON((-178.7858 70.7852, 177.4758 71.2333, 179.7436 71.5733, -178.7858 70.7852))";
 119
 120     bg::strategy::area::spherical
 121         <
 122             ct
 123         > spherical_earth(6373);
 124     bg::read_wkt(poly, geometry);
 125     area = bg::area(geometry, spherical_earth);
 126     // SQL Server gives: 4537.9654419375
 127     // PostGIS gives: 4537.9311668307
 128     // Note: those are Geographic, this test is Spherical
 129     BOOST_CHECK_CLOSE(area, 4506.6389, 0.001);
 130
 131     bg::read_wkt(poly, geometry_geo);
 132     area = bg::area(geometry_geo, area_geographic);
 133     BOOST_CHECK_CLOSE(area, 4537929936.5883484, 0.0001);
 134 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 135     BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.1);
 136 #endif
 137     // Wrangel, more in detail
 138     poly = "POLYGON((-178.568604 71.564148,-178.017548 71.449692,-177.833313 71.3461,\
 139             -177.502838 71.277466 ,-177.439453 71.226929,-177.620026 71.116638,\
 140             -177.9389 71.037491,-178.8186 70.979965,-179.274445 70.907761,\
 141             -180 70.9972,179.678314 70.895538,179.272766 70.888596,\
 142              178.791016 70.7964,178.617737 71.035538,178.872192 71.217484,\
 143              179.530273 71.4383 ,-180 71.535843 ,-179.628601 71.577194,\
 144             -179.305298 71.551361,-179.03421 71.597748,-178.568604 71.564148))";
 145     bg::read_wkt(poly, geometry);
 146     area = bg::area(geometry, spherical_earth);
 147     // SQL Server gives: 7669.10402181435
 148     // PostGIS gives: 7669.55565459832
 149     BOOST_CHECK_CLOSE(area, 7616.523769, 0.001);
 150
 151     bg::read_wkt(poly, geometry_geo);
 152     area = bg::area(geometry_geo, area_geographic);
 153     BOOST_CHECK_CLOSE(area, 7669498457.4130802, 0.0001);
 154
 155     // Check more at the equator
 156     /*
 157     select 1,geography::STGeomFromText('POLYGON((-178.7858 10.7852 , 179.7436 11.5733 , \
 158                         177.4758 11.2333 , -178.7858 10.7852))',4326) .STArea()/1000000.0
 159     union select 2,geography::STGeomFromText('POLYGON((-178.7858 20.7852 , 179.7436 21.5733 ,\
 160                         177.4758 21.2333 , -178.7858 20.7852))',4326) .STArea()/1000000.0
 161     union select 3,geography::STGeomFromText('POLYGON((-178.7858 30.7852 , 179.7436 31.5733 , \
 162                         177.4758 31.2333 , -178.7858 30.7852))',4326) .STArea()/1000000.0
 163     union select 0,geography::STGeomFromText('POLYGON((-178.7858 0.7852 , 179.7436 1.5733 , \
 164                         177.4758 1.2333 , -178.7858 0.7852))',4326) .STArea()/1000000.0
 165     union select 4,geography::STGeomFromText('POLYGON((-178.7858 40.7852 , 179.7436 41.5733 , \
 166                         177.4758 41.2333 , -178.7858 40.7852))',4326) .STArea()/1000000.0
 167     */
 168
 169     poly = "POLYGON((-178.7858 0.7852, 177.4758 1.2333, 179.7436 1.5733, -178.7858 0.7852))";
 170     bg::read_wkt(poly, geometry);
 171     area = bg::area(geometry, spherical_earth);
 172     BOOST_CHECK_CLOSE(area, 14136.09946, 0.001); // SQL Server gives: 14064.1902284513
 173     bg::read_wkt(poly, geometry_geo);
 174     area = bg::area(geometry_geo, area_geographic);
 175     BOOST_CHECK_CLOSE(area, 14064129044.677765, 0.0001);
 176
 177     poly = "POLYGON((-178.7858 10.7852, 177.4758 11.2333, 179.7436 11.5733, -178.7858 10.7852))";
 178     // Geographiclib 1.51 :13696308940.35794067382812
 179     bg::read_wkt(poly, geometry);
 180     area = bg::area(geometry, spherical_earth);
 181     BOOST_CHECK_CLOSE(area, 13760.2456, 0.001); // SQL Server gives: 13697.0941155193
 182     bg::read_wkt(poly, geometry_geo);
 183     area = bg::area(geometry_geo, area_geographic);
 184     BOOST_CHECK_CLOSE(area, 13696308940.342197, 0.0001);
 185
 186     poly = "POLYGON((-178.7858 20.7852, 177.4758 21.2333, 179.7436 21.5733, -178.7858 20.7852))";
 187     bg::read_wkt(poly, geometry);
 188     area = bg::area(geometry, spherical_earth);
 189     BOOST_CHECK_CLOSE(area, 12987.8682, 0.001); // SQL Server gives: 12944.3970990317 -> -39m^2
 190     bg::read_wkt(poly, geometry_geo);
 191     area = bg::area(geometry_geo, area_geographic);
 192     BOOST_CHECK_CLOSE(area, 12943176284.489822, 0.0001);
 193 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 194     BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.00001);
 195 #endif
 196
 197     poly = "POLYGON((-178.7858 30.7852, 177.4758 31.2333, 179.7436 31.5733, -178.7858 30.7852))";
 198     bg::read_wkt(poly, geometry);
 199     area = bg::area(geometry, spherical_earth);
 200     BOOST_CHECK_CLOSE(area, 11856.3935, 0.001); // SQL Server gives: 11838.5338423574 -> -18m^2
 201     bg::read_wkt(poly, geometry_geo);
 202     area = bg::area(geometry_geo, area_geographic);
 203     BOOST_CHECK_CLOSE(area, 11837280445.394035, 0.0001);
 204 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 205     // GEOGRAPHICLIB 1.51: 11837280466.204895
 206     BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.000001);
 207 #endif
 208
 209     poly = "POLYGON((-178.7858 40.7852, 177.4758 41.2333, 179.7436 41.5733, -178.7858 40.7852))";
 210     bg::read_wkt(poly, geometry);
 211     area = bg::area(geometry, spherical_earth);
 212     BOOST_CHECK_CLOSE(area, 10404.627685523914, 0.001);
 213     // SQL Server gives: 10412.0607137119, -> +8m^2
 214     bg::read_wkt(poly, geometry_geo);
 215     area = bg::area(geometry_geo, area_geographic);
 216     BOOST_CHECK_CLOSE(area, 10411098789.387386, 0.0001);
 217 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 218     // GEOGRAPHICLIB 1.51: 10411098790.577026
 219     BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.000001);
 220 #endif
 221
 222     // Concave
 223     poly = "POLYGON((0 40,1 42,0 44,2 43,4 44,3 42,4 40,2 41,0 40))";
 224     bg::read_wkt(poly, geometry);
 225     area = bg::area(geometry, spherical_earth);
 226     BOOST_CHECK_CLOSE(area, 73538.2958, 0.001); // SQL Server gives: 73604.2047689719
 227     bg::read_wkt(poly, geometry_geo);
 228     area = bg::area(geometry_geo, area_geographic);
 229     BOOST_CHECK_CLOSE(area, 73604163095.545319, 0.0001);
 230 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 231     // GEOGRAPHICLIB 1.51: 73604163091.274048
 232     BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.000001);
 233 #endif
 234
 235     // With hole POLYGON((0 40,4 40,4 44,0 44,0 40),(1 41,2 43,3 42,1 41))
 236     poly = "POLYGON((0 40,0 44,4 44,4 40,0 40),(1 41,3 42,2 43,1 41))";
 237     bg::read_wkt(poly, geometry);
 238     area = bg::area(geometry, spherical_earth);
 239     BOOST_CHECK_CLOSE(area, 133233.844876, 0.001); // SQL Server gives: 133353.335
 240     bg::read_wkt(poly, geometry_geo);
 241     area = bg::area(geometry_geo, area_geographic);
 242     BOOST_CHECK_CLOSE(area, 133353077343.26692, 0.0001);
 243 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 244     // GEOGRAPHICLIB 1.51: 147189206350.20142
 245     //BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area, 0.9);
 246 #endif
 247
 248     // mean Earth's radius^2
 249     double r2 = bg::math::sqr(bg::get_radius<0>(bg::srs::sphere<double>()));
 250
 251     // around 0 meridian
 252     {
 253         std::string poly1 = "POLYGON((-10 0,-10 10,0 10,0 0,-10 0))";
 254         std::string poly2 = "POLYGON((0 0,0 10,10 10,10 0,0 0))";
 255         std::string poly3 = "POLYGON((-5 0,-5 10,5 10,5 0,-5 0))";
 256         bg::read_wkt(poly1, geometry);
 257         ct area1 = bg::area(geometry);
 258         bg::read_wkt(poly2, geometry);
 259         ct area2 = bg::area(geometry);
 260         bg::read_wkt(poly3, geometry);
 261         ct area3 = bg::area(geometry);
 262         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 263         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 264         BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1848, 0.001);
 265         //geographic
 266         bg::read_wkt(poly1, geometry_geo);
 267         area1 = bg::area(geometry_geo, area_geographic);
 268         bg::read_wkt(poly2, geometry_geo);
 269         area2 = bg::area(geometry_geo, area_geographic);
 270         bg::read_wkt(poly3, geometry_geo);
 271         area3 = bg::area(geometry_geo, area_geographic);
 272         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 273         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 274         BOOST_CHECK_CLOSE(area1, 1227877191611.3574, 0.001);
 275 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 276         // GEOGRAPHICLIB 1.51: 1227877191609.6267
 277         BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area3, 0.000001);
 278 #endif
 279     }
 280     {
 281         std::string poly1 = "POLYGON((-10 -5,-10 5,0 5,0 -5,-10 -5))";
 282         std::string poly2 = "POLYGON((0 -5,0 5,10 5,10 -5,0 -5))";
 283         std::string poly3 = "POLYGON((-5 -5,-5 5,5 5,5 -5,-5 -5))";
 284         bg::read_wkt(poly1, geometry);
 285         ct area1 = bg::area(geometry);
 286         bg::read_wkt(poly2, geometry);
 287         ct area2 = bg::area(geometry);
 288         bg::read_wkt(poly3, geometry);
 289         ct area3 = bg::area(geometry);
 290         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 291         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 292         BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0261, 0.001);
 293         //geographic
 294         bg::read_wkt(poly1, geometry_geo);
 295         area1 = bg::area(geometry_geo, area_geographic);
 296         bg::read_wkt(poly2, geometry_geo);
 297         area2 = bg::area(geometry_geo, area_geographic);
 298         bg::read_wkt(poly3, geometry_geo);
 299         area3 = bg::area(geometry_geo, area_geographic);
 300         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 301         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 302         BOOST_CHECK_CLOSE(area1, 1232514639151.7168, 0.001);
 303 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 304         // GEOGRAPHICLIB 1.51: 1232514639151.6357
 305         BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area3, 0.000001);
 306 #endif
 307     }
 308     // around 180 meridian
 309     {
 310         std::string poly1 = "POLYGON((-180 0,-180 10,-170 10,-170 0,-180 0))";
 311         std::string poly2 = "POLYGON((175 0,175 10,-175 10,-175 0,175 0))";
 312         std::string poly3 = "POLYGON((170 0,170 10,180 10,180 0,170 0))";
 313         std::string poly4 = "POLYGON((170 0,170 10,-180 10,-180 0,170 0))";
 314         std::string poly5 = "POLYGON((180 0,180 10,-170 10,-170 0,180 0))";
 315         bg::read_wkt(poly1, geometry);
 316         ct area1 = bg::area(geometry);
 317         bg::read_wkt(poly2, geometry);
 318         ct area2 = bg::area(geometry);
 319         bg::read_wkt(poly3, geometry);
 320         ct area3 = bg::area(geometry);
 321         bg::read_wkt(poly4, geometry);
 322         ct area4 = bg::area(geometry);
 323         bg::read_wkt(poly5, geometry);
 324         ct area5 = bg::area(geometry);
 325         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 326         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 327         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 328         BOOST_CHECK_CLOSE(area4, area5, 0.001);
 329         BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1833, 0.001);
 330         //geographic
 331         bg::read_wkt(poly1, geometry_geo);
 332         area1 = bg::area(geometry_geo, area_geographic);
 333         bg::read_wkt(poly2, geometry_geo);
 334         area2 = bg::area(geometry_geo, area_geographic);
 335         bg::read_wkt(poly3, geometry_geo);
 336         area3 = bg::area(geometry_geo, area_geographic);
 337         bg::read_wkt(poly4, geometry_geo);
 338         area4 = bg::area(geometry_geo, area_geographic);
 339         bg::read_wkt(poly5, geometry_geo);
 340         area5 = bg::area(geometry_geo, area_geographic);
 341         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 342         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 343         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 344         BOOST_CHECK_CLOSE(area4, area5, 0.001);
 345         BOOST_CHECK_CLOSE(area1, 1227877191611.3574, 0.001);
 346 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 347         // GEOGRAPHICLIB 1.51: 1227877191609.6267
 348         BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area5, 0.000001);
 349 #endif
 350     }
 351     {
 352         std::string poly1 = "POLYGON((-180 -5,-180 5,-170 5,-170 -5,-180 -5))";
 353         std::string poly2 = "POLYGON((175 -5,175 5,-175 5,-175 -5,175 -5))";
 354         std::string poly3 = "POLYGON((170 -5,170 5,180 5,180 -5,170 -5))";
 355         std::string poly4 = "POLYGON((170 -5,170 5,180 5,180 -5,170 -5))";
 356         std::string poly5 = "POLYGON((180 -5,180 5,-170 5,-170 -5,180 -5))";
 357         bg::read_wkt(poly1, geometry);
 358         ct area1 = bg::area(geometry);
 359         bg::read_wkt(poly2, geometry);
 360         ct area2 = bg::area(geometry);
 361         bg::read_wkt(poly3, geometry);
 362         ct area3 = bg::area(geometry);
 363         bg::read_wkt(poly4, geometry);
 364         ct area4 = bg::area(geometry);
 365         bg::read_wkt(poly5, geometry);
 366         ct area5 = bg::area(geometry);
 367         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 368         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 369         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 370         BOOST_CHECK_CLOSE(area4, area5, 0.001);
 371         BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0247, 0.001);
 372         //geographic
 373         bg::read_wkt(poly1, geometry_geo);
 374         area1 = bg::area(geometry_geo, area_geographic);
 375         bg::read_wkt(poly2, geometry_geo);
 376         area2 = bg::area(geometry_geo, area_geographic);
 377         bg::read_wkt(poly3, geometry_geo);
 378         area3 = bg::area(geometry_geo, area_geographic);
 379         bg::read_wkt(poly4, geometry_geo);
 380         area4 = bg::area(geometry_geo, area_geographic);
 381         bg::read_wkt(poly5, geometry_geo);
 382         area5 = bg::area(geometry_geo, area_geographic);
 383         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 384         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 385         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 386         BOOST_CHECK_CLOSE(area4, area5, 0.001);
 387         BOOST_CHECK_CLOSE(area1, 1232514639151.7168, 0.001);
 388 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
 389         // GEOGRAPHICLIB 1.51: 1232514639151.6357
 390         BOOST_CHECK_CLOSE(geographiclib_area(geometry_geo), area5, 0.000001);
 391 #endif
 392     }
 393     // around poles
 394     {
 395         std::string poly1 = "POLYGON((0 80,-90 80,-180 80,90 80,0 80))";
 396         std::string poly2 = "POLYGON((0 80,-90 80,180 80,90 80,0 80))";
 397         std::string poly3 = "POLYGON((0 -80,90 -80,-180 -80,-90 -80,0 -80))";
 398         std::string poly4 = "POLYGON((0 -80,90 -80,180 -80,-90 -80,0 -80))";
 399         bg::read_wkt(poly1, geometry);
 400         ct area1 = bg::area(geometry);
 401         bg::read_wkt(poly2, geometry);
 402         ct area2 = bg::area(geometry);
 403         bg::read_wkt(poly3, geometry);
 404         ct area3 = bg::area(geometry);
 405         bg::read_wkt(poly4, geometry);
 406         ct area4 = bg::area(geometry);
 407         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 408         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 409         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 410         //geographic
 411         bg::read_wkt(poly1, geometry_geo);
 412         area1 = bg::area(geometry_geo, area_geographic);
 413         bg::read_wkt(poly2, geometry_geo);
 414         area2 = bg::area(geometry_geo, area_geographic);
 415         bg::read_wkt(poly3, geometry_geo);
 416         area3 = bg::area(geometry_geo, area_geographic);
 417         bg::read_wkt(poly4, geometry_geo);
 418         area4 = bg::area(geometry_geo, area_geographic);
 419         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 420         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 421         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 422     }
 423
 424     {
 425         bg::model::ring<pt> aurha; // a'dam-utr-rott.-den haag-a'dam
 426         std::string poly = "POLYGON((4.892 52.373,5.119 52.093,4.479 51.930,\
 427                                      4.23 52.08,4.892 52.373))";
 428         bg::read_wkt(poly, aurha);
 429         /*if (polar)
 430         {
 431             // Create colatitudes (measured from pole)
 432             for (pt& p : aurha)
 433             {
 434                 bg::set<1>(p, ct(90) - bg::get<1>(p));
 435             }
 436             bg::correct(aurha);
 437         }*/
 438         bg::strategy::area::spherical
 439             <
 440                 ct
 441             > area_spherical(6372.795);
 442         area = bg::area(aurha, area_spherical);
 443         BOOST_CHECK_CLOSE(area, 1476.645675, 0.0001);
 444         //geographic
 445         bg::read_wkt(poly, geometry_geo);
 446         area = bg::area(geometry_geo, area_geographic);
 447         BOOST_CHECK_CLOSE(area, 1481555970.0765088, 0.001);
 448
 449         // SQL Server gives: 1481.55595960659
 450         // for select geography::STGeomFromText('POLYGON((4.892 52.373,4.23 52.08,
 451         //      4.479 51.930,5.119 52.093,4.892 52.373))',4326).STArea()/1000000.0
 452     }
 453
 454     {
 455         bg::model::polygon<pt, false> geometry_sph;
 456         std::string wkt = "POLYGON((0 0, 5 0, 5 5, 0 5, 0 0))";
 457         bg::read_wkt(wkt, geometry_sph);
 458
 459         area = bg::area(geometry_sph, bg::strategy::area::spherical<>(6371228.0));
 460         BOOST_CHECK_CLOSE(area, 308932296103.83051, 0.0001);
 461
 462         bg::model::polygon<pt_geo, false> geometry_geo;
 463         bg::read_wkt(wkt, geometry_geo);
 464
 465         area = bg::area(geometry_geo, bg::strategy::area::geographic<>(bg::srs::spheroid<double>(6371228.0, 6371228.0)));
 466         BOOST_CHECK_CLOSE(area, 308932296103.82574, 0.001);
 467     }
 468
 469     {
 470         // small areas: ~20m^2
 471         // see https://github.com/boostorg/geometry/issues/799
 472         // geographiclib 1.50.1 returns: 25.5736
 473
 474         bg::model::polygon<pt> geometry;
 475         std::string wkt ="POLYGON((-0.020000 51.470027,-0.020031 51.470019,-0.020043 51.470000,\
 476                  -0.020031 51.469981,-0.020000 51.469973,-0.019969 51.469981,\
 477                  -0.019957 51.470000,-0.019969 51.470019,-0.020000 51.470027))";
 478         bg::read_wkt(wkt, geometry);
 479
 480         bg::strategy::area::spherical
 481             <
 482                 ct
 483             > area_spherical(6371228.0);
 484         area = bg::area(geometry, area_spherical);
 485         BOOST_CHECK_CLOSE(area, -25.48014, 0.0001);
 486
 487         //geographic
 488
 489         using pt_geo_d = bg::model::point
 490             <
 491                 double, 2, bg::cs::geographic<bg::degree>
 492             >;
 493         using pt_geo_ld = bg::model::point
 494             <
 495                 long double, 2, bg::cs::geographic<bg::degree>
 496             >;
 497
 498         bg::strategy::area::geographic<bg::strategy::andoyer> area_a;
 499         bg::strategy::area::geographic<bg::strategy::thomas> area_t;
 500         bg::strategy::area::geographic<bg::strategy::vincenty> area_v;
 501         bg::strategy::area::geographic<bg::strategy::karney> area_k;
 502
 503         bg::model::polygon<pt_geo_d> geometry_geo_d;
 504
 505         bg::read_wkt(wkt, geometry_geo_d);
 506         area = bg::area(geometry_geo_d, area_a);
 507         BOOST_CHECK_CLOSE(area, -25.47837, 0.001);
 508         area = bg::area(geometry_geo_d, area_t);
 509         BOOST_CHECK_CLOSE(area, -25.57355, 0.001);
 510         area = bg::area(geometry_geo_d, area_v);
 511         BOOST_CHECK_CLOSE(area, -25.55881, 0.001);
 512         area = bg::area(geometry_geo_d, area_k);
 513         BOOST_CHECK_CLOSE(area, -25.57357, 0.001);
 514
 515         bg::model::polygon<pt_geo_ld> geometry_geo_ld;
 516
 517         bg::read_wkt(wkt, geometry_geo_ld);
 518         area = bg::area(geometry_geo_ld, area_a);
 519         BOOST_CHECK_CLOSE(area, -25.57978, 0.4); // -25.478374 with vc-14.1
 520         area = bg::area(geometry_geo_ld, area_t);
 521         BOOST_CHECK_CLOSE(area, -25.57359, 0.001);
 522         area = bg::area(geometry_geo_ld, area_v);
 523         BOOST_CHECK_CLOSE(area, -25.57394, 0.06); // -25.558816 with vc-14.1
 524         area = bg::area(geometry_geo_ld, area_k);
 525         BOOST_CHECK_CLOSE(area, -25.57359, 0.001);
 526     }
 527
 528     {
 529         // very small areas: ~2m^2
 530         // geographiclib 1.50.1 returns: 2.24581
 531
 532         bg::model::polygon<pt> geometry;
 533         std::string wkt ="POLYGON((0.020000 51.470027,0.020005 51.470027,\
 534             0.020005 51.470000,0.020007 51.4699,0.020000 51.470027))";
 535         bg::read_wkt(wkt, geometry);
 536
 537         bg::strategy::area::spherical
 538             <
 539                 ct
 540             > area_spherical(6371228.0);
 541         area = bg::area(geometry, area_spherical);
 542         BOOST_CHECK_CLOSE(area, 2.2375981058307093, 0.001);
 543
 544         //geographic
 545
 546         using pt_geo_d = bg::model::point
 547             <
 548                 double, 2, bg::cs::geographic<bg::degree>
 549             >;
 550         using pt_geo_ld = bg::model::point
 551             <
 552                 long double, 2, bg::cs::geographic<bg::degree>
 553             >;
 554
 555         bg::strategy::area::geographic<bg::strategy::andoyer> area_a;
 556         bg::strategy::area::geographic<bg::strategy::thomas> area_t;
 557         bg::strategy::area::geographic<bg::strategy::vincenty> area_v;
 558         bg::strategy::area::geographic<bg::strategy::karney> area_k;
 559
 560         bg::model::polygon<pt_geo_d> geometry_geo_d;
 561
 562         bg::read_wkt(wkt, geometry_geo_d);
 563         area = bg::area(geometry_geo_d, area_a);
 564         BOOST_CHECK_CLOSE(area, 2.2374430036130444, 0.001);
 565         area = bg::area(geometry_geo_d, area_t);
 566         BOOST_CHECK_CLOSE(area, 2.245814319435655, 0.001);
 567         area = bg::area(geometry_geo_d, area_v);
 568         BOOST_CHECK_CLOSE(area, 2.2374430036130444, 0.001);
 569         area = bg::area(geometry_geo_d, area_k);
 570         BOOST_CHECK_CLOSE(area, 2.2458140167194878, 0.001);
 571
 572         bg::model::polygon<pt_geo_ld> geometry_geo_ld;
 573
 574         bg::read_wkt(wkt, geometry_geo_ld);
 575         area = bg::area(geometry_geo_ld, area_a);
 576         BOOST_CHECK_CLOSE(area, 2.2374430039839437, 0.001);
 577         area = bg::area(geometry_geo_ld, area_t);
 578         BOOST_CHECK_CLOSE(area, 2.2457934740573235, 0.001);
 579         area = bg::area(geometry_geo_ld, area_v);
 580         BOOST_CHECK_CLOSE(area, 2.2374430039839437, 0.001);
 581         area = bg::area(geometry_geo_ld, area_k);
 582         BOOST_CHECK_CLOSE(area, 2.2458050314935392, 0.001);
 583     }
 584
 585 }
 586
 587 int test_main(int, char* [])
 588 {
 589
 590     test_spherical_geo<double>();
 591
 592     return 0;
 593 }