ENH: prepared geometry as additional (cached) attribute on GeometryObject

src/pygeom.c

@@ -19,6 +19,7 @@ PyObject *GeometryObject_FromGEOS(PyTypeObject *type, GEOSGeometry *ptr)
         return NULL;
     } else {
         self->ptr = ptr;
+        self->ptr_prepared = NULL;
         return (PyObject *) self;
     }
 }
@@ -29,11 +30,15 @@ static void GeometryObject_dealloc(GeometryObject *self)
     if (self->ptr != NULL) {
         GEOSGeom_destroy_r(context_handle, self->ptr);
     }
+    if (self->ptr_prepared != NULL) {
+        GEOSPreparedGeom_destroy_r(context_handle, self->ptr_prepared);
+    }
     Py_TYPE(self)->tp_free((PyObject *) self);
 }
 
 static PyMemberDef GeometryObject_members[] = {
     {"_ptr", T_PYSSIZET, offsetof(GeometryObject, ptr), READONLY, "pointer to GEOSGeometry"},
+    {"_ptr_prepared", T_PYSSIZET, offsetof(GeometryObject, ptr_prepared), READONLY, "pointer to PreparedGEOSGeometry"},
     {NULL}  /* Sentinel */
 };
 
@@ -164,6 +169,24 @@ char get_geom(GeometryObject *obj, GEOSGeometry **out) {
     }
 }
 
+/* Get a GEOSPreparedGeometry pointer from a GeometryObject, or NULL if the input is
+Py_None. Returns 0 on error, 1 on success. */
+char get_geom_prepared(GeometryObject *obj, GEOSPreparedGeometry **out) {
+    if (!PyObject_IsInstance((PyObject *) obj, (PyObject *) &GeometryType)) {
+        if ((PyObject *) obj == Py_None) {
+            *out = NULL;
+            return 1;
+        } else {
+            PyErr_Format(PyExc_TypeError, "One of the arguments is of incorrect type. Please provide only Geometry objects.");
+            return 0;
+        }
+    } else {
+        *out = obj->ptr_prepared;
+        return 1;
+    }
+}
+
+
 int
 init_geom_type(PyObject *m)
 {

src/pygeom.h

@@ -8,6 +8,7 @@
 typedef struct {
     PyObject_HEAD
     void *ptr;
+    void *ptr_prepared;
 } GeometryObject;
 
 
@@ -17,6 +18,7 @@ extern PyTypeObject GeometryType;
 extern PyObject *GeometryObject_FromGEOS(PyTypeObject *type, GEOSGeometry *ptr);
 /* Get a GEOSGeometry from a GeometryObject */
 extern char get_geom(GeometryObject *obj, GEOSGeometry **out);
+extern char get_geom_prepared(GeometryObject *obj, GEOSPreparedGeometry **out);
 
 extern int init_geom_type(PyObject *m);
 

src/ufuncs.c

@@ -161,6 +161,44 @@ static void YY_b_func(char **args, npy_intp *dimensions,
 }
 static PyUFuncGenericFunction YY_b_funcs[1] = {&YY_b_func};
 
+
+/* Define the geom, geom -> bool functions (YY_b) prepared */
+static void *intersects_prepared_data[1] = {GEOSIntersects_r};
+typedef char FuncGEOS_YY_b_p(void *context, void *a, void *b);
+static char YY_b_p_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_BOOL};
+static void YY_b_p_func(char **args, npy_intp *dimensions,
+                      npy_intp *steps, void *data)
+{
+    FuncGEOS_YY_b_p *func = (FuncGEOS_YY_b_p *)data;
+    void *context_handle = geos_context[0];
+    GEOSGeometry *in1, *in2;
+    GEOSPreparedGeometry *in1_prepared;
+    char ret;
+
+    BINARY_LOOP {
+        /* get the geometries: return on error */
+        if (!get_geom(*(GeometryObject **)ip1, &in1)) { return; }
+        if (!get_geom(*(GeometryObject **)ip2, &in2)) { return; }
+        if (!get_geom_prepared(*(GeometryObject **)ip1, &in1_prepared)) { return; }
+        if ((in1 == NULL) | (in2 == NULL)) {
+            /* in case of a missing value: return 0 (False) */
+            ret = 0;
+        } else {
+            if (in1_prepared == NULL) {
+                ret = GEOSIntersects_r(context_handle, in1, in2);
+            } else {
+                ret = GEOSPreparedIntersects_r(context_handle, in1_prepared, in2);
+            }
+            /* return for illegal values (trust HandleGEOSError for SetErr) */
+            if (ret == 2) { return; }
+        }
+        *(npy_bool *)op1 = ret;
+    }
+}
+static PyUFuncGenericFunction YY_b_p_funcs[1] = {&YY_b_p_func};
+
+
+
 /* Define the geom -> geom functions (Y_Y) */
 static void *envelope_data[1] = {GEOSEnvelope_r};
 static void *convex_hull_data[1] = {GEOSConvexHull_r};
