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