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, 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>
  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
 407 int test_main(int, char* [])
 408 {
 409
 410     test_spherical_geo<double>();
 411
 412     return 0;
 413 }