4

I am using shapely Python package in order to create polygon from several connected lines. My code worked with a simple polygon, but failed when the polygon to be created had interior ring that touches the exterior ring. I need a solution that works in both cases.

I tried also to flatten the lists containing the lines' coordinates and then construct the polygon from points but this failed in my case because while the lines itself are ordered, some lines may be "reversed" (it should be p1-p2, p2-p3, but sometimes it is p1-p2, p3-p2). That is why I am using MultiLineString() and merge() which worked but not universally.

Here is the code:

from shapely.geometry import Polygon, MultiLineString
from shapely import line_merge

# polygon can be created
t = [[(1, 1), (1, 6)], [(1, 6), (6, 6)], [(6, 6), (6, 1)], [(1, 1), (6, 1)]]
tt = line_merge(MultiLineString(t))
p1 = Polygon(tt)

# polygon cannot be created
f = [[(3, 1), (1, 1)], [(1, 1), (1, 6)], [(1, 6), (6, 6)], [(6, 6), (6, 1)], [(6, 1), (3, 1)], [(3, 1), (5, 3)], [(5, 3), (2, 3)], [(3, 1), (2, 3)]]
ff = line_merge(MultiLineString(f))
p2 = Polygon(ff)

and this is the error:


  File ~\.conda\envs\geo_env2\Lib\site-packages\shapely\geometry\polygon.py:230 in __new__
    shell = LinearRing(shell)

  File ~\.conda\envs\geo_env2\Lib\site-packages\shapely\geometry\polygon.py:93 in __new__
    coordinates = np.array([_coords(o) for o in coordinates])

TypeError: 'MultiLineString' object is not iterable

This is how the two merged lines look like. The first can be converted to polygon but the second cannot:

enter image description here

Python 3.12 Shapely 2.0.2

1
  • 3
    The name for this is an "inversion". ArcGIS models it as a special case of a simple polygon (single ring); OGC rules call it a special case of a two-ring polygon (polygon with "hole" that touches).
    – Vince
    Commented Dec 13, 2023 at 14:55

1 Answer 1

4

As the error states, a MultiLinesString object isn't iterable, so the Polygon constructor cannot loop over the points to create the polygon. You can get around this by extracting the coordinates explicitly with e.g. shapely.get_coordinates.

Using this (flat) list of coordinates, you can create a polygon, but it will be invalid because the boundary touches itself. Using shapely.make_valid can fix this.

This is the script with the changes applied:

from matplotlib import pyplot as plt
from shapely.geometry import Polygon, MultiLineString
from shapely import get_coordinates, line_merge, make_valid
import shapely.plotting as plot

# polygon can be created
t = [[(1, 1), (1, 6)], [(1, 6), (6, 6)], [(6, 6), (6, 1)], [(1, 1), (6, 1)]]
tt = line_merge(MultiLineString(t))
p1 = Polygon(tt)

# polygon cannot be created
f = [[(3, 1), (1, 1)], [(1, 1), (1, 6)], [(1, 6), (6, 6)], [(6, 6), (6, 1)], [(6, 1), (3, 1)], [(3, 1), (5, 3)], [(5, 3), (2, 3)], [(3, 1), (2, 3)]]
ff = line_merge(MultiLineString(f))

# A MultiLineString is not iterable, so extract the coordinates explicitly
p2 = Polygon(get_coordinates(ff))
# The polygon is not valid because the coordinates weren't properly structured.
# Using make_valid() fixes this.
p2 = make_valid(p2)

# Plot result
plot.plot_polygon(p2)
plt.show()

This is the resulting polygon:

enter image description here

3
  • This solution answers my question and works for the example coordinates I provided. Thanks. Unfortunately, when I put it in the context of my full script and run with real data another problem emerged. The problem is when try to construct polygon with holes passing two arguments to Polygon(), the external ring coordinates and a list of one or more rings' coordinates of the nested polygons. Then I get strange results such as multipolygons or even geometryCollection consisting of polygons and lines. But this may require posting another question.
    – ABC
    Commented Dec 14, 2023 at 11:27
  • 1
    Yes, there are different ways/interpretations that can be followed when making an invalid polygon valid, and it is not abnormal that this leads to multipolygons or geometrycollections being formed. Indeed best to post a new question with the specific examples if you don't find a satisfactory solution.
    – Pieter
    Commented Dec 14, 2023 at 17:56
  • 1
    I decided to go round the problem. Instead of creating polygon with hole right from the beginning I created two polygons, surrounding and nested, and then found the difference between them (like "erase" in ArcGIS). For both polygons I used your suggested method (get_coordinates() and make_valid()). As I had to do all this for multiple polygons I created two geodataframes and overlayed them using gdf1.overlay(gdf2 , how='symmetric_difference').
    – ABC
    Commented Dec 15, 2023 at 10:48

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.