Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Segmentation fault when intersecting certain geometries using Shapely #329

Closed
Embi opened this issue Jul 7, 2020 · 10 comments
Closed

Segmentation fault when intersecting certain geometries using Shapely #329

Embi opened this issue Jul 7, 2020 · 10 comments

Comments

@Embi
Copy link

Embi commented Jul 7, 2020

Description

Not sure if this is necessarily a GEOS issue - might be just some Shapely -> GEOS interface issue. When using GEOS 3.7.3 the segfault is not happening. The issue was reported to Shapely back in March shapely/shapely#860

Environment

reproduced on Arch Linux and Ubuntu 20.04
Python 3.8
Shapely 1.7.0 (and also reproduced with 1.6.4.post2)
GEOS 3.8.0
GDAL 3.0.4
PROJ 6.3.1-1

Steps to reproduce

#!/usr/bin/env python3

from shapely.geometry import shape

polygon = shape({
    "type": "Polygon",
    "coordinates": [[
        [1.9280898571014473, 41.26269360649398],
        [1.9280898571014473, 41.267195716664716],
        [1.934056631049079, 41.267195716664716],
        [1.934056631049079, 41.26269360649398],
        [1.9280898571014473, 41.26269360649398]
    ]]
})

collection = shape({
    "type": "GeometryCollection",
    "geometries": [{
        "type": "MultiPolygon",
        "coordinates": [
            [[
                [1.933593750000014, 41.26954950284258],
                [1.9226074218749893, 41.26954950284258],
                [1.9226074218749893, 41.26129149391986],
                [1.933593750000014, 41.26129149391986],
                [1.933593750000014, 41.26954950284258]
            ]],
            [[
                [1.9225859642028806, 41.26033982031611],
                [1.9335937499999998, 41.26033982031611],
                [1.9335937499999998, 41.261130194285094],
                [1.9225859642028806, 41.261130194285094],
                [1.9225859642028806, 41.26033982031611]
            ]]
        ]
    }]
})

# Intersecting this particular collection fails
polygon.intersection(collection)

GDB backtrace:

#0  0x00007ffff6aa4d77 in geos::algorithm::locate::SimplePointInAreaLocator::locatePointInPolygon(geos::geom::Coordinate const&, geos::geom::Polygon const*) () from /usr/lib/libgeos-3.8.0.so
#1  0x00007ffff6ad0057 in geos::geomgraph::EdgeEndStar::getLocation(int, geos::geom::Coordinate const&, std::vector<geos::geomgraph::GeometryGraph*, std::allocator<geos::geomgraph::GeometryGraph*> >*) ()
   from /usr/lib/libgeos-3.8.0.so
