Shapely: Polylabel bug

[benediktbrandt]

Expected behavior and actual behavior.

Holes are not correctly accounted for in shapely's implementation of the polylabel algorithm.

Steps to reproduce the problem.

from shapely.geometry import Polygon
from shapely.algorithms.polylabel import polylabel
poly1=Polygon(shell=[(0,0),(10,0),(10,10),(0,10),(0,0)][::-1], holes=[[(2,2),(6,2),(6,6),(2,6),(2,2)]])
import geopandas
gdf2 = geopandas.GeoDataFrame(geometry=[poly1])
gdf2.plot()
label_location = polylabel(poly1, tolerance=0.05) 
print(label_location.x)
print(label_location.y)

The output of this script is:

6.25
3.75

But the correct position is something close to 7.65625, 7.65625. This can be verified by looking at a plot of this polygon

Maybe I am doing something wrong in the definition of the polygon but I have tried many different ways of instantiating this polygon and I always get the same result.

Operating system

Linux

Shapely version and provenance

Tested with shapely version 1.6 and 1.7a3

[kannes]

Could you please verify that the coordinate that https://github.com/mapbox/polylabel itself suggests is different?

[snorfalorpagus]

Original algorithm in TypeScript:

import polylabel from "polylabel";

let exterior = [[0,0],[10,0],[10,10],[0,10],[0,0]];
let hole = [[2,2],[6,2],[6,6],[2,6],[2,2]];
let polygon = [exterior, hole]

console.log(polylabel(polygon, 0.01));

PS> .\node_modules\.bin\tsc test.ts
PS> node .\test.js
[ 7.65625, 7.65625 ]

Shapely:

>>> from shapely.ops import polylabel
>>> from shapely.geometry import Polygon
>>> exterior = [[0,0],[10,0],[10,10],[0,10],[0,0]]
>>> hole = [[2,2],[6,2],[6,6],[2,6],[2,2]]
>>> p = Polygon(exterior, [hole]) 
>>> polylabel(p, 0.01).coords[:]
[(6.25, 3.75)]

[snorfalorpagus]

You are correct @benediktbrandt. The Shapely implementation doesn't handle the holes correctly, as it only considers the distance to the exterior of the polygon. I think this is a bug this should be fixed to match the behaviour of the original algorithm.

https://github.com/Toblerity/Shapely/blob/cfa66d418fe3afe2fcb48438fcc9156142abaeb8/shapely/algorithms/polylabel.py#L44-L53

The conversion of the polygon exterior to a linestring is a little inefficient here too as it will be done for every cell created.