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