BUG: Fix STRtree memory leaks (#434)

docker/Dockerfile.valgrind

@@ -3,16 +3,31 @@
 # Then run the pytest suite with valgrind enabled:
 # docker run --rm pygeos/valgrind:latest valgrind --show-leak-kinds=definite --log-file=/tmp/valgrind-output python -m pytest -vv --valgrind --valgrind-log=/tmp/valgrind-output > valgrind.log
 
-FROM python:3.6-stretch
 
-RUN apt-get update && apt-get install -y valgrind libgeos-dev --no-install-recommends
+FROM python:3.9-slim-buster
 
-RUN pip install numpy Cython pytest pytest-valgrind
+RUN apt-get update && apt-get install -y build-essential valgrind curl --no-install-recommends
 
-COPY . /code
+RUN pip install cmake numpy Cython pytest pytest-valgrind
 
 WORKDIR /code
 
-RUN python setup.py build_ext --inplace
-
 ENV PYTHONMALLOC malloc
+
+ENV LD_LIBRARY_PATH /usr/local/lib
+
+# Build GEOS
+RUN export GEOS_VERSION=3.10.1 && \
+    mkdir /code/geos && cd /code/geos && \
+    curl -OL http://download.osgeo.org/geos/geos-$GEOS_VERSION.tar.bz2 && \
+    tar xfj geos-$GEOS_VERSION.tar.bz2 && \
+    cd geos-$GEOS_VERSION && mkdir build && cd build && \
+    cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_LIBDIR=lib .. && \
+    cmake --build . -j 4 && \
+    cmake --install . && \
+    cd /code
+
+COPY . /code
+
+# Build pygeos
+RUN python setup.py build_ext --inplace && python setup.py install

pygeos/tests/test_strtree.py

@@ -211,9 +211,70 @@ def test_query_unsupported_predicate(tree):
         tree.query(pygeos.points(1, 1), predicate="disjoint")
 
 
-def test_query_predicate_errors(tree):
+@pytest.mark.parametrize(
+    "predicate,expected",
+    [
+        ("intersects", [0, 1, 2]),
+        ("within", []),
+        ("contains", [1]),
+        ("overlaps", []),
+        ("crosses", []),
+        ("covers", [0, 1, 2]),
+        ("covered_by", []),
+        ("contains_properly", [1]),
+    ],
+)
+def test_query_prepared_inputs(tree, predicate, expected):
+    geom = box(0, 0, 2, 2)
+    pygeos.prepare(geom)
+    assert_array_equal(tree.query(geom, predicate=predicate), expected)
+
+
+@pytest.mark.parametrize(
+    "predicate",
+    [
+        pytest.param(
+            "intersects",
+            marks=pytest.mark.xfail(reason="intersects does not raise exception"),
+        ),
+        pytest.param(
+            "within",
+            marks=pytest.mark.xfail(
+                pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8"
+            ),
+        ),
+        pytest.param(
+            "contains",
+            marks=pytest.mark.xfail(
+                pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8"
+            ),
+        ),
+        "overlaps",
+        "crosses",
+        "touches",
+        pytest.param(
+            "covers",
+            marks=pytest.mark.xfail(
+                pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8"
+            ),
+        ),
+        pytest.param(
+            "covered_by",
+            marks=pytest.mark.xfail(
+                pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8"
+            ),
+        ),
+        pytest.param(
+            "contains_properly",
+            marks=pytest.mark.xfail(
+                pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8"
+            ),
+        ),
+    ],
+)
+def test_query_predicate_errors(tree, predicate):
     with pytest.raises(pygeos.GEOSException):
-        tree.query(pygeos.linestrings([1, 1], [1, float("nan")]), predicate="touches")
+        tree.query(pygeos.linestrings([1, 1], [1, float("nan")]), predicate=predicate)
 
 
 def test_query_tree_with_none():
@@ -1295,6 +1356,13 @@ def test_query_bulk_intersects_polygons(poly_tree, geometry, expected):
     assert_array_equal(poly_tree.query_bulk(geometry, predicate="intersects"), expected)
 
 
+def test_query_bulk_predicate_errors(tree):
+    with pytest.raises(pygeos.GEOSException):
+        tree.query_bulk(
+            [pygeos.linestrings([1, 1], [1, float("nan")])], predicate="touches"
+        )
+
+
 @pytest.mark.skipif(pygeos.geos_version >= (3, 10, 0), reason="GEOS >= 3.10")
 def test_query_bulk_dwithin_geos_version(tree):
     with pytest.raises(UnsupportedGEOSOperation, match="requires GEOS >= 3.10"):

src/strtree.c

@@ -209,6 +209,7 @@ static PyObject* STRtree_new(PyTypeObject* type, PyObject* args, PyObject* kwds)
       return NULL;
     }
     GEOSSTRtree_query_r(ctx, tree, dummy, dummy_query_callback, NULL);
+    GEOSGeom_destroy_r(ctx, dummy);
   }
 
   STRtreeObject* self = (STRtreeObject*)type->tp_alloc(type, 0);
