11

Is it possible to convert a Multipolygon into a Polygon which fills all holes or missing inner areas using Shapely? I have been trying since a while but I can't find it in the documentation. The following image I show an example of a multipolygon with those holes I want to fill, and those squares I want to remove.

Multipolygon with holes

Alberto Bonsanto
  • 17,556
  • 10
  • 64
  • 93
  • 1
    Nice image ("I can see your house from here"). Is it relevant to your question? – Jongware Jan 03 '18 at 17:52
  • @usr2564301 I was trying to show those holes that I pretend to fill. – Alberto Bonsanto Jan 03 '18 at 18:41
  • Mmmm... You mean, the green overlay? So is that map under it part of the problem? As for your problem: it really looks like it's some sort of bitmap. Even if it is not (in origin), it looks like it could trivially be converted to one. Then you could use a bitmap filter. – Jongware Jan 03 '18 at 18:47
  • @usr2564301 yes the green overlay – Alberto Bonsanto Jan 03 '18 at 20:20
  • What are the input data and the expected output? It's not clear if each square is a polygon, or all the big green area is one polygon. It's also not really clear which polygons you want to remove. There is a square that touches the big polygon in two points rather than one, and on the bottom right there are two small polygons that form a ring with the big polygon. Do you want to keep or remove those small polygons? Also, what should happen if there are several polygons with holes? – Georgy Mar 15 '20 at 11:54

2 Answers2

6

An approximate approach could be to:

  1. extract the outer boundary of individual components of the MultiPolygon of interest
  2. expand-shrink each of the outer boundaries in order to fill "holes" which these boundaries approximately encompass, e.g., to handle a "doughnut with a cut"
  3. merge all geometries obtained in previous step

For example:

#!/usr/bin/env python
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import cascaded_union

a = 0.25
delta = 0.49

P = MultiPolygon([
    (
        ((0,0),(0,3),(3,3),(3,2-delta),(2,2-delta),(2,2),(1,2),(1,1),(2,1),(2,1+delta),(3,1+delta),(3,0),(0,0)),
        [((a, a), (1-a,a), (1-a,1-a), (a,1-a), (a,a))]
    )
])


eps = 0.01

omega = cascaded_union([
    Polygon(component.exterior).buffer(eps).buffer(-eps) for component in P
])
for x,y in zip(*omega.exterior.coords.xy):
    print(x, y)

The MultiPolygon P looks like: enter image description here

while the script listed above produces as expected an approximate square with side of length 3, i.e., it fills the hole in the lower-left corner as well as the "empty space" in the center of the MultiPolygon which is rendered equivalent to a hole in the expand-shrink procedure with sufficient high value of parameter eps.

ewcz
  • 12,819
  • 1
  • 25
  • 47
6

If it's enough to fill the holes of some MultiPolygon m, you could do this:

no_holes = MultiPolygon(Polygon(p.exterior) for p in m)

If you also need to fill the holes that arise from touching Polygons inside your MultiPolygon, the following should work:

# Create a polygon `b` that contains `m`
xmin, ymin, xmax, ymax = m.bounds
b = Polygon([(xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin)]).buffer(1, resolution=1)

# Take the complement of `m` in `b`, which is a MultiPolygon.
# Pick the outer polygon and take the complement in `b`.
no_holes = b - min(b - m, key=lambda p: p.bounds)
emulbreh
  • 3,421
  • 23
  • 27