Skip to content

Commit ce7d24a

Browse files
lonvialonvia
authored
Merge pull request #2474 from joto/spherical-length
Add spherical_length() function to geometries in Lua
2 parents e7a33e7 + 61ac715 commit ce7d24a

11 files changed

Lines changed: 103 additions & 6 deletions

‎src/flex-lua-geom.cpp‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -100,6 +100,25 @@ int geom_length(lua_State *lua_state)
100100
return 1;
101101
}
102102

103+
int geom_spherical_length(lua_State *lua_state)
104+
{
105+
auto const *const input_geometry = unpack_geometry(lua_state);
106+
107+
if (input_geometry->srid() != PROJ_LATLONG) {
108+
throw std::runtime_error{"Can only calculate spherical length for "
109+
"geometries in WGS84 (4326) coordinates."};
110+
}
111+
112+
try {
113+
lua_pushnumber(lua_state, geom::spherical_length(*input_geometry));
114+
} catch (...) {
115+
return luaL_error(lua_state,
116+
"Unknown error in 'spherical_length()'.\n");
117+
}
118+
119+
return 1;
120+
}
121+
103122
int geom_centroid(lua_State *lua_state)
104123
{
105124
auto const *const input_geometry = unpack_geometry(lua_state);
@@ -334,6 +353,7 @@ void init_geometry_class(lua_State *lua_state)
334353
{"segmentize", geom_segmentize},
335354
{"simplify", geom_simplify},
336355
{"spherical_area", geom_spherical_area},
356+
{"spherical_length", geom_spherical_length},
337357
{"srid", geom_srid},
338358
{"transform", geom_transform}});
339359
}

‎src/geom-functions.cpp‎

Lines changed: 34 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -367,18 +367,28 @@ double area(geometry_t const &geom)
367367

368368
namespace {
369369

370+
using sph_point = boost::geometry::model::point<
371+
double, 2, boost::geometry::cs::geographic<boost::geometry::degree>>;
372+
370373
double spherical_area(polygon_t const &geom)
371374
{
372-
using sph_point = boost::geometry::model::point<
373-
double, 2, boost::geometry::cs::geographic<boost::geometry::degree>>;
374-
375375
boost::geometry::model::polygon<sph_point> sph_geom;
376376
boost::geometry::convert(geom, sph_geom);
377+
377378
return boost::geometry::area(sph_geom,
378379
boost::geometry::strategy::area::geographic<
379380
boost::geometry::strategy::vincenty>{});
380381
}
381382

383+
double spherical_length(linestring_t const &geom)
384+
{
385+
boost::geometry::model::linestring<sph_point> sph_geom;
386+
boost::geometry::convert(geom, sph_geom);
387+
388+
return static_cast<double>(boost::geometry::length(
389+
sph_geom, boost::geometry::strategy::distance::vincenty<>{}));
390+
}
391+
382392
} // anonymous namespace
383393

384394
double spherical_area(geometry_t const &geom)
@@ -403,6 +413,27 @@ double spherical_area(geometry_t const &geom)
403413
[](auto const & /*input*/) { return 0.0; }}));
404414
}
405415

416+
double spherical_length(geometry_t const &geom)
417+
{
418+
assert(geom.srid() == PROJ_LATLONG);
419+
420+
return geom.visit(overloaded{
421+
[](geom::collection_t const &input) {
422+
return std::accumulate(input.cbegin(), input.cend(), 0.0,
423+
[](double sum, auto const &geom) {
424+
return sum + spherical_length(geom);
425+
});
426+
},
427+
[](geom::linestring_t const &input) { return spherical_length(input); },
428+
[](geom::multilinestring_t const &input) {
429+
return std::accumulate(input.cbegin(), input.cend(), 0.0,
430+
[](double sum, auto const &geom) {
431+
return sum + spherical_length(geom);
432+
});
433+
},
434+
[](auto const & /*input*/) { return 0.0; }});
435+
}
436+
406437
/****************************************************************************/
407438

408439
double length(geometry_t const &geom)

‎src/geom-functions.hpp‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -161,6 +161,17 @@ double area(geometry_t const &geom);
161161
*/
162162
double spherical_area(geometry_t const &geom);
163163

