Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Filter Shapely Polygons By Centroid, Clip, or Something Else?

Tags:

python

shapely

I drew this flower of life by buffering points to polygons. I wanted each overlapping region to be its own polygon, so I used union and polygonize on the lines.

plot1

I have filtered the polygons by area to eliminate sliver polygons, and now I'd like to filter them again and am stuck. I only want to keep the circles that are complete, so the first circle at 0,0 and the first level of surrounding rings (or petals). I want circles like this:

enter image description here

I am wondering if I can filter by centroid location, something like:

complete_polys = [polygon for polygon in filtered_polys if centroid[i].x < 4]
complete_polys = [polygon for polygon in complete_polys_x if centroid[i].x > -4]

Obviously this doesn't work, and I don't even know if it is possible. Perhaps this is the wrong approach entirely, and maybe snap() or clip_by_rect() might be better options?

Thanks in advance for your insight and help.

Here's the code to generate the circles:

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)

x = 0
y = 0
h = 1.73205080757

points = [# center
          Point(x, y),
          #  first ring
          Point((x + 2), y),
          Point((x - 2), y),
          Point((x + 1), (y + h)),
          Point((x - 1), (y + h)),
          Point((x + 1), (y - h)),
          Point((x - 1), (y - h)),
          # second ring
          Point((x + 3), h),
          Point((x - 3), h),
          Point((x + 3), -h),
          Point((x - 3), -h),
          Point((x + 2), (h + h)),
          Point((x - 2), (h + h)),
          Point((x + 2), (-h + -h)),
          Point((x - 2), (-h + -h)),
          Point((x + 4), y),
          Point((x - 4), y),
          Point(x, (h + h)),
          Point(x, (-h + -h)),
          #third ring
          Point((x + 4), (h + h)),
          Point((x - 4), (h + h)),
          Point((x + 4), (-h + -h)),
          Point((x - 4), (-h + -h)),
          Point((x + 1), (h + h + h)),
          Point((x - 1), (h + h + h)),
          Point((x + 1), (-h + -h + -h)),
          Point((x - 1), (-h + -h + -h)),
          Point((x + 5), h),
          Point((x - 5), h),
          Point((x + 5), -h),
          Point((x - 5), -h)]

# buffer points to create circle polygons

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


# unary_union and polygonize to find overlaps

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

# remove tiny sliver polygons
threshold = 0.01
filtered_polys = [polygon for polygon in result_polys if polygon.area > threshold]


print("total polygons = " + str(len(result_polys)))
print("filtered polygons = " + str(len(filtered_polys)))


colors = cm.viridis(np.linspace(0, 1, len(filtered_polys)))


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


plot_polys(filtered_polys, colors)


ax.set_aspect('equal')
plt.show()
like image 694
Doda Avatar asked Oct 16 '25 08:10

Doda


1 Answers

Is this what you wanted?

I used the x^2 + y^2 = r^2 circle formula to filter.

complete_polys = [polygon for polygon in filtered_polys if (polygon.centroid.x**2 + polygon.centroid.y**2 < 4**2)]

plot_polys(complete_polys, colors)

enter image description here

like image 81
11574713 Avatar answered Oct 18 '25 22:10

11574713