#2  0x00007ffff6ad02af in geos::geomgraph::EdgeEndStar::computeLabelling(std::vector<geos::geomgraph::GeometryGraph*, std::allocator<geos::geomgraph::GeometryGraph*> >*) () from /usr/lib/libgeos-3.8.0.so
#3  0x00007ffff6acb8ff in geos::geomgraph::DirectedEdgeStar::computeLabelling(std::vector<geos::geomgraph::GeometryGraph*, std::allocator<geos::geomgraph::GeometryGraph*> >*) () from /usr/lib/libgeos-3.8.0.so
#4  0x00007ffff6b2a196 in geos::operation::overlay::OverlayOp::computeLabelling() () from /usr/lib/libgeos-3.8.0.so
#5  0x00007ffff6b2bba8 in geos::operation::overlay::OverlayOp::computeOverlay(geos::operation::overlay::OverlayOp::OpCode) () from /usr/lib/libgeos-3.8.0.so
#6  0x00007ffff6b2befa in geos::operation::overlay::OverlayOp::getResultGeometry(geos::operation::overlay::OverlayOp::OpCode) () from /usr/lib/libgeos-3.8.0.so
#7  0x00007ffff6b2c39e in geos::operation::overlay::OverlayOp::overlayOp(geos::geom::Geometry const*, geos::geom::Geometry const*, geos::operation::overlay::OverlayOp::OpCode) () from /usr/lib/libgeos-3.8.0.so
#8  0x00007ffff6ab2f67 in std::unique_ptr<geos::geom::Geometry, std::default_delete<geos::geom::Geometry> > geos::geom::BinaryOp<geos::operation::overlay::overlayOp>(geos::geom::Geometry const*, geos::geom::Geometry const*, geos::operation::overlay::overlayOp) () from /usr/lib/libgeos-3.8.0.so
#9  0x00007ffff6ab0844 in geos::geom::Geometry::intersection(geos::geom::Geometry const*) const () from /usr/lib/libgeos-3.8.0.so
#10 0x00007ffff6bda3ba in GEOSIntersection_r () from /usr/lib/libgeos_c.so.1
#11 0x00007ffff715d69a in ffi_call_unix64 () from /usr/lib/libffi.so.6
#12 0x00007ffff715cfb6 in ffi_call () from /usr/lib/libffi.so.6
#13 0x00007ffff71bf5b7 in _ctypes_callproc () from /home/ember/.envs/sk/lib/python3.6/lib-dynload/_ctypes.cpython-36m-x86_64-linux-gnu.so
#14 0x00007ffff71b9658 in ?? () from /home/ember/.envs/sk/lib/python3.6/lib-dynload/_ctypes.cpython-36m-x86_64-linux-gnu.so
#15 0x00007ffff7d214c8 in PyObject_Call (func=0x7ffff6c4a110, args=<optimized out>, kwargs=<optimized out>) at Objects/abstract.c:2261
#16 0x00007ffff7e4c5dc in partial_call (pto=0x7ffff67d7b38, args=<optimized out>, kw=<optimized out>) at ./Modules/_functoolsmodule.c:186
#17 0x00007ffff7d214c8 in PyObject_Call (func=0x7ffff67d7b38, args=<optimized out>, kwargs=<optimized out>) at Objects/abstract.c:2261
#18 0x00007ffff7deac03 in do_call_core (kwdict=0x0, callargs=0x7fffe0641208, func=0x7ffff67d7b38) at Python/ceval.c:5106
#19 _PyEval_EvalFrameDefault (f=<optimized out>, throwflag=<optimized out>) at Python/ceval.c:3404
#20 0x00007ffff7de5ec9 in _PyEval_EvalCodeWithName (_co=_co@entry=0x7ffff7128780, globals=globals@entry=0x7ffff711bbd0, locals=locals@entry=0x0, args=args@entry=0x7fffffffd3c0, argcount=argcount@entry=3, 
    kwnames=kwnames@entry=0x0, kwargs=0x0, kwcount=0, kwstep=2, defs=0x0, defcount=0, kwdefs=0x0, closure=0x0, name=0x7ffff73e0170, qualname=0x7ffff67d7ee0) at Python/ceval.c:4166
#21 0x00007ffff7dedf09 in _PyFunction_FastCallDict (func=func@entry=0x7ffff67ef158, args=args@entry=0x7fffffffd3c0, nargs=nargs@entry=3, kwargs=kwargs@entry=0x0) at Python/ceval.c:5070
#22 0x00007ffff7d216f2 in _PyObject_FastCallDict (func=func@entry=0x7ffff67ef158, args=args@entry=0x7fffffffd3c0, nargs=nargs@entry=3, kwargs=kwargs@entry=0x0) at Objects/abstract.c:2310
#23 0x00007ffff7d2195e in _PyObject_Call_Prepend (func=0x7ffff67ef158, obj=<optimized out>, args=0x7fffe063ad88, kwargs=0x0) at Objects/abstract.c:2373
#24 0x00007ffff7d214c8 in PyObject_Call (func=0x7ffff732d948, args=<optimized out>, kwargs=<optimized out>) at Objects/abstract.c:2261
#25 0x00007ffff7d83861 in slot_tp_call (self=self@entry=0x7ffff67f5160, args=args@entry=0x7fffe063ad88, kwds=kwds@entry=0x0) at Objects/typeobject.c:6207
#26 0x00007ffff7d2164a in _PyObject_FastCallDict (func=0x7ffff67f5160, args=<optimized out>, nargs=2, kwargs=kwargs@entry=0x0) at Objects/abstract.c:2331
#27 0x00007ffff7d21c3a in _PyObject_FastCallKeywords (func=func@entry=0x7ffff67f5160, stack=<optimized out>, nargs=nargs@entry=2, kwnames=kwnames@entry=0x0) at Objects/abstract.c:2496
#28 0x00007ffff7de61f7 in call_function (pp_stack=pp_stack@entry=0x7fffffffd5d8, oparg=<optimized out>, kwnames=kwnames@entry=0x0) at Python/ceval.c:4861
#29 0x00007ffff7de84df in _PyEval_EvalFrameDefault (f=<optimized out>, throwflag=<optimized out>) at Python/ceval.c:3335
#30 0x00007ffff7de5494 in _PyFunction_FastCall (co=<optimized out>, args=<optimized out>, nargs=2, globals=<optimized out>) at Python/ceval.c:4919
#31 0x00007ffff7de6167 in fast_function (func=<optimized out>, stack=<optimized out>, nargs=<optimized out>, kwnames=<optimized out>) at Python/ceval.c:4961
#32 0x00007ffff7de6265 in call_function (pp_stack=pp_stack@entry=0x7fffffffd788, oparg=<optimized out>, kwnames=kwnames@entry=0x0) at Python/ceval.c:4858
#33 0x00007ffff7de84df in _PyEval_EvalFrameDefault (f=<optimized out>, throwflag=<optimized out>) at Python/ceval.c:3335
#34 0x00007ffff7de5ec9 in _PyEval_EvalCodeWithName (_co=_co@entry=0x7ffff7352780, globals=globals@entry=0x7ffff7394510, locals=locals@entry=0x7ffff7394510, args=args@entry=0x0, argcount=argcount@entry=0, 
    kwnames=kwnames@entry=0x0, kwargs=0x0, kwcount=0, kwstep=2, defs=0x0, defcount=0, kwdefs=0x0, closure=0x0, name=0x0, qualname=0x0) at Python/ceval.c:4166
