QGIS API Documentation  3.17.0-Master (3b262f2a79)
qgsgeometrymakevalid.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgeometrymakevalid.cpp
3  --------------------------------------
4  Date : January 2017
5  Copyright : (C) 2017 by Martin Dobias
6  Email : wonder dot sk at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 // The code in this file has been ported from lwgeom library (part of PostGIS).
17 // See lwgeom_geos_clean.c for the original code by Sandro Santilli.
18 //
19 // Ideally one day the implementation will go to GEOS library...
20 
21 #include "qgsgeometry.h"
22 #include "qgsgeos.h"
23 #include "qgslogger.h"
24 
25 #include "qgsgeometrycollection.h"
26 #include "qgsmultipoint.h"
27 #include "qgsmultilinestring.h"
28 #include "qgsmultipolygon.h"
29 #include "qgspolygon.h"
30 
31 #include <memory>
32 
33 
34 // ------------ BuildArea stuff ---------------------------------------------------------------------
35 
36 typedef struct Face_t
37 {
38  const GEOSGeometry *geom = nullptr;
39  GEOSGeometry *env = nullptr;
40  double envarea;
41  struct Face_t *parent; // if this face is an hole of another one, or NULL
42 } Face;
43 
44 static Face *newFace( const GEOSGeometry *g );
45 static void delFace( Face *f );
46 static unsigned int countParens( const Face *f );
47 static void findFaceHoles( Face **faces, int nfaces );
48 
49 static Face *newFace( const GEOSGeometry *g )
50 {
51  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
52 
53  Face *f = new Face;
54  f->geom = g;
55  f->env = GEOSEnvelope_r( handle, f->geom );
56  GEOSArea_r( handle, f->env, &f->envarea );
57  f->parent = nullptr;
58  return f;
59 }
60 
61 static unsigned int countParens( const Face *f )
62 {
63  unsigned int pcount = 0;
64  while ( f->parent )
65  {
66  ++pcount;
67  f = f->parent;
68  }
69  return pcount;
70 }
71 
72 // Destroy the face and release memory associated with it
73 static void delFace( Face *f )
74 {
75  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
76  GEOSGeom_destroy_r( handle, f->env );
77  delete f;
78 }
79 
80 
81 static int compare_by_envarea( const void *g1, const void *g2 )
82 {
83  Face *f1 = *( Face ** )g1;
84  Face *f2 = *( Face ** )g2;
85  double n1 = f1->envarea;
86  double n2 = f2->envarea;
87 
88  if ( n1 < n2 ) return 1;
89  if ( n1 > n2 ) return -1;
90  return 0;
91 }
92 
93 // Find holes of each face
94 static void findFaceHoles( Face **faces, int nfaces )
95 {
96  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
97 
98  // We sort by envelope area so that we know holes are only
99  // after their shells
100  qsort( faces, nfaces, sizeof( Face * ), compare_by_envarea );
101  for ( int i = 0; i < nfaces; ++i )
102  {
103  Face *f = faces[i];
104  int nholes = GEOSGetNumInteriorRings_r( handle, f->geom );
105  for ( int h = 0; h < nholes; ++h )
106  {
107  const GEOSGeometry *hole = GEOSGetInteriorRingN_r( handle, f->geom, h );
108  for ( int j = i + 1; j < nfaces; ++j )
109  {
110  Face *f2 = faces[j];
111  if ( f2->parent )
112  continue; // hole already assigned
113  /* TODO: can be optimized as the ring would have the
114  * same vertices, possibly in different order.
115  * maybe comparing number of points could already be
116  * useful.
117  */
118  if ( GEOSEquals_r( handle, GEOSGetExteriorRing_r( handle, f2->geom ), hole ) )
119  {
120  f2->parent = f;
121  break;
122  }
123  }
124  }
125  }
126 }
127 
128 static GEOSGeometry *collectFacesWithEvenAncestors( Face **faces, int nfaces )
129 {
130  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
131 
132  unsigned int ngeoms = 0;
133  GEOSGeometry **geoms = new GEOSGeometry*[nfaces];
134  for ( int i = 0; i < nfaces; ++i )
135  {
136  Face *f = faces[i];
137  if ( countParens( f ) % 2 )
138  continue; // we skip odd parents geoms
139  geoms[ngeoms++] = GEOSGeom_clone_r( handle, f->geom );
140  }
141 
142  GEOSGeometry *ret = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOLYGON, geoms, ngeoms );
143  delete [] geoms;
144  return ret;
145 }
146 
148 static GEOSGeometry *LWGEOM_GEOS_buildArea( const GEOSGeometry *geom_in, QString &errorMessage )
149 {
150  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
151 
152  GEOSGeometry *tmp = nullptr;
153  GEOSGeometry *shp = nullptr;
154  GEOSGeometry *geos_result = nullptr;
155  int srid = GEOSGetSRID_r( handle, geom_in );
156 
157  GEOSGeometry const *vgeoms[1];
158  vgeoms[0] = geom_in;
159  try
160  {
161  geos_result = GEOSPolygonize_r( handle, vgeoms, 1 );
162  }
163  catch ( GEOSException &e )
164  {
165  errorMessage = QStringLiteral( "GEOSPolygonize(): %1" ).arg( e.what() );
166  return nullptr;
167  }
168 
169  // We should now have a collection
170 #if PARANOIA_LEVEL > 0
171  if ( GEOSGeometryTypeId_r( handle, geos_result ) != COLLECTIONTYPE )
172  {
173  GEOSGeom_destroy_r( handle, geos_result );
174  errorMessage = "Unexpected return from GEOSpolygonize";
175  return nullptr;
176  }
177 #endif
178 
179  int ngeoms = GEOSGetNumGeometries_r( handle, geos_result );
180 
181  // No geometries in collection, early out
182  if ( ngeoms == 0 )
183  {
184  GEOSSetSRID_r( handle, geos_result, srid );
185  return geos_result;
186  }
187 
188  // Return first geometry if we only have one in collection,
189  // to avoid the unnecessary Geometry clone below.
190  if ( ngeoms == 1 )
191  {
192  tmp = ( GEOSGeometry * )GEOSGetGeometryN_r( handle, geos_result, 0 );
193  if ( ! tmp )
194  {
195  GEOSGeom_destroy_r( handle, geos_result );
196  return nullptr;
197  }
198  shp = GEOSGeom_clone_r( handle, tmp );
199  GEOSGeom_destroy_r( handle, geos_result ); // only safe after the clone above
200  GEOSSetSRID_r( handle, shp, srid );
201  return shp;
202  }
203 
204  /*
205  * Polygonizer returns a polygon for each face in the built topology.
206  *
207  * This means that for any face with holes we'll have other faces
208  * representing each hole. We can imagine a parent-child relationship
209  * between these faces.
210  *
211  * In order to maximize the number of visible rings in output we
212  * only use those faces which have an even number of parents.
213  *
214  * Example:
215  *
216  * +---------------+
217  * | L0 | L0 has no parents
218  * | +---------+ |
219  * | | L1 | | L1 is an hole of L0
220  * | | +---+ | |
221  * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0)
222  * | | | | | |
223  * | | +---+ | |
224  * | +---------+ |
225  * | |
226  * +---------------+
227  *
228  * See http://trac.osgeo.org/postgis/ticket/1806
229  *
230  */
231 
232  // Prepare face structures for later analysis
233  Face **geoms = new Face*[ngeoms];
234  for ( int i = 0; i < ngeoms; ++i )
235  geoms[i] = newFace( GEOSGetGeometryN_r( handle, geos_result, i ) );
236 
237  // Find faces representing other faces holes
238  findFaceHoles( geoms, ngeoms );
239 
240  // Build a MultiPolygon composed only by faces with an
241  // even number of ancestors
242  tmp = collectFacesWithEvenAncestors( geoms, ngeoms );
243 
244  // Cleanup face structures
245  for ( int i = 0; i < ngeoms; ++i )
246  delFace( geoms[i] );
247  delete [] geoms;
248 
249  // Faces referenced memory owned by geos_result.
250  // It is safe to destroy geos_result after deleting them.
251  GEOSGeom_destroy_r( handle, geos_result );
252 
253  // Run a single overlay operation to dissolve shared edges
254  shp = GEOSUnionCascaded_r( handle, tmp );
255  if ( !shp )
256  {
257  GEOSGeom_destroy_r( handle, tmp );
258  return nullptr;
259  }
260 
261  GEOSGeom_destroy_r( handle, tmp );
262 
263  GEOSSetSRID_r( handle, shp, srid );
264 
265  return shp;
266 }
268 
269 // Return Nth vertex in GEOSGeometry as a POINT.
270 // May return NULL if the geometry has NO vertexex.
271 static GEOSGeometry *LWGEOM_GEOS_getPointN( const GEOSGeometry *g_in, uint32_t n )
272 {
273  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
274 
275  uint32_t dims;
276  const GEOSCoordSequence *seq_in = nullptr;
277  GEOSCoordSeq seq_out;
278  double val;
279  uint32_t sz;
280  int gn;
281  GEOSGeometry *ret = nullptr;
282 
283  switch ( GEOSGeomTypeId_r( handle, g_in ) )
284  {
285  case GEOS_MULTIPOINT:
286  case GEOS_MULTILINESTRING:
287  case GEOS_MULTIPOLYGON:
288  case GEOS_GEOMETRYCOLLECTION:
289  {
290  for ( gn = 0; gn < GEOSGetNumGeometries_r( handle, g_in ); ++gn )
291  {
292  const GEOSGeometry *g = GEOSGetGeometryN_r( handle, g_in, gn );
293  ret = LWGEOM_GEOS_getPointN( g, n );
294  if ( ret ) return ret;
295  }
296  break;
297  }
298 
299  case GEOS_POLYGON:
300  {
301  ret = LWGEOM_GEOS_getPointN( GEOSGetExteriorRing_r( handle, g_in ), n );
302  if ( ret ) return ret;
303  for ( gn = 0; gn < GEOSGetNumInteriorRings_r( handle, g_in ); ++gn )
304  {
305  const GEOSGeometry *g = GEOSGetInteriorRingN_r( handle, g_in, gn );
306  ret = LWGEOM_GEOS_getPointN( g, n );
307  if ( ret ) return ret;
308  }
309  break;
310  }
311 
312  case GEOS_POINT:
313  case GEOS_LINESTRING:
314  case GEOS_LINEARRING:
315  break;
316 
317  }
318 
319  seq_in = GEOSGeom_getCoordSeq_r( handle, g_in );
320  if ( ! seq_in ) return nullptr;
321  if ( ! GEOSCoordSeq_getSize_r( handle, seq_in, &sz ) ) return nullptr;
322  if ( ! sz ) return nullptr;
323 
324  if ( ! GEOSCoordSeq_getDimensions_r( handle, seq_in, &dims ) ) return nullptr;
325 
326  seq_out = GEOSCoordSeq_create_r( handle, 1, dims );
327  if ( ! seq_out ) return nullptr;
328 
329  if ( ! GEOSCoordSeq_getX_r( handle, seq_in, n, &val ) ) return nullptr;
330  if ( ! GEOSCoordSeq_setX_r( handle, seq_out, n, val ) ) return nullptr;
331  if ( ! GEOSCoordSeq_getY_r( handle, seq_in, n, &val ) ) return nullptr;
332  if ( ! GEOSCoordSeq_setY_r( handle, seq_out, n, val ) ) return nullptr;
333  if ( dims > 2 )
334  {
335  if ( ! GEOSCoordSeq_getZ_r( handle, seq_in, n, &val ) ) return nullptr;
336  if ( ! GEOSCoordSeq_setZ_r( handle, seq_out, n, val ) ) return nullptr;
337  }
338 
339  return GEOSGeom_createPoint_r( handle, seq_out );
340 }
341 
342 
343 static bool lwline_make_geos_friendly( QgsLineString *line );
344 static bool lwpoly_make_geos_friendly( QgsPolygon *poly );
345 static bool lwcollection_make_geos_friendly( QgsGeometryCollection *g );
346 
347 
348 // Ensure the geometry is "structurally" valid (enough for GEOS to accept it)
349 static bool lwgeom_make_geos_friendly( QgsAbstractGeometry *geom )
350 {
351  QgsDebugMsgLevel( QStringLiteral( "lwgeom_make_geos_friendly enter (type %1)" ).arg( geom->wkbType() ), 3 );
352  switch ( QgsWkbTypes::flatType( geom->wkbType() ) )
353  {
354  case QgsWkbTypes::Point:
356  // a point is always valid
357  return true;
358  break;
359 
361  // lines need at least 2 points
362  return lwline_make_geos_friendly( qgsgeometry_cast<QgsLineString *>( geom ) );
363  break;
364 
366  // polygons need all rings closed and with npoints > 3
367  return lwpoly_make_geos_friendly( qgsgeometry_cast<QgsPolygon *>( geom ) );
368  break;
369 
373  return lwcollection_make_geos_friendly( qgsgeometry_cast<QgsGeometryCollection *>( geom ) );
374  break;
375 
376  default:
377  QgsDebugMsg( QStringLiteral( "lwgeom_make_geos_friendly: unsupported input geometry type: %1" ).arg( geom->wkbType() ) );
378  break;
379  }
380  return false;
381 }
382 
383 
384 static bool ring_make_geos_friendly( QgsCurve *ring )
385 {
386  if ( ring->nCoordinates() == 0 )
387  return false;
388 
389  // earlier we allowed in only geometries with straight segments
390  QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( ring );
391  Q_ASSERT_X( linestring, "ring_make_geos_friendly", "ring was not a linestring" );
392 
393  // close the ring if not already closed (2d only)
394 
395  QgsPoint p1 = linestring->startPoint(), p2 = linestring->endPoint();
396  if ( p1.x() != p2.x() || p1.y() != p2.y() )
397  linestring->addVertex( p1 );
398 
399  // must have at least 4 coordinates to be accepted by GEOS
400 
401  while ( linestring->nCoordinates() < 4 )
402  linestring->addVertex( p1 );
403 
404  return true;
405 }
406 
407 // Make sure all rings are closed and have > 3 points.
408 static bool lwpoly_make_geos_friendly( QgsPolygon *poly )
409 {
410  // If the polygon has no rings there's nothing to do
411  // TODO: in qgis representation there always is exterior ring
412  //if ( ! poly->nrings ) return true;
413 
414  // All rings must be closed and have > 3 points
415  for ( int i = 0; i < poly->numInteriorRings(); i++ )
416  {
417  if ( !ring_make_geos_friendly( qgsgeometry_cast<QgsCurve *>( poly->interiorRing( i ) ) ) )
418  return false;
419  }
420 
421  return true;
422 }
423 
424 // Need NO or >1 points. Duplicate first if only one.
425 static bool lwline_make_geos_friendly( QgsLineString *line )
426 {
427  if ( line->numPoints() == 1 ) // 0 is fine, 2 is fine
428  {
429  line->addVertex( line->startPoint() );
430  }
431  return true;
432 }
433 
434 static bool lwcollection_make_geos_friendly( QgsGeometryCollection *g )
435 {
436  for ( int i = 0; i < g->numGeometries(); i++ )
437  {
438  if ( !lwgeom_make_geos_friendly( g->geometryN( i ) ) )
439  return false;
440  }
441 
442  return true;
443 }
444 
445 
446 // Fully node given linework
447 static GEOSGeometry *LWGEOM_GEOS_nodeLines( const GEOSGeometry *lines )
448 {
449  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
450 
451  // Union with first geometry point, obtaining full noding
452  // and dissolving of duplicated repeated points
453  //
454  // TODO: substitute this with UnaryUnion?
455 
456  GEOSGeometry *point = LWGEOM_GEOS_getPointN( lines, 0 );
457  if ( ! point )
458  return nullptr;
459 
460  GEOSGeometry *noded = nullptr;
461  try
462  {
463  noded = GEOSUnion_r( handle, lines, point );
464  }
465  catch ( GEOSException & )
466  {
467  // no need to do anything here - we'll return nullptr anyway
468  }
469  GEOSGeom_destroy_r( handle, point );
470  return noded;
471 }
472 
473 // Will return NULL on error (expect error handler being called by then)
475 static GEOSGeometry *LWGEOM_GEOS_makeValidPolygon( const GEOSGeometry *gin, QString &errorMessage )
476 {
477  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
478 
479  GEOSGeom gout;
480  GEOSGeom geos_bound;
481  GEOSGeom geos_cut_edges, geos_area, collapse_points;
482  GEOSGeometry *vgeoms[3]; // One for area, one for cut-edges
483  unsigned int nvgeoms = 0;
484 
485  Q_ASSERT( GEOSGeomTypeId_r( handle, gin ) == GEOS_POLYGON ||
486  GEOSGeomTypeId_r( handle, gin ) == GEOS_MULTIPOLYGON );
487 
488  geos_bound = GEOSBoundary_r( handle, gin );
489  if ( !geos_bound )
490  return nullptr;
491 
492  // Use noded boundaries as initial "cut" edges
493 
494  geos_cut_edges = LWGEOM_GEOS_nodeLines( geos_bound );
495  if ( !geos_cut_edges )
496  {
497  GEOSGeom_destroy_r( handle, geos_bound );
498  errorMessage = QStringLiteral( "LWGEOM_GEOS_nodeLines() failed" );
499  return nullptr;
500  }
501 
502  // NOTE: the noding process may drop lines collapsing to points.
503  // We want to retrieve any of those
504  {
505  GEOSGeometry *pi = nullptr;
506  GEOSGeometry *po = nullptr;
507 
508  try
509  {
510  pi = GEOSGeom_extractUniquePoints_r( handle, geos_bound );
511  }
512  catch ( GEOSException &e )
513  {
514  GEOSGeom_destroy_r( handle, geos_bound );
515  errorMessage = QStringLiteral( "GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
516  return nullptr;
517  }
518 
519  try
520  {
521  po = GEOSGeom_extractUniquePoints_r( handle, geos_cut_edges );
522  }
523  catch ( GEOSException &e )
524  {
525  GEOSGeom_destroy_r( handle, geos_bound );
526  GEOSGeom_destroy_r( handle, pi );
527  errorMessage = QStringLiteral( "GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
528  return nullptr;
529  }
530 
531  try
532  {
533  collapse_points = GEOSDifference_r( handle, pi, po );
534  }
535  catch ( GEOSException &e )
536  {
537  GEOSGeom_destroy_r( handle, geos_bound );
538  GEOSGeom_destroy_r( handle, pi );
539  GEOSGeom_destroy_r( handle, po );
540  errorMessage = QStringLiteral( "GEOSDifference(): %1" ).arg( e.what() );
541  return nullptr;
542  }
543 
544  GEOSGeom_destroy_r( handle, pi );
545  GEOSGeom_destroy_r( handle, po );
546  }
547  GEOSGeom_destroy_r( handle, geos_bound );
548 
549  // And use an empty geometry as initial "area"
550  try
551  {
552  geos_area = GEOSGeom_createEmptyPolygon_r( handle );
553  }
554  catch ( GEOSException &e )
555  {
556  errorMessage = QStringLiteral( "GEOSGeom_createEmptyPolygon(): %1" ).arg( e.what() );
557  GEOSGeom_destroy_r( handle, geos_cut_edges );
558  return nullptr;
559  }
560 
561  // See if an area can be build with the remaining edges
562  // and if it can, symdifference with the original area.
563  // Iterate this until no more polygons can be created
564  // with left-over edges.
565  while ( GEOSGetNumGeometries_r( handle, geos_cut_edges ) )
566  {
567  GEOSGeometry *new_area = nullptr;
568  GEOSGeometry *new_area_bound = nullptr;
569  GEOSGeometry *symdif = nullptr;
570  GEOSGeometry *new_cut_edges = nullptr;
571 
572  // ASSUMPTION: cut_edges should already be fully noded
573 
574  try
575  {
576  new_area = LWGEOM_GEOS_buildArea( geos_cut_edges, errorMessage );
577  }
578  catch ( GEOSException &e )
579  {
580  GEOSGeom_destroy_r( handle, geos_cut_edges );
581  GEOSGeom_destroy_r( handle, geos_area );
582  errorMessage = QStringLiteral( "LWGEOM_GEOS_buildArea() threw an error: %1" ).arg( e.what() );
583  return nullptr;
584  }
585 
586  if ( GEOSisEmpty_r( handle, new_area ) )
587  {
588  // no more rings can be build with the edges
589  GEOSGeom_destroy_r( handle, new_area );
590  break;
591  }
592 
593  // We succeeded in building a ring !
594 
595  // Save the new ring boundaries first (to compute
596  // further cut edges later)
597  try
598  {
599  new_area_bound = GEOSBoundary_r( handle, new_area );
600  }
601  catch ( GEOSException &e )
602  {
603  // We did check for empty area already so
604  // this must be some other error
605  errorMessage = QStringLiteral( "GEOSBoundary() threw an error: %1" ).arg( e.what() );
606  GEOSGeom_destroy_r( handle, new_area );
607  GEOSGeom_destroy_r( handle, geos_area );
608  return nullptr;
609  }
610 
611  // Now symdiff new and old area
612  try
613  {
614  symdif = GEOSSymDifference_r( handle, geos_area, new_area );
615  }
616  catch ( GEOSException &e )
617  {
618  GEOSGeom_destroy_r( handle, geos_cut_edges );
619  GEOSGeom_destroy_r( handle, new_area );
620  GEOSGeom_destroy_r( handle, new_area_bound );
621  GEOSGeom_destroy_r( handle, geos_area );
622  errorMessage = QStringLiteral( "GEOSSymDifference() threw an error: %1" ).arg( e.what() );
623  return nullptr;
624  }
625 
626  GEOSGeom_destroy_r( handle, geos_area );
627  GEOSGeom_destroy_r( handle, new_area );
628  geos_area = symdif;
629  symdif = nullptr;
630 
631  // Now let's re-set geos_cut_edges with what's left
632  // from the original boundary.
633  // ASSUMPTION: only the previous cut-edges can be
634  // left, so we don't need to reconsider
635  // the whole original boundaries
636  //
637  // NOTE: this is an expensive operation.
638 
639  try
640  {
641  new_cut_edges = GEOSDifference_r( handle, geos_cut_edges, new_area_bound );
642  }
643  catch ( GEOSException &e )
644  {
645  GEOSGeom_destroy_r( handle, geos_cut_edges );
646  GEOSGeom_destroy_r( handle, new_area_bound );
647  GEOSGeom_destroy_r( handle, geos_area );
648  errorMessage = QStringLiteral( "GEOSDifference() threw an error: %1" ).arg( e.what() );
649  return nullptr;
650  }
651  GEOSGeom_destroy_r( handle, geos_cut_edges );
652  GEOSGeom_destroy_r( handle, new_area_bound );
653  geos_cut_edges = new_cut_edges;
654  }
655 
656  if ( !GEOSisEmpty_r( handle, geos_area ) )
657  {
658  vgeoms[nvgeoms++] = geos_area;
659  }
660  else
661  {
662  GEOSGeom_destroy_r( handle, geos_area );
663  }
664 
665  if ( !GEOSisEmpty_r( handle, geos_cut_edges ) )
666  {
667  vgeoms[nvgeoms++] = geos_cut_edges;
668  }
669  else
670  {
671  GEOSGeom_destroy_r( handle, geos_cut_edges );
672  }
673 
674  if ( !GEOSisEmpty_r( handle, collapse_points ) )
675  {
676  vgeoms[nvgeoms++] = collapse_points;
677  }
678  else
679  {
680  GEOSGeom_destroy_r( handle, collapse_points );
681  }
682 
683  if ( 1 == nvgeoms )
684  {
685  // Return cut edges
686  gout = vgeoms[0];
687  }
688  else
689  {
690  // Collect areas and lines (if any line)
691  try
692  {
693  gout = GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms );
694  }
695  catch ( GEOSException &e )
696  {
697  errorMessage = QStringLiteral( "GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
698  // TODO: cleanup!
699  return nullptr;
700  }
701  }
702 
703  return gout;
704 
705 }
707 
708 static GEOSGeometry *LWGEOM_GEOS_makeValidLine( const GEOSGeometry *gin, QString &errorMessage )
709 {
710  Q_UNUSED( errorMessage )
711  return LWGEOM_GEOS_nodeLines( gin );
712 }
713 
714 static GEOSGeometry *LWGEOM_GEOS_makeValidMultiLine( const GEOSGeometry *gin, QString &errorMessage )
715 {
716  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
717 
718  int ngeoms = GEOSGetNumGeometries_r( handle, gin );
719  uint32_t nlines_alloc = ngeoms;
720  QVector<GEOSGeometry *> lines;
721  QVector<GEOSGeometry *> points;
722  lines.reserve( nlines_alloc );
723  points.reserve( ngeoms );
724 
725  for ( int i = 0; i < ngeoms; ++i )
726  {
727  const GEOSGeometry *g = GEOSGetGeometryN_r( handle, gin, i );
728  GEOSGeometry *vg = LWGEOM_GEOS_makeValidLine( g, errorMessage );
729  if ( GEOSisEmpty_r( handle, vg ) )
730  {
731  // we don't care about this one
732  GEOSGeom_destroy_r( handle, vg );
733  }
734  if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_POINT )
735  {
736  points.append( vg );
737  }
738  else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_LINESTRING )
739  {
740  lines.append( vg );
741  }
742  else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_MULTILINESTRING )
743  {
744  int nsubgeoms = GEOSGetNumGeometries_r( handle, vg );
745  nlines_alloc += nsubgeoms;
746  lines.reserve( nlines_alloc );
747  for ( int j = 0; j < nsubgeoms; ++j )
748  {
749  // NOTE: ownership of the cloned geoms will be taken by final collection
750  lines.append( GEOSGeom_clone_r( handle, GEOSGetGeometryN_r( handle, vg, j ) ) );
751  }
752  }
753  else
754  {
755  // NOTE: return from GEOSGeomType will leak but we really don't expect this to happen
756  errorMessage = QStringLiteral( "unexpected geom type returned by LWGEOM_GEOS_makeValid: %1" ).arg( GEOSGeomTypeId_r( handle, vg ) );
757  }
758  }
759 
760  GEOSGeometry *mpoint_out = nullptr;
761  if ( !points.isEmpty() )
762  {
763  if ( points.count() > 1 )
764  {
765  mpoint_out = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOINT, points.data(), points.count() );
766  }
767  else
768  {
769  mpoint_out = points[0];
770  }
771  }
772 
773  GEOSGeometry *mline_out = nullptr;
774  if ( !lines.isEmpty() )
775  {
776  if ( lines.count() > 1 )
777  {
778  mline_out = GEOSGeom_createCollection_r( handle, GEOS_MULTILINESTRING, lines.data(), lines.count() );
779  }
780  else
781  {
782  mline_out = lines[0];
783  }
784  }
785 
786  if ( mline_out && mpoint_out )
787  {
788  GEOSGeometry *collection[2];
789  collection[0] = mline_out;
790  collection[1] = mpoint_out;
791  return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, collection, 2 );
792  }
793  else if ( mline_out )
794  {
795  return mline_out;
796  }
797  else if ( mpoint_out )
798  {
799  return mpoint_out;
800  }
801  else
802  {
803  return nullptr;
804  }
805 }
806 
807 
808 static GEOSGeometry *LWGEOM_GEOS_makeValid( const GEOSGeometry *gin, QString &errorMessage );
809 
810 // Will return NULL on error (expect error handler being called by then)
811 static GEOSGeometry *LWGEOM_GEOS_makeValidCollection( const GEOSGeometry *gin, QString &errorMessage )
812 {
813  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
814 
815  int nvgeoms = GEOSGetNumGeometries_r( handle, gin );
816  if ( nvgeoms == -1 )
817  {
818  errorMessage = QStringLiteral( "GEOSGetNumGeometries: %1" ).arg( QLatin1String( "?" ) );
819  return nullptr;
820  }
821 
822  QVector<GEOSGeometry *> vgeoms( nvgeoms );
823  for ( int i = 0; i < nvgeoms; ++i )
824  {
825  vgeoms[i] = LWGEOM_GEOS_makeValid( GEOSGetGeometryN_r( handle, gin, i ), errorMessage );
826  if ( ! vgeoms[i] )
827  {
828  while ( i-- )
829  GEOSGeom_destroy_r( handle, vgeoms[i] );
830  // we expect lwerror being called already by makeValid
831  return nullptr;
832  }
833  }
834 
835  // Collect areas and lines (if any line)
836  try
837  {
838  return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms.data(), nvgeoms );
839  }
840  catch ( GEOSException &e )
841  {
842  // cleanup and throw
843  for ( int i = 0; i < nvgeoms; ++i )
844  GEOSGeom_destroy_r( handle, vgeoms[i] );
845  errorMessage = QStringLiteral( "GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
846  return nullptr;
847  }
848 }
849 
850 
851 static GEOSGeometry *LWGEOM_GEOS_makeValid( const GEOSGeometry *gin, QString &errorMessage )
852 {
853  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
854 
855  // return what we got so far if already valid
857  try
858  {
859  if ( GEOSisValid_r( handle, gin ) )
860  {
861  // It's valid at this step, return what we have
862  return GEOSGeom_clone_r( handle, gin );
863  }
864  }
865  catch ( GEOSException &e )
866  {
867  // I don't think should ever happen
868  errorMessage = QStringLiteral( "GEOSisValid(): %1" ).arg( e.what() );
869  return nullptr;
870  }
872 
873  // make what we got valid
874 
875  switch ( GEOSGeomTypeId_r( handle, gin ) )
876  {
877  case GEOS_MULTIPOINT:
878  case GEOS_POINT:
879  // points are always valid, but we might have invalid ordinate values
880  QgsDebugMsg( QStringLiteral( "PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up" ) );
881  return nullptr;
882 
883  case GEOS_LINESTRING:
884  return LWGEOM_GEOS_makeValidLine( gin, errorMessage );
885 
886  case GEOS_MULTILINESTRING:
887  return LWGEOM_GEOS_makeValidMultiLine( gin, errorMessage );
888 
889  case GEOS_POLYGON:
890  case GEOS_MULTIPOLYGON:
891  return LWGEOM_GEOS_makeValidPolygon( gin, errorMessage );
892 
893  case GEOS_GEOMETRYCOLLECTION:
894  return LWGEOM_GEOS_makeValidCollection( gin, errorMessage );
895 
896  default:
897  errorMessage = QStringLiteral( "ST_MakeValid: doesn't support geometry type: %1" ).arg( GEOSGeomTypeId_r( handle, gin ) );
898  return nullptr;
899  }
900 }
901 
902 
903 std::unique_ptr< QgsAbstractGeometry > _qgis_lwgeom_make_valid( const QgsAbstractGeometry *lwgeom_in, QString &errorMessage )
904 {
905  //bool is3d = FLAGS_GET_Z(lwgeom_in->flags);
906 
907  // try to convert to GEOS, if impossible, clean that up first
908  // otherwise (adding only duplicates of existing points)
909  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
910 
911  geos::unique_ptr geosgeom = QgsGeos::asGeos( lwgeom_in );
912  if ( !geosgeom )
913  {
914  QgsDebugMsgLevel( QStringLiteral( "Original geom can't be converted to GEOS - will try cleaning that up first" ), 3 );
915 
916  std::unique_ptr<QgsAbstractGeometry> lwgeom_in_clone( lwgeom_in->clone() );
917  if ( !lwgeom_make_geos_friendly( lwgeom_in_clone.get() ) )
918  {
919  QgsDebugMsg( QStringLiteral( "Could not make a valid geometry out of input" ) );
920  }
921 
922  // try again as we did cleanup now
923  // TODO: invoke LWGEOM2GEOS directly with autoclean ?
924  geosgeom = QgsGeos::asGeos( lwgeom_in_clone.get() );
925 
926  if ( ! geosgeom )
927  {
928  errorMessage = QStringLiteral( "Could not convert QGIS geom to GEOS" );
929  return nullptr;
930  }
931  }
932  else
933  {
934  QgsDebugMsgLevel( QStringLiteral( "original geom converted to GEOS" ), 4 );
935  }
936 
937  GEOSGeometry *geosout = LWGEOM_GEOS_makeValid( geosgeom.get(), errorMessage );
938  if ( !geosout )
939  return nullptr;
940 
941  std::unique_ptr< QgsAbstractGeometry > lwgeom_out = QgsGeos::fromGeos( geosout );
942  GEOSGeom_destroy_r( handle, geosout );
943  if ( !lwgeom_out )
944  return nullptr;
945 
946  // force multi-type if we had a multi-type before
947  if ( QgsWkbTypes::isMultiType( lwgeom_in->wkbType() ) && !QgsWkbTypes::isMultiType( lwgeom_out->wkbType() ) )
948  {
949  QgsGeometryCollection *collection = nullptr;
950  switch ( QgsWkbTypes::multiType( lwgeom_out->wkbType() ) )
951  {
953  collection = new QgsMultiPoint();
954  break;
956  collection = new QgsMultiLineString();
957  break;
959  collection = new QgsMultiPolygon();
960  break;
961  default:
962  collection = new QgsGeometryCollection();
963  break;
964  }
965  collection->addGeometry( lwgeom_out.release() ); // takes ownership
966  lwgeom_out.reset( collection );
967  }
968 
969  return lwgeom_out;
970 }
QgsPoint startPoint() const override SIP_HOLDGIL
Returns the starting point of the curve.
double y
Definition: qgspoint.h:42
int numGeometries() const SIP_HOLDGIL
Returns the number of geometries within the collection.
int numInteriorRings() const SIP_HOLDGIL
Returns the number of interior rings contained with the curve polygon.
static Type multiType(Type type) SIP_HOLDGIL
Returns the multi type for a WKB type.
Definition: qgswkbtypes.h:302
Multi point geometry collection.
Definition: qgsmultipoint.h:29
int numPoints() const override SIP_HOLDGIL
Returns the number of points in the curve.
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
Multi line string geometry collection.
static Type flatType(Type type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:702
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:2928
static std::unique_ptr< QgsAbstractGeometry > fromGeos(const GEOSGeometry *geos)
Create a geometry from a GEOSGeometry.
Definition: qgsgeos.cpp:1098
void addVertex(const QgsPoint &pt)
Adds a new vertex to the end of the line string.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
Geometry collection.
GEOSGeometry * env
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:79
struct Face_t Face
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
Abstract base class for all geometries.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
struct Face_t * parent
const QgsCurve * interiorRing(int i) const SIP_HOLDGIL
Retrieves an interior ring from the curve polygon.
#define Q_NOWARN_UNREACHABLE_POP
Definition: qgis.h:800
static geos::unique_ptr asGeos(const QgsGeometry &geometry, double precision=0)
Returns a geos geometry - caller takes ownership of the object (should be deleted with GEOSGeom_destr...
Definition: qgsgeos.cpp:163
Multi polygon geometry collection.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
std::unique_ptr< QgsAbstractGeometry > _qgis_lwgeom_make_valid(const QgsAbstractGeometry *lwgeom_in, QString &errorMessage)
Implementation of QgsGeometry::makeValid(). Not a public API.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
const GEOSGeometry * geom
int nCoordinates() const override SIP_HOLDGIL
Returns the number of nodes contained in the geometry.
static bool isMultiType(Type type) SIP_HOLDGIL
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:832
QgsPoint endPoint() const override SIP_HOLDGIL
Returns the end point of the curve.
#define Q_NOWARN_UNREACHABLE_PUSH
Definition: qgis.h:799
Polygon geometry type.
Definition: qgspolygon.h:33
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
double x
Definition: qgspoint.h:41