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