Intersection of geographic segment that crosses the antimeridian with its envelope

[anddmb]

Dear maintainers,

with boost 1.77, the following program

#include <boost/geometry.hpp>

#include <iostream>

namespace bg = boost::geometry;
namespace bm = bg::model::d2;

using point = bm::point_xy<double, bg::cs::geographic<bg::degree> >;
using linestring = bg::model::linestring<point>;
using multi_linestring = bg::model::multi_linestring<linestring>;
using box = bg::model::box<point>;

point a() { return { -179.9958, -47.2835 }; }
point b() { return {  179.9972, -47.2753 }; }

linestring ls() {return {a(), { b() }};}

int main() {
    multi_linestring result;
    auto const envelope = bg::return_envelope<box>(ls());
    bool const success = bg::intersection (envelope, ls(), result);
    std::cout << "input    : " << bg::wkt(ls()) << std::endl;
    std::cout << "envelope : " << bg::wkt(envelope) << std::endl;
    std::cout << "success  : " << success << std::endl;
    std::cout << "output   : " << bg::wkt(result) << std::endl;
    return 0;
}

yields the output (Compiler Explorer)

input    : LINESTRING(-179.996 -47.2835,179.997 -47.2753)
envelope : POLYGON((179.997 -47.2835,179.997 -47.2753,180.004 -47.2753,180.004 -47.2835,179.997 -47.2835))
success  : 1
output   : MULTILINESTRING((179.997 -47.2753))

The envelope appears to be correct. However, the intersection of the segment with its own envelope only yields a single point.

In the case of the intersection, it seems that boost::geometry assumes that the geodesic between a() and b() does not cross the antimeridian.

While this appears to be a bug, is there any simple workaround other than rotating back/forth before/after the transformation?

Thanks a lot in advance! 😄

Kind regards, Andreas

^{Andreas Dirks [email protected] on behalf of MBition GmbH Provider Information}

[awulkiew]

Thanks for the report!

Yes, this is a bug but it's more severe than you suggested. For set operations of linear geometry vs box the use of cartesian strategy is hardcoded: https://github.com/boostorg/geometry/blob/6b74f7c8a344495384cc5ac4f74bb7a307b7c1d1/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp#L682-L685 https://github.com/boostorg/geometry/blob/6b74f7c8a344495384cc5ac4f74bb7a307b7c1d1/include/boost/geometry/algorithms/detail/overlay/clip_linestring.hpp#L39-L51 AFAIS this set operation was always implemented like that.

Since the operation performs cartesian computation it finds the intersection between a cartesian segment (-179.996 -47.2835,179.997 -47.2753) and cartesian box (179.997 -47.2835, 180.004 -47.2753) which is indeed a single point.

The fix would be the implementation of correct strategies for spherical and geographic coordinate systems.

As for the workaround. Now that you know that the operation is in fact cartesian you could account for that by adapting data and performing the operation several times. Another possibility would be the conversion of the box to a polygon. However then its shape would be different because in Boost.Geometry edges of boxes are defined by meridians and parallels but edges of polygons are defined by geodesics. Or maybe yet another approach would be better in your case. It depends what would you like to do with it?