Skip to content

Commit e01ecdf

Browse files
jotojoto
committed
Use spheroid instead of sphere for spherical_area()
Use spheroid instead of sphere for spherical_area() calculation. This might be a bit more expensive to calculate but is nearer to the default that PostGIS ST_Area returns when called with a geography and no parameter. The old behaviour was using the use_spheroid = false version of ST_Area. See https://postgis.net/docs/manual-3.5/ST_Area.html
1 parent 8b7d01a commit e01ecdf

8 files changed

Lines changed: 21 additions & 12 deletions

‎src/geom-functions.cpp‎

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -369,16 +369,14 @@ namespace {
369369

370370
double spherical_area(polygon_t const &geom)
371371
{
372-
boost::geometry::strategy::area::spherical<> const spherical_earth{
373-
6371008.8};
374-
375372
using sph_point = boost::geometry::model::point<
376-
double, 2,
377-
boost::geometry::cs::spherical_equatorial<boost::geometry::degree>>;
373+
double, 2, boost::geometry::cs::geographic<boost::geometry::degree>>;
378374

379375
boost::geometry::model::polygon<sph_point> sph_geom;
380376
boost::geometry::convert(geom, sph_geom);
381-
return boost::geometry::area(sph_geom, spherical_earth);
377+
return boost::geometry::area(sph_geom,
378+
boost::geometry::strategy::area::geographic<
379+
boost::geometry::strategy::vincenty>{});
382380
}
383381

384382
} // anonymous namespace

‎tests/test-geom-linestrings.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,7 @@ TEST_CASE("line geometry", "[NoDB]")
5555
REQUIRE(dimension(geom) == 1);
5656
REQUIRE(num_geometries(geom) == 1);
5757
REQUIRE(area(geom) == Approx(0.0));
58+
REQUIRE(spherical_area(geom) == Approx(0.0));
5859
REQUIRE(length(geom) == Approx(1.41421));
5960
REQUIRE(geometry_type(geom) == "LINESTRING");
6061
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});
@@ -87,6 +88,7 @@ TEST_CASE("create_linestring from OSM data", "[NoDB]")
8788
REQUIRE(dimension(geom) == 1);
8889
REQUIRE(num_geometries(geom) == 1);
8990
REQUIRE(area(geom) == Approx(0.0));
91+
REQUIRE(spherical_area(geom) == Approx(0.0));
9092
REQUIRE(length(geom) == Approx(1.41421));
9193
REQUIRE(geom.get<geom::linestring_t>() ==
9294
geom::linestring_t{{1, 1}, {2, 2}});

‎tests/test-geom-multilinestrings.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ TEST_CASE("create_multilinestring with single line", "[NoDB]")
3535
REQUIRE(dimension(geom) == 1);
3636
REQUIRE(num_geometries(geom) == 1);
3737
REQUIRE(area(geom) == Approx(0.0));
38+
REQUIRE(spherical_area(geom) == Approx(0.0));
3839
REQUIRE(length(geom) == Approx(1.0));
3940
auto const &ml = geom.get<geom::multilinestring_t>();
4041
REQUIRE(ml.num_geometries() == 1);
@@ -60,6 +61,7 @@ TEST_CASE("create_multilinestring with single line and no force_multi",
6061
REQUIRE(geometry_type(geom) == "LINESTRING");
6162
REQUIRE(num_geometries(geom) == 1);
6263
REQUIRE(area(geom) == Approx(0.0));
64+
REQUIRE(spherical_area(geom) == Approx(0.0));
6365
REQUIRE(length(geom) == Approx(1.0));
6466
auto const &l = geom.get<geom::linestring_t>();
6567
REQUIRE(l.num_geometries() == 1);
@@ -93,6 +95,7 @@ TEST_CASE(
9395
REQUIRE(geometry_type(geom) == "LINESTRING");
9496
REQUIRE(num_geometries(geom) == 1);
9597
REQUIRE(area(geom) == Approx(0.0));
98+
REQUIRE(spherical_area(geom) == Approx(0.0));
9699
REQUIRE(length(geom) == Approx(1.0));
97100
auto const &l = geom.get<geom::linestring_t>();
98101
REQUIRE(l.num_geometries() == 1);

‎tests/test-geom-multipoints.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ TEST_CASE("multipoint_t with a single point", "[NoDB]")
3232
REQUIRE(dimension(geom) == 0);
3333
REQUIRE(num_geometries(geom) == 1);
3434
REQUIRE(area(geom) == Approx(0.0));
35+
REQUIRE(spherical_area(geom) == Approx(0.0));
3536
REQUIRE(length(geom) == Approx(0.0));
3637
REQUIRE(reverse(geom) == geom);
3738
REQUIRE(centroid(geom) == geom::geometry_t{point});
@@ -55,6 +56,7 @@ TEST_CASE("multipoint_t with several points", "[NoDB]")
5556
REQUIRE(geometry_type(geom) == "MULTIPOINT");
5657
REQUIRE(num_geometries(geom) == 3);
5758
REQUIRE(area(geom) == Approx(0.0));
59+
REQUIRE(spherical_area(geom) == Approx(0.0));
5860
REQUIRE(length(geom) == Approx(0.0));
5961
REQUIRE(reverse(geom) == geom);
6062
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{2, 1}});
@@ -92,6 +94,7 @@ TEST_CASE("create_multipoint from OSM data", "[NoDB]")
9294
REQUIRE(c[3] == geom::point_t{3, 1});
9395

