10

I have few polygons or boxes. Ones intersect each other but some are isolated. For example I have free figures:

import shapely.geometry as sg
import shapely.ops as so
import matplotlib.pyplot as plt
 
r1 = sg.Polygon([(0,0),(0,1),(1,1),(1,0),(0,0)])
r2 = sg.box(0.5,0.5,1.5,1.5)
r3 = sg.box(4,4,5,5)

First two are intersected and r3 is on some distance from them. I unite them via cascaded_union:

new_shape = so.cascaded_union([r1, r2, r3])

Then I try to plot it (one united figure of r1 and r2 and one distanced box r3)

xs, ys = new_shape.exterior.xy
fig, axs = plt.subplots()
axs.fill(xs, ys, alpha=0.5, fc='r', ec='none')
plt.show()

Instead of plot I receive an AttributeError: 'MultiPolygon' object has no attribute 'exterior'. Is there a pythonic way to display a multipolygon or to iterate through it and draw its parts?

Kadir Şahbaz
78.6k57 gold badges260 silver badges407 bronze badges
asked Mar 5, 2020 at 15:12

4 Answers 4

20

Shapely Polygon object has attribute exterior. Shapely MultiPolygon object has Polygon object sequence. You should iterate over those polygons. You can do that using attribute geoms of MultiPolygon.

Use this way:

import shapely.geometry as sg
import shapely.ops as so
import matplotlib.pyplot as plt
r1 = sg.Polygon([(0,0),(0,1),(1,1),(1,0),(0,0)])
r2 = sg.box(0.5,0.5,1.5,1.5)
r3 = sg.box(4,4,5,5)
new_shape = so.unary_union([r1, r2, r3])
fig, axs = plt.subplots()
axs.set_aspect('equal', 'datalim')
for geom in new_shape.geoms: 
 xs, ys = geom.exterior.xy 
 axs.fill(xs, ys, alpha=0.5, fc='r', ec='none')
plt.show()

enter image description here

answered Mar 5, 2020 at 15:55
1
  • cascaded_union is deprecated. See here. Commented Nov 27, 2022 at 12:56
13

An alternative, shorter way of plotting using @Kadir Şahbaz's answer:

Edit: Per @petezurich's comment, cascaded_union is superseded by unary_union as of 1.8.0

new_shape = so.unary_union([r1, r2, r3])
# Plot each polygon shape directly
for geom in new_shape.geoms:
 plt.plot(*geom.exterior.xy)
# Set (current) axis to be equal before showing plot
plt.gca().axis("equal")
plt.show()
answered Apr 8, 2020 at 20:33
1
  • 1
    cascaded_union is deprecated. See here. Commented Nov 27, 2022 at 12:55
6

Look at Plot shapefile with islands with matplotlib for example.

As with polygons you can use matplotlib paths and patches and there is a Python module dedicated to plot polygons from shapefiles using these functions Descartes.

new_shape= so.unary_union([r1, r2, r3])
from descartes import PolygonPatch
import matplotlib.pyplot as plt
BLUE = '#6699cc'
GRAY = '#999999'
fig = plt.figure() 
ax = fig.gca() 
ax.add_patch(PolygonPatch(new_shape, fc=GRAY, ec=BLUE, alpha=0.5, zorder=2 ))
ax.axis('scaled')
plt.show()

enter image description here

answered Mar 5, 2020 at 16:26
3

You can use GeoPandas' .plot() function:

import shapely.geometry as sg
import shapely.ops as so
import matplotlib.pyplot as plt
import geopandas as gpd
r1 = sg.Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
r2 = sg.box(0.5, 0.5, 1.5, 1.5)
r3 = sg.box(4, 4, 5, 5)
new_shape = so.cascaded_union([r1, r2, r3])
gpd.GeoSeries(new_shape).plot()
plt.show()
answered Mar 21, 2022 at 11:15
1
  • This is the easiest solution, I'd wager. Commented May 15, 2023 at 8:20

Your Answer

Draft saved
Draft discarded

Sign up or log in

Sign up using Google
Sign up using Email and Password

Post as a guest

Required, but never shown

Post as a guest

Required, but never shown

By clicking "Post Your Answer", you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.