GEOSPointOnSurface_r returns a point that doesn't intersect with the input polygon

[theroggy]

Related to shapely/shapely#1988

For a certain valid polygon, GEOSPointOnSurface_r returns a point that doesn't intersect with the polygon. I seems like a robustness issue... but nonetheless...

This is an image of the polygon (in blue) and the point returned (in red):

Script that reproduces:

import matplotlib.pyplot as plt
from shapely import wkt
from shapely import plotting

poly = wkt.loads(
    "POLYGON ((2653.9 425.9100000000003, 2653.8999999999996 1067.39, 1331.61 1067.3900000000003, 1331.61 1223.5586000000003, 2653.9 1223.5586000000003, 2653.9 425.9100000000003))"
)
print(f"{poly.is_valid=}") # prints True
rp = poly.representative_point()
print(f"{rp=}") # prints False
print(f"{poly.intersects(rp)=}") # prints False

plotting.plot_polygon(poly)
plotting.plot_points(rp, color="red")
plt.show()

[dr-jts]

Here's a geosop reproducer:

geosop -a "POLYGON ((2653.9 425.9100000000003, 2653.8999999999996 1067.39, 1331.61 1067.3900000000003, 1331.61 1223.5586000000003, 2653.9 1223.5586000000003, 2653.9 425.9100000000003))" interiorPoint | geosop -a "POLYGON ((2653.9 425.9100000000003, 2653.8999999999996 1067.39, 1331.61 1067.3900000000003, 1331.61 1223.5586000000003, 2653.9 1223.5586000000003, 2653.9 425.9100000000003))" -b stdin contains

The generated point is POINT (2653.8999999999996 746.6500000000002):

Problem Analysis

There is no floating-point number between 2653.8999999999996 and 2653.9, so the current Interior Point algorithm is doomed to fail for this input, since it uses only one stabbing line.

Workaround

A workaround is to use GEOSMaximumInscribedCircle(). Even a large tolerance (say 1) produces a good result.

Solutions

Perhaps the Interior Point algorithm should check if the point is contained, and if not rerun using MaximumInscribedCircle with a fairly large tolerance. Or, more than one stabbing line can be used (horizontally and vertically).

It will always be possible to create a very thin pathological geometry which fails, however.