38 const GEOSGeometry *
geom =
nullptr;
39 GEOSGeometry *
env =
nullptr;
44 static Face *newFace(
const GEOSGeometry *g );
45 static void delFace(
Face *
f );
46 static unsigned int countParens(
const Face *
f );
47 static void findFaceHoles(
Face **faces,
int nfaces );
49 static Face *newFace(
const GEOSGeometry *g )
55 f->
env = GEOSEnvelope_r( handle, f->
geom );
61 static unsigned int countParens(
const Face *
f )
63 unsigned int pcount = 0;
73 static void delFace(
Face *f )
76 GEOSGeom_destroy_r( handle, f->
env );
81 static int compare_by_envarea(
const void *g1,
const void *g2 )
88 if ( n1 < n2 )
return 1;
89 if ( n1 > n2 )
return -1;
94 static void findFaceHoles(
Face **faces,
int nfaces )
100 qsort( faces, nfaces,
sizeof(
Face * ), compare_by_envarea );
101 for (
int i = 0; i < nfaces; ++i )
104 int nholes = GEOSGetNumInteriorRings_r( handle, f->
geom );
105 for (
int h = 0; h < nholes; ++h )
107 const GEOSGeometry *hole = GEOSGetInteriorRingN_r( handle, f->
geom, h );
108 for (
int j = i + 1; j < nfaces; ++j )
118 if ( GEOSEquals_r( handle, GEOSGetExteriorRing_r( handle, f2->
geom ), hole ) )
128 static GEOSGeometry *collectFacesWithEvenAncestors(
Face **faces,
int nfaces )
132 unsigned int ngeoms = 0;
133 GEOSGeometry **geoms =
new GEOSGeometry*[nfaces];
134 for (
int i = 0; i < nfaces; ++i )
137 if ( countParens( f ) % 2 )
139 geoms[ngeoms++] = GEOSGeom_clone_r( handle, f->
geom );
142 GEOSGeometry *ret = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOLYGON, geoms, ngeoms );
148 static GEOSGeometry *LWGEOM_GEOS_buildArea(
const GEOSGeometry *geom_in, QString &errorMessage )
152 GEOSGeometry *tmp =
nullptr;
153 GEOSGeometry *shp =
nullptr;
154 GEOSGeometry *geos_result =
nullptr;
155 int srid = GEOSGetSRID_r( handle, geom_in );
157 GEOSGeometry
const *vgeoms[1];
161 geos_result = GEOSPolygonize_r( handle, vgeoms, 1 );
163 catch ( GEOSException &e )
165 errorMessage = QStringLiteral(
"GEOSPolygonize(): %1" ).arg( e.what() );
170 #if PARANOIA_LEVEL > 0 171 if ( GEOSGeometryTypeId_r( handle, geos_result ) != COLLECTIONTYPE )
173 GEOSGeom_destroy_r( handle, geos_result );
174 errorMessage =
"Unexpected return from GEOSpolygonize";
179 int ngeoms = GEOSGetNumGeometries_r( handle, geos_result );
184 GEOSSetSRID_r( handle, geos_result, srid );
192 tmp = ( GEOSGeometry * )GEOSGetGeometryN_r( handle, geos_result, 0 );
195 GEOSGeom_destroy_r( handle, geos_result );
198 shp = GEOSGeom_clone_r( handle, tmp );
199 GEOSGeom_destroy_r( handle, geos_result );
200 GEOSSetSRID_r( handle, shp, srid );
234 for (
int i = 0; i < ngeoms; ++i )
235 geoms[i] = newFace( GEOSGetGeometryN_r( handle, geos_result, i ) );
238 findFaceHoles( geoms, ngeoms );
242 tmp = collectFacesWithEvenAncestors( geoms, ngeoms );
245 for (
int i = 0; i < ngeoms; ++i )
251 GEOSGeom_destroy_r( handle, geos_result );
254 shp = GEOSUnionCascaded_r( handle, tmp );
257 GEOSGeom_destroy_r( handle, tmp );
261 GEOSGeom_destroy_r( handle, tmp );
263 GEOSSetSRID_r( handle, shp, srid );
271 static GEOSGeometry *LWGEOM_GEOS_getPointN(
const GEOSGeometry *g_in, uint32_t n )
276 const GEOSCoordSequence *seq_in =
nullptr;
277 GEOSCoordSeq seq_out;
281 GEOSGeometry *ret =
nullptr;
283 switch ( GEOSGeomTypeId_r( handle, g_in ) )
285 case GEOS_MULTIPOINT:
286 case GEOS_MULTILINESTRING:
287 case GEOS_MULTIPOLYGON:
288 case GEOS_GEOMETRYCOLLECTION:
290 for ( gn = 0; gn < GEOSGetNumGeometries_r( handle, g_in ); ++gn )
292 const GEOSGeometry *g = GEOSGetGeometryN_r( handle, g_in, gn );
293 ret = LWGEOM_GEOS_getPointN( g, n );
294 if ( ret )
return ret;
301 ret = LWGEOM_GEOS_getPointN( GEOSGetExteriorRing_r( handle, g_in ), n );
302 if ( ret )
return ret;
303 for ( gn = 0; gn < GEOSGetNumInteriorRings_r( handle, g_in ); ++gn )
305 const GEOSGeometry *g = GEOSGetInteriorRingN_r( handle, g_in, gn );
306 ret = LWGEOM_GEOS_getPointN( g, n );
307 if ( ret )
return ret;
313 case GEOS_LINESTRING:
314 case GEOS_LINEARRING:
319 seq_in = GEOSGeom_getCoordSeq_r( handle, g_in );
320 if ( ! seq_in )
return nullptr;
321 if ( ! GEOSCoordSeq_getSize_r( handle, seq_in, &sz ) )
return nullptr;
322 if ( ! sz )
return nullptr;
324 if ( ! GEOSCoordSeq_getDimensions_r( handle, seq_in, &dims ) )
return nullptr;
326 seq_out = GEOSCoordSeq_create_r( handle, 1, dims );
327 if ( ! seq_out )
return nullptr;
329 if ( ! GEOSCoordSeq_getX_r( handle, seq_in, n, &val ) )
return nullptr;
330 if ( ! GEOSCoordSeq_setX_r( handle, seq_out, n, val ) )
return nullptr;
331 if ( ! GEOSCoordSeq_getY_r( handle, seq_in, n, &val ) )
return nullptr;
332 if ( ! GEOSCoordSeq_setY_r( handle, seq_out, n, val ) )
return nullptr;
335 if ( ! GEOSCoordSeq_getZ_r( handle, seq_in, n, &val ) )
return nullptr;
336 if ( ! GEOSCoordSeq_setZ_r( handle, seq_out, n, val ) )
return nullptr;
339 return GEOSGeom_createPoint_r( handle, seq_out );
343 static bool lwline_make_geos_friendly(
QgsLineString *line );
344 static bool lwpoly_make_geos_friendly(
QgsPolygon *poly );
351 QgsDebugMsg( QString(
"lwgeom_make_geos_friendly enter (type %1)" ).arg( geom->
wkbType() ) );
362 return lwline_make_geos_friendly( qgsgeometry_cast<QgsLineString *>( geom ) );
367 return lwpoly_make_geos_friendly( qgsgeometry_cast<QgsPolygon *>( geom ) );
373 return lwcollection_make_geos_friendly( qgsgeometry_cast<QgsGeometryCollection *>( geom ) );
377 QgsDebugMsg( QString(
"lwgeom_make_geos_friendly: unsupported input geometry type: %1" ).arg( geom->
wkbType() ) );
384 static bool ring_make_geos_friendly(
QgsCurve *ring )
391 Q_ASSERT_X( linestring,
"ring_make_geos_friendly",
"ring was not a linestring" );
396 if ( p1.
x() != p2.x() || p1.
y() != p2.y() )
408 static bool lwpoly_make_geos_friendly(
QgsPolygon *poly )
417 if ( !ring_make_geos_friendly( qgsgeometry_cast<QgsCurve *>( poly->
interiorRing( i ) ) ) )
438 if ( !lwgeom_make_geos_friendly( g->
geometryN( i ) ) )
447 static GEOSGeometry *LWGEOM_GEOS_nodeLines(
const GEOSGeometry *lines )
456 GEOSGeometry *point = LWGEOM_GEOS_getPointN( lines, 0 );
460 GEOSGeometry *noded =
nullptr;
463 noded = GEOSUnion_r( handle, lines, point );
465 catch ( GEOSException & )
469 GEOSGeom_destroy_r( handle, point );
475 static GEOSGeometry *LWGEOM_GEOS_makeValidPolygon(
const GEOSGeometry *gin, QString &errorMessage )
481 GEOSGeom geos_cut_edges, geos_area, collapse_points;
482 GEOSGeometry *vgeoms[3];
483 unsigned int nvgeoms = 0;
485 Q_ASSERT( GEOSGeomTypeId_r( handle, gin ) == GEOS_POLYGON ||
486 GEOSGeomTypeId_r( handle, gin ) == GEOS_MULTIPOLYGON );
488 geos_bound = GEOSBoundary_r( handle, gin );
494 geos_cut_edges = LWGEOM_GEOS_nodeLines( geos_bound );
495 if ( !geos_cut_edges )
497 GEOSGeom_destroy_r( handle, geos_bound );
498 errorMessage = QStringLiteral(
"LWGEOM_GEOS_nodeLines() failed" );
505 GEOSGeometry *pi =
nullptr;
506 GEOSGeometry *po =
nullptr;
510 pi = GEOSGeom_extractUniquePoints_r( handle, geos_bound );
512 catch ( GEOSException &e )
514 GEOSGeom_destroy_r( handle, geos_bound );
515 errorMessage = QStringLiteral(
"GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
521 po = GEOSGeom_extractUniquePoints_r( handle, geos_cut_edges );
523 catch ( GEOSException &e )
525 GEOSGeom_destroy_r( handle, geos_bound );
526 GEOSGeom_destroy_r( handle, pi );
527 errorMessage = QStringLiteral(
"GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
533 collapse_points = GEOSDifference_r( handle, pi, po );
535 catch ( GEOSException &e )
537 GEOSGeom_destroy_r( handle, geos_bound );
538 GEOSGeom_destroy_r( handle, pi );
539 GEOSGeom_destroy_r( handle, po );
540 errorMessage = QStringLiteral(
"GEOSDifference(): %1" ).arg( e.what() );
544 GEOSGeom_destroy_r( handle, pi );
545 GEOSGeom_destroy_r( handle, po );
547 GEOSGeom_destroy_r( handle, geos_bound );
552 geos_area = GEOSGeom_createEmptyPolygon_r( handle );
554 catch ( GEOSException &e )
556 errorMessage = QStringLiteral(
"GEOSGeom_createEmptyPolygon(): %1" ).arg( e.what() );
557 GEOSGeom_destroy_r( handle, geos_cut_edges );
565 while ( GEOSGetNumGeometries_r( handle, geos_cut_edges ) )
567 GEOSGeometry *new_area =
nullptr;
568 GEOSGeometry *new_area_bound =
nullptr;
569 GEOSGeometry *symdif =
nullptr;
570 GEOSGeometry *new_cut_edges =
nullptr;
576 new_area = LWGEOM_GEOS_buildArea( geos_cut_edges, errorMessage );
578 catch ( GEOSException &e )
580 GEOSGeom_destroy_r( handle, geos_cut_edges );
581 GEOSGeom_destroy_r( handle, geos_area );
582 errorMessage = QStringLiteral(
"LWGEOM_GEOS_buildArea() threw an error: %1" ).arg( e.what() );
586 if ( GEOSisEmpty_r( handle, new_area ) )
589 GEOSGeom_destroy_r( handle, new_area );
599 new_area_bound = GEOSBoundary_r( handle, new_area );
601 catch ( GEOSException &e )
605 errorMessage = QStringLiteral(
"GEOSBoundary() threw an error: %1" ).arg( e.what() );
606 GEOSGeom_destroy_r( handle, new_area );
607 GEOSGeom_destroy_r( handle, geos_area );
614 symdif = GEOSSymDifference_r( handle, geos_area, new_area );
616 catch ( GEOSException &e )
618 GEOSGeom_destroy_r( handle, geos_cut_edges );
619 GEOSGeom_destroy_r( handle, new_area );
620 GEOSGeom_destroy_r( handle, new_area_bound );
621 GEOSGeom_destroy_r( handle, geos_area );
622 errorMessage = QStringLiteral(
"GEOSSymDifference() threw an error: %1" ).arg( e.what() );
626 GEOSGeom_destroy_r( handle, geos_area );
627 GEOSGeom_destroy_r( handle, new_area );
641 new_cut_edges = GEOSDifference_r( handle, geos_cut_edges, new_area_bound );
643 catch ( GEOSException &e )
645 GEOSGeom_destroy_r( handle, geos_cut_edges );
646 GEOSGeom_destroy_r( handle, new_area_bound );
647 GEOSGeom_destroy_r( handle, geos_area );
648 errorMessage = QStringLiteral(
"GEOSDifference() threw an error: %1" ).arg( e.what() );
651 GEOSGeom_destroy_r( handle, geos_cut_edges );
652 GEOSGeom_destroy_r( handle, new_area_bound );
653 geos_cut_edges = new_cut_edges;
656 if ( !GEOSisEmpty_r( handle, geos_area ) )
658 vgeoms[nvgeoms++] = geos_area;
662 GEOSGeom_destroy_r( handle, geos_area );
665 if ( !GEOSisEmpty_r( handle, geos_cut_edges ) )
667 vgeoms[nvgeoms++] = geos_cut_edges;
671 GEOSGeom_destroy_r( handle, geos_cut_edges );
674 if ( !GEOSisEmpty_r( handle, collapse_points ) )
676 vgeoms[nvgeoms++] = collapse_points;
680 GEOSGeom_destroy_r( handle, collapse_points );
693 gout = GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms );
695 catch ( GEOSException &e )
697 errorMessage = QStringLiteral(
"GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
708 static GEOSGeometry *LWGEOM_GEOS_makeValidLine(
const GEOSGeometry *gin, QString &errorMessage )
710 Q_UNUSED( errorMessage );
711 return LWGEOM_GEOS_nodeLines( gin );
714 static GEOSGeometry *LWGEOM_GEOS_makeValidMultiLine(
const GEOSGeometry *gin, QString &errorMessage )
718 int ngeoms = GEOSGetNumGeometries_r( handle, gin );
719 uint32_t nlines_alloc = ngeoms;
720 QVector<GEOSGeometry *> lines;
721 QVector<GEOSGeometry *> points;
722 lines.reserve( nlines_alloc );
723 points.reserve( ngeoms );
725 for (
int i = 0; i < ngeoms; ++i )
727 const GEOSGeometry *g = GEOSGetGeometryN_r( handle, gin, i );
728 GEOSGeometry *vg = LWGEOM_GEOS_makeValidLine( g, errorMessage );
729 if ( GEOSisEmpty_r( handle, vg ) )
732 GEOSGeom_destroy_r( handle, vg );
734 if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_POINT )
738 else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_LINESTRING )
742 else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_MULTILINESTRING )
744 int nsubgeoms = GEOSGetNumGeometries_r( handle, vg );
745 nlines_alloc += nsubgeoms;
746 lines.reserve( nlines_alloc );
747 for (
int j = 0; j < nsubgeoms; ++j )
750 lines.append( GEOSGeom_clone_r( handle, GEOSGetGeometryN_r( handle, vg, j ) ) );
756 errorMessage = QStringLiteral(
"unexpected geom type returned by LWGEOM_GEOS_makeValid: %1" ).arg( GEOSGeomTypeId_r( handle, vg ) );
760 GEOSGeometry *mpoint_out =
nullptr;
761 if ( points.count() )
763 if ( points.count() > 1 )
765 mpoint_out = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOINT, points.data(), points.count() );
769 mpoint_out = points[0];
773 GEOSGeometry *mline_out =
nullptr;
776 if ( lines.count() > 1 )
778 mline_out = GEOSGeom_createCollection_r( handle, GEOS_MULTILINESTRING, lines.data(), lines.count() );
782 mline_out = lines[0];
786 if ( mline_out && mpoint_out )
788 GEOSGeometry *collection[2];
789 collection[0] = mline_out;
790 collection[1] = mpoint_out;
791 return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, collection, 2 );
793 else if ( mline_out )
797 else if ( mpoint_out )
808 static GEOSGeometry *LWGEOM_GEOS_makeValid(
const GEOSGeometry *gin, QString &errorMessage );
811 static GEOSGeometry *LWGEOM_GEOS_makeValidCollection(
const GEOSGeometry *gin, QString &errorMessage )
815 int nvgeoms = GEOSGetNumGeometries_r( handle, gin );
818 errorMessage = QStringLiteral(
"GEOSGetNumGeometries: %1" ).arg( QStringLiteral(
"?" ) );
822 QVector<GEOSGeometry *> vgeoms( nvgeoms );
823 for (
int i = 0; i < nvgeoms; ++i )
825 vgeoms[i] = LWGEOM_GEOS_makeValid( GEOSGetGeometryN_r( handle, gin, i ), errorMessage );
829 GEOSGeom_destroy_r( handle, vgeoms[i] );
838 return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms.data(), nvgeoms );
840 catch ( GEOSException &e )
843 for (
int i = 0; i < nvgeoms; ++i )
844 GEOSGeom_destroy_r( handle, vgeoms[i] );
845 errorMessage = QStringLiteral(
"GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
851 static GEOSGeometry *LWGEOM_GEOS_makeValid(
const GEOSGeometry *gin, QString &errorMessage )
859 if ( GEOSisValid_r( handle, gin ) )
862 return GEOSGeom_clone_r( handle, gin );
865 catch ( GEOSException &e )
868 errorMessage = QStringLiteral(
"GEOSisValid(): %1" ).arg( e.what() );
875 switch ( GEOSGeomTypeId_r( handle, gin ) )
877 case GEOS_MULTIPOINT:
880 QgsDebugMsg(
"PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up" );
883 case GEOS_LINESTRING:
884 return LWGEOM_GEOS_makeValidLine( gin, errorMessage );
886 case GEOS_MULTILINESTRING:
887 return LWGEOM_GEOS_makeValidMultiLine( gin, errorMessage );
890 case GEOS_MULTIPOLYGON:
891 return LWGEOM_GEOS_makeValidPolygon( gin, errorMessage );
893 case GEOS_GEOMETRYCOLLECTION:
894 return LWGEOM_GEOS_makeValidCollection( gin, errorMessage );
897 errorMessage = QStringLiteral(
"ST_MakeValid: doesn't support geometry type: %1" ).arg( GEOSGeomTypeId_r( handle, gin ) );
914 QgsDebugMsg(
"Original geom can't be converted to GEOS - will try cleaning that up first" );
916 std::unique_ptr<QgsAbstractGeometry> lwgeom_in_clone( lwgeom_in->
clone() );
917 if ( !lwgeom_make_geos_friendly( lwgeom_in_clone.get() ) )
919 QgsDebugMsg(
"Could not make a valid geometry out of input" );
928 errorMessage = QStringLiteral(
"Could not convert QGIS geom to GEOS" );
937 GEOSGeometry *geosout = LWGEOM_GEOS_makeValid( geosgeom.get(), errorMessage );
941 std::unique_ptr< QgsAbstractGeometry > lwgeom_out =
QgsGeos::fromGeos( geosout );
942 GEOSGeom_destroy_r( handle, geosout );
966 lwgeom_out.reset( collection );
const QgsCurve * interiorRing(int i) const
static Type multiType(Type type)
Returns the multi type for a WKB type.
Multi point geometry collection.
static bool isMultiType(Type type)
Returns true if the WKB type is a multi type.
Multi line string geometry collection.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
QgsPoint endPoint() const override
Returns the end point of the curve.
static GEOSContextHandle_t getGEOSHandler()
int numPoints() const override
Returns the number of points in the curve.
int numInteriorRings() const
static std::unique_ptr< QgsAbstractGeometry > fromGeos(const GEOSGeometry *geos)
Create a geometry from a GEOSGeometry.
void addVertex(const QgsPoint &pt)
Adds a new vertex to the end of the line string.
#define QgsDebugMsgLevel(str, level)
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
static geos::unique_ptr asGeos(const QgsAbstractGeometry *geom, double precision=0)
Abstract base class for curved geometry type.
Abstract base class for all geometries.
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
QgsPoint startPoint() const override
Returns the starting point of the curve.
int numGeometries() const
Returns the number of geometries within the collection.
#define Q_NOWARN_UNREACHABLE_POP
Multi polygon geometry collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
Line string geometry type, with support for z-dimension and m-values.
std::unique_ptr< QgsAbstractGeometry > _qgis_lwgeom_make_valid(const QgsAbstractGeometry *lwgeom_in, QString &errorMessage)
Implementation of QgsGeometry::makeValid(). Not a public API.
const GEOSGeometry * geom
#define Q_NOWARN_UNREACHABLE_PUSH
int nCoordinates() const override
Returns the number of nodes contained in the geometry.
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
static Type flatType(Type type)
Returns the flat type for a WKB type.
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.