Fix geographic intersection calculation in sjoberg_intersection when a segment lies on the equator

[fbriol]

PR Description

Description

This PR fixes a bug in the geographic intersection algorithm where the longitude crossing of the equator is calculated incorrectly when one of the segments lies exactly along the equator.

The root cause is in boost/geometry/formulas/sjoberg_intersection.hpp under sjoberg_geodesic::lon_of_equator_intersection(). The formula calculates the offset of longitude from the start point lonj using - asin_tj_t0j + dLj, but forgets to multiply this offset by sign_lon_diff. As a result, when sign_lon_diff is -1, the offset is added in the wrong direction, causing a completely incorrect intersection longitude.

This incorrect coordinate causes overlay algorithms (like clipping or set operations) to fail to detect that a linestring crosses the boundary of a polygon containing equatorial segments.

Fixes #1482

Proposed Changes

Multiply the longitude offset term by sign_lon_diff in lon_of_equator_intersection():

    CT lon_of_equator_intersection() const
    {
        CT const c0 = 0;
        CT const dLj = d_lambda(c0);
        return lonj + sign_lon_diff * (-asin_tj_t0j + dLj);
    }

Reproduction Example

#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/strategies/geographic/intersection.hpp>
#include <boost/geometry/formulas/sjoberg_intersection.hpp>

namespace bg = boost::geometry;

int main() {
    typedef double CT;
    bg::srs::spheroid<double> spheroid;
    
    // Segment 1 (crossing the equator)
    CT lon_a1 = 3.268218994140625 * bg::math::d2r<double>();
    CT lat_a1 = 0.01789649948477745 * bg::math::d2r<double>();
    CT lon_a2 = 3.3428595066070557 * bg::math::d2r<double>();
    CT lat_a2 = -0.4950112998485565 * bg::math::d2r<double>();
    
    // Segment 2 (lying exactly on the equator)
    CT lon_b1 = 4.0 * bg::math::d2r<double>();
    CT lat_b1 = 0.0 * bg::math::d2r<double>();
    CT lon_b2 = 3.0 * bg::math::d2r<double>();
    CT lat_b2 = 0.0 * bg::math::d2r<double>();
    
    typedef bg::strategy::andoyer::template inverse<double, false, true, false, false, false> inverse_type;
    typename inverse_type::result_type res1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid);
    typename inverse_type::result_type res2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid);
    
    CT lon, lat;
    bool ok = bg::formula::sjoberg_intersection<double, bg::strategy::andoyer::template inverse, 4>
              ::apply(lon_a1, lat_a1, lon_a2, lat_a2, res1.azimuth,
                      lon_b1, lat_b1, lon_b2, lat_b2, res2.azimuth,
                      lon, lat, spheroid);
                      
    std::cout << "sjoberg_intersection::apply: " << (ok ? "true" : "false") << std::endl;
    if (ok) {
        std::cout << "  lon (deg): " << lon * bg::math::r2d<double>() << std::endl;
        std::cout << "  lat (deg): " << lat * bg::math::r2d<double>() << std::endl;
    }
    return 0;
}

• Before Patch: Output is lon (deg): 3.41489 (incorrect, outside segment range).
• After Patch: Output is lon (deg): 3.27082 (correct).

[tinko92]

Thanks, test case looks plausible, and change looks consistent with all other sites of this sign-factor.

[barendgehrels]

Thanks! This is fine, I also checked it locally with a union/intersection having the same problem.
I will add that as another unit test, but that can be another PR.

This is approved and can be merged.

[barendgehrels]

Merging, having 2 approvals
@vissarion will I add it to the release notes or can you take it?