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