164+
/**
165+
* Calculate length of geometry on the spheroid.
166+
* For geometry types other than linestring or multilinestring this will always
167+
* return 0.
168+
*
169+
* \param geom Input geometry.
170+
* \returns Length in m.
171+
* \pre \code geom.srid() == 4326 \endcode
172+
*/
173+
double spherical_length(geometry_t const &geom);
174+
164175
/**
165176
* Split multigeometries into their parts. Non-multi geometries are left
166177
* alone and will end up as the only geometry in the result vector. If the

‎tests/bdd/flex/geometry-linestring.feature‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@ Feature: Creating linestring features from way
1717
{ column = 'mgeom', type = 'multilinestring', projection = 4326 },
1818
{ column = 'xgeom', type = 'multilinestring', projection = 4326 },
1919
{ column = 'npoints', type = 'int' },
20+
{ column = 'length', type = 'real' },
21+
{ column = 'slength', type = 'real' },
2022
})
2123
2224
function osm2pgsql.process_way(object)
@@ -26,6 +28,8 @@ Feature: Creating linestring features from way
2628
mgeom = object:as_multilinestring(),
2729
xgeom = object:as_linestring(),
2830
npoints = object:as_linestring():n_points(),
31+
length = object:as_linestring():length(),
32+
slength = object:as_linestring():spherical_length(),
2933
})
3034
end
3135
end
@@ -34,9 +38,9 @@ Feature: Creating linestring features from way
3438
When running osm2pgsql flex
3539

3640
Then table osm2pgsql_test_lines contains exactly
37-
| way_id | sgeom!geo | mgeom!geo | xgeom!geo | npoints |
38-
| 20 | 1, 2, 3 | [ 1, 2, 3 ] | [ 1, 2, 3 ] | 3 |
39-
| 21 | 4, 5 | [ 4, 5 ] | [ 4, 5 ] | 2 |
41+
| way_id | sgeom!geo | mgeom!geo | xgeom!geo | npoints | length | slength |
42+
| 20 | 1, 2, 3 | [ 1, 2, 3 ] | [ 1, 2, 3 ] | 3 | 0.24142136 | 25718.176 |
43+
| 21 | 4, 5 | [ 4, 5 ] | [ 4, 5 ] | 2 | 0.14142136 | 15235.885 |
4044

4145
Scenario:
4246
Given the grid

‎tests/test-geom-linestrings.cpp‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,7 @@ TEST_CASE("line geometry", "[NoDB]")
5858
REQUIRE(area(geom) == Approx(0.0));
5959
REQUIRE(spherical_area(geom) == Approx(0.0));
6060
REQUIRE(length(geom) == Approx(1.41421));
61+
REQUIRE(spherical_length(geom) == Approx(156876.14940188668).epsilon(0.0000001));
6162
REQUIRE(geometry_type(geom) == "LINESTRING");
6263
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});
6364
REQUIRE(geometry_n(geom, 1) == geom);
@@ -91,6 +92,7 @@ TEST_CASE("create_linestring from OSM data", "[NoDB]")
9192
REQUIRE(area(geom) == Approx(0.0));
9293
REQUIRE(spherical_area(geom) == Approx(0.0));
9394
REQUIRE(length(geom) == Approx(1.41421));
95+
REQUIRE(spherical_length(geom) == Approx(156876.14940188668).epsilon(0.0000001));
9496
REQUIRE(geom.get<geom::linestring_t>() ==
9597
geom::linestring_t{{1, 1}, {2, 2}});
9698
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]")
361363
REQUIRE(l[1] == input.get<geom::linestring_t>()[2]);
362364
}
363365
}
366+
367+
TEST_CASE("long line length - equator", "[NoDB]")
368+
{
369+
geom::geometry_t const geom{geom::linestring_t{{0, 0}, {180, 0}}};
370+
REQUIRE(length(geom) == Approx(180.0));
371+
REQUIRE(spherical_length(geom) == Approx(20003931.458625447).epsilon(0.0000001));
372+
}
373+
374+
TEST_CASE("long line length - to pole", "[NoDB]")
375+
{
376+
geom::geometry_t const geom{geom::linestring_t{{0, -90}, {0, 90}}};
377+
REQUIRE(length(geom) == Approx(180.0));
378+
REQUIRE(spherical_length(geom) == Approx(20003931.458625447).epsilon(0.0000001));
379+
}
380+
381+
TEST_CASE("line length - more points", "[NoDB]")
382+
{
383+
geom::geometry_t const geom{geom::linestring_t{{20, 19.8}, {20.1, 19.8}, {20.2, 19.9}}};
384+
REQUIRE(length(geom) == Approx(0.2414213562373079));
385+
REQUIRE(spherical_length(geom) == Approx(25718.175297824535).epsilon(0.0000001));
386+
}

‎tests/test-geom-multilinestrings.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ TEST_CASE("create_multilinestring with single line", "[NoDB]")
3838
REQUIRE(area(geom) == Approx(0.0));
3939
REQUIRE(spherical_area(geom) == Approx(0.0));
4040
REQUIRE(length(geom) == Approx(1.0));
41+
REQUIRE(spherical_length(geom) == Approx(111302.64933943082));
4142
auto const &ml = geom.get<geom::multilinestring_t>();
4243
REQUIRE(ml.num_geometries() == 1);
4344
REQUIRE(ml[0] == expected);
@@ -65,6 +66,7 @@ TEST_CASE("create_multilinestring with single line and no force_multi",
6566
REQUIRE(area(geom) == Approx(0.0));
6667
REQUIRE(spherical_area(geom) == Approx(0.0));
6768
REQUIRE(length(geom) == Approx(1.0));
69+
REQUIRE(spherical_length(geom) == Approx(111302.64933943082));
6870
auto const &l = geom.get<geom::linestring_t>();
6971
REQUIRE(l.num_geometries() == 1);
7072
REQUIRE(l == expected);
@@ -99,6 +101,7 @@ TEST_CASE(
99101
REQUIRE(area(geom) == Approx(0.0));
100102
REQUIRE(spherical_area(geom) == Approx(0.0));
101103
REQUIRE(length(geom) == Approx(1.0));
104+
REQUIRE(spherical_length(geom) == Approx(111302.64933943082));
102105
auto const &l = geom.get<geom::linestring_t>();
103106
REQUIRE(l.num_geometries() == 1);
104107
REQUIRE(l == expected);

‎tests/test-geom-multipoints.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ TEST_CASE("multipoint_t with a single point", "[NoDB]")
3535
REQUIRE(area(geom) == Approx(0.0));
3636
REQUIRE(spherical_area(geom) == Approx(0.0));
3737
REQUIRE(length(geom) == Approx(0.0));
38+
REQUIRE(spherical_length(geom) == Approx(0.0));
3839
REQUIRE(reverse(geom) == geom);
3940
REQUIRE(centroid(geom) == geom::geometry_t{point});
4041

‎tests/test-geom-multipolygons.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ TEST_CASE("multipolygon geometry with single outer, no inner", "[NoDB]")
3232
REQUIRE(area(geom) == Approx(1.0));
3333
REQUIRE(spherical_area(geom) == Approx(12308778361.469454).epsilon(0.00001));
3434
REQUIRE(length(geom) == Approx(0.0));
35+
REQUIRE(spherical_length(geom) == Approx(0.0));
3536
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{0.5, 0.5}});
3637
REQUIRE(geometry_n(geom, 1) ==
3738
geom::geometry_t{geom::polygon_t{

‎tests/test-geom-null.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ TEST_CASE("null geometry", "[NoDB]")
2424
REQUIRE(area(geom) == Approx(0.0));
2525
REQUIRE(spherical_area(geom) == Approx(0.0));
2626
REQUIRE(length(geom) == Approx(0.0));
27+
REQUIRE(spherical_length(geom) == Approx(0.0));
2728
REQUIRE(geometry_type(geom) == "NULL");
2829
REQUIRE(centroid(geom).is_null());
2930
REQUIRE(geometry_n(geom, 1).is_null());

‎tests/test-geom-points.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,7 @@ TEST_CASE("create_point from OSM data", "[NoDB]")
7171
REQUIRE(area(geom) == Approx(0.0));
7272
REQUIRE(spherical_area(geom) == Approx(0.0));
7373
REQUIRE(length(geom) == Approx(0.0));
74+
REQUIRE(spherical_length(geom) == Approx(0.0));
7475
REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.1, 2.2}});
7576
REQUIRE(geometry_n(geom, 1) == geom);
7677
REQUIRE(reverse(geom) == geom);

0 commit comments

Comments
 (0)