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