Remove Sliver Polygons Using Shapely, Matplotlib in Python-CodePudding

I have code that starts with points, buffers them into polygons, and finds overlapping regions and makes them into their own polygons. I ran into problems trying to color each polygon individually, because there are 7 very, very small polygons near the center of the overlapping circles (near 4, 0). I decided to plot the centroids of the polygons to help me with the coloring, but this is getting weird and complicated and I think I should figure out how to remove these tiny polygons rather than work around them. I have tried following [these directions][1] (dilating and eroding by a small factor [eps]), but I get an error since my "circles" object is a list. Am I doing it wrong?

Is there a way to remove polygons under a certain size?

Here's my code:

import matplotlib.pyplot as plt
from shapely.geometry import Point, LineString
from shapely.ops import unary_union, polygonize

from matplotlib.pyplot import cm
import numpy as np

def plot_coords(coords, color):
    pts = list(coords)
    x, y = zip(*pts)
    # print(color)
    plt.plot(x,y, color='k', linewidth=1)
    plt.fill_between(x, y, facecolor=color)


def plot_polys(polys, color):
    for poly, color in zip(polys, color):
        plot_coords(poly.exterior.coords, color)
# start with points

points = [Point(6, 0),
          Point(5, 1.73205080757),
          Point(3, 1.73205080757),
          Point(2, 0),
          Point(3, -1.73205080757),
          Point(5, -1.73205080757)]

# buffer points to create circle polygons

circles = []
for point in points:
    circles.append(point.buffer(2))

# attempt at dilating and eroding for slivers, gives error
# eps = 0.001
# circles.buffer(eps).buffer(-eps)


# unary_union and polygonize to find overlaps

rings = [LineString(list(pol.exterior.coords)) for pol in circles]
union = unary_union(rings)
result = [geom for geom in polygonize(union)]

print(len(result))


pink = "#ffbeff" 
blue = "#71e2ff"
green = "#98e68e"


colors = [
    pink, blue, pink, # 0, 1, 2
    green, blue, pink, # 3, 4, 5
    pink, pink, pink, # 6, 7, 8
    pink, green, blue, # 9, 10, 11
    pink, green, blue, # 12, 13, 14
    pink, pink, pink, # 15, 16, 17
    pink, pink, pink, # 18, 19, 20
    pink, pink, pink, pink] #21, 22, 23, 24

# using pre-made color palettes
# colors = cm.viridis(np.linspace(0, 1, len(result)))


fig = plt.figure()
ax = fig.add_subplot()
fig.subplots_adjust(top=0.85)

# determine centroids of polygons in result object
for poly in result:
        centroidCoords = poly.centroid.wkt
        print(centroidCoords)

# centroids for all of the polygons generated. 
# You can see the 7 slivers I have commented out 
# that have centroids all very close to (0,4). 

centroidCoords = [Point(7.207241767209244, 0), #0
Point(5.810892076402134, -1.0455444375165859), #1
Point(5.603603552320669, -2.7775889642668727), #2
Point(4.500000000000001, -0.8660254037850003), #3
Point(4.000000000000002, -2.091010554244934),  #4
# Point(4.000081836086183, -0.0024979955436593), #5
# Point(4.001428127882854, -0.0000000000000001), #6
# Point(3.999999999999999, -0.0000000000000016), #7
Point(5.810892076402136, 1.045544437516586),   #8, #5
# Point(4.000081836086183, 0.0024979955436591),  #9
Point(4.500000000000001, 0.8660254037849999),  #10, #6
Point(3.9999999999999987, 2.0910105542449338), #11, #7
Point(5.603603552320668, 2.777588964266872),   #12, #8
Point(4.999999999999998, 0),                   #13, #9
# Point(3.9999181639138155, 0.0024979955436591), #14
Point(3.500000000000001, 0.8660254037850004),  #15, #10
Point(2.1891079235978643, 1.045544437516586),  #16, #11
Point(2.396396447679332, 2.777588964266872),   #17, #12
# Point(3.9985718721171457, -0.0000000000000001),#18
Point(3, 0),                                   #19, #13
Point(2.1891079235978648, -1.045544437516586), #20, #14
Point(0.7927582327907562, 0.0000000000000001), #21, #15
# Point(3.9999181639138155, -0.0024979955436593),#22, #16
Point(3.5, -0.8660254037849997),               #23, #16
Point(2.3963964476793325, -2.777588964266871)] #24, #17


# plot polygons
plot_polys(result, colors)

# plot centroids
for i in range(len(centroidCoords)):
    plt.scatter(centroidCoords[i].x,centroidCoords[i].y, marker="$" str(i) "$")


ax.set_aspect('equal')
plt.show()


  [1]: https://gis.stackexchange.com/questions/120286/removing-small-polygon-gaps-in-shapely-polygon

CodePudding user response：

Polygons in shapely have a area attribute:

threshold = ... # some value
filtered = [polygon for polygon in result if polygon.area > threshold]