In the cython branch in GeoPandas, we had a C implementation of a spatial join algorithm (using the STRTree from GEOS). This was added in geopandas/geopandas#475

Since it's written in C directly against the GEOS C API, it would be nice to include it in this package (primarily for use in geopandas, but it might also be useful in general).

One question is how it would fit in the current API of pygeos.
The version in GeoPandas is currently written as a C function that gets 2 C arrays of geometry objects:

size_vector sjoin(GEOSContextHandle_t handle,
                  GEOSPreparedPredicate predicate,
                  GEOSGeometry **left, size_t nleft,
                  GEOSGeometry **right, size_t nright);

which is then wrapped in cython to provide a python interface with a signature of

def sjoin(left: np.ndarray[1D], right: np.ndarray[1D], predicate_name: str) -> (left_out, right_out)

where left_out and right_out are arrays of indices into left and out, respectively, for all the matches according to the predicate.

I think it is difficult to generalize this to multiple dimensions, as currently all other functions in pygeos do.