1

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()

enter image description here

Please can someone help me out - what I am 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.)

Thanks in advance

4

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()

result

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.