Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
378 views
in Technique[技术] by (71.8m points)

python - Splitting self-intersecting polygon only returned one polygon in Shapely

I am using Python 3.5 64 bit in Windows 7 64 bit, shapely version 1.5.13.

I have the following code that returned me a self-intersecting polygon:

import numpy as np
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt

x = np.array([ 0.38517325,  0.40859912,  0.43296919,  0.4583215 ,  0.4583215 ,
               0.43296919,  0.40859912,  0.38517325,  0.36265506,  0.34100929])
y = np.array([ 62.5       ,  56.17977528,  39.39698492,   0.        ,
               0.        ,  17.34605377,  39.13341671,  60.4180932 ,
               76.02574417,  85.47008547])
polygon = Polygon(np.c_[x, y])
plt.plot(*polygon.exterior.xy)

Self-intersecting polygon

This is correct. Then I tried to obtain the two individual polygons by using buffer(0):

split_polygon = polygon.buffer(0)
plt.plot(*polygon.exterior.xy)
print(type(split_polygon))
plt.fill(*split_polygon.exterior.xy)

Unfortunately, it only returned of the the two polygons:

Only returned one polygon

Could anyone please help? Thanks!

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

The first step is to close the LineString to make a LinearRing, which is what Polygons are made of.

from shapely.geometry import LineString, MultiPolygon
from shapely.ops import polygonize, unary_union

# original data
ls = LineString(np.c_[x, y])
# closed, non-simple
lr = LineString(ls.coords[:] + ls.coords[0:1])
lr.is_simple  # False

However, note that it is non-simple, since the lines cross to make a bow-tie. (The widely used buffer(0) trick usually does not work for fixing bow-ties in my experience). This is unsuitable for a LinearRing, so it needs further work. Make it simple and MultiLineString with unary_union:

mls = unary_union(lr)
mls.geom_type  # MultiLineString'

Then use polygonize to find the Polygons from the linework:

for polygon in polygonize(mls):
    print(polygon)

Or if you want one MultiPolygon geometry:

mp = MultiPolygon(list(polygonize(mls)))

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...