Thanks so much Dan.
I have modified your code slightly to keep the largest polygon:
if not (p1.is_valid):
ext1 = LineString(p1.exterior.coords)
segs = list(segments(ext1))
noded_segments = cascaded_union(segs)
lgstpoly = 0.0
for patch in polygonize(noded_segments):
if patch.area > lgstpoly:
lgstpoly = patch.area
geom = patch.wkt
size = patch.area
minx,miny,maxx,maxy = patch.bounds
regards,
Robert
>>> "Harasty, Daniel J" <dhara...@appcomsci.com> 30/04/2014 9:53 a.m. >>>
I'm using Shapely 1.2.14; this works for me. Let me know if it works in your
general case.
Dan Harasty
'''
Created on Jul 9, 2012
@author: dharasty
'''
from shapely.wkt import loads
from shapely.geometry import LineString, Polygon
from shapely.ops import cascaded_union, polygonize
p1 = loads('POLYGON((0 1, 1 2, 3 0, 5 2 , 6 1, 0 1))')
print "wkt:", p1.wkt
print
ext1 = LineString(p1.exterior.coords)
def segments(linestring):
'''This function takes a LineString, and yields individual LineStrings
for each of its segments.'''
if isinstance(linestring, LineString):
coords = list(linestring.coords)
for n in range(len(coords)-1):
yield LineString([coords[n], coords[n+1]])
segs = list(segments(ext1))
# cascaded_union() happens to correctly "node" the LinearStrings
# (that is, break them where they intersect)
# which is the required input to polygonize(), below
noded_segments = cascaded_union(segs)
for patch in polygonize(noded_segments):
print "wkt: ", patch
print
-----Original Message-----
From: community-boun...@lists.gispython.org
[mailto:community-boun...@lists.gispython.org] On Behalf Of Robert Sanson
Sent: Tuesday, April 29, 2014 5:25 PM
To: gispython.org community projects
Subject: Re: [Community] bowtie problem
Thanks for your suggestion Dan. Your first solution works if there is a single
bowtie and yields a Multipolygon:
from shapely.wkt import loads
from shapely.geometry import Polygon
poly = loads('POLYGON ((0 0, 2 2, 2 0, 0 2, 0 0))')
poly1 = poly.buffer(0)
p2list = list(poly.exterior.coords)
p2list.reverse()
poly2 = Polygon(p2list).buffer(0)
poly3 = poly1.union(poly2)
geom = poly3.wkt
print geom
produces:
MULTIPOLYGON (((1.0 1.0, 2.0 2., 2.0 0.0, 1.0 1.0)), ((0.0 0.0, 0.0 2.0, 1.0
1.0, 0.0 0.0)))
However, it does not generalise if there are more bowties. See the two attached
images. bowtie1.png produces bowtie2.png
I would be interested in your script for your second solution.
Many thanks,
Robert
>>> "Harasty, Daniel J" <dhara...@appcomsci.com<mailto:dhara...@appcomsci.com>>
>>> 30/04/2014 2:11 a.m. >>>
I a hunch that this is the crux of your problem: Since p1 -- as a "bowtie" --
is non-simple, part of it is "inside out" (with respect to the conventional
winding rules). The buffer() operation is -- perhaps by design -- not returning
the buffered version of the "inside out" portion.
This seems to work (and seems to corroborate the above hunch):
p1r = Polygon(list(reversed(p1.exterior.coords))
p2r = p1r.buffer(0)
Since p1r is "inside out" in the "other sense", the buffer() operation seems to
return the "other part" of the bowtie.
Thus the union of p2 and p2r seem to be what you are looking for.
If that doesn't work for you and your complex polygon, I can think of a way
that is more general, but probably at the expense of using more computations.
Here is an outline:
* Split up the complex/non-simple polygon in to its constituent
LinearRings or LineStrings or line segments.
* Use "cascaded_union" to node all the lines (where they intersect)
* Use "polygonize" to stitch back together constituent polygons
That works for me, too, but the code snippet is a bit long; let me know if
you'd like me to post it, and/or send it to you directly.
Cheers,
Dan Harasty
-----Original Message-----
From:
community-boun...@lists.gispython.org<mailto:community-boun...@lists.gispython.org>
[mailto:community-boun...@lists.gispython.org] On Behalf Of Robert Sanson
Sent: Monday, April 28, 2014 9:57 PM
To: gispython.org community projects
Subject: [Community] bowtie problem
Hi All
I am trying to use shapely to fix an invalid polygon with a massive bowtie:
import shapely
from shapely.wkt import loads
p1 = loads('POLYGON((0 0, 2 2, 2 0, 0 2 ,0 0))')
p2 = p1.buffer(0)
p2.wkt returns 'POLYGON ((1 1, 2 2, 2 0, 1 1))'
I am losing half the original polygon.
How can I get a full valid polygon back?
Many thanks,
Robert Sanson
This email and any attachments are confidential and intended solely for the
addressee(s). If you are not the intended recipient, please notify us
immediately and then delete this email from your system.
This message has been scanned for Malware and Viruses by Websense Hosted
Security.
www.websense.com<http://www.websense.com<http://www.websense.com%3chttp:/www.websense.com>>
_______________________________________________
Community mailing list
Community@lists.gispython.org<mailto:Community@lists.gispython.org<mailto:Community@lists.gispython.org%3cmailto:Community@lists.gispython.org>>
http://lists.gispython.org/mailman/listinfo/community
_______________________________________________
Community mailing list
Community@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community