#35 0x00007ffff7de640e in PyEval_EvalCodeEx (_co=_co@entry=0x7ffff7352780, globals=globals@entry=0x7ffff7394510, locals=locals@entry=0x7ffff7394510, args=args@entry=0x0, argcount=argcount@entry=0, 
    kws=kws@entry=0x0, kwcount=0, defs=0x0, defcount=0, kwdefs=0x0, closure=0x0) at Python/ceval.c:4187
#36 0x00007ffff7de643c in PyEval_EvalCode (co=co@entry=0x7ffff7352780, globals=globals@entry=0x7ffff7394510, locals=locals@entry=0x7ffff7394510) at Python/ceval.c:731
#37 0x00007ffff7e10366 in run_mod (mod=mod@entry=0x555555631e90, filename=filename@entry=0x7ffff72c2270, globals=globals@entry=0x7ffff7394510, locals=locals@entry=0x7ffff7394510, 
    flags=flags@entry=0x7fffffffdacc, arena=arena@entry=0x7ffff73ae2e8) at Python/pythonrun.c:1025
#38 0x00007ffff7e12984 in PyRun_FileExFlags (fp=fp@entry=0x55555555a4c0, filename_str=filename_str@entry=0x7ffff732f5f0 "reproduce.py", start=start@entry=257, globals=globals@entry=0x7ffff7394510, 
    locals=locals@entry=0x7ffff7394510, closeit=closeit@entry=1, flags=0x7fffffffdacc) at Python/pythonrun.c:978
#39 0x00007ffff7e12af2 in PyRun_SimpleFileExFlags (fp=fp@entry=0x55555555a4c0, filename=<optimized out>, closeit=closeit@entry=1, flags=flags@entry=0x7fffffffdacc) at Python/pythonrun.c:419
#40 0x00007ffff7e12fa4 in PyRun_AnyFileExFlags (fp=fp@entry=0x55555555a4c0, filename=<optimized out>, closeit=closeit@entry=1, flags=flags@entry=0x7fffffffdacc) at Python/pythonrun.c:81
#41 0x00007ffff7e29d52 in run_file (p_cf=0x7fffffffdacc, filename=0x55555555b6e0 L"reproduce.py", fp=0x55555555a4c0) at Modules/main.c:340
#42 Py_Main (argc=2, argv=0x5555555592a0) at Modules/main.c:810
#43 0x000055555555519f in main ()

@sdvillal
Copy link

sdvillal commented Nov 4, 2020

We are experiencing the same problem with shapely installed in conda environments both in macos and linux. Downgrading to geos 3.7.2 (which also brings shapely down from 1.7.1 to 1.6.4) also solves the problem for us.

