Find holes in a union of rectangles?

Question

I have a number of random rectangles (black) in and around a unit square (red) and need to extract all the polygonal regions inside the unit square that are not covered by any rectangle.

enter image description here

It looks like this can be done with Shapely and I've gotten to the point when I have the union of the rectangles (green) but I'm not sure how to subtract that from the unit square and retrieve a list of polygons.

Here is my code to generate the test data:

import pylab
import random
from matplotlib import pyplot
from shapely.geometry import Point, Polygon
from shapely.ops import cascaded_union
from descartes import PolygonPatch

def make_square(x, y, size1, size2):
    dx = [size1, -size1, -size1, size1, size1]
    dy = [size2, size2, -size2, -size2, size2]
    return [(x+sx, y+sy) for sx, sy in zip(dx, dy)]


pylab.figure()

square = make_square(0.5, 0.5, 1.0, 1.0)
a, b = zip(*square)
pylab.plot(a, b, 'r-')
polygons = []

for i in xrange(10):
    x = random.random()
    y = random.random()
    s1 = random.random()
    s2 = random.random()

    square = make_square(x, y, s1, s2)
    polygons.append(Polygon(square))
    a, b = zip(*square)
    pylab.plot(a, b, 'k-')

u = cascaded_union(polygons)
patch2b = PolygonPatch(u, fc='#00ff00', ec='#00ff00', alpha=0.5, zorder=2)
pylab.gca().add_patch(patch2b)


pylab.show()

Can the empty polygon regions returned by the code simply be a list of rectangles that may be adjacent to each other?. Or do you want a list of polygons defined by the number of points making up that region. For example, the three open space regions in your diagram can be accurately described by seven rectangles. Maybe more, maybe less. It becomes much harder of a problem to describe it as three distinct polygons with an arbitrary number of points. — selbie, Aug 19 '14 at 02:23

score 3 · Accepted Answer · answered Aug 19 '14 at 03:00

Basically, you want to take the difference of your "unioned" polygons with the large square and then polygonize the result to get individual, separate polygons. For example:

#-- Get the region not covered by individual squares.
uncovered_region = Polygon(bigsquare).difference(union)

# In some cases, the result will be a single polygon...
if not isinstance(uncovered_region, MultiPolygon):
    uncovered_region = [uncovered_region]

for poly in polygonize(uncovered_region):
    patch = PolygonPatch(poly, fc='purple', alpha=0.5, zorder=2)
    ax.add_patch(patch)

As a full example, based on yours:

import matplotlib.pyplot as plt
import random
from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import cascaded_union, polygonize
from descartes import PolygonPatch

def make_square(x, y, size1, size2):
    dx = [size1, -size1, -size1, size1, size1]
    dy = [size2, size2, -size2, -size2, size2]
    return [(x+sx, y+sy) for sx, sy in zip(dx, dy)]


fig, ax = plt.subplots()

bigsquare = make_square(0.5, 0.5, 1.0, 1.0)
a, b = zip(*bigsquare)
ax.plot(a, b, 'r-')
polygons = []

for i in xrange(10):
    x = random.random()
    y = random.random()
    s1 = random.random()
    s2 = random.random()

    square = make_square(x, y, s1, s2)
    polygons.append(Polygon(square))
    a, b = zip(*square)
    ax.plot(a, b, 'k-')

union = cascaded_union(polygons)
patch2b = PolygonPatch(union, fc='#00ff00', ec='#00ff00', alpha=0.5, zorder=2)
ax.add_patch(patch2b)

#-- Get the region not covered by individual squares.
uncovered_region = Polygon(bigsquare).difference(union)

# In some cases, the result will be a single polygon...
if not isinstance(uncovered_region, MultiPolygon):
    uncovered_region = [uncovered_region]

for poly in polygonize(uncovered_region):
    print poly
    patch = PolygonPatch(poly, fc='purple', alpha=0.5, zorder=2)
    ax.add_patch(patch)

plt.margins(0.05)
plt.show()

enter image description here

Find holes in a union of rectangles?

1 Answers1

Linked