3

I believe I understand the shapely methodology fairly well at this point. However, I'm struggling with getting my holes to populate to the polygons. I've tried a number of approaches.

What I'm doing is taking a buffer feature that I then want to simplify the edges. So, I take the exterior and the individual interior LineStrings and split them by the points that were used to create the buffer.

All of this works fine. However, I then want to take that exterior and the separate interiors and create a single Polygon feature with the holes. Here's some of the various methods that I have tried:

polygonExterior = Polygon(self.SplitLinesByPoints(inputPolygon.exterior,...)) # returns a LinearRing
polygonInteriors = []
for i in range(len(inputPolygon.interiors)):
    polygonInteriors.append(self.SplitLinesByPoints(inputPolygon.interiors[i],...))

#when I check polygonInteriors, it has 3 valid LinearRings

#Trial 1
newPolygon = Polygon(polygonExteriors, polygonInteriors)
# returns a polygon with the correct exterior and no holes (should be 3 in my test dataset)

#Trial 2
interiorCoords = []
for interior in polygonInteriors:
    interiorCoords.append(list(interior.coords))
newPolygon = Polygon(polygonExteriors, interiorCoords)
# same result as Trial 1

#Trial 3
newPolygon = Polygon(polygonExterior)
for interior in polygonInteriors:
    newPolygon = newPolygon.difference(Polygon(interior))
#same result

#Trial 4 based on first received answer
newPolygon = Polygon(polygonExterior, [[pt for pt in inner.coords] for inner in polygonInteriors])
#same result

# testing output
print(str(len(newPolygon.interiors)))
# gives 0

I output the list of list of tuples that are passed in trial 4 and it looks like this:

[[(3149326.399938697, 13547034.20647637), (3149326.3938959567, 13547034.177421806),...,(3149326.399938697, 13547034.20647637)], [(3149339.754347628, 13546944.667812834), (3149339.6827200023, 13546944.533919774),...,(3149339.754347628, 13546944.667812834)], [(3149334.7216034834, 13547168.84576934), (3149334.729913973, 13547168.803259838),..., (3149334.7216034834, 13547168.84576934)]]

I also just wanted to confirm that the interiors were where they should be. So I created shapefiles for each of the interior circles as well. This is what I expect to get: expectedOutput

As you can see, I have the actual output with the three holes as separate shapefiles.

But, what I get is this:

actual output

That's just the TestInterior2... file that should have the holes.

I'm just really at a loss as to what to try next.

Anyone have any ideas what I'm doing wrong or how to get the holes into the Polygon?

3
  • 1
    You are not sharing all that's necessary here. What is your instance (I see self?) and what does its SplitLinesByPoints method do? Is this a Shapely problem or actually an ArcGIS one? Please reduce your problem to a stackoverflow.com/help/minimal-reproducible-example ! Can you share example data?
    – inc42
    Nov 13, 2019 at 19:24
  • It's a shapely issue, though I think something I'm not loading correctly, not an actual issue with shapely. I was just hoping someone would have some methods that I could use to check or know of other ways to get the holes in the shapely polygons. The SplitLines function is returning a LinearRing and shapely is telling me that the LinearRings are valid so it doesn't seem to matter what the SplitLines is doing. I can also see that the SplitLines is giving me the lines that I expect when I export them to shapefiles. Nov 13, 2019 at 20:00
  • 1
    Please provide the return values from the SplitLinesByPoints for us as reproducable data then. Otherwise it seems impossible to help you.
    – inc42
    Nov 13, 2019 at 20:27

2 Answers 2

2

I think the only issue here is that you need to pass the coordinates to the Polygon constructor, not the list of coordinates. So in trial one, instead of

#Trial 1
newPolygon = Polygon(polygonExteriors, polygonInteriors)

do:

#Trial 1
newPolygon = Polygon(polygonExteriors, [[pt for pt in inner.coords] for inner in polygonInteriors])

As I don't have all of your code, here is an example to demonstrate the solution, which returns the same polygon as inputPolygon:

from shapely.geometry import Polygon

# Example polygon with two holes
inputPolygon = Polygon(((0,0),(10,0),(10,10),(0,10)), (((1,3),(5,3),(5,1),(1,1)), ((9,9),(9,8),(8,8),(8,9))))

polygonExterior = inputPolygon.exterior
polygonInteriors = []
for i in range(len(inputPolygon.interiors)):
    # do the stuff with your polygons
    polygonInteriors.append(inputPolygon.interiors[i])

newPolygon = Polygon(polygonExterior, [[pt for pt in inner.coords] for inner in polygonInteriors])
7
  • I tried that and I still didn't get any interior holes. I'm going to update the OP with screenshots and additional info. Nov 13, 2019 at 13:46
  • 1
    Do you get a valid Polygon with the expected shape when you do Polygon(polygonInteriors[0])?
    – bexi
    Nov 13, 2019 at 16:05
  • I do. All of the interior polygons are valid as separate entities. Nov 13, 2019 at 17:53
  • Can you reproduce my posted example (i.e. do you get a valid polygon including the two holes when you use the example polygon that I posted)?
    – bexi
    Nov 13, 2019 at 18:03
  • Yep. I get this: screenshot Nov 13, 2019 at 18:12
1

I wasn't able to come up with an answer using Shapely.

Instead, I output each of the interior shapes and used the ESRI Erase functions. I'm still not sure why the Difference method didn't work either. But, using the ESRI erase worked and I finally have the output I wanted. It's just a lot more writing temporary shapefiles.

Your Answer

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

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