/
polylabel.py
139 lines (112 loc) · 4.57 KB
/
polylabel.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
from heapq import heappop, heappush
from ..errors import TopologicalError
from ..geometry import Point
class Cell:
"""A `Cell`'s centroid property is a potential solution to finding the pole
of inaccessibility for a given polygon. Rich comparison operators are used
for sorting `Cell` objects in a priority queue based on the potential
maximum distance of any theoretical point within a cell to a given
polygon's exterior boundary.
"""
def __init__(self, x, y, h, polygon):
self.x = x
self.y = y
self.h = h # half of cell size
self.centroid = Point(x, y) # cell centroid, potential solution
# distance from cell centroid to polygon exterior
self.distance = self._dist(polygon)
# max distance to polygon exterior within a cell
self.max_distance = self.distance + h * 1.4142135623730951 # sqrt(2)
# rich comparison operators for sorting in minimum priority queue
def __lt__(self, other):
return self.max_distance > other.max_distance
def __le__(self, other):
return self.max_distance >= other.max_distance
def __eq__(self, other):
return self.max_distance == other.max_distance
def __ne__(self, other):
return self.max_distance != other.max_distance
def __gt__(self, other):
return self.max_distance < other.max_distance
def __ge__(self, other):
return self.max_distance <= other.max_distance
def _dist(self, polygon):
"""Signed distance from Cell centroid to polygon outline. The returned
value is negative if the point is outside of the polygon exterior
boundary.
"""
inside = polygon.contains(self.centroid)
distance = self.centroid.distance(polygon.exterior)
for interior in polygon.interiors:
distance = min(distance, self.centroid.distance(interior))
if inside:
return distance
return -distance
def polylabel(polygon, tolerance=1.0):
"""Finds pole of inaccessibility for a given polygon. Based on
Vladimir Agafonkin's https://github.com/mapbox/polylabel
Parameters
----------
polygon : shapely.geometry.Polygon
tolerance : int or float, optional
`tolerance` represents the highest resolution in units of the
input geometry that will be considered for a solution. (default
value is 1.0).
Returns
-------
shapely.geometry.Point
A point representing the pole of inaccessibility for the given input
polygon.
Raises
------
shapely.errors.TopologicalError
If the input polygon is not a valid geometry.
Example
-------
>>> from shapely import LineString
>>> polygon = LineString([(0, 0), (50, 200), (100, 100), (20, 50),
... (-100, -20), (-150, -200)]).buffer(100)
>>> polylabel(polygon, tolerance=10).wkt
'POINT (59.35615556364569 121.83919629746435)'
"""
if not polygon.is_valid:
raise TopologicalError("Invalid polygon")
minx, miny, maxx, maxy = polygon.bounds
width = maxx - minx
height = maxy - miny
cell_size = min(width, height)
h = cell_size / 2.0
cell_queue = []
# First best cell approximation is one constructed from the centroid
# of the polygon
x, y = polygon.centroid.coords[0]
best_cell = Cell(x, y, 0, polygon)
# Special case for rectangular polygons avoiding floating point error
bbox_cell = Cell(minx + width / 2.0, miny + height / 2, 0, polygon)
if bbox_cell.distance > best_cell.distance:
best_cell = bbox_cell
# build a regular square grid covering the polygon
x = minx
while x < maxx:
y = miny
while y < maxy:
heappush(cell_queue, Cell(x + h, y + h, h, polygon))
y += cell_size
x += cell_size
# minimum priority queue
while cell_queue:
cell = heappop(cell_queue)
# update the best cell if we find a better one
if cell.distance > best_cell.distance:
best_cell = cell
# continue to the next iteration if we can't find a better solution
# based on tolerance
if cell.max_distance - best_cell.distance <= tolerance:
continue
# split the cell into quadrants
h = cell.h / 2.0
heappush(cell_queue, Cell(cell.x - h, cell.y - h, h, polygon))
heappush(cell_queue, Cell(cell.x + h, cell.y - h, h, polygon))
heappush(cell_queue, Cell(cell.x - h, cell.y + h, h, polygon))
heappush(cell_queue, Cell(cell.x + h, cell.y + h, h, polygon))
return best_cell.centroid