Add spherical_length() function to geometries in Lua

src/flex-lua-geom.cpp

@@ -100,6 +100,25 @@ int geom_length(lua_State *lua_state)
     return 1;
 }
 
+int geom_spherical_length(lua_State *lua_state)
+{
+    auto const *const input_geometry = unpack_geometry(lua_state);
+
+    if (input_geometry->srid() != PROJ_LATLONG) {
+        throw std::runtime_error{"Can only calculate spherical length for "
+                                 "geometries in WGS84 (4326) coordinates."};
+    }
+
+    try {
+        lua_pushnumber(lua_state, geom::spherical_length(*input_geometry));
+    } catch (...) {
+        return luaL_error(lua_state,
+                          "Unknown error in 'spherical_length()'.\n");
+    }
+
+    return 1;
+}
+
 int geom_centroid(lua_State *lua_state)
 {
     auto const *const input_geometry = unpack_geometry(lua_state);
@@ -334,6 +353,7 @@ void init_geometry_class(lua_State *lua_state)
          {"segmentize", geom_segmentize},
          {"simplify", geom_simplify},
          {"spherical_area", geom_spherical_area},
+         {"spherical_length", geom_spherical_length},
          {"srid", geom_srid},
          {"transform", geom_transform}});
 }

src/geom-functions.cpp

@@ -367,18 +367,28 @@ double area(geometry_t const &geom)
 
 namespace {
 
+using sph_point = boost::geometry::model::point<
+    double, 2, boost::geometry::cs::geographic<boost::geometry::degree>>;
+
 double spherical_area(polygon_t const &geom)
 {
-    using sph_point = boost::geometry::model::point<
-        double, 2, boost::geometry::cs::geographic<boost::geometry::degree>>;
-
     boost::geometry::model::polygon<sph_point> sph_geom;
     boost::geometry::convert(geom, sph_geom);
+
     return boost::geometry::area(sph_geom,
                                  boost::geometry::strategy::area::geographic<
                                      boost::geometry::strategy::vincenty>{});
 }
 
+double spherical_length(linestring_t const &geom)
+{
+    boost::geometry::model::linestring<sph_point> sph_geom;
+    boost::geometry::convert(geom, sph_geom);
+
+    return static_cast<double>(boost::geometry::length(
+        sph_geom, boost::geometry::strategy::distance::vincenty<>{}));
+}
+
 } // anonymous namespace
 
 double spherical_area(geometry_t const &geom)
@@ -403,6 +413,27 @@ double spherical_area(geometry_t const &geom)
         [](auto const & /*input*/) { return 0.0; }}));
 }
 
+double spherical_length(geometry_t const &geom)
+{
+    assert(geom.srid() == PROJ_LATLONG);
+
+    return geom.visit(overloaded{
+        [](geom::collection_t const &input) {
+            return std::accumulate(input.cbegin(), input.cend(), 0.0,
+                                   [](double sum, auto const &geom) {
+                                       return sum + spherical_length(geom);
+                                   });
+        },
+        [](geom::linestring_t const &input) { return spherical_length(input); },
+        [](geom::multilinestring_t const &input) {
+            return std::accumulate(input.cbegin(), input.cend(), 0.0,
+                                   [](double sum, auto const &geom) {
+                                       return sum + spherical_length(geom);
+                                   });
+        },
+        [](auto const & /*input*/) { return 0.0; }});
+}
+
 /****************************************************************************/
 
 double length(geometry_t const &geom)

src/geom-functions.hpp

@@ -161,6 +161,17 @@ double area(geometry_t const &geom);
  */
 double spherical_area(geometry_t const &geom);
 
+/**
+ * Calculate length of geometry on the spheroid.
+ * For geometry types other than linestring or multilinestring this will always
+ * return 0.
+ *
+ * \param geom Input geometry.
+ * \returns Length in m.
+ * \pre \code geom.srid() == 4326 \endcode
+ */
+double spherical_length(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/bdd/flex/geometry-linestring.feature

@@ -17,6 +17,8 @@ Feature: Creating linestring features from way
                 { column = 'mgeom', type = 'multilinestring', projection = 4326 },
                 { column = 'xgeom', type = 'multilinestring', projection = 4326 },
                 { column = 'npoints', type = 'int' },
+                { column = 'length', type = 'real' },
+                { column = 'slength', type = 'real' },
             })
 
             function osm2pgsql.process_way(object)
@@ -26,6 +28,8 @@ Feature: Creating linestring features from way
                         mgeom = object:as_multilinestring(),
                         xgeom = object:as_linestring(),
                         npoints = object:as_linestring():n_points(),
+                        length = object:as_linestring():length(),
+                        slength = object:as_linestring():spherical_length(),
                     })
                 end
             end
@@ -34,9 +38,9 @@ Feature: Creating linestring features from way
         When running osm2pgsql flex
 
         Then table osm2pgsql_test_lines contains exactly
-            | way_id | sgeom!geo | mgeom!geo   | xgeom!geo   | npoints |
-            | 20     | 1, 2, 3   | [ 1, 2, 3 ] | [ 1, 2, 3 ] | 3       |
-            | 21     | 4, 5      | [ 4, 5 ]    | [ 4, 5 ]    | 2       |
+            | way_id | sgeom!geo | mgeom!geo   | xgeom!geo   | npoints | length     | slength   |
+            | 20     | 1, 2, 3   | [ 1, 2, 3 ] | [ 1, 2, 3 ] | 3       | 0.24142136 | 25718.176 |
+            | 21     | 4, 5      | [ 4, 5 ]    | [ 4, 5 ]    | 2       | 0.14142136 | 15235.885 |
 
     Scenario:
         Given the grid

