Repository: shapely/shapely (language: python)

Issue:
Segfault when cleaning up result from specific intersection_all case
With the following python script:

```python
import shapely
from shapely import LineString
line1 = LineString([(0, 0), (2, 2)])

print(shapely.intersection_all([line1, None, None]))
```

I see a segfault when existing the interpreter or when finishing the script, so on clean-up of the python objects I suppose. 

When running the script under gdb, I see:

```
$ gdb --args python test_shapely_segfault.py 
...
(gdb) run
Starting program: ...
...
LINESTRING (0 0, 2 2)

Thread 1 "python" received signal SIGSEGV, Segmentation fault.
0x0000555555944010 in ?? ()
(gdb) bt
#0  0x0000555555944010 in ?? ()
#1  0x00007ffff732ba1e in GeometryObject_dealloc (self=0x7fffb81a3e80) at src/pygeom.c:48
#2  0x00005555556b0909 in subtype_dealloc (self=<LineString at remote 0x7fffb81a3e80>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/typeobject.c:1351
#3  0x0000555555677e38 in _Py_Dealloc (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/object.c:2209
#4  _Py_DECREF (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Include/object.h:430
#5  _Py_XDECREF (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Include/object.h:497
#6  free_keys_object (keys=0x555555ded7e0) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/dictobject.c:598
#7  dictkeys_decref (dk=0x555555ded7e0) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/dictobject.c:333
#8  dict_dealloc (mp=0x7ffff7596dc0) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/dictobject.c:2026
#9  0x00005555557bd0c0 in _Py_DECREF (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Include/object.h:430
#10 _Py_XDECREF (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Include/object.h:497
#11 module_dealloc (m=0x7ffff759e680) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/moduleobject.c:677
#12 0x00005555556ec90f in _Py_Dealloc (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/object.c:2209
#13 _Py_DECREF (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Include/object.h:430
#14 _Py_XDECREF (op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Include/object.h:497
#15 insertdict (value=<optimized out>, hash=<optimized out>, key='__main__', mp=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/dictobject.c:1123
#16 PyDict_SetItem (value=<optimized out>, key='__main__', op=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/dictobject.c:1573
#17 dict_ass_sub (mp=<optimized out>, v='__main__', w=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Objects/dictobject.c:2173
#18 0x00005555557e0d0c in _PyImport_Cleanup (tstate=0x555555946ba0) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Python/import.c:526
#19 0x00005555557e1106 in Py_FinalizeEx () at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Python/pylifecycle.c:1429
#20 0x00005555557e6e6b in Py_RunMain () at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Modules/main.c:683
#21 0x00005555557e7239 in Py_BytesMain (argc=<optimized out>, argv=<optimized out>) at /home/conda/feedstock_root/build_artifacts/python-split_1639967286578/work/Modules/main.c:1094
#22 0x00007ffff7c59083 in __libc_start_main () from /lib/x86_64-linux-gnu/libc.so.6
#23 0x0000555555760141 in _start ()
(gdb) q
```