Let me know if you want more example geometries to test.

@dr-jts
Copy link
Contributor

dr-jts commented Nov 4, 2020

Technically speaking the overlay operations (including intersection) don't work on GeometryCollections (even ones like this which are isomorphic to a valid MultiPolygon). This is why the case works after doing a union() - that turns the collection into a MultiPolygon.

I'm not sure why this isn't explicitly disallowed via an exception - this is what JTS does. (@strk do you know?)

That said, I'm surprised that this case doesn't in fact work, since it gets as far as it does.

In any case, the overlay code has been switched over to OverlayNG in GEOS 3.9, so this issue may disappear. Although the caution on not supporting GeometryCollections still stands.

@strk
Copy link
Member

strk commented Nov 4, 2020 via email

@pramsey
Copy link
Member

pramsey commented Nov 4, 2020

9 years ago. I'm embarrassed to say I have no recollection of why. There are some tickets around the issue with various clues.

https://trac.osgeo.org/geos/ticket/298
https://trac.osgeo.org/geos/ticket/725
https://lists.osgeo.org/pipermail/geos-devel/2009-January/003880.html

I remain intrigued by the possibility of actually stubbing out the cases where reasonable results can be returned. (Intersects, for example)

@dr-jts
Copy link
Contributor

dr-jts commented Nov 4, 2020

I remain intrigued by the possibility of actually stubbing out the cases where reasonable results can be returned. (Intersects, for example)

This issue is about overlay operations, don't forget. Predicates are another matter (and agreed, it would be nice to extend them to handle collections).

For overlay operations the safest path is to disallow GeometryCollection inputs. Although overlay code can handle some kinds of collection correctly (and OverlayNG is able to handle a wider range of collection inputs which are technically invalid MultiPolygons). But there's no easy way to tell a priori which are safe to run.

Collections which can be handled are homogeneous collections of line, points and non-overlapping polygons (since these are isomorphic to Multi geometries of that type). Handling hetergeneous collections might be possible via some chained operations, assuming a suitable semantic can be defined for each overlay operation (union and intersection are fairly obvious, the others might take some thought).

For now really Shapely should prohibit running overlay ops on collection inputs, using its native error mechanism. But it would be nice to have a guard in GEOS as well (especially since I'm not sure this restriction is well-documented).

@pramsey
Copy link
Member

pramsey commented Nov 4, 2020

I feel like union, intersection and difference are all doable using the chain of

  • first unaryunion each side to remove all duplicate areas
  • then pairwise apply the operation to each component

@dr-jts
Copy link
Contributor

dr-jts commented Nov 4, 2020

Actually JTS has a hack to compute a form of intersection for GeometryCollections. I don't like it, especially since as noted homogeneous collections can be handled in a more straightforward way.

I'm not sure about doing "pairwise apply" - that seems like it could potentially produce a lot of overlapping resultants. But I guess it depends on what sort of semantics are desired. To my mind the closer the semantics to "normal" overlay operations the better (i.e. the result should be non-overlapping and noded). One rationale is that this provides better chainability for overlay operations.

I guess you could make the argument that if the client wanted "normal" semantics it should supply a Multi geometry in the first place. But I think that a lot of clients won't understand this fine distinction, and just use GeometryCollections because their data is sloppily defined.

@dr-jts
Copy link
Contributor

dr-jts commented Nov 5, 2020

On 3.9 this test now works.

@dbaston
Copy link
Member

dbaston commented Feb 1, 2021

I think this can be closed then? Please re-open if I'm mistaken.

@dbaston dbaston closed this as completed Feb 1, 2021
@siodoni
Copy link

siodoni commented Sep 14, 2022

If I may, I would like to share how I solved this problem.

I solved this problem as follows:
First by finding geometries that were invalid.

select * from TABLE a where ST_IsValid(GEOM_COLUMN) = false;

After that, I corrected them, as follows:

update TABLE set GEOM_COLUMN = ST_CollectionExtract(ST_MakeValid(GEOM_COLUMN),3) where ST_IsValid(GEOM_COLUMN) = false;

After that I performed a reindex on the table.

reindex table CONCURRENTLY TABLE;

Hope this is helpful for someone with this kind of problem.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

7 participants