QGIS API Documentation  3.10.0-A Coruña (6c816b4204)
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  QgsDebugMsg( QStringLiteral( "lwgeom_make_geos_friendly enter (type %1)" ).arg( geom->wkbType() ) );
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 thes 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( QStringLiteral( "?" ) );
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  QgsDebugMsg( QStringLiteral( "Original geom can't be converted to GEOS - will try cleaning that up first" ) );
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 }
double y
Definition: qgspoint.h:42
static Type multiType(Type type)
Returns the multi type for a WKB type.
Definition: qgswkbtypes.h:301
Multi point geometry collection.
Definition: qgsmultipoint.h:29
static bool isMultiType(Type type)
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:706
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
Multi line string geometry collection.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
QgsPoint endPoint() const override
Returns the end point of the curve.
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:2825
int numPoints() const override
Returns the number of points in the curve.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
static std::unique_ptr< QgsAbstractGeometry > fromGeos(const GEOSGeometry *geos)
Create a geometry from a GEOSGeometry.
Definition: qgsgeos.cpp:1081
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.
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsPoint startPoint() const override
Returns the starting point of the curve.
int numGeometries() const
Returns the number of geometries within the collection.
struct Face_t * parent
#define Q_NOWARN_UNREACHABLE_POP
Definition: qgis.h:652
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:166
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.
const GEOSGeometry * geom
#define Q_NOWARN_UNREACHABLE_PUSH
Definition: qgis.h:651
Polygon geometry type.
Definition: qgspolygon.h:31
int nCoordinates() const override
Returns the number of nodes contained in the geometry.
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:576
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
double x
Definition: qgspoint.h:41