PostGIS  2.2.8dev-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  * Copyright 2011-2014 Sandro Santilli <strk@keybit.net>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "lwgeom_geos.h"
14 #include "liblwgeom.h"
15 #include "liblwgeom_internal.h"
16 #include "lwgeom_log.h"
17 
18 #include <stdlib.h>
19 
20 LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
21 
22 #undef LWGEOM_PROFILE_BUILDAREA
23 
24 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
26 
27 extern void
28 lwgeom_geos_error(const char *fmt, ...)
29 {
30  va_list ap;
31 
32  va_start(ap, fmt);
33 
34  /* Call the supplied function */
36  {
38  }
39 
40  va_end(ap);
41 }
42 
43 
44 /*
45 ** GEOS <==> PostGIS conversion functions
46 **
47 ** Default conversion creates a GEOS point array, then iterates through the
48 ** PostGIS points, setting each value in the GEOS array one at a time.
49 **
50 */
51 
52 /* Return a POINTARRAY from a GEOSCoordSeq */
53 POINTARRAY *
54 ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, char want3d)
55 {
56  uint32_t dims=2;
57  uint32_t size, i;
58  POINTARRAY *pa;
59  POINT4D point;
60 
61  LWDEBUG(2, "ptarray_fromGEOSCoordSeq called");
62 
63  if ( ! GEOSCoordSeq_getSize(cs, &size) )
64  lwerror("Exception thrown");
65 
66  LWDEBUGF(4, " GEOSCoordSeq size: %d", size);
67 
68  if ( want3d )
69  {
70  if ( ! GEOSCoordSeq_getDimensions(cs, &dims) )
71  lwerror("Exception thrown");
72 
73  LWDEBUGF(4, " GEOSCoordSeq dimensions: %d", dims);
74 
75  /* forget higher dimensions (if any) */
76  if ( dims > 3 ) dims = 3;
77  }
78 
79  LWDEBUGF(4, " output dimensions: %d", dims);
80 
81  pa = ptarray_construct((dims==3), 0, size);
82 
83  for (i=0; i<size; i++)
84  {
85  GEOSCoordSeq_getX(cs, i, &(point.x));
86  GEOSCoordSeq_getY(cs, i, &(point.y));
87  if ( dims >= 3 ) GEOSCoordSeq_getZ(cs, i, &(point.z));
88  ptarray_set_point4d(pa,i,&point);
89  }
90 
91  return pa;
92 }
93 
94 /* Return an LWGEOM from a Geometry */
95 LWGEOM *
96 GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
97 {
98  int type = GEOSGeomTypeId(geom) ;
99  int hasZ;
100  int SRID = GEOSGetSRID(geom);
101 
102  /* GEOS's 0 is equivalent to our unknown as for SRID values */
103  if ( SRID == 0 ) SRID = SRID_UNKNOWN;
104 
105  if ( want3d )
106  {
107  hasZ = GEOSHasZ(geom);
108  if ( ! hasZ )
109  {
110  LWDEBUG(3, "Geometry has no Z, won't provide one");
111 
112  want3d = 0;
113  }
114  }
115 
116 /*
117  if ( GEOSisEmpty(geom) )
118  {
119  return (LWGEOM*)lwcollection_construct_empty(COLLECTIONTYPE, SRID, want3d, 0);
120  }
121 */
122 
123  switch (type)
124  {
125  const GEOSCoordSequence *cs;
126  POINTARRAY *pa, **ppaa;
127  const GEOSGeometry *g;
128  LWGEOM **geoms;
129  uint32_t i, ngeoms;
130 
131  case GEOS_POINT:
132  LWDEBUG(4, "lwgeom_from_geometry: it's a Point");
133  cs = GEOSGeom_getCoordSeq(geom);
134  if ( GEOSisEmpty(geom) )
135  return (LWGEOM*)lwpoint_construct_empty(SRID, want3d, 0);
136  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
137  return (LWGEOM *)lwpoint_construct(SRID, NULL, pa);
138 
139  case GEOS_LINESTRING:
140  case GEOS_LINEARRING:
141  LWDEBUG(4, "lwgeom_from_geometry: it's a LineString or LinearRing");
142  if ( GEOSisEmpty(geom) )
143  return (LWGEOM*)lwline_construct_empty(SRID, want3d, 0);
144 
145  cs = GEOSGeom_getCoordSeq(geom);
146  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
147  return (LWGEOM *)lwline_construct(SRID, NULL, pa);
148 
149  case GEOS_POLYGON:
150  LWDEBUG(4, "lwgeom_from_geometry: it's a Polygon");
151  if ( GEOSisEmpty(geom) )
152  return (LWGEOM*)lwpoly_construct_empty(SRID, want3d, 0);
153  ngeoms = GEOSGetNumInteriorRings(geom);
154  ppaa = lwalloc(sizeof(POINTARRAY *)*(ngeoms+1));
155  g = GEOSGetExteriorRing(geom);
156  cs = GEOSGeom_getCoordSeq(g);
157  ppaa[0] = ptarray_from_GEOSCoordSeq(cs, want3d);
158  for (i=0; i<ngeoms; i++)
159  {
160  g = GEOSGetInteriorRingN(geom, i);
161  cs = GEOSGeom_getCoordSeq(g);
162  ppaa[i+1] = ptarray_from_GEOSCoordSeq(cs,
163  want3d);
164  }
165  return (LWGEOM *)lwpoly_construct(SRID, NULL,
166  ngeoms+1, ppaa);
167 
168  case GEOS_MULTIPOINT:
169  case GEOS_MULTILINESTRING:
170  case GEOS_MULTIPOLYGON:
171  case GEOS_GEOMETRYCOLLECTION:
172  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
173 
174  ngeoms = GEOSGetNumGeometries(geom);
175  geoms = NULL;
176  if ( ngeoms )
177  {
178  geoms = lwalloc(sizeof(LWGEOM *)*ngeoms);
179  for (i=0; i<ngeoms; i++)
180  {
181  g = GEOSGetGeometryN(geom, i);
182  geoms[i] = GEOS2LWGEOM(g, want3d);
183  }
184  }
185  return (LWGEOM *)lwcollection_construct(type,
186  SRID, NULL, ngeoms, geoms);
187 
188  default:
189  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
190  return NULL;
191 
192  }
193 
194 }
195 
196 
197 
198 GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *);
199 
200 
201 GEOSCoordSeq
203 {
204  uint32_t dims = 2;
205  uint32_t i;
206  const POINT3DZ *p3d;
207  const POINT2D *p2d;
208  GEOSCoordSeq sq;
209 
210  if ( FLAGS_GET_Z(pa->flags) )
211  dims = 3;
212 
213  if ( ! (sq = GEOSCoordSeq_create(pa->npoints, dims)) )
214  lwerror("Error creating GEOS Coordinate Sequence");
215 
216  for ( i=0; i < pa->npoints; i++ )
217  {
218  if ( dims == 3 )
219  {
220  p3d = getPoint3dz_cp(pa, i);
221  p2d = (const POINT2D *)p3d;
222  LWDEBUGF(4, "Point: %g,%g,%g", p3d->x, p3d->y, p3d->z);
223  }
224  else
225  {
226  p2d = getPoint2d_cp(pa, i);
227  LWDEBUGF(4, "Point: %g,%g", p2d->x, p2d->y);
228  }
229 
230 #if POSTGIS_GEOS_VERSION < 33
231  /* Make sure we don't pass any infinite values down into GEOS */
232  /* GEOS 3.3+ is supposed to handle this stuff OK */
233  if ( isinf(p2d->x) || isinf(p2d->y) || (dims == 3 && isinf(p3d->z)) )
234  lwerror("Infinite coordinate value found in geometry.");
235  if ( isnan(p2d->x) || isnan(p2d->y) || (dims == 3 && isnan(p3d->z)) )
236  lwerror("NaN coordinate value found in geometry.");
237 #endif
238 
239  GEOSCoordSeq_setX(sq, i, p2d->x);
240  GEOSCoordSeq_setY(sq, i, p2d->y);
241 
242  if ( dims == 3 )
243  GEOSCoordSeq_setZ(sq, i, p3d->z);
244  }
245  return sq;
246 }
247 
248 static GEOSGeometry *
249 ptarray_to_GEOSLinearRing(const POINTARRAY *pa, int autofix)
250 {
251  GEOSCoordSeq sq;
252  GEOSGeom g;
253  POINTARRAY *npa = 0;
254 
255  if ( autofix )
256  {
257  /* check ring for being closed and fix if not */
258  if ( ! ptarray_is_closed_2d(pa) )
259  {
260  npa = ptarray_addPoint(pa, getPoint_internal(pa, 0), FLAGS_NDIMS(pa->flags), pa->npoints);
261  pa = npa;
262  }
263  /* TODO: check ring for having at least 4 vertices */
264 #if 0
265  while ( pa->npoints < 4 )
266  {
267  npa = ptarray_addPoint(npa, getPoint_internal(pa, 0), FLAGS_NDIMS(pa->flags), pa->npoints);
268  }
269 #endif
270  }
271 
272  sq = ptarray_to_GEOSCoordSeq(pa);
273  if ( npa ) ptarray_free(npa);
274  g = GEOSGeom_createLinearRing(sq);
275  return g;
276 }
277 
278 GEOSGeometry *
279 GBOX2GEOS(const GBOX *box)
280 {
281  GEOSGeometry* envelope;
282  GEOSGeometry* ring;
283  GEOSCoordSequence* seq = GEOSCoordSeq_create(5, 2);
284  if (!seq)
285  {
286  return NULL;
287  }
288 
289  GEOSCoordSeq_setX(seq, 0, box->xmin);
290  GEOSCoordSeq_setY(seq, 0, box->ymin);
291 
292  GEOSCoordSeq_setX(seq, 1, box->xmax);
293  GEOSCoordSeq_setY(seq, 1, box->ymin);
294 
295  GEOSCoordSeq_setX(seq, 2, box->xmax);
296  GEOSCoordSeq_setY(seq, 2, box->ymax);
297 
298  GEOSCoordSeq_setX(seq, 3, box->xmin);
299  GEOSCoordSeq_setY(seq, 3, box->ymax);
300 
301  GEOSCoordSeq_setX(seq, 4, box->xmin);
302  GEOSCoordSeq_setY(seq, 4, box->ymin);
303 
304  ring = GEOSGeom_createLinearRing(seq);
305  if (!ring)
306  {
307  GEOSCoordSeq_destroy(seq);
308  return NULL;
309  }
310 
311  envelope = GEOSGeom_createPolygon(ring, NULL, 0);
312  if (!envelope)
313  {
314  GEOSGeom_destroy(ring);
315  return NULL;
316  }
317 
318  return envelope;
319 }
320 
321 GEOSGeometry *
322 LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
323 {
324  GEOSCoordSeq sq;
325  GEOSGeom g, shell;
326  GEOSGeom *geoms = NULL;
327  /*
328  LWGEOM *tmp;
329  */
330  uint32_t ngeoms, i, j;
331  int geostype;
332 #if LWDEBUG_LEVEL >= 4
333  char *wkt;
334 #endif
335 
336  LWDEBUGF(4, "LWGEOM2GEOS got a %s", lwtype_name(lwgeom->type));
337 
338  if (lwgeom_has_arc(lwgeom))
339  {
340  LWGEOM *lwgeom_stroked = lwgeom_stroke(lwgeom, 32);
341  GEOSGeometry *g = LWGEOM2GEOS(lwgeom_stroked, autofix);
342  lwgeom_free(lwgeom_stroked);
343  return g;
344  }
345 
346  switch (lwgeom->type)
347  {
348  LWPOINT *lwp = NULL;
349  LWPOLY *lwpoly = NULL;
350  LWLINE *lwl = NULL;
351  LWCOLLECTION *lwc = NULL;
352 #if POSTGIS_GEOS_VERSION < 33
353  POINTARRAY *pa = NULL;
354 #endif
355 
356  case POINTTYPE:
357  lwp = (LWPOINT *)lwgeom;
358 
359  if ( lwgeom_is_empty(lwgeom) )
360  {
361 #if POSTGIS_GEOS_VERSION < 33
362  pa = ptarray_construct_empty(lwgeom_has_z(lwgeom), lwgeom_has_m(lwgeom), 2);
363  sq = ptarray_to_GEOSCoordSeq(pa);
364  shell = GEOSGeom_createLinearRing(sq);
365  g = GEOSGeom_createPolygon(shell, NULL, 0);
366 #else
367  g = GEOSGeom_createEmptyPolygon();
368 #endif
369  }
370  else
371  {
372  sq = ptarray_to_GEOSCoordSeq(lwp->point);
373  g = GEOSGeom_createPoint(sq);
374  }
375  if ( ! g )
376  {
377  /* lwnotice("Exception in LWGEOM2GEOS"); */
378  return NULL;
379  }
380  break;
381  case LINETYPE:
382  lwl = (LWLINE *)lwgeom;
383  /* TODO: if (autofix) */
384  if ( lwl->points->npoints == 1 ) {
385  /* Duplicate point, to make geos-friendly */
386  lwl->points = ptarray_addPoint(lwl->points,
387  getPoint_internal(lwl->points, 0),
388  FLAGS_NDIMS(lwl->points->flags),
389  lwl->points->npoints);
390  }
391  sq = ptarray_to_GEOSCoordSeq(lwl->points);
392  g = GEOSGeom_createLineString(sq);
393  if ( ! g )
394  {
395  /* lwnotice("Exception in LWGEOM2GEOS"); */
396  return NULL;
397  }
398  break;
399 
400  case POLYGONTYPE:
401  lwpoly = (LWPOLY *)lwgeom;
402  if ( lwgeom_is_empty(lwgeom) )
403  {
404 #if POSTGIS_GEOS_VERSION < 33
405  POINTARRAY *pa = ptarray_construct_empty(lwgeom_has_z(lwgeom), lwgeom_has_m(lwgeom), 2);
406  sq = ptarray_to_GEOSCoordSeq(pa);
407  shell = GEOSGeom_createLinearRing(sq);
408  g = GEOSGeom_createPolygon(shell, NULL, 0);
409 #else
410  g = GEOSGeom_createEmptyPolygon();
411 #endif
412  }
413  else
414  {
415  shell = ptarray_to_GEOSLinearRing(lwpoly->rings[0], autofix);
416  if ( ! shell ) return NULL;
417  /*lwerror("LWGEOM2GEOS: exception during polygon shell conversion"); */
418  ngeoms = lwpoly->nrings-1;
419  if ( ngeoms > 0 )
420  geoms = malloc(sizeof(GEOSGeom)*ngeoms);
421 
422  for (i=1; i<lwpoly->nrings; ++i)
423  {
424  geoms[i-1] = ptarray_to_GEOSLinearRing(lwpoly->rings[i], autofix);
425  if ( ! geoms[i-1] )
426  {
427  --i;
428  while (i) GEOSGeom_destroy(geoms[--i]);
429  free(geoms);
430  GEOSGeom_destroy(shell);
431  return NULL;
432  }
433  /*lwerror("LWGEOM2GEOS: exception during polygon hole conversion"); */
434  }
435  g = GEOSGeom_createPolygon(shell, geoms, ngeoms);
436  if (geoms) free(geoms);
437  }
438  if ( ! g ) return NULL;
439  break;
440  case MULTIPOINTTYPE:
441  case MULTILINETYPE:
442  case MULTIPOLYGONTYPE:
443  case COLLECTIONTYPE:
444  if ( lwgeom->type == MULTIPOINTTYPE )
445  geostype = GEOS_MULTIPOINT;
446  else if ( lwgeom->type == MULTILINETYPE )
447  geostype = GEOS_MULTILINESTRING;
448  else if ( lwgeom->type == MULTIPOLYGONTYPE )
449  geostype = GEOS_MULTIPOLYGON;
450  else
451  geostype = GEOS_GEOMETRYCOLLECTION;
452 
453  lwc = (LWCOLLECTION *)lwgeom;
454 
455  ngeoms = lwc->ngeoms;
456  if ( ngeoms > 0 )
457  geoms = malloc(sizeof(GEOSGeom)*ngeoms);
458 
459  j = 0;
460  for (i=0; i<ngeoms; ++i)
461  {
462  GEOSGeometry* g;
463 
464  if( lwgeom_is_empty(lwc->geoms[i]) )
465  continue;
466 
467  g = LWGEOM2GEOS(lwc->geoms[i], 0);
468  if ( ! g )
469  {
470  while (j) GEOSGeom_destroy(geoms[--j]);
471  free(geoms);
472  return NULL;
473  }
474  geoms[j++] = g;
475  }
476  g = GEOSGeom_createCollection(geostype, geoms, j);
477  if ( geoms ) free(geoms);
478  if ( ! g ) return NULL;
479  break;
480 
481  default:
482  lwerror("Unknown geometry type: %d - %s", lwgeom->type, lwtype_name(lwgeom->type));
483  return NULL;
484  }
485 
486  GEOSSetSRID(g, lwgeom->srid);
487 
488 #if LWDEBUG_LEVEL >= 4
489  wkt = GEOSGeomToWKT(g);
490  LWDEBUGF(4, "LWGEOM2GEOS: GEOSGeom: %s", wkt);
491  free(wkt);
492 #endif
493 
494  return g;
495 }
496 
497 const char*
499 {
500  const char *ver = GEOSversion();
501  return ver;
502 }
503 
504 LWGEOM *
506 {
507  LWGEOM *result ;
508  GEOSGeometry *g1;
509  int is3d ;
510  int srid ;
511 
512  srid = (int)(geom1->srid);
513  is3d = FLAGS_GET_Z(geom1->flags);
514 
515  initGEOS(lwnotice, lwgeom_geos_error);
516 
517  g1 = LWGEOM2GEOS(geom1, 0);
518  if ( 0 == g1 ) /* exception thrown at construction */
519  {
520  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
521  return NULL ;
522  }
523 
524  if ( -1 == GEOSNormalize(g1) )
525  {
526  lwerror("Error in GEOSNormalize: %s", lwgeom_geos_errmsg);
527  return NULL; /* never get here */
528  }
529 
530  GEOSSetSRID(g1, srid); /* needed ? */
531  result = GEOS2LWGEOM(g1, is3d);
532  GEOSGeom_destroy(g1);
533 
534  if (result == NULL)
535  {
536  lwerror("Error performing intersection: GEOS2LWGEOM: %s",
538  return NULL ; /* never get here */
539  }
540 
541  return result ;
542 }
543 
544 LWGEOM *
545 lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
546 {
547  LWGEOM *result ;
548  GEOSGeometry *g1, *g2, *g3 ;
549  int is3d ;
550  int srid ;
551 
552  /* A.Intersection(Empty) == Empty */
553  if ( lwgeom_is_empty(geom2) )
554  return lwgeom_clone_deep(geom2);
555 
556  /* Empty.Intersection(A) == Empty */
557  if ( lwgeom_is_empty(geom1) )
558  return lwgeom_clone_deep(geom1);
559 
560  /* ensure srids are identical */
561  srid = (int)(geom1->srid);
562  error_if_srid_mismatch(srid, (int)(geom2->srid));
563 
564  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
565 
566  initGEOS(lwnotice, lwgeom_geos_error);
567 
568  LWDEBUG(3, "intersection() START");
569 
570  g1 = LWGEOM2GEOS(geom1, 0);
571  if ( 0 == g1 ) /* exception thrown at construction */
572  {
573  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
574  return NULL ;
575  }
576 
577  g2 = LWGEOM2GEOS(geom2, 0);
578  if ( 0 == g2 ) /* exception thrown at construction */
579  {
580  lwerror("Second argument geometry could not be converted to GEOS.");
581  GEOSGeom_destroy(g1);
582  return NULL ;
583  }
584 
585  LWDEBUG(3, " constructed geometrys - calling geos");
586  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
587  LWDEBUGF(3, " g2 = %s", GEOSGeomToWKT(g2));
588  /*LWDEBUGF(3, "g2 is valid = %i",GEOSisvalid(g2)); */
589  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
590 
591  g3 = GEOSIntersection(g1,g2);
592 
593  LWDEBUG(3, " intersection finished");
594 
595  if (g3 == NULL)
596  {
597  GEOSGeom_destroy(g1);
598  GEOSGeom_destroy(g2);
599  lwerror("Error performing intersection: %s",
601  return NULL; /* never get here */
602  }
603 
604  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
605 
606  GEOSSetSRID(g3, srid);
607 
608  result = GEOS2LWGEOM(g3, is3d);
609 
610  if (result == NULL)
611  {
612  GEOSGeom_destroy(g1);
613  GEOSGeom_destroy(g2);
614  GEOSGeom_destroy(g3);
615  lwerror("Error performing intersection: GEOS2LWGEOM: %s",
617  return NULL ; /* never get here */
618  }
619 
620  GEOSGeom_destroy(g1);
621  GEOSGeom_destroy(g2);
622  GEOSGeom_destroy(g3);
623 
624  return result ;
625 }
626 
627 LWGEOM *
629 {
630  LWGEOM *result ;
631  GEOSGeometry *g1, *g3 ;
632  int is3d = FLAGS_GET_Z(geom1->flags);
633  int srid = geom1->srid;
634 
635  /* Empty.Linemerge() == Empty */
636  if ( lwgeom_is_empty(geom1) )
637  return (LWGEOM*)lwcollection_construct_empty( COLLECTIONTYPE, srid, is3d,
638  lwgeom_has_m(geom1) );
639 
640  initGEOS(lwnotice, lwgeom_geos_error);
641 
642  LWDEBUG(3, "linemerge() START");
643 
644  g1 = LWGEOM2GEOS(geom1, 0);
645  if ( 0 == g1 ) /* exception thrown at construction */
646  {
647  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
648  return NULL ;
649  }
650 
651  LWDEBUG(3, " constructed geometrys - calling geos");
652  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
653  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
654 
655  g3 = GEOSLineMerge(g1);
656 
657  LWDEBUG(3, " linemerge finished");
658 
659  if (g3 == NULL)
660  {
661  GEOSGeom_destroy(g1);
662  lwerror("Error performing linemerge: %s",
664  return NULL; /* never get here */
665  }
666 
667  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
668 
669  GEOSSetSRID(g3, srid);
670 
671  result = GEOS2LWGEOM(g3, is3d);
672 
673  if (result == NULL)
674  {
675  GEOSGeom_destroy(g1);
676  GEOSGeom_destroy(g3);
677  lwerror("Error performing linemerge: GEOS2LWGEOM: %s",
679  return NULL ; /* never get here */
680  }
681 
682  GEOSGeom_destroy(g1);
683  GEOSGeom_destroy(g3);
684 
685  return result ;
686 }
687 
688 LWGEOM *
690 {
691  LWGEOM *result ;
692  GEOSGeometry *g1, *g3 ;
693  int is3d = FLAGS_GET_Z(geom1->flags);
694  int srid = geom1->srid;
695 
696  /* Empty.UnaryUnion() == Empty */
697  if ( lwgeom_is_empty(geom1) )
698  return lwgeom_clone_deep(geom1);
699 
700  initGEOS(lwnotice, lwgeom_geos_error);
701 
702  g1 = LWGEOM2GEOS(geom1, 0);
703  if ( 0 == g1 ) /* exception thrown at construction */
704  {
705  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
706  return NULL ;
707  }
708 
709  g3 = GEOSUnaryUnion(g1);
710 
711  if (g3 == NULL)
712  {
713  GEOSGeom_destroy(g1);
714  lwerror("Error performing unaryunion: %s",
716  return NULL; /* never get here */
717  }
718 
719  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
720 
721  GEOSSetSRID(g3, srid);
722 
723  result = GEOS2LWGEOM(g3, is3d);
724 
725  if (result == NULL)
726  {
727  GEOSGeom_destroy(g1);
728  GEOSGeom_destroy(g3);
729  lwerror("Error performing unaryunion: GEOS2LWGEOM: %s",
731  return NULL ; /* never get here */
732  }
733 
734  GEOSGeom_destroy(g1);
735  GEOSGeom_destroy(g3);
736 
737  return result ;
738 }
739 
740 LWGEOM *
741 lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
742 {
743  GEOSGeometry *g1, *g2, *g3;
744  LWGEOM *result;
745  int is3d;
746  int srid;
747 
748  /* A.Difference(Empty) == A */
749  if ( lwgeom_is_empty(geom2) )
750  return lwgeom_clone_deep(geom1);
751 
752  /* Empty.Intersection(A) == Empty */
753  if ( lwgeom_is_empty(geom1) )
754  return lwgeom_clone_deep(geom1);
755 
756  /* ensure srids are identical */
757  srid = (int)(geom1->srid);
758  error_if_srid_mismatch(srid, (int)(geom2->srid));
759 
760  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
761 
762  initGEOS(lwnotice, lwgeom_geos_error);
763 
764  g1 = LWGEOM2GEOS(geom1, 0);
765  if ( 0 == g1 ) /* exception thrown at construction */
766  {
767  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
768  return NULL;
769  }
770 
771  g2 = LWGEOM2GEOS(geom2, 0);
772  if ( 0 == g2 ) /* exception thrown at construction */
773  {
774  GEOSGeom_destroy(g1);
775  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
776  return NULL;
777  }
778 
779  g3 = GEOSDifference(g1,g2);
780 
781  if (g3 == NULL)
782  {
783  GEOSGeom_destroy(g1);
784  GEOSGeom_destroy(g2);
785  lwerror("GEOSDifference: %s", lwgeom_geos_errmsg);
786  return NULL ; /* never get here */
787  }
788 
789  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
790 
791  GEOSSetSRID(g3, srid);
792 
793  result = GEOS2LWGEOM(g3, is3d);
794 
795  if (result == NULL)
796  {
797  GEOSGeom_destroy(g1);
798  GEOSGeom_destroy(g2);
799  GEOSGeom_destroy(g3);
800  lwerror("Error performing difference: GEOS2LWGEOM: %s",
802  return NULL; /* never get here */
803  }
804 
805  GEOSGeom_destroy(g1);
806  GEOSGeom_destroy(g2);
807  GEOSGeom_destroy(g3);
808 
809  /* compressType(result); */
810 
811  return result;
812 }
813 
814 LWGEOM *
815 lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2)
816 {
817  GEOSGeometry *g1, *g2, *g3;
818  LWGEOM *result;
819  int is3d;
820  int srid;
821 
822  /* A.SymDifference(Empty) == A */
823  if ( lwgeom_is_empty(geom2) )
824  return lwgeom_clone_deep(geom1);
825 
826  /* Empty.DymDifference(B) == B */
827  if ( lwgeom_is_empty(geom1) )
828  return lwgeom_clone_deep(geom2);
829 
830  /* ensure srids are identical */
831  srid = (int)(geom1->srid);
832  error_if_srid_mismatch(srid, (int)(geom2->srid));
833 
834  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
835 
836  initGEOS(lwnotice, lwgeom_geos_error);
837 
838  g1 = LWGEOM2GEOS(geom1, 0);
839 
840  if ( 0 == g1 ) /* exception thrown at construction */
841  {
842  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
843  return NULL;
844  }
845 
846  g2 = LWGEOM2GEOS(geom2, 0);
847 
848  if ( 0 == g2 ) /* exception thrown at construction */
849  {
850  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
851  GEOSGeom_destroy(g1);
852  return NULL;
853  }
854 
855  g3 = GEOSSymDifference(g1,g2);
856 
857  if (g3 == NULL)
858  {
859  GEOSGeom_destroy(g1);
860  GEOSGeom_destroy(g2);
861  lwerror("GEOSSymDifference: %s", lwgeom_geos_errmsg);
862  return NULL; /*never get here */
863  }
864 
865  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
866 
867  GEOSSetSRID(g3, srid);
868 
869  result = GEOS2LWGEOM(g3, is3d);
870 
871  if (result == NULL)
872  {
873  GEOSGeom_destroy(g1);
874  GEOSGeom_destroy(g2);
875  GEOSGeom_destroy(g3);
876  lwerror("GEOS symdifference() threw an error (result postgis geometry formation)!");
877  return NULL ; /*never get here */
878  }
879 
880  GEOSGeom_destroy(g1);
881  GEOSGeom_destroy(g2);
882  GEOSGeom_destroy(g3);
883 
884  return result;
885 }
886 
887 LWGEOM*
888 lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
889 {
890  int is3d;
891  int srid;
892  GEOSGeometry *g1, *g2, *g3;
893  LWGEOM *result;
894 
895  LWDEBUG(2, "in geomunion");
896 
897  /* A.Union(empty) == A */
898  if ( lwgeom_is_empty(geom1) )
899  return lwgeom_clone_deep(geom2);
900 
901  /* B.Union(empty) == B */
902  if ( lwgeom_is_empty(geom2) )
903  return lwgeom_clone_deep(geom1);
904 
905 
906  /* ensure srids are identical */
907  srid = (int)(geom1->srid);
908  error_if_srid_mismatch(srid, (int)(geom2->srid));
909 
910  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
911 
912  initGEOS(lwnotice, lwgeom_geos_error);
913 
914  g1 = LWGEOM2GEOS(geom1, 0);
915 
916  if ( 0 == g1 ) /* exception thrown at construction */
917  {
918  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
919  return NULL;
920  }
921 
922  g2 = LWGEOM2GEOS(geom2, 0);
923 
924  if ( 0 == g2 ) /* exception thrown at construction */
925  {
926  GEOSGeom_destroy(g1);
927  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
928  return NULL;
929  }
930 
931  LWDEBUGF(3, "g1=%s", GEOSGeomToWKT(g1));
932  LWDEBUGF(3, "g2=%s", GEOSGeomToWKT(g2));
933 
934  g3 = GEOSUnion(g1,g2);
935 
936  LWDEBUGF(3, "g3=%s", GEOSGeomToWKT(g3));
937 
938  GEOSGeom_destroy(g1);
939  GEOSGeom_destroy(g2);
940 
941  if (g3 == NULL)
942  {
943  lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
944  return NULL; /* never get here */
945  }
946 
947 
948  GEOSSetSRID(g3, srid);
949 
950  result = GEOS2LWGEOM(g3, is3d);
951 
952  GEOSGeom_destroy(g3);
953 
954  if (result == NULL)
955  {
956  lwerror("Error performing union: GEOS2LWGEOM: %s",
958  return NULL; /*never get here */
959  }
960 
961  return result;
962 }
963 
964 LWGEOM *
965 lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1)
966 {
967 #if POSTGIS_GEOS_VERSION < 35
968  lwerror("The GEOS version this postgis binary "
969  "was compiled against (%d) doesn't support "
970  "'GEOSClipByRect' function (3.5.0+ required)",
972  return NULL;
973 #else /* POSTGIS_GEOS_VERSION >= 35 */
974  LWGEOM *result ;
975  GEOSGeometry *g1, *g3 ;
976  int is3d ;
977 
978  /* A.Intersection(Empty) == Empty */
979  if ( lwgeom_is_empty(geom1) )
980  return lwgeom_clone_deep(geom1);
981 
982  is3d = FLAGS_GET_Z(geom1->flags);
983 
984  initGEOS(lwnotice, lwgeom_geos_error);
985 
986  LWDEBUG(3, "clip_by_rect() START");
987 
988  g1 = LWGEOM2GEOS(geom1, 1); /* auto-fix structure */
989  if ( 0 == g1 ) /* exception thrown at construction */
990  {
991  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
992  return NULL ;
993  }
994 
995  LWDEBUG(3, " constructed geometrys - calling geos");
996  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
997  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
998 
999  g3 = GEOSClipByRect(g1,x0,y0,x1,y1);
1000  GEOSGeom_destroy(g1);
1001 
1002  LWDEBUG(3, " clip_by_rect finished");
1003 
1004  if (g3 == NULL)
1005  {
1006  lwnotice("Error performing rectangular clipping: %s", lwgeom_geos_errmsg);
1007  return NULL;
1008  }
1009 
1010  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
1011 
1012  result = GEOS2LWGEOM(g3, is3d);
1013  GEOSGeom_destroy(g3);
1014 
1015  if (result == NULL)
1016  {
1017  lwerror("Error performing intersection: GEOS2LWGEOM: %s", lwgeom_geos_errmsg);
1018  return NULL ; /* never get here */
1019  }
1020 
1021  result->srid = geom1->srid;
1022 
1023  return result ;
1024 #endif /* POSTGIS_GEOS_VERSION >= 35 */
1025 }
1026 
1027 
1028 /* ------------ BuildArea stuff ---------------------------------------------------------------------{ */
1029 
1030 typedef struct Face_t {
1031  const GEOSGeometry* geom;
1032  GEOSGeometry* env;
1033  double envarea;
1034  struct Face_t* parent; /* if this face is an hole of another one, or NULL */
1035 } Face;
1036 
1037 static Face* newFace(const GEOSGeometry* g);
1038 static void delFace(Face* f);
1039 static unsigned int countParens(const Face* f);
1040 static void findFaceHoles(Face** faces, int nfaces);
1041 
1042 static Face*
1043 newFace(const GEOSGeometry* g)
1044 {
1045  Face* f = lwalloc(sizeof(Face));
1046  f->geom = g;
1047  f->env = GEOSEnvelope(f->geom);
1048  GEOSArea(f->env, &f->envarea);
1049  f->parent = NULL;
1050  /* lwnotice("Built Face with area %g and %d holes", f->envarea, GEOSGetNumInteriorRings(f->geom)); */
1051  return f;
1052 }
1053 
1054 static unsigned int
1056 {
1057  unsigned int pcount = 0;
1058  while ( f->parent ) {
1059  ++pcount;
1060  f = f->parent;
1061  }
1062  return pcount;
1063 }
1064 
1065 /* Destroy the face and release memory associated with it */
1066 static void
1068 {
1069  GEOSGeom_destroy(f->env);
1070  lwfree(f);
1071 }
1072 
1073 
1074 static int
1075 compare_by_envarea(const void* g1, const void* g2)
1076 {
1077  Face* f1 = *(Face**)g1;
1078  Face* f2 = *(Face**)g2;
1079  double n1 = f1->envarea;
1080  double n2 = f2->envarea;
1081 
1082  if ( n1 < n2 ) return 1;
1083  if ( n1 > n2 ) return -1;
1084  return 0;
1085 }
1086 
1087 /* Find holes of each face */
1088 static void
1089 findFaceHoles(Face** faces, int nfaces)
1090 {
1091  int i, j, h;
1092 
1093  /* We sort by envelope area so that we know holes are only
1094  * after their shells */
1095  qsort(faces, nfaces, sizeof(Face*), compare_by_envarea);
1096  for (i=0; i<nfaces; ++i) {
1097  Face* f = faces[i];
1098  int nholes = GEOSGetNumInteriorRings(f->geom);
1099  LWDEBUGF(2, "Scanning face %d with env area %g and %d holes", i, f->envarea, nholes);
1100  for (h=0; h<nholes; ++h) {
1101  const GEOSGeometry *hole = GEOSGetInteriorRingN(f->geom, h);
1102  LWDEBUGF(2, "Looking for hole %d/%d of face %d among %d other faces", h+1, nholes, i, nfaces-i-1);
1103  for (j=i+1; j<nfaces; ++j) {
1104  const GEOSGeometry *f2er;
1105  Face* f2 = faces[j];
1106  if ( f2->parent ) continue; /* hole already assigned */
1107  f2er = GEOSGetExteriorRing(f2->geom);
1108  /* TODO: can be optimized as the ring would have the
1109  * same vertices, possibly in different order.
1110  * maybe comparing number of points could already be
1111  * useful.
1112  */
1113  if ( GEOSEquals(f2er, hole) ) {
1114  LWDEBUGF(2, "Hole %d/%d of face %d is face %d", h+1, nholes, i, j);
1115  f2->parent = f;
1116  break;
1117  }
1118  }
1119  }
1120  }
1121 }
1122 
1123 static GEOSGeometry*
1125 {
1126  GEOSGeometry **geoms = lwalloc(sizeof(GEOSGeometry*)*nfaces);
1127  GEOSGeometry *ret;
1128  unsigned int ngeoms = 0;
1129  int i;
1130 
1131  for (i=0; i<nfaces; ++i) {
1132  Face *f = faces[i];
1133  if ( countParens(f) % 2 ) continue; /* we skip odd parents geoms */
1134  geoms[ngeoms++] = GEOSGeom_clone(f->geom);
1135  }
1136 
1137  ret = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, ngeoms);
1138  lwfree(geoms);
1139  return ret;
1140 }
1141 
1142 GEOSGeometry*
1143 LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in)
1144 {
1145  GEOSGeometry *tmp;
1146  GEOSGeometry *geos_result, *shp;
1147  GEOSGeometry const *vgeoms[1];
1148  uint32_t i, ngeoms;
1149  int srid = GEOSGetSRID(geom_in);
1150  Face ** geoms;
1151 
1152  vgeoms[0] = geom_in;
1153 #ifdef LWGEOM_PROFILE_BUILDAREA
1154  lwnotice("Polygonizing");
1155 #endif
1156  geos_result = GEOSPolygonize(vgeoms, 1);
1157 
1158  LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result);
1159 
1160  /* Null return from GEOSpolygonize (an exception) */
1161  if ( ! geos_result ) return 0;
1162 
1163  /*
1164  * We should now have a collection
1165  */
1166 #if PARANOIA_LEVEL > 0
1167  if ( GEOSGeometryTypeId(geos_result) != COLLECTIONTYPE )
1168  {
1169  GEOSGeom_destroy(geos_result);
1170  lwerror("Unexpected return from GEOSpolygonize");
1171  return 0;
1172  }
1173 #endif
1174 
1175  ngeoms = GEOSGetNumGeometries(geos_result);
1176 #ifdef LWGEOM_PROFILE_BUILDAREA
1177  lwnotice("Num geometries from polygonizer: %d", ngeoms);
1178 #endif
1179 
1180 
1181  LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
1182  LWDEBUGF(3, "GEOSpolygonize: polygonized:%s",
1183  lwgeom_to_ewkt(GEOS2LWGEOM(geos_result, 0)));
1184 
1185  /*
1186  * No geometries in collection, early out
1187  */
1188  if ( ngeoms == 0 )
1189  {
1190  GEOSSetSRID(geos_result, srid);
1191  return geos_result;
1192  }
1193 
1194  /*
1195  * Return first geometry if we only have one in collection,
1196  * to avoid the unnecessary Geometry clone below.
1197  */
1198  if ( ngeoms == 1 )
1199  {
1200  tmp = (GEOSGeometry *)GEOSGetGeometryN(geos_result, 0);
1201  if ( ! tmp )
1202  {
1203  GEOSGeom_destroy(geos_result);
1204  return 0; /* exception */
1205  }
1206  shp = GEOSGeom_clone(tmp);
1207  GEOSGeom_destroy(geos_result); /* only safe after the clone above */
1208  GEOSSetSRID(shp, srid);
1209  return shp;
1210  }
1211 
1212  LWDEBUGF(2, "Polygonize returned %d geoms", ngeoms);
1213 
1214  /*
1215  * Polygonizer returns a polygon for each face in the built topology.
1216  *
1217  * This means that for any face with holes we'll have other faces
1218  * representing each hole. We can imagine a parent-child relationship
1219  * between these faces.
1220  *
1221  * In order to maximize the number of visible rings in output we
1222  * only use those faces which have an even number of parents.
1223  *
1224  * Example:
1225  *
1226  * +---------------+
1227  * | L0 | L0 has no parents
1228  * | +---------+ |
1229  * | | L1 | | L1 is an hole of L0
1230  * | | +---+ | |
1231  * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0)
1232  * | | | | | |
1233  * | | +---+ | |
1234  * | +---------+ |
1235  * | |
1236  * +---------------+
1237  *
1238  * See http://trac.osgeo.org/postgis/ticket/1806
1239  *
1240  */
1241 
1242 #ifdef LWGEOM_PROFILE_BUILDAREA
1243  lwnotice("Preparing face structures");
1244 #endif
1245 
1246  /* Prepare face structures for later analysis */
1247  geoms = lwalloc(sizeof(Face**)*ngeoms);
1248  for (i=0; i<ngeoms; ++i)
1249  geoms[i] = newFace(GEOSGetGeometryN(geos_result, i));
1250 
1251 #ifdef LWGEOM_PROFILE_BUILDAREA
1252  lwnotice("Finding face holes");
1253 #endif
1254 
1255  /* Find faces representing other faces holes */
1256  findFaceHoles(geoms, ngeoms);
1257 
1258 #ifdef LWGEOM_PROFILE_BUILDAREA
1259  lwnotice("Colletting even ancestor faces");
1260 #endif
1261 
1262  /* Build a MultiPolygon composed only by faces with an
1263  * even number of ancestors */
1264  tmp = collectFacesWithEvenAncestors(geoms, ngeoms);
1265 
1266 #ifdef LWGEOM_PROFILE_BUILDAREA
1267  lwnotice("Cleaning up");
1268 #endif
1269 
1270  /* Cleanup face structures */
1271  for (i=0; i<ngeoms; ++i) delFace(geoms[i]);
1272  lwfree(geoms);
1273 
1274  /* Faces referenced memory owned by geos_result.
1275  * It is safe to destroy geos_result after deleting them. */
1276  GEOSGeom_destroy(geos_result);
1277 
1278 #ifdef LWGEOM_PROFILE_BUILDAREA
1279  lwnotice("Self-unioning");
1280 #endif
1281 
1282  /* Run a single overlay operation to dissolve shared edges */
1283  shp = GEOSUnionCascaded(tmp);
1284  if ( ! shp )
1285  {
1286  GEOSGeom_destroy(tmp);
1287  return 0; /* exception */
1288  }
1289 
1290 #ifdef LWGEOM_PROFILE_BUILDAREA
1291  lwnotice("Final cleanup");
1292 #endif
1293 
1294  GEOSGeom_destroy(tmp);
1295 
1296  GEOSSetSRID(shp, srid);
1297 
1298  return shp;
1299 }
1300 
1301 LWGEOM*
1303 {
1304  GEOSGeometry* geos_in;
1305  GEOSGeometry* geos_out;
1306  LWGEOM* geom_out;
1307  int SRID = (int)(geom->srid);
1308  int is3d = FLAGS_GET_Z(geom->flags);
1309 
1310  /* Can't build an area from an empty! */
1311  if ( lwgeom_is_empty(geom) )
1312  {
1313  return (LWGEOM*)lwpoly_construct_empty(SRID, is3d, 0);
1314  }
1315 
1316  LWDEBUG(3, "buildarea called");
1317 
1318  LWDEBUGF(3, "ST_BuildArea got geom @ %p", geom);
1319 
1320  initGEOS(lwnotice, lwgeom_geos_error);
1321 
1322  geos_in = LWGEOM2GEOS(geom, 0);
1323 
1324  if ( 0 == geos_in ) /* exception thrown at construction */
1325  {
1326  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1327  return NULL;
1328  }
1329  geos_out = LWGEOM_GEOS_buildArea(geos_in);
1330  GEOSGeom_destroy(geos_in);
1331 
1332  if ( ! geos_out ) /* exception thrown.. */
1333  {
1334  lwerror("LWGEOM_GEOS_buildArea: %s", lwgeom_geos_errmsg);
1335  return NULL;
1336  }
1337 
1338  /* If no geometries are in result collection, return NULL */
1339  if ( GEOSGetNumGeometries(geos_out) == 0 )
1340  {
1341  GEOSGeom_destroy(geos_out);
1342  return NULL;
1343  }
1344 
1345  geom_out = GEOS2LWGEOM(geos_out, is3d);
1346  GEOSGeom_destroy(geos_out);
1347 
1348 #if PARANOIA_LEVEL > 0
1349  if ( geom_out == NULL )
1350  {
1351  lwerror("serialization error");
1352  return NULL;
1353  }
1354 
1355 #endif
1356 
1357  return geom_out;
1358 }
1359 
1360 int
1362 {
1363  GEOSGeometry* geos_in;
1364  int simple;
1365 
1366  /* Empty is always simple */
1367  if ( lwgeom_is_empty(geom) )
1368  {
1369  return 1;
1370  }
1371 
1372  initGEOS(lwnotice, lwgeom_geos_error);
1373 
1374  geos_in = LWGEOM2GEOS(geom, 0);
1375  if ( 0 == geos_in ) /* exception thrown at construction */
1376  {
1377  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1378  return -1;
1379  }
1380  simple = GEOSisSimple(geos_in);
1381  GEOSGeom_destroy(geos_in);
1382 
1383  if ( simple == 2 ) /* exception thrown */
1384  {
1385  lwerror("lwgeom_is_simple: %s", lwgeom_geos_errmsg);
1386  return -1;
1387  }
1388 
1389  return simple ? 1 : 0;
1390 }
1391 
1392 /* ------------ end of BuildArea stuff ---------------------------------------------------------------------} */
1393 
1394 LWGEOM*
1395 lwgeom_geos_noop(const LWGEOM* geom_in)
1396 {
1397  GEOSGeometry *geosgeom;
1398  LWGEOM* geom_out;
1399 
1400  int is3d = FLAGS_GET_Z(geom_in->flags);
1401 
1402  initGEOS(lwnotice, lwgeom_geos_error);
1403  geosgeom = LWGEOM2GEOS(geom_in, 0);
1404  if ( ! geosgeom ) {
1405  lwerror("Geometry could not be converted to GEOS: %s",
1407  return NULL;
1408  }
1409  geom_out = GEOS2LWGEOM(geosgeom, is3d);
1410  GEOSGeom_destroy(geosgeom);
1411  if ( ! geom_out ) {
1412  lwerror("GEOS Geometry could not be converted to LWGEOM: %s",
1414  }
1415  return geom_out;
1416 
1417 }
1418 
1419 LWGEOM*
1420 lwgeom_snap(const LWGEOM* geom1, const LWGEOM* geom2, double tolerance)
1421 {
1422 #if POSTGIS_GEOS_VERSION < 33
1423  lwerror("The GEOS version this lwgeom library "
1424  "was compiled against (%d) doesn't support "
1425  "'Snap' function (3.3.0+ required)",
1427  return NULL;
1428 #else /* POSTGIS_GEOS_VERSION >= 33 */
1429 
1430  int srid, is3d;
1431  GEOSGeometry *g1, *g2, *g3;
1432  LWGEOM* out;
1433 
1434  srid = geom1->srid;
1435  error_if_srid_mismatch(srid, (int)(geom2->srid));
1436 
1437  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
1438 
1439  initGEOS(lwnotice, lwgeom_geos_error);
1440 
1441  g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1, 0);
1442  if ( 0 == g1 ) /* exception thrown at construction */
1443  {
1444  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1445  return NULL;
1446  }
1447 
1448  g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2, 0);
1449  if ( 0 == g2 ) /* exception thrown at construction */
1450  {
1451  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1452  GEOSGeom_destroy(g1);
1453  return NULL;
1454  }
1455 
1456  g3 = GEOSSnap(g1, g2, tolerance);
1457  if (g3 == NULL)
1458  {
1459  GEOSGeom_destroy(g1);
1460  GEOSGeom_destroy(g2);
1461  lwerror("GEOSSnap: %s", lwgeom_geos_errmsg);
1462  return NULL;
1463  }
1464 
1465  GEOSGeom_destroy(g1);
1466  GEOSGeom_destroy(g2);
1467 
1468  GEOSSetSRID(g3, srid);
1469  out = GEOS2LWGEOM(g3, is3d);
1470  if (out == NULL)
1471  {
1472  GEOSGeom_destroy(g3);
1473  lwerror("GEOSSnap() threw an error (result LWGEOM geometry formation)!");
1474  return NULL;
1475  }
1476  GEOSGeom_destroy(g3);
1477 
1478  return out;
1479 
1480 #endif /* POSTGIS_GEOS_VERSION >= 33 */
1481 }
1482 
1483 LWGEOM*
1484 lwgeom_sharedpaths(const LWGEOM* geom1, const LWGEOM* geom2)
1485 {
1486 #if POSTGIS_GEOS_VERSION < 33
1487  lwerror("The GEOS version this postgis binary "
1488  "was compiled against (%d) doesn't support "
1489  "'SharedPaths' function (3.3.0+ required)",
1491  return NULL;
1492 #else /* POSTGIS_GEOS_VERSION >= 33 */
1493  GEOSGeometry *g1, *g2, *g3;
1494  LWGEOM *out;
1495  int is3d, srid;
1496 
1497  srid = geom1->srid;
1498  error_if_srid_mismatch(srid, (int)(geom2->srid));
1499 
1500  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
1501 
1502  initGEOS(lwnotice, lwgeom_geos_error);
1503 
1504  g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1, 0);
1505  if ( 0 == g1 ) /* exception thrown at construction */
1506  {
1507  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1508  return NULL;
1509  }
1510 
1511  g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2, 0);
1512  if ( 0 == g2 ) /* exception thrown at construction */
1513  {
1514  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1515  GEOSGeom_destroy(g1);
1516  return NULL;
1517  }
1518 
1519  g3 = GEOSSharedPaths(g1,g2);
1520 
1521  GEOSGeom_destroy(g1);
1522  GEOSGeom_destroy(g2);
1523 
1524  if (g3 == NULL)
1525  {
1526  lwerror("GEOSSharedPaths: %s", lwgeom_geos_errmsg);
1527  return NULL;
1528  }
1529 
1530  GEOSSetSRID(g3, srid);
1531  out = GEOS2LWGEOM(g3, is3d);
1532  GEOSGeom_destroy(g3);
1533 
1534  if (out == NULL)
1535  {
1536  lwerror("GEOS2LWGEOM threw an error");
1537  return NULL;
1538  }
1539 
1540  return out;
1541 #endif /* POSTGIS_GEOS_VERSION >= 33 */
1542 }
1543 
1544 LWGEOM*
1545 lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
1546 {
1547 #if POSTGIS_GEOS_VERSION < 32
1548  lwerror("lwgeom_offsetcurve: GEOS 3.2 or higher required");
1549 #else
1550  GEOSGeometry *g1, *g3;
1551  LWGEOM *lwgeom_result;
1552  LWGEOM *lwgeom_in = lwline_as_lwgeom(lwline);
1553 
1554  initGEOS(lwnotice, lwgeom_geos_error);
1555 
1556  g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in, 0);
1557  if ( ! g1 )
1558  {
1559  lwerror("lwgeom_offsetcurve: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1560  return NULL;
1561  }
1562 
1563 #if POSTGIS_GEOS_VERSION < 33
1564  /* Size is always positive for GEOSSingleSidedBuffer, and a flag determines left/right */
1565  g3 = GEOSSingleSidedBuffer(g1, size < 0 ? -size : size,
1566  quadsegs, joinStyle, mitreLimit,
1567  size < 0 ? 0 : 1);
1568 #else
1569  g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
1570 #endif
1571  /* Don't need input geometry anymore */
1572  GEOSGeom_destroy(g1);
1573 
1574  if (g3 == NULL)
1575  {
1576  lwerror("GEOSOffsetCurve: %s", lwgeom_geos_errmsg);
1577  return NULL;
1578  }
1579 
1580  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1581 
1582  GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
1583  lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
1584  GEOSGeom_destroy(g3);
1585 
1586  if (lwgeom_result == NULL)
1587  {
1588  lwerror("lwgeom_offsetcurve: GEOS2LWGEOM returned null");
1589  return NULL;
1590  }
1591 
1592  return lwgeom_result;
1593 
1594 #endif /* POSTGIS_GEOS_VERSION < 32 */
1595 }
1596 
1597 LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
1598  int type = GEOSGeomTypeId(geom);
1599  int hasZ;
1600  int SRID = GEOSGetSRID(geom);
1601 
1602  /* GEOS's 0 is equivalent to our unknown as for SRID values */
1603  if ( SRID == 0 ) SRID = SRID_UNKNOWN;
1604 
1605  if ( want3d ) {
1606  hasZ = GEOSHasZ(geom);
1607  if ( ! hasZ ) {
1608  LWDEBUG(3, "Geometry has no Z, won't provide one");
1609  want3d = 0;
1610  }
1611  }
1612 
1613  switch (type) {
1614  LWTRIANGLE **geoms;
1615  uint32_t i, ngeoms;
1616  case GEOS_GEOMETRYCOLLECTION:
1617  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
1618 
1619  ngeoms = GEOSGetNumGeometries(geom);
1620  geoms = NULL;
1621  if ( ngeoms ) {
1622  geoms = lwalloc(ngeoms * sizeof *geoms);
1623  if (!geoms) {
1624  lwerror("lwtin_from_geos: can't allocate geoms");
1625  return NULL;
1626  }
1627  for (i=0; i<ngeoms; i++) {
1628  const GEOSGeometry *poly, *ring;
1629  const GEOSCoordSequence *cs;
1630  POINTARRAY *pa;
1631 
1632  poly = GEOSGetGeometryN(geom, i);
1633  ring = GEOSGetExteriorRing(poly);
1634  cs = GEOSGeom_getCoordSeq(ring);
1635  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
1636 
1637  geoms[i] = lwtriangle_construct(SRID, NULL, pa);
1638  }
1639  }
1640  return (LWTIN *)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM **)geoms);
1641  case GEOS_POLYGON:
1642  case GEOS_MULTIPOINT:
1643  case GEOS_MULTILINESTRING:
1644  case GEOS_MULTIPOLYGON:
1645  case GEOS_LINESTRING:
1646  case GEOS_LINEARRING:
1647  case GEOS_POINT:
1648  lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
1649  break;
1650 
1651  default:
1652  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
1653  return NULL;
1654  }
1655 
1656  /* shouldn't get here */
1657  return NULL;
1658 }
1659 /*
1660  * output = 1 for edges, 2 for TIN, 0 for polygons
1661  */
1662 LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output) {
1663 #if POSTGIS_GEOS_VERSION < 34
1664  lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
1665  return NULL;
1666 #else
1667  GEOSGeometry *g1, *g3;
1668  LWGEOM *lwgeom_result;
1669 
1670  if (output < 0 || output > 2) {
1671  lwerror("lwgeom_delaunay_triangulation: invalid output type specified %d", output);
1672  return NULL;
1673  }
1674 
1675  initGEOS(lwnotice, lwgeom_geos_error);
1676 
1677  g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in, 0);
1678  if ( ! g1 )
1679  {
1680  lwerror("lwgeom_delaunay_triangulation: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1681  return NULL;
1682  }
1683 
1684  /* if output != 1 we want polys */
1685  g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
1686 
1687  /* Don't need input geometry anymore */
1688  GEOSGeom_destroy(g1);
1689 
1690  if (g3 == NULL)
1691  {
1692  lwerror("GEOSDelaunayTriangulation: %s", lwgeom_geos_errmsg);
1693  return NULL;
1694  }
1695 
1696  /* LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); */
1697 
1698  GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
1699 
1700  if (output == 2) {
1701  lwgeom_result = (LWGEOM *)lwtin_from_geos(g3, lwgeom_has_z(lwgeom_in));
1702  } else {
1703  lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
1704  }
1705 
1706  GEOSGeom_destroy(g3);
1707 
1708  if (lwgeom_result == NULL) {
1709  if (output != 2) {
1710  lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
1711  } else {
1712  lwerror("lwgeom_delaunay_triangulation: lwtin_from_geos returned null");
1713  }
1714  return NULL;
1715  }
1716 
1717  return lwgeom_result;
1718 
1719 #endif /* POSTGIS_GEOS_VERSION < 34 */
1720 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:536
double x
Definition: liblwgeom.h:336
#define LINETYPE
Definition: liblwgeom.h:71
LWGEOM * lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)
Definition: lwstroke.c:446
double z
Definition: liblwgeom.h:318
#define POSTGIS_GEOS_VERSION
Definition: sqldefines.h:10
double y
Definition: liblwgeom.h:318
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:51
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
double x
Definition: liblwgeom.h:318
static void delFace(Face *f)
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:61
void lwfree(void *mem)
Definition: lwutil.c:214
int npoints
Definition: liblwgeom.h:355
static int compare_by_envarea(const void *g1, const void *g2)
#define POLYGONTYPE
Definition: liblwgeom.h:72
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:469
uint8_t flags
Definition: liblwgeom.h:381
double xmax
Definition: liblwgeom.h:277
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:330
#define MULTIPOINTTYPE
Definition: liblwgeom.h:73
def fmt
Definition: pixval.py:92
LWGEOM * lwgeom_symdifference(const LWGEOM *geom1, const LWGEOM *geom2)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
GEOSGeometry * GBOX2GEOS(const GBOX *box)
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:822
struct Face_t Face
LWTIN * lwtin_from_geos(const GEOSGeometry *geom, int want3d)
static void findFaceHoles(Face **faces, int nfaces)
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:433
void error_if_srid_mismatch(int srid1, int srid2)
Definition: lwutil.c:341
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:120
POINTARRAY * point
Definition: liblwgeom.h:395
int ptarray_is_closed_2d(const POINTARRAY *pa)
Definition: ptarray.c:694
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:836
int32_t srid
Definition: liblwgeom.h:383
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
double x
Definition: liblwgeom.h:312
GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *)
LWGEOM * lwgeom_linemerge(const LWGEOM *geom1)
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:509
double ymin
Definition: liblwgeom.h:278
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:188
LWGEOM * geom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
double xmin
Definition: liblwgeom.h:276
void lwgeom_geos_error(const char *fmt,...)
LWGEOM * lwgeom_sharedpaths(const LWGEOM *geom1, const LWGEOM *geom2)
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:472
uint8_t flags
Definition: liblwgeom.h:353
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:29
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
LWGEOM * lwgeom_geos_noop(const LWGEOM *geom_in)
Convert an LWGEOM to a GEOS Geometry and convert back – for debug only.
LWGEOM ** geoms
Definition: liblwgeom.h:493
static GEOSGeometry * collectFacesWithEvenAncestors(Face **faces, int nfaces)
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:172
#define TINTYPE
Definition: liblwgeom.h:84
const GEOSGeometry * geom
LWGEOM * lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
static unsigned int countParens(const Face *f)
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1706
LWGEOM * lwgeom_normalize(const LWGEOM *geom1)
LWGEOM * lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
POINTARRAY ** rings
Definition: liblwgeom.h:441
LWGEOM * lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output)
Take vertices of a geometry and build a delaunay triangulation on them.
LWGEOM * lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1)
int nrings
Definition: liblwgeom.h:439
double ymax
Definition: liblwgeom.h:279
int lwgeom_has_arc(const LWGEOM *geom)
Definition: lwstroke.c:41
static Face * newFace(const GEOSGeometry *g)
double y
Definition: liblwgeom.h:312
#define LWGEOM_GEOS_ERRMSG_MAXSIZE
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:124
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
double z
Definition: liblwgeom.h:336
LWTRIANGLE * lwtriangle_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwtriangle.c:27
GEOSGeometry * env
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:75
LWGEOM * lwgeom_unaryunion(const LWGEOM *geom1)
const POINT3DZ * getPoint3dz_cp(const POINTARRAY *pa, int n)
Returns a POINT3DZ pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:486
LWGEOM * lwgeom_snap(const LWGEOM *geom1, const LWGEOM *geom2, double tolerance)
Snap vertices and segments of a geometry to another using a given tolerance.
static GEOSGeometry * ptarray_to_GEOSLinearRing(const POINTARRAY *pa, int autofix)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:70
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:66
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
uint8_t type
Definition: liblwgeom.h:380
type
Definition: ovdump.py:41
struct Face_t * parent
void free(void *)
const char * lwgeom_geos_version()
Return GEOS version string (not to be freed)
void * malloc(YYSIZE_T)
LWGEOM * lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:98
void * lwalloc(size_t size)
Definition: lwutil.c:199
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:1297
double y
Definition: liblwgeom.h:336
#define MULTILINETYPE
Definition: liblwgeom.h:74
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
POINTARRAY * ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, char want3d)
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:136
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:843
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:74
#define COLLECTIONTYPE
Definition: liblwgeom.h:76
This library is the generic geometry handling section of PostGIS.
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
POINTARRAY * points
Definition: liblwgeom.h:406