add sjoin algorithm

geopandas/algos.h

@@ -0,0 +1,23 @@
+#include <geos_c.h>
+#include "kvec.h"
+
+#ifndef GEOPANDAS_ALGOS_C
+#define GEOPANDAS_ALGOS_C
+
+typedef char (*GEOSPredicate)(GEOSContextHandle_t handle, const GEOSGeometry *left, const GEOSGeometry *right);
+typedef char (*GEOSPreparedPredicate)(GEOSContextHandle_t handle, const GEOSPreparedGeometry *left, const GEOSGeometry *right);
+
+typedef struct
+{
+    size_t n, m;
+    size_t *a;
+} size_vector;
+
+void sjoin_callback(void *item, void *vec);
+
+size_vector sjoin(GEOSContextHandle_t handle,
+                  GEOSPreparedPredicate predicate,
+                  GEOSGeometry **left, size_t nleft,
+                  GEOSGeometry **right, size_t nright);
+
+#endif

geopandas/kvec.h

@@ -0,0 +1,90 @@
+/* The MIT License
+
+   Copyright (c) 2008, by Attractive Chaos <attractor@live.co.uk>
+
+   Permission is hereby granted, free of charge, to any person obtaining
+   a copy of this software and associated documentation files (the
+   "Software"), to deal in the Software without restriction, including
+   without limitation the rights to use, copy, modify, merge, publish,
+   distribute, sublicense, and/or sell copies of the Software, and to
+   permit persons to whom the Software is furnished to do so, subject to
+   the following conditions:
+
+   The above copyright notice and this permission notice shall be
+   included in all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+   SOFTWARE.
+*/
+
+/*
+  An example:
+
+#include "kvec.h"
+int main() {
+	kvec_t(int) array;
+	kv_init(array);
+	kv_push(int, array, 10); // append
+	kv_a(int, array, 20) = 5; // dynamic
+	kv_A(array, 20) = 4; // static
+	kv_destroy(array);
+	return 0;
+}
+*/
+
+/*
+  2008-09-22 (0.1.0):
+
+	* The initial version.
+
+*/
+
+#ifndef AC_KVEC_H
+#define AC_KVEC_H
+
+#include <stdlib.h>
+
+#define kv_roundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
+
+#define kvec_t(type) struct { size_t n, m; type *a; }
+#define kv_init(v) ((v).n = (v).m = 0, (v).a = 0)
+#define kv_destroy(v) free((v).a)
+#define kv_A(v, i) ((v).a[(i)])
+#define kv_pop(v) ((v).a[--(v).n])
+#define kv_size(v) ((v).n)
+#define kv_max(v) ((v).m)
+
+#define kv_resize(type, v, s)  ((v).m = (s), (v).a = (type*)realloc((v).a, sizeof(type) * (v).m))
+
+#define kv_copy(type, v1, v0) do {							\
+		if ((v1).m < (v0).n) kv_resize(type, v1, (v0).n);	\
+		(v1).n = (v0).n;									\
+		memcpy((v1).a, (v0).a, sizeof(type) * (v0).n);		\
+	} while (0)												\
+
+#define kv_push(type, v, x) do {									\
+		if ((v).n == (v).m) {										\
+			(v).m = (v).m? (v).m<<1 : 2;							\
+			(v).a = (type*)realloc((v).a, sizeof(type) * (v).m);	\
+		}															\
+		(v).a[(v).n++] = (x);										\
+	} while (0)
+
+#define kv_pushp(type, v) (((v).n == (v).m)?							\
+						   ((v).m = ((v).m? (v).m<<1 : 2),				\
+							(v).a = (type*)realloc((v).a, sizeof(type) * (v).m), 0)	\
+						   : 0), ((v).a + ((v).n++))
+
+#define kv_a(type, v, i) (((v).m <= (size_t)(i)? \
+						  ((v).m = (v).n = (i) + 1, kv_roundup32((v).m), \
+						   (v).a = (type*)realloc((v).a, sizeof(type) * (v).m), 0) \
+						  : (v).n <= (size_t)(i)? (v).n = (i) + 1 \
+						  : 0), (v).a[(i)])
+
+#endif

geopandas/tests/test_vectorized.py

@@ -3,7 +3,7 @@
 import random
 import shapely
 from geopandas.vectorized import (GeometryArray, points_from_xy,
-        from_shapely, serialize, deserialize)
+        from_shapely, serialize, deserialize, cysjoin)
 from shapely.geometry.base import (CAP_STYLE, JOIN_STYLE)
 
 import pytest
@@ -306,6 +306,53 @@ def test_pickle():
     assert T.equals(T2).all()
 
 
+@pytest.mark.parametrize('predicate', [
+    'contains',
+    'covers',
+    'crosses',
+    'disjoint',
+    # 'equals',
+    'intersects',
+    'overlaps',
+    'touches',
+    'within',
+])
+def test_sjoin(predicate):
+    result = cysjoin(T.data, P.data, predicate)
+
+    assert isinstance(result, np.ndarray)
+    assert result.dtype == T.data.dtype
+    n, m = result.shape
+    assert m == 2
+    assert n < (len(T) * len(P))
+
+    for i, j in result:
+        left = triangles[i]
+        right = points[j]
+        assert getattr(left, predicate)(right)
+
+
+@pytest.mark.skip
+def test_bench_sjoin():
+    last = time.time()
+    triangles = [shapely.geometry.Polygon([(random.random(), random.random())
+                                           for i in range(3)])
+                 for _ in range(1000)]
+
+    points = points_from_xy(np.random.random(10000),
+                            np.random.random(10000))
+    print("creation", time.time() - last); last = time.time()
+
+    T = from_shapely(triangles)
+    P = from_shapely(points)
+
+    print("vectorize", time.time() - last); last = time.time()
+
+    result = cysjoin(T.data, P.data, 'intersects')
+
+    print("join", time.time() - last); last = time.time()
+
+
 def test_raise_on_bad_sizes():
     with pytest.raises(ValueError) as info:
         T.contains(P)