From 61ac71521d00d1dc401b14e97a810f11dbcf9d26 Mon Sep 17 00:00:00 2001 From: Jochen Topf Date: Thu, 14 May 2026 22:01:55 +0200 Subject: [PATCH] Add spherical_length() function to geometries in Lua For linestrings and multilinestrings this returns the same value as the PostGIS ST_Length() function on a geography. For other geometry types this always returns 0.0. --- src/flex-lua-geom.cpp | 20 ++++++++++++ src/geom-functions.cpp | 37 ++++++++++++++++++++-- src/geom-functions.hpp | 11 +++++++ tests/bdd/flex/geometry-linestring.feature | 10 ++++-- tests/test-geom-linestrings.cpp | 23 ++++++++++++++ tests/test-geom-multilinestrings.cpp | 3 ++ tests/test-geom-multipoints.cpp | 1 + tests/test-geom-multipolygons.cpp | 1 + tests/test-geom-null.cpp | 1 + tests/test-geom-points.cpp | 1 + tests/test-geom-polygons.cpp | 1 + 11 files changed, 103 insertions(+), 6 deletions(-) diff --git a/src/flex-lua-geom.cpp b/src/flex-lua-geom.cpp index 6bfee5293..cd6290d74 100644 --- a/src/flex-lua-geom.cpp +++ b/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}}); } diff --git a/src/geom-functions.cpp b/src/geom-functions.cpp index 434f88586..d0aaf20d5 100644 --- a/src/geom-functions.cpp +++ b/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>; + double spherical_area(polygon_t const &geom) { - using sph_point = boost::geometry::model::point< - double, 2, boost::geometry::cs::geographic>; - boost::geometry::model::polygon 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_geom; + boost::geometry::convert(geom, sph_geom); + + return static_cast(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) diff --git a/src/geom-functions.hpp b/src/geom-functions.hpp index 6bcae6918..50b11679e 100644 --- a/src/geom-functions.hpp +++ b/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 diff --git a/tests/bdd/flex/geometry-linestring.feature b/tests/bdd/flex/geometry-linestring.feature index 71c4dc972..96c397657 100644 --- a/tests/bdd/flex/geometry-linestring.feature +++ b/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 diff --git a/tests/test-geom-linestrings.cpp b/tests/test-geom-linestrings.cpp index 630af4d7e..4b74bc363 100644 --- a/tests/test-geom-linestrings.cpp +++ b/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{{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()[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)); +} diff --git a/tests/test-geom-multilinestrings.cpp b/tests/test-geom-multilinestrings.cpp index 25a9c2f03..398aef5d5 100644 --- a/tests/test-geom-multilinestrings.cpp +++ b/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(); 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(); 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(); REQUIRE(l.num_geometries() == 1); REQUIRE(l == expected); diff --git a/tests/test-geom-multipoints.cpp b/tests/test-geom-multipoints.cpp index b85d55df8..2a01d668f 100644 --- a/tests/test-geom-multipoints.cpp +++ b/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}); diff --git a/tests/test-geom-multipolygons.cpp b/tests/test-geom-multipolygons.cpp index 3a6a7534f..8ad673291 100644 --- a/tests/test-geom-multipolygons.cpp +++ b/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{ diff --git a/tests/test-geom-null.cpp b/tests/test-geom-null.cpp index 545c3efc4..ebbfeb7bf 100644 --- a/tests/test-geom-null.cpp +++ b/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()); diff --git a/tests/test-geom-points.cpp b/tests/test-geom-points.cpp index 8e6d5de5c..7aaaf03a4 100644 --- a/tests/test-geom-points.cpp +++ b/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); diff --git a/tests/test-geom-polygons.cpp b/tests/test-geom-polygons.cpp index 58622c467..4758fecd6 100644 --- a/tests/test-geom-polygons.cpp +++ b/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);