(from checking some specific examples in https://github.com/shapely/shapely/pull/1501)

cc @caspervdw 


Discussion / hints:
I am not sure how stable / correct this backtrace is, but with some additional operations, I got initially a different backtrace:

```
import shapely
from shapely import LineString
line1 = LineString([(0, 0), (2, 2)])

print(shapely.intersection_all([line1, None, None]))
print(shapely.intersection_all([line1, None]))
print(shapely.intersection_all([line1]))
```
gives

```
Thread 1 "python" received signal SIGSEGV, Segmentation fault.
0x0000000000000000 in ?? ()
(gdb) bt
#0  0x0000000000000000 in ?? ()
#1  0x00007ffff72fa4dc in GEOSisEmpty_r () from /home/joris/miniconda3/envs/shapely-dev/lib/libgeos_c.so.1
#2  0x00007ffff732a636 in wkt_empty_3d_geometry (ctx=ctx@entry=0x555555defdf0, geom=0x555555dcde70, wkt=wkt@entry=0x7fffffffa698) at src/geos.c:316
#3  0x00007ffff73371bd in to_wkt_func (args=<optimized out>, dimensions=<optimized out>, steps=<optimized out>, data=<optimized out>) at src/ufuncs.c:2967
#4  0x00007ffff6cedaa5 in generic_wrapped_legacy_loop () from /home/joris/miniconda3/envs/shapely-dev/lib/python3.9/site-packages/numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so
#5  0x00007ffff6ceffee in try_trivial_single_output_loop () from /home/joris/miniconda3/envs/shapely-dev/lib/python3.9/site-packages/numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so
#6  0x00007ffff6cf6b94 in ufunc_generic_fastcall () from /home/joris/miniconda3/envs/shapely-dev/lib/python3.9/site-packages/numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so
#7  0x00005555556b4a5b in PyVectorcall_Call (kwargs=<optimized out>, 
```
Looking at the reduce C implementation, I think that this happens always when you reduce exactly 1 geometry.

In that case the in and out geometries are different PyObjects that share the same GEOS pointer. Can you confirm that?

Fix patch:
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 396da5a03..e32f2c567 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -756,7 +756,7 @@ static void YY_Y_func_reduce(char** args, npy_intp* dimensions, npy_intp* steps,
   GEOSGeometry *in1 = NULL, *in2 = NULL, *out = NULL;
 
   // Whether to destroy a temporary intermediate value of `out`:
-  char do_destroy = 0;
+  char out_ownership = 0;
 
   GEOS_INIT_THREADS;
 
@@ -780,12 +780,12 @@ static void YY_Y_func_reduce(char** args, npy_intp* dimensions, npy_intp* steps,
         out = func(ctx, in1, in2);
 
         // Discard in1 if it was a temporary intermediate
-        if (do_destroy) {
+        if (out_ownership) {
           GEOSGeom_destroy_r(ctx, in1);
         }
 
         // Mark the newly generated geometry as intermediate. Note: out will become in1.
-        do_destroy = 1;
+        out_ownership = 1;
 
         // Break on error (we do this after discarding in1 to avoid memleaks)
         if (out == NULL) {
@@ -799,7 +799,7 @@ static void YY_Y_func_reduce(char** args, npy_intp* dimensions, npy_intp* steps,
         // Keep in2 as 'outcome' of the operation.
         out = in2;
         // Ensure that it will not be destroyed (it is owned by python)
-        do_destroy = 0;
+        out_ownership = 0;
       }
 
       // 3. (not NULL, NULL); When a None value is encountered after a not-None
@@ -809,6 +809,15 @@ static void YY_Y_func_reduce(char** args, npy_intp* dimensions, npy_intp* steps,
     }
   }
 
+  // In case we do not own the output, make a clone (else we end up with 2 PyObjects
+  // referencing the same GEOS Geometry)
+  if ((errstate == PGERR_SUCCESS) && (!out_ownership)) {
+    out = GEOSGeom_clone_r(ctx, out);
+    if (out == NULL) {
+      errstate = PGERR_GEOS_EXCEPTION;
+    }
+  }
+
   GEOS_FINISH_THREADS;
 
   // fill the numpy array with a single PyObject while holding the GIL


Test patch:
diff --git a/shapely/tests/test_set_operations.py b/shapely/tests/test_set_operations.py
index 5e64f3419..d143d435e 100644
--- a/shapely/tests/test_set_operations.py
+++ b/shapely/tests/test_set_operations.py
@@ -100,6 +100,13 @@ def test_set_operation_reduce_1dim(n, func, related_func):
     assert shapely.equals(actual, expected)
 
 
+@pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
+def test_set_operation_reduce_single_geom(func, related_func):
+    geom = shapely.Point(1, 1)
+    actual = func([geom, None, None])
+    assert shapely.equals(actual, geom)
+
+
 @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
 def test_set_operation_reduce_axis(func, related_func):
     data = [[point] * 2] * 3  # shape = (3, 2)