1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2009-2014 Sandro Santilli <strk@kbt.io>
22 * Copyright 2008 Paul Ramsey <pramsey@cleverelephant.ca>
23 * Copyright 2001-2003 Refractions Research Inc.
24 *
25 **********************************************************************/
26
27
28 #include "../postgis_config.h"
29
30 #include "float.h" /* for DBL_DIG */
31
32 /* PostgreSQL */
33 #include "postgres.h"
34 #include "funcapi.h"
35 #include "utils/array.h"
36 #include "utils/builtins.h"
37 #include "utils/lsyscache.h"
38 #include "utils/numeric.h"
39 #include "access/htup_details.h"
40
41 /* PostGIS */
42 #include "lwgeom_functions_analytic.h" /* for point_in_polygon */
43 #include "lwgeom_geos.h"
44 #include "liblwgeom.h"
45 #include "lwgeom_rtree.h"
46 #include "lwgeom_geos_prepared.h"
47 #include "lwgeom_accum.h"
48
49
50 /* Return NULL on GEOS error
51 *
52 * Prints error message only if it was not for interruption, in which
53 * case we let PostgreSQL deal with the error.
54 */
55 #define HANDLE_GEOS_ERROR(label) \
56 { \
57 if (strstr(lwgeom_geos_errmsg, "InterruptedException")) \
58 ereport(ERROR, \
59 (errcode(ERRCODE_QUERY_CANCELED), errmsg("canceling statement due to user request"))); \
60 else \
61 lwpgerror("%s: %s", (label), lwgeom_geos_errmsg); \
62 PG_RETURN_NULL(); \
63 }
64
65 /*
66 ** Prototypes for SQL-bound functions
67 */
68 Datum relate_full(PG_FUNCTION_ARGS);
69 Datum relate_pattern(PG_FUNCTION_ARGS);
70 Datum disjoint(PG_FUNCTION_ARGS);
71 Datum touches(PG_FUNCTION_ARGS);
72 Datum ST_Intersects(PG_FUNCTION_ARGS);
73 Datum crosses(PG_FUNCTION_ARGS);
74 Datum contains(PG_FUNCTION_ARGS);
75 Datum within(PG_FUNCTION_ARGS);
76 Datum containsproperly(PG_FUNCTION_ARGS);
77 Datum covers(PG_FUNCTION_ARGS);
78 Datum overlaps(PG_FUNCTION_ARGS);
79 Datum isvalid(PG_FUNCTION_ARGS);
80 Datum isvalidreason(PG_FUNCTION_ARGS);
81 Datum isvaliddetail(PG_FUNCTION_ARGS);
82 Datum buffer(PG_FUNCTION_ARGS);
83 Datum ST_Intersection(PG_FUNCTION_ARGS);
84 Datum convexhull(PG_FUNCTION_ARGS);
85 Datum topologypreservesimplify(PG_FUNCTION_ARGS);
86 Datum ST_Difference(PG_FUNCTION_ARGS);
87 Datum boundary(PG_FUNCTION_ARGS);
88 Datum symdifference(PG_FUNCTION_ARGS);
89 Datum ST_Union(PG_FUNCTION_ARGS);
90 Datum issimple(PG_FUNCTION_ARGS);
91 Datum isring(PG_FUNCTION_ARGS);
92 Datum pointonsurface(PG_FUNCTION_ARGS);
93 Datum GEOSnoop(PG_FUNCTION_ARGS);
94 Datum postgis_geos_version(PG_FUNCTION_ARGS);
95 Datum centroid(PG_FUNCTION_ARGS);
96 Datum polygonize_garray(PG_FUNCTION_ARGS);
97 Datum clusterintersecting_garray(PG_FUNCTION_ARGS);
98 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS);
99 Datum linemerge(PG_FUNCTION_ARGS);
100 Datum coveredby(PG_FUNCTION_ARGS);
101 Datum hausdorffdistance(PG_FUNCTION_ARGS);
102 Datum hausdorffdistancedensify(PG_FUNCTION_ARGS);
103 Datum ST_FrechetDistance(PG_FUNCTION_ARGS);
104 Datum ST_UnaryUnion(PG_FUNCTION_ARGS);
105 Datum ST_Equals(PG_FUNCTION_ARGS);
106 Datum ST_BuildArea(PG_FUNCTION_ARGS);
107 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS);
108
109 Datum pgis_union_geometry_array(PG_FUNCTION_ARGS);
110 Datum pgis_geometry_union_finalfn(PG_FUNCTION_ARGS);
111
112 /*
113 ** Prototypes end
114 */
115
116
117 PG_FUNCTION_INFO_V1(postgis_geos_version);
postgis_geos_version(PG_FUNCTION_ARGS)118 Datum postgis_geos_version(PG_FUNCTION_ARGS)
119 {
120 const char *ver = lwgeom_geos_version();
121 text *result = cstring_to_text(ver);
122 PG_RETURN_POINTER(result);
123 }
124
125 static char
is_poly(const GSERIALIZED * g)126 is_poly(const GSERIALIZED* g)
127 {
128 int type = gserialized_get_type(g);
129 return type == POLYGONTYPE || type == MULTIPOLYGONTYPE;
130 }
131
132 static char
is_point(const GSERIALIZED * g)133 is_point(const GSERIALIZED* g)
134 {
135 int type = gserialized_get_type(g);
136 return type == POINTTYPE || type == MULTIPOINTTYPE;
137 }
138
139 /* Utility function that checks a LWPOINT and a GSERIALIZED poly against a cache.
140 * Serialized poly may be a multipart.
141 */
142 static int
pip_short_circuit(RTREE_POLY_CACHE * poly_cache,LWPOINT * point,GSERIALIZED * gpoly)143 pip_short_circuit(RTREE_POLY_CACHE* poly_cache, LWPOINT* point, GSERIALIZED* gpoly)
144 {
145 int result;
146
147 if ( poly_cache && poly_cache->ringIndices )
148 {
149 result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
150 }
151 else
152 {
153 LWGEOM* poly = lwgeom_from_gserialized(gpoly);
154 if ( lwgeom_get_type(poly) == POLYGONTYPE )
155 {
156 result = point_in_polygon(lwgeom_as_lwpoly(poly), point);
157 }
158 else
159 {
160 result = point_in_multipolygon(lwgeom_as_lwmpoly(poly), point);
161 }
162 lwgeom_free(poly);
163 }
164
165 return result;
166 }
167
168 /**
169 * @brief Compute the Hausdorff distance thanks to the corresponding GEOS function
170 * @example ST_HausdorffDistance {@link #hausdorffdistance} - SELECT ST_HausdorffDistance(
171 * 'POLYGON((0 0, 0 2, 1 2, 2 2, 2 0, 0 0))'::geometry,
172 * 'POLYGON((0.5 0.5, 0.5 2.5, 1.5 2.5, 2.5 2.5, 2.5 0.5, 0.5 0.5))'::geometry);
173 */
174
175 PG_FUNCTION_INFO_V1(hausdorffdistance);
hausdorffdistance(PG_FUNCTION_ARGS)176 Datum hausdorffdistance(PG_FUNCTION_ARGS)
177 {
178 GSERIALIZED *geom1;
179 GSERIALIZED *geom2;
180 GEOSGeometry *g1;
181 GEOSGeometry *g2;
182 double result;
183 int retcode;
184
185 geom1 = PG_GETARG_GSERIALIZED_P(0);
186 geom2 = PG_GETARG_GSERIALIZED_P(1);
187
188 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
189 PG_RETURN_NULL();
190
191 initGEOS(lwpgnotice, lwgeom_geos_error);
192
193 g1 = POSTGIS2GEOS(geom1);
194 if (!g1)
195 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
196
197 g2 = POSTGIS2GEOS(geom2);
198 if (!g2)
199 {
200 GEOSGeom_destroy(g1);
201 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
202 }
203
204 retcode = GEOSHausdorffDistance(g1, g2, &result);
205 GEOSGeom_destroy(g1);
206 GEOSGeom_destroy(g2);
207
208 if (retcode == 0) HANDLE_GEOS_ERROR("GEOSHausdorffDistance");
209
210 PG_FREE_IF_COPY(geom1, 0);
211 PG_FREE_IF_COPY(geom2, 1);
212
213 PG_RETURN_FLOAT8(result);
214 }
215
216 /**
217 * @brief Compute the Hausdorff distance with densification thanks to the corresponding GEOS function
218 * @example hausdorffdistancedensify {@link #hausdorffdistancedensify} - SELECT ST_HausdorffDistance(
219 * 'POLYGON((0 0, 0 2, 1 2, 2 2, 2 0, 0 0))'::geometry,
220 * 'POLYGON((0.5 0.5, 0.5 2.5, 1.5 2.5, 2.5 2.5, 2.5 0.5, 0.5 0.5))'::geometry,
221 * 0.5);
222 */
223
224 PG_FUNCTION_INFO_V1(hausdorffdistancedensify);
hausdorffdistancedensify(PG_FUNCTION_ARGS)225 Datum hausdorffdistancedensify(PG_FUNCTION_ARGS)
226 {
227 GSERIALIZED *geom1;
228 GSERIALIZED *geom2;
229 GEOSGeometry *g1;
230 GEOSGeometry *g2;
231 double densifyFrac;
232 double result;
233 int retcode;
234
235 geom1 = PG_GETARG_GSERIALIZED_P(0);
236 geom2 = PG_GETARG_GSERIALIZED_P(1);
237 densifyFrac = PG_GETARG_FLOAT8(2);
238
239 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
240 PG_RETURN_NULL();
241
242 initGEOS(lwpgnotice, lwgeom_geos_error);
243
244 g1 = POSTGIS2GEOS(geom1);
245 if (!g1)
246 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
247
248 g2 = POSTGIS2GEOS(geom2);
249 if (!g2)
250 {
251 GEOSGeom_destroy(g1);
252 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
253 }
254
255 retcode = GEOSHausdorffDistanceDensify(g1, g2, densifyFrac, &result);
256 GEOSGeom_destroy(g1);
257 GEOSGeom_destroy(g2);
258
259 if (retcode == 0) HANDLE_GEOS_ERROR("GEOSHausdorffDistanceDensify");
260
261 PG_FREE_IF_COPY(geom1, 0);
262 PG_FREE_IF_COPY(geom2, 1);
263
264 PG_RETURN_FLOAT8(result);
265 }
266
267 /**
268 * @brief Compute the Frechet distance with optional densification thanks to the corresponding GEOS function
269 * @example ST_FrechetDistance {@link #frechetdistance} - SELECT ST_FrechetDistance(
270 * 'LINESTRING (0 0, 50 200, 100 0, 150 200, 200 0)'::geometry,
271 * 'LINESTRING (0 200, 200 150, 0 100, 200 50, 0 0)'::geometry, 0.5);
272 */
273
274 PG_FUNCTION_INFO_V1(ST_FrechetDistance);
ST_FrechetDistance(PG_FUNCTION_ARGS)275 Datum ST_FrechetDistance(PG_FUNCTION_ARGS)
276 {
277 #if POSTGIS_GEOS_VERSION < 37
278
279 lwpgerror("The GEOS version this PostGIS binary "
280 "was compiled against (%d) doesn't support "
281 "'GEOSFechetDistance' function (3.7.0+ required)",
282 POSTGIS_GEOS_VERSION);
283 PG_RETURN_NULL();
284
285 #else /* POSTGIS_GEOS_VERSION >= 37 */
286 GSERIALIZED *geom1;
287 GSERIALIZED *geom2;
288 GEOSGeometry *g1;
289 GEOSGeometry *g2;
290 double densifyFrac;
291 double result;
292 int retcode;
293
294 geom1 = PG_GETARG_GSERIALIZED_P(0);
295 geom2 = PG_GETARG_GSERIALIZED_P(1);
296 densifyFrac = PG_GETARG_FLOAT8(2);
297
298 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
299 PG_RETURN_NULL();
300
301 initGEOS(lwpgnotice, lwgeom_geos_error);
302
303 g1 = POSTGIS2GEOS(geom1);
304 if (!g1)
305 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
306
307 g2 = POSTGIS2GEOS(geom2);
308 if (!g2)
309 {
310 GEOSGeom_destroy(g1);
311 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
312 }
313
314 if (densifyFrac <= 0.0)
315 {
316 retcode = GEOSFrechetDistance(g1, g2, &result);
317 }
318 else
319 {
320 retcode = GEOSFrechetDistanceDensify(g1, g2, densifyFrac, &result);
321 }
322
323 GEOSGeom_destroy(g1);
324 GEOSGeom_destroy(g2);
325
326 if (retcode == 0) HANDLE_GEOS_ERROR("GEOSFrechetDistance");
327
328 PG_FREE_IF_COPY(geom1, 0);
329 PG_FREE_IF_COPY(geom2, 1);
330
331 PG_RETURN_FLOAT8(result);
332
333 #endif /* POSTGIS_GEOS_VERSION >= 37 */
334 }
335
336 /**
337 * @brief This is the final function for GeomUnion
338 * aggregate. Will have as input an array of Geometries.
339 * Will iteratively call GEOSUnion on the GEOS-converted
340 * versions of them and return PGIS-converted version back.
341 * Changing combination order *might* speed up performance.
342 */
343 PG_FUNCTION_INFO_V1(pgis_union_geometry_array);
pgis_union_geometry_array(PG_FUNCTION_ARGS)344 Datum pgis_union_geometry_array(PG_FUNCTION_ARGS)
345 {
346 ArrayType *array;
347
348 ArrayIterator iterator;
349 Datum value;
350 bool isnull;
351
352 int is3d = LW_FALSE, gotsrid = LW_FALSE;
353 int nelems = 0, geoms_size = 0, curgeom = 0, count = 0;
354
355 GSERIALIZED *gser_out = NULL;
356
357 GEOSGeometry *g = NULL;
358 GEOSGeometry *g_union = NULL;
359 GEOSGeometry **geoms = NULL;
360
361 int32_t srid = SRID_UNKNOWN;
362
363 int empty_type = 0;
364
365 /* Null array, null geometry (should be empty?) */
366 if ( PG_ARGISNULL(0) )
367 PG_RETURN_NULL();
368
369 array = PG_GETARG_ARRAYTYPE_P(0);
370 nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
371
372 /* Empty array? Null return */
373 if ( nelems == 0 ) PG_RETURN_NULL();
374
375 /* Quick scan for nulls */
376 iterator = array_create_iterator(array, 0, NULL);
377 while (array_iterate(iterator, &value, &isnull))
378 {
379 /* Skip null array items */
380 if (isnull) continue;
381 count++;
382 }
383 array_free_iterator(iterator);
384
385
386 /* All-nulls? Return null */
387 if ( count == 0 )
388 PG_RETURN_NULL();
389
390 /* One geom, good geom? Return it */
391 if ( count == 1 && nelems == 1 )
392 {
393 #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)
394 #endif
395 PG_RETURN_POINTER((GSERIALIZED *)(ARR_DATA_PTR(array)));
396 #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)
397 #endif
398 }
399
400 /* Ok, we really need GEOS now ;) */
401 initGEOS(lwpgnotice, lwgeom_geos_error);
402
403 /*
404 ** Collect the non-empty inputs and stuff them into a GEOS collection
405 */
406 geoms_size = nelems;
407 geoms = palloc(sizeof(GEOSGeometry*) * geoms_size);
408
409 /*
410 ** We need to convert the array of GSERIALIZED into a GEOS collection.
411 ** First make an array of GEOS geometries.
412 */
413 iterator = array_create_iterator(array, 0, NULL);
414 while (array_iterate(iterator, &value, &isnull))
415 {
416 GSERIALIZED *gser_in;
417
418 /* Skip null array items */
419 if (isnull) continue;
420 gser_in = (GSERIALIZED *)DatumGetPointer(value);
421
422 /* Check for SRID mismatch in array elements */
423 if ( gotsrid )
424 gserialized_error_if_srid_mismatch_reference(gser_in, srid, __func__);
425 else
426 {
427 /* Initialize SRID/dimensions info */
428 srid = gserialized_get_srid(gser_in);
429 is3d = gserialized_has_z(gser_in);
430 gotsrid = 1;
431 }
432
433 /* Don't include empties in the union */
434 if ( gserialized_is_empty(gser_in) )
435 {
436 int gser_type = gserialized_get_type(gser_in);
437 if (gser_type > empty_type)
438 {
439 empty_type = gser_type;
440 POSTGIS_DEBUGF(4, "empty_type = %d gser_type = %d", empty_type, gser_type);
441 }
442 }
443 else
444 {
445 g = POSTGIS2GEOS(gser_in);
446
447 /* Uh oh! Exception thrown at construction... */
448 if ( ! g )
449 {
450 HANDLE_GEOS_ERROR(
451 "One of the geometries in the set "
452 "could not be converted to GEOS");
453 }
454
455 /* Ensure we have enough space in our storage array */
456 if ( curgeom == geoms_size )
457 {
458 geoms_size *= 2;
459 geoms = repalloc( geoms, sizeof(GEOSGeometry*) * geoms_size );
460 }
461
462 geoms[curgeom] = g;
463 curgeom++;
464 }
465
466 }
467 array_free_iterator(iterator);
468
469 /*
470 ** Take our GEOS geometries and turn them into a GEOS collection,
471 ** then pass that into cascaded union.
472 */
473 if (curgeom > 0)
474 {
475 g = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, curgeom);
476 if (!g) HANDLE_GEOS_ERROR("Could not create GEOS COLLECTION from geometry array");
477
478 g_union = GEOSUnaryUnion(g);
479 GEOSGeom_destroy(g);
480 if (!g_union) HANDLE_GEOS_ERROR("GEOSUnaryUnion");
481
482 GEOSSetSRID(g_union, srid);
483 gser_out = GEOS2POSTGIS(g_union, is3d);
484 GEOSGeom_destroy(g_union);
485 }
486 /* No real geometries in our array, any empties? */
487 else
488 {
489 /* If it was only empties, we'll return the largest type number */
490 if ( empty_type > 0 )
491 {
492 PG_RETURN_POINTER(geometry_serialize(lwgeom_construct_empty(empty_type, srid, is3d, 0)));
493 }
494 /* Nothing but NULL, returns NULL */
495 else
496 {
497 PG_RETURN_NULL();
498 }
499 }
500
501 if ( ! gser_out )
502 {
503 /* Union returned a NULL geometry */
504 PG_RETURN_NULL();
505 }
506
507 PG_RETURN_POINTER(gser_out);
508 }
509
510
511 PG_FUNCTION_INFO_V1(pgis_geometry_union_finalfn);
pgis_geometry_union_finalfn(PG_FUNCTION_ARGS)512 Datum pgis_geometry_union_finalfn(PG_FUNCTION_ARGS)
513 {
514 CollectionBuildState *state;
515 ListCell *l;
516 LWGEOM **geoms;
517 GSERIALIZED *gser_out;
518 size_t ngeoms = 0;
519 int empty_type = 0;
520 bool first = true;
521 int32_t srid = SRID_UNKNOWN;
522 int has_z = LW_FALSE;
523
524 if (PG_ARGISNULL(0))
525 PG_RETURN_NULL(); /* returns null iff no input values */
526
527 state = (CollectionBuildState *)PG_GETARG_POINTER(0);
528 geoms = palloc(list_length(state->geoms) * sizeof(LWGEOM*));
529
530 /* Read contents of list into an array of only non-null values */
531 foreach (l, state->geoms)
532 {
533 LWGEOM *geom = (LWGEOM*)(lfirst(l));
534 if (geom)
535 {
536 if (!lwgeom_is_empty(geom))
537 {
538 geoms[ngeoms++] = geom;
539 if (first)
540 {
541 srid = lwgeom_get_srid(geom);
542 has_z = lwgeom_has_z(geom);
543 first = false;
544 }
545 }
546 else
547 {
548 int type = lwgeom_get_type(geom);
549 empty_type = type > empty_type ? type : empty_type;
550 srid = (srid != SRID_UNKNOWN ? srid : lwgeom_get_srid(geom));
551 }
552 }
553 }
554
555 /*
556 ** Take our array of LWGEOM* and turn it into a GEOS collection,
557 ** then pass that into cascaded union.
558 */
559 if (ngeoms > 0)
560 {
561 GEOSGeometry *g = NULL;
562 GEOSGeometry *g_union = NULL;
563 LWCOLLECTION* col = lwcollection_construct(COLLECTIONTYPE, srid, NULL, ngeoms, geoms);
564
565 initGEOS(lwpgnotice, lwgeom_geos_error);
566 g = LWGEOM2GEOS((LWGEOM*)col, LW_FALSE);
567 if (!g)
568 HANDLE_GEOS_ERROR("Could not create GEOS COLLECTION from geometry array");
569
570 g_union = GEOSUnaryUnion(g);
571 GEOSGeom_destroy(g);
572 if (!g_union)
573 HANDLE_GEOS_ERROR("GEOSUnaryUnion");
574
575 GEOSSetSRID(g_union, srid);
576 gser_out = GEOS2POSTGIS(g_union, has_z);
577 GEOSGeom_destroy(g_union);
578 }
579 /* No real geometries in our array, any empties? */
580 else
581 {
582 /* If it was only empties, we'll return the largest type number */
583 if (empty_type > 0)
584 PG_RETURN_POINTER(
585 geometry_serialize(lwgeom_construct_empty(empty_type, srid, has_z, 0)));
586
587 /* Nothing but NULL, returns NULL */
588 else
589 PG_RETURN_NULL();
590 }
591
592 if (!gser_out)
593 {
594 /* Union returned a NULL geometry */
595 PG_RETURN_NULL();
596 }
597
598 PG_RETURN_POINTER(gser_out);
599 }
600
601
602
603 /**
604 * @example ST_UnaryUnion {@link #geomunion} SELECT ST_UnaryUnion(
605 * 'POLYGON((0 0, 10 0, 0 10, 10 10, 0 0))'
606 * );
607 *
608 */
609 PG_FUNCTION_INFO_V1(ST_UnaryUnion);
ST_UnaryUnion(PG_FUNCTION_ARGS)610 Datum ST_UnaryUnion(PG_FUNCTION_ARGS)
611 {
612 GSERIALIZED *geom1;
613 GSERIALIZED *result;
614 LWGEOM *lwgeom1, *lwresult ;
615
616 geom1 = PG_GETARG_GSERIALIZED_P(0);
617
618 lwgeom1 = lwgeom_from_gserialized(geom1) ;
619
620 lwresult = lwgeom_unaryunion(lwgeom1);
621 result = geometry_serialize(lwresult) ;
622
623 lwgeom_free(lwgeom1) ;
624 lwgeom_free(lwresult) ;
625
626 PG_FREE_IF_COPY(geom1, 0);
627
628 PG_RETURN_POINTER(result);
629 }
630
631 PG_FUNCTION_INFO_V1(ST_Union);
ST_Union(PG_FUNCTION_ARGS)632 Datum ST_Union(PG_FUNCTION_ARGS)
633 {
634 GSERIALIZED *geom1;
635 GSERIALIZED *geom2;
636 GSERIALIZED *result;
637 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
638
639 geom1 = PG_GETARG_GSERIALIZED_P(0);
640 geom2 = PG_GETARG_GSERIALIZED_P(1);
641
642 lwgeom1 = lwgeom_from_gserialized(geom1);
643 lwgeom2 = lwgeom_from_gserialized(geom2);
644
645 lwresult = lwgeom_union(lwgeom1, lwgeom2);
646 result = geometry_serialize(lwresult);
647
648 lwgeom_free(lwgeom1);
649 lwgeom_free(lwgeom2);
650 lwgeom_free(lwresult);
651
652 PG_FREE_IF_COPY(geom1, 0);
653 PG_FREE_IF_COPY(geom2, 1);
654
655 PG_RETURN_POINTER(result);
656 }
657
658 /**
659 * @example symdifference {@link #symdifference} - SELECT ST_SymDifference(
660 * 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
661 * 'POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
662 */
663 PG_FUNCTION_INFO_V1(symdifference);
symdifference(PG_FUNCTION_ARGS)664 Datum symdifference(PG_FUNCTION_ARGS)
665 {
666 GSERIALIZED *geom1;
667 GSERIALIZED *geom2;
668 GSERIALIZED *result;
669 LWGEOM *lwgeom1, *lwgeom2, *lwresult ;
670
671 geom1 = PG_GETARG_GSERIALIZED_P(0);
672 geom2 = PG_GETARG_GSERIALIZED_P(1);
673
674 lwgeom1 = lwgeom_from_gserialized(geom1) ;
675 lwgeom2 = lwgeom_from_gserialized(geom2) ;
676
677 lwresult = lwgeom_symdifference(lwgeom1, lwgeom2) ;
678 result = geometry_serialize(lwresult) ;
679
680 lwgeom_free(lwgeom1) ;
681 lwgeom_free(lwgeom2) ;
682 lwgeom_free(lwresult) ;
683
684 PG_FREE_IF_COPY(geom1, 0);
685 PG_FREE_IF_COPY(geom2, 1);
686
687 PG_RETURN_POINTER(result);
688 }
689
690
691 PG_FUNCTION_INFO_V1(boundary);
boundary(PG_FUNCTION_ARGS)692 Datum boundary(PG_FUNCTION_ARGS)
693 {
694 GSERIALIZED *geom1;
695 GEOSGeometry *g1, *g3;
696 GSERIALIZED *result;
697 LWGEOM *lwgeom;
698 int32_t srid;
699
700 geom1 = PG_GETARG_GSERIALIZED_P(0);
701
702 /* Empty.Boundary() == Empty */
703 if ( gserialized_is_empty(geom1) )
704 PG_RETURN_POINTER(geom1);
705
706 srid = gserialized_get_srid(geom1);
707
708 lwgeom = lwgeom_from_gserialized(geom1);
709 if ( ! lwgeom ) {
710 lwpgerror("POSTGIS2GEOS: unable to deserialize input");
711 PG_RETURN_NULL();
712 }
713
714 /* GEOS doesn't do triangle type, so we special case that here */
715 if (lwgeom->type == TRIANGLETYPE)
716 {
717 lwgeom->type = LINETYPE;
718 result = geometry_serialize(lwgeom);
719 lwgeom_free(lwgeom);
720 PG_RETURN_POINTER(result);
721 }
722
723 initGEOS(lwpgnotice, lwgeom_geos_error);
724
725 g1 = LWGEOM2GEOS(lwgeom, 0);
726 lwgeom_free(lwgeom);
727
728 if (!g1)
729 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
730
731 g3 = GEOSBoundary(g1);
732
733 if (!g3)
734 {
735 GEOSGeom_destroy(g1);
736 HANDLE_GEOS_ERROR("GEOSBoundary");
737 }
738
739 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
740
741 GEOSSetSRID(g3, srid);
742
743 result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
744
745 if (!result)
746 {
747 GEOSGeom_destroy(g1);
748 GEOSGeom_destroy(g3);
749 elog(NOTICE,
750 "GEOS2POSTGIS threw an error (result postgis geometry "
751 "formation)!");
752 PG_RETURN_NULL(); /* never get here */
753 }
754
755 GEOSGeom_destroy(g1);
756 GEOSGeom_destroy(g3);
757
758 PG_FREE_IF_COPY(geom1, 0);
759
760 PG_RETURN_POINTER(result);
761 }
762
763 PG_FUNCTION_INFO_V1(convexhull);
convexhull(PG_FUNCTION_ARGS)764 Datum convexhull(PG_FUNCTION_ARGS)
765 {
766 GSERIALIZED *geom1;
767 GEOSGeometry *g1, *g3;
768 GSERIALIZED *result;
769 LWGEOM *lwout;
770 int32_t srid;
771 GBOX bbox;
772
773 geom1 = PG_GETARG_GSERIALIZED_P(0);
774
775 /* Empty.ConvexHull() == Empty */
776 if ( gserialized_is_empty(geom1) )
777 PG_RETURN_POINTER(geom1);
778
779 srid = gserialized_get_srid(geom1);
780
781 initGEOS(lwpgnotice, lwgeom_geos_error);
782
783 g1 = POSTGIS2GEOS(geom1);
784
785 if (!g1)
786 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
787
788 g3 = GEOSConvexHull(g1);
789 GEOSGeom_destroy(g1);
790
791 if (!g3) HANDLE_GEOS_ERROR("GEOSConvexHull");
792
793 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
794
795 GEOSSetSRID(g3, srid);
796
797 lwout = GEOS2LWGEOM(g3, gserialized_has_z(geom1));
798 GEOSGeom_destroy(g3);
799
800 if (!lwout)
801 {
802 elog(ERROR,
803 "convexhull() failed to convert GEOS geometry to LWGEOM");
804 PG_RETURN_NULL(); /* never get here */
805 }
806
807 /* Copy input bbox if any */
808 if ( gserialized_get_gbox_p(geom1, &bbox) )
809 {
810 /* Force the box to have the same dimensionality as the lwgeom */
811 bbox.flags = lwout->flags;
812 lwout->bbox = gbox_copy(&bbox);
813 }
814
815 result = geometry_serialize(lwout);
816 lwgeom_free(lwout);
817
818 if (!result)
819 {
820 elog(ERROR,"GEOS convexhull() threw an error (result postgis geometry formation)!");
821 PG_RETURN_NULL(); /* never get here */
822 }
823
824 PG_FREE_IF_COPY(geom1, 0);
825 PG_RETURN_POINTER(result);
826 }
827
828 PG_FUNCTION_INFO_V1(topologypreservesimplify);
topologypreservesimplify(PG_FUNCTION_ARGS)829 Datum topologypreservesimplify(PG_FUNCTION_ARGS)
830 {
831 GSERIALIZED *geom1;
832 double tolerance;
833 GEOSGeometry *g1, *g3;
834 GSERIALIZED *result;
835 uint32_t type;
836
837 geom1 = PG_GETARG_GSERIALIZED_P(0);
838 tolerance = PG_GETARG_FLOAT8(1);
839
840 /* Empty.Simplify() == Empty */
841 type = gserialized_get_type(geom1);
842 if ( gserialized_is_empty(geom1) || type == TINTYPE || type == TRIANGLETYPE )
843 PG_RETURN_POINTER(geom1);
844
845 initGEOS(lwpgnotice, lwgeom_geos_error);
846
847 g1 = POSTGIS2GEOS(geom1);
848 if (!g1)
849 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
850
851 g3 = GEOSTopologyPreserveSimplify(g1,tolerance);
852 GEOSGeom_destroy(g1);
853
854 if (!g3) HANDLE_GEOS_ERROR("GEOSTopologyPreserveSimplify");
855
856 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
857
858 GEOSSetSRID(g3, gserialized_get_srid(geom1));
859
860 result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
861 GEOSGeom_destroy(g3);
862
863 if (!result)
864 {
865 elog(ERROR,"GEOS topologypreservesimplify() threw an error (result postgis geometry formation)!");
866 PG_RETURN_NULL(); /* never get here */
867 }
868
869 PG_FREE_IF_COPY(geom1, 0);
870 PG_RETURN_POINTER(result);
871 }
872
873 PG_FUNCTION_INFO_V1(buffer);
buffer(PG_FUNCTION_ARGS)874 Datum buffer(PG_FUNCTION_ARGS)
875 {
876 GEOSBufferParams *bufferparams;
877 GEOSGeometry *g1, *g3 = NULL;
878 GSERIALIZED *result;
879 LWGEOM *lwg;
880 int quadsegs = 8; /* the default */
881 int singleside = 0; /* the default */
882 enum
883 {
884 ENDCAP_ROUND = 1,
885 ENDCAP_FLAT = 2,
886 ENDCAP_SQUARE = 3
887 };
888 enum
889 {
890 JOIN_ROUND = 1,
891 JOIN_MITRE = 2,
892 JOIN_BEVEL = 3
893 };
894 const double DEFAULT_MITRE_LIMIT = 5.0;
895 const int DEFAULT_ENDCAP_STYLE = ENDCAP_ROUND;
896 const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
897 double mitreLimit = DEFAULT_MITRE_LIMIT;
898 int endCapStyle = DEFAULT_ENDCAP_STYLE;
899 int joinStyle = DEFAULT_JOIN_STYLE;
900
901 GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
902 double size = PG_GETARG_FLOAT8(1);
903 text *params_text;
904
905 if (PG_NARGS() > 2)
906 {
907 params_text = PG_GETARG_TEXT_P(2);
908 }
909 else
910 {
911 params_text = palloc(VARHDRSZ);
912 SET_VARSIZE(params_text, 0);
913 }
914
915 /* Empty.Buffer() == Empty[polygon] */
916 if ( gserialized_is_empty(geom1) )
917 {
918 lwg = lwpoly_as_lwgeom(lwpoly_construct_empty(
919 gserialized_get_srid(geom1),
920 0, 0)); // buffer wouldn't give back z or m anyway
921 PG_RETURN_POINTER(geometry_serialize(lwg));
922 }
923
924 initGEOS(lwpgnotice, lwgeom_geos_error);
925
926 g1 = POSTGIS2GEOS(geom1);
927 if (!g1)
928 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
929
930
931 if (VARSIZE_ANY_EXHDR(params_text) > 0)
932 {
933 char *param;
934 char *params = text_to_cstring(params_text);
935
936 for (param=params; ; param=NULL)
937 {
938 char *key, *val;
939 param = strtok(param, " ");
940 if (!param) break;
941 POSTGIS_DEBUGF(3, "Param: %s", param);
942
943 key = param;
944 val = strchr(key, '=');
945 if (!val || *(val + 1) == '\0')
946 {
947 lwpgerror("Missing value for buffer parameter %s", key);
948 break;
949 }
950 *val = '\0';
951 ++val;
952
953 POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
954
955 if ( !strcmp(key, "endcap") )
956 {
957 /* Supported end cap styles:
958 * "round", "flat", "square"
959 */
960 if ( !strcmp(val, "round") )
961 {
962 endCapStyle = ENDCAP_ROUND;
963 }
964 else if ( !strcmp(val, "flat") ||
965 !strcmp(val, "butt") )
966 {
967 endCapStyle = ENDCAP_FLAT;
968 }
969 else if ( !strcmp(val, "square") )
970 {
971 endCapStyle = ENDCAP_SQUARE;
972 }
973 else
974 {
975 lwpgerror("Invalid buffer end cap "
976 "style: %s (accept: "
977 "'round', 'flat', 'butt' "
978 "or 'square'"
979 ")", val);
980 break;
981 }
982
983 }
984 else if ( !strcmp(key, "join") )
985 {
986 if ( !strcmp(val, "round") )
987 {
988 joinStyle = JOIN_ROUND;
989 }
990 else if ( !strcmp(val, "mitre") ||
991 !strcmp(val, "miter") )
992 {
993 joinStyle = JOIN_MITRE;
994 }
995 else if ( !strcmp(val, "bevel") )
996 {
997 joinStyle = JOIN_BEVEL;
998 }
999 else
1000 {
1001 lwpgerror("Invalid buffer end cap "
1002 "style: %s (accept: "
1003 "'round', 'mitre', 'miter' "
1004 " or 'bevel'"
1005 ")", val);
1006 break;
1007 }
1008 }
1009 else if ( !strcmp(key, "mitre_limit") ||
1010 !strcmp(key, "miter_limit") )
1011 {
1012 /* mitreLimit is a float */
1013 mitreLimit = atof(val);
1014 }
1015 else if ( !strcmp(key, "quad_segs") )
1016 {
1017 /* quadrant segments is an int */
1018 quadsegs = atoi(val);
1019 }
1020 else if ( !strcmp(key, "side") )
1021 {
1022 if ( !strcmp(val, "both") )
1023 {
1024 singleside = 0;
1025 }
1026 else if ( !strcmp(val, "left") )
1027 {
1028 singleside = 1;
1029 }
1030 else if ( !strcmp(val, "right") )
1031 {
1032 singleside = 1;
1033 size *= -1;
1034 }
1035 else
1036 {
1037 lwpgerror("Invalid side parameter: %s (accept: 'right', 'left', 'both')", val);
1038 break;
1039 }
1040 }
1041 else
1042 {
1043 lwpgerror(
1044 "Invalid buffer parameter: %s (accept: 'endcap', 'join', 'mitre_limit', 'miter_limit', 'quad_segs' and 'side')",
1045 key);
1046 break;
1047 }
1048 }
1049 pfree(params); /* was pstrduped */
1050 }
1051
1052
1053 POSTGIS_DEBUGF(3, "endCap:%d joinStyle:%d mitreLimit:%g",
1054 endCapStyle, joinStyle, mitreLimit);
1055
1056 bufferparams = GEOSBufferParams_create();
1057 if (bufferparams)
1058 {
1059 if (GEOSBufferParams_setEndCapStyle(bufferparams, endCapStyle) &&
1060 GEOSBufferParams_setJoinStyle(bufferparams, joinStyle) &&
1061 GEOSBufferParams_setMitreLimit(bufferparams, mitreLimit) &&
1062 GEOSBufferParams_setQuadrantSegments(bufferparams, quadsegs) &&
1063 GEOSBufferParams_setSingleSided(bufferparams, singleside))
1064 {
1065 g3 = GEOSBufferWithParams(g1, bufferparams, size);
1066 }
1067 else
1068 {
1069 lwpgerror("Error setting buffer parameters.");
1070 }
1071 GEOSBufferParams_destroy(bufferparams);
1072 }
1073 else
1074 {
1075 lwpgerror("Error setting buffer parameters.");
1076 }
1077
1078 GEOSGeom_destroy(g1);
1079
1080 if (!g3) HANDLE_GEOS_ERROR("GEOSBuffer");
1081
1082 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1083
1084 GEOSSetSRID(g3, gserialized_get_srid(geom1));
1085
1086 result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
1087 GEOSGeom_destroy(g3);
1088
1089 if (!result)
1090 {
1091 elog(ERROR,"GEOS buffer() threw an error (result postgis geometry formation)!");
1092 PG_RETURN_NULL(); /* never get here */
1093 }
1094
1095 PG_FREE_IF_COPY(geom1, 0);
1096 PG_RETURN_POINTER(result);
1097 }
1098
1099 /*
1100 * Generate a field of random points within the area of a
1101 * polygon or multipolygon. Throws an error for other geometry
1102 * types.
1103 */
1104 Datum ST_GeneratePoints(PG_FUNCTION_ARGS);
1105 PG_FUNCTION_INFO_V1(ST_GeneratePoints);
ST_GeneratePoints(PG_FUNCTION_ARGS)1106 Datum ST_GeneratePoints(PG_FUNCTION_ARGS)
1107 {
1108 GSERIALIZED *gser_input;
1109 GSERIALIZED *gser_result;
1110 LWGEOM *lwgeom_input;
1111 LWGEOM *lwgeom_result;
1112 int32 npoints;
1113 int32 seed = 0;
1114
1115 gser_input = PG_GETARG_GSERIALIZED_P(0);
1116 npoints = PG_GETARG_INT32(1);
1117
1118 if (npoints < 0)
1119 PG_RETURN_NULL();
1120
1121 if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
1122 {
1123 seed = PG_GETARG_INT32(2);
1124 if (seed < 1)
1125 {
1126 lwpgerror("ST_GeneratePoints: seed must be greater than zero");
1127 PG_RETURN_NULL();
1128 }
1129 }
1130
1131 /* Types get checked in the code, we'll keep things small out there */
1132 lwgeom_input = lwgeom_from_gserialized(gser_input);
1133 lwgeom_result = (LWGEOM*)lwgeom_to_points(lwgeom_input, npoints, seed);
1134 lwgeom_free(lwgeom_input);
1135 PG_FREE_IF_COPY(gser_input, 0);
1136
1137 /* Return null as null */
1138 if (!lwgeom_result)
1139 PG_RETURN_NULL();
1140
1141 /* Serialize and return */
1142 gser_result = geometry_serialize(lwgeom_result);
1143 lwgeom_free(lwgeom_result);
1144 PG_RETURN_POINTER(gser_result);
1145 }
1146
1147
1148 /*
1149 * Compute at offset curve to a line
1150 */
1151 Datum ST_OffsetCurve(PG_FUNCTION_ARGS);
1152 PG_FUNCTION_INFO_V1(ST_OffsetCurve);
ST_OffsetCurve(PG_FUNCTION_ARGS)1153 Datum ST_OffsetCurve(PG_FUNCTION_ARGS)
1154 {
1155 GSERIALIZED *gser_input;
1156 GSERIALIZED *gser_result;
1157 LWGEOM *lwgeom_input;
1158 LWGEOM *lwgeom_result;
1159 double size;
1160 int quadsegs = 8; /* the default */
1161 int nargs;
1162
1163 enum
1164 {
1165 JOIN_ROUND = 1,
1166 JOIN_MITRE = 2,
1167 JOIN_BEVEL = 3
1168 };
1169
1170 static const double DEFAULT_MITRE_LIMIT = 5.0;
1171 static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
1172 double mitreLimit = DEFAULT_MITRE_LIMIT;
1173 int joinStyle = DEFAULT_JOIN_STYLE;
1174 char *param = NULL;
1175 char *paramstr = NULL;
1176
1177 /* Read SQL arguments */
1178 nargs = PG_NARGS();
1179 gser_input = PG_GETARG_GSERIALIZED_P(0);
1180 size = PG_GETARG_FLOAT8(1);
1181
1182 /* For distance == 0, just return the input. */
1183 if (size == 0) PG_RETURN_POINTER(gser_input);
1184
1185 /* Read the lwgeom, check for errors */
1186 lwgeom_input = lwgeom_from_gserialized(gser_input);
1187 if ( ! lwgeom_input )
1188 lwpgerror("ST_OffsetCurve: lwgeom_from_gserialized returned NULL");
1189
1190 /* For empty inputs, just echo them back */
1191 if ( lwgeom_is_empty(lwgeom_input) )
1192 PG_RETURN_POINTER(gser_input);
1193
1194 /* Process the optional arguments */
1195 if ( nargs > 2 )
1196 {
1197 text *wkttext = PG_GETARG_TEXT_P(2);
1198 paramstr = text_to_cstring(wkttext);
1199
1200 POSTGIS_DEBUGF(3, "paramstr: %s", paramstr);
1201
1202 for ( param=paramstr; ; param=NULL )
1203 {
1204 char *key, *val;
1205 param = strtok(param, " ");
1206 if (!param) break;
1207 POSTGIS_DEBUGF(3, "Param: %s", param);
1208
1209 key = param;
1210 val = strchr(key, '=');
1211 if (!val || *(val + 1) == '\0')
1212 {
1213 lwpgerror("ST_OffsetCurve: Missing value for buffer parameter %s", key);
1214 break;
1215 }
1216 *val = '\0';
1217 ++val;
1218
1219 POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
1220
1221 if ( !strcmp(key, "join") )
1222 {
1223 if ( !strcmp(val, "round") )
1224 {
1225 joinStyle = JOIN_ROUND;
1226 }
1227 else if ( !(strcmp(val, "mitre") && strcmp(val, "miter")) )
1228 {
1229 joinStyle = JOIN_MITRE;
1230 }
1231 else if ( ! strcmp(val, "bevel") )
1232 {
1233 joinStyle = JOIN_BEVEL;
1234 }
1235 else
1236 {
1237 lwpgerror(
1238 "Invalid buffer end cap style: %s (accept: 'round', 'mitre', 'miter' or 'bevel')",
1239 val);
1240 break;
1241 }
1242 }
1243 else if ( !strcmp(key, "mitre_limit") ||
1244 !strcmp(key, "miter_limit") )
1245 {
1246 /* mitreLimit is a float */
1247 mitreLimit = atof(val);
1248 }
1249 else if ( !strcmp(key, "quad_segs") )
1250 {
1251 /* quadrant segments is an int */
1252 quadsegs = atoi(val);
1253 }
1254 else
1255 {
1256 lwpgerror(
1257 "Invalid buffer parameter: %s (accept: 'join', 'mitre_limit', 'miter_limit and 'quad_segs')",
1258 key);
1259 break;
1260 }
1261 }
1262 POSTGIS_DEBUGF(3, "joinStyle:%d mitreLimit:%g", joinStyle, mitreLimit);
1263 pfree(paramstr); /* alloc'ed in text_to_cstring */
1264 }
1265
1266 lwgeom_result = lwgeom_offsetcurve(lwgeom_input, size, quadsegs, joinStyle, mitreLimit);
1267
1268 if (!lwgeom_result)
1269 lwpgerror("ST_OffsetCurve: lwgeom_offsetcurve returned NULL");
1270
1271 gser_result = geometry_serialize(lwgeom_result);
1272 lwgeom_free(lwgeom_input);
1273 lwgeom_free(lwgeom_result);
1274 PG_RETURN_POINTER(gser_result);
1275 }
1276
1277 PG_FUNCTION_INFO_V1(ST_Intersection);
ST_Intersection(PG_FUNCTION_ARGS)1278 Datum ST_Intersection(PG_FUNCTION_ARGS)
1279 {
1280 GSERIALIZED *geom1;
1281 GSERIALIZED *geom2;
1282 GSERIALIZED *result;
1283 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
1284
1285 geom1 = PG_GETARG_GSERIALIZED_P(0);
1286 geom2 = PG_GETARG_GSERIALIZED_P(1);
1287
1288 lwgeom1 = lwgeom_from_gserialized(geom1);
1289 lwgeom2 = lwgeom_from_gserialized(geom2);
1290
1291 lwresult = lwgeom_intersection(lwgeom1, lwgeom2);
1292 result = geometry_serialize(lwresult);
1293
1294 lwgeom_free(lwgeom1);
1295 lwgeom_free(lwgeom2);
1296 lwgeom_free(lwresult);
1297
1298 PG_FREE_IF_COPY(geom1, 0);
1299 PG_FREE_IF_COPY(geom2, 1);
1300
1301 PG_RETURN_POINTER(result);
1302 }
1303
1304 PG_FUNCTION_INFO_V1(ST_Difference);
ST_Difference(PG_FUNCTION_ARGS)1305 Datum ST_Difference(PG_FUNCTION_ARGS)
1306 {
1307 GSERIALIZED *geom1;
1308 GSERIALIZED *geom2;
1309 GSERIALIZED *result;
1310 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
1311
1312 geom1 = PG_GETARG_GSERIALIZED_P(0);
1313 geom2 = PG_GETARG_GSERIALIZED_P(1);
1314
1315 lwgeom1 = lwgeom_from_gserialized(geom1);
1316 lwgeom2 = lwgeom_from_gserialized(geom2);
1317
1318 lwresult = lwgeom_difference(lwgeom1, lwgeom2);
1319 result = geometry_serialize(lwresult);
1320
1321 lwgeom_free(lwgeom1);
1322 lwgeom_free(lwgeom2);
1323 lwgeom_free(lwresult);
1324
1325 PG_FREE_IF_COPY(geom1, 0);
1326 PG_FREE_IF_COPY(geom2, 1);
1327
1328 PG_RETURN_POINTER(result);
1329 }
1330
1331 /**
1332 @example pointonsurface - {@link #pointonsurface} SELECT ST_PointOnSurface('POLYGON((0 0, 10 0, 10 10, 0 10, 0
1333 0))');
1334 */
1335 PG_FUNCTION_INFO_V1(pointonsurface);
pointonsurface(PG_FUNCTION_ARGS)1336 Datum pointonsurface(PG_FUNCTION_ARGS)
1337 {
1338 GSERIALIZED *geom, *result;
1339 LWGEOM *lwgeom, *lwresult;
1340
1341 geom = PG_GETARG_GSERIALIZED_P(0);
1342
1343 lwgeom = lwgeom_from_gserialized(geom);
1344 lwresult = lwgeom_pointonsurface(lwgeom);
1345 lwgeom_free(lwgeom);
1346 PG_FREE_IF_COPY(geom, 0);
1347
1348 if (!lwresult) PG_RETURN_NULL();
1349
1350 result = geometry_serialize(lwresult);
1351 lwgeom_free(lwresult);
1352 PG_RETURN_POINTER(result);
1353 }
1354
1355 PG_FUNCTION_INFO_V1(centroid);
centroid(PG_FUNCTION_ARGS)1356 Datum centroid(PG_FUNCTION_ARGS)
1357 {
1358 GSERIALIZED *geom, *result;
1359 LWGEOM *lwgeom, *lwresult;
1360
1361 geom = PG_GETARG_GSERIALIZED_P(0);
1362
1363 lwgeom = lwgeom_from_gserialized(geom);
1364 lwresult = lwgeom_centroid(lwgeom);
1365 lwgeom_free(lwgeom);
1366 PG_FREE_IF_COPY(geom, 0);
1367
1368 if (!lwresult) PG_RETURN_NULL();
1369
1370 result = geometry_serialize(lwresult);
1371 lwgeom_free(lwresult);
1372 PG_RETURN_POINTER(result);
1373 }
1374
1375 Datum ST_ClipByBox2d(PG_FUNCTION_ARGS);
1376 PG_FUNCTION_INFO_V1(ST_ClipByBox2d);
ST_ClipByBox2d(PG_FUNCTION_ARGS)1377 Datum ST_ClipByBox2d(PG_FUNCTION_ARGS)
1378 {
1379 GSERIALIZED *geom1;
1380 GSERIALIZED *result;
1381 LWGEOM *lwgeom1, *lwresult ;
1382 const GBOX *bbox1;
1383 GBOX *bbox2;
1384
1385 geom1 = PG_GETARG_GSERIALIZED_P(0);
1386 lwgeom1 = lwgeom_from_gserialized(geom1) ;
1387
1388 bbox1 = lwgeom_get_bbox(lwgeom1);
1389 if ( ! bbox1 )
1390 {
1391 /* empty clips to empty, no matter rect */
1392 lwgeom_free(lwgeom1);
1393 PG_RETURN_POINTER(geom1);
1394 }
1395
1396 /* WARNING: this is really a BOX2DF, use only xmin and ymin fields */
1397 bbox2 = (GBOX *)PG_GETARG_POINTER(1);
1398 bbox2->flags = 0;
1399
1400 /* If bbox1 outside of bbox2, return empty */
1401 if ( ! gbox_overlaps_2d(bbox1, bbox2) )
1402 {
1403 lwresult = lwgeom_construct_empty(lwgeom1->type, lwgeom1->srid, 0, 0);
1404 lwgeom_free(lwgeom1);
1405 PG_FREE_IF_COPY(geom1, 0);
1406 result = geometry_serialize(lwresult) ;
1407 lwgeom_free(lwresult) ;
1408 PG_RETURN_POINTER(result);
1409 }
1410
1411 /* if bbox1 is covered by bbox2, return lwgeom1 */
1412 if ( gbox_contains_2d(bbox2, bbox1) )
1413 {
1414 lwgeom_free(lwgeom1);
1415 PG_RETURN_POINTER(geom1);
1416 }
1417
1418 lwresult = lwgeom_clip_by_rect(lwgeom1, bbox2->xmin, bbox2->ymin,
1419 bbox2->xmax, bbox2->ymax);
1420
1421 lwgeom_free(lwgeom1);
1422 PG_FREE_IF_COPY(geom1, 0);
1423
1424 if (!lwresult) PG_RETURN_NULL();
1425
1426 result = geometry_serialize(lwresult) ;
1427 lwgeom_free(lwresult) ;
1428 PG_RETURN_POINTER(result);
1429 }
1430
1431
1432 /*---------------------------------------------*/
1433
1434 PG_FUNCTION_INFO_V1(isvalid);
isvalid(PG_FUNCTION_ARGS)1435 Datum isvalid(PG_FUNCTION_ARGS)
1436 {
1437 GSERIALIZED *geom1;
1438 LWGEOM *lwgeom;
1439 char result;
1440 GEOSGeom g1;
1441
1442 geom1 = PG_GETARG_GSERIALIZED_P(0);
1443
1444 /* Empty.IsValid() == TRUE */
1445 if ( gserialized_is_empty(geom1) )
1446 PG_RETURN_BOOL(true);
1447
1448 initGEOS(lwpgnotice, lwgeom_geos_error);
1449
1450 lwgeom = lwgeom_from_gserialized(geom1);
1451 if ( ! lwgeom )
1452 {
1453 lwpgerror("unable to deserialize input");
1454 }
1455 g1 = LWGEOM2GEOS(lwgeom, 0);
1456 lwgeom_free(lwgeom);
1457
1458 if ( ! g1 )
1459 {
1460 /* should we drop the following
1461 * notice now that we have ST_isValidReason ?
1462 */
1463 lwpgnotice("%s", lwgeom_geos_errmsg);
1464 PG_RETURN_BOOL(false);
1465 }
1466
1467 result = GEOSisValid(g1);
1468 GEOSGeom_destroy(g1);
1469
1470 if (result == 2)
1471 {
1472 elog(ERROR,"GEOS isvalid() threw an error!");
1473 PG_RETURN_NULL(); /*never get here */
1474 }
1475
1476 PG_FREE_IF_COPY(geom1, 0);
1477 PG_RETURN_BOOL(result);
1478 }
1479
1480 PG_FUNCTION_INFO_V1(isvalidreason);
isvalidreason(PG_FUNCTION_ARGS)1481 Datum isvalidreason(PG_FUNCTION_ARGS)
1482 {
1483 GSERIALIZED *geom = NULL;
1484 char *reason_str = NULL;
1485 text *result = NULL;
1486 const GEOSGeometry *g1 = NULL;
1487
1488 geom = PG_GETARG_GSERIALIZED_P(0);
1489
1490 initGEOS(lwpgnotice, lwgeom_geos_error);
1491
1492 g1 = POSTGIS2GEOS(geom);
1493 if ( g1 )
1494 {
1495 reason_str = GEOSisValidReason(g1);
1496 GEOSGeom_destroy((GEOSGeometry *)g1);
1497 if (!reason_str) HANDLE_GEOS_ERROR("GEOSisValidReason");
1498 result = cstring_to_text(reason_str);
1499 GEOSFree(reason_str);
1500 }
1501 else
1502 {
1503 result = cstring_to_text(lwgeom_geos_errmsg);
1504 }
1505
1506 PG_FREE_IF_COPY(geom, 0);
1507 PG_RETURN_POINTER(result);
1508 }
1509
1510 PG_FUNCTION_INFO_V1(isvaliddetail);
isvaliddetail(PG_FUNCTION_ARGS)1511 Datum isvaliddetail(PG_FUNCTION_ARGS)
1512 {
1513 GSERIALIZED *geom = NULL;
1514 const GEOSGeometry *g1 = NULL;
1515 char *values[3]; /* valid bool, reason text, location geometry */
1516 char *geos_reason = NULL;
1517 char *reason = NULL;
1518 GEOSGeometry *geos_location = NULL;
1519 LWGEOM *location = NULL;
1520 char valid = 0;
1521 HeapTupleHeader result;
1522 TupleDesc tupdesc;
1523 HeapTuple tuple;
1524 AttInMetadata *attinmeta;
1525 int flags = 0;
1526
1527 /*
1528 * Build a tuple description for a
1529 * valid_detail tuple
1530 */
1531 get_call_result_type(fcinfo, 0, &tupdesc);
1532 BlessTupleDesc(tupdesc);
1533
1534 /*
1535 * generate attribute metadata needed later to produce
1536 * tuples from raw C strings
1537 */
1538 attinmeta = TupleDescGetAttInMetadata(tupdesc);
1539
1540 geom = PG_GETARG_GSERIALIZED_P(0);
1541 flags = PG_GETARG_INT32(1);
1542
1543 initGEOS(lwpgnotice, lwgeom_geos_error);
1544
1545 g1 = POSTGIS2GEOS(geom);
1546
1547 if ( g1 )
1548 {
1549 valid = GEOSisValidDetail(g1, flags, &geos_reason, &geos_location);
1550 GEOSGeom_destroy((GEOSGeometry *)g1);
1551 if ( geos_reason )
1552 {
1553 reason = pstrdup(geos_reason);
1554 GEOSFree(geos_reason);
1555 }
1556 if ( geos_location )
1557 {
1558 location = GEOS2LWGEOM(geos_location, GEOSHasZ(geos_location));
1559 GEOSGeom_destroy(geos_location);
1560 }
1561
1562 if (valid == 2)
1563 {
1564 /* NOTE: should only happen on OOM or similar */
1565 lwpgerror("GEOS isvaliddetail() threw an exception!");
1566 PG_RETURN_NULL(); /* never gets here */
1567 }
1568 }
1569 else
1570 {
1571 /* TODO: check lwgeom_geos_errmsg for validity error */
1572 reason = pstrdup(lwgeom_geos_errmsg);
1573 }
1574
1575 /* the boolean validity */
1576 values[0] = valid ? "t" : "f";
1577
1578 /* the reason */
1579 values[1] = reason;
1580
1581 /* the location */
1582 values[2] = location ? lwgeom_to_hexwkb(location, WKB_EXTENDED, 0) : 0;
1583
1584 tuple = BuildTupleFromCStrings(attinmeta, values);
1585 result = (HeapTupleHeader) palloc(tuple->t_len);
1586 memcpy(result, tuple->t_data, tuple->t_len);
1587 heap_freetuple(tuple);
1588
1589 PG_RETURN_HEAPTUPLEHEADER(result);
1590 }
1591
1592 /**
1593 * overlaps(GSERIALIZED g1,GSERIALIZED g2)
1594 * @param g1
1595 * @param g2
1596 * @return if GEOS::g1->overlaps(g2) returns true
1597 * @throw an error (elog(ERROR,...)) if GEOS throws an error
1598 */
1599 PG_FUNCTION_INFO_V1(overlaps);
overlaps(PG_FUNCTION_ARGS)1600 Datum overlaps(PG_FUNCTION_ARGS)
1601 {
1602 GSERIALIZED *geom1;
1603 GSERIALIZED *geom2;
1604 GEOSGeometry *g1, *g2;
1605 char result;
1606 GBOX box1, box2;
1607
1608 geom1 = PG_GETARG_GSERIALIZED_P(0);
1609 geom2 = PG_GETARG_GSERIALIZED_P(1);
1610 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1611
1612 /* A.Overlaps(Empty) == FALSE */
1613 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1614 PG_RETURN_BOOL(false);
1615
1616 /*
1617 * short-circuit 1: if geom2 bounding box does not overlap
1618 * geom1 bounding box we can return FALSE.
1619 */
1620 if ( gserialized_get_gbox_p(geom1, &box1) &&
1621 gserialized_get_gbox_p(geom2, &box2) )
1622 {
1623 if ( ! gbox_overlaps_2d(&box1, &box2) )
1624 {
1625 PG_RETURN_BOOL(false);
1626 }
1627 }
1628
1629 initGEOS(lwpgnotice, lwgeom_geos_error);
1630
1631 g1 = POSTGIS2GEOS(geom1);
1632 if (!g1)
1633 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1634
1635 g2 = POSTGIS2GEOS(geom2);
1636
1637 if (!g2)
1638 {
1639 GEOSGeom_destroy(g1);
1640 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1641 }
1642
1643 result = GEOSOverlaps(g1,g2);
1644
1645 GEOSGeom_destroy(g1);
1646 GEOSGeom_destroy(g2);
1647 if (result == 2) HANDLE_GEOS_ERROR("GEOSOverlaps");
1648
1649 PG_FREE_IF_COPY(geom1, 0);
1650 PG_FREE_IF_COPY(geom2, 1);
1651
1652 PG_RETURN_BOOL(result);
1653 }
1654
1655
1656 PG_FUNCTION_INFO_V1(contains);
contains(PG_FUNCTION_ARGS)1657 Datum contains(PG_FUNCTION_ARGS)
1658 {
1659 GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
1660 GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
1661 int result;
1662 GEOSGeometry *g1, *g2;
1663 GBOX box1, box2;
1664 PrepGeomCache *prep_cache;
1665 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1666
1667 /* A.Contains(Empty) == FALSE */
1668 if (gserialized_is_empty(geom1) || gserialized_is_empty(geom2))
1669 PG_RETURN_BOOL(false);
1670
1671 POSTGIS_DEBUG(3, "contains called.");
1672
1673 /*
1674 ** short-circuit 1: if geom2 bounding box is not completely inside
1675 ** geom1 bounding box we can return FALSE.
1676 */
1677 if (gserialized_get_gbox_p(geom1, &box1) &&
1678 gserialized_get_gbox_p(geom2, &box2))
1679 {
1680 if (!gbox_contains_2d(&box1, &box2))
1681 PG_RETURN_BOOL(false);
1682 }
1683
1684 /*
1685 ** short-circuit 2: if geom2 is a point and geom1 is a polygon
1686 ** call the point-in-polygon function.
1687 */
1688 if (is_poly(geom1) && is_point(geom2))
1689 {
1690 GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
1691 GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
1692 RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
1693 int retval;
1694
1695 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
1696 if (gserialized_get_type(gpoint) == POINTTYPE)
1697 {
1698 LWGEOM* point = lwgeom_from_gserialized(gpoint);
1699 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
1700 lwgeom_free(point);
1701
1702 retval = (pip_result == 1); /* completely inside */
1703 }
1704 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
1705 {
1706 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
1707 uint32_t i;
1708 int found_completely_inside = LW_FALSE;
1709
1710 retval = LW_TRUE;
1711 for (i = 0; i < mpoint->ngeoms; i++)
1712 {
1713 /* We need to find at least one point that's completely inside the
1714 * polygons (pip_result == 1). As long as we have one point that's
1715 * completely inside, we can have as many as we want on the boundary
1716 * itself. (pip_result == 0)
1717 */
1718 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
1719 if (pip_result == 1)
1720 found_completely_inside = LW_TRUE;
1721
1722 if (pip_result == -1) /* completely outside */
1723 {
1724 retval = LW_FALSE;
1725 break;
1726 }
1727 }
1728
1729 retval = retval && found_completely_inside;
1730 lwmpoint_free(mpoint);
1731 }
1732 else
1733 {
1734 /* Never get here */
1735 elog(ERROR,"Type isn't point or multipoint!");
1736 PG_RETURN_BOOL(false);
1737 }
1738
1739 return retval > 0;
1740 }
1741 else
1742 {
1743 POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
1744 }
1745
1746 initGEOS(lwpgnotice, lwgeom_geos_error);
1747
1748 prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
1749
1750 if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
1751 {
1752 g1 = POSTGIS2GEOS(geom2);
1753 if (!g1)
1754 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
1755
1756 POSTGIS_DEBUG(4, "containsPrepared: cache is live, running preparedcontains");
1757 result = GEOSPreparedContains( prep_cache->prepared_geom, g1);
1758 GEOSGeom_destroy(g1);
1759 }
1760 else
1761 {
1762 g1 = POSTGIS2GEOS(geom1);
1763 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1764 g2 = POSTGIS2GEOS(geom2);
1765 if (!g2)
1766 {
1767 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1768 GEOSGeom_destroy(g1);
1769 }
1770 POSTGIS_DEBUG(4, "containsPrepared: cache is not ready, running standard contains");
1771 result = GEOSContains( g1, g2);
1772 GEOSGeom_destroy(g1);
1773 GEOSGeom_destroy(g2);
1774 }
1775
1776 if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
1777
1778 PG_FREE_IF_COPY(geom1, 0);
1779 PG_FREE_IF_COPY(geom2, 1);
1780 PG_RETURN_BOOL(result > 0);
1781 }
1782
1783
1784 PG_FUNCTION_INFO_V1(within);
within(PG_FUNCTION_ARGS)1785 Datum within(PG_FUNCTION_ARGS)
1786 {
1787 PG_RETURN_DATUM(CallerFInfoFunctionCall2(contains, fcinfo->flinfo, InvalidOid,
1788 PG_GETARG_DATUM(1), PG_GETARG_DATUM(0)));
1789 }
1790
1791
1792
1793 PG_FUNCTION_INFO_V1(containsproperly);
containsproperly(PG_FUNCTION_ARGS)1794 Datum containsproperly(PG_FUNCTION_ARGS)
1795 {
1796 GSERIALIZED * geom1;
1797 GSERIALIZED * geom2;
1798 char result;
1799 GBOX box1, box2;
1800 PrepGeomCache * prep_cache;
1801
1802 geom1 = PG_GETARG_GSERIALIZED_P(0);
1803 geom2 = PG_GETARG_GSERIALIZED_P(1);
1804 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1805
1806 /* A.ContainsProperly(Empty) == FALSE */
1807 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1808 PG_RETURN_BOOL(false);
1809
1810 /*
1811 * short-circuit: if geom2 bounding box is not completely inside
1812 * geom1 bounding box we can return FALSE.
1813 */
1814 if ( gserialized_get_gbox_p(geom1, &box1) &&
1815 gserialized_get_gbox_p(geom2, &box2) )
1816 {
1817 if ( ! gbox_contains_2d(&box1, &box2) )
1818 PG_RETURN_BOOL(false);
1819 }
1820
1821 initGEOS(lwpgnotice, lwgeom_geos_error);
1822
1823 prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
1824
1825 if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
1826 {
1827 GEOSGeometry *g = POSTGIS2GEOS(geom2);
1828 if (!g) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1829 result = GEOSPreparedContainsProperly( prep_cache->prepared_geom, g);
1830 GEOSGeom_destroy(g);
1831 }
1832 else
1833 {
1834 GEOSGeometry *g2;
1835 GEOSGeometry *g1;
1836
1837 g1 = POSTGIS2GEOS(geom1);
1838 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1839 g2 = POSTGIS2GEOS(geom2);
1840 if (!g2)
1841 {
1842 GEOSGeom_destroy(g1);
1843 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1844 }
1845 result = GEOSRelatePattern( g1, g2, "T**FF*FF*" );
1846 GEOSGeom_destroy(g1);
1847 GEOSGeom_destroy(g2);
1848 }
1849
1850 if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
1851
1852 PG_FREE_IF_COPY(geom1, 0);
1853 PG_FREE_IF_COPY(geom2, 1);
1854
1855 PG_RETURN_BOOL(result);
1856 }
1857
1858 /*
1859 * Described at
1860 * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
1861 */
1862 PG_FUNCTION_INFO_V1(covers);
covers(PG_FUNCTION_ARGS)1863 Datum covers(PG_FUNCTION_ARGS)
1864 {
1865 GSERIALIZED *geom1;
1866 GSERIALIZED *geom2;
1867 int result;
1868 GBOX box1, box2;
1869 PrepGeomCache *prep_cache;
1870
1871 geom1 = PG_GETARG_GSERIALIZED_P(0);
1872 geom2 = PG_GETARG_GSERIALIZED_P(1);
1873
1874 /* A.Covers(Empty) == FALSE */
1875 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1876 PG_RETURN_BOOL(false);
1877
1878 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1879
1880 /*
1881 * short-circuit 1: if geom2 bounding box is not completely inside
1882 * geom1 bounding box we can return FALSE.
1883 */
1884 if ( gserialized_get_gbox_p(geom1, &box1) &&
1885 gserialized_get_gbox_p(geom2, &box2) )
1886 {
1887 if ( ! gbox_contains_2d(&box1, &box2) )
1888 {
1889 PG_RETURN_BOOL(false);
1890 }
1891 }
1892 /*
1893 * short-circuit 2: if geom2 is a point and geom1 is a polygon
1894 * call the point-in-polygon function.
1895 */
1896 if (is_poly(geom1) && is_point(geom2))
1897 {
1898 GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
1899 GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
1900 RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
1901 int retval;
1902
1903 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
1904 if (gserialized_get_type(gpoint) == POINTTYPE)
1905 {
1906 LWGEOM* point = lwgeom_from_gserialized(gpoint);
1907 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
1908 lwgeom_free(point);
1909
1910 retval = (pip_result != -1); /* not outside */
1911 }
1912 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
1913 {
1914 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
1915 uint32_t i;
1916
1917 retval = LW_TRUE;
1918 for (i = 0; i < mpoint->ngeoms; i++)
1919 {
1920 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
1921 if (pip_result == -1)
1922 {
1923 retval = LW_FALSE;
1924 break;
1925 }
1926 }
1927
1928 lwmpoint_free(mpoint);
1929 }
1930 else
1931 {
1932 /* Never get here */
1933 elog(ERROR,"Type isn't point or multipoint!");
1934 PG_RETURN_NULL();
1935 }
1936
1937 PG_FREE_IF_COPY(geom1, 0);
1938 PG_FREE_IF_COPY(geom2, 1);
1939 PG_RETURN_BOOL(retval);
1940 }
1941 else
1942 {
1943 POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
1944 }
1945
1946 initGEOS(lwpgnotice, lwgeom_geos_error);
1947
1948 prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
1949
1950 if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
1951 {
1952 GEOSGeometry *g1 = POSTGIS2GEOS(geom2);
1953 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1954 result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
1955 GEOSGeom_destroy(g1);
1956 }
1957 else
1958 {
1959 GEOSGeometry *g1;
1960 GEOSGeometry *g2;
1961
1962 g1 = POSTGIS2GEOS(geom1);
1963 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1964 g2 = POSTGIS2GEOS(geom2);
1965 if (!g2)
1966 {
1967 GEOSGeom_destroy(g1);
1968 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1969 }
1970 result = GEOSRelatePattern( g1, g2, "******FF*" );
1971 GEOSGeom_destroy(g1);
1972 GEOSGeom_destroy(g2);
1973 }
1974
1975 if (result == 2) HANDLE_GEOS_ERROR("GEOSCovers");
1976
1977 PG_FREE_IF_COPY(geom1, 0);
1978 PG_FREE_IF_COPY(geom2, 1);
1979
1980 PG_RETURN_BOOL(result);
1981
1982 }
1983
1984
1985 /**
1986 * ST_Within(A, B) => ST_Contains(B, A) so we just delegate this calculation to the
1987 * Contains implementation.
1988 PG_FUNCTION_INFO_V1(within);
1989 Datum within(PG_FUNCTION_ARGS)
1990 */
1991
1992 /*
1993 * Described at:
1994 * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
1995 */
1996 PG_FUNCTION_INFO_V1(coveredby);
coveredby(PG_FUNCTION_ARGS)1997 Datum coveredby(PG_FUNCTION_ARGS)
1998 {
1999 GSERIALIZED *geom1;
2000 GSERIALIZED *geom2;
2001 GEOSGeometry *g1, *g2;
2002 int result;
2003 GBOX box1, box2;
2004 char *patt = "**F**F***";
2005
2006 geom1 = PG_GETARG_GSERIALIZED_P(0);
2007 geom2 = PG_GETARG_GSERIALIZED_P(1);
2008 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2009
2010 /* A.CoveredBy(Empty) == FALSE */
2011 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2012 PG_RETURN_BOOL(false);
2013
2014 /*
2015 * short-circuit 1: if geom1 bounding box is not completely inside
2016 * geom2 bounding box we can return FALSE.
2017 */
2018 if ( gserialized_get_gbox_p(geom1, &box1) &&
2019 gserialized_get_gbox_p(geom2, &box2) )
2020 {
2021 if ( ! gbox_contains_2d(&box2, &box1) )
2022 {
2023 PG_RETURN_BOOL(false);
2024 }
2025
2026 POSTGIS_DEBUG(3, "bounding box short-circuit missed.");
2027 }
2028 /*
2029 * short-circuit 2: if geom1 is a point and geom2 is a polygon
2030 * call the point-in-polygon function.
2031 */
2032 if (is_point(geom1) && is_poly(geom2))
2033 {
2034 GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
2035 GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
2036 RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
2037 int retval;
2038
2039 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2040 if (gserialized_get_type(gpoint) == POINTTYPE)
2041 {
2042 LWGEOM* point = lwgeom_from_gserialized(gpoint);
2043 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2044 lwgeom_free(point);
2045
2046 retval = (pip_result != -1); /* not outside */
2047 }
2048 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2049 {
2050 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2051 uint32_t i;
2052
2053 retval = LW_TRUE;
2054 for (i = 0; i < mpoint->ngeoms; i++)
2055 {
2056 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2057 if (pip_result == -1)
2058 {
2059 retval = LW_FALSE;
2060 break;
2061 }
2062 }
2063
2064 lwmpoint_free(mpoint);
2065 }
2066 else
2067 {
2068 /* Never get here */
2069 elog(ERROR,"Type isn't point or multipoint!");
2070 PG_RETURN_NULL();
2071 }
2072
2073 PG_FREE_IF_COPY(geom1, 0);
2074 PG_FREE_IF_COPY(geom2, 1);
2075 PG_RETURN_BOOL(retval);
2076 }
2077 else
2078 {
2079 POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2080 }
2081
2082 initGEOS(lwpgnotice, lwgeom_geos_error);
2083
2084 g1 = POSTGIS2GEOS(geom1);
2085
2086 if (!g1)
2087 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2088
2089 g2 = POSTGIS2GEOS(geom2);
2090
2091 if (!g2)
2092 {
2093 GEOSGeom_destroy(g1);
2094 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2095 }
2096
2097 result = GEOSRelatePattern(g1,g2,patt);
2098
2099 GEOSGeom_destroy(g1);
2100 GEOSGeom_destroy(g2);
2101
2102 if (result == 2) HANDLE_GEOS_ERROR("GEOSCoveredBy");
2103
2104 PG_FREE_IF_COPY(geom1, 0);
2105 PG_FREE_IF_COPY(geom2, 1);
2106
2107 PG_RETURN_BOOL(result);
2108 }
2109
2110 PG_FUNCTION_INFO_V1(crosses);
crosses(PG_FUNCTION_ARGS)2111 Datum crosses(PG_FUNCTION_ARGS)
2112 {
2113 GSERIALIZED *geom1;
2114 GSERIALIZED *geom2;
2115 GEOSGeometry *g1, *g2;
2116 int result;
2117 GBOX box1, box2;
2118
2119 geom1 = PG_GETARG_GSERIALIZED_P(0);
2120 geom2 = PG_GETARG_GSERIALIZED_P(1);
2121 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2122
2123 /* A.Crosses(Empty) == FALSE */
2124 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2125 PG_RETURN_BOOL(false);
2126
2127 /*
2128 * short-circuit 1: if geom2 bounding box does not overlap
2129 * geom1 bounding box we can return FALSE.
2130 */
2131 if ( gserialized_get_gbox_p(geom1, &box1) &&
2132 gserialized_get_gbox_p(geom2, &box2) )
2133 {
2134 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2135 {
2136 PG_RETURN_BOOL(false);
2137 }
2138 }
2139
2140 initGEOS(lwpgnotice, lwgeom_geos_error);
2141
2142 g1 = POSTGIS2GEOS(geom1);
2143 if (!g1)
2144 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2145
2146 g2 = POSTGIS2GEOS(geom2);
2147 if (!g2)
2148 {
2149 GEOSGeom_destroy(g1);
2150 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2151 }
2152
2153 result = GEOSCrosses(g1,g2);
2154
2155 GEOSGeom_destroy(g1);
2156 GEOSGeom_destroy(g2);
2157
2158 if (result == 2) HANDLE_GEOS_ERROR("GEOSCrosses");
2159
2160 PG_FREE_IF_COPY(geom1, 0);
2161 PG_FREE_IF_COPY(geom2, 1);
2162
2163 PG_RETURN_BOOL(result);
2164 }
2165
2166 PG_FUNCTION_INFO_V1(ST_Intersects);
ST_Intersects(PG_FUNCTION_ARGS)2167 Datum ST_Intersects(PG_FUNCTION_ARGS)
2168 {
2169 GSERIALIZED *geom1;
2170 GSERIALIZED *geom2;
2171 int result;
2172 GBOX box1, box2;
2173 PrepGeomCache *prep_cache;
2174
2175 geom1 = PG_GETARG_GSERIALIZED_P(0);
2176 geom2 = PG_GETARG_GSERIALIZED_P(1);
2177 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2178
2179 /* A.Intersects(Empty) == FALSE */
2180 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2181 PG_RETURN_BOOL(false);
2182
2183 /*
2184 * short-circuit 1: if geom2 bounding box does not overlap
2185 * geom1 bounding box we can return FALSE.
2186 */
2187 if ( gserialized_get_gbox_p(geom1, &box1) &&
2188 gserialized_get_gbox_p(geom2, &box2) )
2189 {
2190 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2191 PG_RETURN_BOOL(false);
2192 }
2193
2194 /*
2195 * short-circuit 2: if the geoms are a point and a polygon,
2196 * call the point_outside_polygon function.
2197 */
2198 if ((is_point(geom1) && is_poly(geom2)) || (is_poly(geom1) && is_point(geom2)))
2199 {
2200 GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
2201 GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
2202 RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
2203 int retval;
2204
2205 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2206 if (gserialized_get_type(gpoint) == POINTTYPE)
2207 {
2208 LWGEOM* point = lwgeom_from_gserialized(gpoint);
2209 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2210 lwgeom_free(point);
2211
2212 retval = (pip_result != -1); /* not outside */
2213 }
2214 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2215 {
2216 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2217 uint32_t i;
2218
2219 retval = LW_FALSE;
2220 for (i = 0; i < mpoint->ngeoms; i++)
2221 {
2222 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2223 if (pip_result != -1) /* not outside */
2224 {
2225 retval = LW_TRUE;
2226 break;
2227 }
2228 }
2229
2230 lwmpoint_free(mpoint);
2231 }
2232 else
2233 {
2234 /* Never get here */
2235 elog(ERROR,"Type isn't point or multipoint!");
2236 PG_RETURN_NULL();
2237 }
2238
2239 PG_FREE_IF_COPY(geom1, 0);
2240 PG_FREE_IF_COPY(geom2, 1);
2241 PG_RETURN_BOOL(retval);
2242 }
2243
2244 initGEOS(lwpgnotice, lwgeom_geos_error);
2245 prep_cache = GetPrepGeomCache( fcinfo, geom1, geom2 );
2246
2247 if ( prep_cache && prep_cache->prepared_geom )
2248 {
2249 if ( prep_cache->gcache.argnum == 1 )
2250 {
2251 GEOSGeometry *g = POSTGIS2GEOS(geom2);
2252 if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2253 result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2254 GEOSGeom_destroy(g);
2255 }
2256 else
2257 {
2258 GEOSGeometry *g = POSTGIS2GEOS(geom1);
2259 if (!g)
2260 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2261 result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2262 GEOSGeom_destroy(g);
2263 }
2264 }
2265 else
2266 {
2267 GEOSGeometry *g1;
2268 GEOSGeometry *g2;
2269 g1 = POSTGIS2GEOS(geom1);
2270 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2271 g2 = POSTGIS2GEOS(geom2);
2272 if (!g2)
2273 {
2274 GEOSGeom_destroy(g1);
2275 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2276 }
2277 result = GEOSIntersects( g1, g2);
2278 GEOSGeom_destroy(g1);
2279 GEOSGeom_destroy(g2);
2280 }
2281
2282 if (result == 2) HANDLE_GEOS_ERROR("GEOSIntersects");
2283
2284 PG_FREE_IF_COPY(geom1, 0);
2285 PG_FREE_IF_COPY(geom2, 1);
2286
2287 PG_RETURN_BOOL(result);
2288 }
2289
2290
2291 PG_FUNCTION_INFO_V1(touches);
touches(PG_FUNCTION_ARGS)2292 Datum touches(PG_FUNCTION_ARGS)
2293 {
2294 GSERIALIZED *geom1;
2295 GSERIALIZED *geom2;
2296 GEOSGeometry *g1, *g2;
2297 char result;
2298 GBOX box1, box2;
2299
2300 geom1 = PG_GETARG_GSERIALIZED_P(0);
2301 geom2 = PG_GETARG_GSERIALIZED_P(1);
2302 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2303
2304 /* A.Touches(Empty) == FALSE */
2305 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2306 PG_RETURN_BOOL(false);
2307
2308 /*
2309 * short-circuit 1: if geom2 bounding box does not overlap
2310 * geom1 bounding box we can return FALSE.
2311 */
2312 if ( gserialized_get_gbox_p(geom1, &box1) &&
2313 gserialized_get_gbox_p(geom2, &box2) )
2314 {
2315 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2316 {
2317 PG_RETURN_BOOL(false);
2318 }
2319 }
2320
2321 initGEOS(lwpgnotice, lwgeom_geos_error);
2322
2323 g1 = POSTGIS2GEOS(geom1 );
2324 if (!g1)
2325 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2326
2327 g2 = POSTGIS2GEOS(geom2 );
2328 if (!g2)
2329 {
2330 GEOSGeom_destroy(g1);
2331 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2332 }
2333
2334 result = GEOSTouches(g1,g2);
2335
2336 GEOSGeom_destroy(g1);
2337 GEOSGeom_destroy(g2);
2338
2339 if (result == 2) HANDLE_GEOS_ERROR("GEOSTouches");
2340
2341 PG_FREE_IF_COPY(geom1, 0);
2342 PG_FREE_IF_COPY(geom2, 1);
2343
2344 PG_RETURN_BOOL(result);
2345 }
2346
2347
2348 PG_FUNCTION_INFO_V1(disjoint);
disjoint(PG_FUNCTION_ARGS)2349 Datum disjoint(PG_FUNCTION_ARGS)
2350 {
2351 GSERIALIZED *geom1;
2352 GSERIALIZED *geom2;
2353 GEOSGeometry *g1, *g2;
2354 char result;
2355 GBOX box1, box2;
2356
2357 geom1 = PG_GETARG_GSERIALIZED_P(0);
2358 geom2 = PG_GETARG_GSERIALIZED_P(1);
2359 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2360
2361 /* A.Disjoint(Empty) == TRUE */
2362 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2363 PG_RETURN_BOOL(true);
2364
2365 /*
2366 * short-circuit 1: if geom2 bounding box does not overlap
2367 * geom1 bounding box we can return TRUE.
2368 */
2369 if ( gserialized_get_gbox_p(geom1, &box1) &&
2370 gserialized_get_gbox_p(geom2, &box2) )
2371 {
2372 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2373 {
2374 PG_RETURN_BOOL(true);
2375 }
2376 }
2377
2378 initGEOS(lwpgnotice, lwgeom_geos_error);
2379
2380 g1 = POSTGIS2GEOS(geom1);
2381 if (!g1)
2382 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2383
2384 g2 = POSTGIS2GEOS(geom2);
2385 if (!g2)
2386 {
2387 GEOSGeom_destroy(g1);
2388 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2389 }
2390
2391 result = GEOSDisjoint(g1,g2);
2392
2393 GEOSGeom_destroy(g1);
2394 GEOSGeom_destroy(g2);
2395
2396 if (result == 2) HANDLE_GEOS_ERROR("GEOSDisjoint");
2397
2398 PG_FREE_IF_COPY(geom1, 0);
2399 PG_FREE_IF_COPY(geom2, 1);
2400
2401 PG_RETURN_BOOL(result);
2402 }
2403
2404
2405 PG_FUNCTION_INFO_V1(relate_pattern);
relate_pattern(PG_FUNCTION_ARGS)2406 Datum relate_pattern(PG_FUNCTION_ARGS)
2407 {
2408 GSERIALIZED *geom1;
2409 GSERIALIZED *geom2;
2410 char *patt;
2411 char result;
2412 GEOSGeometry *g1, *g2;
2413 size_t i;
2414
2415 geom1 = PG_GETARG_GSERIALIZED_P(0);
2416 geom2 = PG_GETARG_GSERIALIZED_P(1);
2417 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2418
2419 /* TODO handle empty */
2420
2421 initGEOS(lwpgnotice, lwgeom_geos_error);
2422
2423 g1 = POSTGIS2GEOS(geom1);
2424 if (!g1)
2425 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2426 g2 = POSTGIS2GEOS(geom2);
2427 if (!g2)
2428 {
2429 GEOSGeom_destroy(g1);
2430 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2431 }
2432
2433 patt = DatumGetCString(DirectFunctionCall1(textout,
2434 PointerGetDatum(PG_GETARG_DATUM(2))));
2435
2436 /*
2437 ** Need to make sure 't' and 'f' are upper-case before handing to GEOS
2438 */
2439 for ( i = 0; i < strlen(patt); i++ )
2440 {
2441 if ( patt[i] == 't' ) patt[i] = 'T';
2442 if ( patt[i] == 'f' ) patt[i] = 'F';
2443 }
2444
2445 result = GEOSRelatePattern(g1,g2,patt);
2446 GEOSGeom_destroy(g1);
2447 GEOSGeom_destroy(g2);
2448 pfree(patt);
2449
2450 if (result == 2) HANDLE_GEOS_ERROR("GEOSRelatePattern");
2451
2452 PG_FREE_IF_COPY(geom1, 0);
2453 PG_FREE_IF_COPY(geom2, 1);
2454
2455 PG_RETURN_BOOL(result);
2456 }
2457
2458
2459
2460 PG_FUNCTION_INFO_V1(relate_full);
relate_full(PG_FUNCTION_ARGS)2461 Datum relate_full(PG_FUNCTION_ARGS)
2462 {
2463 GSERIALIZED *geom1;
2464 GSERIALIZED *geom2;
2465 GEOSGeometry *g1, *g2;
2466 char *relate_str;
2467 text *result;
2468 int bnr = GEOSRELATE_BNR_OGC;
2469
2470 /* TODO handle empty */
2471 geom1 = PG_GETARG_GSERIALIZED_P(0);
2472 geom2 = PG_GETARG_GSERIALIZED_P(1);
2473 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2474
2475 if ( PG_NARGS() > 2 )
2476 bnr = PG_GETARG_INT32(2);
2477
2478 initGEOS(lwpgnotice, lwgeom_geos_error);
2479
2480 g1 = POSTGIS2GEOS(geom1 );
2481 if (!g1)
2482 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2483 g2 = POSTGIS2GEOS(geom2 );
2484 if (!g2)
2485 {
2486 GEOSGeom_destroy(g1);
2487 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2488 }
2489
2490 POSTGIS_DEBUG(3, "constructed geometries ");
2491
2492 POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
2493 POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
2494
2495 relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
2496
2497 GEOSGeom_destroy(g1);
2498 GEOSGeom_destroy(g2);
2499
2500 if (!relate_str) HANDLE_GEOS_ERROR("GEOSRelate");
2501
2502 result = cstring_to_text(relate_str);
2503 GEOSFree(relate_str);
2504
2505 PG_FREE_IF_COPY(geom1, 0);
2506 PG_FREE_IF_COPY(geom2, 1);
2507
2508 PG_RETURN_TEXT_P(result);
2509 }
2510
2511
2512 PG_FUNCTION_INFO_V1(ST_Equals);
ST_Equals(PG_FUNCTION_ARGS)2513 Datum ST_Equals(PG_FUNCTION_ARGS)
2514 {
2515 GSERIALIZED *geom1;
2516 GSERIALIZED *geom2;
2517 GEOSGeometry *g1, *g2;
2518 char result;
2519 GBOX box1, box2;
2520
2521 geom1 = PG_GETARG_GSERIALIZED_P(0);
2522 geom2 = PG_GETARG_GSERIALIZED_P(1);
2523 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2524
2525 /* Empty == Empty */
2526 if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
2527 PG_RETURN_BOOL(true);
2528
2529 /*
2530 * short-circuit: If geom1 and geom2 do not have the same bounding box
2531 * we can return FALSE.
2532 */
2533 if ( gserialized_get_gbox_p(geom1, &box1) &&
2534 gserialized_get_gbox_p(geom2, &box2) )
2535 {
2536 if ( gbox_same_2d_float(&box1, &box2) == LW_FALSE )
2537 {
2538 PG_RETURN_BOOL(false);
2539 }
2540 }
2541
2542 /*
2543 * short-circuit: if geom1 and geom2 are binary-equivalent, we can return
2544 * TRUE. This is much faster than doing the comparison using GEOS.
2545 */
2546 if (VARSIZE(geom1) == VARSIZE(geom2) && !memcmp(geom1, geom2, VARSIZE(geom1))) {
2547 PG_RETURN_BOOL(true);
2548 }
2549
2550 initGEOS(lwpgnotice, lwgeom_geos_error);
2551
2552 g1 = POSTGIS2GEOS(geom1);
2553
2554 if (!g1)
2555 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2556
2557 g2 = POSTGIS2GEOS(geom2);
2558
2559 if (!g2)
2560 {
2561 GEOSGeom_destroy(g1);
2562 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2563 }
2564
2565 result = GEOSEquals(g1,g2);
2566
2567 GEOSGeom_destroy(g1);
2568 GEOSGeom_destroy(g2);
2569
2570 if (result == 2) HANDLE_GEOS_ERROR("GEOSEquals");
2571
2572 PG_FREE_IF_COPY(geom1, 0);
2573 PG_FREE_IF_COPY(geom2, 1);
2574
2575 PG_RETURN_BOOL(result);
2576 }
2577
2578 PG_FUNCTION_INFO_V1(issimple);
issimple(PG_FUNCTION_ARGS)2579 Datum issimple(PG_FUNCTION_ARGS)
2580 {
2581 GSERIALIZED *geom;
2582 LWGEOM *lwgeom_in;
2583 int result;
2584
2585 POSTGIS_DEBUG(2, "issimple called");
2586
2587 geom = PG_GETARG_GSERIALIZED_P(0);
2588
2589 if ( gserialized_is_empty(geom) )
2590 PG_RETURN_BOOL(true);
2591
2592 lwgeom_in = lwgeom_from_gserialized(geom);
2593 result = lwgeom_is_simple(lwgeom_in);
2594 lwgeom_free(lwgeom_in) ;
2595 PG_FREE_IF_COPY(geom, 0);
2596
2597 if (result == -1) {
2598 PG_RETURN_NULL(); /*never get here */
2599 }
2600
2601 PG_RETURN_BOOL(result);
2602 }
2603
2604 PG_FUNCTION_INFO_V1(isring);
isring(PG_FUNCTION_ARGS)2605 Datum isring(PG_FUNCTION_ARGS)
2606 {
2607 GSERIALIZED *geom;
2608 GEOSGeometry *g1;
2609 int result;
2610
2611 geom = PG_GETARG_GSERIALIZED_P(0);
2612
2613 /* Empty things can't close */
2614 if ( gserialized_is_empty(geom) )
2615 PG_RETURN_BOOL(false);
2616
2617 initGEOS(lwpgnotice, lwgeom_geos_error);
2618
2619 g1 = POSTGIS2GEOS(geom);
2620 if (!g1)
2621 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2622
2623 if ( GEOSGeomTypeId(g1) != GEOS_LINESTRING )
2624 {
2625 GEOSGeom_destroy(g1);
2626 elog(ERROR, "ST_IsRing() should only be called on a linear feature");
2627 }
2628
2629 result = GEOSisRing(g1);
2630 GEOSGeom_destroy(g1);
2631
2632 if (result == 2) HANDLE_GEOS_ERROR("GEOSisRing");
2633
2634 PG_FREE_IF_COPY(geom, 0);
2635 PG_RETURN_BOOL(result);
2636 }
2637
2638 GSERIALIZED *
GEOS2POSTGIS(GEOSGeom geom,char want3d)2639 GEOS2POSTGIS(GEOSGeom geom, char want3d)
2640 {
2641 LWGEOM *lwgeom;
2642 GSERIALIZED *result;
2643
2644 lwgeom = GEOS2LWGEOM(geom, want3d);
2645 if ( ! lwgeom )
2646 {
2647 lwpgerror("%s: GEOS2LWGEOM returned NULL", __func__);
2648 return NULL;
2649 }
2650
2651 POSTGIS_DEBUGF(4, "%s: GEOS2LWGEOM returned a %s", __func__, lwgeom_summary(lwgeom, 0));
2652
2653 if (lwgeom_needs_bbox(lwgeom)) lwgeom_add_bbox(lwgeom);
2654
2655 result = geometry_serialize(lwgeom);
2656 lwgeom_free(lwgeom);
2657
2658 return result;
2659 }
2660
2661 /*-----=POSTGIS2GEOS= */
2662
2663
2664 GEOSGeometry *
POSTGIS2GEOS(GSERIALIZED * pglwgeom)2665 POSTGIS2GEOS(GSERIALIZED *pglwgeom)
2666 {
2667 GEOSGeometry *ret;
2668 LWGEOM *lwgeom = lwgeom_from_gserialized(pglwgeom);
2669 if ( ! lwgeom )
2670 {
2671 lwpgerror("POSTGIS2GEOS: unable to deserialize input");
2672 return NULL;
2673 }
2674 ret = LWGEOM2GEOS(lwgeom, 0);
2675 lwgeom_free(lwgeom);
2676
2677 return ret;
2678 }
2679
array_nelems_not_null(ArrayType * array)2680 uint32_t array_nelems_not_null(ArrayType* array) {
2681 ArrayIterator iterator;
2682 Datum value;
2683 bool isnull;
2684 uint32_t nelems_not_null = 0;
2685 iterator = array_create_iterator(array, 0, NULL);
2686 while(array_iterate(iterator, &value, &isnull) )
2687 if (!isnull)
2688 nelems_not_null++;
2689
2690 array_free_iterator(iterator);
2691
2692 return nelems_not_null;
2693 }
2694
2695 /* ARRAY2LWGEOM: Converts the non-null elements of a Postgres array into a LWGEOM* array */
ARRAY2LWGEOM(ArrayType * array,uint32_t nelems,int * is3d,int * srid)2696 LWGEOM** ARRAY2LWGEOM(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2697 {
2698 ArrayIterator iterator;
2699 Datum value;
2700 bool isnull;
2701 bool gotsrid = false;
2702 uint32_t i = 0;
2703
2704 LWGEOM** lw_geoms = palloc(nelems * sizeof(LWGEOM*));
2705
2706 iterator = array_create_iterator(array, 0, NULL);
2707
2708 while (array_iterate(iterator, &value, &isnull))
2709 {
2710 GSERIALIZED *geom = (GSERIALIZED *)DatumGetPointer(value);
2711
2712 if (isnull)
2713 continue;
2714
2715 *is3d = *is3d || gserialized_has_z(geom);
2716
2717 lw_geoms[i] = lwgeom_from_gserialized(geom);
2718 if (!lw_geoms[i]) /* error in creation */
2719 {
2720 lwpgerror("Geometry deserializing geometry");
2721 return NULL;
2722 }
2723 if (!gotsrid)
2724 {
2725 gotsrid = true;
2726 *srid = gserialized_get_srid(geom);
2727 }
2728 else
2729 gserialized_error_if_srid_mismatch_reference(geom, *srid, __func__);
2730
2731 i++;
2732 }
2733
2734 return lw_geoms;
2735 }
2736
2737 /* ARRAY2GEOS: Converts the non-null elements of a Postgres array into a GEOSGeometry* array */
ARRAY2GEOS(ArrayType * array,uint32_t nelems,int * is3d,int * srid)2738 GEOSGeometry** ARRAY2GEOS(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2739 {
2740 ArrayIterator iterator;
2741 Datum value;
2742 bool isnull;
2743 bool gotsrid = false;
2744 uint32_t i = 0;
2745
2746 GEOSGeometry** geos_geoms = palloc(nelems * sizeof(GEOSGeometry*));
2747
2748 iterator = array_create_iterator(array, 0, NULL);
2749
2750 while(array_iterate(iterator, &value, &isnull))
2751 {
2752 GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2753
2754 if (isnull)
2755 continue;
2756
2757 *is3d = *is3d || gserialized_has_z(geom);
2758
2759 geos_geoms[i] = POSTGIS2GEOS(geom);
2760 if (!geos_geoms[i])
2761 {
2762 uint32_t j;
2763 lwpgerror("Geometry could not be converted to GEOS");
2764
2765 for (j = 0; j < i; j++) {
2766 GEOSGeom_destroy(geos_geoms[j]);
2767 }
2768 return NULL;
2769 }
2770
2771 if (!gotsrid)
2772 {
2773 *srid = gserialized_get_srid(geom);
2774 gotsrid = true;
2775 }
2776 else if (*srid != gserialized_get_srid(geom))
2777 {
2778 uint32_t j;
2779 for (j = 0; j <= i; j++) {
2780 GEOSGeom_destroy(geos_geoms[j]);
2781 }
2782 gserialized_error_if_srid_mismatch_reference(geom, *srid, __func__);
2783 return NULL;
2784 }
2785
2786 i++;
2787 }
2788
2789 array_free_iterator(iterator);
2790 return geos_geoms;
2791 }
2792
2793 PG_FUNCTION_INFO_V1(GEOSnoop);
GEOSnoop(PG_FUNCTION_ARGS)2794 Datum GEOSnoop(PG_FUNCTION_ARGS)
2795 {
2796 GSERIALIZED *geom;
2797 GEOSGeometry *geosgeom;
2798 GSERIALIZED *lwgeom_result;
2799
2800 initGEOS(lwpgnotice, lwgeom_geos_error);
2801
2802 geom = PG_GETARG_GSERIALIZED_P(0);
2803 geosgeom = POSTGIS2GEOS(geom);
2804 if ( ! geosgeom ) PG_RETURN_NULL();
2805
2806 lwgeom_result = GEOS2POSTGIS(geosgeom, gserialized_has_z(geom));
2807 GEOSGeom_destroy(geosgeom);
2808
2809 PG_FREE_IF_COPY(geom, 0);
2810
2811 PG_RETURN_POINTER(lwgeom_result);
2812 }
2813
2814 PG_FUNCTION_INFO_V1(polygonize_garray);
polygonize_garray(PG_FUNCTION_ARGS)2815 Datum polygonize_garray(PG_FUNCTION_ARGS)
2816 {
2817 ArrayType *array;
2818 int is3d = 0;
2819 uint32 nelems, i;
2820 GSERIALIZED *result;
2821 GEOSGeometry *geos_result;
2822 const GEOSGeometry **vgeoms;
2823 int32_t srid = SRID_UNKNOWN;
2824 #if POSTGIS_DEBUG_LEVEL >= 3
2825 static int call=1;
2826 #endif
2827
2828 #if POSTGIS_DEBUG_LEVEL >= 3
2829 call++;
2830 #endif
2831
2832 if (PG_ARGISNULL(0))
2833 PG_RETURN_NULL();
2834
2835 array = PG_GETARG_ARRAYTYPE_P(0);
2836 nelems = array_nelems_not_null(array);
2837
2838 if (nelems == 0)
2839 PG_RETURN_NULL();
2840
2841 POSTGIS_DEBUGF(3, "polygonize_garray: number of non-null elements: %d", nelems);
2842
2843 /* Ok, we really need geos now ;) */
2844 initGEOS(lwpgnotice, lwgeom_geos_error);
2845
2846 vgeoms = (const GEOSGeometry**) ARRAY2GEOS(array, nelems, &is3d, &srid);
2847
2848 POSTGIS_DEBUG(3, "polygonize_garray: invoking GEOSpolygonize");
2849
2850 geos_result = GEOSPolygonize(vgeoms, nelems);
2851
2852 POSTGIS_DEBUG(3, "polygonize_garray: GEOSpolygonize returned");
2853
2854 for (i=0; i<nelems; ++i) GEOSGeom_destroy((GEOSGeometry *)vgeoms[i]);
2855 pfree(vgeoms);
2856
2857 if ( ! geos_result ) PG_RETURN_NULL();
2858
2859 GEOSSetSRID(geos_result, srid);
2860 result = GEOS2POSTGIS(geos_result, is3d);
2861 GEOSGeom_destroy(geos_result);
2862 if (!result)
2863 {
2864 elog(ERROR, "%s returned an error", __func__);
2865 PG_RETURN_NULL(); /*never get here */
2866 }
2867
2868 PG_RETURN_POINTER(result);
2869 }
2870
2871
2872 PG_FUNCTION_INFO_V1(clusterintersecting_garray);
clusterintersecting_garray(PG_FUNCTION_ARGS)2873 Datum clusterintersecting_garray(PG_FUNCTION_ARGS)
2874 {
2875 Datum* result_array_data;
2876 ArrayType *array, *result;
2877 int is3d = 0;
2878 uint32 nelems, nclusters, i;
2879 GEOSGeometry **geos_inputs, **geos_results;
2880 int32_t srid = SRID_UNKNOWN;
2881
2882 /* Parameters used to construct a result array */
2883 int16 elmlen;
2884 bool elmbyval;
2885 char elmalign;
2886
2887 /* Null array, null geometry (should be empty?) */
2888 if (PG_ARGISNULL(0))
2889 PG_RETURN_NULL();
2890
2891 array = PG_GETARG_ARRAYTYPE_P(0);
2892 nelems = array_nelems_not_null(array);
2893
2894 POSTGIS_DEBUGF(3, "clusterintersecting_garray: number of non-null elements: %d", nelems);
2895
2896 if ( nelems == 0 ) PG_RETURN_NULL();
2897
2898 /* TODO short-circuit for one element? */
2899
2900 /* Ok, we really need geos now ;) */
2901 initGEOS(lwpgnotice, lwgeom_geos_error);
2902
2903 geos_inputs = ARRAY2GEOS(array, nelems, &is3d, &srid);
2904 if(!geos_inputs)
2905 {
2906 PG_RETURN_NULL();
2907 }
2908
2909 if (cluster_intersecting(geos_inputs, nelems, &geos_results, &nclusters) != LW_SUCCESS)
2910 {
2911 elog(ERROR, "clusterintersecting: Error performing clustering");
2912 PG_RETURN_NULL();
2913 }
2914 pfree(geos_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
2915
2916 if (!geos_results) PG_RETURN_NULL();
2917
2918 result_array_data = palloc(nclusters * sizeof(Datum));
2919 for (i=0; i<nclusters; ++i)
2920 {
2921 result_array_data[i] = PointerGetDatum(GEOS2POSTGIS(geos_results[i], is3d));
2922 GEOSGeom_destroy(geos_results[i]);
2923 }
2924 lwfree(geos_results);
2925
2926 get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
2927 result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
2928
2929 if (!result)
2930 {
2931 elog(ERROR, "clusterintersecting: Error constructing return-array");
2932 PG_RETURN_NULL();
2933 }
2934
2935 PG_RETURN_POINTER(result);
2936 }
2937
2938 PG_FUNCTION_INFO_V1(cluster_within_distance_garray);
cluster_within_distance_garray(PG_FUNCTION_ARGS)2939 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS)
2940 {
2941 Datum* result_array_data;
2942 ArrayType *array, *result;
2943 int is3d = 0;
2944 uint32 nelems, nclusters, i;
2945 LWGEOM** lw_inputs;
2946 LWGEOM** lw_results;
2947 double tolerance;
2948 int32_t srid = SRID_UNKNOWN;
2949
2950 /* Parameters used to construct a result array */
2951 int16 elmlen;
2952 bool elmbyval;
2953 char elmalign;
2954
2955 /* Null array, null geometry (should be empty?) */
2956 if (PG_ARGISNULL(0))
2957 PG_RETURN_NULL();
2958
2959 array = PG_GETARG_ARRAYTYPE_P(0);
2960
2961 tolerance = PG_GETARG_FLOAT8(1);
2962 if (tolerance < 0)
2963 {
2964 lwpgerror("Tolerance must be a positive number.");
2965 PG_RETURN_NULL();
2966 }
2967
2968 nelems = array_nelems_not_null(array);
2969
2970 POSTGIS_DEBUGF(3, "cluster_within_distance_garray: number of non-null elements: %d", nelems);
2971
2972 if ( nelems == 0 ) PG_RETURN_NULL();
2973
2974 /* TODO short-circuit for one element? */
2975
2976 /* Ok, we really need geos now ;) */
2977 initGEOS(lwpgnotice, lwgeom_geos_error);
2978
2979 lw_inputs = ARRAY2LWGEOM(array, nelems, &is3d, &srid);
2980 if (!lw_inputs)
2981 {
2982 PG_RETURN_NULL();
2983 }
2984
2985 if (cluster_within_distance(lw_inputs, nelems, tolerance, &lw_results, &nclusters) != LW_SUCCESS)
2986 {
2987 elog(ERROR, "cluster_within: Error performing clustering");
2988 PG_RETURN_NULL();
2989 }
2990 pfree(lw_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
2991
2992 if (!lw_results) PG_RETURN_NULL();
2993
2994 result_array_data = palloc(nclusters * sizeof(Datum));
2995 for (i=0; i<nclusters; ++i)
2996 {
2997 result_array_data[i] = PointerGetDatum(geometry_serialize(lw_results[i]));
2998 lwgeom_free(lw_results[i]);
2999 }
3000 lwfree(lw_results);
3001
3002 get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3003 result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3004
3005 if (!result)
3006 {
3007 elog(ERROR, "clusterwithin: Error constructing return-array");
3008 PG_RETURN_NULL();
3009 }
3010
3011 PG_RETURN_POINTER(result);
3012 }
3013
3014 PG_FUNCTION_INFO_V1(linemerge);
linemerge(PG_FUNCTION_ARGS)3015 Datum linemerge(PG_FUNCTION_ARGS)
3016 {
3017 GSERIALIZED *geom1;
3018 GSERIALIZED *result;
3019 LWGEOM *lwgeom1, *lwresult ;
3020
3021 geom1 = PG_GETARG_GSERIALIZED_P(0);
3022
3023
3024 lwgeom1 = lwgeom_from_gserialized(geom1) ;
3025
3026 lwresult = lwgeom_linemerge(lwgeom1);
3027 result = geometry_serialize(lwresult) ;
3028
3029 lwgeom_free(lwgeom1) ;
3030 lwgeom_free(lwresult) ;
3031
3032 PG_FREE_IF_COPY(geom1, 0);
3033
3034 PG_RETURN_POINTER(result);
3035 }
3036
3037 /*
3038 * Take a geometry and return an areal geometry
3039 * (Polygon or MultiPolygon).
3040 * Actually a wrapper around GEOSpolygonize,
3041 * transforming the resulting collection into
3042 * a valid polygon Geometry.
3043 */
3044 PG_FUNCTION_INFO_V1(ST_BuildArea);
ST_BuildArea(PG_FUNCTION_ARGS)3045 Datum ST_BuildArea(PG_FUNCTION_ARGS)
3046 {
3047 GSERIALIZED *result;
3048 GSERIALIZED *geom;
3049 LWGEOM *lwgeom_in, *lwgeom_out;
3050
3051 geom = PG_GETARG_GSERIALIZED_P(0);
3052 lwgeom_in = lwgeom_from_gserialized(geom);
3053
3054 lwgeom_out = lwgeom_buildarea(lwgeom_in);
3055 lwgeom_free(lwgeom_in) ;
3056
3057 if ( ! lwgeom_out )
3058 {
3059 PG_FREE_IF_COPY(geom, 0);
3060 PG_RETURN_NULL();
3061 }
3062
3063 result = geometry_serialize(lwgeom_out) ;
3064 lwgeom_free(lwgeom_out) ;
3065
3066 PG_FREE_IF_COPY(geom, 0);
3067 PG_RETURN_POINTER(result);
3068 }
3069
3070 /*
3071 * Take the vertices of a geometry and builds
3072 * Delaunay triangles around them.
3073 */
3074 PG_FUNCTION_INFO_V1(ST_DelaunayTriangles);
ST_DelaunayTriangles(PG_FUNCTION_ARGS)3075 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS)
3076 {
3077 GSERIALIZED *result;
3078 GSERIALIZED *geom;
3079 LWGEOM *lwgeom_in, *lwgeom_out;
3080 double tolerance = 0.0;
3081 int flags = 0;
3082
3083 geom = PG_GETARG_GSERIALIZED_P(0);
3084 tolerance = PG_GETARG_FLOAT8(1);
3085 flags = PG_GETARG_INT32(2);
3086
3087 lwgeom_in = lwgeom_from_gserialized(geom);
3088 lwgeom_out = lwgeom_delaunay_triangulation(lwgeom_in, tolerance, flags);
3089 lwgeom_free(lwgeom_in) ;
3090
3091 if ( ! lwgeom_out )
3092 {
3093 PG_FREE_IF_COPY(geom, 0);
3094 PG_RETURN_NULL();
3095 }
3096
3097 result = geometry_serialize(lwgeom_out) ;
3098 lwgeom_free(lwgeom_out) ;
3099
3100 PG_FREE_IF_COPY(geom, 0);
3101 PG_RETURN_POINTER(result);
3102 }
3103
3104 /*
3105 * ST_Snap
3106 *
3107 * Snap a geometry to another with a given tolerance
3108 */
3109 Datum ST_Snap(PG_FUNCTION_ARGS);
3110 PG_FUNCTION_INFO_V1(ST_Snap);
ST_Snap(PG_FUNCTION_ARGS)3111 Datum ST_Snap(PG_FUNCTION_ARGS)
3112 {
3113 GSERIALIZED *geom1, *geom2, *result;
3114 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
3115 double tolerance;
3116
3117 geom1 = PG_GETARG_GSERIALIZED_P(0);
3118 geom2 = PG_GETARG_GSERIALIZED_P(1);
3119 tolerance = PG_GETARG_FLOAT8(2);
3120
3121 lwgeom1 = lwgeom_from_gserialized(geom1);
3122 lwgeom2 = lwgeom_from_gserialized(geom2);
3123
3124 lwresult = lwgeom_snap(lwgeom1, lwgeom2, tolerance);
3125 lwgeom_free(lwgeom1);
3126 lwgeom_free(lwgeom2);
3127 PG_FREE_IF_COPY(geom1, 0);
3128 PG_FREE_IF_COPY(geom2, 1);
3129
3130 result = geometry_serialize(lwresult);
3131 lwgeom_free(lwresult);
3132
3133 PG_RETURN_POINTER(result);
3134 }
3135
3136 /*
3137 * ST_Split
3138 *
3139 * Split polygon by line, line by line, line by point.
3140 * Returns at most components as a collection.
3141 * First element of the collection is always the part which
3142 * remains after the cut, while the second element is the
3143 * part which has been cut out. We arbitrarily take the part
3144 * on the *right* of cut lines as the part which has been cut out.
3145 * For a line cut by a point the part which remains is the one
3146 * from start of the line to the cut point.
3147 *
3148 *
3149 * Author: Sandro Santilli <strk@kbt.io>
3150 *
3151 * Work done for Faunalia (http://www.faunalia.it) with fundings
3152 * from Regione Toscana - Sistema Informativo per il Governo
3153 * del Territorio e dell'Ambiente (RT-SIGTA).
3154 *
3155 * Thanks to the PostGIS community for sharing poly/line ideas [1]
3156 *
3157 * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
3158 *
3159 */
3160 Datum ST_Split(PG_FUNCTION_ARGS);
3161 PG_FUNCTION_INFO_V1(ST_Split);
ST_Split(PG_FUNCTION_ARGS)3162 Datum ST_Split(PG_FUNCTION_ARGS)
3163 {
3164 GSERIALIZED *in, *blade_in, *out;
3165 LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out;
3166
3167 in = PG_GETARG_GSERIALIZED_P(0);
3168 blade_in = PG_GETARG_GSERIALIZED_P(1);
3169 gserialized_error_if_srid_mismatch(in, blade_in, __func__);
3170
3171 lwgeom_in = lwgeom_from_gserialized(in);
3172 lwblade_in = lwgeom_from_gserialized(blade_in);
3173
3174 lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
3175 lwgeom_free(lwgeom_in);
3176 lwgeom_free(lwblade_in);
3177
3178 if ( ! lwgeom_out )
3179 {
3180 PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3181 PG_FREE_IF_COPY(blade_in, 1);
3182 PG_RETURN_NULL();
3183 }
3184
3185 out = geometry_serialize(lwgeom_out);
3186 lwgeom_free(lwgeom_out);
3187 PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3188 PG_FREE_IF_COPY(blade_in, 1);
3189
3190 PG_RETURN_POINTER(out);
3191 }
3192
3193 /**********************************************************************
3194 *
3195 * ST_SharedPaths
3196 *
3197 * Return the set of paths shared between two linear geometries,
3198 * and their direction (same or opposite).
3199 *
3200 * Developed by Sandro Santilli (strk@kbt.io) for Faunalia
3201 * (http://www.faunalia.it) with funding from Regione Toscana - Sistema
3202 * Informativo per la Gestione del Territorio e dell' Ambiente
3203 * [RT-SIGTA]". For the project: "Sviluppo strumenti software per il
3204 * trattamento di dati geografici basati su QuantumGIS e Postgis (CIG
3205 * 0494241492)"
3206 *
3207 **********************************************************************/
3208 Datum ST_SharedPaths(PG_FUNCTION_ARGS);
3209 PG_FUNCTION_INFO_V1(ST_SharedPaths);
ST_SharedPaths(PG_FUNCTION_ARGS)3210 Datum ST_SharedPaths(PG_FUNCTION_ARGS)
3211 {
3212 GSERIALIZED *geom1, *geom2, *out;
3213 LWGEOM *g1, *g2, *lwgeom_out;
3214
3215 geom1 = PG_GETARG_GSERIALIZED_P(0);
3216 geom2 = PG_GETARG_GSERIALIZED_P(1);
3217
3218 g1 = lwgeom_from_gserialized(geom1);
3219 g2 = lwgeom_from_gserialized(geom2);
3220
3221 lwgeom_out = lwgeom_sharedpaths(g1, g2);
3222 lwgeom_free(g1);
3223 lwgeom_free(g2);
3224
3225 if ( ! lwgeom_out )
3226 {
3227 PG_FREE_IF_COPY(geom1, 0);
3228 PG_FREE_IF_COPY(geom2, 1);
3229 PG_RETURN_NULL();
3230 }
3231
3232 out = geometry_serialize(lwgeom_out);
3233 lwgeom_free(lwgeom_out);
3234
3235 PG_FREE_IF_COPY(geom1, 0);
3236 PG_FREE_IF_COPY(geom2, 1);
3237 PG_RETURN_POINTER(out);
3238 }
3239
3240
3241 /**********************************************************************
3242 *
3243 * ST_Node
3244 *
3245 * Fully node a set of lines using the least possible nodes while
3246 * preserving all of the input ones.
3247 *
3248 **********************************************************************/
3249 Datum ST_Node(PG_FUNCTION_ARGS);
3250 PG_FUNCTION_INFO_V1(ST_Node);
ST_Node(PG_FUNCTION_ARGS)3251 Datum ST_Node(PG_FUNCTION_ARGS)
3252 {
3253 GSERIALIZED *geom1, *out;
3254 LWGEOM *g1, *lwgeom_out;
3255
3256 geom1 = PG_GETARG_GSERIALIZED_P(0);
3257
3258 g1 = lwgeom_from_gserialized(geom1);
3259
3260 lwgeom_out = lwgeom_node(g1);
3261 lwgeom_free(g1);
3262
3263 if ( ! lwgeom_out )
3264 {
3265 PG_FREE_IF_COPY(geom1, 0);
3266 PG_RETURN_NULL();
3267 }
3268
3269 out = geometry_serialize(lwgeom_out);
3270 lwgeom_free(lwgeom_out);
3271
3272 PG_FREE_IF_COPY(geom1, 0);
3273 PG_RETURN_POINTER(out);
3274 }
3275
3276 /******************************************
3277 *
3278 * ST_Voronoi
3279 *
3280 * Returns a Voronoi diagram constructed
3281 * from the points of the input geometry.
3282 *
3283 ******************************************/
3284 Datum ST_Voronoi(PG_FUNCTION_ARGS);
3285 PG_FUNCTION_INFO_V1(ST_Voronoi);
ST_Voronoi(PG_FUNCTION_ARGS)3286 Datum ST_Voronoi(PG_FUNCTION_ARGS)
3287 {
3288 GSERIALIZED* input;
3289 GSERIALIZED* clip;
3290 GSERIALIZED* result;
3291 LWGEOM* lwgeom_input;
3292 LWGEOM* lwgeom_result;
3293 double tolerance;
3294 GBOX clip_envelope;
3295 int custom_clip_envelope;
3296 int return_polygons;
3297
3298 /* Return NULL on NULL geometry */
3299 if (PG_ARGISNULL(0))
3300 PG_RETURN_NULL();
3301
3302 /* Read our tolerance value */
3303 if (PG_ARGISNULL(2))
3304 {
3305 lwpgerror("Tolerance must be a positive number.");
3306 PG_RETURN_NULL();
3307 }
3308
3309 tolerance = PG_GETARG_FLOAT8(2);
3310
3311 if (tolerance < 0)
3312 {
3313 lwpgerror("Tolerance must be a positive number.");
3314 PG_RETURN_NULL();
3315 }
3316
3317 /* Are we returning lines or polygons? */
3318 if (PG_ARGISNULL(3))
3319 {
3320 lwpgerror("return_polygons must be true or false.");
3321 PG_RETURN_NULL();
3322 }
3323 return_polygons = PG_GETARG_BOOL(3);
3324
3325 /* Read our clipping envelope, if applicable. */
3326 custom_clip_envelope = !PG_ARGISNULL(1);
3327 if (custom_clip_envelope) {
3328 clip = PG_GETARG_GSERIALIZED_P(1);
3329 if (!gserialized_get_gbox_p(clip, &clip_envelope))
3330 {
3331 lwpgerror("Could not determine envelope of clipping geometry.");
3332 PG_FREE_IF_COPY(clip, 1);
3333 PG_RETURN_NULL();
3334 }
3335 PG_FREE_IF_COPY(clip, 1);
3336 }
3337
3338 /* Read our input geometry */
3339 input = PG_GETARG_GSERIALIZED_P(0);
3340
3341 lwgeom_input = lwgeom_from_gserialized(input);
3342
3343 if(!lwgeom_input)
3344 {
3345 lwpgerror("Could not read input geometry.");
3346 PG_FREE_IF_COPY(input, 0);
3347 PG_RETURN_NULL();
3348 }
3349
3350 lwgeom_result = lwgeom_voronoi_diagram(lwgeom_input, custom_clip_envelope ? &clip_envelope : NULL, tolerance, !return_polygons);
3351 lwgeom_free(lwgeom_input);
3352
3353 if (!lwgeom_result)
3354 {
3355 lwpgerror("Error computing Voronoi diagram.");
3356 PG_FREE_IF_COPY(input, 0);
3357 PG_RETURN_NULL();
3358 }
3359
3360 result = geometry_serialize(lwgeom_result);
3361 lwgeom_free(lwgeom_result);
3362
3363 PG_FREE_IF_COPY(input, 0);
3364 PG_RETURN_POINTER(result);
3365 }
3366
3367 /******************************************
3368 *
3369 * ST_MinimumClearance
3370 *
3371 * Returns the minimum clearance of a geometry.
3372 *
3373 ******************************************/
3374 Datum ST_MinimumClearance(PG_FUNCTION_ARGS);
3375 PG_FUNCTION_INFO_V1(ST_MinimumClearance);
ST_MinimumClearance(PG_FUNCTION_ARGS)3376 Datum ST_MinimumClearance(PG_FUNCTION_ARGS)
3377 {
3378 GSERIALIZED* input;
3379 GEOSGeometry* input_geos;
3380 int error;
3381 double result;
3382
3383 initGEOS(lwpgnotice, lwgeom_geos_error);
3384
3385 input = PG_GETARG_GSERIALIZED_P(0);
3386 input_geos = POSTGIS2GEOS(input);
3387 if (!input_geos)
3388 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3389
3390 error = GEOSMinimumClearance(input_geos, &result);
3391 GEOSGeom_destroy(input_geos);
3392 if (error) HANDLE_GEOS_ERROR("Error computing minimum clearance");
3393
3394 PG_FREE_IF_COPY(input, 0);
3395 PG_RETURN_FLOAT8(result);
3396 }
3397
3398 /******************************************
3399 *
3400 * ST_MinimumClearanceLine
3401 *
3402 * Returns the minimum clearance line of a geometry.
3403 *
3404 ******************************************/
3405 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS);
3406 PG_FUNCTION_INFO_V1(ST_MinimumClearanceLine);
ST_MinimumClearanceLine(PG_FUNCTION_ARGS)3407 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS)
3408 {
3409 GSERIALIZED* input;
3410 GSERIALIZED* result;
3411 GEOSGeometry* input_geos;
3412 GEOSGeometry* result_geos;
3413 int32_t srid;
3414
3415 initGEOS(lwpgnotice, lwgeom_geos_error);
3416
3417 input = PG_GETARG_GSERIALIZED_P(0);
3418 srid = gserialized_get_srid(input);
3419 input_geos = POSTGIS2GEOS(input);
3420 if (!input_geos)
3421 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3422
3423 result_geos = GEOSMinimumClearanceLine(input_geos);
3424 GEOSGeom_destroy(input_geos);
3425 if (!result_geos)
3426 HANDLE_GEOS_ERROR("Error computing minimum clearance");
3427
3428 GEOSSetSRID(result_geos, srid);
3429 result = GEOS2POSTGIS(result_geos, LW_FALSE);
3430 GEOSGeom_destroy(result_geos);
3431
3432 PG_FREE_IF_COPY(input, 0);
3433 PG_RETURN_POINTER(result);
3434 }
3435
3436 /******************************************
3437 *
3438 * ST_OrientedEnvelope
3439 *
3440 ******************************************/
3441 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS);
3442 PG_FUNCTION_INFO_V1(ST_OrientedEnvelope);
ST_OrientedEnvelope(PG_FUNCTION_ARGS)3443 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS)
3444 {
3445 GSERIALIZED* input;
3446 GSERIALIZED* result;
3447 GEOSGeometry* input_geos;
3448 GEOSGeometry* result_geos;
3449 int32_t srid;
3450
3451 initGEOS(lwpgnotice, lwgeom_geos_error);
3452
3453 input = PG_GETARG_GSERIALIZED_P(0);
3454 srid = gserialized_get_srid(input);
3455 input_geos = POSTGIS2GEOS(input);
3456 if (!input_geos)
3457 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3458
3459 result_geos = GEOSMinimumRotatedRectangle(input_geos);
3460 GEOSGeom_destroy(input_geos);
3461 if (!result_geos)
3462 HANDLE_GEOS_ERROR("Error computing oriented envelope");
3463
3464 GEOSSetSRID(result_geos, srid);
3465 result = GEOS2POSTGIS(result_geos, LW_FALSE);
3466 GEOSGeom_destroy(result_geos);
3467
3468 PG_FREE_IF_COPY(input, 0);
3469 PG_RETURN_POINTER(result);
3470 }
3471