Add spherical_area() function

To calculate the "real" area of a (multi)polygon. This is still an
approximation. Results are similar to what ST_Area(geography, false)
calculates in PostGIS. Accessible from Lua.

flex-config/geometries.lua

@@ -34,7 +34,10 @@ tables.polygons = osm2pgsql.define_area_table('polygons', {
     -- if we have multiple geometry columns and some of them can be valid
     -- and others not.
     { column = 'geom', type = 'geometry', projection = 4326 },
+    -- In this column we'll put the area calculated in Mercator coordinates
     { column = 'area', type = 'real' },
+    -- In this column we'll put the true area calculated on the spheroid
+    { column = 'spherical_area', type = 'real' },
 })
 
 tables.boundaries = osm2pgsql.define_relation_table('boundaries', {
@@ -135,7 +138,9 @@ function osm2pgsql.process_way(object)
             geom = geom,
             -- Calculate the area in Mercator projection and store in the
             -- area column
-            area = geom:transform(3857):area()
+            area = geom:transform(3857):area(),
+            -- Also calculate "real" area in spheroid
+            spherical_area = geom:spherical_area()
         })
     else
         -- We want to split long lines into smaller segments. We can use
@@ -192,7 +197,9 @@ function osm2pgsql.process_relation(object)
             geom = geom,
             -- Calculate the area in Mercator projection and store in the
             -- area column
-            area = geom:transform(3857):area()
+            area = geom:transform(3857):area(),
+            -- Also calculate "real" area in spheroid
+            spherical_area = geom:spherical_area()
         })
     end
 end

src/flex-lua-geom.cpp

@@ -66,6 +66,24 @@ static int geom_area(lua_State *lua_state)
     return 1;
 }
 
+static int geom_spherical_area(lua_State *lua_state)
+{
+    auto const *const input_geometry = unpack_geometry(lua_state);
+
+    if (input_geometry->srid() != 4326) {
+        throw std::runtime_error{"Can only calculate spherical area for "
+                                 "geometries in WGS84 (4326) coordinates."};
+    }
+
+    try {
+        lua_pushnumber(lua_state, geom::spherical_area(*input_geometry));
+    } catch (...) {
+        return luaL_error(lua_state, "Unknown error in 'spherical_area()'.\n");
+    }
+
+    return 1;
+}
+
 static int geom_length(lua_State *lua_state)
 {
     auto const *const input_geometry = unpack_geometry(lua_state);
@@ -286,6 +304,7 @@ void init_geometry_class(lua_State *lua_state)
                         geom_pole_of_inaccessibility);
     luaX_add_table_func(lua_state, "segmentize", geom_segmentize);
     luaX_add_table_func(lua_state, "simplify", geom_simplify);
+    luaX_add_table_func(lua_state, "spherical_area", geom_spherical_area);
     luaX_add_table_func(lua_state, "srid", geom_srid);
     luaX_add_table_func(lua_state, "transform", geom_transform);
 

src/geom-functions.cpp

@@ -364,6 +364,44 @@ double area(geometry_t const &geom)
 
 /****************************************************************************/
 
