Plotting shapely polygon with holes does not plot all holes

Question

I'm trying to plot a filled shapely.geometry.Polygon polygon with holes in. I have tried 2 methods, both give the same result. Below is the code I have been using to debug

from shapely.geometry import Polygon
from descartes import PolygonPatch
import geopandas as gpd
import matplotlib.pyplot as plt


# Polygon exterior:
p = [[20,767],[54,744],[107,707],
 [190,654],[265,609],[363,548],
 [462,484],[514,447],[603,389],
 [682,337],[726,310],[757,290],
 [786,277],[820,259],[843,249],
 [881,231],[921,215],[975,197],
 [1048,174],[1089,163],[1141,152],
 [1212,137],[1270,121],[1271,64],
 [1207,78],[1163,89],[1096,103],
 [1048,115],[1001,129],[949,144],
 [905,157],[874,170],[830,187],
 [781,208],[730,236],[696,255],
 [652,282],[606,306],[561,340],
 [512,370],[478,393],[436,418],
 [385,453],[330,490],[285,521],
 [229,566],[183,603],[123,652],
 [70,698],[13,749]]

# Define interior "holes":
interiors = {}
interiors[0] = [[290,543],[301,560],[393,501],[377,482]]
interiors[1] = [[507,392],[549,363],[553,367],[572,352],[588,372],[522,415]]
interiors[2] = [[599,340],[636,316],[648,334],[612,357]]
interiors[3] = [[714,262],[727,284],[821,238],[811,215]]
interiors[4] = [[850,218],[935,185],[937,187],[850,221]]
interiors[5] = [[959,159],[1066,129],[1071,146],[966,177]]
interiors[6] = [[1119,133],[1175,122],[1178,123],[1119,134]]
interiors[7] = [[1211,102],[1266,91],[1267,97],[1212,108],[1211,102]]

i_p = {k: Polygon(v) for k, v in interiors.items()}


zone = Polygon(p, [zone.exterior.coords for zone in i_p.values() \
                    if zone.within(Polygon(p)) is True])

I now have a Polygon called zone which I think should have the correct exterior and interior. Working in a jupyter notebook I am able to run a cell with just the object name

zone

As hoped for the polygon displays correctly. However I now want to plot this polygon in all sorts of other places, so I look online and try to use descartes.PolygonPatch

fig, ax = plt.subplots()
patch = PolygonPatch(zone)
ax.add_patch(patch)
ax.set_xlim(0, 1300)
ax.set_ylim(0, 777)
plt.show()

The output is now incorrect, two of the interior areas are filled in. Checking that a point that lies in one of these interior holes is not within the zone polygon with zone.contains(<point>) returns a False (as desired). However running the same method on the PolygonPatch created, patch, returned True. Each interior returned True to is_closed().

I then tried geopandas and returned the same

p = gpd.GeoSeries(zone)
p.plot()
plt.show()

What am I doing wrong?

I have also spent hours trying to find what method is running to produce the first plot when the object zone is "run" in interactive python, as I figure that could give me some insight.)

Answer 1

The reason is that exterior and interiors have a canonical form which assumes counter-clockwise coordinates for exterior and clockwise for interiors. You have to normalize your result to see it correctly. Shapely does not have normalize function yet (which is available in PyGEOS) but doing the buffer(0) does the trick.

fig, ax = plt.subplots()
patch = PolygonPatch(zone.buffer(0))
ax.add_patch(patch)
ax.set_xlim(0, 1300)
ax.set_ylim(0, 777)
plt.show()

Answer 2

As previously answered - the interior boundary must be clockwise if the exterior is counterclockwise (need to has opposite rotational direction). Of course it cannot cross itself. This is because matplotlib treats it like a single polygon that has very thin line connecting hole with outer space, but it is not plot as contour.

The simplest way is to invert the list of points for wrongly plotted interior curves:

for index, curve in interiors.items():
    if is_the_other_way_around(curve):
        interiors[index] = curve[::-1]

if you need to check it automatically you can use this method: https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order