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