PostGIS  2.5.0dev-r@@SVN_REVISION@@
liblwgeom/lwgeom_geos.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright 2011-2014 Sandro Santilli <strk@kbt.io>
22  * Copyright 2017-2018 Darafei Praliaskouski <me@komzpa.net>
23  *
24  **********************************************************************/
25 
26 #include "lwgeom_geos.h"
27 #include "liblwgeom.h"
28 #include "liblwgeom_internal.h"
29 #include "lwgeom_log.h"
30 
31 #include <stdlib.h>
32 #include <time.h>
33 
34 LWTIN* lwtin_from_geos(const GEOSGeometry* geom, uint8_t want3d);
35 
36 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
38 
39 extern void
40 lwgeom_geos_error(const char* fmt, ...)
41 {
42  va_list ap;
43  va_start(ap, fmt);
44 
45  /* Call the supplied function */
48 
49  va_end(ap);
50 }
51 
52 /*
53 ** GEOS <==> PostGIS conversion functions
54 **
55 ** Default conversion creates a GEOS point array, then iterates through the
56 ** PostGIS points, setting each value in the GEOS array one at a time.
57 **
58 */
59 
60 /* Return a POINTARRAY from a GEOSCoordSeq */
62 ptarray_from_GEOSCoordSeq(const GEOSCoordSequence* cs, uint8_t want3d)
63 {
64  uint32_t dims = 2;
65  uint32_t size, i;
66  POINTARRAY* pa;
67  POINT4D point;
68 
69  LWDEBUG(2, "ptarray_fromGEOSCoordSeq called");
70 
71  if (!GEOSCoordSeq_getSize(cs, &size)) lwerror("Exception thrown");
72 
73  LWDEBUGF(4, " GEOSCoordSeq size: %d", size);
74 
75  if (want3d)
76  {
77  if (!GEOSCoordSeq_getDimensions(cs, &dims)) lwerror("Exception thrown");
78 
79  LWDEBUGF(4, " GEOSCoordSeq dimensions: %d", dims);
80 
81  /* forget higher dimensions (if any) */
82  if (dims > 3) dims = 3;
83  }
84 
85  LWDEBUGF(4, " output dimensions: %d", dims);
86 
87  pa = ptarray_construct((dims == 3), 0, size);
88 
89  for (i = 0; i < size; i++)
90  {
91  GEOSCoordSeq_getX(cs, i, &(point.x));
92  GEOSCoordSeq_getY(cs, i, &(point.y));
93  if (dims >= 3) GEOSCoordSeq_getZ(cs, i, &(point.z));
94  ptarray_set_point4d(pa, i, &point);
95  }
96 
97  return pa;
98 }
99 
100 /* Return an LWGEOM from a Geometry */
101 LWGEOM*
102 GEOS2LWGEOM(const GEOSGeometry* geom, uint8_t want3d)
103 {
104  int type = GEOSGeomTypeId(geom);
105  int SRID = GEOSGetSRID(geom);
106 
107  /* GEOS's 0 is equivalent to our unknown as for SRID values */
108  if (SRID == 0) SRID = SRID_UNKNOWN;
109 
110  if (want3d && !GEOSHasZ(geom))
111  {
112  LWDEBUG(3, "Geometry has no Z, won't provide one");
113  want3d = 0;
114  }
115 
116  switch (type)
117  {
118  const GEOSCoordSequence* cs;
119  POINTARRAY *pa, **ppaa;
120  const GEOSGeometry* g;
121  LWGEOM** geoms;
122  uint32_t i, ngeoms;
123 
124  case GEOS_POINT:
125  LWDEBUG(4, "lwgeom_from_geometry: it's a Point");
126  cs = GEOSGeom_getCoordSeq(geom);
127  if (GEOSisEmpty(geom)) return (LWGEOM*)lwpoint_construct_empty(SRID, want3d, 0);
128  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
129  return (LWGEOM*)lwpoint_construct(SRID, NULL, pa);
130 
131  case GEOS_LINESTRING:
132  case GEOS_LINEARRING:
133  LWDEBUG(4, "lwgeom_from_geometry: it's a LineString or LinearRing");
134  if (GEOSisEmpty(geom)) return (LWGEOM*)lwline_construct_empty(SRID, want3d, 0);
135 
136  cs = GEOSGeom_getCoordSeq(geom);
137  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
138  return (LWGEOM*)lwline_construct(SRID, NULL, pa);
139 
140  case GEOS_POLYGON:
141  LWDEBUG(4, "lwgeom_from_geometry: it's a Polygon");
142  if (GEOSisEmpty(geom)) return (LWGEOM*)lwpoly_construct_empty(SRID, want3d, 0);
143  ngeoms = GEOSGetNumInteriorRings(geom);
144  ppaa = lwalloc(sizeof(POINTARRAY*) * (ngeoms + 1));
145  g = GEOSGetExteriorRing(geom);
146  cs = GEOSGeom_getCoordSeq(g);
147  ppaa[0] = ptarray_from_GEOSCoordSeq(cs, want3d);
148  for (i = 0; i < ngeoms; i++)
149  {
150  g = GEOSGetInteriorRingN(geom, i);
151  cs = GEOSGeom_getCoordSeq(g);
152  ppaa[i + 1] = ptarray_from_GEOSCoordSeq(cs, want3d);
153  }
154  return (LWGEOM*)lwpoly_construct(SRID, NULL, ngeoms + 1, ppaa);
155 
156  case GEOS_MULTIPOINT:
157  case GEOS_MULTILINESTRING:
158  case GEOS_MULTIPOLYGON:
159  case GEOS_GEOMETRYCOLLECTION:
160  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
161 
162  ngeoms = GEOSGetNumGeometries(geom);
163  geoms = NULL;
164  if (ngeoms)
165  {
166  geoms = lwalloc(sizeof(LWGEOM*) * ngeoms);
167  for (i = 0; i < ngeoms; i++)
168  {
169  g = GEOSGetGeometryN(geom, i);
170  geoms[i] = GEOS2LWGEOM(g, want3d);
171  }
172  }
173  return (LWGEOM*)lwcollection_construct(type, SRID, NULL, ngeoms, geoms);
174 
175  default:
176  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
177  return NULL;
178  }
179 }
180 
181 GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY*, uint8_t fix_ring);
182 
183 GEOSCoordSeq
185 {
186  uint32_t dims = 2;
187  uint32_t i;
188  int append_points = 0;
189  const POINT3DZ* p3d;
190  const POINT2D* p2d;
191  GEOSCoordSeq sq;
192 
193  if (FLAGS_GET_Z(pa->flags)) dims = 3;
194 
195  if (fix_ring)
196  {
197  if (pa->npoints < 1)
198  {
199  lwerror("ptarray_to_GEOSCoordSeq called with fix_ring and 0 vertices in ring, cannot fix");
200  return NULL;
201  }
202  else
203  {
204  if (pa->npoints < 4) append_points = 4 - pa->npoints;
205  if (!ptarray_is_closed_2d(pa) && append_points == 0) append_points = 1;
206  }
207  }
208 
209  if (!(sq = GEOSCoordSeq_create(pa->npoints + append_points, dims)))
210  {
211  lwerror("Error creating GEOS Coordinate Sequence");
212  return NULL;
213  }
214 
215  for (i = 0; i < pa->npoints; i++)
216  {
217  if (dims == 3)
218  {
219  p3d = getPoint3dz_cp(pa, i);
220  p2d = (const POINT2D*)p3d;
221  LWDEBUGF(4, "Point: %g,%g,%g", p3d->x, p3d->y, p3d->z);
222  }
223  else
224  {
225  p2d = getPoint2d_cp(pa, i);
226  LWDEBUGF(4, "Point: %g,%g", p2d->x, p2d->y);
227  }
228 
229  GEOSCoordSeq_setX(sq, i, p2d->x);
230  GEOSCoordSeq_setY(sq, i, p2d->y);
231 
232  if (dims == 3) GEOSCoordSeq_setZ(sq, i, p3d->z);
233  }
234 
235  if (append_points)
236  {
237  if (dims == 3)
238  {
239  p3d = getPoint3dz_cp(pa, 0);
240  p2d = (const POINT2D*)p3d;
241  }
242  else
243  p2d = getPoint2d_cp(pa, 0);
244  for (i = pa->npoints; i < pa->npoints + append_points; i++)
245  {
246  GEOSCoordSeq_setX(sq, i, p2d->x);
247  GEOSCoordSeq_setY(sq, i, p2d->y);
248 
249  if (dims == 3) GEOSCoordSeq_setZ(sq, i, p3d->z);
250  }
251  }
252 
253  return sq;
254 }
255 
256 static inline GEOSGeometry*
258 {
259  GEOSCoordSeq sq;
260  GEOSGeom g;
261  sq = ptarray_to_GEOSCoordSeq(pa, autofix);
262  g = GEOSGeom_createLinearRing(sq);
263  return g;
264 }
265 
266 GEOSGeometry*
267 GBOX2GEOS(const GBOX* box)
268 {
269  GEOSGeometry* envelope;
270  GEOSGeometry* ring;
271  GEOSCoordSequence* seq = GEOSCoordSeq_create(5, 2);
272  if (!seq) return NULL;
273 
274  GEOSCoordSeq_setX(seq, 0, box->xmin);
275  GEOSCoordSeq_setY(seq, 0, box->ymin);
276 
277  GEOSCoordSeq_setX(seq, 1, box->xmax);
278  GEOSCoordSeq_setY(seq, 1, box->ymin);
279 
280  GEOSCoordSeq_setX(seq, 2, box->xmax);
281  GEOSCoordSeq_setY(seq, 2, box->ymax);
282 
283  GEOSCoordSeq_setX(seq, 3, box->xmin);
284  GEOSCoordSeq_setY(seq, 3, box->ymax);
285 
286  GEOSCoordSeq_setX(seq, 4, box->xmin);
287  GEOSCoordSeq_setY(seq, 4, box->ymin);
288 
289  ring = GEOSGeom_createLinearRing(seq);
290  if (!ring)
291  {
292  GEOSCoordSeq_destroy(seq);
293  return NULL;
294  }
295 
296  envelope = GEOSGeom_createPolygon(ring, NULL, 0);
297  if (!envelope)
298  {
299  GEOSGeom_destroy(ring);
300  return NULL;
301  }
302 
303  return envelope;
304 }
305 
306 GEOSGeometry*
307 LWGEOM2GEOS(const LWGEOM* lwgeom, uint8_t autofix)
308 {
309  GEOSCoordSeq sq;
310  GEOSGeom g, shell;
311  GEOSGeom* geoms = NULL;
312  uint32_t ngeoms, i, j;
313  int geostype;
314 #if LWDEBUG_LEVEL >= 4
315  char* wkt;
316 #endif
317 
318  if (autofix)
319  {
320  /* cross fingers and try without autofix, maybe it'll work? */
321  g = LWGEOM2GEOS(lwgeom, LW_FALSE);
322  if (g) return g;
323  }
324 
325  LWDEBUGF(4, "LWGEOM2GEOS got a %s", lwtype_name(lwgeom->type));
326 
327  if (lwgeom_has_arc(lwgeom))
328  {
329  LWGEOM* lwgeom_stroked = lwgeom_stroke(lwgeom, 32);
330  GEOSGeometry* g = LWGEOM2GEOS(lwgeom_stroked, autofix);
331  lwgeom_free(lwgeom_stroked);
332  return g;
333  }
334 
335  LWPOINT* lwp = NULL;
336  LWPOLY* lwpoly = NULL;
337  LWLINE* lwl = NULL;
338  LWCOLLECTION* lwc = NULL;
339 
340  switch (lwgeom->type)
341  {
342  case POINTTYPE:
343  lwp = (LWPOINT*)lwgeom;
344 
345  if (lwgeom_is_empty(lwgeom))
346  g = GEOSGeom_createEmptyPolygon();
347  else
348  {
349  sq = ptarray_to_GEOSCoordSeq(lwp->point, 0);
350  g = GEOSGeom_createPoint(sq);
351  }
352  if (!g) return NULL;
353  break;
354 
355  case LINETYPE:
356  lwl = (LWLINE*)lwgeom;
357  /* TODO: if (autofix) */
358  if (lwl->points->npoints == 1)
359  {
360  /* Duplicate point, to make geos-friendly */
361  lwl->points = ptarray_addPoint(lwl->points,
362  getPoint_internal(lwl->points, 0),
363  FLAGS_NDIMS(lwl->points->flags),
364  lwl->points->npoints);
365  }
366  sq = ptarray_to_GEOSCoordSeq(lwl->points, 0);
367  g = GEOSGeom_createLineString(sq);
368  if (!g) return NULL;
369  break;
370 
371  case POLYGONTYPE:
372  lwpoly = (LWPOLY*)lwgeom;
373  if (lwgeom_is_empty(lwgeom))
374  g = GEOSGeom_createEmptyPolygon();
375  else
376  {
377  shell = ptarray_to_GEOSLinearRing(lwpoly->rings[0], autofix);
378  if (!shell) return NULL;
379  ngeoms = lwpoly->nrings - 1;
380  if (ngeoms > 0) geoms = malloc(sizeof(GEOSGeom) * ngeoms);
381 
382  for (i = 1; i < lwpoly->nrings; i++)
383  {
384  geoms[i - 1] = ptarray_to_GEOSLinearRing(lwpoly->rings[i], autofix);
385  if (!geoms[i - 1])
386  {
387  uint32_t k;
388  for (k = 0; k < i - 1; k++)
389  GEOSGeom_destroy(geoms[k]);
390  free(geoms);
391  GEOSGeom_destroy(shell);
392  return NULL;
393  }
394  }
395  g = GEOSGeom_createPolygon(shell, geoms, ngeoms);
396  if (geoms) free(geoms);
397  }
398  if (!g) return NULL;
399  break;
400  case MULTIPOINTTYPE:
401  case MULTILINETYPE:
402  case MULTIPOLYGONTYPE:
403  case COLLECTIONTYPE:
404  if (lwgeom->type == MULTIPOINTTYPE)
405  geostype = GEOS_MULTIPOINT;
406  else if (lwgeom->type == MULTILINETYPE)
407  geostype = GEOS_MULTILINESTRING;
408  else if (lwgeom->type == MULTIPOLYGONTYPE)
409  geostype = GEOS_MULTIPOLYGON;
410  else
411  geostype = GEOS_GEOMETRYCOLLECTION;
412 
413  lwc = (LWCOLLECTION*)lwgeom;
414 
415  ngeoms = lwc->ngeoms;
416  if (ngeoms > 0) geoms = malloc(sizeof(GEOSGeom) * ngeoms);
417 
418  j = 0;
419  for (i = 0; i < ngeoms; ++i)
420  {
421  GEOSGeometry* g;
422 
423  if (lwgeom_is_empty(lwc->geoms[i])) continue;
424 
425  g = LWGEOM2GEOS(lwc->geoms[i], 0);
426  if (!g)
427  {
428  uint32_t k;
429  for (k = 0; k < j; k++)
430  GEOSGeom_destroy(geoms[k]);
431  free(geoms);
432  return NULL;
433  }
434  geoms[j++] = g;
435  }
436  g = GEOSGeom_createCollection(geostype, geoms, j);
437  if (ngeoms > 0) free(geoms);
438  if (!g) return NULL;
439  break;
440 
441  default:
442  lwerror("Unknown geometry type: %d - %s", lwgeom->type, lwtype_name(lwgeom->type));
443  return NULL;
444  }
445 
446  GEOSSetSRID(g, lwgeom->srid);
447 
448 #if LWDEBUG_LEVEL >= 4
449  wkt = GEOSGeomToWKT(g);
450  LWDEBUGF(4, "LWGEOM2GEOS: GEOSGeom: %s", wkt);
451  free(wkt);
452 #endif
453 
454  return g;
455 }
456 
457 GEOSGeometry*
458 make_geos_point(double x, double y)
459 {
460  GEOSCoordSequence* seq = GEOSCoordSeq_create(1, 2);
461  GEOSGeometry* geom = NULL;
462 
463  if (!seq) return NULL;
464 
465  GEOSCoordSeq_setX(seq, 0, x);
466  GEOSCoordSeq_setY(seq, 0, y);
467 
468  geom = GEOSGeom_createPoint(seq);
469  if (!geom) GEOSCoordSeq_destroy(seq);
470  return geom;
471 }
472 
473 GEOSGeometry*
474 make_geos_segment(double x1, double y1, double x2, double y2)
475 {
476  GEOSCoordSequence* seq = GEOSCoordSeq_create(2, 2);
477  GEOSGeometry* geom = NULL;
478 
479  if (!seq) return NULL;
480 
481  GEOSCoordSeq_setX(seq, 0, x1);
482  GEOSCoordSeq_setY(seq, 0, y1);
483  GEOSCoordSeq_setX(seq, 1, x2);
484  GEOSCoordSeq_setY(seq, 1, y2);
485 
486  geom = GEOSGeom_createLineString(seq);
487  if (!geom) GEOSCoordSeq_destroy(seq);
488  return geom;
489 }
490 
491 const char*
493 {
494  const char* ver = GEOSversion();
495  return ver;
496 }
497 
498 inline static int32_t
499 get_result_srid(const LWGEOM* geom1, const LWGEOM* geom2, const char* funcname)
500 {
501  if (!geom1)
502  {
503  lwerror("%s: First argument is null pointer", funcname);
504  return SRID_INVALID;
505  }
506  if (geom2 && (geom1->srid != geom2->srid))
507  {
508  lwerror("%s: Operation on mixed SRID geometries (%d != %d)", funcname, geom1->srid, geom2->srid);
509  return SRID_INVALID;
510  }
511  if (geom1->srid > SRID_MAXIMUM)
512  {
513  lwerror("%s: SRID is more than maximum (%d > %d)", funcname, geom1->srid, SRID_USER_MAXIMUM);
514  return SRID_INVALID;
515  }
516  return geom1->srid;
517 }
518 
519 /* Input decoder and sanity checker for GEOS wrappers */
520 inline static uint8_t
521 input_lwgeom_to_geos(GEOSGeometry** g, const LWGEOM* geom, const char* funcname)
522 {
523  *g = LWGEOM2GEOS(geom, 1);
524  if (!*g)
525  {
526  lwerror("%s: argument geometry could not be converted to GEOS: %s", funcname, lwgeom_geos_errmsg);
527  return LW_FALSE;
528  }
529  return LW_TRUE;
530 }
531 
532 /* Output encoder and sanity checker for GEOS wrappers */
533 inline static uint8_t
534 output_geos_as_lwgeom(GEOSGeometry** g, LWGEOM** geom, const int32_t srid, const uint8_t is3d, const char* funcname)
535 {
536  GEOSSetSRID(*g, srid);
537  *geom = GEOS2LWGEOM(*g, is3d);
538  if (!*geom)
539  {
540  lwerror("%s: result geometry could not be converted from GEOS: %s", funcname, lwgeom_geos_errmsg);
541  return LW_FALSE;
542  }
543  return LW_TRUE;
544 }
545 
546 /* Output encoder and sanity checker for GEOS wrappers */
547 inline static LWGEOM*
548 geos_clean_and_fail(GEOSGeometry* g1, GEOSGeometry* g2, GEOSGeometry* g3, const char* funcname)
549 {
550  if (g1) GEOSGeom_destroy(g1);
551  if (g2) GEOSGeom_destroy(g2);
552  if (g3) GEOSGeom_destroy(g3);
553  lwerror("%s: GEOS Error: %s", funcname, lwgeom_geos_errmsg);
554  return NULL;
555 }
556 
557 /* Output encoder and sanity checker for GEOS wrappers */
558 inline static void
559 geos_clean(GEOSGeometry* g1, GEOSGeometry* g2, GEOSGeometry* g3)
560 {
561  if (g1) GEOSGeom_destroy(g1);
562  if (g2) GEOSGeom_destroy(g2);
563  if (g3) GEOSGeom_destroy(g3);
564  return;
565 }
566 
567 LWGEOM*
569 {
570  LWGEOM* result;
571  int32_t srid = get_result_srid(geom, NULL, __func__);
572  uint8_t is3d = FLAGS_GET_Z(geom->flags);
573  GEOSGeometry* g;
574 
575  if (srid == SRID_INVALID) return NULL;
576 
577  initGEOS(lwnotice, lwgeom_geos_error);
578 
579  if (!input_lwgeom_to_geos(&g, geom, __func__)) return NULL;
580 
581  if (GEOSNormalize(g) == -1) return geos_clean_and_fail(g, NULL, NULL, __func__);
582 
583  output_geos_as_lwgeom(&g, &result, srid, is3d, __func__);
584  GEOSGeom_destroy(g);
585  return result;
586 }
587 
588 LWGEOM*
589 lwgeom_intersection(const LWGEOM* geom1, const LWGEOM* geom2)
590 {
591  LWGEOM* result;
592  int32_t srid = get_result_srid(geom1, geom2, __func__);
593  uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
594  GEOSGeometry* g1;
595  GEOSGeometry* g2;
596  GEOSGeometry* g3;
597 
598  if (srid == SRID_INVALID) return NULL;
599 
600  /* A.Intersection(Empty) == Empty */
601  if (lwgeom_is_empty(geom2)) return lwgeom_clone_deep(geom2); /* match empty type? */
602 
603  /* Empty.Intersection(A) == Empty */
604  if (lwgeom_is_empty(geom1)) return lwgeom_clone_deep(geom1); /* match empty type? */
605 
606  initGEOS(lwnotice, lwgeom_geos_error);
607 
608  if (!input_lwgeom_to_geos(&g1, geom1, __func__)) return NULL;
609  if (!input_lwgeom_to_geos(&g2, geom2, __func__)) return geos_clean_and_fail(g1, NULL, NULL, __func__);
610 
611  g3 = GEOSIntersection(g1, g2);
612 
613  if (!g3)
614  {
615  GEOSGeometry *g1v, *g2v;
616  lwnotice("%s: GEOS Error: %s", __func__, lwgeom_geos_errmsg);
617 
618  if (!GEOSisValid(g1) || !GEOSisValid(g2))
619  {
620  lwnotice(
621  "Your geometry dataset is not valid per OGC Specification. "
622  "Please fix it with manual review of entries that are not ST_IsValid(geom). "
623  "Retrying GEOS operation with ST_MakeValid of your input.");
624  g1v = LWGEOM_GEOS_makeValid(g1);
625  g2v = LWGEOM_GEOS_makeValid(g2);
626  g3 = GEOSIntersection(g1v, g2v);
627  geos_clean(g1v, g2v, NULL);
628  }
629  }
630 
631  if (!g3) return geos_clean_and_fail(g1, g2, NULL, __func__);
632 
633  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
634  return geos_clean_and_fail(g1, g2, g3, __func__);
635 
636  geos_clean(g1, g2, g3);
637  return result;
638 }
639 
640 LWGEOM*
642 {
643  LWGEOM* result;
644  int32_t srid = get_result_srid(geom, NULL, __func__);
645  uint8_t is3d = FLAGS_GET_Z(geom->flags);
646  GEOSGeometry* g1;
647  GEOSGeometry* g3;
648 
649  if (srid == SRID_INVALID) return NULL;
650 
651  /* Empty.Linemerge() == Empty */
652  if (lwgeom_is_empty(geom)) return lwgeom_clone_deep(geom); /* match empty type to linestring? */
653 
654  initGEOS(lwnotice, lwgeom_geos_error);
655 
656  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
657 
658  g3 = GEOSLineMerge(g1);
659 
660  if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
661 
662  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
663  return geos_clean_and_fail(g1, NULL, g3, __func__);
664 
665  geos_clean(g1, NULL, g3);
666 
667  return result;
668 }
669 
670 LWGEOM*
672 {
673  LWGEOM* result;
674  int32_t srid = get_result_srid(geom, NULL, __func__);
675  uint8_t is3d = FLAGS_GET_Z(geom->flags);
676  GEOSGeometry* g1;
677  GEOSGeometry* g3;
678 
679  if (srid == SRID_INVALID) return NULL;
680 
681  /* Empty.UnaryUnion() == Empty */
682  if (lwgeom_is_empty(geom)) return lwgeom_clone_deep(geom);
683 
684  initGEOS(lwnotice, lwgeom_geos_error);
685 
686  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
687 
688  g3 = GEOSUnaryUnion(g1);
689 
690  if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
691 
692  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
693  return geos_clean_and_fail(g1, NULL, g3, __func__);
694 
695  geos_clean(g1, NULL, g3);
696 
697  return result;
698 }
699 
700 LWGEOM*
701 lwgeom_difference(const LWGEOM* geom1, const LWGEOM* geom2)
702 {
703  LWGEOM* result;
704  int32_t srid = get_result_srid(geom1, geom2, __func__);
705  uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
706  GEOSGeometry *g1, *g2, *g3;
707 
708  if (srid == SRID_INVALID) return NULL;
709 
710  /* A.Intersection(Empty) == Empty */
711  if (lwgeom_is_empty(geom2)) return lwgeom_clone_deep(geom1); /* match empty type? */
712 
713  /* Empty.Intersection(A) == Empty */
714  if (lwgeom_is_empty(geom1)) return lwgeom_clone_deep(geom1); /* match empty type? */
715 
716  initGEOS(lwnotice, lwgeom_geos_error);
717 
718  if (!input_lwgeom_to_geos(&g1, geom1, __func__)) return NULL;
719  if (!input_lwgeom_to_geos(&g2, geom2, __func__)) return geos_clean_and_fail(g1, NULL, NULL, __func__);
720 
721  g3 = GEOSDifference(g1, g2);
722 
723  if (!g3)
724  {
725  GEOSGeometry *g1v, *g2v;
726  lwnotice("%s: GEOS Error: %s", __func__, lwgeom_geos_errmsg);
727 
728  if (!GEOSisValid(g1) || !GEOSisValid(g2))
729  {
730  lwnotice(
731  "Your geometry dataset is not valid per OGC Specification. "
732  "Please fix it with manual review of entries that are not ST_IsValid(geom). "
733  "Retrying GEOS operation with ST_MakeValid of your input.");
734  g1v = LWGEOM_GEOS_makeValid(g1);
735  g2v = LWGEOM_GEOS_makeValid(g2);
736  g3 = GEOSDifference(g1v, g2v);
737  geos_clean(g1v, g2v, NULL);
738  }
739  }
740 
741  if (!g3) return geos_clean_and_fail(g1, g2, NULL, __func__);
742 
743  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
744  return geos_clean_and_fail(g1, g2, g3, __func__);
745 
746  geos_clean(g1, g2, g3);
747  return result;
748 }
749 
750 LWGEOM*
751 lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2)
752 {
753  LWGEOM* result;
754  int32_t srid = get_result_srid(geom1, geom2, __func__);
755  uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
756  GEOSGeometry *g1, *g2, *g3;
757 
758  if (srid == SRID_INVALID) return NULL;
759 
760  /* A.SymDifference(Empty) == A */
761  if (lwgeom_is_empty(geom2)) return lwgeom_clone_deep(geom1);
762 
763  /* Empty.DymDifference(B) == B */
764  if (lwgeom_is_empty(geom1)) return lwgeom_clone_deep(geom2);
765 
766  initGEOS(lwnotice, lwgeom_geos_error);
767 
768  if (!input_lwgeom_to_geos(&g1, geom1, __func__)) return NULL;
769  if (!input_lwgeom_to_geos(&g2, geom2, __func__)) return geos_clean_and_fail(g1, NULL, NULL, __func__);
770 
771  g3 = GEOSSymDifference(g1, g2);
772 
773  if (!g3)
774  {
775  GEOSGeometry *g1v, *g2v;
776  lwnotice("%s: GEOS Error: %s", __func__, lwgeom_geos_errmsg);
777 
778  if (!GEOSisValid(g1) || !GEOSisValid(g2))
779  {
780  lwnotice(
781  "Your geometry dataset is not valid per OGC Specification. "
782  "Please fix it with manual review of entries that are not ST_IsValid(geom). "
783  "Retrying GEOS operation with ST_MakeValid of your input.");
784  g1v = LWGEOM_GEOS_makeValid(g1);
785  g2v = LWGEOM_GEOS_makeValid(g2);
786  g3 = GEOSSymDifference(g1v, g2v);
787  geos_clean(g1v, g2v, NULL);
788  }
789  }
790 
791  if (!g3) return geos_clean_and_fail(g1, g2, NULL, __func__);
792 
793  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
794  return geos_clean_and_fail(g1, g2, g3, __func__);
795 
796  geos_clean(g1, g2, g3);
797  return result;
798 }
799 
800 LWGEOM*
802 {
803  LWGEOM* result;
804  int32_t srid = get_result_srid(geom, NULL, __func__);
805  uint8_t is3d = FLAGS_GET_Z(geom->flags);
806  GEOSGeometry *g1, *g3;
807 
808  if (srid == SRID_INVALID) return NULL;
809 
810  if (lwgeom_is_empty(geom))
811  {
812  LWPOINT* lwp = lwpoint_construct_empty(srid, is3d, lwgeom_has_m(geom));
813  return lwpoint_as_lwgeom(lwp);
814  }
815 
816  initGEOS(lwnotice, lwgeom_geos_error);
817 
818  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
819 
820  g3 = GEOSGetCentroid(g1);
821 
822  if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
823 
824  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
825  return geos_clean_and_fail(g1, NULL, g3, __func__);
826 
827  geos_clean(g1, NULL, g3);
828 
829  return result;
830 }
831 
832 LWGEOM*
833 lwgeom_union(const LWGEOM* geom1, const LWGEOM* geom2)
834 {
835  LWGEOM* result;
836  int32_t srid = get_result_srid(geom1, geom2, __func__);
837  uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
838  GEOSGeometry *g1, *g2, *g3;
839 
840  if (srid == SRID_INVALID) return NULL;
841 
842  /* A.Union(empty) == A */
843  if (lwgeom_is_empty(geom1)) return lwgeom_clone_deep(geom2);
844 
845  /* B.Union(empty) == B */
846  if (lwgeom_is_empty(geom2)) return lwgeom_clone_deep(geom1);
847 
848  initGEOS(lwnotice, lwgeom_geos_error);
849 
850  if (!input_lwgeom_to_geos(&g1, geom1, __func__)) return NULL;
851  if (!input_lwgeom_to_geos(&g2, geom2, __func__)) return geos_clean_and_fail(g1, NULL, NULL, __func__);
852 
853  g3 = GEOSUnion(g1, g2);
854 
855  if (!g3)
856  {
857  GEOSGeometry *g1v, *g2v;
858  lwnotice("%s: GEOS Error: %s", __func__, lwgeom_geos_errmsg);
859 
860  if (!GEOSisValid(g1) || !GEOSisValid(g2))
861  {
862  lwnotice(
863  "Your geometry dataset is not valid per OGC Specification. "
864  "Please fix it with manual review of entries that are not ST_IsValid(geom). "
865  "Retrying GEOS operation with ST_MakeValid of your input.");
866  g1v = LWGEOM_GEOS_makeValid(g1);
867  g2v = LWGEOM_GEOS_makeValid(g2);
868  g3 = GEOSUnion(g1v, g2v);
869  geos_clean(g1v, g2v, NULL);
870  }
871  }
872 
873  if (!g3) return geos_clean_and_fail(g1, g2, NULL, __func__);
874 
875  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
876  return geos_clean_and_fail(g1, g2, g3, __func__);
877 
878  geos_clean(g1, g2, g3);
879  return result;
880 }
881 
882 LWGEOM *
883 lwgeom_clip_by_rect(const LWGEOM *geom1, double x1, double y1, double x2, double y2)
884 {
885  LWGEOM *result;
886  LWGEOM *tmp;
887  LWGEOM *envelope = (LWGEOM*)lwpoly_construct_envelope(geom1->srid, x1, y1, x2, y2);
888 
889  result = lwgeom_intersection(geom1, envelope);
890  lwgeom_free(envelope);
891 
892  if (!result) return NULL;
893 
894  /* clipping should not produce lower dimension objects */
895  if (
896  /* input has exact dimensionality, isn't a generic collection */
897  geom1->type != COLLECTIONTYPE &&
898  /* output may have different things inside */
899  result->type == COLLECTIONTYPE)
900  {
902  lwfree(result);
903  result = tmp;
904  if (!result) return NULL;
905  }
906 
907  /* clean up stray points on geometry boundary */
908  lwgeom_simplify_in_place(result, 0.0, LW_TRUE);
909 
910  return result;
911 
912 #if 0 /* POSTGIS_GEOS_VERSION >= 35, enable only after bugs in geos are fixed */
913  int32_t srid = get_result_srid(geom, NULL, __func__);
914  uint8_t is3d = FLAGS_GET_Z(geom->flags);
915  GEOSGeometry *g1, *g3;
916 
917  if (srid == SRID_INVALID) return NULL;
918 
919  /* A.Intersection(Empty) == Empty */
921 
922  initGEOS(lwnotice, lwgeom_geos_error);
923 
924  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
925 
926  g3 = GEOSClipByRect(g1, x1, y1, x2, y2);
927 
928  if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
929 
930  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
931  return geos_clean_and_fail(g1, NULL, g3, __func__);
932 
933  geos_clean(g1, NULL, g3);
934 
935  return result;
936 #endif /* POSTGIS_GEOS_VERSION >= 35 */
937 }
938 
939 /* ------------ BuildArea stuff ---------------------------------------------------------------------{ */
940 
941 typedef struct Face_t
942 {
943  const GEOSGeometry* geom;
944  GEOSGeometry* env;
945  double envarea;
946  struct Face_t* parent; /* if this face is an hole of another one, or NULL */
947 } Face;
948 
949 static Face* newFace(const GEOSGeometry* g);
950 static void delFace(Face* f);
951 static unsigned int countParens(const Face* f);
952 static void findFaceHoles(Face** faces, int nfaces);
953 
954 static Face*
955 newFace(const GEOSGeometry* g)
956 {
957  Face* f = lwalloc(sizeof(Face));
958  f->geom = g;
959  f->env = GEOSEnvelope(f->geom);
960  GEOSArea(f->env, &f->envarea);
961  f->parent = NULL;
962  return f;
963 }
964 
965 static unsigned int
966 countParens(const Face* f)
967 {
968  unsigned int pcount = 0;
969  while (f->parent)
970  {
971  ++pcount;
972  f = f->parent;
973  }
974  return pcount;
975 }
976 
977 /* Destroy the face and release memory associated with it */
978 static void
980 {
981  GEOSGeom_destroy(f->env);
982  lwfree(f);
983 }
984 
985 static int
986 compare_by_envarea(const void* g1, const void* g2)
987 {
988  Face* f1 = *(Face**)g1;
989  Face* f2 = *(Face**)g2;
990  double n1 = f1->envarea;
991  double n2 = f2->envarea;
992 
993  if (n1 < n2) return 1;
994  if (n1 > n2) return -1;
995  return 0;
996 }
997 
998 /* Find holes of each face */
999 static void
1000 findFaceHoles(Face** faces, int nfaces)
1001 {
1002  int i, j, h;
1003 
1004  /* We sort by envelope area so that we know holes are only after their shells */
1005  qsort(faces, nfaces, sizeof(Face*), compare_by_envarea);
1006  for (i = 0; i < nfaces; ++i)
1007  {
1008  Face* f = faces[i];
1009  int nholes = GEOSGetNumInteriorRings(f->geom);
1010  LWDEBUGF(2, "Scanning face %d with env area %g and %d holes", i, f->envarea, nholes);
1011  for (h = 0; h < nholes; ++h)
1012  {
1013  const GEOSGeometry* hole = GEOSGetInteriorRingN(f->geom, h);
1014  LWDEBUGF(2,
1015  "Looking for hole %d/%d of face %d among %d other faces",
1016  h + 1,
1017  nholes,
1018  i,
1019  nfaces - i - 1);
1020  for (j = i + 1; j < nfaces; ++j)
1021  {
1022  const GEOSGeometry* f2er;
1023  Face* f2 = faces[j];
1024  if (f2->parent) continue; /* hole already assigned */
1025  f2er = GEOSGetExteriorRing(f2->geom);
1026  /* TODO: can be optimized as the ring would have the same vertices, possibly in
1027  * different order. Maybe comparing number of points could already be useful. */
1028  if (GEOSEquals(f2er, hole))
1029  {
1030  LWDEBUGF(2, "Hole %d/%d of face %d is face %d", h + 1, nholes, i, j);
1031  f2->parent = f;
1032  break;
1033  }
1034  }
1035  }
1036  }
1037 }
1038 
1039 static GEOSGeometry*
1041 {
1042  GEOSGeometry** geoms = lwalloc(sizeof(GEOSGeometry*) * nfaces);
1043  GEOSGeometry* ret;
1044  unsigned int ngeoms = 0;
1045  int i;
1046 
1047  for (i = 0; i < nfaces; ++i)
1048  {
1049  Face* f = faces[i];
1050  if (countParens(f) % 2) continue; /* we skip odd parents geoms */
1051  geoms[ngeoms++] = GEOSGeom_clone(f->geom);
1052  }
1053 
1054  ret = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, ngeoms);
1055  lwfree(geoms);
1056  return ret;
1057 }
1058 
1059 GEOSGeometry*
1060 LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in)
1061 {
1062  GEOSGeometry* tmp;
1063  GEOSGeometry *geos_result, *shp;
1064  GEOSGeometry const* vgeoms[1];
1065  uint32_t i, ngeoms;
1066  int srid = GEOSGetSRID(geom_in);
1067  Face** geoms;
1068 #if POSTGIS_DEBUG_LEVEL >= 3
1069  LWGEOM *geos_geom;
1070  char *geom_ewkt;
1071 #endif
1072 
1073  vgeoms[0] = geom_in;
1074  geos_result = GEOSPolygonize(vgeoms, 1);
1075 
1076  LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result);
1077 
1078  /* Null return from GEOSpolygonize (an exception) */
1079  if (!geos_result) return 0;
1080 
1081  /* We should now have a collection */
1082 #if PARANOIA_LEVEL > 0
1083  if (GEOSGeomTypeId(geos_result) != COLLECTIONTYPE)
1084  {
1085  GEOSGeom_destroy(geos_result);
1086  lwerror("%s [%d] Unexpected return from GEOSpolygonize", __FILE__, __LINE__);
1087  return 0;
1088  }
1089 #endif
1090 
1091  ngeoms = GEOSGetNumGeometries(geos_result);
1092 
1093 #if POSTGIS_DEBUG_LEVEL >= 3
1094  LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
1095  geos_geom = GEOS2LWGEOM(geos_result, 0);
1096  geom_ewkt = lwgeom_to_ewkt(geos_geom);
1097  LWDEBUGF(3, "GEOSpolygonize: polygonized:%s", geom_ewkt);
1098  lwgeom_free(geos_geom);
1099  lwfree(geom_ewkt);
1100 #endif
1101 
1102  /* No geometries in collection, early out */
1103  if (ngeoms == 0)
1104  {
1105  GEOSSetSRID(geos_result, srid);
1106  return geos_result;
1107  }
1108 
1109  /* Return first geometry if we only have one in collection, to avoid the unnecessary Geometry clone below. */
1110  if (ngeoms == 1)
1111  {
1112  tmp = (GEOSGeometry*)GEOSGetGeometryN(geos_result, 0);
1113  if (!tmp)
1114  {
1115  GEOSGeom_destroy(geos_result);
1116  return 0; /* exception */
1117  }
1118  shp = GEOSGeom_clone(tmp);
1119  GEOSGeom_destroy(geos_result); /* only safe after the clone above */
1120  GEOSSetSRID(shp, srid);
1121  return shp;
1122  }
1123 
1124  LWDEBUGF(2, "Polygonize returned %d geoms", ngeoms);
1125 
1126  /*
1127  * Polygonizer returns a polygon for each face in the built topology.
1128  *
1129  * This means that for any face with holes we'll have other faces representing each hole. We can imagine a
1130  * parent-child relationship between these faces.
1131  *
1132  * In order to maximize the number of visible rings in output we only use those faces which have an even number
1133  * of parents.
1134  *
1135  * Example:
1136  *
1137  * +---------------+
1138  * | L0 | L0 has no parents
1139  * | +---------+ |
1140  * | | L1 | | L1 is an hole of L0
1141  * | | +---+ | |
1142  * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0)
1143  * | | | | | |
1144  * | | +---+ | |
1145  * | +---------+ |
1146  * | |
1147  * +---------------+
1148  *
1149  * See http://trac.osgeo.org/postgis/ticket/1806
1150  *
1151  */
1152 
1153  /* Prepare face structures for later analysis */
1154  geoms = lwalloc(sizeof(Face**) * ngeoms);
1155  for (i = 0; i < ngeoms; ++i)
1156  geoms[i] = newFace(GEOSGetGeometryN(geos_result, i));
1157 
1158  /* Find faces representing other faces holes */
1159  findFaceHoles(geoms, ngeoms);
1160 
1161  /* Build a MultiPolygon composed only by faces with an even number of ancestors */
1162  tmp = collectFacesWithEvenAncestors(geoms, ngeoms);
1163 
1164  /* Cleanup face structures */
1165  for (i = 0; i < ngeoms; ++i)
1166  delFace(geoms[i]);
1167  lwfree(geoms);
1168 
1169  /* Faces referenced memory owned by geos_result. It is safe to destroy geos_result after deleting them. */
1170  GEOSGeom_destroy(geos_result);
1171 
1172  /* Run a single overlay operation to dissolve shared edges */
1173  shp = GEOSUnionCascaded(tmp);
1174  if (!shp)
1175  {
1176  GEOSGeom_destroy(tmp);
1177  return 0; /* exception */
1178  }
1179 
1180  GEOSGeom_destroy(tmp);
1181 
1182  GEOSSetSRID(shp, srid);
1183 
1184  return shp;
1185 }
1186 
1187 LWGEOM*
1189 {
1190  LWGEOM* result;
1191  int32_t srid = get_result_srid(geom, NULL, __func__);
1192  uint8_t is3d = FLAGS_GET_Z(geom->flags);
1193  GEOSGeometry *g1, *g3;
1194 
1195  if (srid == SRID_INVALID) return NULL;
1196 
1197  /* Can't build an area from an empty! */
1198  if (lwgeom_is_empty(geom)) return (LWGEOM*)lwpoly_construct_empty(srid, is3d, 0);
1199 
1200  initGEOS(lwnotice, lwgeom_geos_error);
1201 
1202  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
1203 
1204  g3 = LWGEOM_GEOS_buildArea(g1);
1205 
1206  if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
1207 
1208  /* If no geometries are in result collection, return NULL */
1209  if (GEOSGetNumGeometries(g3) == 0) return geos_clean_and_fail(g1, NULL, NULL, __func__);
1210 
1211  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
1212  return geos_clean_and_fail(g1, NULL, g3, __func__);
1213 
1214  geos_clean(g1, NULL, g3);
1215 
1216  return result;
1217 }
1218 
1219 /* ------------ end of BuildArea stuff ---------------------------------------------------------------------} */
1220 
1221 int
1223 {
1224  GEOSGeometry* g;
1225  int simple;
1226 
1227  /* Empty is always simple */
1228  if (lwgeom_is_empty(geom)) return LW_TRUE;
1229 
1230  initGEOS(lwnotice, lwgeom_geos_error);
1231 
1232  if (!input_lwgeom_to_geos(&g, geom, __func__)) return -1;
1233 
1234  simple = GEOSisSimple(g);
1235  GEOSGeom_destroy(g);
1236 
1237  if (simple == 2) /* exception thrown */
1238  {
1239  lwerror("lwgeom_is_simple: %s", lwgeom_geos_errmsg);
1240  return -1;
1241  }
1242 
1243  return simple ? LW_TRUE : LW_FALSE;
1244 }
1245 
1246 LWGEOM*
1248 {
1249  LWGEOM* result;
1250  int32_t srid = get_result_srid(geom, NULL, __func__);
1251  uint8_t is3d = FLAGS_GET_Z(geom->flags);
1252  GEOSGeometry* g;
1253 
1254  if (srid == SRID_INVALID) return NULL;
1255 
1256  initGEOS(lwnotice, lwgeom_geos_error);
1257 
1258  if (!input_lwgeom_to_geos(&g, geom, __func__)) return NULL;
1259 
1260  if (!g) return geos_clean_and_fail(g, NULL, NULL, __func__);
1261 
1262  if (!output_geos_as_lwgeom(&g, &result, srid, is3d, __func__))
1263  return geos_clean_and_fail(g, NULL, NULL, __func__);
1264 
1265  geos_clean(g, NULL, NULL);
1266 
1267  return result;
1268 }
1269 
1270 LWGEOM*
1271 lwgeom_snap(const LWGEOM* geom1, const LWGEOM* geom2, double tolerance)
1272 {
1273  LWGEOM* result;
1274  int32_t srid = get_result_srid(geom1, geom2, __func__);
1275  uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
1276  GEOSGeometry *g1, *g2, *g3;
1277 
1278  if (srid == SRID_INVALID) return NULL;
1279 
1280  initGEOS(lwnotice, lwgeom_geos_error);
1281 
1282  if (!input_lwgeom_to_geos(&g1, geom1, __func__)) return NULL;
1283  if (!input_lwgeom_to_geos(&g2, geom2, __func__)) return geos_clean_and_fail(g1, NULL, NULL, __func__);
1284 
1285  g3 = GEOSSnap(g1, g2, tolerance);
1286 
1287  if (!g3) return geos_clean_and_fail(g1, g2, NULL, __func__);
1288 
1289  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
1290  return geos_clean_and_fail(g1, g2, g3, __func__);
1291 
1292  geos_clean(g1, g2, g3);
1293  return result;
1294 }
1295 
1296 LWGEOM*
1297 lwgeom_sharedpaths(const LWGEOM* geom1, const LWGEOM* geom2)
1298 {
1299  LWGEOM* result;
1300  int32_t srid = get_result_srid(geom1, geom2, __func__);
1301  uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
1302  GEOSGeometry *g1, *g2, *g3;
1303 
1304  if (srid == SRID_INVALID) return NULL;
1305 
1306  initGEOS(lwnotice, lwgeom_geos_error);
1307 
1308  if (!input_lwgeom_to_geos(&g1, geom1, __func__)) return NULL;
1309  if (!input_lwgeom_to_geos(&g2, geom2, __func__)) return geos_clean_and_fail(g1, NULL, NULL, __func__);
1310 
1311  g3 = GEOSSharedPaths(g1, g2);
1312 
1313  if (!g3) return geos_clean_and_fail(g1, g2, NULL, __func__);
1314 
1315  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
1316  return geos_clean_and_fail(g1, g2, g3, __func__);
1317 
1318  geos_clean(g1, g2, g3);
1319  return result;
1320 }
1321 
1322 static LWGEOM *
1323 lwline_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
1324 {
1325  LWGEOM* result;
1326  LWGEOM* geom = lwline_as_lwgeom(lwline);
1327  int32_t srid = get_result_srid(geom, NULL, __func__);
1328  uint8_t is3d = FLAGS_GET_Z(geom->flags);
1329  GEOSGeometry *g1, *g3;
1330 
1331  if (srid == SRID_INVALID) return NULL;
1332 
1333  initGEOS(lwnotice, lwgeom_geos_error);
1334 
1335  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
1336 
1337  g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
1338 
1339  if (!g3)
1340  {
1341  geos_clean(g1, NULL, NULL);
1342  return NULL;
1343  }
1344 
1345  if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
1346  return geos_clean_and_fail(g1, NULL, g3, __func__);
1347 
1348  geos_clean(g1, NULL, g3);
1349  return result;
1350 }
1351 
1352 static LWGEOM *
1353 lwcollection_offsetcurve(const LWCOLLECTION *col, double size, int quadsegs, int joinStyle, double mitreLimit)
1354 {
1355  const LWGEOM *geom = lwcollection_as_lwgeom(col);
1356  int32_t srid = get_result_srid(geom, NULL, __func__);
1357  uint8_t is3d = FLAGS_GET_Z(col->flags);
1358  LWCOLLECTION *result;
1359  LWGEOM *tmp;
1360  uint32_t i;
1361  if (srid == SRID_INVALID) return NULL;
1362 
1363  result = lwcollection_construct_empty(MULTILINETYPE, srid, is3d, LW_FALSE);
1364 
1365  for (i = 0; i < col->ngeoms; i++)
1366  {
1367  tmp = lwgeom_offsetcurve(col->geoms[i], size, quadsegs, joinStyle, mitreLimit);
1368 
1369  if (!tmp)
1370  {
1371  lwcollection_free(result);
1372  return NULL;
1373  }
1374 
1375  if (!lwgeom_is_empty(tmp))
1376  {
1377  if (lwgeom_is_collection(tmp))
1379  else
1380  result = lwcollection_add_lwgeom(result, tmp);
1381 
1382  if (!result)
1383  {
1384  lwgeom_free(tmp);
1385  return NULL;
1386  }
1387  }
1388  }
1389 
1390  if (result->ngeoms == 1)
1391  {
1392  tmp = result->geoms[0];
1393  lwcollection_release(result);
1394  return tmp;
1395  }
1396  else
1397  return lwcollection_as_lwgeom(result);
1398 }
1399 
1400 LWGEOM*
1401 lwgeom_offsetcurve(const LWGEOM* geom, double size, int quadsegs, int joinStyle, double mitreLimit)
1402 {
1403  int32_t srid = get_result_srid(geom, NULL, __func__);
1404  LWGEOM *result = NULL;
1405  LWGEOM *noded = NULL;
1406  if (srid == SRID_INVALID) return NULL;
1407 
1408  if (lwgeom_dimension(geom) != 1)
1409  {
1410  lwerror("%s: input is not linear", __func__, lwtype_name(geom->type));
1411  return NULL;
1412  }
1413 
1414  while (!result)
1415  {
1416  switch (geom->type)
1417  {
1418  case LINETYPE:
1419  result = lwline_offsetcurve(lwgeom_as_lwline(geom), size, quadsegs, joinStyle, mitreLimit);
1420  break;
1421  case COLLECTIONTYPE:
1422  case MULTILINETYPE:
1423  result = lwcollection_offsetcurve(lwgeom_as_lwcollection(geom), size, quadsegs, joinStyle, mitreLimit);
1424  break;
1425  default:
1426  lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type));
1427  return NULL;
1428  }
1429 
1430  if (result)
1431  return result;
1432  else if (!noded)
1433  {
1434  noded = lwgeom_node(geom);
1435  if (!noded)
1436  {
1437  lwfree(noded);
1438  lwerror("lwgeom_offsetcurve: cannot node input");
1439  return NULL;
1440  }
1441  geom = noded;
1442  }
1443  else
1444  {
1445  lwerror("lwgeom_offsetcurve: noded geometry cannot be offset");
1446  return NULL;
1447  }
1448  }
1449  return result;
1450 }
1451 
1452 LWMPOINT*
1453 lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints)
1454 {
1455  double area, bbox_area, bbox_width, bbox_height;
1456  GBOX bbox;
1457  const LWGEOM* lwgeom = (LWGEOM*)lwpoly;
1458  uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
1459  double sample_cell_size;
1460  uint32_t i, j, n;
1461  uint32_t iterations = 0;
1462  uint32_t npoints_generated = 0;
1463  uint32_t npoints_tested = 0;
1464  GEOSGeometry* g;
1465  const GEOSPreparedGeometry* gprep;
1466  GEOSGeometry* gpt;
1467  GEOSCoordSequence* gseq;
1468  LWMPOINT* mpt;
1469  int srid = lwgeom_get_srid(lwgeom);
1470  int done = 0;
1471  int* cells;
1472  const size_t size = 2 * sizeof(int);
1473  char tmp[2 * sizeof(int)];
1474  const size_t stride = 2 * sizeof(int);
1475 
1476  if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
1477  {
1478  lwerror("%s: only polygons supported", __func__);
1479  return NULL;
1480  }
1481 
1482  if (npoints == 0 || lwgeom_is_empty(lwgeom)) return NULL;
1483 
1484  if (!lwpoly->bbox)
1485  lwgeom_calculate_gbox(lwgeom, &bbox);
1486  else
1487  bbox = *(lwpoly->bbox);
1488 
1489  area = lwpoly_area(lwpoly);
1490  bbox_width = bbox.xmax - bbox.xmin;
1491  bbox_height = bbox.ymax - bbox.ymin;
1492  bbox_area = bbox_width * bbox_height;
1493 
1494  if (area == 0.0 || bbox_area == 0.0)
1495  {
1496  lwerror("%s: zero area input polygon, TBD", __func__);
1497  return NULL;
1498  }
1499 
1500  /* Gross up our test set a bit to increase odds of getting coverage in one pass */
1501  sample_npoints = npoints * bbox_area / area;
1502 
1503  /* We're going to generate points using a sample grid as described
1504  * http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html to try and get a more uniform
1505  * "random" set of points. So we have to figure out how to stick a grid into our box */
1506  sample_sqrt = lround(sqrt(sample_npoints));
1507  if (sample_sqrt == 0) sample_sqrt = 1;
1508 
1509  /* Calculate the grids we're going to randomize within */
1510  if (bbox_width > bbox_height)
1511  {
1512  sample_width = sample_sqrt;
1513  sample_height = ceil((double)sample_npoints / (double)sample_width);
1514  sample_cell_size = bbox_width / sample_width;
1515  }
1516  else
1517  {
1518  sample_height = sample_sqrt;
1519  sample_width = ceil((double)sample_npoints / (double)sample_height);
1520  sample_cell_size = bbox_height / sample_height;
1521  }
1522 
1523  /* Prepare the polygon for fast true/false testing */
1524  initGEOS(lwnotice, lwgeom_geos_error);
1525  g = (GEOSGeometry*)LWGEOM2GEOS(lwgeom, 0);
1526  if (!g)
1527  {
1528  lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
1529  return NULL;
1530  }
1531  gprep = GEOSPrepare(g);
1532 
1533  /* Get an empty multi-point ready to return */
1534  mpt = lwmpoint_construct_empty(srid, 0, 0);
1535 
1536  /* Init random number generator */
1537  srand(time(NULL));
1538 
1539  /* Now we fill in an array of cells, and then shuffle that array, */
1540  /* so we can visit the cells in random order to avoid visual ugliness */
1541  /* caused by visiting them sequentially */
1542  cells = lwalloc(2 * sizeof(int) * sample_height * sample_width);
1543  for (i = 0; i < sample_width; i++)
1544  {
1545  for (j = 0; j < sample_height; j++)
1546  {
1547  cells[2 * (i * sample_height + j)] = i;
1548  cells[2 * (i * sample_height + j) + 1] = j;
1549  }
1550  }
1551 
1552  /* shuffle */
1553  n = sample_height * sample_width;
1554  if (n > 1)
1555  {
1556  for (i = 0; i < n - 1; ++i)
1557  {
1558  size_t rnd = (size_t)rand();
1559  size_t j = i + rnd / (RAND_MAX / (n - i) + 1);
1560 
1561  memcpy(tmp, (char *)cells + j * stride, size);
1562  memcpy((char *)cells + j * stride, (char *)cells + i * stride, size);
1563  memcpy((char *)cells + i * stride, tmp, size);
1564  }
1565  }
1566 
1567  /* Start testing points */
1568  while (npoints_generated < npoints)
1569  {
1570  iterations++;
1571  for (i = 0; i < sample_width * sample_height; i++)
1572  {
1573  int contains = 0;
1574  double y = bbox.ymin + cells[2 * i] * sample_cell_size;
1575  double x = bbox.xmin + cells[2 * i + 1] * sample_cell_size;
1576  x += rand() * sample_cell_size / RAND_MAX;
1577  y += rand() * sample_cell_size / RAND_MAX;
1578  if (x >= bbox.xmax || y >= bbox.ymax) continue;
1579 
1580  gseq = GEOSCoordSeq_create(1, 2);
1581  GEOSCoordSeq_setX(gseq, 0, x);
1582  GEOSCoordSeq_setY(gseq, 0, y);
1583  gpt = GEOSGeom_createPoint(gseq);
1584 
1585  contains = GEOSPreparedIntersects(gprep, gpt);
1586 
1587  GEOSGeom_destroy(gpt);
1588 
1589  if (contains == 2)
1590  {
1591  GEOSPreparedGeom_destroy(gprep);
1592  GEOSGeom_destroy(g);
1593  lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
1594  return NULL;
1595  }
1596  if (contains)
1597  {
1598  npoints_generated++;
1599  mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
1600  if (npoints_generated == npoints)
1601  {
1602  done = 1;
1603  break;
1604  }
1605  }
1606 
1607  /* Short-circuit check for ctrl-c occasionally */
1608  npoints_tested++;
1609  if (npoints_tested % 10000 == 0)
1610  LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
1611 
1612  if (done) break;
1613  }
1614  if (done || iterations > 100) break;
1615  }
1616 
1617  GEOSPreparedGeom_destroy(gprep);
1618  GEOSGeom_destroy(g);
1619  lwfree(cells);
1620 
1621  return mpt;
1622 }
1623 
1624 /* Allocate points to sub-geometries by area, then call lwgeom_poly_to_points and bundle up final result in a single
1625  * multipoint. */
1626 LWMPOINT*
1627 lwmpoly_to_points(const LWMPOLY* lwmpoly, uint32_t npoints)
1628 {
1629  const LWGEOM* lwgeom = (LWGEOM*)lwmpoly;
1630  double area;
1631  uint32_t i;
1632  LWMPOINT* mpt = NULL;
1633 
1634  if (lwgeom_get_type(lwgeom) != MULTIPOLYGONTYPE)
1635  {
1636  lwerror("%s: only multipolygons supported", __func__);
1637  return NULL;
1638  }
1639  if (npoints == 0 || lwgeom_is_empty(lwgeom)) return NULL;
1640 
1641  area = lwgeom_area(lwgeom);
1642 
1643  for (i = 0; i < lwmpoly->ngeoms; i++)
1644  {
1645  double sub_area = lwpoly_area(lwmpoly->geoms[i]);
1646  int sub_npoints = lround(npoints * sub_area / area);
1647  if (sub_npoints > 0)
1648  {
1649  LWMPOINT* sub_mpt = lwpoly_to_points(lwmpoly->geoms[i], sub_npoints);
1650  if (!mpt)
1651  mpt = sub_mpt;
1652  else
1653  {
1654  uint32_t j;
1655  for (j = 0; j < sub_mpt->ngeoms; j++)
1656  mpt = lwmpoint_add_lwpoint(mpt, sub_mpt->geoms[j]);
1657  /* Just free the shell, leave the underlying lwpoints alone, as they are now owned by
1658  * the returning multipoint */
1659  lwfree(sub_mpt->geoms);
1660  lwgeom_release((LWGEOM*)sub_mpt);
1661  }
1662  }
1663  }
1664  return mpt;
1665 }
1666 
1667 LWMPOINT*
1668 lwgeom_to_points(const LWGEOM* lwgeom, uint32_t npoints)
1669 {
1670  switch (lwgeom_get_type(lwgeom))
1671  {
1672  case MULTIPOLYGONTYPE:
1673  return lwmpoly_to_points((LWMPOLY*)lwgeom, npoints);
1674  case POLYGONTYPE:
1675  return lwpoly_to_points((LWPOLY*)lwgeom, npoints);
1676  default:
1677  lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(lwgeom_get_type(lwgeom)));
1678  return NULL;
1679  }
1680 }
1681 
1682 LWTIN*
1683 lwtin_from_geos(const GEOSGeometry* geom, uint8_t want3d)
1684 {
1685  int type = GEOSGeomTypeId(geom);
1686  int SRID = GEOSGetSRID(geom);
1687 
1688  /* GEOS's 0 is equivalent to our unknown as for SRID values */
1689  if (SRID == 0) SRID = SRID_UNKNOWN;
1690 
1691  if (want3d && !GEOSHasZ(geom))
1692  {
1693  LWDEBUG(3, "Geometry has no Z, won't provide one");
1694  want3d = 0;
1695  }
1696 
1697  switch (type)
1698  {
1699  LWTRIANGLE** geoms;
1700  uint32_t i, ngeoms;
1701  case GEOS_GEOMETRYCOLLECTION:
1702  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
1703 
1704  ngeoms = GEOSGetNumGeometries(geom);
1705  geoms = NULL;
1706  if (ngeoms)
1707  {
1708  geoms = lwalloc(ngeoms * sizeof *geoms);
1709  if (!geoms)
1710  {
1711  lwerror("lwtin_from_geos: can't allocate geoms");
1712  return NULL;
1713  }
1714  for (i = 0; i < ngeoms; i++)
1715  {
1716  const GEOSGeometry *poly, *ring;
1717  const GEOSCoordSequence* cs;
1718  POINTARRAY* pa;
1719 
1720  poly = GEOSGetGeometryN(geom, i);
1721  ring = GEOSGetExteriorRing(poly);
1722  cs = GEOSGeom_getCoordSeq(ring);
1723  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
1724 
1725  geoms[i] = lwtriangle_construct(SRID, NULL, pa);
1726  }
1727  }
1728  return (LWTIN*)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM**)geoms);
1729  case GEOS_POLYGON:
1730  case GEOS_MULTIPOINT:
1731  case GEOS_MULTILINESTRING:
1732  case GEOS_MULTIPOLYGON:
1733  case GEOS_LINESTRING:
1734  case GEOS_LINEARRING:
1735  case GEOS_POINT:
1736  lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
1737  break;
1738 
1739  default:
1740  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
1741  return NULL;
1742  }
1743 
1744  /* shouldn't get here */
1745  return NULL;
1746 }
1747 /*
1748  * output = 1 for edges, 2 for TIN, 0 for polygons
1749  */
1750 LWGEOM*
1751 lwgeom_delaunay_triangulation(const LWGEOM* geom, double tolerance, int32_t output)
1752 {
1753 #if POSTGIS_GEOS_VERSION < 34
1754  lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
1755  return NULL;
1756 #else
1757  LWGEOM* result;
1758  int32_t srid = get_result_srid(geom, NULL, __func__);
1759  uint8_t is3d = FLAGS_GET_Z(geom->flags);
1760  GEOSGeometry *g1, *g3;
1761 
1762  if (output < 0 || output > 2)
1763  {
1764  lwerror("%s: invalid output type specified %d", __func__, output);
1765  return NULL;
1766  }
1767 
1768  if (srid == SRID_INVALID) return NULL;
1769 
1770  initGEOS(lwnotice, lwgeom_geos_error);
1771 
1772  if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
1773 
1774  /* if output != 1 we want polys */
1775  g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
1776 
1777  if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
1778 
1779  if (output == 2)
1780  {
1781  result = (LWGEOM*)lwtin_from_geos(g3, is3d);
1782  if (!result)
1783  {
1784  geos_clean(g1, NULL, g3);
1785  lwerror("%s: cannot convert output geometry", __func__);
1786  return NULL;
1787  }
1788  lwgeom_set_srid(result, srid);
1789  }
1790  else if (!output_geos_as_lwgeom(&g3, &result, srid, is3d, __func__))
1791  return geos_clean_and_fail(g1, NULL, g3, __func__);
1792 
1793  geos_clean(g1, NULL, g3);
1794  return result;
1795 #endif /* POSTGIS_GEOS_VERSION < 34 */
1796 }
1797 
1798 #if POSTGIS_GEOS_VERSION < 35
1799 LWGEOM*
1800 lwgeom_voronoi_diagram(const LWGEOM* g, const GBOX* env, double tolerance, int output_edges)
1801 {
1802  lwerror("lwgeom_voronoi_diagram: GEOS 3.5 or higher required");
1803  return NULL;
1804 }
1805 #else /* POSTGIS_GEOS_VERSION >= 35 */
1806 static GEOSCoordSequence*
1808 {
1809  uint32_t i = 0;
1810  uint8_t num_dims = 2;
1811  LWPOINTITERATOR* it;
1812  GEOSCoordSequence* coords;
1813  POINT4D tmp;
1814 
1815  coords = GEOSCoordSeq_create(num_points, num_dims);
1816  if (!coords) return NULL;
1817 
1818  it = lwpointiterator_create(g);
1819  while (lwpointiterator_next(it, &tmp))
1820  {
1821  if (i >= num_points)
1822  {
1823  lwerror("Incorrect num_points provided to lwgeom_get_geos_coordseq_2d");
1824  GEOSCoordSeq_destroy(coords);
1826  return NULL;
1827  }
1828 
1829  if (!GEOSCoordSeq_setX(coords, i, tmp.x) || !GEOSCoordSeq_setY(coords, i, tmp.y))
1830  {
1831  GEOSCoordSeq_destroy(coords);
1833  return NULL;
1834  }
1835  i++;
1836  }
1838 
1839  return coords;
1840 }
1841 
1842 LWGEOM*
1843 lwgeom_voronoi_diagram(const LWGEOM* g, const GBOX* env, double tolerance, int output_edges)
1844 {
1845  uint32_t num_points = lwgeom_count_vertices(g);
1846  LWGEOM* lwgeom_result;
1847  char is_3d = LW_FALSE;
1848  int srid = lwgeom_get_srid(g);
1849  GEOSCoordSequence* coords;
1850  GEOSGeometry* geos_geom;
1851  GEOSGeometry* geos_env = NULL;
1852  GEOSGeometry* geos_result;
1853 
1854  if (num_points < 2)
1855  {
1857  return lwcollection_as_lwgeom(empty);
1858  }
1859 
1860  initGEOS(lwnotice, lwgeom_geos_error);
1861 
1862  /* Instead of using the standard LWGEOM2GEOS transformer, we read the vertices of the LWGEOM directly and put
1863  * them into a single GEOS CoordinateSeq that can be used to define a LineString. This allows us to process
1864  * geometry types that may not be supported by GEOS, and reduces the memory requirements in cases of many
1865  * geometries with few points (such as LWMPOINT).*/
1866  coords = lwgeom_get_geos_coordseq_2d(g, num_points);
1867  if (!coords) return NULL;
1868 
1869  geos_geom = GEOSGeom_createLineString(coords);
1870  if (!geos_geom)
1871  {
1872  GEOSCoordSeq_destroy(coords);
1873  return NULL;
1874  }
1875 
1876  if (env) geos_env = GBOX2GEOS(env);
1877 
1878  geos_result = GEOSVoronoiDiagram(geos_geom, geos_env, tolerance, output_edges);
1879 
1880  GEOSGeom_destroy(geos_geom);
1881  if (env) GEOSGeom_destroy(geos_env);
1882 
1883  if (!geos_result)
1884  {
1885  lwerror("GEOSVoronoiDiagram: %s", lwgeom_geos_errmsg);
1886  return NULL;
1887  }
1888 
1889  lwgeom_result = GEOS2LWGEOM(geos_result, is_3d);
1890  GEOSGeom_destroy(geos_result);
1891 
1892  lwgeom_set_srid(lwgeom_result, srid);
1893 
1894  return lwgeom_result;
1895 }
1896 #endif /* POSTGIS_GEOS_VERSION >= 35 */
LWGEOM * lwgeom_unaryunion(const LWGEOM *geom)
double x
Definition: liblwgeom.h:354
#define LINETYPE
Definition: liblwgeom.h:85
POINTARRAY * ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, uint8_t want3d)
LWGEOM * lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)
Definition: lwstroke.c:715
LWGEOM * lwgeom_voronoi_diagram(const LWGEOM *g, const GBOX *env, double tolerance, int output_edges)
Take vertices of a geometry and build the Voronoi diagram.
double z
Definition: liblwgeom.h:336
double y
Definition: liblwgeom.h:336
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
int lwgeom_is_simple(const LWGEOM *geom)
LWLINE * lwline_construct_empty(int srid, char hasz, char hasm)
Definition: lwline.c:64
LWGEOM * lwgeom_clip_by_rect(const LWGEOM *geom1, double x1, double y1, double x2, double y2)
GEOSGeometry * make_geos_segment(double x1, double y1, double x2, double y2)
uint32_t ngeoms
Definition: liblwgeom.h:470
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1085
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, uint8_t want3d)
int lwpointiterator_next(LWPOINTITERATOR *s, POINT4D *p)
Attempts to assign the next point in the iterator to p, and advances the iterator to the next point...
Definition: lwiterator.c:212
LWCOLLECTION * lwcollection_extract(LWCOLLECTION *col, int type)
Takes a potentially heterogeneous collection and returns a homogeneous collection consisting only of ...
Definition: lwcollection.c:386
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:43
double x
Definition: liblwgeom.h:336
static void delFace(Face *f)
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:923
uint32_t ngeoms
Definition: liblwgeom.h:496
LWGEOM * lwgeom_linemerge(const LWGEOM *geom)
static int compare_by_envarea(const void *g1, const void *g2)
#define POLYGONTYPE
Definition: liblwgeom.h:86
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
Datum area(PG_FUNCTION_ARGS)
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:556
uint8_t flags
Definition: liblwgeom.h:399
double xmax
Definition: liblwgeom.h:295
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
static LWGEOM * geos_clean_and_fail(GEOSGeometry *g1, GEOSGeometry *g2, GEOSGeometry *g3, const char *funcname)
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
def fmt
Definition: pixval.py:92
LWGEOM * lwgeom_symdifference(const LWGEOM *geom1, const LWGEOM *geom2)
GEOSGeometry * LWGEOM_GEOS_makeValid(const GEOSGeometry *)
Datum contains(PG_FUNCTION_ARGS)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:425
static GEOSGeometry * ptarray_to_GEOSLinearRing(const POINTARRAY *pa, uint8_t autofix)
static LWGEOM * lwline_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
GEOSGeometry * GBOX2GEOS(const GBOX *box)
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:916
#define LW_ON_INTERRUPT(x)
struct Face_t Face
static void findFaceHoles(Face **faces, int nfaces)
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:520
GBOX * bbox
Definition: liblwgeom.h:455
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:151
LWGEOM * lwgeom_centroid(const LWGEOM *geom)
uint32_t ngeoms
Definition: liblwgeom.h:509
POINTARRAY * point
Definition: liblwgeom.h:413
int ptarray_is_closed_2d(const POINTARRAY *pa)
Definition: ptarray.c:688
static LWGEOM * lwcollection_offsetcurve(const LWCOLLECTION *col, double size, int quadsegs, int joinStyle, double mitreLimit)
LWPOINTITERATOR * lwpointiterator_create(const LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM*.
Definition: lwiterator.c:244
int32_t srid
Definition: liblwgeom.h:401
GEOSGeometry * make_geos_point(double x, double y)
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
static GEOSCoordSequence * lwgeom_get_geos_coordseq_2d(const LWGEOM *g, uint32_t num_points)
LWCOLLECTION * lwcollection_concat_in_place(LWCOLLECTION *col1, const LWCOLLECTION *col2)
Appends all geometries from col2 to col1 in place.
Definition: lwcollection.c:240
uint32_t nrings
Definition: liblwgeom.h:457
static uint8_t output_geos_as_lwgeom(GEOSGeometry **g, LWGEOM **geom, const int32_t srid, const uint8_t is3d, const char *funcname)
LWTIN * lwtin_from_geos(const GEOSGeometry *geom, uint8_t want3d)
void lwgeom_simplify_in_place(LWGEOM *igeom, double dist, int preserve_collapsed)
Definition: lwgeom.c:1748
unsigned int uint32_t
Definition: uthash.h:78
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:746
double x
Definition: liblwgeom.h:330
uint8_t flags
Definition: liblwgeom.h:506
#define SRID_USER_MAXIMUM
Maximum valid SRID value for the user We reserve 1000 values for internal use.
Definition: liblwgeom.h:184
LWMPOINT * lwpoly_to_points(const LWPOLY *lwpoly, uint32_t npoints)
LWPOLY * lwpoly_construct_envelope(int srid, double x1, double y1, double x2, double y2)
Definition: lwpoly.c:98
LWGEOM * lwgeom_buildarea(const LWGEOM *geom)
Take a geometry and return an areal geometry (Polygon or MultiPolygon).
POINTARRAY * ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
Add a point in a pointarray.
Definition: ptarray.c:503
const POINT3DZ * getPoint3dz_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT3DZ pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:378
double ymin
Definition: liblwgeom.h:296
void lwpointiterator_destroy(LWPOINTITERATOR *s)
Free all memory associated with the iterator.
Definition: lwiterator.c:269
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
LWGEOM * geom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:330
double xmin
Definition: liblwgeom.h:294
void lwgeom_geos_error(const char *fmt,...)
LWGEOM * lwgeom_node(const LWGEOM *lwgeom_in)
#define LW_FALSE
Definition: liblwgeom.h:76
LWGEOM * lwgeom_sharedpaths(const LWGEOM *geom1, const LWGEOM *geom2)
static void geos_clean(GEOSGeometry *g1, GEOSGeometry *g2, GEOSGeometry *g3)
uint8_t flags
Definition: liblwgeom.h:371
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:43
LWMPOINT * lwgeom_to_points(const LWGEOM *lwgeom, uint32_t npoints)
LWPOLY ** geoms
Definition: liblwgeom.h:498
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
LWGEOM ** geoms
Definition: liblwgeom.h:511
static GEOSGeometry * collectFacesWithEvenAncestors(Face **faces, int nfaces)
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
#define TINTYPE
Definition: liblwgeom.h:98
const GEOSGeometry * geom
LWGEOM * lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
static unsigned int countParens(const Face *f)
LWGEOM * lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
#define SRID_INVALID
Definition: liblwgeom.h:191
POINTARRAY ** rings
Definition: liblwgeom.h:459
LWPOINT ** geoms
Definition: liblwgeom.h:472
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
Definition: lwpoly.c:441
LWMPOINT * lwmpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoint.c:39
LWGEOM * lwgeom_delaunay_triangulation(const LWGEOM *geom, double tolerance, int32_t output)
Take vertices of a geometry and build a delaunay triangulation on them.
double ymax
Definition: liblwgeom.h:297
int lwgeom_has_arc(const LWGEOM *geom)
Definition: lwstroke.c:55
static Face * newFace(const GEOSGeometry *g)
double y
Definition: liblwgeom.h:330
#define LWGEOM_GEOS_ERRMSG_MAXSIZE
LWGEOM * lwgeom_geos_noop(const LWGEOM *geom)
Convert an LWGEOM to a GEOS Geometry and convert back – for debug only.
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:139
double z
Definition: liblwgeom.h:354
LWTRIANGLE * lwtriangle_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwtriangle.c:40
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
GEOSGeometry * env
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1235
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
LWMPOINT * lwmpoly_to_points(const LWMPOLY *lwmpoly, uint32_t npoints)
void lwcollection_release(LWCOLLECTION *lwcollection)
Definition: lwcollection.c:36
static uint8_t input_lwgeom_to_geos(GEOSGeometry **g, const LWGEOM *geom, const char *funcname)
LWGEOM * lwgeom_snap(const LWGEOM *geom1, const LWGEOM *geom2, double tolerance)
Snap vertices and segments of a geometry to another using a given tolerance.
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:224
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1857
int lwgeom_dimension(const LWGEOM *geom)
For an LWGEOM, returns 0 for points, 1 for lines, 2 for polygons, 3 for volume, and the max dimension...
Definition: lwgeom.c:1287
uint8_t * getPoint_internal(const POINTARRAY *pa, uint32_t n)
Definition: ptarray.c:1743
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:161
void lwgeom_release(LWGEOM *lwgeom)
Free the containing LWGEOM and the associated BOX.
Definition: lwgeom.c:459
uint8_t type
Definition: liblwgeom.h:398
type
Definition: ovdump.py:41
struct Face_t * parent
void free(void *)
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:356
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:335
const char * lwgeom_geos_version()
Return GEOS version string (not to be freed)
void * malloc(YYSIZE_T)
void lwgeom_set_srid(LWGEOM *geom, int srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
static int32_t get_result_srid(const LWGEOM *geom1, const LWGEOM *geom2, const char *funcname)
#define SRID_MAXIMUM
Maximum allowed SRID value in serialized geometry.
Definition: liblwgeom.h:178
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:129
void * lwalloc(size_t size)
Definition: lwutil.c:229
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1393
double y
Definition: liblwgeom.h:354
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
LWGEOM * lwgeom_normalize(const LWGEOM *geom)
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
unsigned char uint8_t
Definition: uthash.h:79
LWGEOM * lwgeom_offsetcurve(const LWGEOM *geom, double size, int quadsegs, int joinStyle, double mitreLimit)
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:937
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
This library is the generic geometry handling section of PostGIS.
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *, uint8_t fix_ring)
POINTARRAY * points
Definition: liblwgeom.h:424
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:300
uint32_t npoints
Definition: liblwgeom.h:373