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