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