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