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