Pole of Inaccessibility via Shapely's PolyLabel - Contained within the given polygon?

[diehl]

I recently wrote code to render a visualization with representative points marked inside a collection of polygons. I used Shapely's polylabel function to select representative points inside the polygons. Perusing the results, I started to notice examples where the points produced by polylabel fell outside of their corresponding polygons. Am I misinterpreting the documentation by assuming the point should always fall inside the polygon? Maybe there's an issue with how I'm setting the tolerance parameter? It's not clear how to interpret that parameter from the documentation. Any pointers would be greatly appreciated!

Update: Looking at the original algorithm description, it states "a fast algorithm for finding polygon pole of inaccessibility, the most distant internal point from the polygon outline." Based on that statement, it appears the behavior I'm seeing is unexpected, unless the approximation algorithm is not guaranteed to produce an internal point.

[jorisvandenbossche]

My interpretation of the documentation is also that the center point should always be inside the polygon. In GEOS, this is called Maximum Inscribed Circle, and quoting its docs:

Constructs the "maximum inscribed circle" (MIC) for a polygonal geometry,
up to a specified tolerance.
The MIC is determined by a point in the interior of the area
which has the farthest distance from the area boundary, along with a boundary point at that distance.
In the context of geography the center of the MIC is known as the
"pole of inaccessibility". A cartographic use case is to determine a suitable point
to place a map label within a polygon.

Now, currently we have a custom implementation in Python, and are not using the GEOS function (there is a draft PR for doing that: #1307).

Do you have an example of such a polygon where the result is not inside the polygon?

[diehl]

@jorisvandenbossche as promised, here's a geojson file with the 190 geometries for which the pole of inaccessibility is not contained within each corresponding polygon: https://gist.github.com/diehl/ed12dfeb531b3f9a1e19ba61573e6acc please let me know if you need anything else!

[diehl]

@jorisvandenbossche do you have a sense of when the PR might be integrated into the library? I might need to find an interim solution depending on the timing.

[jorisvandenbossche]

Checking the first geometry of the GeoJSON file you provided, I can reproduce that the returned point is not inside the polygon. Comparing it with the implementation based on GEOS in #1307, I got the same result though, but by doing that noticed that the default value for tolerance is:

tolerance represents the highest resolution in units of the
input geometry that will be considered for a solution. (default
value is 1.0).

So it's an absolute tolerance in the units of the geometry, not some relative number. And your coordinates are in degrees, so the value of 1 is huge. Passing a much small tolerance value ensures the result is actually inside the polygon.
In #1307, the new shapely.maximum_inscribed_circle function that I am adding therefore also has a non-default tolerance, but by default calculated from the extent of the geometry (to exactly avoid this mistake).

Illustration:

import json
import shapely

with open("poi_violations.geojson") as f:
    poi = json.load(f)
polygon = shapely.geometry.shape(poi["features"][0]["geometry"])

from shapely.algorithms.polylabel import polylabel
from shapely.plotting import plot_polygon, plot_points

p1 = polylabel(polygon)
p2 = polylabel(polygon, tolerance=1e-4)

plot_polygon(polygon)
plot_points(p1, color="C1")
plot_points(p2, color="C2")

The reason for the wrong result is that the iterative algoritms starts with the center of the bounding box of the geometry, and with a tolerance of 1.0 it already stops directly at that point, I assume.

[diehl]

@jorisvandenbossche Brilliant!! I was wondering about the impact of tolerance since the Shapely documentation said nothing about how that is defined. In a moment of synchronicity today, a kind soul pointed me to this Shapely issue from a little over a year ago where you pointed out the same issue. Maybe a near-term fix is to add the tolerance definition from above to the documentation? I will circle back to my data and transform the coordinates to a system where I can specify the tolerance in meters which should hopefully resolve my issue. Thank you for digging into this!!

[diehl]

This indeed did the trick. I've got the results I need now. Thanks again @jorisvandenbossche!