How to get a list of every Point inside a MultiPolygon using Shapely

Question

I have the following MultiPolygon:

MULTIPOLYGON (
(
(10.8849956 49.8901705, 10.8849507 49.8902499, 10.884969 49.8902588, 10.8851033 49.8903298, 10.8851183 49.8903132, 10.88512882654868 49.8903054, 10.8851246 49.8903054, 10.8851246 49.8902754, 10.8851546 49.8902754, 10.8851546 49.89028643275958, 10.8853289 49.8901612, 10.885421 49.8901035, 10.8854414638889 49.8900896, 10.8854205 49.8900896, 10.8854205 49.8900596, 10.8854505 49.8900596, 10.8854505 49.89008346226415, 10.885527 49.8900315, 10.885519 49.8899952, 10.8854851 49.8899903, 10.8853164 49.8899957, 10.8852419 49.8899981, 10.8851711 49.8899919, 10.8851165 49.8899814, 10.8850728 49.8899652, 10.8850692 49.8899713, 10.8849925 49.8900275, 10.8850251 49.890083, 10.8850275 49.8901159, 10.8850185 49.8901733, 10.8849956 49.8901705),
(10.8852028 49.8901715, 10.8852328 49.8901715, 10.8852328 49.8902015, 10.8852028 49.8902015, 10.8852028 49.8901715),
(10.8852889 49.8900884, 10.8853146 49.8900884, 10.8853146 49.8901184, 10.8853078 49.8901184, 10.8853078 49.8901337, 10.8852808 49.8901337, 10.8852808 49.8901463, 10.8852508 49.8901463, 10.8852508 49.8901346, 10.8852239 49.8901346, 10.8852239 49.8901046, 10.8852305 49.8901046, 10.8852305 49.8900815, 10.8852589 49.8900815, 10.8852589 49.8900812, 10.8852889 49.8900812, 10.8852889 49.8900884),
(10.8851133 49.890201, 10.8851433 49.890201, 10.8851433 49.890231, 10.8851133 49.890231, 10.8851133 49.890201),
(10.8849849 49.8902202, 10.8850149 49.8902202, 10.8850149 49.8902502, 10.8849849 49.8902502, 10.8849849 49.8902202)
),

(
(10.8852605 49.8901112, 10.8852605 49.8901115, 10.8852539 49.8901115, 10.8852539 49.8901163, 10.8852778 49.8901163, 10.8852778 49.8901112, 10.8852605 49.8901112)
)
)

How do I get a flat list containing every single (unique) point? The list items do not have to be Shapely Points, instead they could also be tuples. I don't really understand how to iterate over this structure.

Answer 1

Unfortunately, Shapely doesn't provide the functionality to extract immediately all points from a MultiPolygon object. Instead, you would have to first, iterate over individual polygons of a MultiPolygon, and second, extract individual points of each Polygon.

One could come up with different ways to approach the problem. For example, if you know that none of your polygons have holes, you could simply do the following:

points = []
for polygon in multipolygon:
    points.extend(polygon.exterior.coords[:-1])

Note the [:-1] which prevents duplicating the first vertex. You can remove it if you'd like to have a cleaner syntax and don't care about having one duplicate point for each polygon.
This can also be written in one line using list comprehension with two loops:

points = [point for polygon in multipolygon for point in polygon.exterior.coords[:-1]]

or with the help of itertools.chain.from_iterable:

from itertools import chain

points = list(chain.from_iterable(polygon.exterior.coords[:-1] for polygon in multipolygon))

---

In general, when the polygons can contain holes, we could, for example, write the following function to extract coordinates from the interior rings:

def to_coords(multipolygon):
    for polygon in multipolygon:
        yield from polygon.exterior.coords[:-1]
        yield from chain.from_iterable(interior.coords[:-1] for interior in polygon.interiors)

Example of usage:

mp = MultiPolygon([Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
                   Polygon([(2, 0), (3, 0), (3, 1), (2, 1)], 
                           holes=[[(2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75)]])])

