CODEWITHSUNDEEP
pythonpolygonregionshapely
Matthie456
Matthie456

Reputation: 81

python region growing with polygons performance

I'm trying to create regions of polygons on the condition that they touch. In my example I have an example dataset with 382 polygons that need to be grouped together (but the full dataset contains 6355 polygons). (I would show a picture, but I don't have enough reputation to do that..)

I though of doing this brute force, but of course that takes very long and is not very optimal.

def groupBuildings(blds):
    # blds is a list with shapely polygons
    groups = []
    for bld in blds:
        group = []
        group.append(bld)
        for other in blds:
            for any in group:
                if any != other and any.intersects(other):
                    group.append(other)
        groups.append(group)
    return groups

I learned about region growing and thought that that would be a possible solution, but still the performance is terrible. I've implemented this in the following way:

def groupBuildings(blds):
    # blds is a list with shapely polygons
    others = blds
    groups = []
    while blds != []:
        done = []
        group = []
        first = blds.pop(0)
        done.append(first)
        group.append(first)
        for other in others:
            if (other in blds) and first.touches(other):
                group.append(other)
                blds.remove(other)

return groups

But I think the problem here is that I don't have any nearest neighbors, so I still have to iterate over every building twice.

So my question is: are nearest neighbors essential for region growing? Or is there another way of doing this efficiently?

Upvotes: 1

Views: 453

Answers (1)

Logan Byers
Logan Byers

Reputation: 1573

You will be best served using shapely.ops.cascaded_union() (docs here).

from shapely.geometry import Point, Polygon, MultiPolygon
from shapely.ops import cascaded_union
import numpy as np
polygons = [Point(200*x,200*y).buffer(b) for x,y,b in np.random.random((6000,3))]
multi = MultiPolygon(polygons)
unioned = cascaded_union(multi)

%%timeit
unioned = cascaded_union(multi) 
# 2.8 seconds for me

Upvotes: 1

Related Questions