9496
REQUIRE(area(geom) == Approx(0.0));
97+
REQUIRE(spherical_area(geom) == Approx(0.0));
9598
REQUIRE(length(geom) == Approx(0.0));
9699
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{2, 1}});
97100
}
@@ -118,6 +121,7 @@ TEST_CASE("create_multipoint from OSM data with only a single point", "[NoDB]")
118121
REQUIRE(num_geometries(geom) == 1);
119122
REQUIRE(geom.get<geom::point_t>() == geom::point_t{1, 0});
120123
REQUIRE(area(geom) == Approx(0.0));
124+
REQUIRE(spherical_area(geom) == Approx(0.0));
121125
REQUIRE(length(geom) == Approx(0.0));
122126
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1, 0}});
123127
}

‎tests/test-geom-multipolygons.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ TEST_CASE("multipolygon geometry with single outer, no inner", "[NoDB]")
2929
REQUIRE(dimension(geom) == 2);
3030
REQUIRE(num_geometries(geom) == 1);
3131
REQUIRE(area(geom) == Approx(1.0));
32-
REQUIRE(spherical_area(geom) == Approx(12364031798.5));
32+
REQUIRE(spherical_area(geom) == Approx(12308778361.469454).epsilon(0.00001));
3333
REQUIRE(length(geom) == Approx(0.0));
3434
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
3535
REQUIRE(geometry_n(geom, 1) ==
@@ -59,7 +59,7 @@ TEST_CASE("multipolygon geometry with two polygons", "[NoDB]")
5959
REQUIRE(dimension(geom) == 2);
6060
REQUIRE(num_geometries(geom) == 2);
6161
REQUIRE(area(geom) == Approx(9.0));
62-
REQUIRE(spherical_area(geom) == Approx(111106540105.7));
62+
REQUIRE(spherical_area(geom) == Approx(110615268622.783).epsilon(0.00001));
6363
REQUIRE(length(geom) == Approx(0.0));
6464
}
6565

‎tests/test-geom-null.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@ TEST_CASE("null geometry", "[NoDB]")
2020

2121
REQUIRE(dimension(geom) == 0);
2222
REQUIRE(num_geometries(geom) == 0);
23-
REQUIRE(area(geom) == 0.0);
23+
REQUIRE(area(geom) == Approx(0.0));
24+
REQUIRE(spherical_area(geom) == Approx(0.0));
2425
REQUIRE(length(geom) == Approx(0.0));
2526
REQUIRE(geometry_type(geom) == "NULL");
2627
REQUIRE(centroid(geom).is_null());

‎tests/test-geom-points.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,7 @@ TEST_CASE("create_point from OSM data", "[NoDB]")
6868
REQUIRE(dimension(geom) == 0);
6969
REQUIRE(num_geometries(geom) == 1);
7070
REQUIRE(area(geom) == Approx(0.0));
71+
REQUIRE(spherical_area(geom) == Approx(0.0));
7172
REQUIRE(length(geom) == Approx(0.0));
7273
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.1, 2.2}});
7374
REQUIRE(geometry_n(geom, 1) == geom);

‎tests/test-geom-polygons.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ TEST_CASE("polygon geometry without inner", "[NoDB]")
2525
REQUIRE(dimension(geom) == 2);
2626
REQUIRE(num_geometries(geom) == 1);
2727
REQUIRE(area(geom) == Approx(1.0));
28-
REQUIRE(spherical_area(geom) == Approx(12364031798.5));
28+
REQUIRE(spherical_area(geom) == Approx(12308778361.469454).epsilon(0.00001));
2929
REQUIRE(length(geom) == Approx(0.0));
3030
REQUIRE(geometry_type(geom) == "POLYGON");
3131
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
@@ -41,7 +41,7 @@ TEST_CASE("polygon geometry without inner (reverse)", "[NoDB]")
4141
REQUIRE(dimension(geom) == 2);
4242
REQUIRE(num_geometries(geom) == 1);
4343
REQUIRE(area(geom) == Approx(1.0));
44-
REQUIRE(spherical_area(geom) == Approx(12364031798.5));
44+
REQUIRE(spherical_area(geom) == Approx(12308778361.469454).epsilon(0.00001));
4545
REQUIRE(length(geom) == Approx(0.0));
4646
REQUIRE(geometry_type(geom) == "POLYGON");
4747
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
@@ -65,7 +65,7 @@ TEST_CASE("geom::polygon_t", "[NoDB]")
6565
REQUIRE(dimension(geom) == 2);
6666
REQUIRE(num_geometries(geom) == 1);
6767
REQUIRE(area(geom) == Approx(8.0));
68-
REQUIRE(spherical_area(geom) == Approx(98893356298.4));
68+
REQUIRE(spherical_area(geom) == Approx(98452667625.52686).epsilon(0.00001));
6969
REQUIRE(length(geom) == Approx(0.0));
7070
REQUIRE(geometry_type(geom) == "POLYGON");
7171
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});

0 commit comments

Comments
 (0)