On Mar 8, 2010, at 6:23 PM, Ben Hesketh wrote:
Hi,
I'm not sure if I should be asking this here or on the libgeos list,
please point me in the right direction if I'm in the wrong place.
I have a multipolygon that is made up of valid individual polygons
that touch, but don't intersect. When I call is_valid it returns
false. This is the multipolygon in question.
MULTIPOLYGON (((1.0000000000000000 1.0000000000000000,
1.0000000000000000 3.0000000000000000, 1.7500000000000000
3.2500000000000000, 2.5000000000000000 2.5000000000000000,
1.0000000000000000 1.0000000000000000)), ((2.5000000000000000
3.5000000000000000, 1.7500000000000000 3.2500000000000000,
1.0000000000000000 4.0000000000000000, 2.5000000000000000
3.5000000000000000)), ((2.5000000000000000 3.5000000000000000,
4.0000000000000000 4.0000000000000000, 3.2500000000000000
3.2500000000000000, 2.5000000000000000 3.5000000000000000)),
((2.5000000000000000 2.5000000000000000, 3.2500000000000000
3.2500000000000000, 4.0000000000000000 3.0000000000000000,
4.0000000000000000 1.0000000000000000, 2.5000000000000000
2.5000000000000000)), ((1.7500000000000000 3.2500000000000000,
2.5000000000000000 3.5000000000000000, 3.2500000000000000
3.2500000000000000, 2.5000000000000000 2.5000000000000000,
1.7500000000000000 3.2500000000000000)))
What does it mean for a MultiPolygon to be valid? I can't find
anything in the docs that explains.
Thanks,
Ben
--
Co-founder, Compass Engine
http://compassengine.com
I've always found that GEOS is best documented by the PostGIS project:
http://postgis.refractions.net/documentation/manual-svn/ch04.html#OGC_Validity
The first four polygons in that collection make up a valid multi (I
checked with shapely.wkt.loads). They touch, but at exactly one point
only between each other. The problem is the 5th polygon, which touches
all the others along lines (an infinite number of points). I used
matplotlib and Shapely 1.2b3 to visualize the situation:
>>> from shapely.wkt import loads
>>> g = loads('MULTIPOLYGON (((1.0000000000000000
1.0000000000000000 ...')
>>> import pylab
>>> x, y = g.geoms[0].exterior.xy
>>> pylab.plot(x, y)
[<matplotlib.lines.Line2D object at 0x1d8ef10>]
>>> x, y = g.geoms[1].exterior.xy
>>> pylab.plot(x, y, 'r')
[<matplotlib.lines.Line2D object at 0x16447870>]
...
Though technically invalid, many operations will still be supported on
that multipolygon.
>>> w = Polygon([(-1.5, 2), (-1.5, 3.5), (3.5, 3.5), (1.5, 3.5)])
>>> g.intersection(w)
<shapely.geometry.multilinestring.MultiLineString object at 0x16451650>
>>> g.intersection(w).wkt
'MULTILINESTRING ((1.5000000000000000 3.5000000000000000,
2.5000000000000000 3.5000000000000000), (2.5000000000000000
3.5000000000000000, 3.5000000000000000 3.5000000000000000))'
--
Sean
_______________________________________________
Community mailing list
[email protected]
http://lists.gispython.org/mailman/listinfo/community
