oriented_envelope divide by zero warning when applied to an axes-aligned rectangle

[DOSull]

Reporting this here so it doesn't fall between the cracks.

I found this on shapely (see https://github.com/shapely/shapely/issues/2215) where it shows up as a warning calling minimum_rotated_rectangle on a rectangle parallel to the axes. As reported there:

Python 3.13.1 | packaged by conda-forge | (main, Jan 13 2025, 09:45:31) [Clang 18.1.8 ] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from shapely.geometry import *
>>> Polygon([Point(-100,-100),Point(-100,100),Point(100,100),Point(100,-100)]).minimum_rotated_rectangle
/opt/miniconda3/envs/shapely/lib/python3.13/site-packages/shapely/constructive.py:996: RuntimeWarning: divide by zero encountered in oriented_envelope
  return lib.oriented_envelope(geometry, **kwargs)
<POLYGON ((-100 -100, -100 100, 100 100, 100 -100, -100 -100))>

Apparently it's an Apple Silicon problem only. It's only a warning, and the reported result appears correct. So not urgent, but kind of annoying!

I thought I'd see if it was happening on the CLI using geosop but that became a bit of a rabbit hole. I assume the function in question is minAreaRectangle. In any case it doesn't run into the problem at the command line:

geosop -a "POLYGON ((-100 -100, -100 100, 100 100, 100 -100, -100 -100))" minAreaRectangle 
POLYGON ((-100 -100, -100 100, 100 100, 100 -100, -100 -100))

Which is curious...

[mwtoews]

Until recently (#1228), geosop didn't catch/show floating point exceptions. But if you are able to build and test from main, try the same command in -v verbose mode to see if it shows any other messages. I didn't see any message from Linux, so it could be tested on Apple Silicon.

The algorithm in question is mostly in https://github.com/libgeos/geos/blob/main/src/algorithm/MinimumAreaRectangle.cpp which doesn't show any obvious source of division. Note that pin-pointing the exact place where a FPE occurs can be time-consuming, as I don't know how to register a trace.

[dr-jts]

Indeed this does show a FPE in geosop on Apple silicon:

geosop -v -a "POLYGON ((-100 -100, -100 100, 100 100, 100 -100, -100 -100))" minAreaRectangle 
Input A: WKT literal
Read 1 geometries, 5 vertices  -- 1,478 usec
[ 1] minAreaRectangle: A[1] Polygon( 5 )  -> Polygon( 5 )  --  519 usec
POLYGON ((-100 -100, -100 100, 100 100, 100 -100, -100 -100))
Operation raised floating-point environment flag(s): FE_DIVBYZERO FE_INVALID
Ran 1 minAreaRectangle ops ( 5 vertices)  -- 519 usec    (GEOS 3.14.0dev)

Now, how to find out where that's occurring...

[DOSull]

If I had to guess I'd say it's somewhere in the convex hull stuff that is invoked? But I am not even not an expert on C++ or GEOS, so I'll leave figuring out where it's happening to the experts!

[johnthagen]

These warnings hit us when upgrading from shapely 2.0.7 to 2.1.0. We run our unit tests with -Werror this results in a failure. These warnings are appearing on macOS ARM64 platform. We don't see them in our Linux/x86 CI test output.

[tool.pytest.ini_options]
filterwarnings = [
    "error",
]

A minimal reproducer:

from shapely.geometry.polygon import Polygon

coords = ((0.0, 10.0), (0.0, 11.0), (3.0, 11.0), (3.0, 10.0))
polygon = Polygon(coords)
rect = minimum_rotated_rectangle(polygon)

Results in this warning:

geometry = <POLYGON ((0 10, 0 11, 3 11, 3 10, 0 10))>, kwargs = {}

    @multithreading_enabled
    def _oriented_envelope_geos(geometry, **kwargs):
>       return lib.oriented_envelope(geometry, **kwargs)
E       RuntimeWarning: divide by zero encountered in oriented_envelope

Filtering this warning:

[tool.pytest.ini_options]
filterwarnings = [
    "error",
    "ignore:divide by zero encountered in oriented_envelope:RuntimeWarning:shapely",
]

Reveals another warning:

geometry = <POLYGON ((0 10, 0 11, 3 11, 3 10, 0 10))>, kwargs = {}

    @multithreading_enabled
    def _oriented_envelope_geos(geometry, **kwargs):
>       return lib.oriented_envelope(geometry, **kwargs)
E       RuntimeWarning: invalid value encountered in oriented_envelope

To fully suppress this warning in a targeted way, we had to:

[tool.pytest.ini_options]
filterwarnings = [
    "error",
    "ignore:divide by zero encountered in oriented_envelope:RuntimeWarning:shapely",
    "ignore:invalid value encountered in oriented_envelope:RuntimeWarning:shapely",
]

[pramsey]

In term of raising an FE_DIVBYZERO, I found five tests that did so,

./bin/test_geos_unit geos::algorithm::Intersection 3

./bin/test_geos_unit geos::math::DD 3
./bin/test_geos_unit geos::math::DD 15

./bin/test_geos_unit geos::algorithm::CircularArcIntersector 36
./bin/test_geos_unit geos::algorithm::CircularArcIntersector 41

Fixing CGAlgorithmsDD::intersection should dispense with the first one, but that still leaves a couple that are more fiddly. The DD one is at the core, and raises the question of whether DD should also represent Inf as well as NaN, and if it should be working around double flags. The CircularArcIntersector has the error is being thrown during the calculation of CircularArcs::getCenterImpl in the calculation of Hx and Hy because the value of e is zero, coming from CircularArcIntersectorTest<36> which implies maybe fixing higher up the call stack, although patching at the level of DD would quiet the signal.

[pramsey]

Actually my fix only fixed the one unit test geos::algorithm::Intersection 3. The minAreaRectangle example remains at large. Putting in a trap, I found that the exception happens in Rectangle::createLineForStandardEquation at line 109 ... but only if the compile is a Release build, a Debug build finds no exception at all! I do not see how we can hunt down this flag if it does not manifest where we can fix it.