points = list(to_coords(mp))
# [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0),
#  (2.0, 0.0), (3.0, 0.0), (3.0, 1.0), (2.0, 1.0), 
#  (2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75)]

---

One could go even furter and generalize this for any input geometry (Python ≥3.7):

from functools import singledispatch
from itertools import chain
from typing import (List, 
                    Tuple,
                    TypeVar)

from shapely.geometry import (GeometryCollection,
                              LinearRing,
                              LineString,
                              Point,
                              Polygon)
from shapely.geometry.base import (BaseGeometry,
                                   BaseMultipartGeometry)

Geometry = TypeVar('Geometry', bound=BaseGeometry)


@singledispatch
def to_coords(geometry: Geometry) -> List[Tuple[float, float]]:
    """Returns a list of unique vertices of a given geometry object."""
    raise NotImplementedError(f"Unsupported Geometry {type(geometry)}")


@to_coords.register
def _(geometry: Point):
    return [(geometry.x, geometry.y)]


@to_coords.register
def _(geometry: LineString):
    return list(geometry.coords)


@to_coords.register
def _(geometry: LinearRing):
    return list(geometry.coords[:-1])


@to_coords.register
def _(geometry: BaseMultipartGeometry):
    return list(set(chain.from_iterable(map(to_coords, geometry))))


@to_coords.register
def _(geometry: Polygon):
    return to_coords(GeometryCollection([geometry.exterior, *geometry.interiors]))

Example of usage:

from shapely.geometry import (MultiLineString,
                              MultiPoint,
                              MultiPolygon)

geometry_objects = [Point(0, 0),
                    LineString([(0, 0), (1, 1)]),
                    LinearRing([(0, 0), (1, 0), (1, 1)]),
                    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)], 
                            holes=[[(0.25, 0.25), (0.75, 0.25), (0.75, 0.75), (0.25, 0.75)]]),
                    MultiPoint([(0, 0), (1, 1)]),
                    MultiLineString([LineString([(0, 0), (1, 1)]), LineString([(2, 0), (3, 1)])]),
                    MultiPolygon([Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
                                  Polygon([(2, 0), (3, 0), (3, 1), (2, 1)], 
                                          holes=[[(2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75)]])]),
                    GeometryCollection([Point(0, 0), LineString([(0, 0), (1, 1)])])]

for geometry in geometry_objects:
    print(f"For {geometry.wkt}\nwe got:\n"
          f"{to_coords(geometry)}\n")

Output:

For POINT (0 0)
we got:
[(0.0, 0.0)]

For LINESTRING (0 0, 1 1)
we got:
[(0.0, 0.0), (1.0, 1.0)]

For LINEARRING (0 0, 1 0, 1 1, 0 0)
we got:
[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0)]

For POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0), (0.25 0.25, 0.75 0.25, 0.75 0.75, 0.25 0.75, 0.25 0.25))
we got:
[(0.0, 1.0), (0.0, 0.0), (0.25, 0.25), (0.75, 0.25), (0.75, 0.75), (0.25, 0.75), (1.0, 0.0), (1.0, 1.0)]

For MULTIPOINT (0 0, 1 1)
we got:
[(0.0, 0.0), (1.0, 1.0)]

For MULTILINESTRING ((0 0, 1 1), (2 0, 3 1))
we got:
[(2.0, 0.0), (0.0, 0.0), (3.0, 1.0), (1.0, 1.0)]

For MULTIPOLYGON (((0 0, 1 0, 1 1, 0 1, 0 0)), ((2 0, 3 0, 3 1, 2 1, 2 0), (2.25 0.25, 2.75 0.25, 2.75 0.75, 2.25 0.75, 2.25 0.25)))
we got:
[(0.0, 1.0), (0.0, 0.0), (3.0, 0.0), (3.0, 1.0), (2.0, 1.0), (2.0, 0.0), (2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75), (1.0, 0.0), (1.0, 1.0)]

For GEOMETRYCOLLECTION (POINT (0 0), LINESTRING (0 0, 1 1))
we got:
[(0.0, 0.0), (1.0, 1.0)]