Repository: pygeos/pygeos

Issue:
Union_all skips NaG, while intersection_all does not
Compare:

```
line_1 = Geometry("LINESTRING(0 0, 2 2)")
line_2 = Geometry("LINESTRING(2 2, 3 3)")
union_all([[line_1, line_2, NaG]], axis=1).tolist()
[<pygeos.Geometry MULTILINESTRING ((0 0, 2 2), (2 2, 3 3))>]
intersection_all([[line_1, line_2, NaG]], axis=1).tolist()
[<pygeos.NaG>]
```

This is because union_all makes use of the UnaryUnion function from GEOS, while intersection_all uses the ufunc.reduce feature.

I guess the neatest would be adding a `skipna` kwarg.

Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 3e7be0d..4ad2275 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -11,6 +11,7 @@ Version 0.9 (unreleased)
 * Addition of a ``is_ccw()`` function for GEOS >= 3.7 (#201)
 * Addition of a ``minimum_clearance`` function for GEOS >= 3.6.0 (#223)
 * Addition of a ``offset_curve`` function (#229)
+* Addition of a ``relate_pattern`` function (#245)
 * Added support for pickling to ``Geometry`` objects (#190)
 * Limited the length of geometry repr to 80 characters (#189)
 * Argument in ``line_interpolate_point`` and ``line_locate_point``
@@ -26,10 +27,17 @@ Version 0.9 (unreleased)
 * API change: geometry and counting functions (``get_num_coordinates``,
   ``get_num_geometries``, ``get_num_interior_rings``, ``get_num_points``) now return 0
   for ``None`` input values instead of -1 (#218).
+* API change: ``intersection_all`` and ``symmetric_difference_all`` now ignore None values
+  instead of returning None if any value is None (#249).
+* API change: ``union_all`` now returns None (instead of ``GEOMETRYCOLLECTION EMPTY``) if
+  all input values are None (#249).
 * Fixed internal GEOS error code detection for ``get_dimensions`` and ``get_srid`` (#218).
 * Addition of ``prepare`` function that generates a GEOS prepared geometry which is stored on
   the Geometry object itself. All binary predicates (except ``equals``) make use of this (#92).
-
+* Use previously prepared geometries within ``STRtree`` ``query`` and ``query_bulk``
+  functions if available (#246).
+* Ensure that ``python setup.py clean`` removes all previously Cythonized and compiled
+  files (#239).
 
 Version 0.8 (2020-09-06)
 ------------------------
diff --git a/pygeos/_geometry.pyx b/pygeos/_geometry.pyx
index f045927..090f256 100644
--- a/pygeos/_geometry.pyx
+++ b/pygeos/_geometry.pyx
@@ -1,3 +1,5 @@
+# distutils: define_macros=GEOS_USE_ONLY_R_API
+
 from cpython cimport PyObject
 cimport cython
 
@@ -28,17 +30,13 @@ def get_parts(object[:] array):
     cdef Py_ssize_t part_idx = 0
     cdef Py_ssize_t idx = 0
     cdef GEOSGeometry *geom = NULL
-    cdef GEOSGeometry *part = NULL
+    cdef const GEOSGeometry *part = NULL
 
     counts = pygeos.get_num_geometries(array)
+    cdef Py_ssize_t count = counts.sum()
 
-    # None elements in array return -1 for count, so
-    # they must be filtered out before calculating total count
-    cdef Py_ssize_t count = counts[counts>0].sum()
-
-    if count <= 0:
-        # return immeidately if there are no geometries to return
-        # count is negative when the only entries in array are None
+    if count == 0:
+        # return immediately if there are no geometries to return
         return (
             np.empty(shape=(0, ), dtype=np.object),
             np.empty(shape=(0, ), dtype=np.intp)
@@ -70,7 +68,9 @@ def get_parts(object[:] array):
 
                 # clone the geometry to keep it separate from the inputs
                 part = GEOSGeom_clone_r(geos_handle, part)
-                parts_view[idx] = PyGEOS_CreateGeometry(part, geos_handle)
+                # cast part back to <GEOSGeometry> to discard const qualifier
+                # pending issue #227
+                parts_view[idx] = PyGEOS_CreateGeometry(<GEOSGeometry *>part, geos_handle)
 
                 idx += 1
 
diff --git a/pygeos/_geos.pyx b/pygeos/_geos.pyx
index 18138cb..2804ea7 100644
--- a/pygeos/_geos.pyx
+++ b/pygeos/_geos.pyx
@@ -1,3 +1,5 @@
+# distutils: define_macros=GEOS_USE_ONLY_R_API
+
 cdef class get_geos_handle:
     '''This class provides a context manager that wraps the GEOS context handle.
 
diff --git a/pygeos/predicates.py b/pygeos/predicates.py
index fc97532..cc8f721 100644
--- a/pygeos/predicates.py
+++ b/pygeos/predicates.py
@@ -28,6 +28,7 @@ __all__ = [
     "within",
     "equals_exact",
     "relate",
+    "relate_pattern",
 ]
 
 
@@ -57,7 +58,7 @@ def is_ccw(geometry, **kwargs):
     Note that there are no checks on whether lines are actually closed and
     not self-intersecting, while this is a requirement for is_ccw. The recommended
     usage of this function for linestrings is ``is_ccw(g) & is_simple(g)`` and for
-    linearrings ``is_ccw(g) & is_valid(g)``. 
+    linearrings ``is_ccw(g) & is_valid(g)``.
 
     Parameters
     ----------
@@ -768,3 +769,33 @@ def relate(a, b, **kwargs):
     'F0FFFF102'
     """
     return lib.relate(a, b, **kwargs)
+
+
+@multithreading_enabled
+def relate_pattern(a, b, pattern, **kwargs):
+    """
+    Returns True if the DE-9IM string code for the relationship between
+    the geometries satisfies the pattern, else False.
+
+    This function compares the DE-9IM code string for two geometries
+    against a specified pattern. If the string matches the pattern then
+    ``True`` is returned, otherwise ``False``. The pattern specified can
+    be an exact match (``0``, ``1`` or ``2``), a boolean match
+    (uppercase ``T`` or ``F``), or a wildcard (``*``). For example,
+    the pattern for the `within` predicate is ``'T*F**F***'``.
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+    pattern : string
+
+    Examples
+    --------
+    >>> point = Geometry("POINT (0.5 0.5)")
+    >>> square = Geometry("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))")
+    >>> relate(point, square)
+    '0FFFFF212'
+    >>> relate_pattern(point, square, "T*F**F***")
+    True
+    """
+    return lib.relate_pattern(a, b, pattern, **kwargs)
diff --git a/pygeos/set_operations.py b/pygeos/set_operations.py
index 0a335bd..eaad1ae 100644
--- a/pygeos/set_operations.py
+++ b/pygeos/set_operations.py
@@ -61,6 +61,9 @@ def intersection(a, b, **kwargs):
 def intersection_all(geometries, axis=0, **kwargs):
     """Returns the intersection of multiple geometries.
 
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None, None is returned.
+
     Parameters
     ----------
     geometries : array_like
@@ -81,7 +84,7 @@ def intersection_all(geometries, axis=0, **kwargs):
     >>> intersection_all([line_1, line_2])
     <pygeos.Geometry LINESTRING (1 1, 2 2)>
     >>> intersection_all([[line_1, line_2, None]], axis=1).tolist()
-    [None]
+    [<pygeos.Geometry LINESTRING (1 1, 2 2)>]
     """
     return lib.intersection.reduce(geometries, axis=axis, **kwargs)
 
@@ -111,6 +114,9 @@ def symmetric_difference(a, b, **kwargs):
 def symmetric_difference_all(geometries, axis=0, **kwargs):
     """Returns the symmetric difference of multiple geometries.
 
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None, None is returned.
+
     Parameters
     ----------
     geometries : array_like
@@ -131,7 +137,7 @@ def symmetric_difference_all(geometries, axis=0, **kwargs):
     >>> symmetric_difference_all([line_1, line_2])
     <pygeos.Geometry MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))>
     >>> symmetric_difference_all([[line_1, line_2, None]], axis=1).tolist()
-    [None]
+    [<pygeos.Geometry MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))>]
     """
     return lib.symmetric_difference.reduce(geometries, axis=axis, **kwargs)
 
@@ -162,6 +168,9 @@ def union(a, b, **kwargs):
 def union_all(geometries, axis=0, **kwargs):
     """Returns the union of multiple geometries.
 
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None, None is returned.
+
     Parameters
     ----------
     geometries : array_like
@@ -196,7 +205,15 @@ def union_all(geometries, axis=0, **kwargs):
         )
     # create_collection acts on the inner axis
     collections = lib.create_collection(geometries, GeometryType.GEOMETRYCOLLECTION)
-    return lib.unary_union(collections, **kwargs)
+    result = lib.unary_union(collections, **kwargs)
+    # for consistency with other _all functions, we replace GEOMETRY COLLECTION EMPTY
+    # if the original collection had no geometries
+    only_none = lib.get_num_geometries(collections) == 0
+    if np.isscalar(only_none):
+        return result if not only_none else None
+    else:
+        result[only_none] = None
+        return result
 
 
 @requires_geos("3.8.0")
diff --git a/setup.py b/setup.py
index 7b8e2a9..66335a1 100644
--- a/setup.py
+++ b/setup.py
@@ -1,5 +1,6 @@
 import builtins
 import os
+from pathlib import Path
 import subprocess
 import sys
 from distutils.version import LooseVersion
@@ -110,12 +111,14 @@ def get_geos_paths():
     }
 
 
-# Add numpy include dirs without importing numpy on module level.
-# derived from scikit-hep:
-# https://github.com/scikit-hep/root_numpy/pull/292
 class build_ext(_build_ext):
     def finalize_options(self):
         _build_ext.finalize_options(self)
+
+        # Add numpy include dirs without importing numpy on module level.
+        # derived from scikit-hep:
+        # https://github.com/scikit-hep/root_numpy/pull/292
+
         # Prevent numpy from thinking it is still in its setup process:
         try:
             del builtins.__NUMPY_SETUP__
@@ -129,7 +132,15 @@ class build_ext(_build_ext):
 
 ext_modules = []
 
-if "clean" not in sys.argv and "sdist" not in sys.argv:
+if "clean" in sys.argv:
+    # delete any previously Cythonized or compiled files in pygeos
+    p = Path("pygeos")
+    for pattern in ["*.c", "*.so", "*.pyd"]:
+        for filename in p.glob(pattern):
+            print("removing '{}'".format(filename))
+            filename.unlink()
+
+elif "sdist" not in sys.argv:
     ext_options = get_geos_paths()
 
     ext_modules = [
diff --git a/src/c_api.c b/src/c_api.c
index 59b3e69..00719fe 100644
--- a/src/c_api.c
+++ b/src/c_api.c
@@ -21,6 +21,6 @@ extern PyObject* PyGEOS_CreateGeometry(GEOSGeometry *ptr, GEOSContextHandle_t ct
     return GeometryObject_FromGEOS(ptr, ctx);
 }
 
-extern char PyGEOS_GetGEOSGeometry(GeometryObject *obj, GEOSGeometry **out) {
-    return get_geom(obj, out);
+extern char PyGEOS_GetGEOSGeometry(PyObject *obj, GEOSGeometry **out) {
+    return get_geom((GeometryObject*)obj, out);
 }
diff --git a/src/c_api.h b/src/c_api.h
index 61833a1..ced48d2 100644
--- a/src/c_api.h
+++ b/src/c_api.h
@@ -30,7 +30,7 @@
 /* char PyGEOS_GetGEOSGeometry(GeometryObject *obj, GEOSGeometry **out) */
 #define PyGEOS_GetGEOSGeometry_NUM 1
 #define PyGEOS_GetGEOSGeometry_RETURN char
-#define PyGEOS_GetGEOSGeometry_PROTO (GeometryObject * obj, GEOSGeometry * *out)
+#define PyGEOS_GetGEOSGeometry_PROTO (PyObject * obj, GEOSGeometry * *out)
 
 /* Total number of C API pointers */
 #define PyGEOS_API_num_pointers 2
diff --git a/src/geos.h b/src/geos.h
index 7fcfc97..1212de9 100644
--- a/src/geos.h
+++ b/src/geos.h
@@ -4,7 +4,9 @@
 #include <Python.h>
 
 /* To avoid accidental use of non reentrant GEOS API. */
+#ifndef GEOS_USE_ONLY_R_API
 #define GEOS_USE_ONLY_R_API
+#endif
 
 // wrap geos.h import to silence geos gcc warnings
 #pragma GCC diagnostic push
diff --git a/src/strtree.c b/src/strtree.c
index 932c137..8748b6d 100644
--- a/src/strtree.c
+++ b/src/strtree.c
@@ -210,29 +210,45 @@ void query_callback(void* item, void* user_data) {
  * Parameters
  * ----------
  * predicate_func: pointer to a prepared predicate function, e.g.,
- * GEOSPreparedIntersects_r geom: input geometry to prepare and test against each geometry
- * in the tree specified by in_indexes. tree_geometries: pointer to ndarray of all
- * geometries in the tree in_indexes: dynamic vector of indexes of tree geometries that
- * have overlapping envelopes with envelope of input geometry. out_indexes: dynamic vector
- * of indexes of tree geometries that meet predicate function. count: pointer to an
- * integer where the number of geometries that met the predicate will be written
+ *   GEOSPreparedIntersects_r
+ *
+ * geom: input geometry to prepare and test against each geometry in the tree specified by
+ *   in_indexes.
+ *
+ * prepared_geom: input prepared geometry, only if previously created.  If NULL, geom
+ *   will be prepared instead.
+ *
+ * tree_geometries: pointer to ndarray of all geometries in the tree
+ *
+ * in_indexes: dynamic vector of indexes of tree geometries that have overlapping
+ *   envelopes with envelope of input geometry.
+ *
+ * out_indexes: dynamic vector of indexes of tree geometries that meet predicate function.
+ *
+ * count: pointer to an integer where the number of geometries that met the predicate will
+ *   be written.
  *
  * Returns PGERR_GEOS_EXCEPTION if an error was encountered or PGERR_SUCCESS otherwise
  * */
 
 static char evaluate_predicate(void* context, FuncGEOS_YpY_b* predicate_func,
-                               GEOSGeometry* geom, pg_geom_obj_vec* tree_geometries,
-                               npy_intp_vec* in_indexes, npy_intp_vec* out_indexes,
-                               npy_intp* count) {
+                               GEOSGeometry* geom, GEOSPreparedGeometry* prepared_geom,
+                               pg_geom_obj_vec* tree_geometries, npy_intp_vec* in_indexes,
+                               npy_intp_vec* out_indexes, npy_intp* count) {
   GeometryObject* pg_geom;
   GEOSGeometry* target_geom;
-  const GEOSPreparedGeometry* prepared_geom;
+  const GEOSPreparedGeometry* prepared_geom_tmp;
   npy_intp i, size, index;
 
-  // Create prepared geometry
-  prepared_geom = GEOSPrepare_r(context, geom);
   if (prepared_geom == NULL) {
-    return PGERR_GEOS_EXCEPTION;
+    // geom was not previously prepared, prepare it now
+    prepared_geom_tmp = GEOSPrepare_r(context, geom);
+    if (prepared_geom_tmp == NULL) {
+      return PGERR_GEOS_EXCEPTION;
+    }
+  } else {
+    // cast to const only needed until larger refactor of all geom pointers to const
+    prepared_geom_tmp = (const GEOSPreparedGeometry*)prepared_geom;
   }
 
   size = kv_size(*in_indexes);
@@ -249,13 +265,17 @@ static char evaluate_predicate(void* context, FuncGEOS_YpY_b* predicate_func,
     get_geom((GeometryObject*)pg_geom, &target_geom);
 
     // keep the index value if it passes the predicate
-    if (predicate_func(context, prepared_geom, target_geom)) {
+    if (predicate_func(context, prepared_geom_tmp, target_geom)) {
       kv_push(npy_intp, *out_indexes, index);
       (*count)++;
     }
   }
 
-  GEOSPreparedGeom_destroy_r(context, prepared_geom);
+  if (prepared_geom == NULL) {
+    // only if we created prepared_geom_tmp here, destroy it
+    GEOSPreparedGeom_destroy_r(context, prepared_geom_tmp);
+    prepared_geom_tmp = NULL;
+  }
 
   return PGERR_SUCCESS;
 }
@@ -272,9 +292,10 @@ static char evaluate_predicate(void* context, FuncGEOS_YpY_b* predicate_func,
  * */
 
 static PyObject* STRtree_query(STRtreeObject* self, PyObject* args) {
-  GeometryObject* geometry;
+  GeometryObject* geometry = NULL;
   int predicate_id = 0;  // default no predicate
-  GEOSGeometry* geom;
+  GEOSGeometry* geom = NULL;
+  GEOSPreparedGeometry* prepared_geom = NULL;
   npy_intp_vec query_indexes,
       predicate_indexes;  // Resizable array for matches for each geometry
   npy_intp count;
@@ -290,7 +311,7 @@ static PyObject* STRtree_query(STRtreeObject* self, PyObject* args) {
     return NULL;
   }
 
-  if (!get_geom(geometry, &geom)) {
+  if (!get_geom_with_prepared(geometry, &geom, &prepared_geom)) {
     PyErr_SetString(PyExc_TypeError, "Invalid geometry");
     return NULL;
   }
@@ -326,8 +347,8 @@ static PyObject* STRtree_query(STRtreeObject* self, PyObject* args) {
   }
 
   kv_init(predicate_indexes);
-  errstate = evaluate_predicate(ctx, predicate_func, geom, &self->_geoms, &query_indexes,
-                                &predicate_indexes, &count);
+  errstate = evaluate_predicate(ctx, predicate_func, geom, prepared_geom, &self->_geoms,
+                                &query_indexes, &predicate_indexes, &count);
   if (errstate != PGERR_SUCCESS) {
     // error performing predicate
     kv_destroy(query_indexes);
@@ -361,9 +382,10 @@ static PyObject* STRtree_query(STRtreeObject* self, PyObject* args) {
 static PyObject* STRtree_query_bulk(STRtreeObject* self, PyObject* args) {
   PyObject* arr;
   PyArrayObject* pg_geoms;
-  GeometryObject* pg_geom;
+  GeometryObject* pg_geom = NULL;
   int predicate_id = 0;  // default no predicate
-  GEOSGeometry* geom;
+  GEOSGeometry* geom = NULL;
+  GEOSPreparedGeometry* prepared_geom = NULL;
   npy_intp_vec query_indexes, src_indexes, target_indexes;
   npy_intp i, j, n, size;
   FuncGEOS_YpY_b* predicate_func = NULL;
@@ -416,7 +438,7 @@ static PyObject* STRtree_query_bulk(STRtreeObject* self, PyObject* args) {
   for (i = 0; i < n; i++) {
     // get pygeos geometry from input geometry array
     pg_geom = *(GeometryObject**)PyArray_GETPTR1(pg_geoms, i);
-    if (!get_geom(pg_geom, &geom)) {
+    if (!get_geom_with_prepared(pg_geom, &geom, &prepared_geom)) {
       errstate = PGERR_NOT_A_GEOMETRY;
       break;
     }
@@ -442,8 +464,9 @@ static PyObject* STRtree_query_bulk(STRtreeObject* self, PyObject* args) {
       }
     } else {
       // this pushes directly onto target_indexes
-      errstate = evaluate_predicate(ctx, predicate_func, geom, &self->_geoms,
-                                    &query_indexes, &target_indexes, &size);
+      errstate =
+          evaluate_predicate(ctx, predicate_func, geom, prepared_geom, &self->_geoms,
+                             &query_indexes, &target_indexes, &size);
 
       if (errstate != PGERR_SUCCESS) {
         kv_destroy(query_indexes);
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 2bc3b05..f820d87 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -717,46 +717,57 @@ static void YY_Y_func_reduce(char** args, npy_intp* dimensions, npy_intp* steps,
   FuncGEOS_YY_Y* func = (FuncGEOS_YY_Y*)data;
   GEOSGeometry *in1 = NULL, *in2 = NULL, *out = NULL;
 
+  // Whether to destroy a temporary intermediate value of `out`:
+  char do_destroy = 0;
+
   GEOS_INIT_THREADS;
 
   if (!get_geom(*(GeometryObject**)args[0], &out)) {
     errstate = PGERR_NOT_A_GEOMETRY;
-  } else if (out != NULL) {
+  } else {
     BINARY_LOOP {
-      // Cleanup previous in1
-      if (i > 1) {
-        // If i == 0, in1 is undefined
-        // If i == 1, in1 is the first input, which is owned by python
-        GEOSGeom_destroy_r(ctx, in1);
-      }
-      // This is the main reduce logic: in1 becomes previous out
+      // Get the geometry inputs; in1 from previous iteration, in2 from array
       in1 = out;
-      // Get the other geometry (as normal)
       if (!get_geom(*(GeometryObject**)ip2, &in2)) {
         errstate = PGERR_NOT_A_GEOMETRY;
         break;
       }
-      if (in2 == NULL) {
-        // We break the loop as the answer will remain NULL (None)
-        out = NULL;
-        break;
-      } else {
+
+      /* Either (or both) in1 and in2 could be NULL (Python: None).
+       * Reduction operations should skip None values. We have 4 possible combinations:
+       */
+
+      // 1. (not NULL, not NULL); run the GEOS function
+      if ((in1 != NULL) && (in2 != NULL)) {
         out = func(ctx, in1, in2);
+  
+        // Discard in1 if it was a temporary intermediate
+        if (do_destroy) {
+          GEOSGeom_destroy_r(ctx, in1);
+        }
+
+        // Mark the newly generated geometry as intermediate. Note: out will become in1. 
+        do_destroy = 1;
+
+        // Break on error (we do this after discarding in1 to avoid memleaks)
         if (out == NULL) {
           errstate = PGERR_GEOS_EXCEPTION;
           break;
         }
       }
-    }
-    // We need to cleanup the intermediate geometry stored in in1.
-    // However, in some cases the intermediate geometry equals the first input, which is
-    // 'owned' by python. Destoying that would lead to a segfault when the python object
-    // is dereferenced.
-    // We check for that situation explicitly from the ufunc input (misusing the variable
-    // 'in2').
-    get_geom(*(GeometryObject**)args[0], &in2);
-    if (in1 != in2) {
-      GEOSGeom_destroy_r(ctx, in1);
+
+      // 2. (NULL, not NULL); When the first element of the reduction axis is None
+      else if ((in1 == NULL) && (in2 != NULL)) {
+        // Keep in2 as 'outcome' of the operation.
+        out = in2;
+        // Ensure that it will not be destroyed (it is owned by python)
+        do_destroy = 0;
+      }
+
+      // 3. (not NULL, NULL); When a None value is encountered after a not-None
+      //    Don't do `out = in1`, as that is already the case.
+      // 4. (NULL, NULL); When we have not yet encountered any not-None
+      //    Do nothing; out will remain NULL
     }
   }
 
@@ -1535,6 +1546,66 @@ finish:
 }
 static PyUFuncGenericFunction relate_funcs[1] = {&relate_func};
 
+static char relate_pattern_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT, NPY_BOOL};
+static void relate_pattern_func(char** args, npy_intp* dimensions, npy_intp* steps,
+                                void* data) {
+  GEOSGeometry *in1 = NULL, *in2 = NULL;
+  const char* pattern = NULL;
+  npy_bool ret;
+
+  /* get the pattern argument (only deal with scalar for now) */
+  char* ip3 = args[2];
+  npy_intp is3 = steps[2];
+
+  if (is3 != 0) {
+    PyErr_Format(PyExc_ValueError, "pattern keyword only supports scalar argument");
+    return;
+  }
+  PyObject* in3 = *(PyObject**)ip3;
+  if (PyUnicode_Check(in3)) {
+    pattern = PyUnicode_AsUTF8(in3);
+    if (pattern == NULL) {
+      /* error happened in PyUnicode_AsUTF8, error already set by Python */
+      return;
+    }
+  } else {
+    PyErr_Format(PyExc_TypeError, "pattern keyword expected string, got %s",
+                 Py_TYPE(in3)->tp_name);
+    return;
+  }
+
+  GEOS_INIT_THREADS;
+
+  TERNARY_LOOP {
+    /* get the geometries: return on error */
+    if (!get_geom(*(GeometryObject**)ip1, &in1)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      goto finish;
+    }
+    if (!get_geom(*(GeometryObject**)ip2, &in2)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      goto finish;
+    }
+    /* ip3 is already handled above */
+
+    if ((in1 == NULL) | (in2 == NULL)) {
+      /* in case of a missing value: return 0 (False) */
+      ret = 0;
+    } else {
+      ret = GEOSRelatePattern_r(ctx, in1, in2, pattern);
+      if (ret == 2) {
+        errstate = PGERR_GEOS_EXCEPTION;
+        goto finish;
+      }
+    }
+    *(npy_bool*)op1 = ret;
+  }
+
+finish:
+  GEOS_FINISH_THREADS;
+}
+static PyUFuncGenericFunction relate_pattern_funcs[1] = {&relate_pattern_func};
+
 /* define double -> geometry construction functions */
 static char points_dtypes[2] = {NPY_DOUBLE, NPY_OBJECT};
 static void points_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
@@ -2387,6 +2458,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_CUSTOM(voronoi_polygons, 4);
   DEFINE_CUSTOM(is_valid_reason, 1);
   DEFINE_CUSTOM(relate, 2);
+  DEFINE_CUSTOM(relate_pattern, 3);
   DEFINE_GENERALIZED(points, 1, "(d)->()");
   DEFINE_GENERALIZED(linestrings, 1, "(i, d)->()");
   DEFINE_GENERALIZED(linearrings, 1, "(i, d)->()");


Test patch:
diff --git a/pygeos/test/test_predicates.py b/pygeos/test/test_predicates.py
index 2148282..4093881 100644
--- a/pygeos/test/test_predicates.py
+++ b/pygeos/test/test_predicates.py
@@ -3,7 +3,7 @@ import pygeos
 from pygeos import Geometry
 import numpy as np
 
-from .common import point, all_types, polygon, geometry_collection
+from .common import point, all_types, line_string, polygon, geometry_collection, empty
 
 UNARY_PREDICATES = (
     pygeos.is_empty,
@@ -111,6 +111,43 @@ def test_relate_none(g1, g2):
     assert pygeos.relate(g1, g2) is None
 
 
+def test_relate_pattern():
+    line_string = pygeos.linestrings([(0, 0), (1, 0), (1, 1)])
+    polygon = pygeos.box(0, 0, 2, 2)
+    assert pygeos.relate(line_string, polygon) == "11F00F212"
+    assert pygeos.relate_pattern(line_string, polygon, "11F00F212")
+    assert pygeos.relate_pattern(line_string, polygon, "*********")
+    assert not pygeos.relate_pattern(line_string, polygon, "F********")
+
+
+def test_relate_pattern_empty():
+    assert pygeos.relate_pattern(empty, empty, "*" * 9).item() is True
+
+
+@pytest.mark.parametrize("g1, g2", [(point, None), (None, point), (None, None)])
+def test_relate_pattern_none(g1, g2):
+    assert pygeos.relate_pattern(g1, g2, "*" * 9).item() is False
+
+
+def test_relate_pattern_incorrect_length():
+    with pytest.raises(pygeos.GEOSException, match="Should be length 9"):
+        pygeos.relate_pattern(point, polygon, "**")
+
+    with pytest.raises(pygeos.GEOSException, match="Should be length 9"):
+        pygeos.relate_pattern(point, polygon, "**********")
+
+
+@pytest.mark.parametrize("pattern", [b"*********", 10, None])
+def test_relate_pattern_non_string(pattern):
+    with pytest.raises(TypeError, match="expected string"):
+        pygeos.relate_pattern(point, polygon, pattern)
+
+
+def test_relate_pattern_non_scalar():
+    with pytest.raises(ValueError, match="only supports scalar"):
+        pygeos.relate_pattern([point] *2, polygon, ["*********"] * 2)
+
+
 @pytest.mark.skipif(pygeos.geos_version < (3, 7, 0), reason="GEOS < 3.7")
 @pytest.mark.parametrize("geom, expected", [
     (Geometry("LINEARRING (0 0, 0 1, 1 1, 0 0)"), False),
diff --git a/pygeos/test/test_set_operations.py b/pygeos/test/test_set_operations.py
index 886a50a..097937c 100644
--- a/pygeos/test/test_set_operations.py
+++ b/pygeos/test/test_set_operations.py
@@ -65,6 +65,46 @@ def test_set_operation_reduce_axis(func, related_func):
     assert actual.shape == (3,)
 
 
+@pytest.mark.parametrize("none_position", range(3))
+@pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
+def test_set_operation_reduce_one_none(func, related_func, none_position):
+    # API change: before, intersection_all and symmetric_difference_all returned
+    # None if any input geometry was None.
+    # The new behaviour is to ignore None values.
+    test_data = reduce_test_data[:2]
+    test_data.insert(none_position, None)
+    actual = func(test_data)
+    expected = related_func(reduce_test_data[0], reduce_test_data[1])
+    assert pygeos.equals(actual, expected)
+
+
+@pytest.mark.parametrize("none_position", range(3))
+@pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
+def test_set_operation_reduce_two_none(func, related_func, none_position):
+    test_data = reduce_test_data[:2]
+    test_data.insert(none_position, None)
+    test_data.insert(none_position, None)
+    actual = func(test_data)
+    expected = related_func(reduce_test_data[0], reduce_test_data[1])
+    assert pygeos.equals(actual, expected)
+
+
+@pytest.mark.parametrize("n", range(1, 3))
+@pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
+def test_set_operation_reduce_all_none(n, func, related_func):
+    # API change: before, union_all([None]) yielded EMPTY GEOMETRYCOLLECTION
+    # The new behaviour is that it returns None if all inputs are None.
+    assert func([None] * n) is None
+
+
+@pytest.mark.parametrize("n", range(1, 3))
+@pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
+def test_set_operation_reduce_all_none_arr(n, func, related_func):
+    # API change: before, union_all([None]) yielded EMPTY GEOMETRYCOLLECTION
+    # The new behaviour is that it returns None if all inputs are None.
+    assert func([[None] * n] * 2, axis=1).tolist() == [None, None]
+
+
 @pytest.mark.skipif(pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8")
 @pytest.mark.parametrize("n", range(1, 4))
 def test_coverage_union_reduce_1dim(n):
@@ -73,11 +113,7 @@ def test_coverage_union_reduce_1dim(n):
       1. It expects only non-overlapping polygons
       2. It expects GEOS 3.8.0+
     """
-    test_data = [
-        pygeos.box(0, 0, 1, 1),
-        pygeos.box(1, 0, 2, 1),
-        pygeos.box(2, 0, 3, 1),
-    ]
+    test_data = [pygeos.box(0, 0, 1, 1), pygeos.box(1, 0, 2, 1), pygeos.box(2, 0, 3, 1)]
     actual = pygeos.coverage_union_all(test_data[:n])
     # perform the reduction in a python loop and compare
     expected = test_data[0]
@@ -105,18 +141,20 @@ def test_coverage_union_overlapping_inputs():
     polygon = Geometry("POLYGON ((1 1, 1 0, 0 0, 0 1, 1 1))")
 
     # Overlapping polygons raise an error
-    with pytest.raises(pygeos.GEOSException, match="CoverageUnion cannot process incorrectly noded inputs."):
-        pygeos.coverage_union(polygon, Geometry("POLYGON ((1 0, 0.9 1, 2 1, 2 0, 1 0))"))
+    with pytest.raises(
+        pygeos.GEOSException,
+        match="CoverageUnion cannot process incorrectly noded inputs.",
+    ):
+        pygeos.coverage_union(
+            polygon, Geometry("POLYGON ((1 0, 0.9 1, 2 1, 2 0, 1 0))")
+        )
 
 
 @pytest.mark.skipif(pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8")
 @pytest.mark.parametrize(
     "geom_1, geom_2",
     # All possible polygon, non_polygon combinations
-    [
-        [polygon, non_polygon]
-        for non_polygon in non_polygon_types
-    ] +
+    [[polygon, non_polygon] for non_polygon in non_polygon_types] +
     # All possible non_polygon, non_polygon combinations
     [
         [non_polygon_1, non_polygon_2]
@@ -126,5 +164,7 @@ def test_coverage_union_overlapping_inputs():
 )
 def test_coverage_union_non_polygon_inputs(geom_1, geom_2):
     # Non polygon geometries raise an error
-    with pytest.raises(pygeos.GEOSException, match="Unhandled geometry type in CoverageUnion."):
+    with pytest.raises(
+        pygeos.GEOSException, match="Unhandled geometry type in CoverageUnion."
+    ):
         pygeos.coverage_union(geom_1, geom_2)
diff --git a/pygeos/test/test_strtree.py b/pygeos/test/test_strtree.py
index b59739b..1677f0b 100644
--- a/pygeos/test/test_strtree.py
+++ b/pygeos/test/test_strtree.py
@@ -217,6 +217,15 @@ def test_query_tree_with_none():
 
 ### predicate == 'intersects'
 
+
+def test_query_with_prepared(tree):
+    geom = box(0, 0, 1, 1)
+    expected = tree.query(geom, predicate="intersects")
+
+    pygeos.prepare(geom)
+    assert_array_equal(expected, tree.query(geom, predicate="intersects"))
+
+
 # TEMPORARY xfail: MultiPoint intersects with prepared geometries does not work
 # properly on GEOS 3.5.x; it was fixed in 3.6+
 @pytest.mark.parametrize(
@@ -1045,6 +1054,20 @@ def test_query_bulk_invalid_predicate(tree):
 
 ### predicate == 'intersects'
 
+
+def test_query_bulk_with_prepared(tree):
+    geom = np.array([box(0, 0, 1, 1), box(3, 3, 5, 5)])
+    expected = tree.query_bulk(geom, predicate="intersects")
+
+    # test with array of partially prepared geometries
+    pygeos.prepare(geom[0])
+    assert_array_equal(expected, tree.query_bulk(geom, predicate="intersects"))
+
+    # test with fully prepared geometries
+    pygeos.prepare(geom)
+    assert_array_equal(expected, tree.query_bulk(geom, predicate="intersects"))
+
+
 # TEMPORARY xfail: MultiPoint intersects with prepared geometries does not work
 # properly on GEOS 3.5.x; it was fixed in 3.6+
 @pytest.mark.parametrize(