[geos-devel] Re: PolygonBuilder::findShell assertion 'shellcount <=
1' failed
kyle cronan
kyle at pbx.org
Wed Feb 2 21:52:21 EST 2011
Hi,
I posted the message below back in august, but have only recently had
a chance to look into this further. Here is a minimal test case that
triggers the assertion failure using Shapely. Below that I have
included details on the debugging I did and a patch that fixes the
issue for me (no idea if it's correct though).
from shapely.wkt import loads
p1 = loads('MULTIPOLYGON (((60.0000000000000000 6.5105151320986412,
44.0044859469790026 11.6931320480208569, 0.0000000000000000
25.9507790663861222, 0.0000000000000000 26.8608278557796467,
0.0000000000000000 29.8387923019253307, 60.0000000000000000
10.3985283676378408, 60.0000000000000000 7.8021345594223774,
60.0000000000000000 6.6570998796460161, 60.0000000000000000
6.5105151320986412)), ((43.3161197496508308 0.0000000000000000,
0.0000000000000000 0.0000000000000000, 0.0000000000000000
14.0346133423735822, 0.0000000000000000 17.9226661292310787,
0.0000000000000000 21.5874865260243638, 34.0258524396557860
6.8981402622972743, 55.3161197496508308 0.0000000000000000,
50.0044466166182886 0.0000000000000000, 43.3161197496508308
0.0000000000000000)), ((13.4455725323347899 36.0000000000000000,
60.0000000000000000 36.0000000000000000, 60.0000000000000000
16.7944518298098018, 60.0000000000000000 16.3644011555093201,
60.0000000000000000 14.0439960304547569, 2.9187843276549756
36.0000000000000000, 11.8945390820010992 36.0000000000000000,
13.4455725323347899 36.0000000000000000)))')
p2 = loads('POLYGON ((50.0044466166182886 0.0000000000000000,
0.0000000000000000 21.5874865260243638, 0.0000000000000000
35.7392139719321804, 60.0000000000000000 13.1838946818537934,
60.0000000000000000 0.0000000000000000, 50.0044466166182886
0.0000000000000000))')
p1.union(p2)
Debug details: I've got a python program that uses GEOS via Shapely.
I'm getting an assertion failure within
geos::operation::overlay::PolygonBuilder::findShell, where it checks
that the number of shells is <= 1. I changed the assert to just
return NULL, so that my program will continue past this occasional
error (that's why you see a breakpoint set on "return NULL" instead of
abort(), below). Here is the backtrace:
Breakpoint 1, geos::operation::overlay::PolygonBuilder::findShell (
this=0x7fffffffb830, minEdgeRings=0x209d540) at PolygonBuilder.cpp:263
263 return NULL;
(gdb) bt
#0 geos::operation::overlay::PolygonBuilder::findShell (this=0x7fffffffb830,
minEdgeRings=0x209d540) at PolygonBuilder.cpp:263
#1 0x00007ffff3135f0d in geos::operation::overlay::PolygonBuilder::buildMinimal
EdgeRings (this=0x7fffffffb830, maxEdgeRings=0x20b6c50,
newShellList=0x7fffffffb838, freeHoleList=0x7fffffffb690)
at PolygonBuilder.cpp:213
#2 0x00007ffff3135b30 in geos::operation::overlay::PolygonBuilder::add (
this=0x7fffffffb830, dirEdges=0x7fffffffb740, nodes=0x7fffffffb720)
at PolygonBuilder.cpp:127
#3 0x00007ffff3135985 in geos::operation::overlay::PolygonBuilder::add (
this=0x7fffffffb830, graph=0x7fffffffb9e8) at PolygonBuilder.cpp:100
#4 0x00007ffff3130d65 in geos::operation::overlay::OverlayOp::computeOverlay (
this=0x7fffffffb930, opCode=geos::operation::overlay::OverlayOp::opUNION)
at OverlayOp.cpp:756
#5 0x00007ffff312f508 in
geos::operation::overlay::OverlayOp::getResultGeometry
(this=0x7fffffffb930,
funcCode=geos::operation::overlay::OverlayOp::opUNION)
at OverlayOp.cpp:188
#6 0x00007ffff312e845 in geos::operation::overlay::OverlayOp::overlayOp (
geom0=0x2094490, geom1=0x207ee90,
opCode=geos::operation::overlay::OverlayOp::opUNION) at OverlayOp.cpp:94
#7 0x00007ffff341a28f in geos::operation::overlay::overlayOp::operator() (
this=0x7fffffffbb00, g0=0x2094490, g1=0x207ee90)
at ../source/headers/geos/operation/overlay/OverlayOp.h:344
#8 0x00007ffff341ac59 in
geos::geom::BinaryOp<geos::operation::overlay::overlayOp>
(g0=0x2094490, g1=0x207ee90, _Op=...)
at ../source/headers/geos/geom/BinaryOp.h:215
#9 0x00007ffff34133a5 in GEOSUnion_r (extHandle=0xa13fa0, g1=0x2094490,
g2=0x207ee90) at geos_ts_c.cpp:1675
#10 0x00007ffff363ce4c in ffi_call_unix64 ()
from /usr/lib/python2.6/lib-dynload/_ctypes.so
#11 0x00007ffff363cbd4 in ffi_call ()
from /usr/lib/python2.6/lib-dynload/_ctypes.so
#12 0x00007ffff3637664 in _CallProc ()
from /usr/lib/python2.6/lib-dynload/_ctypes.so
#13 0x00007ffff362ef33 in ?? () from /usr/lib/python2.6/lib-dynload/_ctypes.so
#14 0x000000000041c9d7 in PyObject_Call ()
...
I figured I could create a minimal test case by getting the WKT for
these geometries and then going back and trying to execute just this
one union operation:
(gdb) frame 9
#9 0x00007ffff34133a5 in GEOSUnion_r (extHandle=0xa13fa0, g1=0x2094490,
g2=0x207ee90) at geos_ts_c.cpp:1675
1675 GeomAutoPtr g3 = BinaryOp(g1, g2,
overlayOp(OverlayOp::opUNION));
(gdb) printf "%s", GEOSGeomToWKT_r(extHandle, g1)
MULTIPOLYGON (((60.0000 ... [see test case]
(gdb) printf "%s", GEOSGeomToWKT_r(extHandle, g2)
POLYGON ((50.0044 ... [see test case]
I looked at these in a visualizer. g1 is made up of three distinct
polygons, each of which is simple. g2 is a polygon that intersects
with two of these. So the union should be a multipolygon with 2
components.
I tried doing this union with Shapely, and with my modified GEOS where
it returns NULL instead of calling abort(), I got the correct result.
But I don't know how this code works, so maybe that's not a proper
fix.
Here is a patch that applies against the latest nightly snapshot. I'm
running with this modification now, so I'd be very interested to hear
if this might cause incorrect results in other circumstances.
--- geos-20110201/src/operation/overlay/PolygonBuilder.cpp
2011-02-02 04:20:04.000000000 +0000
+++ geos-20110201.new/src/operation/overlay/PolygonBuilder.cpp
2011-02-03 02:47:10.796078297 +0000
@@ -258,7 +258,7 @@
#endif
}
}
- assert(shellCount <= 1); // found two shells in MinimalEdgeRing list
+ if (shellCount > 1) return NULL; // found two shells in
MinimalEdgeRing list
return shell;
}
Thanks,
Kyle Cronan
On Mon, Aug 23, 2010 at 3:52 PM, kyle cronan <kyle at pbx.org> wrote:
> Hi everyone,
>
> I ran into this assertion failure while developing a python script
> that uses Shapely 1.2.3 with GEOS 3.2.2. My program makes it most of
> the way through a large dataset and then gives this message:
>
> python: PolygonBuilder.cpp:261: geos::geomgraph::EdgeRing*
> geos::operation::overlay::PolygonBuilder::findShell(std::vector<geos::operation::overlay::MinimalEdgeRing*,
> std::allocator<geos::operation::overlay::MinimalEdgeRing*> >*):
> Assertion `shellCount <= 1' failed.
>
> I figured I'd try to isolate some code that can recreate the problem,
> and see if I can get the same behavior with some test code in C, but
> first I just wanted to ask what this failure is all about. And what
> is a minimal edge ring? Any pointers that could help me understand
> how to debug the problem would be appreciated.
>
> Thanks,
> Kyle Cronan
> <kyle at pbx.org>
>
More information about the geos-devel
mailing list