@@ -272,6 +273,7 @@ static char evaluate_predicate(void* context, FuncGEOS_YpY_b* predicate_func,
                                GEOSGeometry* geom, GEOSPreparedGeometry* prepared_geom,
                                tree_geom_vec_t* in_geoms, tree_geom_vec_t* out_geoms,
                                npy_intp* count) {
+  char errstate = PGERR_SUCCESS;
   GeometryObject* pg_geom;
   GeometryObject** pg_geom_loc;  // address of geometry in tree geometries (_geoms)
   GEOSGeometry* target_geom;
@@ -306,7 +308,9 @@ static char evaluate_predicate(void* context, FuncGEOS_YpY_b* predicate_func,
     // keep the geometry if it passes the predicate
     predicate_result = predicate_func(context, prepared_geom_tmp, target_geom);
     if (predicate_result == 2) {
-      return PGERR_GEOS_EXCEPTION;
+      // error evaluating predicate; break and cleanup prepared geometry below
+      errstate = PGERR_GEOS_EXCEPTION;
+      break;
     } else if (predicate_result == 1) {
       kv_push(GeometryObject**, *out_geoms, pg_geom_loc);
       (*count)++;
@@ -319,7 +323,7 @@ static char evaluate_predicate(void* context, FuncGEOS_YpY_b* predicate_func,
     prepared_geom_tmp = NULL;
   }
 
-  return PGERR_SUCCESS;
+  return errstate;
 }
 
 /* Query the tree based on input geometry and predicate function.
@@ -528,8 +532,6 @@ static PyObject* STRtree_query_bulk(STRtreeObject* self, PyObject* args) {
 
       if (errstate != PGERR_SUCCESS) {
         kv_destroy(query_geoms);
-        kv_destroy(src_indexes);
-        kv_destroy(target_geoms);
         break;
       }
 
@@ -544,6 +546,8 @@ static PyObject* STRtree_query_bulk(STRtreeObject* self, PyObject* args) {
   GEOS_FINISH_THREADS;
 
   if (errstate != PGERR_SUCCESS) {
+    kv_destroy(src_indexes);
+    kv_destroy(target_geoms);
     return NULL;
   }
 
@@ -555,6 +559,8 @@ static PyObject* STRtree_query_bulk(STRtreeObject* self, PyObject* args) {
   result = (PyArrayObject*)PyArray_SimpleNew(2, dims, NPY_INTP);
   if (result == NULL) {
     PyErr_SetString(PyExc_RuntimeError, "could not allocate numpy array");
+    kv_destroy(src_indexes);
+    kv_destroy(target_geoms);
     return NULL;
   }
 
@@ -799,6 +805,8 @@ static PyObject* STRtree_nearest(STRtreeObject* self, PyObject* arr) {
   result = (PyArrayObject*)PyArray_SimpleNew(2, index_dims, NPY_INTP);
   if (result == NULL) {
     PyErr_SetString(PyExc_RuntimeError, "could not allocate numpy array");
+    kv_destroy(src_indexes);
+    kv_destroy(nearest_indexes);
     return NULL;
   }
 
@@ -934,19 +942,13 @@ static PyObject* STRtree_nearest_all(STRtreeObject* self, PyObject* args) {
       // if max_distance is defined, prescreen geometries using simple bbox expansion
       if (get_bounds(ctx, geom, &xmin, &ymin, &xmax, &ymax) == 0) {
         errstate = PGERR_GEOS_EXCEPTION;
-        kv_destroy(src_indexes);
-        kv_destroy(nearest_geoms);
-        kv_destroy(nearest_dist);
         break;
       }
 
       envelope = create_box(ctx, xmin - max_distance, ymin - max_distance,
                             xmax + max_distance, ymax + max_distance, 1);
       if (envelope == NULL) {
         errstate = PGERR_GEOS_EXCEPTION;
-        kv_destroy(src_indexes);
-        kv_destroy(nearest_geoms);
-        kv_destroy(nearest_dist);
         break;
       }
 
@@ -962,6 +964,11 @@ static PyObject* STRtree_nearest_all(STRtreeObject* self, PyObject* args) {
       }
     }
 
+    if (errstate == PGERR_GEOS_EXCEPTION) {
+      // break outer loop
+      break;
+    }
+
     // reset loop-dependent values of userdata
     kv_init(dist_pairs);
     userdata.min_distance = DBL_MAX;
@@ -973,9 +980,6 @@ static PyObject* STRtree_nearest_all(STRtreeObject* self, PyObject* args) {
     if (nearest_result == NULL) {
       errstate = PGERR_GEOS_EXCEPTION;
       kv_destroy(dist_pairs);
-      kv_destroy(src_indexes);
-      kv_destroy(nearest_geoms);
-      kv_destroy(nearest_dist);
       break;
     }
 
@@ -999,6 +1003,9 @@ static PyObject* STRtree_nearest_all(STRtreeObject* self, PyObject* args) {
   GEOS_FINISH_THREADS;
 
   if (errstate != PGERR_SUCCESS) {
+    kv_destroy(src_indexes);
+    kv_destroy(nearest_geoms);
+    kv_destroy(nearest_dist);
     return NULL;
   }
 
@@ -1009,6 +1016,9 @@ static PyObject* STRtree_nearest_all(STRtreeObject* self, PyObject* args) {
   result_indexes = (PyArrayObject*)PyArray_SimpleNew(2, index_dims, NPY_INTP);
   if (result_indexes == NULL) {
     PyErr_SetString(PyExc_RuntimeError, "could not allocate numpy array");
+    kv_destroy(src_indexes);
+    kv_destroy(nearest_geoms);
+    kv_destroy(nearest_dist);
     return NULL;
   }
 
@@ -1017,6 +1027,9 @@ static PyObject* STRtree_nearest_all(STRtreeObject* self, PyObject* args) {
   result_distances = (PyArrayObject*)PyArray_SimpleNew(1, distance_dims, NPY_DOUBLE);
   if (result_distances == NULL) {
     PyErr_SetString(PyExc_RuntimeError, "could not allocate numpy array");
+    kv_destroy(src_indexes);
+    kv_destroy(nearest_geoms);
+    kv_destroy(nearest_dist);
     return NULL;
   }