+static double spherical_area(polygon_t const &geom)
+{
+    boost::geometry::strategy::area::spherical<> const spherical_earth{
+        6371008.8};
+
+    using sph_point = boost::geometry::model::point<
+        double, 2,
+        boost::geometry::cs::spherical_equatorial<boost::geometry::degree>>;
+
+    boost::geometry::model::polygon<sph_point> sph_geom;
+    boost::geometry::convert(geom, sph_geom);
+    return boost::geometry::area(sph_geom, spherical_earth);
+}
+
+double spherical_area(geometry_t const &geom)
+{
+    assert(geom.srid() == 4326);
+
+    return std::abs(geom.visit(overloaded{
+        [&](geom::nullgeom_t const & /*input*/) { return 0.0; },
+        [&](geom::collection_t const &input) {
+            return std::accumulate(input.cbegin(), input.cend(), 0.0,
+                                   [](double sum, auto const &geom) {
+                                       return sum + spherical_area(geom);
+                                   });
+        },
+        [&](geom::polygon_t const &input) { return spherical_area(input); },
+        [&](geom::multipolygon_t const &input) {
+            return std::accumulate(input.cbegin(), input.cend(), 0.0,
+                                   [&](double sum, auto const &geom) {
+                                       return sum + spherical_area(geom);
+                                   });
+        },
+        [&](auto const & /*input*/) { return 0.0; }}));
+}
+
+/****************************************************************************/
+
 double length(geometry_t const &geom)
 {
     return geom.visit(overloaded{

src/geom-functions.hpp

@@ -147,6 +147,17 @@ geometry_t segmentize(geometry_t const &input, double max_segment_length);
  */
 double area(geometry_t const &geom);
 
+/**
+ * Calculate area of geometry on the spheroid.
+ * For geometry types other than polygon or multipolygon this will always
+ * return 0.
+ *
+ * \param geom Input geometry.
+ * \returns Area in m².
+ * \pre \code geom.srid() == 4326 \endcode
+ */
+double spherical_area(geometry_t const &geom);
+
 /**
  * Split multigeometries into their parts. Non-multi geometries are left
  * alone and will end up as the only geometry in the result vector. If the

tests/test-geom-multipolygons.cpp

@@ -21,13 +21,15 @@ TEST_CASE("multipolygon geometry with single outer, no inner", "[NoDB]")
     geom::geometry_t geom{geom::multipolygon_t{}};
     auto &mp = geom.get<geom::multipolygon_t>();
 
+    // MULTIPOLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)))
     mp.add_geometry(
         geom::polygon_t{geom::ring_t{{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}}});
 
     REQUIRE(geometry_type(geom) == "MULTIPOLYGON");
     REQUIRE(dimension(geom) == 2);
     REQUIRE(num_geometries(geom) == 1);
     REQUIRE(area(geom) == Approx(1.0));
+    REQUIRE(spherical_area(geom) == Approx(12364031798.5));
     REQUIRE(length(geom) == Approx(0.0));
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
     REQUIRE(geometry_n(geom, 1) ==
@@ -40,6 +42,7 @@ TEST_CASE("multipolygon geometry with two polygons", "[NoDB]")
     geom::geometry_t geom{geom::multipolygon_t{}};
     auto &mp = geom.get<geom::multipolygon_t>();
 
+    // MULTIPOLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)), ((2 2, 2 5, 5 5, 5 2, 2 2), (3 3, 4 3, 4 4, 3 4, 3 3)))
     mp.add_geometry(
         geom::polygon_t{geom::ring_t{{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}}});
 
@@ -56,6 +59,7 @@ TEST_CASE("multipolygon geometry with two polygons", "[NoDB]")
     REQUIRE(dimension(geom) == 2);
     REQUIRE(num_geometries(geom) == 2);
     REQUIRE(area(geom) == Approx(9.0));
+    REQUIRE(spherical_area(geom) == Approx(111106540105.7));
     REQUIRE(length(geom) == Approx(0.0));
 }
 

tests/test-geom-polygons.cpp

@@ -18,12 +18,14 @@
 
 TEST_CASE("polygon geometry without inner", "[NoDB]")
 {
+    // POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
     geom::geometry_t const geom{
         geom::polygon_t{geom::ring_t{{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}}}};
 
     REQUIRE(dimension(geom) == 2);
     REQUIRE(num_geometries(geom) == 1);
     REQUIRE(area(geom) == Approx(1.0));
+    REQUIRE(spherical_area(geom) == Approx(12364031798.5));
     REQUIRE(length(geom) == Approx(0.0));
     REQUIRE(geometry_type(geom) == "POLYGON");
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
@@ -32,12 +34,14 @@ TEST_CASE("polygon geometry without inner", "[NoDB]")
 
 TEST_CASE("polygon geometry without inner (reverse)", "[NoDB]")
 {
+    // POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))
     geom::geometry_t const geom{
         geom::polygon_t{geom::ring_t{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}}}};
 
     REQUIRE(dimension(geom) == 2);
     REQUIRE(num_geometries(geom) == 1);
     REQUIRE(area(geom) == Approx(1.0));
+    REQUIRE(spherical_area(geom) == Approx(12364031798.5));
     REQUIRE(length(geom) == Approx(0.0));
     REQUIRE(geometry_type(geom) == "POLYGON");
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
@@ -48,6 +52,8 @@ TEST_CASE("geom::polygon_t", "[NoDB]")
     geom::polygon_t polygon;
 
     REQUIRE(polygon.outer().empty());
+
+    // POLYGON((0 0, 0 3, 3 3, 3 0, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))
     polygon.outer() = geom::ring_t{{0, 0}, {0, 3}, {3, 3}, {3, 0}, {0, 0}};
     polygon.inners().emplace_back(
         geom::ring_t{{1, 1}, {2, 1}, {2, 2}, {1, 2}, {1, 1}});
@@ -59,6 +65,7 @@ TEST_CASE("geom::polygon_t", "[NoDB]")
     REQUIRE(dimension(geom) == 2);
     REQUIRE(num_geometries(geom) == 1);
     REQUIRE(area(geom) == Approx(8.0));
+    REQUIRE(spherical_area(geom) == Approx(98893356298.4));
     REQUIRE(length(geom) == Approx(0.0));
     REQUIRE(geometry_type(geom) == "POLYGON");
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});