1 #define PY_SSIZE_T_CLEAN
2 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
3
4 #include <Python.h>
5 #include <math.h>
6
7 #define NO_IMPORT_ARRAY
8 #define NO_IMPORT_UFUNC
9 #define PY_ARRAY_UNIQUE_SYMBOL pygeos_ARRAY_API
10 #define PY_UFUNC_UNIQUE_SYMBOL pygeos_UFUNC_API
11 #include <numpy/ndarraytypes.h>
12 #include <numpy/npy_3kcompat.h>
13 #include <numpy/ufuncobject.h>
14
15 #include "fast_loop_macros.h"
16 #include "geos.h"
17 #include "pygeom.h"
18
19 #define OUTPUT_Y \
20 PyObject* ret = GeometryObject_FromGEOS(ret_ptr, ctx); \
21 PyObject** out = (PyObject**)op1; \
22 Py_XDECREF(*out); \
23 *out = ret
24
25 #define OUTPUT_Y_I(I, RET_PTR) \
26 PyObject* ret##I = GeometryObject_FromGEOS(RET_PTR, ctx); \
27 PyObject** out##I = (PyObject**)op##I; \
28 Py_XDECREF(*out##I); \
29 *out##I = ret##I
30
31 // Fail if inputs output multiple times on the same place in memory. That would
32 // lead to segfaults as the same GEOSGeometry would be 'owned' by multiple PyObjects.
33 #define CHECK_NO_INPLACE_OUTPUT(N) \
34 if ((steps[N] == 0) && (dimensions[0] > 1)) { \
35 PyErr_Format(PyExc_NotImplementedError, \
36 "Zero-strided output detected. Ufunc mode with args[0]=%p, " \
37 "args[N]=%p, steps[0]=%ld, steps[N]=%ld, dimensions[0]=%ld.", \
38 args[0], args[N], steps[0], steps[N], dimensions[0]); \
39 return; \
40 }
41
42 #define CHECK_ALLOC(ARR) \
43 if (ARR == NULL) { \
44 PyErr_SetString(PyExc_MemoryError, "Could not allocate memory"); \
45 return; \
46 }
47
geom_arr_to_npy(GEOSGeometry ** array,char * ptr,npy_intp stride,npy_intp count)48 static void geom_arr_to_npy(GEOSGeometry** array, char* ptr, npy_intp stride,
49 npy_intp count) {
50 npy_intp i;
51 PyObject* ret;
52 PyObject** out;
53
54 GEOS_INIT;
55
56 for (i = 0; i < count; i++, ptr += stride) {
57 ret = GeometryObject_FromGEOS(array[i], ctx);
58 out = (PyObject**)ptr;
59 Py_XDECREF(*out);
60 *out = ret;
61 }
62
63 GEOS_FINISH;
64 }
65
66 /* Define the geom -> bool functions (Y_b) */
67 static void* is_empty_data[1] = {GEOSisEmpty_r};
68 /* the GEOSisSimple_r function fails on geometrycollections */
GEOSisSimpleAllTypes_r(void * context,void * geom)69 static char GEOSisSimpleAllTypes_r(void* context, void* geom) {
70 int type = GEOSGeomTypeId_r(context, geom);
71 if (type == -1) {
72 return 2; // Predicates use a return value of 2 for errors
73 } else if (type == 7) {
74 return 0;
75 } else {
76 return GEOSisSimple_r(context, geom);
77 }
78 }
79 static void* is_simple_data[1] = {GEOSisSimpleAllTypes_r};
80 static void* is_ring_data[1] = {GEOSisRing_r};
81 static void* has_z_data[1] = {GEOSHasZ_r};
82 /* the GEOSisClosed_r function fails on non-linestrings */
GEOSisClosedAllTypes_r(void * context,void * geom)83 static char GEOSisClosedAllTypes_r(void* context, void* geom) {
84 int type = GEOSGeomTypeId_r(context, geom);
85 if (type == -1) {
86 return 2; // Predicates use a return value of 2 for errors
87 } else if ((type == 1) || (type == 5)) {
88 return GEOSisClosed_r(context, geom);
89 } else {
90 return 0;
91 }
92 }
93 static void* is_closed_data[1] = {GEOSisClosedAllTypes_r};
94 static void* is_valid_data[1] = {GEOSisValid_r};
95
96 #if GEOS_SINCE_3_7_0
GEOSGeom_isCCW_r(void * context,void * geom)97 static char GEOSGeom_isCCW_r(void* context, void* geom) {
98 const GEOSCoordSequence* coord_seq;
99 char is_ccw = 2; // return value of 2 means GEOSException
100 int i;
101
102 // Return False for non-linear geometries
103 i = GEOSGeomTypeId_r(context, geom);
104 if (i == -1) {
105 return 2;
106 }
107 if ((i != GEOS_LINEARRING) && (i != GEOS_LINESTRING)) {
108 return 0;
109 }
110
111 // Return False for lines with fewer than 4 points
112 i = GEOSGeomGetNumPoints_r(context, geom);
113 if (i == -1) {
114 return 2;
115 }
116 if (i < 4) {
117 return 0;
118 }
119
120 // Get the coordinatesequence and call isCCW()
121 coord_seq = GEOSGeom_getCoordSeq_r(context, geom);
122 if (coord_seq == NULL) {
123 return 2;
124 }
125 if (!GEOSCoordSeq_isCCW_r(context, coord_seq, &is_ccw)) {
126 return 2;
127 }
128 return is_ccw;
129 }
130 static void* is_ccw_data[1] = {GEOSGeom_isCCW_r};
131 #endif
132
133 typedef char FuncGEOS_Y_b(void* context, void* a);
134 static char Y_b_dtypes[2] = {NPY_OBJECT, NPY_BOOL};
Y_b_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)135 static void Y_b_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
136 FuncGEOS_Y_b* func = (FuncGEOS_Y_b*)data;
137 GEOSGeometry* in1 = NULL;
138 char ret;
139
140 GEOS_INIT_THREADS;
141
142 UNARY_LOOP {
143 /* get the geometry; return on error */
144 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
145 errstate = PGERR_NOT_A_GEOMETRY;
146 goto finish;
147 }
148 if (in1 == NULL) {
149 /* in case of a missing value: return 0 (False) */
150 ret = 0;
151 } else {
152 /* call the GEOS function */
153 ret = func(ctx, in1);
154 /* finish for illegal values */
155 if (ret == 2) {
156 errstate = PGERR_GEOS_EXCEPTION;
157 goto finish;
158 }
159 }
160 *(npy_bool*)op1 = ret;
161 }
162
163 finish:
164 GEOS_FINISH_THREADS;
165 }
166 static PyUFuncGenericFunction Y_b_funcs[1] = {&Y_b_func};
167
168 /* Define the object -> bool functions (O_b) which do not raise on non-geom objects*/
IsMissing(void * context,PyObject * obj)169 static char IsMissing(void* context, PyObject* obj) {
170 GEOSGeometry* g = NULL;
171 if (!get_geom((GeometryObject*)obj, &g)) {
172 return 0;
173 };
174 return g == NULL; // get_geom sets g to NULL for None input
175 }
176 static void* is_missing_data[1] = {IsMissing};
IsGeometry(void * context,PyObject * obj)177 static char IsGeometry(void* context, PyObject* obj) {
178 GEOSGeometry* g = NULL;
179 if (!get_geom((GeometryObject*)obj, &g)) {
180 return 0;
181 }
182 return g != NULL;
183 }
184 static void* is_geometry_data[1] = {IsGeometry};
IsValidInput(void * context,PyObject * obj)185 static char IsValidInput(void* context, PyObject* obj) {
186 GEOSGeometry* g = NULL;
187 return get_geom((GeometryObject*)obj, &g);
188 }
189 static void* is_valid_input_data[1] = {IsValidInput};
190 typedef char FuncGEOS_O_b(void* context, PyObject* obj);
191 static char O_b_dtypes[2] = {NPY_OBJECT, NPY_BOOL};
O_b_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)192 static void O_b_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
193 FuncGEOS_O_b* func = (FuncGEOS_O_b*)data;
194 GEOS_INIT_THREADS;
195 UNARY_LOOP { *(npy_bool*)op1 = func(ctx, *(PyObject**)ip1); }
196 GEOS_FINISH_THREADS;
197 }
198 static PyUFuncGenericFunction O_b_funcs[1] = {&O_b_func};
199
200 /* Define the geom, geom -> bool functions (YY_b) */
201 static void* equals_data[1] = {GEOSEquals_r};
202 typedef char FuncGEOS_YY_b(void* context, void* a, void* b);
203 static char YY_b_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_BOOL};
YY_b_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)204 static void YY_b_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
205 FuncGEOS_YY_b* func = (FuncGEOS_YY_b*)data;
206 GEOSGeometry *in1 = NULL, *in2 = NULL;
207 char ret;
208
209 GEOS_INIT_THREADS;
210
211 BINARY_LOOP {
212 /* get the geometries: return on error */
213 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
214 errstate = PGERR_NOT_A_GEOMETRY;
215 goto finish;
216 }
217 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
218 errstate = PGERR_NOT_A_GEOMETRY;
219 goto finish;
220 }
221 if ((in1 == NULL) || (in2 == NULL)) {
222 /* in case of a missing value: return 0 (False) */
223 ret = 0;
224 } else {
225 /* call the GEOS function */
226 ret = func(ctx, in1, in2);
227 /* return for illegal values */
228 if (ret == 2) {
229 errstate = PGERR_GEOS_EXCEPTION;
230 goto finish;
231 }
232 }
233 *(npy_bool*)op1 = ret;
234 }
235
236 finish:
237 GEOS_FINISH_THREADS;
238 }
239 static PyUFuncGenericFunction YY_b_funcs[1] = {&YY_b_func};
240
241 /* Define the geom, geom -> bool functions (YY_b) prepared */
242 static void* contains_func_tuple[2] = {GEOSContains_r, GEOSPreparedContains_r};
243 static void* contains_data[1] = {contains_func_tuple};
GEOSContainsProperly(void * context,void * g1,void * g2)244 static char GEOSContainsProperly(void* context, void* g1, void* g2) {
245 const GEOSPreparedGeometry* prepared_geom_tmp = NULL;
246 char ret;
247
248 prepared_geom_tmp = GEOSPrepare_r(context, g1);
249 if (prepared_geom_tmp == NULL) {
250 return 2;
251 }
252 ret = GEOSPreparedContainsProperly_r(context, prepared_geom_tmp, g2);
253 GEOSPreparedGeom_destroy_r(context, prepared_geom_tmp);
254 return ret;
255 }
256 static void* contains_properly_func_tuple[2] = {GEOSContainsProperly,
257 GEOSPreparedContainsProperly_r};
258 static void* contains_properly_data[1] = {contains_properly_func_tuple};
259 static void* covered_by_func_tuple[2] = {GEOSCoveredBy_r, GEOSPreparedCoveredBy_r};
260 static void* covered_by_data[1] = {covered_by_func_tuple};
261 static void* covers_func_tuple[2] = {GEOSCovers_r, GEOSPreparedCovers_r};
262 static void* covers_data[1] = {covers_func_tuple};
263 static void* crosses_func_tuple[2] = {GEOSCrosses_r, GEOSPreparedCrosses_r};
264 static void* crosses_data[1] = {crosses_func_tuple};
265 static void* disjoint_func_tuple[2] = {GEOSDisjoint_r, GEOSPreparedDisjoint_r};
266 static void* disjoint_data[1] = {disjoint_func_tuple};
267 static void* intersects_func_tuple[2] = {GEOSIntersects_r, GEOSPreparedIntersects_r};
268 static void* intersects_data[1] = {intersects_func_tuple};
269 static void* overlaps_func_tuple[2] = {GEOSOverlaps_r, GEOSPreparedOverlaps_r};
270 static void* overlaps_data[1] = {overlaps_func_tuple};
271 static void* touches_func_tuple[2] = {GEOSTouches_r, GEOSPreparedTouches_r};
272 static void* touches_data[1] = {touches_func_tuple};
273 static void* within_func_tuple[2] = {GEOSWithin_r, GEOSPreparedWithin_r};
274 static void* within_data[1] = {within_func_tuple};
275 static char YY_b_p_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_BOOL};
YY_b_p_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)276 static void YY_b_p_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
277 FuncGEOS_YY_b* func = ((FuncGEOS_YY_b**)data)[0];
278 FuncGEOS_YY_b* func_prepared = ((FuncGEOS_YY_b**)data)[1];
279
280 GEOSGeometry *in1 = NULL, *in2 = NULL;
281 GEOSPreparedGeometry* in1_prepared = NULL;
282 char ret;
283
284 GEOS_INIT_THREADS;
285
286 BINARY_LOOP {
287 /* get the geometries: return on error */
288 if (!get_geom_with_prepared(*(GeometryObject**)ip1, &in1, &in1_prepared)) {
289 errstate = PGERR_NOT_A_GEOMETRY;
290 goto finish;
291 }
292 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
293 errstate = PGERR_NOT_A_GEOMETRY;
294 goto finish;
295 }
296 if ((in1 == NULL) || (in2 == NULL)) {
297 /* in case of a missing value: return 0 (False) */
298 ret = 0;
299 } else {
300 if (in1_prepared == NULL) {
301 /* call the GEOS function */
302 ret = func(ctx, in1, in2);
303 } else {
304 /* call the prepared GEOS function */
305 ret = func_prepared(ctx, in1_prepared, in2);
306 }
307 /* return for illegal values */
308 if (ret == 2) {
309 errstate = PGERR_GEOS_EXCEPTION;
310 goto finish;
311 }
312 }
313 *(npy_bool*)op1 = ret;
314 }
315
316 finish:
317
318 GEOS_FINISH_THREADS;
319 }
320 static PyUFuncGenericFunction YY_b_p_funcs[1] = {&YY_b_p_func};
321
322 static char is_prepared_dtypes[2] = {NPY_OBJECT, NPY_BOOL};
is_prepared_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)323 static void is_prepared_func(char** args, npy_intp* dimensions, npy_intp* steps,
324 void* data) {
325 GEOSGeometry* in1 = NULL;
326 GEOSPreparedGeometry* in1_prepared = NULL;
327
328 GEOS_INIT_THREADS;
329
330 UNARY_LOOP {
331 /* get the geometry: return on error */
332 if (!get_geom_with_prepared(*(GeometryObject**)ip1, &in1, &in1_prepared)) {
333 errstate = PGERR_NOT_A_GEOMETRY;
334 break;
335 }
336 *(npy_bool*)op1 = (in1_prepared != NULL);
337 }
338
339 GEOS_FINISH_THREADS;
340 }
341 static PyUFuncGenericFunction is_prepared_funcs[1] = {&is_prepared_func};
342
343 /* Define the geom -> geom functions (Y_Y) */
344 static void* envelope_data[1] = {GEOSEnvelope_r};
345 static void* convex_hull_data[1] = {GEOSConvexHull_r};
GEOSBoundaryAllTypes_r(void * context,void * geom)346 static void* GEOSBoundaryAllTypes_r(void* context, void* geom) {
347 char typ = GEOSGeomTypeId_r(context, geom);
348 if (typ == 7) {
349 /* return None for geometrycollections */
350 return NULL;
351 } else {
352 return GEOSBoundary_r(context, geom);
353 }
354 }
355 static void* boundary_data[1] = {GEOSBoundaryAllTypes_r};
356 static void* unary_union_data[1] = {GEOSUnaryUnion_r};
357 static void* point_on_surface_data[1] = {GEOSPointOnSurface_r};
358 static void* centroid_data[1] = {GEOSGetCentroid_r};
359 static void* line_merge_data[1] = {GEOSLineMerge_r};
360 static void* extract_unique_points_data[1] = {GEOSGeom_extractUniquePoints_r};
GetExteriorRing(void * context,void * geom)361 static void* GetExteriorRing(void* context, void* geom) {
362 char typ = GEOSGeomTypeId_r(context, geom);
363 if (typ != 3) {
364 return NULL;
365 }
366 void* ret = (void*)GEOSGetExteriorRing_r(context, geom);
367 /* Create a copy of the obtained geometry */
368 if (ret != NULL) {
369 ret = GEOSGeom_clone_r(context, ret);
370 }
371 return ret;
372 }
373 static void* get_exterior_ring_data[1] = {GetExteriorRing};
374 /* the normalize funcion acts inplace */
GEOSNormalize_r_with_clone(void * context,void * geom)375 static void* GEOSNormalize_r_with_clone(void* context, void* geom) {
376 int ret;
377 void* new_geom = GEOSGeom_clone_r(context, geom);
378 if (new_geom == NULL) {
379 return NULL;
380 }
381 ret = GEOSNormalize_r(context, new_geom);
382 if (ret == -1) {
383 GEOSGeom_destroy_r(context, new_geom);
384 return NULL;
385 }
386 return new_geom;
387 }
388 static void* normalize_data[1] = {GEOSNormalize_r_with_clone};
389 #if GEOS_SINCE_3_8_0
390 static void* build_area_data[1] = {GEOSBuildArea_r};
391 static void* make_valid_data[1] = {GEOSMakeValid_r};
392 static void* coverage_union_data[1] = {GEOSCoverageUnion_r};
GEOSMinimumBoundingCircleWithReturn(void * context,void * geom)393 static void* GEOSMinimumBoundingCircleWithReturn(void* context, void* geom) {
394 GEOSGeometry* center = NULL;
395 double radius;
396 GEOSGeometry* ret = GEOSMinimumBoundingCircle_r(context, geom, &radius, ¢er);
397 if (ret == NULL) {
398 return NULL;
399 }
400 GEOSGeom_destroy_r(context, center);
401 return ret;
402 }
403 static void* minimum_bounding_circle_data[1] = {GEOSMinimumBoundingCircleWithReturn};
404 #endif
405 #if GEOS_SINCE_3_7_0
406 static void* reverse_data[1] = {GEOSReverse_r};
407 #endif
408 #if GEOS_SINCE_3_6_0
409 static void* oriented_envelope_data[1] = {GEOSMinimumRotatedRectangle_r};
410 #endif
411 typedef void* FuncGEOS_Y_Y(void* context, void* a);
412 static char Y_Y_dtypes[2] = {NPY_OBJECT, NPY_OBJECT};
Y_Y_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)413 static void Y_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
414 FuncGEOS_Y_Y* func = (FuncGEOS_Y_Y*)data;
415 GEOSGeometry* in1 = NULL;
416 GEOSGeometry** geom_arr;
417
418 CHECK_NO_INPLACE_OUTPUT(1);
419
420 // allocate a temporary array to store output GEOSGeometry objects
421 geom_arr = malloc(sizeof(void*) * dimensions[0]);
422 CHECK_ALLOC(geom_arr);
423
424 GEOS_INIT_THREADS;
425
426 UNARY_LOOP {
427 // get the geometry: return on error
428 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
429 errstate = PGERR_NOT_A_GEOMETRY;
430 destroy_geom_arr(ctx, geom_arr, i - 1);
431 break;
432 }
433 if (in1 == NULL) {
434 // in case of a missing value: return NULL (None)
435 geom_arr[i] = NULL;
436 } else {
437 geom_arr[i] = func(ctx, in1);
438 // NULL means: exception, but for some functions it may also indicate a
439 // "missing value" (None) (GetExteriorRing, GEOSBoundaryAllTypes_r)
440 // So: check the last_error before setting error state
441 if ((geom_arr[i] == NULL) && (last_error[0] != 0)) {
442 errstate = PGERR_GEOS_EXCEPTION;
443 destroy_geom_arr(ctx, geom_arr, i - 1);
444 break;
445 }
446 }
447 }
448
449 GEOS_FINISH_THREADS;
450
451 // fill the numpy array with PyObjects while holding the GIL
452 if (errstate == PGERR_SUCCESS) {
453 geom_arr_to_npy(geom_arr, args[1], steps[1], dimensions[0]);
454 }
455 free(geom_arr);
456 }
457 static PyUFuncGenericFunction Y_Y_funcs[1] = {&Y_Y_func};
458
459 /* Define the geom -> no return value functions (Y) */
PrepareGeometryObject(void * ctx,GeometryObject * geom)460 static char PrepareGeometryObject(void* ctx, GeometryObject* geom) {
461 if (geom->ptr_prepared == NULL) {
462 geom->ptr_prepared = (GEOSPreparedGeometry*)GEOSPrepare_r(ctx, geom->ptr);
463 if (geom->ptr_prepared == NULL) {
464 return PGERR_GEOS_EXCEPTION;
465 }
466 }
467 return PGERR_SUCCESS;
468 }
DestroyPreparedGeometryObject(void * ctx,GeometryObject * geom)469 static char DestroyPreparedGeometryObject(void* ctx, GeometryObject* geom) {
470 if (geom->ptr_prepared != NULL) {
471 GEOSPreparedGeom_destroy_r(ctx, geom->ptr_prepared);
472 geom->ptr_prepared = NULL;
473 }
474 return PGERR_SUCCESS;
475 }
476
477 static void* prepare_data[1] = {PrepareGeometryObject};
478 static void* destroy_prepared_data[1] = {DestroyPreparedGeometryObject};
479 typedef char FuncPyGEOS_Y(void* ctx, GeometryObject* geom);
480 static char Y_dtypes[1] = {NPY_OBJECT};
Y_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)481 static void Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
482 FuncPyGEOS_Y* func = (FuncPyGEOS_Y*)data;
483 GEOSGeometry* in1 = NULL;
484 GeometryObject* geom_obj = NULL;
485
486 GEOS_INIT;
487
488 NO_OUTPUT_LOOP {
489 geom_obj = *(GeometryObject**)ip1;
490 if (!get_geom(geom_obj, &in1)) {
491 errstate = PGERR_GEOS_EXCEPTION;
492 goto finish;
493 }
494 if (in1 != NULL) {
495 errstate = func(ctx, geom_obj);
496 if (errstate != PGERR_SUCCESS) {
497 goto finish;
498 }
499 }
500 }
501
502 finish:
503
504 GEOS_FINISH;
505 }
506 static PyUFuncGenericFunction Y_funcs[1] = {&Y_func};
507
508 /* Define the geom, double -> geom functions (Yd_Y) */
GEOSInterpolateProtectEmpty_r(void * context,void * geom,double d)509 static void* GEOSInterpolateProtectEmpty_r(void* context, void* geom, double d) {
510 char errstate = geos_interpolate_checker(context, geom);
511 if (errstate == PGERR_SUCCESS) {
512 return GEOSInterpolate_r(context, geom, d);
513 } else if (errstate == PGERR_EMPTY_GEOMETRY) {
514 return GEOSGeom_createEmptyPoint_r(context);
515 } else {
516 return NULL;
517 }
518 }
519 static void* line_interpolate_point_data[1] = {GEOSInterpolateProtectEmpty_r};
GEOSInterpolateNormalizedProtectEmpty_r(void * context,void * geom,double d)520 static void* GEOSInterpolateNormalizedProtectEmpty_r(void* context, void* geom,
521 double d) {
522 char errstate = geos_interpolate_checker(context, geom);
523 if (errstate == PGERR_SUCCESS) {
524 return GEOSInterpolateNormalized_r(context, geom, d);
525 } else if (errstate == PGERR_EMPTY_GEOMETRY) {
526 return GEOSGeom_createEmptyPoint_r(context);
527 } else {
528 return NULL;
529 }
530 }
531 static void* line_interpolate_point_normalized_data[1] = {
532 GEOSInterpolateNormalizedProtectEmpty_r};
533
534 static void* simplify_data[1] = {GEOSSimplify_r};
535 static void* simplify_preserve_topology_data[1] = {GEOSTopologyPreserveSimplify_r};
536
537 #if GEOS_SINCE_3_9_0
538 static void* unary_union_prec_data[1] = {GEOSUnaryUnionPrec_r};
539 #endif
540
541 #if GEOS_SINCE_3_10_0
542 static void* segmentize_data[1] = {GEOSDensify_r};
543 #endif
544
545 typedef void* FuncGEOS_Yd_Y(void* context, void* a, double b);
546 static char Yd_Y_dtypes[3] = {NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT};
Yd_Y_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)547 static void Yd_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
548 FuncGEOS_Yd_Y* func = (FuncGEOS_Yd_Y*)data;
549 GEOSGeometry* in1 = NULL;
550 GEOSGeometry** geom_arr;
551
552 CHECK_NO_INPLACE_OUTPUT(2);
553
554 // allocate a temporary array to store output GEOSGeometry objects
555 geom_arr = malloc(sizeof(void*) * dimensions[0]);
556 CHECK_ALLOC(geom_arr);
557
558 GEOS_INIT_THREADS;
559
560 BINARY_LOOP {
561 // get the geometry: return on error
562 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
563 errstate = PGERR_NOT_A_GEOMETRY;
564 destroy_geom_arr(ctx, geom_arr, i - 1);
565 break;
566 }
567 double in2 = *(double*)ip2;
568 if ((in1 == NULL) || (npy_isnan(in2))) {
569 // in case of a missing value: return NULL (None)
570 geom_arr[i] = NULL;
571 } else {
572 geom_arr[i] = func(ctx, in1, in2);
573 if (geom_arr[i] == NULL) {
574 // Interpolate functions return NULL on PGERR_GEOMETRY_TYPE and on
575 // PGERR_GEOS_EXCEPTION. Distinguish these by the state of last_error.
576 errstate = last_error[0] == 0 ? PGERR_GEOMETRY_TYPE : PGERR_GEOS_EXCEPTION;
577 destroy_geom_arr(ctx, geom_arr, i - 1);
578 break;
579 }
580 }
581 }
582
583 GEOS_FINISH_THREADS;
584
585 // fill the numpy array with PyObjects while holding the GIL
586 if (errstate == PGERR_SUCCESS) {
587 geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
588 }
589 free(geom_arr);
590 }
591 static PyUFuncGenericFunction Yd_Y_funcs[1] = {&Yd_Y_func};
592
593 /* Define the geom, int -> geom functions (Yi_Y) */
594 /* We add bound and type checking to the various indexing functions */
GetPointN(void * context,void * geom,int n)595 static void* GetPointN(void* context, void* geom, int n) {
596 char typ = GEOSGeomTypeId_r(context, geom);
597 int size, i;
598 if ((typ != 1) && (typ != 2)) {
599 return NULL;
600 }
601 size = GEOSGeomGetNumPoints_r(context, geom);
602 if (size == -1) {
603 return NULL;
604 }
605 if (n < 0) {
606 /* Negative indexing: we get it for free */
607 i = size + n;
608 } else {
609 i = n;
610 }
611 if ((i < 0) || (i >= size)) {
612 /* Important, could give segfaults else */
613 return NULL;
614 }
615 return GEOSGeomGetPointN_r(context, geom, i);
616 }
617 static void* get_point_data[1] = {GetPointN};
GetInteriorRingN(void * context,void * geom,int n)618 static void* GetInteriorRingN(void* context, void* geom, int n) {
619 char typ = GEOSGeomTypeId_r(context, geom);
620 int size, i;
621 if (typ != 3) {
622 return NULL;
623 }
624 size = GEOSGetNumInteriorRings_r(context, geom);
625 if (size == -1) {
626 return NULL;
627 }
628 if (n < 0) {
629 /* Negative indexing: we get it for free */
630 i = size + n;
631 } else {
632 i = n;
633 }
634 if ((i < 0) || (i >= size)) {
635 /* Important, could give segfaults else */
636 return NULL;
637 }
638 void* ret = (void*)GEOSGetInteriorRingN_r(context, geom, i);
639 /* Create a copy of the obtained geometry */
640 if (ret != NULL) {
641 ret = GEOSGeom_clone_r(context, ret);
642 }
643 return ret;
644 }
645 static void* get_interior_ring_data[1] = {GetInteriorRingN};
GetGeometryN(void * context,void * geom,int n)646 static void* GetGeometryN(void* context, void* geom, int n) {
647 int size, i;
648 size = GEOSGetNumGeometries_r(context, geom);
649 if (size == -1) {
650 return NULL;
651 }
652 if (n < 0) {
653 /* Negative indexing: we get it for free */
654 i = size + n;
655 } else {
656 i = n;
657 }
658 if ((i < 0) || (i >= size)) {
659 /* Important, could give segfaults else */
660 return NULL;
661 }
662 void* ret = (void*)GEOSGetGeometryN_r(context, geom, i);
663 /* Create a copy of the obtained geometry */
664 if (ret != NULL) {
665 ret = GEOSGeom_clone_r(context, ret);
666 }
667 return ret;
668 }
669 static void* get_geometry_data[1] = {GetGeometryN};
670 /* the set srid funcion acts inplace */
GEOSSetSRID_r_with_clone(void * context,void * geom,int srid)671 static void* GEOSSetSRID_r_with_clone(void* context, void* geom, int srid) {
672 void* ret = GEOSGeom_clone_r(context, geom);
673 if (ret == NULL) {
674 return NULL;
675 }
676 GEOSSetSRID_r(context, ret, srid);
677 return ret;
678 }
679 static void* set_srid_data[1] = {GEOSSetSRID_r_with_clone};
680 typedef void* FuncGEOS_Yi_Y(void* context, void* a, int b);
681 static char Yi_Y_dtypes[3] = {NPY_OBJECT, NPY_INT, NPY_OBJECT};
Yi_Y_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)682 static void Yi_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
683 FuncGEOS_Yi_Y* func = (FuncGEOS_Yi_Y*)data;
684 GEOSGeometry* in1 = NULL;
685 GEOSGeometry** geom_arr;
686
687 CHECK_NO_INPLACE_OUTPUT(2);
688
689 // allocate a temporary array to store output GEOSGeometry objects
690 geom_arr = malloc(sizeof(void*) * dimensions[0]);
691 CHECK_ALLOC(geom_arr);
692
693 GEOS_INIT_THREADS;
694
695 BINARY_LOOP {
696 // get the geometry: return on error
697 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
698 errstate = PGERR_NOT_A_GEOMETRY;
699 destroy_geom_arr(ctx, geom_arr, i - 1);
700 break;
701 }
702 int in2 = *(int*)ip2;
703 if (in1 == NULL) {
704 // in case of a missing value: return NULL (None)
705 geom_arr[i] = NULL;
706 } else {
707 geom_arr[i] = func(ctx, in1, in2);
708 // NULL means: exception, but for some functions it may also indicate a
709 // "missing value" (None) (GetPointN, GetInteriorRingN, GetGeometryN)
710 // So: check the last_error before setting error state
711 if ((geom_arr[i] == NULL) && (last_error[0] != 0)) {
712 errstate = PGERR_GEOS_EXCEPTION;
713 destroy_geom_arr(ctx, geom_arr, i - 1);
714 break;
715 }
716 }
717 }
718
719 GEOS_FINISH_THREADS;
720
721 // fill the numpy array with PyObjects while holding the GIL
722 if (errstate == PGERR_SUCCESS) {
723 geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
724 }
725 free(geom_arr);
726 }
727 static PyUFuncGenericFunction Yi_Y_funcs[1] = {&Yi_Y_func};
728
729 /* Define the geom, geom -> geom functions (YY_Y) */
730 static void* intersection_data[1] = {GEOSIntersection_r};
731 static void* difference_data[1] = {GEOSDifference_r};
732 static void* symmetric_difference_data[1] = {GEOSSymDifference_r};
733 static void* union_data[1] = {GEOSUnion_r};
734 static void* shared_paths_data[1] = {GEOSSharedPaths_r};
735 typedef void* FuncGEOS_YY_Y(void* context, void* a, void* b);
736 static char YY_Y_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
737
738 /* There are two inner loop functions for the YY_Y. This is because this
739 * pattern allows for ufunc.reduce.
740 * A reduce operation requires special attention, because we output the
741 * result in a temporary array. NumPy expects that the output array is
742 * filled during the inner loop, so that it can take the first input of
743 * the reduction operation to be the output array. Easiest to show by
744 * example:
745
746 * function called: out = intersection.reduce(in)
747 * initialization by numpy: out[0] = in[0]; in1 = out; in2[:] = in[:]
748 * first loop: out[0] = func(in1[0], in2[1]) [ = func(out[0], in2[1]) ]
749 * second loop: out[0] = func(in1[0], in2[2]) [ = func(out[0], in2[2]) ]
750 */
YY_Y_func_reduce(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)751 static void YY_Y_func_reduce(char** args, npy_intp* dimensions, npy_intp* steps,
752 void* data) {
753 FuncGEOS_YY_Y* func = (FuncGEOS_YY_Y*)data;
754 GEOSGeometry *in1 = NULL, *in2 = NULL, *out = NULL;
755
756 // Whether to destroy a temporary intermediate value of `out`:
757 char do_destroy = 0;
758
759 GEOS_INIT_THREADS;
760
761 if (!get_geom(*(GeometryObject**)args[0], &out)) {
762 errstate = PGERR_NOT_A_GEOMETRY;
763 } else {
764 BINARY_LOOP {
765 // Get the geometry inputs; in1 from previous iteration, in2 from array
766 in1 = out;
767 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
768 errstate = PGERR_NOT_A_GEOMETRY;
769 break;
770 }
771
772 /* Either (or both) in1 and in2 could be NULL (Python: None).
773 * Reduction operations should skip None values. We have 4 possible combinations:
774 */
775
776 // 1. (not NULL, not NULL); run the GEOS function
777 if ((in1 != NULL) && (in2 != NULL)) {
778 out = func(ctx, in1, in2);
779
780 // Discard in1 if it was a temporary intermediate
781 if (do_destroy) {
782 GEOSGeom_destroy_r(ctx, in1);
783 }
784
785 // Mark the newly generated geometry as intermediate. Note: out will become in1.
786 do_destroy = 1;
787
788 // Break on error (we do this after discarding in1 to avoid memleaks)
789 if (out == NULL) {
790 errstate = PGERR_GEOS_EXCEPTION;
791 break;
792 }
793 }
794
795 // 2. (NULL, not NULL); When the first element of the reduction axis is None
796 else if ((in1 == NULL) && (in2 != NULL)) {
797 // Keep in2 as 'outcome' of the operation.
798 out = in2;
799 // Ensure that it will not be destroyed (it is owned by python)
800 do_destroy = 0;
801 }
802
803 // 3. (not NULL, NULL); When a None value is encountered after a not-None
804 // Don't do `out = in1`, as that is already the case.
805 // 4. (NULL, NULL); When we have not yet encountered any not-None
806 // Do nothing; out will remain NULL
807 }
808 }
809
810 GEOS_FINISH_THREADS;
811
812 // fill the numpy array with a single PyObject while holding the GIL
813 if (errstate == PGERR_SUCCESS) {
814 geom_arr_to_npy(&out, args[2], steps[2], 1);
815 }
816 }
817
YY_Y_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)818 static void YY_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
819 // A reduce is characterized by multiple iterations (dimension[0] > 1) that
820 // are output on the same place in memory (steps[2] == 0).
821 if ((steps[2] == 0) && (dimensions[0] > 1)) {
822 if (args[0] == args[2]) {
823 YY_Y_func_reduce(args, dimensions, steps, data);
824 return;
825 } else {
826 // Fail if inputs do not have the expected structure.
827 PyErr_Format(PyExc_NotImplementedError,
828 "Unknown ufunc mode with args=[%p, %p, %p], steps=[%ld, %ld, %ld], "
829 "dimensions=[%ld].",
830 args[0], args[1], args[2], steps[0], steps[1], steps[2],
831 dimensions[0]);
832 return;
833 }
834 }
835
836 FuncGEOS_YY_Y* func = (FuncGEOS_YY_Y*)data;
837 GEOSGeometry *in1 = NULL, *in2 = NULL;
838 GEOSGeometry** geom_arr;
839
840 // allocate a temporary array to store output GEOSGeometry objects
841 geom_arr = malloc(sizeof(void*) * dimensions[0]);
842 CHECK_ALLOC(geom_arr);
843
844 GEOS_INIT_THREADS;
845
846 BINARY_LOOP {
847 // get the geometries: return on error
848 if (!get_geom(*(GeometryObject**)ip1, &in1) ||
849 !get_geom(*(GeometryObject**)ip2, &in2)) {
850 errstate = PGERR_NOT_A_GEOMETRY;
851 destroy_geom_arr(ctx, geom_arr, i - 1);
852 break;
853 }
854 if ((in1 == NULL) || (in2 == NULL)) {
855 // in case of a missing value: return NULL (None)
856 geom_arr[i] = NULL;
857 } else {
858 geom_arr[i] = func(ctx, in1, in2);
859 if (geom_arr[i] == NULL) {
860 errstate = PGERR_GEOS_EXCEPTION;
861 destroy_geom_arr(ctx, geom_arr, i - 1);
862 break;
863 }
864 }
865 }
866
867 GEOS_FINISH_THREADS;
868
869 // fill the numpy array with PyObjects while holding the GIL
870 if (errstate == PGERR_SUCCESS) {
871 geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
872 }
873 free(geom_arr);
874 }
875 static PyUFuncGenericFunction YY_Y_funcs[1] = {&YY_Y_func};
876
877 /* Define the geom -> double functions (Y_d) */
GetX(void * context,void * a,double * b)878 static int GetX(void* context, void* a, double* b) {
879 char typ = GEOSGeomTypeId_r(context, a);
880 if (typ != 0) {
881 *(double*)b = NPY_NAN;
882 return 1;
883 } else {
884 return GEOSGeomGetX_r(context, a, b);
885 }
886 }
887 static void* get_x_data[1] = {GetX};
GetY(void * context,void * a,double * b)888 static int GetY(void* context, void* a, double* b) {
889 char typ = GEOSGeomTypeId_r(context, a);
890 if (typ != 0) {
891 *(double*)b = NPY_NAN;
892 return 1;
893 } else {
894 return GEOSGeomGetY_r(context, a, b);
895 }
896 }
897 static void* get_y_data[1] = {GetY};
898 #if GEOS_SINCE_3_7_0
GetZ(void * context,void * a,double * b)899 static int GetZ(void* context, void* a, double* b) {
900 char typ = GEOSGeomTypeId_r(context, a);
901 if (typ != 0) {
902 *(double*)b = NPY_NAN;
903 return 1;
904 } else {
905 return GEOSGeomGetZ_r(context, a, b);
906 }
907 }
908 static void* get_z_data[1] = {GetZ};
909 #endif
910 static void* area_data[1] = {GEOSArea_r};
911 static void* length_data[1] = {GEOSLength_r};
912
913 #if GEOS_SINCE_3_6_0
GetPrecision(void * context,void * a,double * b)914 static int GetPrecision(void* context, void* a, double* b) {
915 // GEOS returns -1 on error; 0 indicates double precision; > 0 indicates a precision
916 // grid size was set for this geometry.
917 double out = GEOSGeom_getPrecision_r(context, a);
918 if (out == -1) {
919 return 0;
920 }
921 *(double*)b = out;
922 return 1;
923 }
924 static void* get_precision_data[1] = {GetPrecision};
MinimumClearance(void * context,void * a,double * b)925 static int MinimumClearance(void* context, void* a, double* b) {
926 // GEOSMinimumClearance deviates from the pattern of returning 0 on exception and 1 on
927 // success for functions that return an int (it follows pattern for boolean functions
928 // returning char 0/1 and 2 on exception)
929 int retcode = GEOSMinimumClearance_r(context, a, b);
930 if (retcode == 2) {
931 return 0;
932 } else {
933 return 1;
934 }
935 }
936 static void* minimum_clearance_data[1] = {MinimumClearance};
937 #endif
938 #if GEOS_SINCE_3_8_0
GEOSMinimumBoundingRadius(void * context,GEOSGeometry * geom,double * radius)939 static int GEOSMinimumBoundingRadius(void* context, GEOSGeometry* geom, double* radius) {
940 GEOSGeometry* center = NULL;
941 GEOSGeometry* ret = GEOSMinimumBoundingCircle_r(context, geom, radius, ¢er);
942 if (ret == NULL) {
943 return 0; // exception code
944 }
945 GEOSGeom_destroy_r(context, center);
946 GEOSGeom_destroy_r(context, ret);
947 return 1; // success code
948 }
949 static void* minimum_bounding_radius_data[1] = {GEOSMinimumBoundingRadius};
950 #endif
951 typedef int FuncGEOS_Y_d(void* context, void* a, double* b);
952 static char Y_d_dtypes[2] = {NPY_OBJECT, NPY_DOUBLE};
Y_d_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)953 static void Y_d_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
954 FuncGEOS_Y_d* func = (FuncGEOS_Y_d*)data;
955 GEOSGeometry* in1 = NULL;
956
957 GEOS_INIT_THREADS;
958
959 UNARY_LOOP {
960 /* get the geometry: return on error */
961 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
962 errstate = PGERR_NOT_A_GEOMETRY;
963 goto finish;
964 }
965 if (in1 == NULL) {
966 *(double*)op1 = NPY_NAN;
967 } else {
968 /* let the GEOS function set op1; return on error */
969 if (func(ctx, in1, (npy_double*)op1) == 0) {
970 errstate = PGERR_GEOS_EXCEPTION;
971 goto finish;
972 }
973 }
974 }
975
976 finish:
977 GEOS_FINISH_THREADS;
978 }
979 static PyUFuncGenericFunction Y_d_funcs[1] = {&Y_d_func};
980
981 /* Define the geom -> int functions (Y_i) */
982 /* data values are GEOS func, GEOS error code, return value when input is None */
983 static void* get_type_id_func_tuple[3] = {GEOSGeomTypeId_r, (void*)-1, (void*)-1};
984 static void* get_type_id_data[1] = {get_type_id_func_tuple};
985
986 static void* get_dimensions_func_tuple[3] = {GEOSGeom_getDimensions_r, (void*)0,
987 (void*)-1};
988 static void* get_dimensions_data[1] = {get_dimensions_func_tuple};
989
990 static void* get_coordinate_dimension_func_tuple[3] = {GEOSGeom_getCoordinateDimension_r,
991 (void*)-1, (void*)-1};
992 static void* get_coordinate_dimension_data[1] = {get_coordinate_dimension_func_tuple};
993
994 static void* get_srid_func_tuple[3] = {GEOSGetSRID_r, (void*)0, (void*)-1};
995 static void* get_srid_data[1] = {get_srid_func_tuple};
996
GetNumPoints(void * context,void * geom,int n)997 static int GetNumPoints(void* context, void* geom, int n) {
998 char typ = GEOSGeomTypeId_r(context, geom);
999 if ((typ == 1) || (typ == 2)) { /* Linestring & Linearring */
1000 return GEOSGeomGetNumPoints_r(context, geom);
1001 } else {
1002 return 0;
1003 }
1004 }
1005 static void* get_num_points_func_tuple[3] = {GetNumPoints, (void*)-1, (void*)0};
1006 static void* get_num_points_data[1] = {get_num_points_func_tuple};
1007
GetNumInteriorRings(void * context,void * geom,int n)1008 static int GetNumInteriorRings(void* context, void* geom, int n) {
1009 char typ = GEOSGeomTypeId_r(context, geom);
1010 if (typ == 3) { /* Polygon */
1011 return GEOSGetNumInteriorRings_r(context, geom);
1012 } else {
1013 return 0;
1014 }
1015 }
1016 static void* get_num_interior_rings_func_tuple[3] = {GetNumInteriorRings, (void*)-1,
1017 (void*)0};
1018 static void* get_num_interior_rings_data[1] = {get_num_interior_rings_func_tuple};
1019
1020 static void* get_num_geometries_func_tuple[3] = {GEOSGetNumGeometries_r, (void*)-1,
1021 (void*)0};
1022 static void* get_num_geometries_data[1] = {get_num_geometries_func_tuple};
1023
1024 static void* get_num_coordinates_func_tuple[3] = {GEOSGetNumCoordinates_r, (void*)-1,
1025 (void*)0};
1026 static void* get_num_coordinates_data[1] = {get_num_coordinates_func_tuple};
1027
1028 typedef int FuncGEOS_Y_i(void* context, void* a);
1029 static char Y_i_dtypes[2] = {NPY_OBJECT, NPY_INT};
Y_i_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1030 static void Y_i_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1031 FuncGEOS_Y_i* func = ((FuncGEOS_Y_i**)data)[0];
1032 int errcode = (int)((int**)data)[1];
1033 int none_value = (int)((int**)data)[2];
1034
1035 GEOSGeometry* in1 = NULL;
1036 int result;
1037
1038 GEOS_INIT_THREADS;
1039
1040 UNARY_LOOP {
1041 /* get the geometry: return on error */
1042 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1043 errstate = PGERR_NOT_A_GEOMETRY;
1044 goto finish;
1045 }
1046 if (in1 == NULL) {
1047 /* None results in 0 for counting functions, -1 otherwise */
1048 *(npy_int*)op1 = none_value;
1049 } else {
1050 result = func(ctx, in1);
1051 // Check last_error if the result equals errcode.
1052 // Otherwise we can't be sure if it is an exception
1053 if ((result == errcode) && (last_error[0] != 0)) {
1054 errstate = PGERR_GEOS_EXCEPTION;
1055 goto finish;
1056 }
1057 *(npy_int*)op1 = result;
1058 }
1059 }
1060
1061 finish:
1062 GEOS_FINISH_THREADS;
1063 }
1064 static PyUFuncGenericFunction Y_i_funcs[1] = {&Y_i_func};
1065
1066 /* Define the geom, geom -> double functions (YY_d) */
1067 static void* distance_data[1] = {GEOSDistance_r};
1068 static void* hausdorff_distance_data[1] = {GEOSHausdorffDistance_r};
1069 #if GEOS_SINCE_3_7_0
GEOSFrechetDistanceWrapped_r(void * context,void * a,void * b,double * c)1070 static int GEOSFrechetDistanceWrapped_r(void* context, void* a, void* b, double* c) {
1071 /* Handle empty geometries (they give segfaults) */
1072 if (GEOSisEmpty_r(context, a) || GEOSisEmpty_r(context, b)) {
1073 *c = NPY_NAN;
1074 return 1;
1075 }
1076 return GEOSFrechetDistance_r(context, a, b, c);
1077 }
1078 static void* frechet_distance_data[1] = {GEOSFrechetDistanceWrapped_r};
1079 #endif
1080 /* Project and ProjectNormalize don't return error codes. wrap them. */
GEOSProjectWrapped_r(void * context,void * a,void * b,double * c)1081 static int GEOSProjectWrapped_r(void* context, void* a, void* b, double* c) {
1082 /* Handle empty points (they give segfaults (for b) or give exception (for a)) */
1083 if (GEOSisEmpty_r(context, a) || GEOSisEmpty_r(context, b)) {
1084 *c = NPY_NAN;
1085 } else {
1086 *c = GEOSProject_r(context, a, b);
1087 }
1088 if (*c == -1.0) {
1089 return 0;
1090 } else {
1091 return 1;
1092 }
1093 }
1094 static void* line_locate_point_data[1] = {GEOSProjectWrapped_r};
GEOSProjectNormalizedWrapped_r(void * context,void * a,void * b,double * c)1095 static int GEOSProjectNormalizedWrapped_r(void* context, void* a, void* b, double* c) {
1096 double length;
1097 double distance;
1098
1099 /* Handle empty points (they give segfaults (for b) or give exception (for a)) */
1100 if (GEOSisEmpty_r(context, a) || GEOSisEmpty_r(context, b)) {
1101 *c = NPY_NAN;
1102 } else {
1103 /* Use custom implementation of GEOSProjectNormalized to overcome bug in
1104 older GEOS versions (https://trac.osgeo.org/geos/ticket/1058) */
1105 if (GEOSLength_r(context, a, &length) != 1) {
1106 return 0;
1107 };
1108 distance = GEOSProject_r(context, a, b);
1109 if (distance == -1.0) {
1110 return 0;
1111 } else {
1112 *c = distance / length;
1113 }
1114 }
1115 return 1;
1116 }
1117 static void* line_locate_point_normalized_data[1] = {GEOSProjectNormalizedWrapped_r};
1118 typedef int FuncGEOS_YY_d(void* context, void* a, void* b, double* c);
1119 static char YY_d_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE};
YY_d_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1120 static void YY_d_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1121 FuncGEOS_YY_d* func = (FuncGEOS_YY_d*)data;
1122 GEOSGeometry *in1 = NULL, *in2 = NULL;
1123
1124 GEOS_INIT_THREADS;
1125
1126 BINARY_LOOP {
1127 /* get the geometries: return on error */
1128 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1129 errstate = PGERR_NOT_A_GEOMETRY;
1130 goto finish;
1131 }
1132 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
1133 errstate = PGERR_NOT_A_GEOMETRY;
1134 goto finish;
1135 }
1136 if ((in1 == NULL) || (in2 == NULL)) {
1137 /* in case of a missing value: return NaN */
1138 *(double*)op1 = NPY_NAN;
1139 } else {
1140 /* let the GEOS function set op1; return on error */
1141 if (func(ctx, in1, in2, (double*)op1) == 0) {
1142 errstate = PGERR_GEOS_EXCEPTION;
1143 goto finish;
1144 }
1145 /* incase the outcome is 0.0, check the inputs for emptyness */
1146 if (*op1 == 0.0) {
1147 if (GEOSisEmpty_r(ctx, in1) || GEOSisEmpty_r(ctx, in2)) {
1148 *(double*)op1 = NPY_NAN;
1149 }
1150 }
1151 }
1152 }
1153
1154 finish:
1155 GEOS_FINISH_THREADS;
1156 }
1157 static PyUFuncGenericFunction YY_d_funcs[1] = {&YY_d_func};
1158
1159 /* Define the geom, geom, double -> double functions (YYd_d) */
1160 static void* hausdorff_distance_densify_data[1] = {GEOSHausdorffDistanceDensify_r};
1161 #if GEOS_SINCE_3_7_0
1162 static void* frechet_distance_densify_data[1] = {GEOSFrechetDistanceDensify_r};
1163 #endif
1164 typedef int FuncGEOS_YYd_d(void* context, void* a, void* b, double c, double* d);
1165 static char YYd_d_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_DOUBLE};
YYd_d_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1166 static void YYd_d_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1167 FuncGEOS_YYd_d* func = (FuncGEOS_YYd_d*)data;
1168 GEOSGeometry *in1 = NULL, *in2 = NULL;
1169
1170 GEOS_INIT_THREADS;
1171
1172 TERNARY_LOOP {
1173 /* get the geometries: return on error */
1174 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1175 errstate = PGERR_NOT_A_GEOMETRY;
1176 goto finish;
1177 }
1178 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
1179 errstate = PGERR_NOT_A_GEOMETRY;
1180 goto finish;
1181 }
1182 double in3 = *(double*)ip3;
1183 if ((in1 == NULL) || (in2 == NULL) || npy_isnan(in3) || GEOSisEmpty_r(ctx, in1) ||
1184 GEOSisEmpty_r(ctx, in2)) {
1185 *(double*)op1 = NPY_NAN;
1186 } else {
1187 /* let the GEOS function set op1; return on error */
1188 if (func(ctx, in1, in2, in3, (double*)op1) == 0) {
1189 errstate = PGERR_GEOS_EXCEPTION;
1190 goto finish;
1191 }
1192 }
1193 }
1194
1195 finish:
1196 GEOS_FINISH_THREADS;
1197 }
1198 static PyUFuncGenericFunction YYd_d_funcs[1] = {&YYd_d_func};
1199
1200 #if GEOS_SINCE_3_9_0
1201
1202 /* Define the geom, geom, double -> geom functions (YYd_Y) */
1203 static void* intersection_prec_data[1] = {GEOSIntersectionPrec_r};
1204 static void* difference_prec_data[1] = {GEOSDifferencePrec_r};
1205 static void* symmetric_difference_prec_data[1] = {GEOSSymDifferencePrec_r};
1206 static void* union_prec_data[1] = {GEOSUnionPrec_r};
1207 typedef void* FuncGEOS_YYd_Y(void* context, void* a, void* b, double c);
1208 static char YYd_Y_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT};
1209
YYd_Y_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1210 static void YYd_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1211 FuncGEOS_YYd_Y* func = (FuncGEOS_YYd_Y*)data;
1212 GEOSGeometry *in1 = NULL, *in2 = NULL;
1213 GEOSGeometry** geom_arr;
1214
1215 // allocate a temporary array to store output GEOSGeometry objects
1216 geom_arr = malloc(sizeof(void*) * dimensions[0]);
1217 CHECK_ALLOC(geom_arr);
1218
1219 GEOS_INIT_THREADS;
1220
1221 TERNARY_LOOP {
1222 // get the geometries: return on error
1223 if (!get_geom(*(GeometryObject**)ip1, &in1) ||
1224 !get_geom(*(GeometryObject**)ip2, &in2)) {
1225 errstate = PGERR_NOT_A_GEOMETRY;
1226 destroy_geom_arr(ctx, geom_arr, i - 1);
1227 break;
1228 }
1229 double in3 = *(double*)ip3;
1230 if ((in1 == NULL) || (in2 == NULL) || npy_isnan(in3)) {
1231 // in case of a missing value: return NULL (None)
1232 geom_arr[i] = NULL;
1233 } else {
1234 geom_arr[i] = func(ctx, in1, in2, in3);
1235 if (geom_arr[i] == NULL) {
1236 errstate = PGERR_GEOS_EXCEPTION;
1237 destroy_geom_arr(ctx, geom_arr, i - 1);
1238 break;
1239 }
1240 }
1241 }
1242
1243 GEOS_FINISH_THREADS;
1244
1245 // fill the numpy array with PyObjects while holding the GIL
1246 if (errstate == PGERR_SUCCESS) {
1247 geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]);
1248 }
1249 free(geom_arr);
1250 }
1251 static PyUFuncGenericFunction YYd_Y_funcs[1] = {&YYd_Y_func};
1252 #endif
1253
1254 /* Define functions with unique call signatures */
1255 static char box_dtypes[6] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
1256 NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
box_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1257 static void box_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1258 char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4];
1259 npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4];
1260 npy_intp n = dimensions[0];
1261 npy_intp i;
1262 GEOSGeometry** geom_arr;
1263
1264 CHECK_NO_INPLACE_OUTPUT(5);
1265
1266 // allocate a temporary array to store output GEOSGeometry objects
1267 geom_arr = malloc(sizeof(void*) * n);
1268 CHECK_ALLOC(geom_arr);
1269
1270 GEOS_INIT_THREADS;
1271
1272 for (i = 0; i < n; i++, ip1 += is1, ip2 += is2, ip3 += is3, ip4 += is4, ip5 += is5) {
1273 geom_arr[i] = create_box(ctx, *(double*)ip1, *(double*)ip2, *(double*)ip3,
1274 *(double*)ip4, *(char*)ip5);
1275 if (geom_arr[i] == NULL) {
1276 // result will be NULL for any nan coordinates, which is OK;
1277 // otherwise raise an error
1278 if (!(npy_isnan(*(double*)ip1) || npy_isnan(*(double*)ip2) ||
1279 npy_isnan(*(double*)ip3) || npy_isnan(*(double*)ip4))) {
1280 errstate = PGERR_GEOS_EXCEPTION;
1281 destroy_geom_arr(ctx, geom_arr, i - 1);
1282 break;
1283 }
1284 }
1285 }
1286
1287 GEOS_FINISH_THREADS;
1288
1289 // fill the numpy array with PyObjects while holding the GIL
1290 if (errstate == PGERR_SUCCESS) {
1291 geom_arr_to_npy(geom_arr, args[5], steps[5], dimensions[0]);
1292 }
1293 free(geom_arr);
1294 }
1295
1296 static PyUFuncGenericFunction box_funcs[1] = {&box_func};
1297
1298 static void* null_data[1] = {NULL};
buffer_inner(void * ctx,GEOSBufferParams * params,void * ip1,void * ip2,GEOSGeometry ** geom_arr,npy_intp i)1299 static char buffer_inner(void* ctx, GEOSBufferParams* params, void* ip1, void* ip2,
1300 GEOSGeometry** geom_arr, npy_intp i) {
1301 GEOSGeometry* in1 = NULL;
1302
1303 /* get the geometry: return on error */
1304 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1305 return PGERR_NOT_A_GEOMETRY;
1306 }
1307 double in2 = *(double*)ip2;
1308 /* handle NULL geometries or NaN buffer width */
1309 if ((in1 == NULL) || npy_isnan(in2)) {
1310 geom_arr[i] = NULL;
1311 } else {
1312 geom_arr[i] = GEOSBufferWithParams_r(ctx, in1, params, in2);
1313 if (geom_arr[i] == NULL) {
1314 return PGERR_GEOS_EXCEPTION;
1315 }
1316 }
1317 return PGERR_SUCCESS;
1318 }
1319
1320 static char buffer_dtypes[8] = {NPY_OBJECT, NPY_DOUBLE, NPY_INT, NPY_INT,
1321 NPY_INT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
buffer_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1322 static void buffer_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1323 char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4],
1324 *ip6 = args[5], *ip7 = args[6];
1325 npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4],
1326 is6 = steps[5], is7 = steps[6];
1327 npy_intp n = dimensions[0];
1328 npy_intp i;
1329 GEOSGeometry** geom_arr;
1330
1331 CHECK_NO_INPLACE_OUTPUT(7);
1332
1333 if ((is3 != 0) || (is4 != 0) || (is5 != 0) || (is6 != 0) || (is7 != 0)) {
1334 PyErr_Format(PyExc_ValueError, "Buffer function called with non-scalar parameters");
1335 return;
1336 }
1337
1338 // allocate a temporary array to store output GEOSGeometry objects
1339 geom_arr = malloc(sizeof(void*) * n);
1340 CHECK_ALLOC(geom_arr);
1341
1342 GEOS_INIT_THREADS;
1343
1344 GEOSBufferParams* params = GEOSBufferParams_create_r(ctx);
1345 if (params != 0) {
1346 if (!GEOSBufferParams_setQuadrantSegments_r(ctx, params, *(int*)ip3)) {
1347 errstate = PGERR_GEOS_EXCEPTION;
1348 }
1349 if (!GEOSBufferParams_setEndCapStyle_r(ctx, params, *(int*)ip4)) {
1350 errstate = PGERR_GEOS_EXCEPTION;
1351 }
1352 if (!GEOSBufferParams_setJoinStyle_r(ctx, params, *(int*)ip5)) {
1353 errstate = PGERR_GEOS_EXCEPTION;
1354 }
1355 if (!GEOSBufferParams_setMitreLimit_r(ctx, params, *(double*)ip6)) {
1356 errstate = PGERR_GEOS_EXCEPTION;
1357 }
1358 if (!GEOSBufferParams_setSingleSided_r(ctx, params, *(npy_bool*)ip7)) {
1359 errstate = PGERR_GEOS_EXCEPTION;
1360 }
1361 } else {
1362 errstate = PGERR_GEOS_EXCEPTION;
1363 }
1364
1365 if (errstate == PGERR_SUCCESS) {
1366 for (i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
1367 errstate = buffer_inner(ctx, params, ip1, ip2, geom_arr, i);
1368 if (errstate != PGERR_SUCCESS) {
1369 destroy_geom_arr(ctx, geom_arr, i - 1);
1370 break;
1371 }
1372 }
1373 }
1374
1375 if (params != 0) {
1376 GEOSBufferParams_destroy_r(ctx, params);
1377 }
1378
1379 GEOS_FINISH_THREADS;
1380
1381 // fill the numpy array with PyObjects while holding the GIL
1382 if (errstate == PGERR_SUCCESS) {
1383 geom_arr_to_npy(geom_arr, args[7], steps[7], dimensions[0]);
1384 }
1385 free(geom_arr);
1386 }
1387 static PyUFuncGenericFunction buffer_funcs[1] = {&buffer_func};
1388
1389 static char offset_curve_dtypes[6] = {NPY_OBJECT, NPY_DOUBLE, NPY_INT,
1390 NPY_INT, NPY_DOUBLE, NPY_OBJECT};
offset_curve_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1391 static void offset_curve_func(char** args, npy_intp* dimensions, npy_intp* steps,
1392 void* data) {
1393 char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4];
1394 npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4];
1395 npy_intp n = dimensions[0];
1396 npy_intp i;
1397 GEOSGeometry** geom_arr;
1398 GEOSGeometry* in1 = NULL;
1399
1400 CHECK_NO_INPLACE_OUTPUT(5);
1401
1402 if ((is3 != 0) || (is4 != 0) || (is5 != 0)) {
1403 PyErr_Format(PyExc_ValueError,
1404 "Offset curve function called with non-scalar parameters");
1405 return;
1406 }
1407
1408 double width;
1409 int quadsegs = *(int*)ip3;
1410 int joinStyle = *(int*)ip4;
1411 double mitreLimit = *(double*)ip5;
1412
1413 // allocate a temporary array to store output GEOSGeometry objects
1414 geom_arr = malloc(sizeof(void*) * n);
1415 CHECK_ALLOC(geom_arr);
1416
1417 GEOS_INIT_THREADS;
1418
1419 for (i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
1420 /* get the geometry: return on error */
1421 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1422 errstate = PGERR_NOT_A_GEOMETRY;
1423 destroy_geom_arr(ctx, geom_arr, i - 1);
1424 break;
1425 }
1426
1427 width = *(double*)ip2;
1428 if ((in1 == NULL) || npy_isnan(width)) {
1429 // in case of a missing value: return NULL (None)
1430 geom_arr[i] = NULL;
1431 } else {
1432 geom_arr[i] = GEOSOffsetCurve_r(ctx, in1, width, quadsegs, joinStyle, mitreLimit);
1433 if (geom_arr[i] == NULL) {
1434 errstate = PGERR_GEOS_EXCEPTION;
1435 destroy_geom_arr(ctx, geom_arr, i - 1);
1436 break;
1437 }
1438 }
1439 }
1440
1441 GEOS_FINISH_THREADS;
1442
1443 // fill the numpy array with PyObjects while holding the GIL
1444 if (errstate == PGERR_SUCCESS) {
1445 geom_arr_to_npy(geom_arr, args[5], steps[5], dimensions[0]);
1446 }
1447 free(geom_arr);
1448 }
1449 static PyUFuncGenericFunction offset_curve_funcs[1] = {&offset_curve_func};
1450
1451 static char snap_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT};
snap_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1452 static void snap_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1453 GEOSGeometry *in1 = NULL, *in2 = NULL;
1454 GEOSGeometry** geom_arr;
1455
1456 CHECK_NO_INPLACE_OUTPUT(3);
1457
1458 // allocate a temporary array to store output GEOSGeometry objects
1459 geom_arr = malloc(sizeof(void*) * dimensions[0]);
1460 CHECK_ALLOC(geom_arr);
1461
1462 GEOS_INIT_THREADS;
1463
1464 TERNARY_LOOP {
1465 /* get the geometries: return on error */
1466 if (!get_geom(*(GeometryObject**)ip1, &in1) ||
1467 !get_geom(*(GeometryObject**)ip2, &in2)) {
1468 errstate = PGERR_NOT_A_GEOMETRY;
1469 destroy_geom_arr(ctx, geom_arr, i - 1);
1470 break;
1471 }
1472 double in3 = *(double*)ip3;
1473 if ((in1 == NULL) || (in2 == NULL) || npy_isnan(in3)) {
1474 // in case of a missing value: return NULL (None)
1475 geom_arr[i] = NULL;
1476 } else {
1477 geom_arr[i] = GEOSSnap_r(ctx, in1, in2, in3);
1478 if (geom_arr[i] == NULL) {
1479 errstate = PGERR_GEOS_EXCEPTION;
1480 destroy_geom_arr(ctx, geom_arr, i - 1);
1481 break;
1482 }
1483 }
1484 }
1485
1486 GEOS_FINISH_THREADS;
1487
1488 // fill the numpy array with PyObjects while holding the GIL
1489 if (errstate == PGERR_SUCCESS) {
1490 geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]);
1491 }
1492 free(geom_arr);
1493 }
1494 static PyUFuncGenericFunction snap_funcs[1] = {&snap_func};
1495
1496 static char clip_by_rect_dtypes[6] = {NPY_OBJECT, NPY_DOUBLE, NPY_DOUBLE,
1497 NPY_DOUBLE, NPY_DOUBLE, NPY_OBJECT};
clip_by_rect_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1498 static void clip_by_rect_func(char** args, npy_intp* dimensions, npy_intp* steps,
1499 void* data) {
1500 char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4];
1501 npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4];
1502 npy_intp n = dimensions[0];
1503 npy_intp i;
1504 GEOSGeometry** geom_arr;
1505 GEOSGeometry* in1 = NULL;
1506
1507 CHECK_NO_INPLACE_OUTPUT(5);
1508
1509 if ((is2 != 0) || (is3 != 0) || (is4 != 0) || (is5 != 0)) {
1510 PyErr_Format(PyExc_ValueError,
1511 "clip_by_rect function called with non-scalar parameters");
1512 return;
1513 }
1514
1515 double xmin = *(double*)ip2;
1516 double ymin = *(double*)ip3;
1517 double xmax = *(double*)ip4;
1518 double ymax = *(double*)ip5;
1519
1520 // allocate a temporary array to store output GEOSGeometry objects
1521 geom_arr = malloc(sizeof(void*) * n);
1522 CHECK_ALLOC(geom_arr);
1523
1524 GEOS_INIT_THREADS;
1525
1526 for (i = 0; i < n; i++, ip1 += is1) {
1527 /* get the geometry: return on error */
1528 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1529 errstate = PGERR_NOT_A_GEOMETRY;
1530 destroy_geom_arr(ctx, geom_arr, i - 1);
1531 break;
1532 }
1533
1534 if (in1 == NULL) {
1535 // in case of a missing value: return NULL (None)
1536 geom_arr[i] = NULL;
1537 } else {
1538 geom_arr[i] = GEOSClipByRect_r(ctx, in1, xmin, ymin, xmax, ymax);
1539 if (geom_arr[i] == NULL) {
1540 errstate = PGERR_GEOS_EXCEPTION;
1541 destroy_geom_arr(ctx, geom_arr, i - 1);
1542 break;
1543 }
1544 }
1545 }
1546
1547 GEOS_FINISH_THREADS;
1548
1549 // fill the numpy array with PyObjects while holding the GIL
1550 if (errstate == PGERR_SUCCESS) {
1551 geom_arr_to_npy(geom_arr, args[5], steps[5], dimensions[0]);
1552 }
1553 free(geom_arr);
1554 }
1555 static PyUFuncGenericFunction clip_by_rect_funcs[1] = {&clip_by_rect_func};
1556
1557 static char equals_exact_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_BOOL};
equals_exact_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1558 static void equals_exact_func(char** args, npy_intp* dimensions, npy_intp* steps,
1559 void* data) {
1560 GEOSGeometry *in1 = NULL, *in2 = NULL;
1561 npy_bool ret;
1562
1563 GEOS_INIT_THREADS;
1564
1565 TERNARY_LOOP {
1566 /* get the geometries: return on error */
1567 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1568 errstate = PGERR_NOT_A_GEOMETRY;
1569 goto finish;
1570 }
1571 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
1572 errstate = PGERR_NOT_A_GEOMETRY;
1573 goto finish;
1574 }
1575 double in3 = *(double*)ip3;
1576 if ((in1 == NULL) || (in2 == NULL) || npy_isnan(in3)) {
1577 /* return 0 (False) for missing values */
1578 ret = 0;
1579 } else {
1580 ret = GEOSEqualsExact_r(ctx, in1, in2, in3);
1581 if ((ret != 0) && (ret != 1)) {
1582 errstate = PGERR_GEOS_EXCEPTION;
1583 goto finish;
1584 }
1585 }
1586 *(npy_bool*)op1 = ret;
1587 }
1588
1589 finish:
1590 GEOS_FINISH_THREADS;
1591 }
1592 static PyUFuncGenericFunction equals_exact_funcs[1] = {&equals_exact_func};
1593
1594 static char delaunay_triangles_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
delaunay_triangles_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1595 static void delaunay_triangles_func(char** args, npy_intp* dimensions, npy_intp* steps,
1596 void* data) {
1597 GEOSGeometry* in1 = NULL;
1598 GEOSGeometry** geom_arr;
1599
1600 CHECK_NO_INPLACE_OUTPUT(3);
1601
1602 // allocate a temporary array to store output GEOSGeometry objects
1603 geom_arr = malloc(sizeof(void*) * dimensions[0]);
1604 CHECK_ALLOC(geom_arr);
1605
1606 GEOS_INIT_THREADS;
1607
1608 TERNARY_LOOP {
1609 // get the geometry: return on error
1610 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1611 errstate = PGERR_NOT_A_GEOMETRY;
1612 destroy_geom_arr(ctx, geom_arr, i - 1);
1613 break;
1614 }
1615 double in2 = *(double*)ip2;
1616 npy_bool in3 = *(npy_bool*)ip3;
1617 if ((in1 == NULL) || npy_isnan(in2)) {
1618 // in case of a missing value: return NULL (None)
1619 geom_arr[i] = NULL;
1620 } else {
1621 geom_arr[i] = GEOSDelaunayTriangulation_r(ctx, in1, in2, (int)in3);
1622 if (geom_arr[i] == NULL) {
1623 errstate = PGERR_GEOS_EXCEPTION;
1624 destroy_geom_arr(ctx, geom_arr, i - 1);
1625 break;
1626 }
1627 }
1628 }
1629
1630 GEOS_FINISH_THREADS;
1631
1632 // fill the numpy array with PyObjects while holding the GIL
1633 if (errstate == PGERR_SUCCESS) {
1634 geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]);
1635 }
1636 free(geom_arr);
1637 }
1638 static PyUFuncGenericFunction delaunay_triangles_funcs[1] = {&delaunay_triangles_func};
1639
1640 static char voronoi_polygons_dtypes[5] = {NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT, NPY_BOOL,
1641 NPY_OBJECT};
voronoi_polygons_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1642 static void voronoi_polygons_func(char** args, npy_intp* dimensions, npy_intp* steps,
1643 void* data) {
1644 GEOSGeometry *in1 = NULL, *in3 = NULL;
1645 GEOSGeometry** geom_arr;
1646
1647 CHECK_NO_INPLACE_OUTPUT(4);
1648
1649 // allocate a temporary array to store output GEOSGeometry objects
1650 geom_arr = malloc(sizeof(void*) * dimensions[0]);
1651 CHECK_ALLOC(geom_arr);
1652
1653 GEOS_INIT_THREADS;
1654
1655 QUATERNARY_LOOP {
1656 // get the geometry: return on error
1657 if (!get_geom(*(GeometryObject**)ip1, &in1) ||
1658 !get_geom(*(GeometryObject**)ip3, &in3)) {
1659 errstate = PGERR_NOT_A_GEOMETRY;
1660 destroy_geom_arr(ctx, geom_arr, i - 1);
1661 break;
1662 }
1663 double in2 = *(double*)ip2;
1664 npy_bool in4 = *(npy_bool*)ip4;
1665 if ((in1 == NULL) || npy_isnan(in2)) {
1666 /* propagate NULL geometries; in3 = NULL is actually supported */
1667 geom_arr[i] = NULL;
1668 } else {
1669 geom_arr[i] = GEOSVoronoiDiagram_r(ctx, in1, in3, in2, (int)in4);
1670 if (geom_arr[i] == NULL) {
1671 errstate = PGERR_GEOS_EXCEPTION;
1672 destroy_geom_arr(ctx, geom_arr, i - 1);
1673 break;
1674 }
1675 }
1676 }
1677
1678 GEOS_FINISH_THREADS;
1679
1680 // fill the numpy array with PyObjects while holding the GIL
1681 if (errstate == PGERR_SUCCESS) {
1682 geom_arr_to_npy(geom_arr, args[4], steps[4], dimensions[0]);
1683 }
1684 free(geom_arr);
1685 }
1686 static PyUFuncGenericFunction voronoi_polygons_funcs[1] = {&voronoi_polygons_func};
1687
1688 static char is_valid_reason_dtypes[2] = {NPY_OBJECT, NPY_OBJECT};
is_valid_reason_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1689 static void is_valid_reason_func(char** args, npy_intp* dimensions, npy_intp* steps,
1690 void* data) {
1691 char* reason;
1692 GEOSGeometry* in1 = NULL;
1693
1694 GEOS_INIT;
1695
1696 UNARY_LOOP {
1697 PyObject** out = (PyObject**)op1;
1698 /* get the geometry return on error */
1699 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1700 errstate = PGERR_NOT_A_GEOMETRY;
1701 goto finish;
1702 }
1703 if (in1 == NULL) {
1704 /* Missing geometries give None */
1705 Py_XDECREF(*out);
1706 Py_INCREF(Py_None);
1707 *out = Py_None;
1708 } else {
1709 reason = GEOSisValidReason_r(ctx, in1);
1710 if (reason == NULL) {
1711 errstate = PGERR_GEOS_EXCEPTION;
1712 goto finish;
1713 }
1714 /* convert to python string and set to out */
1715 Py_XDECREF(*out);
1716 *out = PyUnicode_FromString(reason);
1717 GEOSFree_r(ctx, reason);
1718 }
1719 }
1720
1721 finish:
1722 GEOS_FINISH;
1723 }
1724 static PyUFuncGenericFunction is_valid_reason_funcs[1] = {&is_valid_reason_func};
1725
1726 static char relate_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
relate_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1727 static void relate_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
1728 char* pattern;
1729 GEOSGeometry *in1 = NULL, *in2 = NULL;
1730
1731 GEOS_INIT;
1732
1733 BINARY_LOOP {
1734 PyObject** out = (PyObject**)op1;
1735 /* get the geometries: return on error */
1736 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1737 errstate = PGERR_NOT_A_GEOMETRY;
1738 goto finish;
1739 }
1740 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
1741 errstate = PGERR_NOT_A_GEOMETRY;
1742 goto finish;
1743 }
1744 if ((in1 == NULL) || (in2 == NULL)) {
1745 /* Missing geometries give None */
1746 Py_XDECREF(*out);
1747 Py_INCREF(Py_None);
1748 *out = Py_None;
1749 } else {
1750 pattern = GEOSRelate_r(ctx, in1, in2);
1751 if (pattern == NULL) {
1752 errstate = PGERR_GEOS_EXCEPTION;
1753 goto finish;
1754 }
1755 /* convert to python string and set to out */
1756 Py_XDECREF(*out);
1757 *out = PyUnicode_FromString(pattern);
1758 GEOSFree_r(ctx, pattern);
1759 }
1760 }
1761
1762 finish:
1763 GEOS_FINISH;
1764 }
1765 static PyUFuncGenericFunction relate_funcs[1] = {&relate_func};
1766
1767 static char relate_pattern_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT, NPY_BOOL};
relate_pattern_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1768 static void relate_pattern_func(char** args, npy_intp* dimensions, npy_intp* steps,
1769 void* data) {
1770 GEOSGeometry *in1 = NULL, *in2 = NULL;
1771 const char* pattern = NULL;
1772 npy_bool ret;
1773
1774 /* get the pattern argument (only deal with scalar for now) */
1775 char* ip3 = args[2];
1776 npy_intp is3 = steps[2];
1777
1778 if (is3 != 0) {
1779 PyErr_Format(PyExc_ValueError, "pattern keyword only supports scalar argument");
1780 return;
1781 }
1782 PyObject* in3 = *(PyObject**)ip3;
1783 if (PyUnicode_Check(in3)) {
1784 pattern = PyUnicode_AsUTF8(in3);
1785 if (pattern == NULL) {
1786 /* error happened in PyUnicode_AsUTF8, error already set by Python */
1787 return;
1788 }
1789 } else {
1790 PyErr_Format(PyExc_TypeError, "pattern keyword expected string, got %s",
1791 Py_TYPE(in3)->tp_name);
1792 return;
1793 }
1794
1795 GEOS_INIT_THREADS;
1796
1797 TERNARY_LOOP {
1798 /* get the geometries: return on error */
1799 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
1800 errstate = PGERR_NOT_A_GEOMETRY;
1801 goto finish;
1802 }
1803 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
1804 errstate = PGERR_NOT_A_GEOMETRY;
1805 goto finish;
1806 }
1807 /* ip3 is already handled above */
1808
1809 if ((in1 == NULL) || (in2 == NULL)) {
1810 /* in case of a missing value: return 0 (False) */
1811 ret = 0;
1812 } else {
1813 ret = GEOSRelatePattern_r(ctx, in1, in2, pattern);
1814 if (ret == 2) {
1815 errstate = PGERR_GEOS_EXCEPTION;
1816 goto finish;
1817 }
1818 }
1819 *(npy_bool*)op1 = ret;
1820 }
1821
1822 finish:
1823 GEOS_FINISH_THREADS;
1824 }
1825 static PyUFuncGenericFunction relate_pattern_funcs[1] = {&relate_pattern_func};
1826
1827 static char polygonize_dtypes[2] = {NPY_OBJECT, NPY_OBJECT};
polygonize_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1828 static void polygonize_func(char** args, npy_intp* dimensions, npy_intp* steps,
1829 void* data) {
1830 GEOSGeometry* geom = NULL;
1831 unsigned int n_geoms;
1832
1833 GEOS_INIT;
1834
1835 GEOSGeometry** geoms = malloc(sizeof(void*) * dimensions[1]);
1836 if (geoms == NULL) {
1837 errstate = PGERR_NO_MALLOC;
1838 goto finish;
1839 }
1840
1841 SINGLE_COREDIM_LOOP_OUTER {
1842 n_geoms = 0;
1843 SINGLE_COREDIM_LOOP_INNER {
1844 if (!get_geom(*(GeometryObject**)cp1, &geom)) {
1845 errstate = PGERR_NOT_A_GEOMETRY;
1846 goto finish;
1847 }
1848 if (geom == NULL) {
1849 continue;
1850 }
1851 geoms[n_geoms] = geom;
1852 n_geoms++;
1853 }
1854
1855 GEOSGeometry* ret_ptr = GEOSPolygonize_r(ctx, geoms, n_geoms);
1856 if (ret_ptr == NULL) {
1857 errstate = PGERR_GEOS_EXCEPTION;
1858 goto finish;
1859 }
1860 OUTPUT_Y;
1861 }
1862
1863 finish:
1864 if (geoms != NULL) {
1865 free(geoms);
1866 }
1867 GEOS_FINISH;
1868 }
1869 static PyUFuncGenericFunction polygonize_funcs[1] = {&polygonize_func};
1870
1871 static char polygonize_full_dtypes[5] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT, NPY_OBJECT,
1872 NPY_OBJECT};
polygonize_full_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1873 static void polygonize_full_func(char** args, npy_intp* dimensions, npy_intp* steps,
1874 void* data) {
1875 GEOSGeometry* geom = NULL;
1876 GEOSGeometry* geom_copy = NULL;
1877 unsigned int n_geoms;
1878
1879 GEOSGeometry* collection = NULL;
1880 GEOSGeometry* cuts = NULL;
1881 GEOSGeometry* dangles = NULL;
1882 GEOSGeometry* invalidRings = NULL;
1883
1884 GEOS_INIT;
1885
1886 GEOSGeometry** geoms = malloc(sizeof(void*) * dimensions[1]);
1887 if (geoms == NULL) {
1888 errstate = PGERR_NO_MALLOC;
1889 goto finish;
1890 }
1891
1892 SINGLE_COREDIM_LOOP_OUTER_NOUT4 {
1893 n_geoms = 0;
1894 SINGLE_COREDIM_LOOP_INNER {
1895 if (!get_geom(*(GeometryObject**)cp1, &geom)) {
1896 errstate = PGERR_NOT_A_GEOMETRY;
1897 goto finish;
1898 }
1899 if (geom == NULL) {
1900 continue;
1901 }
1902 // need to copy the input geometries, because the Collection takes ownership
1903 geom_copy = GEOSGeom_clone_r(ctx, geom);
1904 if (geom_copy == NULL) {
1905 // if something went wrong before creating the collection, destroy previously
1906 // cloned geoms
1907 for (i = 0; i < n_geoms; i++) {
1908 GEOSGeom_destroy_r(ctx, geoms[i]);
1909 }
1910 errstate = PGERR_GEOS_EXCEPTION;
1911 goto finish;
1912 }
1913 geoms[n_geoms] = geom_copy;
1914 n_geoms++;
1915 }
1916 collection =
1917 GEOSGeom_createCollection_r(ctx, GEOS_GEOMETRYCOLLECTION, geoms, n_geoms);
1918 if (collection == NULL) {
1919 errstate = PGERR_GEOS_EXCEPTION;
1920 goto finish;
1921 }
1922
1923 GEOSGeometry* ret_ptr =
1924 GEOSPolygonize_full_r(ctx, collection, &cuts, &dangles, &invalidRings);
1925 if (ret_ptr == NULL) {
1926 errstate = PGERR_GEOS_EXCEPTION;
1927 goto finish;
1928 }
1929 OUTPUT_Y_I(1, ret_ptr);
1930 OUTPUT_Y_I(2, cuts);
1931 OUTPUT_Y_I(3, dangles);
1932 OUTPUT_Y_I(4, invalidRings);
1933 GEOSGeom_destroy_r(ctx, collection);
1934 collection = NULL;
1935 }
1936
1937 finish:
1938 if (collection != NULL) {
1939 GEOSGeom_destroy_r(ctx, collection);
1940 }
1941 if (geoms != NULL) {
1942 free(geoms);
1943 }
1944 GEOS_FINISH;
1945 }
1946 static PyUFuncGenericFunction polygonize_full_funcs[1] = {&polygonize_full_func};
1947
1948 static char shortest_line_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
shortest_line_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)1949 static void shortest_line_func(char** args, npy_intp* dimensions, npy_intp* steps,
1950 void* data) {
1951 GEOSGeometry* in1 = NULL;
1952 GEOSGeometry* in2 = NULL;
1953 GEOSPreparedGeometry* in1_prepared = NULL;
1954 GEOSGeometry** geom_arr;
1955 GEOSCoordSequence* coord_seq = NULL;
1956
1957 CHECK_NO_INPLACE_OUTPUT(2);
1958
1959 // allocate a temporary array to store output GEOSGeometry objects
1960 geom_arr = malloc(sizeof(void*) * dimensions[0]);
1961 CHECK_ALLOC(geom_arr);
1962
1963 GEOS_INIT_THREADS;
1964
1965 BINARY_LOOP {
1966 /* get the geometries: return on error */
1967 if (!get_geom_with_prepared(*(GeometryObject**)ip1, &in1, &in1_prepared)) {
1968 errstate = PGERR_NOT_A_GEOMETRY;
1969 destroy_geom_arr(ctx, geom_arr, i - 1);
1970 break;
1971 }
1972 if (!get_geom(*(GeometryObject**)ip2, &in2)) {
1973 errstate = PGERR_NOT_A_GEOMETRY;
1974 destroy_geom_arr(ctx, geom_arr, i - 1);
1975 break;
1976 }
1977
1978 if ((in1 == NULL) || (in2 == NULL) || GEOSisEmpty_r(ctx, in1) ||
1979 GEOSisEmpty_r(ctx, in2)) {
1980 // in case of a missing value or empty geometry: return NULL (None)
1981 // GEOSNearestPoints_r returns NULL for empty geometries
1982 // but this is not distinguishable from an actual error, so we handle this ourselves
1983 geom_arr[i] = NULL;
1984 continue;
1985 }
1986 #if GEOS_SINCE_3_9_0
1987 if (in1_prepared != NULL) {
1988 coord_seq = GEOSPreparedNearestPoints_r(ctx, in1_prepared, in2);
1989 } else {
1990 coord_seq = GEOSNearestPoints_r(ctx, in1, in2);
1991 }
1992 #else
1993 coord_seq = GEOSNearestPoints_r(ctx, in1, in2);
1994 #endif
1995 if (coord_seq == NULL) {
1996 errstate = PGERR_GEOS_EXCEPTION;
1997 destroy_geom_arr(ctx, geom_arr, i - 1);
1998 break;
1999 }
2000 geom_arr[i] = GEOSGeom_createLineString_r(ctx, coord_seq);
2001 // Note: coordinate sequence is owned by linestring; if linestring fails to
2002 // construct, it will automatically clean up the coordinate sequence
2003 if (geom_arr[i] == NULL) {
2004 errstate = PGERR_GEOS_EXCEPTION;
2005 destroy_geom_arr(ctx, geom_arr, i - 1);
2006 break;
2007 }
2008 }
2009
2010 GEOS_FINISH_THREADS;
2011
2012 // fill the numpy array with PyObjects while holding the GIL
2013 if (errstate == PGERR_SUCCESS) {
2014 geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
2015 }
2016 free(geom_arr);
2017 }
2018 static PyUFuncGenericFunction shortest_line_funcs[1] = {&shortest_line_func};
2019
2020 #if GEOS_SINCE_3_6_0
2021 static char set_precision_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
set_precision_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2022 static void set_precision_func(char** args, npy_intp* dimensions, npy_intp* steps,
2023 void* data) {
2024 GEOSGeometry* in1 = NULL;
2025 GEOSGeometry** geom_arr;
2026
2027 CHECK_NO_INPLACE_OUTPUT(3);
2028
2029 // allocate a temporary array to store output GEOSGeometry objects
2030 geom_arr = malloc(sizeof(void*) * dimensions[0]);
2031 CHECK_ALLOC(geom_arr);
2032
2033 GEOS_INIT_THREADS;
2034
2035 QUATERNARY_LOOP {
2036 // get the geometry: return on error
2037 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
2038 errstate = PGERR_NOT_A_GEOMETRY;
2039 destroy_geom_arr(ctx, geom_arr, i - 1);
2040 break;
2041 }
2042 // grid size
2043 double in2 = *(double*)ip2;
2044 // preserve topology
2045 npy_bool in3 = *(npy_bool*)ip3;
2046 // flags:
2047 // GEOS_PREC_NO_TOPO (1<<0): if set, do not try to preserve topology
2048 // GEOS_PREC_KEEP_COLLAPSED (1<<1): Not used because uncollapsed geometries are
2049 // invalid and will not be retained in GEOS >= 3.9 anyway.
2050 int flags = in3 ? 0 : GEOS_PREC_NO_TOPO;
2051
2052 if ((in1 == NULL) || npy_isnan(in2)) {
2053 // in case of a missing value: return NULL (None)
2054 geom_arr[i] = NULL;
2055 } else {
2056 geom_arr[i] = GEOSGeom_setPrecision_r(ctx, in1, in2, flags);
2057 if (geom_arr[i] == NULL) {
2058 errstate = PGERR_GEOS_EXCEPTION;
2059 destroy_geom_arr(ctx, geom_arr, i - 1);
2060 break;
2061 }
2062 }
2063 }
2064
2065 GEOS_FINISH_THREADS;
2066
2067 // fill the numpy array with PyObjects while holding the GIL
2068 if (errstate == PGERR_SUCCESS) {
2069 geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]);
2070 }
2071 free(geom_arr);
2072 }
2073
2074 static PyUFuncGenericFunction set_precision_funcs[1] = {&set_precision_func};
2075 #endif
2076
2077 /* define double -> geometry construction functions */
2078 static char points_dtypes[2] = {NPY_DOUBLE, NPY_OBJECT};
points_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2079 static void points_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
2080 GEOSCoordSequence* coord_seq = NULL;
2081 GEOSGeometry** geom_arr;
2082
2083 // allocate a temporary array to store output GEOSGeometry objects
2084 geom_arr = malloc(sizeof(void*) * dimensions[0]);
2085 CHECK_ALLOC(geom_arr);
2086
2087 GEOS_INIT_THREADS;
2088
2089 SINGLE_COREDIM_LOOP_OUTER {
2090 coord_seq = GEOSCoordSeq_create_r(ctx, 1, n_c1);
2091 if (coord_seq == NULL) {
2092 errstate = PGERR_GEOS_EXCEPTION;
2093 destroy_geom_arr(ctx, geom_arr, i - 1);
2094 goto finish;
2095 }
2096 SINGLE_COREDIM_LOOP_INNER {
2097 if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, 0, i_c1, *(double*)cp1)) {
2098 errstate = PGERR_GEOS_EXCEPTION;
2099 GEOSCoordSeq_destroy_r(ctx, coord_seq);
2100 destroy_geom_arr(ctx, geom_arr, i - 1);
2101 goto finish;
2102 }
2103 }
2104 geom_arr[i] = GEOSGeom_createPoint_r(ctx, coord_seq);
2105 // Note: coordinate sequence is owned by point; if point fails to construct, it will
2106 // automatically clean up the coordinate sequence
2107 if (geom_arr[i] == NULL) {
2108 errstate = PGERR_GEOS_EXCEPTION;
2109 destroy_geom_arr(ctx, geom_arr, i - 1);
2110 goto finish;
2111 }
2112 }
2113
2114 finish:
2115 GEOS_FINISH_THREADS;
2116 // fill the numpy array with PyObjects while holding the GIL
2117 if (errstate == PGERR_SUCCESS) {
2118 geom_arr_to_npy(geom_arr, args[1], steps[1], dimensions[0]);
2119 }
2120 free(geom_arr);
2121 }
2122 static PyUFuncGenericFunction points_funcs[1] = {&points_func};
2123
2124 static char linestrings_dtypes[2] = {NPY_DOUBLE, NPY_OBJECT};
linestrings_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2125 static void linestrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
2126 void* data) {
2127 GEOSCoordSequence* coord_seq = NULL;
2128 GEOSGeometry** geom_arr;
2129
2130 // allocate a temporary array to store output GEOSGeometry objects
2131 geom_arr = malloc(sizeof(void*) * dimensions[0]);
2132 CHECK_ALLOC(geom_arr);
2133
2134 GEOS_INIT_THREADS;
2135
2136 DOUBLE_COREDIM_LOOP_OUTER {
2137 coord_seq = GEOSCoordSeq_create_r(ctx, n_c1, n_c2);
2138 if (coord_seq == NULL) {
2139 errstate = PGERR_GEOS_EXCEPTION;
2140 destroy_geom_arr(ctx, geom_arr, i - 1);
2141 goto finish;
2142 }
2143 DOUBLE_COREDIM_LOOP_INNER_1 {
2144 DOUBLE_COREDIM_LOOP_INNER_2 {
2145 if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, i_c1, i_c2, *(double*)cp2)) {
2146 errstate = PGERR_GEOS_EXCEPTION;
2147 GEOSCoordSeq_destroy_r(ctx, coord_seq);
2148 destroy_geom_arr(ctx, geom_arr, i - 1);
2149 goto finish;
2150 }
2151 }
2152 }
2153 geom_arr[i] = GEOSGeom_createLineString_r(ctx, coord_seq);
2154 // Note: coordinate sequence is owned by linestring; if linestring fails to construct,
2155 // it will automatically clean up the coordinate sequence
2156 if (geom_arr[i] == NULL) {
2157 errstate = PGERR_GEOS_EXCEPTION;
2158 destroy_geom_arr(ctx, geom_arr, i - 1);
2159 goto finish;
2160 }
2161 }
2162
2163 finish:
2164 GEOS_FINISH_THREADS;
2165
2166 // fill the numpy array with PyObjects while holding the GIL
2167 if (errstate == PGERR_SUCCESS) {
2168 geom_arr_to_npy(geom_arr, args[1], steps[1], dimensions[0]);
2169 }
2170 free(geom_arr);
2171 }
2172 static PyUFuncGenericFunction linestrings_funcs[1] = {&linestrings_func};
2173
2174 static char linearrings_dtypes[2] = {NPY_DOUBLE, NPY_OBJECT};
linearrings_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2175 static void linearrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
2176 void* data) {
2177 GEOSCoordSequence* coord_seq = NULL;
2178 GEOSGeometry** geom_arr;
2179 char ring_closure = 0;
2180 double first_coord, last_coord;
2181
2182 // allocate a temporary array to store output GEOSGeometry objects
2183 geom_arr = malloc(sizeof(void*) * dimensions[0]);
2184 CHECK_ALLOC(geom_arr);
2185
2186 GEOS_INIT_THREADS;
2187
2188 DOUBLE_COREDIM_LOOP_OUTER {
2189 /* check if first and last coords are equal; duplicate if necessary */
2190 ring_closure = 0;
2191 DOUBLE_COREDIM_LOOP_INNER_2 {
2192 first_coord = *(double*)(ip1 + i_c2 * cs2);
2193 last_coord = *(double*)(ip1 + (n_c1 - 1) * cs1 + i_c2 * cs2);
2194 if (first_coord != last_coord) {
2195 ring_closure = 1;
2196 break;
2197 }
2198 }
2199 /* fill the coordinate sequence */
2200 coord_seq = GEOSCoordSeq_create_r(ctx, n_c1 + ring_closure, n_c2);
2201 if (coord_seq == NULL) {
2202 errstate = PGERR_GEOS_EXCEPTION;
2203 destroy_geom_arr(ctx, geom_arr, i - 1);
2204 goto finish;
2205 }
2206 DOUBLE_COREDIM_LOOP_INNER_1 {
2207 DOUBLE_COREDIM_LOOP_INNER_2 {
2208 if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, i_c1, i_c2, *(double*)cp2)) {
2209 errstate = PGERR_GEOS_EXCEPTION;
2210 GEOSCoordSeq_destroy_r(ctx, coord_seq);
2211 destroy_geom_arr(ctx, geom_arr, i - 1);
2212 goto finish;
2213 }
2214 }
2215 }
2216 /* add the closing coordinate if necessary */
2217 if (ring_closure) {
2218 DOUBLE_COREDIM_LOOP_INNER_2 {
2219 first_coord = *(double*)(ip1 + i_c2 * cs2);
2220 if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, n_c1, i_c2, first_coord)) {
2221 errstate = PGERR_GEOS_EXCEPTION;
2222 GEOSCoordSeq_destroy_r(ctx, coord_seq);
2223 destroy_geom_arr(ctx, geom_arr, i - 1);
2224 goto finish;
2225 }
2226 }
2227 }
2228 geom_arr[i] = GEOSGeom_createLinearRing_r(ctx, coord_seq);
2229 // Note: coordinate sequence is owned by linearring; if linearring fails to construct,
2230 // it will automatically clean up the coordinate sequence
2231 if (geom_arr[i] == NULL) {
2232 errstate = PGERR_GEOS_EXCEPTION;
2233 destroy_geom_arr(ctx, geom_arr, i - 1);
2234 goto finish;
2235 }
2236 }
2237
2238 finish:
2239 GEOS_FINISH_THREADS;
2240
2241 // fill the numpy array with PyObjects while holding the GIL
2242 if (errstate == PGERR_SUCCESS) {
2243 geom_arr_to_npy(geom_arr, args[1], steps[1], dimensions[0]);
2244 }
2245 free(geom_arr);
2246 }
2247 static PyUFuncGenericFunction linearrings_funcs[1] = {&linearrings_func};
2248
2249 static char polygons_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
polygons_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2250 static void polygons_func(char** args, npy_intp* dimensions, npy_intp* steps,
2251 void* data) {
2252 GEOSGeometry *hole, *shell, *hole_copy, *shell_copy;
2253 GEOSGeometry **holes, **geom_arr;
2254 int geom_type;
2255 int n_holes;
2256
2257 // allocate a temporary array to store output GEOSGeometry objects
2258 geom_arr = malloc(sizeof(void*) * dimensions[0]);
2259 CHECK_ALLOC(geom_arr)
2260
2261 // allocate a temporary array to store holes
2262 holes = malloc(sizeof(void*) * dimensions[1]);
2263 CHECK_ALLOC(holes)
2264
2265 GEOS_INIT_THREADS;
2266
2267 BINARY_SINGLE_COREDIM_LOOP_OUTER {
2268 if (!get_geom(*(GeometryObject**)ip1, &shell)) {
2269 errstate = PGERR_NOT_A_GEOMETRY;
2270 destroy_geom_arr(ctx, geom_arr, i - 1);
2271 break;
2272 }
2273 if (shell == NULL) {
2274 // output empty polygon if shell is None (ignoring holes)
2275 geom_arr[i] = GEOSGeom_createEmptyPolygon_r(ctx);
2276 if (geom_arr[i] == NULL) {
2277 errstate = PGERR_GEOS_EXCEPTION;
2278 destroy_geom_arr(ctx, geom_arr, i - 1);
2279 break;
2280 };
2281 continue;
2282 }
2283 geom_type = GEOSGeomTypeId_r(ctx, shell);
2284 // Pre-emptively check the geometry type (https://trac.osgeo.org/geos/ticket/1111)
2285 if (geom_type == -1) {
2286 errstate = PGERR_GEOS_EXCEPTION;
2287 destroy_geom_arr(ctx, geom_arr, i - 1);
2288 break;
2289 } else if (geom_type != 2) {
2290 errstate = PGERR_GEOMETRY_TYPE;
2291 destroy_geom_arr(ctx, geom_arr, i - 1);
2292 break;
2293 }
2294 n_holes = 0;
2295 cp1 = ip2;
2296 BINARY_SINGLE_COREDIM_LOOP_INNER {
2297 if (!get_geom(*(GeometryObject**)cp1, &hole)) {
2298 errstate = PGERR_NOT_A_GEOMETRY;
2299 destroy_geom_arr(ctx, geom_arr, i - 1);
2300 destroy_geom_arr(ctx, holes, n_holes - 1);
2301 goto finish;
2302 }
2303 if (hole == NULL) {
2304 continue;
2305 }
2306 // Pre-emptively check the geometry type (https://trac.osgeo.org/geos/ticket/1111)
2307 geom_type = GEOSGeomTypeId_r(ctx, hole);
2308 if (geom_type == -1) {
2309 errstate = PGERR_GEOS_EXCEPTION;
2310 destroy_geom_arr(ctx, geom_arr, i - 1);
2311 destroy_geom_arr(ctx, holes, n_holes - 1);
2312 goto finish;
2313 } else if (geom_type != 2) {
2314 errstate = PGERR_GEOMETRY_TYPE;
2315 destroy_geom_arr(ctx, geom_arr, i - 1);
2316 destroy_geom_arr(ctx, holes, n_holes - 1);
2317 goto finish;
2318 }
2319 hole_copy = GEOSGeom_clone_r(ctx, hole);
2320 if (hole_copy == NULL) {
2321 errstate = PGERR_GEOS_EXCEPTION;
2322 destroy_geom_arr(ctx, geom_arr, i - 1);
2323 destroy_geom_arr(ctx, holes, n_holes - 1);
2324 goto finish;
2325 }
2326 holes[n_holes] = hole_copy;
2327 n_holes++;
2328 }
2329 shell_copy = GEOSGeom_clone_r(ctx, shell);
2330 if (shell_copy == NULL) {
2331 errstate = PGERR_GEOS_EXCEPTION;
2332 destroy_geom_arr(ctx, geom_arr, i - 1);
2333 destroy_geom_arr(ctx, holes, n_holes - 1);
2334 break;
2335 }
2336 geom_arr[i] = GEOSGeom_createPolygon_r(ctx, shell_copy, holes, n_holes);
2337 if (geom_arr[i] == NULL) {
2338 // We will have a memory leak now (https://trac.osgeo.org/geos/ticket/1111)
2339 // but we have covered all known cases that GEOS would error by pre-emptively
2340 // checking if all inputs are linearrings.
2341 errstate = PGERR_GEOS_EXCEPTION;
2342 destroy_geom_arr(ctx, geom_arr, i - 1);
2343 break;
2344 };
2345 }
2346
2347 finish:
2348 GEOS_FINISH_THREADS;
2349
2350 // fill the numpy array with PyObjects while holding the GIL
2351 if (errstate == PGERR_SUCCESS) {
2352 geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
2353 }
2354 free(geom_arr);
2355 if (holes != NULL) {
2356 free(holes);
2357 }
2358 }
2359 static PyUFuncGenericFunction polygons_funcs[1] = {&polygons_func};
2360
2361 static char create_collection_dtypes[3] = {NPY_OBJECT, NPY_INT, NPY_OBJECT};
create_collection_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2362 static void create_collection_func(char** args, npy_intp* dimensions, npy_intp* steps,
2363 void* data) {
2364 GEOSGeometry *g, *g_copy;
2365 int n_geoms, type, actual_type, expected_type, alt_expected_type;
2366 GEOSGeometry **temp_geoms, **geom_arr;
2367
2368 // allocate a temporary array to store output GEOSGeometry objects
2369 geom_arr = malloc(sizeof(void*) * dimensions[0]);
2370 CHECK_ALLOC(geom_arr)
2371
2372 // allocate a temporary array to store geometries to put in a collection
2373 temp_geoms = malloc(sizeof(void*) * dimensions[1]);
2374 CHECK_ALLOC(temp_geoms)
2375
2376 GEOS_INIT_THREADS;
2377
2378 BINARY_SINGLE_COREDIM_LOOP_OUTER {
2379 type = *(int*)ip2;
2380 switch (type) {
2381 case GEOS_MULTIPOINT:
2382 expected_type = GEOS_POINT;
2383 alt_expected_type = -1;
2384 break;
2385 case GEOS_MULTILINESTRING:
2386 expected_type = GEOS_LINESTRING;
2387 alt_expected_type = GEOS_LINEARRING;
2388 break;
2389 case GEOS_MULTIPOLYGON:
2390 expected_type = GEOS_POLYGON;
2391 alt_expected_type = -1;
2392 break;
2393 case GEOS_GEOMETRYCOLLECTION:
2394 expected_type = -1;
2395 alt_expected_type = -1;
2396 break;
2397 default:
2398 errstate = PGERR_GEOMETRY_TYPE;
2399 destroy_geom_arr(ctx, geom_arr, i - 1);
2400 goto finish;
2401 }
2402 n_geoms = 0;
2403 cp1 = ip1;
2404 BINARY_SINGLE_COREDIM_LOOP_INNER {
2405 if (!get_geom(*(GeometryObject**)cp1, &g)) {
2406 errstate = PGERR_NOT_A_GEOMETRY;
2407 destroy_geom_arr(ctx, geom_arr, i - 1);
2408 destroy_geom_arr(ctx, temp_geoms, n_geoms - 1);
2409 goto finish;
2410 }
2411 if (g == NULL) {
2412 continue;
2413 }
2414 if (expected_type != -1) {
2415 actual_type = GEOSGeomTypeId_r(ctx, g);
2416 if (actual_type == -1) {
2417 errstate = PGERR_GEOS_EXCEPTION;
2418 destroy_geom_arr(ctx, geom_arr, i - 1);
2419 destroy_geom_arr(ctx, temp_geoms, n_geoms - 1);
2420 goto finish;
2421 }
2422 if ((actual_type != expected_type) && (actual_type != alt_expected_type)) {
2423 errstate = PGERR_GEOMETRY_TYPE;
2424 destroy_geom_arr(ctx, geom_arr, i - 1);
2425 destroy_geom_arr(ctx, temp_geoms, n_geoms - 1);
2426 goto finish;
2427 }
2428 }
2429 g_copy = GEOSGeom_clone_r(ctx, g);
2430 if (g_copy == NULL) {
2431 errstate = PGERR_GEOS_EXCEPTION;
2432 destroy_geom_arr(ctx, geom_arr, i - 1);
2433 destroy_geom_arr(ctx, temp_geoms, n_geoms - 1);
2434 goto finish;
2435 }
2436 temp_geoms[n_geoms] = g_copy;
2437 n_geoms++;
2438 }
2439 geom_arr[i] = GEOSGeom_createCollection_r(ctx, type, temp_geoms, n_geoms);
2440 if (geom_arr[i] == NULL) {
2441 // We may have a memory leak now (https://trac.osgeo.org/geos/ticket/1111)
2442 // but we have covered all known cases that GEOS would error by pre-emptively
2443 // checking if all inputs are the correct geometry types.
2444 errstate = PGERR_GEOS_EXCEPTION;
2445 destroy_geom_arr(ctx, geom_arr, i - 1);
2446 break;
2447 };
2448 }
2449
2450 finish:
2451 GEOS_FINISH_THREADS;
2452
2453 // fill the numpy array with PyObjects while holding the GIL
2454 if (errstate == PGERR_SUCCESS) {
2455 geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
2456 }
2457 free(geom_arr);
2458 if (temp_geoms != NULL) {
2459 free(temp_geoms);
2460 }
2461 }
2462 static PyUFuncGenericFunction create_collection_funcs[1] = {&create_collection_func};
2463
2464 static char bounds_dtypes[2] = {NPY_OBJECT, NPY_DOUBLE};
bounds_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2465 static void bounds_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
2466 GEOSGeometry *envelope = NULL, *in1;
2467 const GEOSGeometry* ring;
2468 const GEOSCoordSequence* coord_seq;
2469 int size;
2470 char *ip1 = args[0], *op1 = args[1];
2471 double *x1, *y1, *x2, *y2;
2472
2473 GEOS_INIT_THREADS;
2474
2475 npy_intp is1 = steps[0], os1 = steps[1], cs1 = steps[2];
2476 npy_intp n = dimensions[0], i;
2477 for (i = 0; i < n; i++, ip1 += is1, op1 += os1) {
2478 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
2479 errstate = PGERR_NOT_A_GEOMETRY;
2480 goto finish;
2481 }
2482
2483 /* get the 4 (pointers to) the bbox values from the "core stride 1" (cs1) */
2484 x1 = (double*)(op1);
2485 y1 = (double*)(op1 + cs1);
2486 x2 = (double*)(op1 + 2 * cs1);
2487 y2 = (double*)(op1 + 3 * cs1);
2488
2489 if (in1 == NULL) { /* no geometry => bbox becomes (nan, nan, nan, nan) */
2490 *x1 = *y1 = *x2 = *y2 = NPY_NAN;
2491 } else {
2492 /* construct the envelope */
2493 envelope = GEOSEnvelope_r(ctx, in1);
2494 if (envelope == NULL) {
2495 errstate = PGERR_GEOS_EXCEPTION;
2496 goto finish;
2497 }
2498 size = GEOSGetNumCoordinates_r(ctx, envelope);
2499
2500 /* get the bbox depending on the number of coordinates in the envelope */
2501 if (size == 0) { /* Envelope is empty */
2502 *x1 = *y1 = *x2 = *y2 = NPY_NAN;
2503 } else if (size == 1) { /* Envelope is a point */
2504 if (!GEOSGeomGetX_r(ctx, envelope, x1)) {
2505 errstate = PGERR_GEOS_EXCEPTION;
2506 goto finish;
2507 }
2508 if (!GEOSGeomGetY_r(ctx, envelope, y1)) {
2509 errstate = PGERR_GEOS_EXCEPTION;
2510 goto finish;
2511 }
2512 *x2 = *x1;
2513 *y2 = *y1;
2514 } else if (size == 5) { /* Envelope is a box */
2515 ring = GEOSGetExteriorRing_r(ctx, envelope);
2516 if (ring == NULL) {
2517 errstate = PGERR_GEOS_EXCEPTION;
2518 goto finish;
2519 }
2520 coord_seq = GEOSGeom_getCoordSeq_r(ctx, ring);
2521 if (coord_seq == NULL) {
2522 errstate = PGERR_GEOS_EXCEPTION;
2523 goto finish;
2524 }
2525 if (!GEOSCoordSeq_getX_r(ctx, coord_seq, 0, x1)) {
2526 errstate = PGERR_GEOS_EXCEPTION;
2527 goto finish;
2528 }
2529 if (!GEOSCoordSeq_getY_r(ctx, coord_seq, 0, y1)) {
2530 errstate = PGERR_GEOS_EXCEPTION;
2531 goto finish;
2532 }
2533 if (!GEOSCoordSeq_getX_r(ctx, coord_seq, 2, x2)) {
2534 errstate = PGERR_GEOS_EXCEPTION;
2535 goto finish;
2536 }
2537 if (!GEOSCoordSeq_getY_r(ctx, coord_seq, 2, y2)) {
2538 errstate = PGERR_GEOS_EXCEPTION;
2539 goto finish;
2540 }
2541 } else {
2542 errstate = PGERR_GEOMETRY_TYPE;
2543 goto finish;
2544 }
2545 GEOSGeom_destroy_r(ctx, envelope);
2546 envelope = NULL;
2547 }
2548 }
2549
2550 finish:
2551 if (envelope != NULL) {
2552 GEOSGeom_destroy_r(ctx, envelope);
2553 }
2554 GEOS_FINISH_THREADS;
2555 }
2556 static PyUFuncGenericFunction bounds_funcs[1] = {&bounds_func};
2557
2558 /* Define the object -> geom functions (O_Y) */
2559
2560 static char from_wkb_dtypes[3] = {NPY_OBJECT, NPY_UINT8, NPY_OBJECT};
from_wkb_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2561 static void from_wkb_func(char** args, npy_intp* dimensions, npy_intp* steps,
2562 void* data) {
2563 char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];
2564 npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];
2565 PyObject* in1;
2566 npy_uint8 on_invalid = *(npy_uint8*)ip2;
2567 npy_intp n = dimensions[0];
2568 npy_intp i;
2569 GEOSWKBReader* reader;
2570 unsigned char* wkb;
2571 GEOSGeometry* ret_ptr;
2572 Py_ssize_t size;
2573 char is_hex;
2574
2575 if ((is2 != 0)) {
2576 PyErr_Format(PyExc_ValueError, "from_wkb function called with non-scalar parameters");
2577 return;
2578 }
2579
2580 GEOS_INIT;
2581
2582 /* Create the WKB reader */
2583 reader = GEOSWKBReader_create_r(ctx);
2584 if (reader == NULL) {
2585 errstate = PGERR_GEOS_EXCEPTION;
2586 goto finish;
2587 }
2588
2589 for (i = 0; i < n; i++, ip1 += is1, op1 += os1) {
2590 /* ip1 is pointer to array element PyObject* */
2591 in1 = *(PyObject**)ip1;
2592
2593 if (in1 == Py_None) {
2594 /* None in the input propagates to the output */
2595 ret_ptr = NULL;
2596 } else {
2597 /* Cast the PyObject (only bytes) to char* */
2598 if (PyBytes_Check(in1)) {
2599 size = PyBytes_Size(in1);
2600 wkb = (unsigned char*)PyBytes_AsString(in1);
2601 if (wkb == NULL) {
2602 errstate = PGERR_GEOS_EXCEPTION;
2603 goto finish;
2604 }
2605 } else if (PyUnicode_Check(in1)) {
2606 wkb = (unsigned char*)PyUnicode_AsUTF8AndSize(in1, &size);
2607 if (wkb == NULL) {
2608 errstate = PGERR_GEOS_EXCEPTION;
2609 goto finish;
2610 }
2611 } else {
2612 PyErr_Format(PyExc_TypeError, "Expected bytes, got %s", Py_TYPE(in1)->tp_name);
2613 goto finish;
2614 }
2615
2616 /* Check if this is a HEX WKB */
2617 if (size != 0) {
2618 is_hex = ((wkb[0] == 48) || (wkb[0] == 49));
2619 } else {
2620 is_hex = 0;
2621 }
2622
2623 /* Read the WKB */
2624 if (is_hex) {
2625 ret_ptr = GEOSWKBReader_readHEX_r(ctx, reader, wkb, size);
2626 } else {
2627 ret_ptr = GEOSWKBReader_read_r(ctx, reader, wkb, size);
2628 }
2629 if (ret_ptr == NULL) {
2630 if (on_invalid == 2) {
2631 // raise exception
2632 errstate = PGERR_GEOS_EXCEPTION;
2633 goto finish;
2634 } else if (on_invalid == 1) {
2635 // raise warning, return None
2636 errstate = PGWARN_INVALID_WKB;
2637 }
2638 // else: return None, no warning
2639 }
2640 }
2641 OUTPUT_Y;
2642 }
2643
2644 finish:
2645 GEOSWKBReader_destroy_r(ctx, reader);
2646 GEOS_FINISH;
2647 }
2648 static PyUFuncGenericFunction from_wkb_funcs[1] = {&from_wkb_func};
2649
2650 static char from_wkt_dtypes[3] = {NPY_OBJECT, NPY_UINT8, NPY_OBJECT};
from_wkt_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2651 static void from_wkt_func(char** args, npy_intp* dimensions, npy_intp* steps,
2652 void* data) {
2653 char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];
2654 npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];
2655 PyObject* in1;
2656 npy_uint8 on_invalid = *(npy_uint8*)ip2;
2657 npy_intp n = dimensions[0];
2658 npy_intp i;
2659 GEOSGeometry* ret_ptr;
2660 GEOSWKTReader* reader;
2661 const char* wkt;
2662
2663 if ((is2 != 0)) {
2664 PyErr_Format(PyExc_ValueError, "from_wkt function called with non-scalar parameters");
2665 return;
2666 }
2667
2668 GEOS_INIT;
2669
2670 /* Create the WKT reader */
2671 reader = GEOSWKTReader_create_r(ctx);
2672 if (reader == NULL) {
2673 errstate = PGERR_GEOS_EXCEPTION;
2674 goto finish;
2675 }
2676
2677 for (i = 0; i < n; i++, ip1 += is1, op1 += os1) {
2678 /* ip1 is pointer to array element PyObject* */
2679 in1 = *(PyObject**)ip1;
2680
2681 if (in1 == Py_None) {
2682 /* None in the input propagates to the output */
2683 ret_ptr = NULL;
2684 } else {
2685 /* Cast the PyObject (bytes or str) to char* */
2686 if (PyBytes_Check(in1)) {
2687 wkt = PyBytes_AsString(in1);
2688 if (wkt == NULL) {
2689 errstate = PGERR_GEOS_EXCEPTION;
2690 goto finish;
2691 }
2692 } else if (PyUnicode_Check(in1)) {
2693 wkt = PyUnicode_AsUTF8(in1);
2694 if (wkt == NULL) {
2695 errstate = PGERR_GEOS_EXCEPTION;
2696 goto finish;
2697 }
2698 } else {
2699 PyErr_Format(PyExc_TypeError, "Expected bytes, got %s", Py_TYPE(in1)->tp_name);
2700 goto finish;
2701 }
2702
2703 /* Read the WKT */
2704 ret_ptr = GEOSWKTReader_read_r(ctx, reader, wkt);
2705 if (ret_ptr == NULL) {
2706 if (on_invalid == 2) {
2707 // raise exception
2708 errstate = PGERR_GEOS_EXCEPTION;
2709 goto finish;
2710 } else if (on_invalid == 1) {
2711 // raise warning, return None
2712 errstate = PGWARN_INVALID_WKT;
2713 }
2714 // else: return None, no warning
2715 }
2716 }
2717 OUTPUT_Y;
2718 }
2719
2720 finish:
2721 GEOSWKTReader_destroy_r(ctx, reader);
2722 GEOS_FINISH;
2723 }
2724 static PyUFuncGenericFunction from_wkt_funcs[1] = {&from_wkt_func};
2725
2726 static char from_shapely_dtypes[2] = {NPY_OBJECT, NPY_OBJECT};
from_shapely_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2727 static void from_shapely_func(char** args, npy_intp* dimensions, npy_intp* steps,
2728 void* data) {
2729 GEOSGeometry *in_ptr, *ret_ptr;
2730 PyObject *in1, *attr;
2731
2732 GEOS_INIT;
2733
2734 UNARY_LOOP {
2735 /* ip1 is pointer to array element PyObject* */
2736 in1 = *(PyObject**)ip1;
2737
2738 if (in1 == Py_None) {
2739 /* None in the input propagates to the output */
2740 ret_ptr = NULL;
2741 } else {
2742 /* Check if we have a prepared geometry */
2743 if (PyObject_HasAttrString(in1, "context")) {
2744 in1 = PyObject_GetAttrString(in1, "context");
2745 }
2746 /* Get the __geom__ attribute */
2747 attr = PyObject_GetAttrString(in1, "__geom__");
2748 if (attr == NULL) {
2749 /* Raise if __geom__ does not exist */
2750 PyErr_Format(PyExc_TypeError, "Expected a shapely object or None, got %s",
2751 Py_TYPE(in1)->tp_name);
2752 GEOS_FINISH;
2753 return;
2754 } else if (!PyLong_Check(attr)) {
2755 /* Raise if __geom__ is of incorrect type */
2756 PyErr_Format(PyExc_TypeError, "Expected int for the __geom__ attribute, got %s",
2757 Py_TYPE(attr)->tp_name);
2758 Py_XDECREF(attr);
2759 GEOS_FINISH;
2760 return;
2761 }
2762 /* Convert it to a GEOSGeometry pointer */
2763 in_ptr = PyLong_AsVoidPtr(attr);
2764 Py_XDECREF(attr);
2765 /* Clone the geometry and finish */
2766 ret_ptr = GEOSGeom_clone_r(ctx, in_ptr);
2767 if (ret_ptr == NULL) {
2768 errstate = PGERR_GEOS_EXCEPTION;
2769 goto finish;
2770 }
2771 }
2772 OUTPUT_Y;
2773 }
2774
2775 finish:
2776 GEOS_FINISH;
2777 }
2778 static PyUFuncGenericFunction from_shapely_funcs[1] = {&from_shapely_func};
2779
2780 static char to_wkb_dtypes[6] = {NPY_OBJECT, NPY_BOOL, NPY_INT,
2781 NPY_INT, NPY_BOOL, NPY_OBJECT};
to_wkb_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2782 static void to_wkb_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
2783 char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4],
2784 *op1 = args[5];
2785 npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4],
2786 os1 = steps[5];
2787 npy_intp n = dimensions[0];
2788 npy_intp i;
2789
2790 GEOSGeometry *in1, *temp_geom;
2791 GEOSWKBWriter* writer;
2792 unsigned char* wkb;
2793 size_t size;
2794 #if !GEOS_SINCE_3_10_0
2795 char has_empty;
2796 #endif // !GEOS_SINCE_3_10_0
2797
2798 if ((is2 != 0) || (is3 != 0) || (is4 != 0) || (is5 != 0)) {
2799 PyErr_Format(PyExc_ValueError, "to_wkb function called with non-scalar parameters");
2800 return;
2801 }
2802
2803 GEOS_INIT;
2804
2805 /* Create the WKB writer */
2806 writer = GEOSWKBWriter_create_r(ctx);
2807 if (writer == NULL) {
2808 errstate = PGERR_GEOS_EXCEPTION;
2809 goto finish;
2810 }
2811
2812 char hex = *(npy_bool*)ip2;
2813 GEOSWKBWriter_setOutputDimension_r(ctx, writer, *(int*)ip3);
2814 int byte_order = *(int*)ip4;
2815 if (byte_order != -1) {
2816 GEOSWKBWriter_setByteOrder_r(ctx, writer, *(int*)ip4);
2817 }
2818 GEOSWKBWriter_setIncludeSRID_r(ctx, writer, *(npy_bool*)ip5);
2819
2820 // Check if the above functions caused a GEOS exception
2821 if (last_error[0] != 0) {
2822 errstate = PGERR_GEOS_EXCEPTION;
2823 goto finish;
2824 }
2825
2826 for (i = 0; i < n; i++, ip1 += is1, op1 += os1) {
2827 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
2828 errstate = PGERR_NOT_A_GEOMETRY;
2829 goto finish;
2830 }
2831 PyObject** out = (PyObject**)op1;
2832
2833 if (in1 == NULL) {
2834 Py_XDECREF(*out);
2835 Py_INCREF(Py_None);
2836 *out = Py_None;
2837 } else {
2838 #if !GEOS_SINCE_3_10_0
2839 // WKB Does not allow empty points in GEOS<3.10.
2840 // We check for that and patch the POINT EMPTY if necessary
2841 has_empty = has_point_empty(ctx, in1);
2842 if (has_empty == 2) {
2843 errstate = PGERR_GEOS_EXCEPTION;
2844 goto finish;
2845 }
2846 if (has_empty) {
2847 temp_geom = point_empty_to_nan_all_geoms(ctx, in1);
2848 } else {
2849 temp_geom = in1;
2850 }
2851 #else
2852 temp_geom = in1;
2853 #endif // !GEOS_SINCE_3_10_0
2854 if (hex) {
2855 wkb = GEOSWKBWriter_writeHEX_r(ctx, writer, temp_geom, &size);
2856 } else {
2857 wkb = GEOSWKBWriter_write_r(ctx, writer, temp_geom, &size);
2858 }
2859 #if !GEOS_SINCE_3_10_0
2860 // Destroy the temp_geom if it was patched (POINT EMPTY patch)
2861 if (has_empty) {
2862 GEOSGeom_destroy_r(ctx, temp_geom);
2863 }
2864 #endif // !GEOS_SINCE_3_10_0
2865 if (wkb == NULL) {
2866 errstate = PGERR_GEOS_EXCEPTION;
2867 goto finish;
2868 }
2869 Py_XDECREF(*out);
2870 if (hex) {
2871 *out = PyUnicode_FromStringAndSize((char*)wkb, size);
2872 } else {
2873 *out = PyBytes_FromStringAndSize((char*)wkb, size);
2874 }
2875 GEOSFree_r(ctx, wkb);
2876 }
2877 }
2878
2879 finish:
2880 GEOSWKBWriter_destroy_r(ctx, writer);
2881 GEOS_FINISH;
2882 }
2883 static PyUFuncGenericFunction to_wkb_funcs[1] = {&to_wkb_func};
2884
2885 static char to_wkt_dtypes[6] = {NPY_OBJECT, NPY_INT, NPY_BOOL,
2886 NPY_INT, NPY_BOOL, NPY_OBJECT};
to_wkt_func(char ** args,npy_intp * dimensions,npy_intp * steps,void * data)2887 static void to_wkt_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
2888 char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4],
2889 *op1 = args[5];
2890 npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4],
2891 os1 = steps[5];
2892 npy_intp n = dimensions[0];
2893 npy_intp i;
2894
2895 GEOSGeometry* in1;
2896 GEOSWKTWriter* writer;
2897 char* wkt;
2898
2899 if ((is2 != 0) || (is3 != 0) || (is4 != 0) || (is5 != 0)) {
2900 PyErr_Format(PyExc_ValueError, "to_wkt function called with non-scalar parameters");
2901 return;
2902 }
2903
2904 GEOS_INIT;
2905
2906 /* Create the WKT writer */
2907 writer = GEOSWKTWriter_create_r(ctx);
2908 if (writer == NULL) {
2909 errstate = PGERR_GEOS_EXCEPTION;
2910 goto finish;
2911 }
2912 GEOSWKTWriter_setRoundingPrecision_r(ctx, writer, *(int*)ip2);
2913 GEOSWKTWriter_setTrim_r(ctx, writer, *(npy_bool*)ip3);
2914 GEOSWKTWriter_setOutputDimension_r(ctx, writer, *(int*)ip4);
2915 GEOSWKTWriter_setOld3D_r(ctx, writer, *(npy_bool*)ip5);
2916
2917 // Check if the above functions caused a GEOS exception
2918 if (last_error[0] != 0) {
2919 errstate = PGERR_GEOS_EXCEPTION;
2920 goto finish;
2921 }
2922
2923 for (i = 0; i < n; i++, ip1 += is1, op1 += os1) {
2924 if (!get_geom(*(GeometryObject**)ip1, &in1)) {
2925 errstate = PGERR_NOT_A_GEOMETRY;
2926 goto finish;
2927 }
2928 PyObject** out = (PyObject**)op1;
2929
2930 if (in1 == NULL) {
2931 Py_XDECREF(*out);
2932 Py_INCREF(Py_None);
2933 *out = Py_None;
2934 } else {
2935 errstate = check_to_wkt_compatible(ctx, in1);
2936 if (errstate != PGERR_SUCCESS) {
2937 goto finish;
2938 }
2939 wkt = GEOSWKTWriter_write_r(ctx, writer, in1);
2940 if (wkt == NULL) {
2941 errstate = PGERR_GEOS_EXCEPTION;
2942 goto finish;
2943 }
2944 Py_XDECREF(*out);
2945 *out = PyUnicode_FromString(wkt);
2946 GEOSFree_r(ctx, wkt);
2947 }
2948 }
2949
2950 finish:
2951 GEOSWKTWriter_destroy_r(ctx, writer);
2952 GEOS_FINISH;
2953 }
2954 static PyUFuncGenericFunction to_wkt_funcs[1] = {&to_wkt_func};
2955
2956 /*
2957 TODO polygonizer functions
2958 TODO prepared geometry predicate functions
2959 TODO relate functions
2960 */
2961
2962 #define DEFINE_Y_b(NAME) \
2963 ufunc = PyUFunc_FromFuncAndData(Y_b_funcs, NAME##_data, Y_b_dtypes, 1, 1, 1, \
2964 PyUFunc_None, #NAME, NULL, 0); \
2965 PyDict_SetItemString(d, #NAME, ufunc)
2966
2967 #define DEFINE_O_b(NAME) \
2968 ufunc = PyUFunc_FromFuncAndData(O_b_funcs, NAME##_data, O_b_dtypes, 1, 1, 1, \
2969 PyUFunc_None, #NAME, NULL, 0); \
2970 PyDict_SetItemString(d, #NAME, ufunc)
2971
2972 #define DEFINE_YY_b(NAME) \
2973 ufunc = PyUFunc_FromFuncAndData(YY_b_funcs, NAME##_data, YY_b_dtypes, 1, 2, 1, \
2974 PyUFunc_None, #NAME, "", 0); \
2975 PyDict_SetItemString(d, #NAME, ufunc)
2976
2977 #define DEFINE_YY_b_p(NAME) \
2978 ufunc = PyUFunc_FromFuncAndData(YY_b_p_funcs, NAME##_data, YY_b_p_dtypes, 1, 2, 1, \
2979 PyUFunc_None, #NAME, "", 0); \
2980 PyDict_SetItemString(d, #NAME, ufunc)
2981
2982 #define DEFINE_Y_Y(NAME) \
2983 ufunc = PyUFunc_FromFuncAndData(Y_Y_funcs, NAME##_data, Y_Y_dtypes, 1, 1, 1, \
2984 PyUFunc_None, #NAME, "", 0); \
2985 PyDict_SetItemString(d, #NAME, ufunc)
2986
2987 #define DEFINE_Y(NAME) \
2988 ufunc = PyUFunc_FromFuncAndData(Y_funcs, NAME##_data, Y_dtypes, 1, 1, 0, PyUFunc_None, \
2989 #NAME, "", 0); \
2990 PyDict_SetItemString(d, #NAME, ufunc)
2991
2992 #define DEFINE_Yi_Y(NAME) \
2993 ufunc = PyUFunc_FromFuncAndData(Yi_Y_funcs, NAME##_data, Yi_Y_dtypes, 1, 2, 1, \
2994 PyUFunc_None, #NAME, "", 0); \
2995 PyDict_SetItemString(d, #NAME, ufunc)
2996
2997 #define DEFINE_Yd_Y(NAME) \
2998 ufunc = PyUFunc_FromFuncAndData(Yd_Y_funcs, NAME##_data, Yd_Y_dtypes, 1, 2, 1, \
2999 PyUFunc_None, #NAME, "", 0); \
3000 PyDict_SetItemString(d, #NAME, ufunc)
3001
3002 #define DEFINE_YY_Y(NAME) \
3003 ufunc = PyUFunc_FromFuncAndData(YY_Y_funcs, NAME##_data, YY_Y_dtypes, 1, 2, 1, \
3004 PyUFunc_None, #NAME, "", 0); \
3005 PyDict_SetItemString(d, #NAME, ufunc)
3006
3007 #define DEFINE_YY_Y_REORDERABLE(NAME) \
3008 ufunc = PyUFunc_FromFuncAndData(YY_Y_funcs, NAME##_data, YY_Y_dtypes, 1, 2, 1, \
3009 PyUFunc_ReorderableNone, #NAME, "", 0); \
3010 PyDict_SetItemString(d, #NAME, ufunc)
3011
3012 #define DEFINE_Y_d(NAME) \
3013 ufunc = PyUFunc_FromFuncAndData(Y_d_funcs, NAME##_data, Y_d_dtypes, 1, 1, 1, \
3014 PyUFunc_None, #NAME, "", 0); \
3015 PyDict_SetItemString(d, #NAME, ufunc)
3016
3017 #define DEFINE_Y_B(NAME) \
3018 ufunc = PyUFunc_FromFuncAndData(Y_B_funcs, NAME##_data, Y_B_dtypes, 1, 1, 1, \
3019 PyUFunc_None, #NAME, "", 0); \
3020 PyDict_SetItemString(d, #NAME, ufunc)
3021
3022 #define DEFINE_Y_i(NAME) \
3023 ufunc = PyUFunc_FromFuncAndData(Y_i_funcs, NAME##_data, Y_i_dtypes, 1, 1, 1, \
3024 PyUFunc_None, #NAME, "", 0); \
3025 PyDict_SetItemString(d, #NAME, ufunc)
3026
3027 #define DEFINE_YY_d(NAME) \
3028 ufunc = PyUFunc_FromFuncAndData(YY_d_funcs, NAME##_data, YY_d_dtypes, 1, 2, 1, \
3029 PyUFunc_None, #NAME, "", 0); \
3030 PyDict_SetItemString(d, #NAME, ufunc)
3031
3032 #define DEFINE_YYd_d(NAME) \
3033 ufunc = PyUFunc_FromFuncAndData(YYd_d_funcs, NAME##_data, YYd_d_dtypes, 1, 3, 1, \
3034 PyUFunc_None, #NAME, "", 0); \
3035 PyDict_SetItemString(d, #NAME, ufunc)
3036
3037 #define DEFINE_YYd_Y(NAME) \
3038 ufunc = PyUFunc_FromFuncAndData(YYd_Y_funcs, NAME##_data, YYd_Y_dtypes, 1, 3, 1, \
3039 PyUFunc_None, #NAME, "", 0); \
3040 PyDict_SetItemString(d, #NAME, ufunc)
3041
3042 #define DEFINE_CUSTOM(NAME, N_IN) \
3043 ufunc = PyUFunc_FromFuncAndData(NAME##_funcs, null_data, NAME##_dtypes, 1, N_IN, 1, \
3044 PyUFunc_None, #NAME, "", 0); \
3045 PyDict_SetItemString(d, #NAME, ufunc)
3046
3047 #define DEFINE_GENERALIZED(NAME, N_IN, SIGNATURE) \
3048 ufunc = PyUFunc_FromFuncAndDataAndSignature(NAME##_funcs, null_data, NAME##_dtypes, 1, \
3049 N_IN, 1, PyUFunc_None, #NAME, "", 0, \
3050 SIGNATURE); \
3051 PyDict_SetItemString(d, #NAME, ufunc)
3052
3053 #define DEFINE_GENERALIZED_NOUT4(NAME, N_IN, SIGNATURE) \
3054 ufunc = PyUFunc_FromFuncAndDataAndSignature(NAME##_funcs, null_data, NAME##_dtypes, 1, \
3055 N_IN, 4, PyUFunc_None, #NAME, "", 0, \
3056 SIGNATURE); \
3057 PyDict_SetItemString(d, #NAME, ufunc)
3058
init_ufuncs(PyObject * m,PyObject * d)3059 int init_ufuncs(PyObject* m, PyObject* d) {
3060 PyObject* ufunc;
3061
3062 DEFINE_Y_b(is_empty);
3063 DEFINE_Y_b(is_simple);
3064 DEFINE_Y_b(is_geometry);
3065 DEFINE_Y_b(is_ring);
3066 DEFINE_Y_b(has_z);
3067 DEFINE_Y_b(is_closed);
3068 DEFINE_Y_b(is_valid);
3069
3070 DEFINE_O_b(is_geometry);
3071 DEFINE_O_b(is_missing);
3072 DEFINE_O_b(is_valid_input);
3073
3074 DEFINE_YY_b_p(disjoint);
3075 DEFINE_YY_b_p(touches);
3076 DEFINE_YY_b_p(intersects);
3077 DEFINE_YY_b_p(crosses);
3078 DEFINE_YY_b_p(within);
3079 DEFINE_YY_b_p(contains);
3080 DEFINE_YY_b_p(contains_properly);
3081 DEFINE_YY_b_p(overlaps);
3082 DEFINE_YY_b(equals);
3083 DEFINE_YY_b_p(covers);
3084 DEFINE_YY_b_p(covered_by);
3085 DEFINE_CUSTOM(is_prepared, 1);
3086
3087 DEFINE_Y_Y(envelope);
3088 DEFINE_Y_Y(convex_hull);
3089 DEFINE_Y_Y(boundary);
3090 DEFINE_Y_Y(unary_union);
3091 DEFINE_Y_Y(point_on_surface);
3092 DEFINE_Y_Y(centroid);
3093 DEFINE_Y_Y(line_merge);
3094 DEFINE_Y_Y(extract_unique_points);
3095 DEFINE_Y_Y(get_exterior_ring);
3096 DEFINE_Y_Y(normalize);
3097
3098 DEFINE_Y(prepare);
3099 DEFINE_Y(destroy_prepared);
3100
3101 DEFINE_Yi_Y(get_point);
3102 DEFINE_Yi_Y(get_interior_ring);
3103 DEFINE_Yi_Y(get_geometry);
3104 DEFINE_Yi_Y(set_srid);
3105
3106 DEFINE_Yd_Y(line_interpolate_point);
3107 DEFINE_Yd_Y(line_interpolate_point_normalized);
3108 DEFINE_Yd_Y(simplify);
3109 DEFINE_Yd_Y(simplify_preserve_topology);
3110
3111 // 'REORDERABLE' sets PyUFunc_ReorderableNone instead of PyUFunc_None as 'identity',
3112 // meaning that the order of arguments does not matter. This enables reduction over
3113 // multiple (or all) axes.
3114 DEFINE_YY_Y_REORDERABLE(intersection);
3115 DEFINE_YY_Y(difference);
3116 DEFINE_YY_Y_REORDERABLE(symmetric_difference);
3117 DEFINE_YY_Y_REORDERABLE(union);
3118 DEFINE_YY_Y(shared_paths);
3119
3120 DEFINE_Y_d(get_x);
3121 DEFINE_Y_d(get_y);
3122 DEFINE_Y_d(area);
3123 DEFINE_Y_d(length);
3124
3125 DEFINE_Y_i(get_type_id);
3126 DEFINE_Y_i(get_dimensions);
3127 DEFINE_Y_i(get_coordinate_dimension);
3128 DEFINE_Y_i(get_srid);
3129 DEFINE_Y_i(get_num_points);
3130 DEFINE_Y_i(get_num_interior_rings);
3131 DEFINE_Y_i(get_num_geometries);
3132 DEFINE_Y_i(get_num_coordinates);
3133
3134 DEFINE_YY_d(distance);
3135 DEFINE_YY_d(hausdorff_distance);
3136 DEFINE_YY_d(line_locate_point);
3137 DEFINE_YY_d(line_locate_point_normalized);
3138
3139 DEFINE_YYd_d(hausdorff_distance_densify);
3140
3141 DEFINE_CUSTOM(box, 5);
3142 DEFINE_CUSTOM(buffer, 7);
3143 DEFINE_CUSTOM(offset_curve, 5);
3144 DEFINE_CUSTOM(snap, 3);
3145 DEFINE_CUSTOM(clip_by_rect, 5);
3146 DEFINE_CUSTOM(equals_exact, 3);
3147
3148 DEFINE_CUSTOM(delaunay_triangles, 3);
3149 DEFINE_CUSTOM(voronoi_polygons, 4);
3150 DEFINE_CUSTOM(is_valid_reason, 1);
3151 DEFINE_CUSTOM(relate, 2);
3152 DEFINE_CUSTOM(relate_pattern, 3);
3153 DEFINE_GENERALIZED(polygonize, 1, "(d)->()");
3154 DEFINE_GENERALIZED_NOUT4(polygonize_full, 1, "(d)->(),(),(),()");
3155 DEFINE_CUSTOM(shortest_line, 2);
3156
3157 DEFINE_GENERALIZED(points, 1, "(d)->()");
3158 DEFINE_GENERALIZED(linestrings, 1, "(i, d)->()");
3159 DEFINE_GENERALIZED(linearrings, 1, "(i, d)->()");
3160 DEFINE_GENERALIZED(bounds, 1, "()->(n)");
3161 DEFINE_GENERALIZED(polygons, 2, "(),(i)->()");
3162 DEFINE_GENERALIZED(create_collection, 2, "(i),()->()");
3163
3164 DEFINE_CUSTOM(from_wkb, 2);
3165 DEFINE_CUSTOM(from_wkt, 2);
3166 DEFINE_CUSTOM(to_wkb, 5);
3167 DEFINE_CUSTOM(to_wkt, 5);
3168 DEFINE_CUSTOM(from_shapely, 1);
3169
3170 #if GEOS_SINCE_3_6_0
3171 DEFINE_Y_d(minimum_clearance);
3172 DEFINE_Y_d(get_precision);
3173 DEFINE_Y_Y(oriented_envelope);
3174 DEFINE_CUSTOM(set_precision, 3);
3175 #endif
3176
3177 #if GEOS_SINCE_3_7_0
3178 DEFINE_Y_b(is_ccw);
3179 DEFINE_Y_d(get_z);
3180 DEFINE_Y_Y(reverse);
3181 DEFINE_YY_d(frechet_distance);
3182 DEFINE_YYd_d(frechet_distance_densify);
3183 #endif
3184
3185 #if GEOS_SINCE_3_8_0
3186 DEFINE_Y_Y(make_valid);
3187 DEFINE_Y_Y(build_area);
3188 DEFINE_Y_Y(coverage_union);
3189 DEFINE_Y_Y(minimum_bounding_circle);
3190 DEFINE_Y_d(minimum_bounding_radius);
3191 #endif
3192
3193 #if GEOS_SINCE_3_9_0
3194 DEFINE_YYd_Y(difference_prec);
3195 DEFINE_YYd_Y(intersection_prec);
3196 DEFINE_YYd_Y(symmetric_difference_prec);
3197 DEFINE_YYd_Y(union_prec);
3198 DEFINE_Yd_Y(unary_union_prec);
3199 #endif
3200
3201 #if GEOS_SINCE_3_10_0
3202 DEFINE_Yd_Y(segmentize);
3203 #endif
3204
3205 Py_DECREF(ufunc);
3206 return 0;
3207 }
3208