tests/test-geom-linestrings.cpp

@@ -58,6 +58,7 @@ TEST_CASE("line geometry", "[NoDB]")
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(1.41421));
+    REQUIRE(spherical_length(geom) == Approx(156876.14940188668).epsilon(0.0000001));
     REQUIRE(geometry_type(geom) == "LINESTRING");
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});
     REQUIRE(geometry_n(geom, 1) == geom);
@@ -91,6 +92,7 @@ TEST_CASE("create_linestring from OSM data", "[NoDB]")
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(1.41421));
+    REQUIRE(spherical_length(geom) == Approx(156876.14940188668).epsilon(0.0000001));
     REQUIRE(geom.get<geom::linestring_t>() ==
             geom::linestring_t{{1, 1}, {2, 2}});
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});
@@ -361,3 +363,24 @@ TEST_CASE("geom::simplify of straight line", "[NoDB]")
         REQUIRE(l[1] == input.get<geom::linestring_t>()[2]);
     }
 }
+
+TEST_CASE("long line length - equator", "[NoDB]")
+{
+    geom::geometry_t const geom{geom::linestring_t{{0, 0}, {180, 0}}};
+    REQUIRE(length(geom) == Approx(180.0));
+    REQUIRE(spherical_length(geom) == Approx(20003931.458625447).epsilon(0.0000001));
+}
+
+TEST_CASE("long line length - to pole", "[NoDB]")
+{
+    geom::geometry_t const geom{geom::linestring_t{{0, -90}, {0, 90}}};
+    REQUIRE(length(geom) == Approx(180.0));
+    REQUIRE(spherical_length(geom) == Approx(20003931.458625447).epsilon(0.0000001));
+}
+
+TEST_CASE("line length - more points", "[NoDB]")
+{
+    geom::geometry_t const geom{geom::linestring_t{{20, 19.8}, {20.1, 19.8}, {20.2, 19.9}}};
+    REQUIRE(length(geom) == Approx(0.2414213562373079));
+    REQUIRE(spherical_length(geom) == Approx(25718.175297824535).epsilon(0.0000001));
+}

tests/test-geom-multilinestrings.cpp

@@ -38,6 +38,7 @@ TEST_CASE("create_multilinestring with single line", "[NoDB]")
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(1.0));
+    REQUIRE(spherical_length(geom) == Approx(111302.64933943082));
     auto const &ml = geom.get<geom::multilinestring_t>();
     REQUIRE(ml.num_geometries() == 1);
     REQUIRE(ml[0] == expected);
@@ -65,6 +66,7 @@ TEST_CASE("create_multilinestring with single line and no force_multi",
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(1.0));
+    REQUIRE(spherical_length(geom) == Approx(111302.64933943082));
     auto const &l = geom.get<geom::linestring_t>();
     REQUIRE(l.num_geometries() == 1);
     REQUIRE(l == expected);
@@ -99,6 +101,7 @@ TEST_CASE(
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(1.0));
+    REQUIRE(spherical_length(geom) == Approx(111302.64933943082));
     auto const &l = geom.get<geom::linestring_t>();
     REQUIRE(l.num_geometries() == 1);
     REQUIRE(l == expected);

tests/test-geom-multipoints.cpp

@@ -35,6 +35,7 @@ TEST_CASE("multipoint_t with a single point", "[NoDB]")
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(0.0));
+    REQUIRE(spherical_length(geom) == Approx(0.0));
     REQUIRE(reverse(geom) == geom);
     REQUIRE(centroid(geom) == geom::geometry_t{point});
 

tests/test-geom-multipolygons.cpp

@@ -32,6 +32,7 @@ TEST_CASE("multipolygon geometry with single outer, no inner", "[NoDB]")
     REQUIRE(area(geom) == Approx(1.0));
     REQUIRE(spherical_area(geom) == Approx(12308778361.469454).epsilon(0.00001));
     REQUIRE(length(geom) == Approx(0.0));
+    REQUIRE(spherical_length(geom) == Approx(0.0));
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
     REQUIRE(geometry_n(geom, 1) ==
             geom::geometry_t{geom::polygon_t{

tests/test-geom-null.cpp

@@ -24,6 +24,7 @@ TEST_CASE("null geometry", "[NoDB]")
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(0.0));
+    REQUIRE(spherical_length(geom) == Approx(0.0));
     REQUIRE(geometry_type(geom) == "NULL");
     REQUIRE(centroid(geom).is_null());
     REQUIRE(geometry_n(geom, 1).is_null());

tests/test-geom-points.cpp

@@ -71,6 +71,7 @@ TEST_CASE("create_point from OSM data", "[NoDB]")
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(spherical_area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(0.0));
+    REQUIRE(spherical_length(geom) == Approx(0.0));
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.1, 2.2}});
     REQUIRE(geometry_n(geom, 1) == geom);
     REQUIRE(reverse(geom) == geom);

tests/test-geom-polygons.cpp

@@ -28,6 +28,7 @@ TEST_CASE("polygon geometry without inner", "[NoDB]")
     REQUIRE(area(geom) == Approx(1.0));
     REQUIRE(spherical_area(geom) == Approx(12308778361.469454).epsilon(0.00001));
     REQUIRE(length(geom) == Approx(0.0));
+    REQUIRE(spherical_length(geom) == Approx(0.0));
     REQUIRE(geometry_type(geom) == "POLYGON");
     REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
     REQUIRE(geometry_n(geom, 1) == geom);