PostGIS  2.1.10dev-r@@SVN_REVISION@@
liblwgeom/lwgeom_geos_clean.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: lwgeom_geos.c 5258 2010-02-17 21:02:49Z strk $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  *
7  * Copyright 2009-2010 Sandro Santilli <strk@keybit.net>
8  *
9  * This is free software; you can redistribute and/or modify it under
10  * the terms of the GNU General Public Licence. See the COPYING file.
11  *
12  **********************************************************************
13  *
14  * ST_MakeValid
15  *
16  * Attempts to make an invalid geometries valid w/out losing
17  * points.
18  *
19  * Polygons may become lines or points or a collection of
20  * polygons lines and points (collapsed ring cases).
21  *
22  * Author: Sandro Santilli <strk@keybit.net>
23  *
24  * Work done for Faunalia (http://www.faunalia.it) with fundings
25  * from Regione Toscana - Sistema Informativo per il Governo
26  * del Territorio e dell'Ambiente (RT-SIGTA).
27  *
28  * Thanks to Dr. Horst Duester for previous work on a plpgsql version
29  * of the cleanup logic [1]
30  *
31  * Thanks to Andrea Peri for recommandations on constraints.
32  *
33  * [1] http://www.sogis1.so.ch/sogis/dl/postgis/cleanGeometry.sql
34  *
35  *
36  **********************************************************************/
37 
38 #include "liblwgeom.h"
39 #include "lwgeom_geos.h"
40 #include "liblwgeom_internal.h"
41 #include "lwgeom_log.h"
42 
43 #include <string.h>
44 #include <stdlib.h>
45 #include <assert.h>
46 
47 /* #define POSTGIS_DEBUG_LEVEL 4 */
48 #undef LWGEOM_PROFILE_MAKEVALID
49 
50 
51 /*
52  * Return Nth vertex in GEOSGeometry as a POINT.
53  * May return NULL if the geometry has NO vertexex.
54  */
55 GEOSGeometry* LWGEOM_GEOS_getPointN(const GEOSGeometry*, uint32_t);
56 GEOSGeometry*
57 LWGEOM_GEOS_getPointN(const GEOSGeometry* g_in, uint32_t n)
58 {
59  uint32_t dims;
60  const GEOSCoordSequence* seq_in;
61  GEOSCoordSeq seq_out;
62  double val;
63  uint32_t sz;
64  int gn;
65  GEOSGeometry* ret;
66 
67  switch ( GEOSGeomTypeId(g_in) )
68  {
69  case GEOS_MULTIPOINT:
70  case GEOS_MULTILINESTRING:
71  case GEOS_MULTIPOLYGON:
72  case GEOS_GEOMETRYCOLLECTION:
73  {
74  for (gn=0; gn<GEOSGetNumGeometries(g_in); ++gn)
75  {
76  const GEOSGeometry* g = GEOSGetGeometryN(g_in, gn);
77  ret = LWGEOM_GEOS_getPointN(g,n);
78  if ( ret ) return ret;
79  }
80  break;
81  }
82 
83  case GEOS_POLYGON:
84  {
85  ret = LWGEOM_GEOS_getPointN(GEOSGetExteriorRing(g_in), n);
86  if ( ret ) return ret;
87  for (gn=0; gn<GEOSGetNumInteriorRings(g_in); ++gn)
88  {
89  const GEOSGeometry* g = GEOSGetInteriorRingN(g_in, gn);
90  ret = LWGEOM_GEOS_getPointN(g, n);
91  if ( ret ) return ret;
92  }
93  break;
94  }
95 
96  case GEOS_POINT:
97  case GEOS_LINESTRING:
98  case GEOS_LINEARRING:
99  break;
100 
101  }
102 
103  seq_in = GEOSGeom_getCoordSeq(g_in);
104  if ( ! seq_in ) return NULL;
105  if ( ! GEOSCoordSeq_getSize(seq_in, &sz) ) return NULL;
106  if ( ! sz ) return NULL;
107 
108  if ( ! GEOSCoordSeq_getDimensions(seq_in, &dims) ) return NULL;
109 
110  seq_out = GEOSCoordSeq_create(1, dims);
111  if ( ! seq_out ) return NULL;
112 
113  if ( ! GEOSCoordSeq_getX(seq_in, n, &val) ) return NULL;
114  if ( ! GEOSCoordSeq_setX(seq_out, n, val) ) return NULL;
115  if ( ! GEOSCoordSeq_getY(seq_in, n, &val) ) return NULL;
116  if ( ! GEOSCoordSeq_setY(seq_out, n, val) ) return NULL;
117  if ( dims > 2 )
118  {
119  if ( ! GEOSCoordSeq_getZ(seq_in, n, &val) ) return NULL;
120  if ( ! GEOSCoordSeq_setZ(seq_out, n, val) ) return NULL;
121  }
122 
123  return GEOSGeom_createPoint(seq_out);
124 }
125 
126 
127 
132 
133 /*
134  * Ensure the geometry is "structurally" valid
135  * (enough for GEOS to accept it)
136  * May return the input untouched (if already valid).
137  * May return geometries of lower dimension (on collapses)
138  */
139 static LWGEOM *
141 {
142  LWDEBUGF(2, "lwgeom_make_geos_friendly enter (type %d)", geom->type);
143  switch (geom->type)
144  {
145  case POINTTYPE:
146  case MULTIPOINTTYPE:
147  /* a point is always valid */
148  return geom;
149  break;
150 
151  case LINETYPE:
152  /* lines need at least 2 points */
153  return lwline_make_geos_friendly((LWLINE *)geom);
154  break;
155 
156  case POLYGONTYPE:
157  /* polygons need all rings closed and with npoints > 3 */
158  return lwpoly_make_geos_friendly((LWPOLY *)geom);
159  break;
160 
161  case MULTILINETYPE:
162  case MULTIPOLYGONTYPE:
163  case COLLECTIONTYPE:
165  break;
166 
167  case CIRCSTRINGTYPE:
168  case COMPOUNDTYPE:
169  case CURVEPOLYTYPE:
170  case MULTISURFACETYPE:
171  case MULTICURVETYPE:
172  default:
173  lwerror("lwgeom_make_geos_friendly: unsupported input geometry type: %s (%d)", lwtype_name(geom->type), geom->type);
174  break;
175  }
176  return 0;
177 }
178 
179 /*
180  * Close the point array, if not already closed in 2d.
181  * Returns the input if already closed in 2d, or a newly
182  * constructed POINTARRAY.
183  * TODO: move in ptarray.c
184  */
186 POINTARRAY*
188 {
189  POINTARRAY* newring;
190 
191  /* close the ring if not already closed (2d only) */
192  if ( ! ptarray_is_closed_2d(ring) )
193  {
194  /* close it up */
195  newring = ptarray_addPoint(ring,
196  getPoint_internal(ring, 0),
197  FLAGS_NDIMS(ring->flags),
198  ring->npoints);
199  ring = newring;
200  }
201  return ring;
202 }
203 
204 /* May return the same input or a new one (never zero) */
205 POINTARRAY*
207 {
208  POINTARRAY* closedring;
209 
210  /* close the ring if not already closed (2d only) */
211  closedring = ptarray_close2d(ring);
212  if (closedring != ring )
213  {
214  ptarray_free(ring); /* should we do this ? */
215  ring = closedring;
216  }
217 
218  /* return 0 for collapsed ring (after closeup) */
219 
220  while ( ring->npoints < 4 )
221  {
222  LWDEBUGF(4, "ring has %d points, adding another", ring->npoints);
223  /* let's add another... */
224  ring = ptarray_addPoint(ring,
225  getPoint_internal(ring, 0),
226  FLAGS_NDIMS(ring->flags),
227  ring->npoints);
228  }
229 
230 
231  return ring;
232 }
233 
234 /* Make sure all rings are closed and have > 3 points.
235  * May return the input untouched.
236  */
237 LWGEOM *
239 {
240  LWGEOM* ret;
241  POINTARRAY **new_rings;
242  int i;
243 
244  /* If the polygon has no rings there's nothing to do */
245  if ( ! poly->nrings ) return (LWGEOM*)poly;
246 
247  /* Allocate enough pointers for all rings */
248  new_rings = lwalloc(sizeof(POINTARRAY*)*poly->nrings);
249 
250  /* All rings must be closed and have > 3 points */
251  for (i=0; i<poly->nrings; i++)
252  {
253  POINTARRAY* ring_in = poly->rings[i];
254  POINTARRAY* ring_out = ring_make_geos_friendly(ring_in);
255 
256  if ( ring_in != ring_out )
257  {
258  LWDEBUGF(3, "lwpoly_make_geos_friendly: ring %d cleaned, now has %d points", i, ring_out->npoints);
259  /* this may come right from
260  * the binary representation lands
261  */
262  /*ptarray_free(ring_in); */
263  }
264  else
265  {
266  LWDEBUGF(3, "lwpoly_make_geos_friendly: ring %d untouched", i);
267  }
268 
269  assert ( ring_out );
270  new_rings[i] = ring_out;
271  }
272 
273  lwfree(poly->rings);
274  poly->rings = new_rings;
275  ret = (LWGEOM*)poly;
276 
277  return ret;
278 }
279 
280 /* Need NO or >1 points. Duplicate first if only one. */
281 LWGEOM *
283 {
284  LWGEOM *ret;
285 
286  if (line->points->npoints == 1) /* 0 is fine, 2 is fine */
287  {
288 #if 1
289  /* Duplicate point */
290  line->points = ptarray_addPoint(line->points,
291  getPoint_internal(line->points, 0),
292  FLAGS_NDIMS(line->points->flags),
293  line->points->npoints);
294  ret = (LWGEOM*)line;
295 #else
296  /* Turn into a point */
297  ret = (LWGEOM*)lwpoint_construct(line->srid, 0, line->points);
298 #endif
299  return ret;
300  }
301  else
302  {
303  return (LWGEOM*)line;
304  /* return lwline_clone(line); */
305  }
306 }
307 
308 LWGEOM *
310 {
311  LWGEOM **new_geoms;
312  uint32_t i, new_ngeoms=0;
313  LWCOLLECTION *ret;
314 
315  /* enough space for all components */
316  new_geoms = lwalloc(sizeof(LWGEOM *)*g->ngeoms);
317 
318  ret = lwalloc(sizeof(LWCOLLECTION));
319  memcpy(ret, g, sizeof(LWCOLLECTION));
320  ret->maxgeoms = g->ngeoms;
321 
322  for (i=0; i<g->ngeoms; i++)
323  {
324  LWGEOM* newg = lwgeom_make_geos_friendly(g->geoms[i]);
325  if ( newg ) new_geoms[new_ngeoms++] = newg;
326  }
327 
328  ret->bbox = NULL; /* recompute later... */
329 
330  ret->ngeoms = new_ngeoms;
331  if ( new_ngeoms )
332  {
333  ret->geoms = new_geoms;
334  }
335  else
336  {
337  free(new_geoms);
338  ret->geoms = NULL;
339  ret->maxgeoms = 0;
340  }
341 
342  return (LWGEOM*)ret;
343 }
344 
345 /*
346  * Fully node given linework
347  */
348 static GEOSGeometry*
349 LWGEOM_GEOS_nodeLines(const GEOSGeometry* lines)
350 {
351  GEOSGeometry* noded;
352  GEOSGeometry* point;
353 
354  /*
355  * Union with first geometry point, obtaining full noding
356  * and dissolving of duplicated repeated points
357  *
358  * TODO: substitute this with UnaryUnion?
359  */
360 
361  point = LWGEOM_GEOS_getPointN(lines, 0);
362  if ( ! point ) return NULL;
363 
364  LWDEBUGF(3,
365  "Boundary point: %s",
366  lwgeom_to_ewkt(GEOS2LWGEOM(point, 0)));
367 
368  noded = GEOSUnion(lines, point);
369  if ( NULL == noded )
370  {
371  GEOSGeom_destroy(point);
372  return NULL;
373  }
374 
375  GEOSGeom_destroy(point);
376 
377  LWDEBUGF(3,
378  "LWGEOM_GEOS_nodeLines: in[%s] out[%s]",
379  lwgeom_to_ewkt(GEOS2LWGEOM(lines, 0)),
380  lwgeom_to_ewkt(GEOS2LWGEOM(noded, 0)));
381 
382  return noded;
383 }
384 
385 #if POSTGIS_GEOS_VERSION >= 33
386 /*
387  * We expect initGEOS being called already.
388  * Will return NULL on error (expect error handler being called by then)
389  *
390  */
391 static GEOSGeometry*
392 LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
393 {
394  GEOSGeom gout;
395  GEOSGeom geos_bound;
396  GEOSGeom geos_cut_edges, geos_area, collapse_points;
397  GEOSGeometry *vgeoms[3]; /* One for area, one for cut-edges */
398  unsigned int nvgeoms=0;
399 
400  assert (GEOSGeomTypeId(gin) == GEOS_POLYGON ||
401  GEOSGeomTypeId(gin) == GEOS_MULTIPOLYGON);
402 
403  geos_bound = GEOSBoundary(gin);
404  if ( NULL == geos_bound )
405  {
406  return NULL;
407  }
408 
409  LWDEBUGF(3,
410  "Boundaries: %s",
411  lwgeom_to_ewkt(GEOS2LWGEOM(geos_bound, 0)));
412 
413  /* Use noded boundaries as initial "cut" edges */
414 
415 #ifdef LWGEOM_PROFILE_MAKEVALID
416  lwnotice("ST_MakeValid: noding lines");
417 #endif
418 
419 
420  geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound);
421  if ( NULL == geos_cut_edges )
422  {
423  GEOSGeom_destroy(geos_bound);
424  lwnotice("LWGEOM_GEOS_nodeLines(): %s", lwgeom_geos_errmsg);
425  return NULL;
426  }
427 
428  /* NOTE: the noding process may drop lines collapsing to points.
429  * We want to retrive any of those */
430  {
431  GEOSGeometry* pi;
432  GEOSGeometry* po;
433 
434 #ifdef LWGEOM_PROFILE_MAKEVALID
435  lwnotice("ST_MakeValid: extracting unique points from bounds");
436 #endif
437 
438  pi = GEOSGeom_extractUniquePoints(geos_bound);
439  if ( NULL == pi )
440  {
441  GEOSGeom_destroy(geos_bound);
442  lwnotice("GEOSGeom_extractUniquePoints(): %s",
444  return NULL;
445  }
446 
447  LWDEBUGF(3,
448  "Boundaries input points %s",
449  lwgeom_to_ewkt(GEOS2LWGEOM(pi, 0)));
450 
451 #ifdef LWGEOM_PROFILE_MAKEVALID
452  lwnotice("ST_MakeValid: extracting unique points from cut_edges");
453 #endif
454 
455  po = GEOSGeom_extractUniquePoints(geos_cut_edges);
456  if ( NULL == po )
457  {
458  GEOSGeom_destroy(geos_bound);
459  GEOSGeom_destroy(pi);
460  lwnotice("GEOSGeom_extractUniquePoints(): %s",
462  return NULL;
463  }
464 
465  LWDEBUGF(3,
466  "Boundaries output points %s",
467  lwgeom_to_ewkt(GEOS2LWGEOM(po, 0)));
468 
469 #ifdef LWGEOM_PROFILE_MAKEVALID
470  lwnotice("ST_MakeValid: find collapse points");
471 #endif
472 
473  collapse_points = GEOSDifference(pi, po);
474  if ( NULL == collapse_points )
475  {
476  GEOSGeom_destroy(geos_bound);
477  GEOSGeom_destroy(pi);
478  GEOSGeom_destroy(po);
479  lwnotice("GEOSDifference(): %s", lwgeom_geos_errmsg);
480  return NULL;
481  }
482 
483  LWDEBUGF(3,
484  "Collapse points: %s",
485  lwgeom_to_ewkt(GEOS2LWGEOM(collapse_points, 0)));
486 
487 #ifdef LWGEOM_PROFILE_MAKEVALID
488  lwnotice("ST_MakeValid: cleanup(1)");
489 #endif
490 
491  GEOSGeom_destroy(pi);
492  GEOSGeom_destroy(po);
493  }
494  GEOSGeom_destroy(geos_bound);
495 
496  LWDEBUGF(3,
497  "Noded Boundaries: %s",
498  lwgeom_to_ewkt(GEOS2LWGEOM(geos_cut_edges, 0)));
499 
500  /* And use an empty geometry as initial "area" */
501  geos_area = GEOSGeom_createEmptyPolygon();
502  if ( ! geos_area )
503  {
504  lwnotice("GEOSGeom_createEmptyPolygon(): %s", lwgeom_geos_errmsg);
505  GEOSGeom_destroy(geos_cut_edges);
506  return NULL;
507  }
508 
509  /*
510  * See if an area can be build with the remaining edges
511  * and if it can, symdifference with the original area.
512  * Iterate this until no more polygons can be created
513  * with left-over edges.
514  */
515  while (GEOSGetNumGeometries(geos_cut_edges))
516  {
517  GEOSGeometry* new_area=0;
518  GEOSGeometry* new_area_bound=0;
519  GEOSGeometry* symdif=0;
520  GEOSGeometry* new_cut_edges=0;
521 
522 #ifdef LWGEOM_PROFILE_MAKEVALID
523  lwnotice("ST_MakeValid: building area from %d edges", GEOSGetNumGeometries(geos_cut_edges));
524 #endif
525 
526  /*
527  * ASSUMPTION: cut_edges should already be fully noded
528  */
529 
530  new_area = LWGEOM_GEOS_buildArea(geos_cut_edges);
531  if ( ! new_area ) /* must be an exception */
532  {
533  GEOSGeom_destroy(geos_cut_edges);
534  GEOSGeom_destroy(geos_area);
535  lwnotice("LWGEOM_GEOS_buildArea() threw an error: %s",
537  return NULL;
538  }
539 
540  if ( GEOSisEmpty(new_area) )
541  {
542  /* no more rings can be build with thes edges */
543  GEOSGeom_destroy(new_area);
544  break;
545  }
546 
547  /*
548  * We succeeded in building a ring !
549  */
550 
551 #ifdef LWGEOM_PROFILE_MAKEVALID
552  lwnotice("ST_MakeValid: ring built with %d cut edges, saving boundaries", GEOSGetNumGeometries(geos_cut_edges));
553 #endif
554 
555  /*
556  * Save the new ring boundaries first (to compute
557  * further cut edges later)
558  */
559  new_area_bound = GEOSBoundary(new_area);
560  if ( ! new_area_bound )
561  {
562  /* We did check for empty area already so
563  * this must be some other error */
564  lwnotice("GEOSBoundary('%s') threw an error: %s",
565  lwgeom_to_ewkt(GEOS2LWGEOM(new_area, 0)),
567  GEOSGeom_destroy(new_area);
568  GEOSGeom_destroy(geos_area);
569  return NULL;
570  }
571 
572 #ifdef LWGEOM_PROFILE_MAKEVALID
573  lwnotice("ST_MakeValid: running SymDifference with new area");
574 #endif
575 
576  /*
577  * Now symdif new and old area
578  */
579  symdif = GEOSSymDifference(geos_area, new_area);
580  if ( ! symdif ) /* must be an exception */
581  {
582  GEOSGeom_destroy(geos_cut_edges);
583  GEOSGeom_destroy(new_area);
584  GEOSGeom_destroy(new_area_bound);
585  GEOSGeom_destroy(geos_area);
586  lwnotice("GEOSSymDifference() threw an error: %s",
588  return NULL;
589  }
590 
591  GEOSGeom_destroy(geos_area);
592  GEOSGeom_destroy(new_area);
593  geos_area = symdif;
594  symdif = 0;
595 
596  /*
597  * Now let's re-set geos_cut_edges with what's left
598  * from the original boundary.
599  * ASSUMPTION: only the previous cut-edges can be
600  * left, so we don't need to reconsider
601  * the whole original boundaries
602  *
603  * NOTE: this is an expensive operation.
604  *
605  */
606 
607 #ifdef LWGEOM_PROFILE_MAKEVALID
608  lwnotice("ST_MakeValid: computing new cut_edges (GEOSDifference)");
609 #endif
610 
611  new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound);
612  GEOSGeom_destroy(new_area_bound);
613  if ( ! new_cut_edges ) /* an exception ? */
614  {
615  /* cleanup and throw */
616  GEOSGeom_destroy(geos_cut_edges);
617  GEOSGeom_destroy(geos_area);
618  /* TODO: Shouldn't this be an lwerror ? */
619  lwnotice("GEOSDifference() threw an error: %s",
621  return NULL;
622  }
623  GEOSGeom_destroy(geos_cut_edges);
624  geos_cut_edges = new_cut_edges;
625  }
626 
627 #ifdef LWGEOM_PROFILE_MAKEVALID
628  lwnotice("ST_MakeValid: final checks");
629 #endif
630 
631  if ( ! GEOSisEmpty(geos_area) )
632  {
633  vgeoms[nvgeoms++] = geos_area;
634  }
635  else
636  {
637  GEOSGeom_destroy(geos_area);
638  }
639 
640  if ( ! GEOSisEmpty(geos_cut_edges) )
641  {
642  vgeoms[nvgeoms++] = geos_cut_edges;
643  }
644  else
645  {
646  GEOSGeom_destroy(geos_cut_edges);
647  }
648 
649  if ( ! GEOSisEmpty(collapse_points) )
650  {
651  vgeoms[nvgeoms++] = collapse_points;
652  }
653  else
654  {
655  GEOSGeom_destroy(collapse_points);
656  }
657 
658  if ( 1 == nvgeoms )
659  {
660  /* Return cut edges */
661  gout = vgeoms[0];
662  }
663  else
664  {
665  /* Collect areas and lines (if any line) */
666  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
667  if ( ! gout ) /* an exception again */
668  {
669  /* cleanup and throw */
670  /* TODO: Shouldn't this be an lwerror ? */
671  lwnotice("GEOSGeom_createCollection() threw an error: %s",
673  /* TODO: cleanup! */
674  return NULL;
675  }
676  }
677 
678  return gout;
679 
680 }
681 
682 static GEOSGeometry*
683 LWGEOM_GEOS_makeValidLine(const GEOSGeometry* gin)
684 {
685  GEOSGeometry* noded;
686  noded = LWGEOM_GEOS_nodeLines(gin);
687  return noded;
688 }
689 
690 static GEOSGeometry*
691 LWGEOM_GEOS_makeValidMultiLine(const GEOSGeometry* gin)
692 {
693  GEOSGeometry** lines;
694  GEOSGeometry** points;
695  GEOSGeometry* mline_out=0;
696  GEOSGeometry* mpoint_out=0;
697  GEOSGeometry* gout=0;
698  uint32_t nlines=0, nlines_alloc;
699  uint32_t npoints=0;
700  uint32_t ngeoms=0, nsubgeoms;
701  uint32_t i, j;
702 
703  ngeoms = GEOSGetNumGeometries(gin);
704 
705  nlines_alloc = ngeoms;
706  lines = lwalloc(sizeof(GEOSGeometry*)*nlines_alloc);
707  points = lwalloc(sizeof(GEOSGeometry*)*ngeoms);
708 
709  for (i=0; i<ngeoms; ++i)
710  {
711  const GEOSGeometry* g = GEOSGetGeometryN(gin, i);
712  GEOSGeometry* vg;
714  if ( GEOSisEmpty(vg) )
715  {
716  /* we don't care about this one */
717  GEOSGeom_destroy(vg);
718  }
719  if ( GEOSGeomTypeId(vg) == GEOS_POINT )
720  {
721  points[npoints++] = vg;
722  }
723  else if ( GEOSGeomTypeId(vg) == GEOS_LINESTRING )
724  {
725  lines[nlines++] = vg;
726  }
727  else if ( GEOSGeomTypeId(vg) == GEOS_MULTILINESTRING )
728  {
729  nsubgeoms=GEOSGetNumGeometries(vg);
730  nlines_alloc += nsubgeoms;
731  lines = lwrealloc(lines, sizeof(GEOSGeometry*)*nlines_alloc);
732  for (j=0; j<nsubgeoms; ++j)
733  {
734  const GEOSGeometry* gc = GEOSGetGeometryN(vg, j);
735  /* NOTE: ownership of the cloned geoms will be
736  * taken by final collection */
737  lines[nlines++] = GEOSGeom_clone(gc);
738  }
739  }
740  else
741  {
742  /* NOTE: return from GEOSGeomType will leak
743  * but we really don't expect this to happen */
744  lwerror("unexpected geom type returned "
745  "by LWGEOM_GEOS_makeValid: %s",
746  GEOSGeomType(vg));
747  }
748  }
749 
750  if ( npoints )
751  {
752  if ( npoints > 1 )
753  {
754  mpoint_out = GEOSGeom_createCollection(GEOS_MULTIPOINT,
755  points, npoints);
756  }
757  else
758  {
759  mpoint_out = points[0];
760  }
761  }
762 
763  if ( nlines )
764  {
765  if ( nlines > 1 )
766  {
767  mline_out = GEOSGeom_createCollection(
768  GEOS_MULTILINESTRING, lines, nlines);
769  }
770  else
771  {
772  mline_out = lines[0];
773  }
774  }
775 
776  lwfree(lines);
777 
778  if ( mline_out && mpoint_out )
779  {
780  points[0] = mline_out;
781  points[1] = mpoint_out;
782  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION,
783  points, 2);
784  }
785  else if ( mline_out )
786  {
787  gout = mline_out;
788  }
789  else if ( mpoint_out )
790  {
791  gout = mpoint_out;
792  }
793 
794  lwfree(points);
795 
796  return gout;
797 }
798 
799 static GEOSGeometry* LWGEOM_GEOS_makeValid(const GEOSGeometry*);
800 
801 /*
802  * We expect initGEOS being called already.
803  * Will return NULL on error (expect error handler being called by then)
804  */
805 static GEOSGeometry*
806 LWGEOM_GEOS_makeValidCollection(const GEOSGeometry* gin)
807 {
808  int nvgeoms;
809  GEOSGeometry **vgeoms;
810  GEOSGeom gout;
811  unsigned int i;
812 
813  nvgeoms = GEOSGetNumGeometries(gin);
814  if ( nvgeoms == -1 ) {
815  lwerror("GEOSGetNumGeometries: %s", lwgeom_geos_errmsg);
816  return 0;
817  }
818 
819  vgeoms = lwalloc( sizeof(GEOSGeometry*) * nvgeoms );
820  if ( ! vgeoms ) {
821  lwerror("LWGEOM_GEOS_makeValidCollection: out of memory");
822  return 0;
823  }
824 
825  for ( i=0; i<nvgeoms; ++i ) {
826  vgeoms[i] = LWGEOM_GEOS_makeValid( GEOSGetGeometryN(gin, i) );
827  if ( ! vgeoms[i] ) {
828  while (i--) GEOSGeom_destroy(vgeoms[i]);
829  lwfree(vgeoms);
830  /* we expect lwerror being called already by makeValid */
831  return NULL;
832  }
833  }
834 
835  /* Collect areas and lines (if any line) */
836  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
837  if ( ! gout ) /* an exception again */
838  {
839  /* cleanup and throw */
840  for ( i=0; i<nvgeoms; ++i ) GEOSGeom_destroy(vgeoms[i]);
841  lwfree(vgeoms);
842  lwerror("GEOSGeom_createCollection() threw an error: %s",
844  return NULL;
845  }
846  lwfree(vgeoms);
847 
848  return gout;
849 
850 }
851 
852 
853 static GEOSGeometry*
854 LWGEOM_GEOS_makeValid(const GEOSGeometry* gin)
855 {
856  GEOSGeometry* gout;
857  char ret_char;
858 
859  /*
860  * Step 2: return what we got so far if already valid
861  */
862 
863  ret_char = GEOSisValid(gin);
864  if ( ret_char == 2 )
865  {
866  /* I don't think should ever happen */
867  lwerror("GEOSisValid(): %s", lwgeom_geos_errmsg);
868  return NULL;
869  }
870  else if ( ret_char )
871  {
872  LWDEBUGF(3,
873  "Geometry [%s] is valid. ",
874  lwgeom_to_ewkt(GEOS2LWGEOM(gin, 0)));
875 
876  /* It's valid at this step, return what we have */
877  return GEOSGeom_clone(gin);
878  }
879 
880  LWDEBUGF(3,
881  "Geometry [%s] is still not valid: %s. "
882  "Will try to clean up further.",
884 
885 
886 
887  /*
888  * Step 3 : make what we got valid
889  */
890 
891  switch (GEOSGeomTypeId(gin))
892  {
893  case GEOS_MULTIPOINT:
894  case GEOS_POINT:
895  /* points are always valid, but we might have invalid ordinate values */
896  lwnotice("PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up");
897  return NULL;
898  break;
899 
900  case GEOS_LINESTRING:
901  gout = LWGEOM_GEOS_makeValidLine(gin);
902  if ( ! gout ) /* an exception or something */
903  {
904  /* cleanup and throw */
906  return NULL;
907  }
908  break; /* we've done */
909 
910  case GEOS_MULTILINESTRING:
911  gout = LWGEOM_GEOS_makeValidMultiLine(gin);
912  if ( ! gout ) /* an exception or something */
913  {
914  /* cleanup and throw */
916  return NULL;
917  }
918  break; /* we've done */
919 
920  case GEOS_POLYGON:
921  case GEOS_MULTIPOLYGON:
922  {
923  gout = LWGEOM_GEOS_makeValidPolygon(gin);
924  if ( ! gout ) /* an exception or something */
925  {
926  /* cleanup and throw */
928  return NULL;
929  }
930  break; /* we've done */
931  }
932 
933  case GEOS_GEOMETRYCOLLECTION:
934  {
936  if ( ! gout ) /* an exception or something */
937  {
938  /* cleanup and throw */
940  return NULL;
941  }
942  break; /* we've done */
943  }
944 
945  default:
946  {
947  char* typname = GEOSGeomType(gin);
948  lwnotice("ST_MakeValid: doesn't support geometry type: %s",
949  typname);
950  GEOSFree(typname);
951  return NULL;
952  break;
953  }
954  }
955 
956  /*
957  * Now check if every point of input is also found
958  * in output, or abort by returning NULL
959  *
960  * Input geometry was lwgeom_in
961  */
962  {
963  const int paranoia = 2;
964  /* TODO: check if the result is valid */
965  if (paranoia)
966  {
967  int loss;
968  GEOSGeometry *pi, *po, *pd;
969 
970  /* TODO: handle some errors here...
971  * Lack of exceptions is annoying indeed,
972  * I'm getting old --strk;
973  */
974  pi = GEOSGeom_extractUniquePoints(gin);
975  po = GEOSGeom_extractUniquePoints(gout);
976  pd = GEOSDifference(pi, po); /* input points - output points */
977  GEOSGeom_destroy(pi);
978  GEOSGeom_destroy(po);
979  loss = !GEOSisEmpty(pd);
980  GEOSGeom_destroy(pd);
981  if ( loss )
982  {
983  lwnotice("Vertices lost in LWGEOM_GEOS_makeValid");
984  /* return NULL */
985  }
986  }
987  }
988 
989 
990  return gout;
991 }
992 
993 /* Exported. Uses GEOS internally */
994 LWGEOM*
996 {
997  int is3d;
998  GEOSGeom geosgeom;
999  GEOSGeometry* geosout;
1000  LWGEOM *lwgeom_out;
1001 
1002  is3d = FLAGS_GET_Z(lwgeom_in->flags);
1003 
1004  /*
1005  * Step 1 : try to convert to GEOS, if impossible, clean that up first
1006  * otherwise (adding only duplicates of existing points)
1007  */
1008 
1010 
1011  lwgeom_out = lwgeom_in;
1012  geosgeom = LWGEOM2GEOS(lwgeom_out);
1013  if ( ! geosgeom )
1014  {
1015  LWDEBUGF(4,
1016  "Original geom can't be converted to GEOS (%s)"
1017  " - will try cleaning that up first",
1019 
1020 
1021  lwgeom_out = lwgeom_make_geos_friendly(lwgeom_out);
1022  if ( ! lwgeom_out )
1023  {
1024  lwerror("Could not make a valid geometry out of input");
1025  }
1026 
1027  /* try again as we did cleanup now */
1028  geosgeom = LWGEOM2GEOS(lwgeom_out);
1029  if ( ! geosgeom )
1030  {
1031  lwerror("Couldn't convert POSTGIS geom to GEOS: %s",
1033  return NULL;
1034  }
1035 
1036  }
1037  else
1038  {
1039  LWDEBUG(4, "original geom converted to GEOS");
1040  lwgeom_out = lwgeom_in;
1041  }
1042 
1043  geosout = LWGEOM_GEOS_makeValid(geosgeom);
1044  GEOSGeom_destroy(geosgeom);
1045  if ( ! geosout )
1046  {
1047  return NULL;
1048  }
1049 
1050  lwgeom_out = GEOS2LWGEOM(geosout, is3d);
1051  GEOSGeom_destroy(geosout);
1052 
1053  if ( lwgeom_is_collection(lwgeom_in) && ! lwgeom_is_collection(lwgeom_out) )
1054  {{
1055  LWGEOM **ogeoms = lwalloc(sizeof(LWGEOM*));
1056  LWGEOM *ogeom;
1057  LWDEBUG(3, "lwgeom_make_valid: forcing multi");
1058  /* NOTE: this is safe because lwgeom_out is surely not lwgeom_in or
1059  * otherwise we couldn't have a collection and a non-collection */
1060  assert(lwgeom_in != lwgeom_out);
1061  ogeoms[0] = lwgeom_out;
1062  ogeom = (LWGEOM *)lwcollection_construct(MULTITYPE[lwgeom_out->type],
1063  lwgeom_out->srid, lwgeom_out->bbox, 1, ogeoms);
1064  lwgeom_out->bbox = NULL;
1065  lwgeom_out = ogeom;
1066  }}
1067 
1068  lwgeom_out->srid = lwgeom_in->srid;
1069  return lwgeom_out;
1070 }
1071 
1072 #endif /* POSTGIS_GEOS_VERSION >= 33 */
1073 
#define LINETYPE
Definition: liblwgeom.h:61
POINTARRAY * ptarray_close2d(POINTARRAY *ring)
GBOX * bbox
Definition: liblwgeom.h:354
LWGEOM * lwgeom_make_valid(LWGEOM *lwgeom_in)
Attempts to make an invalid geometries valid w/out losing points.
#define MULTICURVETYPE
Definition: liblwgeom.h:70
static GEOSGeometry * LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry *gin)
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:947
static LWGEOM * lwgeom_make_geos_friendly(LWGEOM *geom)
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
void lwfree(void *mem)
Definition: lwutil.c:190
int npoints
Definition: liblwgeom.h:327
#define POLYGONTYPE
Definition: liblwgeom.h:62
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:425
uint8_t flags
Definition: liblwgeom.h:353
#define CURVEPOLYTYPE
Definition: liblwgeom.h:69
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:315
#define COMPOUNDTYPE
Definition: liblwgeom.h:68
#define MULTIPOINTTYPE
Definition: liblwgeom.h:63
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
GBOX * bbox
Definition: liblwgeom.h:461
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
static GEOSGeometry * LWGEOM_GEOS_makeValid(const GEOSGeometry *)
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
LWGEOM * lwcollection_make_geos_friendly(LWCOLLECTION *g)
int32_t srid
Definition: liblwgeom.h:377
int ptarray_is_closed_2d(const POINTARRAY *pa)
Definition: ptarray.c:672
int32_t srid
Definition: liblwgeom.h:355
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
POINTARRAY * ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
Add a point in a pointarray.
Definition: ptarray.c:481
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
void lwgeom_geos_error(const char *fmt,...)
uint8_t flags
Definition: liblwgeom.h:325
LWGEOM ** geoms
Definition: liblwgeom.h:465
const GEOSGeometry * geom
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1645
POINTARRAY ** rings
Definition: liblwgeom.h:413
int nrings
Definition: liblwgeom.h:411
static GEOSGeometry * LWGEOM_GEOS_makeValidLine(const GEOSGeometry *gin)
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
uint8_t MULTITYPE[NUMTYPES]
Look-up for the correct MULTI* type promotion for singleton types.
Definition: lwgeom.c:264
#define MULTISURFACETYPE
Definition: liblwgeom.h:71
POINTARRAY * ring_make_geos_friendly(POINTARRAY *ring)
GEOSGeometry * LWGEOM_GEOS_getPointN(const GEOSGeometry *, uint32_t)
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:183
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
uint8_t type
Definition: liblwgeom.h:352
LWGEOM * lwpoly_make_geos_friendly(LWPOLY *poly)
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:67
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:96
void * lwalloc(size_t size)
Definition: lwutil.c:175
LWGEOM * lwline_make_geos_friendly(LWLINE *line)
#define MULTILINETYPE
Definition: liblwgeom.h:64
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:118
static GEOSGeometry * LWGEOM_GEOS_makeValidMultiLine(const GEOSGeometry *gin)
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
This library is the generic geometry handling section of PostGIS.
static GEOSGeometry * LWGEOM_GEOS_nodeLines(const GEOSGeometry *lines)
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
static GEOSGeometry * LWGEOM_GEOS_makeValidCollection(const GEOSGeometry *gin)
POINTARRAY * points
Definition: liblwgeom.h:378