QGIS API Documentation 3.27.0-Master (f261cc1f8b)
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
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
37typedef 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
45static Face *newFace( const GEOSGeometry *g );
46static void delFace( Face *f );
47static unsigned int countParens( const Face *f );
48static void findFaceHoles( Face **faces, int nfaces );
49
50static 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
62static 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
74static void delFace( Face *f )
75{
76 GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
77 GEOSGeom_destroy_r( handle, f->env );
78 delete f;
79}
80
81
82static 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
95static 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
129static 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
149static 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.
272static 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
344static bool lwline_make_geos_friendly( QgsLineString *line );
345static bool lwpoly_make_geos_friendly( QgsPolygon *poly );
346static bool lwcollection_make_geos_friendly( QgsGeometryCollection *g );
347
348
349// Ensure the geometry is "structurally" valid (enough for GEOS to accept it)
350static 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 {
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
385static 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.
409static 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.
426static 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
435static 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
448static 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)
476static 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
709static GEOSGeometry *LWGEOM_GEOS_makeValidLine( const GEOSGeometry *gin, QString &errorMessage )
710{
711 Q_UNUSED( errorMessage )
712 return LWGEOM_GEOS_nodeLines( gin );
713}
714
715static 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
809static GEOSGeometry *LWGEOM_GEOS_makeValid( const GEOSGeometry *gin, QString &errorMessage );
810
811// Will return NULL on error (expect error handler being called by then)
812static 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
852static 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
904std::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.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
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.
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
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:45
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:2952
#define Q_NOWARN_UNREACHABLE_POP
Definition: qgis.h:2953
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
#define QgsDebugMsg(str)
Definition: qgslogger.h:38