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