@@ -223,6 +261,40 @@ static void Y_Y_func(char **args, npy_intp *dimensions,
 }
 static PyUFuncGenericFunction Y_Y_funcs[1] = {&Y_Y_func};
 
+
+/* Define the geom -> prepared geom functions (Y_Yp) */
+static void *prepare_data[1] = {GEOSPrepare_r};
+typedef void *FuncGEOS_Y_Yp(void *context, void *a);
+static char Y_Yp_dtypes[2] = {NPY_OBJECT, NPY_OBJECT};
+static void Y_Yp_func(char **args, npy_intp *dimensions,
+                      npy_intp *steps, void *data)
+{
+    FuncGEOS_Y_Yp *func = (FuncGEOS_Y_Yp *)data;
+    void *context_handle = geos_context[0];
+    GEOSGeometry *in1;
+    GEOSPreparedGeometry *ret_ptr;
+
+    UNARY_LOOP {
+        /* get the geometry: return on error */
+        if (!get_geom(*(GeometryObject **)ip1, &in1)) { return; }
+
+        if (in1 == NULL) {
+            ret_ptr = NULL;
+        } else {
+            ret_ptr = func(context_handle, in1);
+            /* trust that GEOS calls HandleGEOSError on error */
+        }
+    GeometryObject *geom_obj = *(GeometryObject **)ip1;
+    geom_obj->ptr_prepared = ret_ptr;
+    /* TODO check if nout can be 0 for ufuncs (because in principle we don't need to return something here) */
+    PyObject **out = (PyObject **)op1;
+    Py_XDECREF(*out);
+    *out = geom_obj;
+    }
+}
+static PyUFuncGenericFunction Y_Yp_funcs[1] = {&Y_Yp_func};
+
+
 /* Define the geom, double -> geom functions (Yd_Y) */
 
 /* GEOS < 3.8 gives segfault for empty linestrings, this is fixed since
@@ -1309,10 +1381,18 @@ TODO relate functions
     ufunc = PyUFunc_FromFuncAndData(YY_b_funcs, NAME ##_data, YY_b_dtypes, 1, 2, 1, PyUFunc_None, # NAME, "", 0);\
     PyDict_SetItemString(d, # NAME, ufunc)
 
+#define DEFINE_YY_b_p(NAME)\
+    ufunc = PyUFunc_FromFuncAndData(YY_b_p_funcs, NAME ##_data, YY_b_p_dtypes, 1, 2, 1, PyUFunc_None, # NAME, "", 0);\
+    PyDict_SetItemString(d, # NAME, ufunc)
+
 #define DEFINE_Y_Y(NAME)\
     ufunc = PyUFunc_FromFuncAndData(Y_Y_funcs, NAME ##_data, Y_Y_dtypes, 1, 1, 1, PyUFunc_None, # NAME, "", 0);\
     PyDict_SetItemString(d, # NAME, ufunc)
 
+#define DEFINE_Y_Yp(NAME)\
+    ufunc = PyUFunc_FromFuncAndData(Y_Yp_funcs, NAME##_data, Y_Yp_dtypes, 1, 1, 1, PyUFunc_None, #NAME, "", 0);\
+    PyDict_SetItemString(d, #NAME, ufunc)
+
 #define DEFINE_Yi_Y(NAME)\
     ufunc = PyUFunc_FromFuncAndData(Yi_Y_funcs, NAME ##_data, Yi_Y_dtypes, 1, 2, 1, PyUFunc_None, # NAME, "", 0);\
     PyDict_SetItemString(d, # NAME, ufunc)
@@ -1376,6 +1456,8 @@ int init_ufuncs(PyObject *m, PyObject *d)
     DEFINE_YY_b (covers);
     DEFINE_YY_b (covered_by);
 
+    DEFINE_YY_b_p (intersects_prepared);
+
     DEFINE_Y_Y (envelope);
     DEFINE_Y_Y (convex_hull);
     DEFINE_Y_Y (boundary);
@@ -1386,6 +1468,8 @@ int init_ufuncs(PyObject *m, PyObject *d)
     DEFINE_Y_Y (extract_unique_points);
     DEFINE_Y_Y (get_exterior_ring);
 
+    DEFINE_Y_Yp(prepare);
+
     DEFINE_Yi_Y (get_point);
     DEFINE_Yi_Y (get_interior_ring);
     DEFINE_Yi_Y (get_geometry);