33#define DEFAULT_QUADRANT_SEGMENTS 8
35#define CATCH_GEOS(r) \
36 catch (GEOSException &) \
41#define CATCH_GEOS_WITH_ERRMSG(r) \
42 catch (GEOSException &e) \
46 *errorMsg = e.what(); \
53static void throwGEOSException(
const char *fmt, ... )
59 vsnprintf( buffer,
sizeof buffer, fmt, ap );
62 QString message = QString::fromUtf8( buffer );
72 throw GEOSException( message );
80 throw GEOSException( message );
85static void printGEOSNotice(
const char *fmt, ... )
92 vsnprintf( buffer,
sizeof buffer, fmt, ap );
103#if defined(USE_THREAD_LOCAL) && !defined(Q_OS_WIN)
106QThreadStorage< QgsGeosContext * > QgsGeosContext::sGeosContext;
112#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=5 )
113 mContext = GEOS_init_r();
114 GEOSContext_setNoticeHandler_r( mContext, printGEOSNotice );
115 GEOSContext_setErrorHandler_r( mContext, throwGEOSException );
117 mContext = initGEOS_r( printGEOSNotice, throwGEOSException );
123#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=5 )
124 GEOS_finish_r( mContext );
126 finishGEOS_r( mContext );
132#if defined(USE_THREAD_LOCAL) && !defined(Q_OS_WIN)
133 return sGeosContext.mContext;
135 GEOSContextHandle_t gContext =
nullptr;
136 if ( sGeosContext.hasLocalData() )
138 gContext = sGeosContext.localData()->mContext;
143 gContext = sGeosContext.localData()->mContext;
153void geos::GeosDeleter::operator()(
GEOSGeometry *geom )
const
158void geos::GeosDeleter::operator()(
const GEOSPreparedGeometry *geom )
const
163void geos::GeosDeleter::operator()( GEOSBufferParams *params )
const
168void geos::GeosDeleter::operator()( GEOSCoordSequence *sequence )
const
207#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<10
209 throw QgsNotSupportedException( QObject::tr(
"The structured method to make geometries valid requires a QGIS build based on GEOS 3.10 or later" ) );
212 throw QgsNotSupportedException( QObject::tr(
"The keep collapsed option for making geometries valid requires a QGIS build based on GEOS 3.10 or later" ) );
213 geos::unique_ptr
geos;
216 geos.reset( GEOSMakeValid_r( context, mGeos.get() ) );
221 GEOSMakeValidParams *params = GEOSMakeValidParams_create_r( context );
225 GEOSMakeValidParams_setMethod_r( context, params, GEOS_MAKE_VALID_LINEWORK );
229 GEOSMakeValidParams_setMethod_r( context, params, GEOS_MAKE_VALID_STRUCTURE );
233 GEOSMakeValidParams_setKeepCollapsed_r( context,
235 keepCollapsed ? 1 : 0 );
237 geos::unique_ptr
geos;
240 geos.reset( GEOSMakeValidWithParams_r( context, mGeos.get(), params ) );
241 GEOSMakeValidParams_destroy_r( context, params );
243 catch ( GEOSException &e )
247 *errorMsg = e.what();
249 GEOSMakeValidParams_destroy_r( context, params );
278 std::unique_ptr< QgsAbstractGeometry > geom =
fromGeos( newPart );
285 mGeosPrepared.reset();
319 return overlay( geom, OverlayIntersection, errorMsg, parameters ).release();
324 return overlay( geom, OverlayDifference, errorMsg, parameters ).release();
339 catch ( GEOSException &e )
341 logError( QStringLiteral(
"GEOS" ), e.what() );
344 *errorMsg = e.what();
353 int partType = GEOSGeomTypeId_r( context, currentPart );
356 if ( partType == GEOS_POINT )
367 if ( partType == GEOS_MULTILINESTRING || partType == GEOS_MULTIPOLYGON || partType == GEOS_GEOMETRYCOLLECTION )
369 int partCount = GEOSGetNumGeometries_r( context, currentPart );
370 for (
int i = 0; i < partCount; ++i )
372 subdivideRecursive( GEOSGetGeometryN_r( context, currentPart, i ), maxNodes, depth, parts, clipRect, gridSize );
383 int vertexCount = GEOSGetNumCoordinates_r( context, currentPart );
384 if ( vertexCount == 0 )
388 else if ( vertexCount < maxNodes )
395 double width = clipRect.
width();
396 double height = clipRect.
height();
399 if ( width > height )
412 halfClipRect1.
setYMinimum( halfClipRect1.
yMinimum() - std::numeric_limits<double>::epsilon() );
413 halfClipRect2.
setYMinimum( halfClipRect2.
yMinimum() - std::numeric_limits<double>::epsilon() );
414 halfClipRect1.
setYMaximum( halfClipRect1.
yMaximum() + std::numeric_limits<double>::epsilon() );
415 halfClipRect2.
setYMaximum( halfClipRect2.
yMaximum() + std::numeric_limits<double>::epsilon() );
419 halfClipRect1.
setXMinimum( halfClipRect1.
xMinimum() - std::numeric_limits<double>::epsilon() );
420 halfClipRect2.
setXMinimum( halfClipRect2.
xMinimum() - std::numeric_limits<double>::epsilon() );
421 halfClipRect1.
setXMaximum( halfClipRect1.
xMaximum() + std::numeric_limits<double>::epsilon() );
422 halfClipRect2.
setXMaximum( halfClipRect2.
xMaximum() + std::numeric_limits<double>::epsilon() );
425 geos::unique_ptr clipPart1( GEOSClipByRect_r( context, currentPart, halfClipRect1.
xMinimum(), halfClipRect1.
yMinimum(), halfClipRect1.
xMaximum(), halfClipRect1.
yMaximum() ) );
426 geos::unique_ptr clipPart2( GEOSClipByRect_r( context, currentPart, halfClipRect2.
xMinimum(), halfClipRect2.
yMinimum(), halfClipRect2.
xMaximum(), halfClipRect2.
yMaximum() ) );
434 clipPart1.reset( GEOSIntersectionPrec_r( context, mGeos.get(), clipPart1.get(), gridSize ) );
436 subdivideRecursive( clipPart1.get(), maxNodes, depth, parts, halfClipRect1, gridSize );
442 clipPart2.reset( GEOSIntersectionPrec_r( context, mGeos.get(), clipPart2.get(), gridSize ) );
444 subdivideRecursive( clipPart2.get(), maxNodes, depth, parts, halfClipRect2, gridSize );
456 maxNodes = std::max( maxNodes, 8 );
465 return std::move( parts );
470 return overlay( geom, OverlayUnion, errorMsg, parameters ).release();
475 std::vector<geos::unique_ptr> geosGeometries;
476 geosGeometries.reserve( geomList.size() );
482 geosGeometries.emplace_back(
asGeos( g, mPrecision ) );
486 geos::unique_ptr geomUnion;
489 geos::unique_ptr geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geosGeometries );
492 geomUnion.reset( GEOSUnaryUnionPrec_r( context, geomCollection.get(), parameters.
gridSize() ) );
496 geomUnion.reset( GEOSUnaryUnion_r( context, geomCollection.get() ) );
501 std::unique_ptr< QgsAbstractGeometry > result =
fromGeos( geomUnion.get() );
502 return result.release();
507 std::vector<geos::unique_ptr> geosGeometries;
508 geosGeometries.reserve( geomList.size() );
514 geosGeometries.emplace_back(
asGeos( g.constGet(), mPrecision ) );
518 geos::unique_ptr geomUnion;
521 geos::unique_ptr geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geosGeometries );
525 geomUnion.reset( GEOSUnaryUnionPrec_r( context, geomCollection.get(), parameters.
gridSize() ) );
529 geomUnion.reset( GEOSUnaryUnion_r( context, geomCollection.get() ) );
535 std::unique_ptr< QgsAbstractGeometry > result =
fromGeos( geomUnion.get() );
536 return result.release();
541 return overlay( geom, OverlaySymDifference, errorMsg, parameters ).release();
544static bool isZVerticalLine(
const QgsAbstractGeometry *geom,
double tolerance = 4 * std::numeric_limits<double>::epsilon() )
554 bool isVertical =
true;
555 if (
const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( geom ) )
557 const int nrPoints = line->numPoints();
565 const double sqrTolerance = tolerance * tolerance;
566 const double *lineX = line->xData();
567 const double *lineY = line->yData();
568 for (
int iVert = nrPoints - 1, jVert = 0; jVert < nrPoints; iVert = jVert++ )
589 geos::unique_ptr otherGeosGeom;
598 otherGeosGeom =
asGeos( &firstPoint, mPrecision );
602 otherGeosGeom =
asGeos( geom, mPrecision );
605 if ( !otherGeosGeom )
615 GEOSPreparedDistance_r( context, mGeosPrepared.get(), otherGeosGeom.get(), &
distance );
619 GEOSDistance_r( context, mGeos.get(), otherGeosGeom.get(), &
distance );
635 geos::unique_ptr point = createGeosPointXY( x, y,
false, 0,
false, 0, 2, 0 );
644 GEOSPreparedDistance_r( context, mGeosPrepared.get(), point.get(), &
distance );
648 GEOSDistance_r( context, mGeos.get(), point.get(), &
distance );
663 geos::unique_ptr otherGeosGeom;
672 otherGeosGeom =
asGeos( &firstPoint );
676 otherGeosGeom =
asGeos( geom, mPrecision );
679 if ( !otherGeosGeom )
694#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
695 return GEOSPreparedDistanceWithin_r( context, mGeosPrepared.get(), otherGeosGeom.get(), maxdist );
697 GEOSPreparedDistance_r( context, mGeosPrepared.get(), otherGeosGeom.get(), &
distance );
702#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
703 return GEOSDistanceWithin_r( context, mGeos.get(), otherGeosGeom.get(), maxdist );
705 GEOSDistance_r( context, mGeos.get(), otherGeosGeom.get(), &
distance );
720#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
723 geos::unique_ptr point = createGeosPointXY( x, y,
false, 0,
false, 0, 2, 0 );
729#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
730 return GEOSPreparedContainsXY_r( context, mGeosPrepared.get(), x, y ) == 1;
732 return GEOSPreparedContains_r( context, mGeosPrepared.get(), point.get() ) == 1;
736#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
737 geos::unique_ptr point = createGeosPointXY( x, y,
false, 0,
false, 0, 2, 0 );
742 result = ( GEOSContains_r( context, mGeos.get(), point.get() ) == 1 );
744 catch ( GEOSException &e )
746 logError( QStringLiteral(
"GEOS" ), e.what() );
749 *errorMsg = e.what();
765 geos::unique_ptr otherGeosGeom(
asGeos( geom, mPrecision ) );
766 if ( !otherGeosGeom )
788 geos::unique_ptr otherGeosGeom(
asGeos( geom, mPrecision ) );
789 if ( !otherGeosGeom )
811 geos::unique_ptr otherGeosGeom(
asGeos( geom, mPrecision ) );
812 if ( !otherGeosGeom )
834 geos::unique_ptr otherGeosGeom(
asGeos( geom, mPrecision ) );
835 if ( !otherGeosGeom )
851 if ( !mGeos || !geom )
856#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
864 return GEOSPreparedIntersectsXY_r(
QgsGeosContext::get(), mGeosPrepared.get(), point->x(), point->y() ) == 1;
866 catch ( GEOSException &e )
868 logError( QStringLiteral(
"GEOS" ), e.what() );
871 *errorMsg = e.what();
879 return relation( geom, RelationIntersects, errorMsg );
884 return relation( geom, RelationTouches, errorMsg );
889 return relation( geom, RelationCrosses, errorMsg );
894 return relation( geom, RelationWithin, errorMsg );
899 return relation( geom, RelationOverlaps, errorMsg );
904 if ( !mGeos || !geom )
909#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
917 return GEOSPreparedContainsXY_r(
QgsGeosContext::get(), mGeosPrepared.get(), point->x(), point->y() ) == 1;
919 catch ( GEOSException &e )
921 logError( QStringLiteral(
"GEOS" ), e.what() );
924 *errorMsg = e.what();
932 return relation( geom, RelationContains, errorMsg );
937 return relation( geom, RelationDisjoint, errorMsg );
947 geos::unique_ptr geosGeom(
asGeos( geom, mPrecision ) );
957 char *r = GEOSRelate_r( context, mGeos.get(), geosGeom.get() );
960 result = QString( r );
961 GEOSFree_r( context, r );
964 catch ( GEOSException &e )
966 logError( QStringLiteral(
"GEOS" ), e.what() );
969 *errorMsg = e.what();
978 if ( !mGeos || !geom )
983 geos::unique_ptr geosGeom(
asGeos( geom, mPrecision ) );
993 result = ( GEOSRelatePattern_r( context, mGeos.get(), geosGeom.get(), pattern.toLocal8Bit().constData() ) == 1 );
995 catch ( GEOSException &e )
997 logError( QStringLiteral(
"GEOS" ), e.what() );
1000 *errorMsg = e.what();
1041 QVector<QgsGeometry> &newGeometries,
1044 QString *errorMsg,
bool skipIntersectionCheck )
const
1060 if ( !GEOSisValid_r( context, mGeos.get() ) )
1068 newGeometries.clear();
1069 geos::unique_ptr splitLineGeos;
1075 splitLineGeos = createGeosLinestring( &splitLine, mPrecision );
1079 splitLineGeos = createGeosPointXY( splitLine.
xAt( 0 ), splitLine.
yAt( 0 ),
false, 0,
false, 0, 2, mPrecision );
1086 if ( !GEOSisValid_r( context, splitLineGeos.get() ) || !GEOSisSimple_r( context, splitLineGeos.get() ) )
1094 if ( !topologicalTestPointsSplit( splitLineGeos.get(), topologyTestPoints ) )
1103 returnCode = splitLinearGeometry( splitLineGeos.get(), newGeometries, skipIntersectionCheck );
1107 returnCode = splitPolygonGeometry( splitLineGeos.get(), newGeometries, skipIntersectionCheck );
1136 geos::unique_ptr intersectionGeom( GEOSIntersection_r( context, mGeos.get(), splitLine ) );
1137 if ( !intersectionGeom )
1140 bool simple =
false;
1141 int nIntersectGeoms = 1;
1142 if ( GEOSGeomTypeId_r( context, intersectionGeom.get() ) == GEOS_LINESTRING
1143 || GEOSGeomTypeId_r( context, intersectionGeom.get() ) == GEOS_POINT )
1147 nIntersectGeoms = GEOSGetNumGeometries_r( context, intersectionGeom.get() );
1149 for (
int i = 0; i < nIntersectGeoms; ++i )
1153 currentIntersectGeom = intersectionGeom.get();
1155 currentIntersectGeom = GEOSGetGeometryN_r( context, intersectionGeom.get(), i );
1157 const GEOSCoordSequence *lineSequence = GEOSGeom_getCoordSeq_r( context, currentIntersectGeom );
1158 unsigned int sequenceSize = 0;
1160 if ( GEOSCoordSeq_getSize_r( context, lineSequence, &sequenceSize ) != 0 )
1162 for (
unsigned int i = 0; i < sequenceSize; ++i )
1164 if ( GEOSCoordSeq_getXYZ_r( context, lineSequence, i, &x, &y, &z ) )
1166 testPoints.push_back(
QgsPoint( x, y, z ) );
1180 int type = GEOSGeomTypeId_r( context, mGeos.get() );
1182 std::unique_ptr< QgsMultiCurve > multiCurve;
1183 if ( type == GEOS_MULTILINESTRING )
1185 multiCurve.reset( qgsgeometry_cast<QgsMultiCurve *>(
mGeometry->
clone() ) );
1187 else if ( type == GEOS_LINESTRING )
1205 std::unique_ptr< QgsMultiPoint > splitPoints;
1207 std::unique_ptr< QgsAbstractGeometry > splitGeom(
fromGeos( GEOSsplitPoint ) );
1209 if ( qgsgeometry_cast<QgsMultiPoint *>( splitGeom.get() ) )
1211 splitPoints.reset( qgsgeometry_cast<QgsMultiPoint *>( splitGeom.release() ) );
1213 else if ( qgsgeometry_cast<QgsPoint *>( splitGeom.get() ) )
1215 splitPoints = std::make_unique< QgsMultiPoint >();
1216 if ( qgsgeometry_cast<QgsPoint *>( splitGeom.get() ) )
1218 splitPoints->addGeometry( qgsgeometry_cast<QgsPoint *>( splitGeom.release() ) );
1226 for (
int geometryIndex = 0; geometryIndex < multiCurve->numGeometries(); ++geometryIndex )
1228 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( multiCurve->geometryN( geometryIndex ) );
1231 const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( multiCurve->geometryN( geometryIndex ) );
1239 QMap< int, QVector< QPair< double, QgsPoint > > >pointMap;
1240 for (
int splitPointIndex = 0; splitPointIndex < splitPoints->numGeometries(); ++splitPointIndex )
1242 const QgsPoint *intersectionPoint = splitPoints->pointN( splitPointIndex );
1248 line->
closestSegment( *intersectionPoint, segmentPoint2D, nextVertex );
1264 const QPair< double, QgsPoint > pair = qMakePair(
distance, *correctSegmentPoint.get() );
1265 if ( pointMap.contains( nextVertex.
vertex - 1 ) )
1266 pointMap[ nextVertex.
vertex - 1 ].append( pair );
1268 pointMap[ nextVertex.
vertex - 1 ] = QVector< QPair< double, QgsPoint > >() << pair;
1273 for (
auto &p : pointMap )
1275 std::sort( p.begin(), p.end(), [](
const QPair< double, QgsPoint > &a,
const QPair< double, QgsPoint > &b ) { return a.first < b.first; } );
1282 for (
int vertexIndex = 0; vertexIndex < nVertices; ++vertexIndex )
1286 if ( pointMap.contains( vertexIndex ) )
1289 for (
int k = 0; k < pointMap[ vertexIndex ].size(); ++k )
1291 splitPoint = pointMap[ vertexIndex ][k].second;
1292 if ( splitPoint == currentPoint )
1298 else if ( splitPoint == line->
pointN( vertexIndex + 1 ) )
1317 return asGeos( &lines, mPrecision );
1322 Q_UNUSED( skipIntersectionCheck )
1331 geos::unique_ptr intersectGeom( GEOSIntersection_r( context, splitLine, mGeos.get() ) );
1332 if ( !intersectGeom || GEOSisEmpty_r( context, intersectGeom.get() ) )
1336 const int linearIntersect = GEOSRelatePattern_r( context, mGeos.get(), splitLine,
"1********" );
1337 if ( linearIntersect > 0 )
1340 geos::unique_ptr splitGeom = linePointDifference( intersectGeom.get() );
1345 std::vector<geos::unique_ptr> lineGeoms;
1347 const int splitType = GEOSGeomTypeId_r( context, splitGeom.get() );
1348 if ( splitType == GEOS_MULTILINESTRING )
1350 const int nGeoms = GEOSGetNumGeometries_r( context, splitGeom.get() );
1351 lineGeoms.reserve( nGeoms );
1352 for (
int i = 0; i < nGeoms; ++i )
1353 lineGeoms.emplace_back( GEOSGeom_clone_r( context, GEOSGetGeometryN_r( context, splitGeom.get(), i ) ) );
1358 lineGeoms.emplace_back( GEOSGeom_clone_r( context, splitGeom.get() ) );
1361 mergeGeometriesMultiTypeSplit( lineGeoms );
1363 for ( geos::unique_ptr &lineGeom : lineGeoms )
1381 if ( !mGeosPrepared )
1387 if ( !skipIntersectionCheck && !GEOSPreparedIntersects_r( context, mGeosPrepared.get(), splitLine ) )
1391 geos::unique_ptr nodedGeometry = nodeGeometries( splitLine, mGeos.get() );
1392 if ( !nodedGeometry )
1396 geos::unique_ptr polygons( GEOSPolygonize_r( context, &noded, 1 ) );
1401 const int numberOfGeometriesPolygon = numberOfGeometries( polygons.get() );
1402 if ( numberOfGeometriesPolygon == 0 )
1409 std::vector<geos::unique_ptr> testedGeometries;
1414 for (
int i = 0; i < numberOfGeometriesPolygon; i++ )
1416 const GEOSGeometry *polygon = GEOSGetGeometryN_r( context, polygons.get(), i );
1418 geos::unique_ptr
pointOnSurface( GEOSPointOnSurface_r( context, polygon ) );
1420 testedGeometries.emplace_back( GEOSGeom_clone_r( context, polygon ) );
1423 const size_t nGeometriesThis = numberOfGeometries( mGeos.get() );
1424 if ( testedGeometries.empty() || testedGeometries.size() == nGeometriesThis )
1434 mergeGeometriesMultiTypeSplit( testedGeometries );
1437 for ( i = 0; i < testedGeometries.size() && GEOSisValid_r( context, testedGeometries[i].get() ); ++i )
1440 if ( i < testedGeometries.size() )
1445 for ( geos::unique_ptr &testedGeometry : testedGeometries )
1455 if ( !splitLine || !geom )
1458 geos::unique_ptr geometryBoundary;
1460 if ( GEOSGeomTypeId_r( context, geom ) == GEOS_POLYGON || GEOSGeomTypeId_r( context, geom ) == GEOS_MULTIPOLYGON )
1461 geometryBoundary.reset( GEOSBoundary_r( context, geom ) );
1463 geometryBoundary.reset( GEOSGeom_clone_r( context, geom ) );
1465 geos::unique_ptr splitLineClone( GEOSGeom_clone_r( context, splitLine ) );
1466 geos::unique_ptr unionGeometry( GEOSUnion_r( context, splitLineClone.get(), geometryBoundary.get() ) );
1468 return unionGeometry;
1471int QgsGeos::mergeGeometriesMultiTypeSplit( std::vector<geos::unique_ptr> &splitResult )
const
1478 int type = GEOSGeomTypeId_r( context, mGeos.get() );
1479 if ( type != GEOS_GEOMETRYCOLLECTION &&
1480 type != GEOS_MULTILINESTRING &&
1481 type != GEOS_MULTIPOLYGON &&
1482 type != GEOS_MULTIPOINT )
1486 std::vector<geos::unique_ptr> unionGeom;
1488 std::vector<geos::unique_ptr> newSplitResult;
1490 for (
size_t i = 0; i < splitResult.size(); ++i )
1493 bool isPart =
false;
1494 for (
int j = 0; j < GEOSGetNumGeometries_r( context, mGeos.get() ); j++ )
1496 if ( GEOSEquals_r( context, splitResult[i].get(), GEOSGetGeometryN_r( context, mGeos.get(), j ) ) )
1505 unionGeom.emplace_back( std::move( splitResult[i] ) );
1509 std::vector<geos::unique_ptr> geomVector;
1510 geomVector.emplace_back( std::move( splitResult[i] ) );
1512 if ( type == GEOS_MULTILINESTRING )
1513 newSplitResult.emplace_back( createGeosCollection( GEOS_MULTILINESTRING, geomVector ) );
1514 else if ( type == GEOS_MULTIPOLYGON )
1515 newSplitResult.emplace_back( createGeosCollection( GEOS_MULTIPOLYGON, geomVector ) );
1519 splitResult = std::move( newSplitResult );
1522 if ( !unionGeom.empty() )
1524 if ( type == GEOS_MULTILINESTRING )
1525 splitResult.emplace_back( createGeosCollection( GEOS_MULTILINESTRING, unionGeom ) );
1526 else if ( type == GEOS_MULTIPOLYGON )
1527 splitResult.emplace_back( createGeosCollection( GEOS_MULTIPOLYGON, unionGeom ) );
1533geos::unique_ptr QgsGeos::createGeosCollection(
int typeId, std::vector<geos::unique_ptr> &geoms )
1535 std::vector<GEOSGeometry *> geomarr;
1536 geomarr.reserve( geoms.size() );
1539 for ( geos::unique_ptr &geomUniquePtr : geoms )
1541 if ( geomUniquePtr )
1543 if ( !GEOSisEmpty_r( context, geomUniquePtr.get() ) )
1547 geomarr.emplace_back( geomUniquePtr.release() );
1551 geos::unique_ptr geomRes;
1555 geomRes.reset( GEOSGeom_createCollection_r( context, typeId, geomarr.data(), geomarr.size() ) );
1557 catch ( GEOSException & )
1561 GEOSGeom_destroy_r( context, geom );
1576 int nCoordDims = GEOSGeom_getCoordinateDimension_r( context,
geos );
1577 int nDims = GEOSGeom_getDimensions_r( context,
geos );
1578 bool hasZ = ( nCoordDims == 3 );
1579 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1581 switch ( GEOSGeomTypeId_r( context,
geos ) )
1585 if ( GEOSisEmpty_r( context,
geos ) )
1588 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context,
geos );
1589 unsigned int nPoints = 0;
1590 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1591 return nPoints > 0 ? std::unique_ptr<QgsAbstractGeometry>(
coordSeqPoint( cs, 0, hasZ, hasM ).clone() ) :
nullptr;
1593 case GEOS_LINESTRING:
1595 return sequenceToLinestring(
geos, hasZ, hasM );
1601 case GEOS_MULTIPOINT:
1603 auto multiPoint = std::make_unique<QgsMultiPoint>();
1604 int nParts = GEOSGetNumGeometries_r( context,
geos );
1605 multiPoint->reserve( nParts );
1606 for (
int i = 0; i < nParts; ++i )
1608 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context, GEOSGetGeometryN_r( context,
geos, i ) );
1611 unsigned int nPoints = 0;
1612 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1614 multiPoint->addGeometry(
coordSeqPoint( cs, 0, hasZ, hasM ).clone() );
1617 return std::move( multiPoint );
1619 case GEOS_MULTILINESTRING:
1621 auto multiLineString = std::make_unique<QgsMultiLineString>();
1622 int nParts = GEOSGetNumGeometries_r( context,
geos );
1623 multiLineString->reserve( nParts );
1624 for (
int i = 0; i < nParts; ++i )
1626 std::unique_ptr< QgsLineString >line( sequenceToLinestring( GEOSGetGeometryN_r( context,
geos, i ), hasZ, hasM ) );
1629 multiLineString->addGeometry( line.release() );
1632 return std::move( multiLineString );
1634 case GEOS_MULTIPOLYGON:
1636 auto multiPolygon = std::make_unique<QgsMultiPolygon>();
1638 int nParts = GEOSGetNumGeometries_r( context,
geos );
1639 multiPolygon->reserve( nParts );
1640 for (
int i = 0; i < nParts; ++i )
1642 std::unique_ptr< QgsPolygon > poly =
fromGeosPolygon( GEOSGetGeometryN_r( context,
geos, i ) );
1645 multiPolygon->addGeometry( poly.release() );
1648 return std::move( multiPolygon );
1650 case GEOS_GEOMETRYCOLLECTION:
1652 auto geomCollection = std::make_unique<QgsGeometryCollection>();
1653 int nParts = GEOSGetNumGeometries_r( context,
geos );
1654 geomCollection->reserve( nParts );
1655 for (
int i = 0; i < nParts; ++i )
1657 std::unique_ptr< QgsAbstractGeometry > geom(
fromGeos( GEOSGetGeometryN_r( context,
geos, i ) ) );
1660 geomCollection->addGeometry( geom.release() );
1663 return std::move( geomCollection );
1672 if ( GEOSGeomTypeId_r( context,
geos ) != GEOS_POLYGON )
1677 int nCoordDims = GEOSGeom_getCoordinateDimension_r( context,
geos );
1678 int nDims = GEOSGeom_getDimensions_r( context,
geos );
1679 bool hasZ = ( nCoordDims == 3 );
1680 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1682 auto polygon = std::make_unique<QgsPolygon>();
1687 polygon->setExteriorRing( sequenceToLinestring( ring, hasZ, hasM ).release() );
1690 QVector<QgsCurve *> interiorRings;
1691 const int ringCount = GEOSGetNumInteriorRings_r( context,
geos );
1692 interiorRings.reserve( ringCount );
1693 for (
int i = 0; i < ringCount; ++i )
1695 ring = GEOSGetInteriorRingN_r( context,
geos, i );
1698 interiorRings.push_back( sequenceToLinestring( ring, hasZ, hasM ).release() );
1701 polygon->setInteriorRings( interiorRings );
1706std::unique_ptr<QgsLineString> QgsGeos::sequenceToLinestring(
const GEOSGeometry *
geos,
bool hasZ,
bool hasM )
1709 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context,
geos );
1711 unsigned int nPoints;
1712 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1714 QVector< double > xOut( nPoints );
1715 QVector< double > yOut( nPoints );
1716 QVector< double > zOut;
1718 zOut.resize( nPoints );
1719 QVector< double > mOut;
1721 mOut.resize( nPoints );
1723 double *x = xOut.data();
1724 double *y = yOut.data();
1725 double *z = zOut.data();
1726 double *m = mOut.data();
1728#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
1729 GEOSCoordSeq_copyToArrays_r( context, cs, x, y, hasZ ? z : nullptr, hasM ? m : nullptr );
1731 for (
unsigned int i = 0; i < nPoints; ++i )
1734 GEOSCoordSeq_getXYZ_r( context, cs, i, x++, y++, z++ );
1736 GEOSCoordSeq_getXY_r( context, cs, i, x++, y++ );
1739 GEOSCoordSeq_getOrdinate_r( context, cs, i, 3, m++ );
1743 auto line = std::make_unique<QgsLineString>( xOut, yOut, zOut, mOut );
1753 int geometryType = GEOSGeomTypeId_r( context, g );
1754 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
1755 || geometryType == GEOS_POLYGON )
1759 return GEOSGetNumGeometries_r( context, g );
1775 GEOSCoordSeq_getXYZ_r( context, cs, i, &x, &y, &z );
1777 GEOSCoordSeq_getXY_r( context, cs, i, &x, &y );
1780 GEOSCoordSeq_getOrdinate_r( context, cs, i, 3, &m );
1816 int geosType = GEOS_GEOMETRYCOLLECTION;
1823 geosType = GEOS_MULTIPOINT;
1827 geosType = GEOS_MULTILINESTRING;
1831 geosType = GEOS_MULTIPOLYGON;
1846 std::vector<geos::unique_ptr> geomVector;
1847 geomVector.reserve(
c->numGeometries() );
1848 for (
int i = 0; i <
c->numGeometries(); ++i )
1855 geomVector.emplace_back( std::move( geosGeom ) );
1857 return createGeosCollection( geosType, geomVector );
1864 const QgsPolyhedralSurface *polyhedralSurface = qgsgeometry_cast<const QgsPolyhedralSurface *>( geom );
1865 if ( !polyhedralSurface )
1868 std::vector<geos::unique_ptr> geomVector;
1869 geomVector.reserve( polyhedralSurface->
numPatches() );
1870 for (
int i = 0; i < polyhedralSurface->
numPatches(); ++i )
1872 geos::unique_ptr geosPolygon = createGeosPolygon( polyhedralSurface->
patchN( i ),
precision );
1877 geomVector.emplace_back( std::move( geosPolygon ) );
1880 return createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
1887 return createGeosPoint(
static_cast<const QgsPoint *
>( geom ), coordDims,
precision, flags );
1905 if ( !mGeos || !geom )
1910 geos::unique_ptr geosGeom(
asGeos( geom, mPrecision ) );
1916 const double gridSize = parameters.
gridSize();
1921 geos::unique_ptr opGeom;
1924 case OverlayIntersection:
1927 opGeom.reset( GEOSIntersectionPrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1931 opGeom.reset( GEOSIntersection_r( context, mGeos.get(), geosGeom.get() ) );
1935 case OverlayDifference:
1938 opGeom.reset( GEOSDifferencePrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1942 opGeom.reset( GEOSDifference_r( context, mGeos.get(), geosGeom.get() ) );
1948 geos::unique_ptr unionGeometry;
1951 unionGeometry.reset( GEOSUnionPrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1955 unionGeometry.reset( GEOSUnion_r( context, mGeos.get(), geosGeom.get() ) );
1958 if ( unionGeometry && GEOSGeomTypeId_r( context, unionGeometry.get() ) == GEOS_MULTILINESTRING )
1960 geos::unique_ptr mergedLines( GEOSLineMerge_r( context, unionGeometry.get() ) );
1963 unionGeometry = std::move( mergedLines );
1967 opGeom = std::move( unionGeometry );
1971 case OverlaySymDifference:
1974 opGeom.reset( GEOSSymDifferencePrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1978 opGeom.reset( GEOSSymDifference_r( context, mGeos.get(), geosGeom.get() ) );
1984 catch ( GEOSException &e )
1986 logError( QStringLiteral(
"GEOS" ), e.what() );
1989 *errorMsg = e.what();
1995bool QgsGeos::relation(
const QgsAbstractGeometry *geom, Relation r, QString *errorMsg )
const
1997 if ( !mGeos || !geom )
2002 geos::unique_ptr geosGeom(
asGeos( geom, mPrecision ) );
2009 bool result =
false;
2012 if ( mGeosPrepared )
2016 case RelationIntersects:
2017 result = ( GEOSPreparedIntersects_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2019 case RelationTouches:
2020 result = ( GEOSPreparedTouches_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2022 case RelationCrosses:
2023 result = ( GEOSPreparedCrosses_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2025 case RelationWithin:
2026 result = ( GEOSPreparedWithin_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2028 case RelationContains:
2029 result = ( GEOSPreparedContains_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2031 case RelationDisjoint:
2032 result = ( GEOSPreparedDisjoint_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2034 case RelationOverlaps:
2035 result = ( GEOSPreparedOverlaps_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2043 case RelationIntersects:
2044 result = ( GEOSIntersects_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2046 case RelationTouches:
2047 result = ( GEOSTouches_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2049 case RelationCrosses:
2050 result = ( GEOSCrosses_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2052 case RelationWithin:
2053 result = ( GEOSWithin_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2055 case RelationContains:
2056 result = ( GEOSContains_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2058 case RelationDisjoint:
2059 result = ( GEOSDisjoint_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2061 case RelationOverlaps:
2062 result = ( GEOSOverlaps_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2066 catch ( GEOSException &e )
2068 logError( QStringLiteral(
"GEOS" ), e.what() );
2071 *errorMsg = e.what();
2086 geos::unique_ptr
geos;
2097 geos::unique_ptr
geos =
buffer( mGeos.get(),
distance, segments, endCapStyle, joinStyle, miterLimit, errorMsg );
2108 geos::unique_ptr
geos;
2111 geos.reset( GEOSBufferWithStyle_r(
QgsGeosContext::get(), geometry,
distance, segments,
static_cast< int >( endCapStyle ),
static_cast< int >( joinStyle ), miterLimit ) );
2123 geos::unique_ptr
geos;
2138 geos::unique_ptr
geos;
2154 geos::unique_ptr
geos;
2161 geos.reset( GEOSGetCentroid_r( context, mGeos.get() ) );
2166 GEOSGeomGetX_r( context,
geos.get(), &x );
2167 GEOSGeomGetY_r( context,
geos.get(), &y );
2180 geos::unique_ptr
geos;
2200 geos::unique_ptr
geos;
2203 geos.reset( GEOSPointOnSurface_r( context, mGeos.get() ) );
2205 if ( !
geos || GEOSisEmpty_r( context,
geos.get() ) != 0 )
2210 GEOSGeomGetX_r( context,
geos.get(), &x );
2211 GEOSGeomGetY_r( context,
geos.get(), &y );
2228 std::unique_ptr< QgsAbstractGeometry > cHullGeom =
fromGeos( cHull.get() );
2229 return cHullGeom.release();
2234std::unique_ptr< QgsAbstractGeometry >
QgsGeos::concaveHull(
double targetPercent,
bool allowHoles, QString *errorMsg )
const
2236#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
2238 ( void )targetPercent;
2240 throw QgsNotSupportedException( QObject::tr(
"Calculating concaveHull requires a QGIS build based on GEOS 3.11 or later" ) );
2251 return concaveHullGeom;
2259#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<12
2261 ( void )invalidEdges;
2263 throw QgsNotSupportedException( QObject::tr(
"Validating coverages requires a QGIS build based on GEOS 3.12 or later" ) );
2268 *errorMsg = QStringLiteral(
"Input geometry was not set" );
2276 const int result = GEOSCoverageIsValid_r( context, mGeos.get(), gapWidth, invalidEdges ? &invalidEdgesGeos : nullptr );
2277 if ( invalidEdges && invalidEdgesGeos )
2279 *invalidEdges =
fromGeos( invalidEdgesGeos );
2281 if ( invalidEdgesGeos )
2283 GEOSGeom_destroy_r( context, invalidEdgesGeos );
2284 invalidEdgesGeos =
nullptr;
2304#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<12
2306 ( void )preserveBoundary;
2308 throw QgsNotSupportedException( QObject::tr(
"Simplifying coverages requires a QGIS build based on GEOS 3.12 or later" ) );
2313 *errorMsg = QStringLiteral(
"Input geometry was not set" );
2319 geos::unique_ptr simplified( GEOSCoverageSimplifyVW_r(
QgsGeosContext::get(), mGeos.get(), tolerance, preserveBoundary ? 1 : 0 ) );
2320 std::unique_ptr< QgsAbstractGeometry > simplifiedGeom =
fromGeos( simplified.get() );
2321 return simplifiedGeom;
2332 *errorMsg = QStringLiteral(
"Input geometry was not set" );
2339 std::unique_ptr< QgsAbstractGeometry > result =
fromGeos( unioned.get() );
2350 *errorMsg = QObject::tr(
"QGIS geometry cannot be converted to a GEOS geometry",
"GEOS Error" );
2359 char res = GEOSisValidDetail_r( context, mGeos.get(), allowSelfTouchingHoles ? GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE : 0, &r, &g1 );
2360 const bool invalid = res != 1;
2365 error = QString( r );
2366 GEOSFree_r( context, r );
2369 if ( invalid && errorMsg )
2372 static const std::map< QString, QString > sTranslatedErrors
2374 { QStringLiteral(
"topology validation error" ), QObject::tr(
"Topology validation error",
"GEOS Error" ) },
2375 { QStringLiteral(
"repeated point" ), QObject::tr(
"Repeated point",
"GEOS Error" ) },
2376 { QStringLiteral(
"hole lies outside shell" ), QObject::tr(
"Hole lies outside shell",
"GEOS Error" ) },
2377 { QStringLiteral(
"holes are nested" ), QObject::tr(
"Holes are nested",
"GEOS Error" ) },
2378 { QStringLiteral(
"interior is disconnected" ), QObject::tr(
"Interior is disconnected",
"GEOS Error" ) },
2379 { QStringLiteral(
"self-intersection" ), QObject::tr(
"Self-intersection",
"GEOS Error" ) },
2380 { QStringLiteral(
"ring self-intersection" ), QObject::tr(
"Ring self-intersection",
"GEOS Error" ) },
2381 { QStringLiteral(
"nested shells" ), QObject::tr(
"Nested shells",
"GEOS Error" ) },
2382 { QStringLiteral(
"duplicate rings" ), QObject::tr(
"Duplicate rings",
"GEOS Error" ) },
2383 { QStringLiteral(
"too few points in geometry component" ), QObject::tr(
"Too few points in geometry component",
"GEOS Error" ) },
2384 { QStringLiteral(
"invalid coordinate" ), QObject::tr(
"Invalid coordinate",
"GEOS Error" ) },
2385 { QStringLiteral(
"ring is not closed" ), QObject::tr(
"Ring is not closed",
"GEOS Error" ) },
2388 const auto translatedError = sTranslatedErrors.find( error.toLower() );
2389 if ( translatedError != sTranslatedErrors.end() )
2390 *errorMsg = translatedError->second;
2394 if ( g1 && errorLoc )
2400 GEOSGeom_destroy_r( context, g1 );
2410 if ( !mGeos || !geom )
2417 geos::unique_ptr geosGeom(
asGeos( geom, mPrecision ) );
2456GEOSCoordSequence *QgsGeos::createCoordinateSequence(
const QgsCurve *curve,
double precision,
bool forceClose )
2460 std::unique_ptr< QgsLineString > segmentized;
2461 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( curve );
2466 line = segmentized.get();
2473 GEOSCoordSequence *coordSeq =
nullptr;
2475 const int numPoints = line->
numPoints();
2477 const bool hasZ = line->
is3D();
2479#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
2482 if ( !forceClose || ( line->
pointN( 0 ) == line->
pointN( numPoints - 1 ) ) )
2487 coordSeq = GEOSCoordSeq_copyFromArrays_r( context, line->
xData(), line->
yData(), line->
zData(),
nullptr, numPoints );
2490 QgsDebugError( QStringLiteral(
"GEOS Exception: Could not create coordinate sequence for %1 points" ).arg( numPoints ) );
2498 QVector< double > x = line->
xVector();
2499 if ( numPoints > 0 )
2500 x.append( x.at( 0 ) );
2501 QVector< double > y = line->
yVector();
2502 if ( numPoints > 0 )
2503 y.append( y.at( 0 ) );
2504 QVector< double > z = line->
zVector();
2505 if ( hasZ && numPoints > 0 )
2506 z.append( z.at( 0 ) );
2509 coordSeq = GEOSCoordSeq_copyFromArrays_r( context, x.constData(), y.constData(), !hasZ ?
nullptr : z.constData(), nullptr, numPoints + 1 );
2512 QgsDebugError( QStringLiteral(
"GEOS Exception: Could not create closed coordinate sequence for %1 points" ).arg( numPoints + 1 ) );
2523 const bool hasM =
false;
2534 int numOutPoints = numPoints;
2535 if ( forceClose && ( line->
pointN( 0 ) != line->
pointN( numPoints - 1 ) ) )
2542 coordSeq = GEOSCoordSeq_create_r( context, numOutPoints, coordDims );
2545 QgsDebugError( QStringLiteral(
"GEOS Exception: Could not create coordinate sequence for %1 points in %2 dimensions" ).arg( numPoints ).arg( coordDims ) );
2549 const double *xData = line->
xData();
2550 const double *yData = line->
yData();
2551 const double *zData = hasZ ? line->
zData() :
nullptr;
2552 const double *mData = hasM ? line->
mData() :
nullptr;
2556 for (
int i = 0; i < numOutPoints; ++i )
2558 if ( i >= numPoints )
2561 xData = line->
xData();
2562 yData = line->
yData();
2563 zData = hasZ ? line->
zData() :
nullptr;
2564 mData = hasM ? line->
mData() :
nullptr;
2576 GEOSCoordSeq_setOrdinate_r( context, coordSeq, i, 3, *mData++ );
2582 for (
int i = 0; i < numOutPoints; ++i )
2584 if ( i >= numPoints )
2587 xData = line->
xData();
2588 yData = line->
yData();
2589 zData = hasZ ? line->
zData() :
nullptr;
2590 mData = hasM ? line->
mData() :
nullptr;
2594 GEOSCoordSeq_setXYZ_r( context, coordSeq, i, *xData++, *yData++, *zData++ );
2598 GEOSCoordSeq_setXY_r( context, coordSeq, i, *xData++, *yData++ );
2602 GEOSCoordSeq_setOrdinate_r( context, coordSeq, i, 3, *mData++ );
2614 const QgsPoint *pt = qgsgeometry_cast<const QgsPoint *>( point );
2621geos::unique_ptr QgsGeos::createGeosPointXY(
double x,
double y,
bool hasZ,
double z,
bool hasM,
double m,
int coordDims,
double precision,
Qgis::GeosCreationFlags )
2626 geos::unique_ptr geosPoint;
2630 if ( coordDims == 2 )
2636 geosPoint.reset( GEOSGeom_createPointFromXY_r( context, x, y ) );
2640 GEOSCoordSequence *coordSeq = GEOSCoordSeq_create_r( context, 1, coordDims );
2643 QgsDebugError( QStringLiteral(
"GEOS Exception: Could not create coordinate sequence for point with %1 dimensions" ).arg( coordDims ) );
2648 GEOSCoordSeq_setX_r( context, coordSeq, 0, std::round( x /
precision ) *
precision );
2649 GEOSCoordSeq_setY_r( context, coordSeq, 0, std::round( y /
precision ) *
precision );
2652 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 2, std::round( z /
precision ) *
precision );
2657 GEOSCoordSeq_setX_r( context, coordSeq, 0, x );
2658 GEOSCoordSeq_setY_r( context, coordSeq, 0, y );
2661 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 2, z );
2667 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 3, m );
2670 geosPoint.reset( GEOSGeom_createPoint_r( context, coordSeq ) );
2678 const QgsCurve *
c = qgsgeometry_cast<const QgsCurve *>( curve );
2682 GEOSCoordSequence *coordSeq = createCoordinateSequence(
c,
precision );
2686 geos::unique_ptr geosGeom;
2697 const QgsCurvePolygon *polygon = qgsgeometry_cast<const QgsCurvePolygon *>( poly );
2702 if ( !exteriorRing )
2708 geos::unique_ptr geosPolygon;
2711 geos::unique_ptr exteriorRingGeos( GEOSGeom_createLinearRing_r( context, createCoordinateSequence( exteriorRing,
precision,
true ) ) );
2717 for (
int i = 0; i < nInteriorRings; ++i )
2720 if ( !interiorRing->
isEmpty() )
2728 nHoles = nInteriorRings;
2736 for (
int i = 0; i < nInteriorRings; ++i )
2741 holes[i] = GEOSGeom_createLinearRing_r( context, createCoordinateSequence( interiorRing,
precision,
true ) );
2744 geosPolygon.reset( GEOSGeom_createPolygon_r( context, exteriorRingGeos.release(), holes, nHoles ) );
2757 geos::unique_ptr offset;
2764 offset.reset( GEOSOffsetCurve_r(
QgsGeosContext::get(), geometry,
distance, segments,
static_cast< int >( joinStyle ), miterLimit ) );
2772 geos::unique_ptr res =
offsetCurve( mGeos.get(),
distance, segments, joinStyle, miterLimit, errorMsg );
2776 return fromGeos( res.get() ).release();
2786 geos::unique_ptr
geos;
2790 geos::buffer_params_unique_ptr bp( GEOSBufferParams_create_r( context ) );
2791 GEOSBufferParams_setSingleSided_r( context, bp.get(), 1 );
2792 GEOSBufferParams_setQuadrantSegments_r( context, bp.get(), segments );
2793 GEOSBufferParams_setJoinStyle_r( context, bp.get(),
static_cast< int >( joinStyle ) );
2794 GEOSBufferParams_setMitreLimit_r( context, bp.get(), miterLimit );
2800 geos.reset( GEOSBufferWithParams_r( context, mGeos.get(), bp.get(),
distance ) );
2813 geos::unique_ptr
geos;
2829 geos::unique_ptr
geos;
2832 geos::unique_ptr boundaryGeos;
2834 boundaryGeos =
asGeos( boundary );
2849 geos::unique_ptr
geos;
2862 return std::numeric_limits< double >::quiet_NaN();
2865 geos::unique_ptr
geos;
2870 return std::numeric_limits< double >::quiet_NaN();
2883 geos::unique_ptr
geos;
2899 geos::unique_ptr
geos;
2910 if ( !mGeos || !other )
2915 geos::unique_ptr
geos;
2918 geos::unique_ptr otherGeos =
asGeos( other );
2942 geos::unique_ptr reshapeLineGeos = createGeosLinestring( &reshapeWithLine, mPrecision );
2946 int numGeoms = GEOSGetNumGeometries_r( context, mGeos.get() );
2947 if ( numGeoms == -1 )
2956 bool isMultiGeom =
false;
2957 int geosTypeId = GEOSGeomTypeId_r( context, mGeos.get() );
2958 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
2965 geos::unique_ptr reshapedGeometry;
2968 reshapedGeometry = reshapeLine( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2972 reshapedGeometry = reshapePolygon( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2977 std::unique_ptr< QgsAbstractGeometry > reshapeResult =
fromGeos( reshapedGeometry.get() );
2978 return reshapeResult;
2985 bool reshapeTookPlace =
false;
2987 geos::unique_ptr currentReshapeGeometry;
2990 for (
int i = 0; i < numGeoms; ++i )
2993 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( context, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2995 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( context, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2997 if ( currentReshapeGeometry )
2999 newGeoms[i] = currentReshapeGeometry.release();
3000 reshapeTookPlace =
true;
3004 newGeoms[i] = GEOSGeom_clone_r( context, GEOSGetGeometryN_r( context, mGeos.get(), i ) );
3008 geos::unique_ptr newMultiGeom;
3011 newMultiGeom.reset( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, newGeoms, numGeoms ) );
3015 newMultiGeom.reset( GEOSGeom_createCollection_r( context, GEOS_MULTIPOLYGON, newGeoms, numGeoms ) );
3019 if ( !newMultiGeom )
3025 if ( reshapeTookPlace )
3029 std::unique_ptr< QgsAbstractGeometry > reshapedMultiGeom =
fromGeos( newMultiGeom.get() );
3030 return reshapedMultiGeom;
3053 if ( GEOSGeomTypeId_r( context, mGeos.get() ) != GEOS_MULTILINESTRING )
3056 geos::unique_ptr
geos;
3059 geos.reset( GEOSLineMerge_r( context, mGeos.get() ) );
3072 geos::unique_ptr otherGeom(
asGeos( other.
constGet(), mPrecision ) );
3083 geos::coord_sequence_unique_ptr nearestCoord;
3084 if ( mGeosPrepared )
3086 nearestCoord.reset( GEOSPreparedNearestPoints_r( context, mGeosPrepared.get(), otherGeom.get() ) );
3090 nearestCoord.reset( GEOSNearestPoints_r( context, mGeos.get(), otherGeom.get() ) );
3093 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 0, &nx );
3094 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 0, &ny );
3096 catch ( GEOSException &e )
3098 logError( QStringLiteral(
"GEOS" ), e.what() );
3101 *errorMsg = e.what();
3106 return std::make_unique< QgsPoint >( nx, ny );
3111 if ( !mGeos || other.
isEmpty() )
3121 if ( !other || other->
isEmpty() )
3124 geos::unique_ptr otherGeom(
asGeos( other, mPrecision ) );
3137 geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( context, mGeos.get(), otherGeom.get() ) );
3139 if ( !nearestCoord )
3142 *errorMsg = QStringLiteral(
"GEOS returned no nearest points" );
3146 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 0, &nx1 );
3147 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 0, &ny1 );
3148 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 1, &nx2 );
3149 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 1, &ny2 );
3151 catch ( GEOSException &e )
3153 logError( QStringLiteral(
"GEOS" ), e.what() );
3156 *errorMsg = e.what();
3161 auto line = std::make_unique< QgsLineString >();
3174 geos::unique_ptr otherGeom(
asGeos( &point, mPrecision ) );
3185 catch ( GEOSException &e )
3187 logError( QStringLiteral(
"GEOS" ), e.what() );
3190 *errorMsg = e.what();
3205 geos::unique_ptr point = createGeosPointXY( x, y,
false, 0,
false, 0, 2, 0 );
3214 catch ( GEOSException &e )
3216 logError( QStringLiteral(
"GEOS" ), e.what() );
3219 *errorMsg = e.what();
3233 geos::unique_ptr l =
asGeos( g );
3236 lineGeosGeometries[validLines] = l.release();
3244 geos::unique_ptr result( GEOSPolygonize_r( context, lineGeosGeometries, validLines ) );
3245 for (
int i = 0; i < validLines; ++i )
3247 GEOSGeom_destroy_r( context, lineGeosGeometries[i] );
3249 delete[] lineGeosGeometries;
3252 catch ( GEOSException &e )
3256 *errorMsg = e.what();
3258 for (
int i = 0; i < validLines; ++i )
3260 GEOSGeom_destroy_r( context, lineGeosGeometries[i] );
3262 delete[] lineGeosGeometries;
3274 geos::unique_ptr extentGeosGeom;
3277 extentGeosGeom =
asGeos( extent, mPrecision );
3278 if ( !extentGeosGeom )
3284 geos::unique_ptr
geos;
3288 geos.reset( GEOSVoronoiDiagram_r( context, mGeos.get(), extentGeosGeom.get(), tolerance, edgesOnly ) );
3290 if ( !
geos || GEOSisEmpty_r( context,
geos.get() ) != 0 )
3308 geos::unique_ptr
geos;
3311 geos.reset( GEOSDelaunayTriangulation_r( context, mGeos.get(), tolerance, edgesOnly ) );
3313 if ( !
geos || GEOSisEmpty_r( context,
geos.get() ) != 0 )
3325#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
3327 throw QgsNotSupportedException( QObject::tr(
"Calculating constrainedDelaunayTriangulation requires a QGIS build based on GEOS 3.11 or later" ) );
3334 geos::unique_ptr
geos;
3338 geos.reset( GEOSConstrainedDelaunayTriangulation_r( context, mGeos.get() ) );
3340 if ( !
geos || GEOSisEmpty_r( context,
geos.get() ) != 0 )
3345 std::unique_ptr< QgsAbstractGeometry > res =
fromGeos(
geos.get() );
3346 if (
const QgsGeometryCollection *collection = qgsgeometry_cast< const QgsGeometryCollection * >( res.get() ) )
3348 return std::unique_ptr< QgsAbstractGeometry >( collection->extractPartsByType(
Qgis::WkbType::Polygon,
true ) );
3360static bool _linestringEndpoints(
const GEOSGeometry *linestring,
double &x1,
double &y1,
double &x2,
double &y2 )
3363 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( context, linestring );
3367 unsigned int coordSeqSize;
3368 if ( GEOSCoordSeq_getSize_r( context, coordSeq, &coordSeqSize ) == 0 )
3371 if ( coordSeqSize < 2 )
3374 GEOSCoordSeq_getX_r( context, coordSeq, 0, &x1 );
3375 GEOSCoordSeq_getY_r( context, coordSeq, 0, &y1 );
3376 GEOSCoordSeq_getX_r( context, coordSeq, coordSeqSize - 1, &x2 );
3377 GEOSCoordSeq_getY_r( context, coordSeq, coordSeqSize - 1, &y2 );
3385 double x1, y1, x2, y2;
3386 if ( !_linestringEndpoints( line1, x1, y1, x2, y2 ) )
3389 double rx1, ry1, rx2, ry2;
3390 if ( !_linestringEndpoints( line2, rx1, ry1, rx2, ry2 ) )
3393 bool intersectionAtOrigLineEndpoint =
3394 ( intersectionPoint.
x() == x1 && intersectionPoint.
y() == y1 ) !=
3395 ( intersectionPoint.
x() == x2 && intersectionPoint.
y() == y2 );
3396 bool intersectionAtReshapeLineEndpoint =
3397 ( intersectionPoint.
x() == rx1 && intersectionPoint.
y() == ry1 ) ||
3398 ( intersectionPoint.
x() == rx2 && intersectionPoint.
y() == ry2 );
3402 if ( intersectionAtOrigLineEndpoint && intersectionAtReshapeLineEndpoint )
3404 geos::unique_ptr g1( GEOSGeom_clone_r( context, line1 ) );
3405 geos::unique_ptr g2( GEOSGeom_clone_r( context, line2 ) );
3406 GEOSGeometry *geoms[2] = { g1.release(), g2.release() };
3407 geos::unique_ptr multiGeom( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, geoms, 2 ) );
3408 geos::unique_ptr res( GEOSLineMerge_r( context, multiGeom.get() ) );
3412 double x1res, y1res, x2res, y2res;
3413 if ( !_linestringEndpoints( res.get(), x1res, y1res, x2res, y2res ) )
3415 if ( ( x1res == x2 && y1res == y2 ) || ( x2res == x1 && y2res == y1 ) )
3416 res.reset( GEOSReverse_r( context, res.get() ) );
3427 if ( !line || !reshapeLineGeos )
3430 bool atLeastTwoIntersections =
false;
3431 bool oneIntersection =
false;
3438 geos::unique_ptr intersectGeom( GEOSIntersection_r( context, line, reshapeLineGeos ) );
3439 if ( intersectGeom )
3441 const int geomType = GEOSGeomTypeId_r( context, intersectGeom.get() );
3442 atLeastTwoIntersections = ( geomType == GEOS_MULTIPOINT && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 1 )
3443 || ( geomType == GEOS_GEOMETRYCOLLECTION && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 0 )
3444 || ( geomType == GEOS_MULTILINESTRING && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 0 );
3446 if ( GEOSGeomTypeId_r( context, intersectGeom.get() ) == GEOS_POINT )
3448 const GEOSCoordSequence *intersectionCoordSeq = GEOSGeom_getCoordSeq_r( context, intersectGeom.get() );
3450 GEOSCoordSeq_getX_r( context, intersectionCoordSeq, 0, &xi );
3451 GEOSCoordSeq_getY_r( context, intersectionCoordSeq, 0, &yi );
3452 oneIntersection =
true;
3457 catch ( GEOSException & )
3459 atLeastTwoIntersections =
false;
3463 if ( oneIntersection )
3464 return _mergeLinestrings( line, reshapeLineGeos, oneIntersectionPoint );
3466 if ( !atLeastTwoIntersections )
3470 double x1, y1, x2, y2;
3471 if ( !_linestringEndpoints( line, x1, y1, x2, y2 ) )
3474 geos::unique_ptr beginLineVertex = createGeosPointXY( x1, y1,
false, 0,
false, 0, 2,
precision );
3475 geos::unique_ptr endLineVertex = createGeosPointXY( x2, y2,
false, 0,
false, 0, 2,
precision );
3477 bool isRing =
false;
3478 if ( GEOSGeomTypeId_r( context, line ) == GEOS_LINEARRING
3479 || GEOSEquals_r( context, beginLineVertex.get(), endLineVertex.get() ) == 1 )
3483 geos::unique_ptr nodedGeometry = nodeGeometries( reshapeLineGeos, line );
3484 if ( !nodedGeometry )
3490 geos::unique_ptr mergedLines( GEOSLineMerge_r( context, nodedGeometry.get() ) );
3496 int numMergedLines = GEOSGetNumGeometries_r( context, mergedLines.get() );
3497 if ( numMergedLines < 2 )
3499 if ( numMergedLines == 1 )
3501 geos::unique_ptr result( GEOSGeom_clone_r( context, reshapeLineGeos ) );
3508 QVector<GEOSGeometry *> resultLineParts;
3509 QVector<GEOSGeometry *> probableParts;
3511 for (
int i = 0; i < numMergedLines; ++i )
3513 const GEOSGeometry *currentGeom = GEOSGetGeometryN_r( context, mergedLines.get(), i );
3516 bool alreadyAdded =
false;
3518 double bufferDistance = std::pow( 10.0L, geomDigits( currentGeom ) - 11 );
3519 for (
const GEOSGeometry *other : std::as_const( resultLineParts ) )
3521 GEOSHausdorffDistance_r( context, currentGeom, other, &
distance );
3524 alreadyAdded =
true;
3531 const GEOSCoordSequence *currentCoordSeq = GEOSGeom_getCoordSeq_r( context, currentGeom );
3532 unsigned int currentCoordSeqSize;
3533 GEOSCoordSeq_getSize_r( context, currentCoordSeq, ¤tCoordSeqSize );
3534 if ( currentCoordSeqSize < 2 )
3538 double xBegin, xEnd, yBegin, yEnd;
3539 GEOSCoordSeq_getX_r( context, currentCoordSeq, 0, &xBegin );
3540 GEOSCoordSeq_getY_r( context, currentCoordSeq, 0, &yBegin );
3541 GEOSCoordSeq_getX_r( context, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
3542 GEOSCoordSeq_getY_r( context, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
3543 geos::unique_ptr beginCurrentGeomVertex = createGeosPointXY( xBegin, yBegin,
false, 0,
false, 0, 2,
precision );
3544 geos::unique_ptr endCurrentGeomVertex = createGeosPointXY( xEnd, yEnd,
false, 0,
false, 0, 2,
precision );
3547 int nEndpointsOnOriginalLine = 0;
3548 if ( pointContainedInLine( beginCurrentGeomVertex.get(), line ) == 1 )
3549 nEndpointsOnOriginalLine += 1;
3551 if ( pointContainedInLine( endCurrentGeomVertex.get(), line ) == 1 )
3552 nEndpointsOnOriginalLine += 1;
3555 int nEndpointsSameAsOriginalLine = 0;
3556 if ( GEOSEquals_r( context, beginCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3557 || GEOSEquals_r( context, beginCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3558 nEndpointsSameAsOriginalLine += 1;
3560 if ( GEOSEquals_r( context, endCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3561 || GEOSEquals_r( context, endCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3562 nEndpointsSameAsOriginalLine += 1;
3565 bool currentGeomOverlapsOriginalGeom =
false;
3566 bool currentGeomOverlapsReshapeLine =
false;
3567 if ( lineContainedInLine( currentGeom, line ) == 1 )
3568 currentGeomOverlapsOriginalGeom =
true;
3570 if ( lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
3571 currentGeomOverlapsReshapeLine =
true;
3574 if ( !isRing && nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3576 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3579 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3581 probableParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3583 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3585 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3587 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3589 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3591 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
3593 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3598 if ( isRing && !probableParts.isEmpty() )
3600 geos::unique_ptr maxGeom;
3602 double maxLength = -std::numeric_limits<double>::max();
3603 double currentLength = 0;
3604 for (
int i = 0; i < probableParts.size(); ++i )
3606 currentGeom = probableParts.at( i );
3607 GEOSLength_r( context, currentGeom, ¤tLength );
3608 if ( currentLength > maxLength )
3610 maxLength = currentLength;
3611 maxGeom.reset( currentGeom );
3615 GEOSGeom_destroy_r( context, currentGeom );
3618 resultLineParts.push_back( maxGeom.release() );
3621 geos::unique_ptr result;
3622 if ( resultLineParts.empty() )
3625 if ( resultLineParts.size() == 1 )
3627 result.reset( resultLineParts[0] );
3632 for (
int i = 0; i < resultLineParts.size(); ++i )
3634 lineArray[i] = resultLineParts[i];
3638 geos::unique_ptr multiLineGeom( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ) );
3639 delete [] lineArray;
3642 result.reset( GEOSLineMerge_r( context, multiLineGeom.get() ) );
3646 if ( GEOSGeomTypeId_r( context, result.get() ) != GEOS_LINESTRING )
3652 bool reverseLine =
false;
3656 char isResultCCW = 0, isOriginCCW = 0;
3657 if ( GEOSCoordSeq_isCCW_r( context, GEOSGeom_getCoordSeq_r( context, result.get() ), &isResultCCW ) &&
3658 GEOSCoordSeq_isCCW_r( context, GEOSGeom_getCoordSeq_r( context, line ), &isOriginCCW )
3662 reverseLine = ( isOriginCCW == 1 && isResultCCW == 0 ) || ( isOriginCCW == 0 && isResultCCW == 1 );
3668 double x1res, y1res, x2res, y2res;
3669 if ( !_linestringEndpoints( result.get(), x1res, y1res, x2res, y2res ) )
3671 geos::unique_ptr beginResultLineVertex = createGeosPointXY( x1res, y1res,
false, 0,
false, 0, 2,
precision );
3672 geos::unique_ptr endResultLineVertex = createGeosPointXY( x2res, y2res,
false, 0,
false, 0, 2,
precision );
3673 reverseLine = GEOSEquals_r( context, beginLineVertex.get(), endResultLineVertex.get() ) == 1 || GEOSEquals_r( context, endLineVertex.get(), beginResultLineVertex.get() ) == 1;
3676 result.reset( GEOSReverse_r( context, result.get() ) );
3684 int nIntersections = 0;
3685 int lastIntersectingRing = -2;
3689 int nRings = GEOSGetNumInteriorRings_r( context, polygon );
3694 const GEOSGeometry *outerRing = GEOSGetExteriorRing_r( context, polygon );
3695 if ( GEOSIntersects_r( context, outerRing, reshapeLineGeos ) == 1 )
3698 lastIntersectingRing = -1;
3699 lastIntersectingGeom = outerRing;
3707 for (
int i = 0; i < nRings; ++i )
3709 innerRings[i] = GEOSGetInteriorRingN_r( context, polygon, i );
3710 if ( GEOSIntersects_r( context, innerRings[i], reshapeLineGeos ) == 1 )
3713 lastIntersectingRing = i;
3714 lastIntersectingGeom = innerRings[i];
3718 catch ( GEOSException & )
3723 if ( nIntersections != 1 )
3725 delete [] innerRings;
3730 geos::unique_ptr reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos,
precision );
3731 if ( !reshapeResult )
3733 delete [] innerRings;
3739 const GEOSCoordSequence *reshapeSequence = GEOSGeom_getCoordSeq_r( context, reshapeResult.get() );
3740 GEOSCoordSequence *newCoordSequence = GEOSCoordSeq_clone_r( context, reshapeSequence );
3742 reshapeResult.reset();
3744 newRing = GEOSGeom_createLinearRing_r( context, newCoordSequence );
3747 delete [] innerRings;
3752 if ( lastIntersectingRing == -1 )
3753 newOuterRing = newRing;
3755 newOuterRing = GEOSGeom_clone_r( context, outerRing );
3758 QVector<GEOSGeometry *> ringList;
3761 GEOSGeometry *outerRingPoly = GEOSGeom_createPolygon_r( context, GEOSGeom_clone_r( context, newOuterRing ),
nullptr, 0 );
3762 if ( outerRingPoly )
3764 ringList.reserve( nRings );
3766 for (
int i = 0; i < nRings; ++i )
3768 if ( lastIntersectingRing == i )
3769 currentRing = newRing;
3771 currentRing = GEOSGeom_clone_r( context, innerRings[i] );
3774 if ( GEOSContains_r( context, outerRingPoly, currentRing ) == 1 )
3775 ringList.push_back( currentRing );
3777 GEOSGeom_destroy_r( context, currentRing );
3780 GEOSGeom_destroy_r( context, outerRingPoly );
3784 for (
int i = 0; i < ringList.size(); ++i )
3785 newInnerRings[i] = ringList.at( i );
3787 delete [] innerRings;
3789 geos::unique_ptr reshapedPolygon( GEOSGeom_createPolygon_r( context, newOuterRing, newInnerRings, ringList.size() ) );
3790 delete[] newInnerRings;
3792 return reshapedPolygon;
3797 if ( !line1 || !line2 )
3802 double bufferDistance = std::pow( 10.0L, geomDigits( line2 ) - 11 );
3809 geos::unique_ptr intersectionGeom( GEOSIntersection_r( context, bufferGeom.get(), line1 ) );
3812 double intersectGeomLength;
3815 GEOSLength_r( context, intersectionGeom.get(), &intersectGeomLength );
3816 GEOSLength_r( context, line1, &line1Length );
3818 double intersectRatio = line1Length / intersectGeomLength;
3819 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
3827 if ( !point || !line )
3830 double bufferDistance = std::pow( 10.0L, geomDigits( line ) - 11 );
3833 geos::unique_ptr lineBuffer( GEOSBuffer_r( context, line, bufferDistance, 8 ) );
3837 bool contained =
false;
3838 if ( GEOSContains_r( context, lineBuffer.get(), point ) == 1 )
3847 geos::unique_ptr bbox( GEOSEnvelope_r( context, geom ) );
3851 const GEOSGeometry *bBoxRing = GEOSGetExteriorRing_r( context, bbox.get() );
3855 const GEOSCoordSequence *bBoxCoordSeq = GEOSGeom_getCoordSeq_r( context, bBoxRing );
3857 if ( !bBoxCoordSeq )
3860 unsigned int nCoords = 0;
3861 if ( !GEOSCoordSeq_getSize_r( context, bBoxCoordSeq, &nCoords ) )
3865 for (
unsigned int i = 0; i < nCoords - 1; ++i )
3868 GEOSCoordSeq_getX_r( context, bBoxCoordSeq, i, &t );
3871 digits = std::ceil( std::log10( std::fabs( t ) ) );
3872 if ( digits > maxDigits )
3875 GEOSCoordSeq_getY_r( context, bBoxCoordSeq, i, &t );
3876 digits = std::ceil( std::log10( std::fabs( t ) ) );
3877 if ( digits > maxDigits )
The Qgis class provides global constants for use throughout the application.
BufferSide
Side of line to buffer.
@ Right
Buffer to right of line.
GeometryOperationResult
Success or failure of a geometry operation.
@ AddPartNotMultiGeometry
The source geometry is not multi.
@ InvalidBaseGeometry
The base geometry on which the operation is done is invalid or empty.
@ RejectOnInvalidSubGeometry
Don't allow geometries with invalid sub-geometries to be created.
@ SkipEmptyInteriorRings
Skip any empty polygon interior ring.
QFlags< GeosCreationFlag > GeosCreationFlags
Geos geometry creation behavior flags.
JoinStyle
Join styles for buffers.
EndCapStyle
End cap styles for buffers.
CoverageValidityResult
Coverage validity results.
@ Valid
Coverage is valid.
@ Invalid
Coverage is invalid. Invalidity includes polygons that overlap, that have gaps smaller than the gap w...
@ Error
An exception occurred while determining validity.
MakeValidMethod
Algorithms to use when repairing invalid geometries.
@ Linework
Combines all rings into a set of noded lines and then extracts valid polygons from that linework.
@ Structure
Structured method, first makes all rings valid and then merges shells and subtracts holes from shells...
WkbType
The WKB type describes the number of dimensions a geometry has.
@ LineStringZM
LineStringZM.
@ GeometryCollection
GeometryCollection.
@ LineStringZ
LineStringZ.
@ PolyhedralSurface
PolyhedralSurface.
Abstract base class for all geometries.
virtual const QgsAbstractGeometry * simplifiedTypeRef() const
Returns a reference to the simplest lossless representation of this geometry, e.g.
bool isMeasure() const
Returns true if the geometry contains m values.
virtual QgsRectangle boundingBox() const
Returns the minimal bounding box for the geometry.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual QgsPoint vertexAt(QgsVertexId id) const =0
Returns the point corresponding to a specified vertex id.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
virtual bool isEmpty() const
Returns true if the geometry is empty.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
Curve polygon geometry type.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
Abstract base class for curved geometry type.
virtual QgsLineString * curveToLine(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve.
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
static Qgis::GeometryOperationResult addPart(QgsAbstractGeometry *geometry, std::unique_ptr< QgsAbstractGeometry > part)
Add a part to multi type geometry.
A geometry engine is a low-level representation of a QgsAbstractGeometry object, optimised for use wi...
const QgsAbstractGeometry * mGeometry
EngineOperationResult
Success or failure of a geometry operation.
@ NothingHappened
Nothing happened, without any error.
@ InvalidBaseGeometry
The geometry on which the operation occurs is not valid.
@ InvalidInput
The input is not valid.
@ NodedGeometryError
Error occurred while creating a noded geometry.
@ EngineError
Error occurred in the geometry engine.
@ SplitCannotSplitPoint
Points cannot be split.
@ Success
Operation succeeded.
void logError(const QString &engineName, const QString &message) const
Logs an error message encountered during an operation.
static std::unique_ptr< QgsGeometryCollection > createCollectionOfType(Qgis::WkbType type)
Returns a new geometry collection matching a specified WKB type.
Encapsulates parameters under which a geometry operation is performed.
double gridSize() const
Returns the grid size which will be used to snap vertices of a geometry.
static double sqrDistance2D(double x1, double y1, double x2, double y2)
Returns the squared 2D distance between (x1, y1) and (x2, y2).
A geometry is the spatial representation of a feature.
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Used to create and store a proj context object, correctly freeing the context upon destruction.
static GEOSContextHandle_t get()
Returns a thread local instance of a GEOS context, safe for use in the current thread.
Does vector analysis using the GEOS library and handles import, export, and exception handling.
double frechetDistanceDensify(const QgsAbstractGeometry *geometry, double densifyFraction, QString *errorMsg=nullptr) const
Returns the Fréchet distance between this geometry and another geometry, restricted to discrete point...
double hausdorffDistanceDensify(const QgsAbstractGeometry *geometry, double densifyFraction, QString *errorMsg=nullptr) const
Returns the Hausdorff distance between this geometry and another geometry.
bool touches(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom touches this.
static geos::unique_ptr offsetCurve(const GEOSGeometry *geometry, double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg=nullptr)
Directly calculates the offset curve for a GEOS geometry object and returns a GEOS geometry result.
QgsAbstractGeometry * symDifference(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters ¶meters=QgsGeometryParameters()) const override
Calculate the symmetric difference of this and geom.
bool isValid(QString *errorMsg=nullptr, bool allowSelfTouchingHoles=false, QgsGeometry *errorLoc=nullptr) const override
Returns true if the geometry is valid.
std::unique_ptr< QgsAbstractGeometry > subdivide(int maxNodes, QString *errorMsg=nullptr, const QgsGeometryParameters ¶meters=QgsGeometryParameters()) const
Subdivides the geometry.
std::unique_ptr< QgsAbstractGeometry > mergeLines(QString *errorMsg=nullptr) const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
QgsAbstractGeometry * intersection(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters ¶meters=QgsGeometryParameters()) const override
Calculate the intersection of this and geom.
std::unique_ptr< QgsAbstractGeometry > constrainedDelaunayTriangulation(QString *errorMsg=nullptr) const
Returns a constrained Delaunay triangulation for the vertices of the geometry.
QgsAbstractGeometry * buffer(double distance, int segments, QString *errorMsg=nullptr) const override
std::unique_ptr< QgsAbstractGeometry > closestPoint(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the closest point on the geometry to the other geometry.
std::unique_ptr< QgsAbstractGeometry > reshapeGeometry(const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg=nullptr) const
Reshapes the geometry using a line.
std::unique_ptr< QgsAbstractGeometry > maximumInscribedCircle(double tolerance, QString *errorMsg=nullptr) const
Returns the maximum inscribed circle.
static geos::unique_ptr asGeos(const QgsGeometry &geometry, double precision=0, Qgis::GeosCreationFlags flags=Qgis::GeosCreationFlags())
Returns a geos geometry - caller takes ownership of the object (should be deleted with GEOSGeom_destr...
QgsAbstractGeometry * difference(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters ¶meters=QgsGeometryParameters()) const override
Calculate the difference of this and geom.
EngineOperationResult splitGeometry(const QgsLineString &splitLine, QVector< QgsGeometry > &newGeometries, bool topological, QgsPointSequence &topologyTestPoints, QString *errorMsg=nullptr, bool skipIntersectionCheck=false) const override
Splits this geometry according to a given line.
std::unique_ptr< QgsAbstractGeometry > sharedPaths(const QgsAbstractGeometry *other, QString *errorMsg=nullptr) const
Find paths shared between the two given lineal geometries (this and other).
std::unique_ptr< QgsAbstractGeometry > clip(const QgsRectangle &rectangle, QString *errorMsg=nullptr) const
Performs a fast, non-robust intersection between the geometry and a rectangle.
std::unique_ptr< QgsAbstractGeometry > node(QString *errorMsg=nullptr) const
Returns a (Multi)LineString representing the fully noded version of a collection of linestrings.
double minimumClearance(QString *errorMsg=nullptr) const
Computes the minimum clearance of a geometry.
std::unique_ptr< QgsAbstractGeometry > singleSidedBuffer(double distance, int segments, Qgis::BufferSide side, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg=nullptr) const
Returns a single sided buffer for a geometry.
bool disjoint(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom is disjoint from this.
std::unique_ptr< QgsAbstractGeometry > makeValid(Qgis::MakeValidMethod method=Qgis::MakeValidMethod::Linework, bool keepCollapsed=false, QString *errorMsg=nullptr) const
Repairs the geometry using GEOS make valid routine.
std::unique_ptr< QgsAbstractGeometry > voronoiDiagram(const QgsAbstractGeometry *extent=nullptr, double tolerance=0.0, bool edgesOnly=false, QString *errorMsg=nullptr) const
Creates a Voronoi diagram for the nodes contained within the geometry.
QgsAbstractGeometry * simplify(double tolerance, QString *errorMsg=nullptr) const override
std::unique_ptr< QgsAbstractGeometry > unionCoverage(QString *errorMsg=nullptr) const
Optimized union algorithm for polygonal inputs that are correctly noded and do not overlap.
static std::unique_ptr< QgsPolygon > fromGeosPolygon(const GEOSGeometry *geos)
double hausdorffDistance(const QgsAbstractGeometry *geometry, QString *errorMsg=nullptr) const
Returns the Hausdorff distance between this geometry and another geometry.
std::unique_ptr< QgsAbstractGeometry > minimumClearanceLine(QString *errorMsg=nullptr) const
Returns a LineString whose endpoints define the minimum clearance of a geometry.
QgsAbstractGeometry * envelope(QString *errorMsg=nullptr) const override
QString relate(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Returns the Dimensional Extended 9 Intersection Model (DE-9IM) representation of the relationship bet...
bool crosses(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom crosses this.
QgsAbstractGeometry * convexHull(QString *errorMsg=nullptr) const override
Calculate the convex hull of this.
double distance(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculates the distance between this and geom.
std::unique_ptr< QgsAbstractGeometry > minimumWidth(QString *errorMsg=nullptr) const
Returns a linestring geometry which represents the minimum diameter of the geometry.
bool isSimple(QString *errorMsg=nullptr) const override
Determines whether the geometry is simple (according to OGC definition).
QgsGeos(const QgsAbstractGeometry *geometry, double precision=0, Qgis::GeosCreationFlags flags=Qgis::GeosCreationFlag::SkipEmptyInteriorRings)
GEOS geometry engine constructor.
QgsAbstractGeometry * combine(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters ¶meters=QgsGeometryParameters()) const override
Calculate the combination of this and geom.
bool within(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom is within this.
std::unique_ptr< QgsAbstractGeometry > shortestLine(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the shortest line joining this geometry to the other geometry.
std::unique_ptr< QgsAbstractGeometry > concaveHull(double targetPercent, bool allowHoles=false, QString *errorMsg=nullptr) const
Returns a possibly concave geometry that encloses the input geometry.
void prepareGeometry() override
Prepares the geometry, so that subsequent calls to spatial relation methods are much faster.
static std::unique_ptr< QgsAbstractGeometry > fromGeos(const GEOSGeometry *geos)
Create a geometry from a GEOSGeometry.
QgsAbstractGeometry * interpolate(double distance, QString *errorMsg=nullptr) const override
bool distanceWithin(const QgsAbstractGeometry *geom, double maxdistance, QString *errorMsg=nullptr) const override
Checks if geom is within maxdistance distance from this geometry.
bool isEqual(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if this is equal to geom.
bool contains(double x, double y, QString *errorMsg=nullptr) const
Returns true if the geometry contains the point at (x, y).
static QgsGeometry polygonize(const QVector< const QgsAbstractGeometry * > &geometries, QString *errorMsg=nullptr)
Creates a GeometryCollection geometry containing possible polygons formed from the constituent linewo...
double frechetDistance(const QgsAbstractGeometry *geometry, QString *errorMsg=nullptr) const
Returns the Fréchet distance between this geometry and another geometry, restricted to discrete point...
bool relatePattern(const QgsAbstractGeometry *geom, const QString &pattern, QString *errorMsg=nullptr) const override
Tests whether two geometries are related by a specified Dimensional Extended 9 Intersection Model (DE...
QgsPoint * centroid(QString *errorMsg=nullptr) const override
Calculates the centroid of this.
double lineLocatePoint(const QgsPoint &point, QString *errorMsg=nullptr) const
Returns a distance representing the location along this linestring of the closest point on this lines...
std::unique_ptr< QgsAbstractGeometry > simplifyCoverageVW(double tolerance, bool preserveBoundary, QString *errorMsg=nullptr) const
Operates on a coverage (represented as a list of polygonal geometry with exactly matching edge geomet...
bool isEmpty(QString *errorMsg=nullptr) const override
static Qgis::GeometryOperationResult addPart(QgsGeometry &geometry, GEOSGeometry *newPart)
Adds a new island polygon to a multipolygon feature.
bool intersects(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom intersects this.
void geometryChanged() override
Should be called whenever the geometry associated with the engine has been modified and the engine mu...
bool overlaps(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom overlaps this.
double area(QString *errorMsg=nullptr) const override
std::unique_ptr< QgsAbstractGeometry > delaunayTriangulation(double tolerance=0.0, bool edgesOnly=false, QString *errorMsg=nullptr) const
Returns the Delaunay triangulation for the vertices of the geometry.
double length(QString *errorMsg=nullptr) const override
std::unique_ptr< QgsAbstractGeometry > largestEmptyCircle(double tolerance, const QgsAbstractGeometry *boundary=nullptr, QString *errorMsg=nullptr) const
Constructs the Largest Empty Circle for a set of obstacle geometries, up to a specified tolerance.
Qgis::CoverageValidityResult validateCoverage(double gapWidth, std::unique_ptr< QgsAbstractGeometry > *invalidEdges, QString *errorMsg=nullptr) const
Analyze a coverage (represented as a collection of polygonal geometry with exactly matching edge geom...
static QgsPoint coordSeqPoint(const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM)
QgsPoint * pointOnSurface(QString *errorMsg=nullptr) const override
Calculate a point that is guaranteed to be on the surface of this.
static QgsGeometry geometryFromGeos(GEOSGeometry *geos)
Creates a new QgsGeometry object, feeding in a geometry in GEOS format.
Line string geometry type, with support for z-dimension and m-values.
const double * yData() const
Returns a const pointer to the y vertex data.
const double * xData() const
Returns a const pointer to the x vertex data.
QVector< double > xVector() const
Returns the x vertex values as a vector.
const double * zData() const
Returns a const pointer to the z vertex data, or nullptr if the linestring does not have z values.
int numPoints() const override
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
QVector< double > yVector() const
Returns the y vertex values as a vector.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
QVector< double > zVector() const
Returns the z vertex values as a vector.
double closestSegment(const QgsPoint &pt, QgsPoint &segmentPt, QgsVertexId &vertexAfter, int *leftOf=nullptr, double epsilon=4 *std::numeric_limits< double >::epsilon()) const override
Searches for the closest segment of the geometry to a given point.
const double * mData() const
Returns a const pointer to the m vertex data, or nullptr if the linestring does not have m values.
void addVertex(const QgsPoint &pt)
Adds a new vertex to the end of the line string.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
Multi curve geometry collection.
bool addGeometry(QgsAbstractGeometry *g) override
Adds a geometry and takes ownership. Returns true in case of success.
Custom exception class which is raised when an operation is not supported.
A class to represent a 2D point.
Point geometry type, with support for z-dimension and m-values.
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Polyhedral surface geometry type.
int numPatches() const
Returns the number of patches contained with the polyhedral surface.
const QgsPolygon * patchN(int i) const
Retrieves a patch from the polyhedral surface.
A rectangle specified with double values.
void setYMinimum(double y)
Set the minimum y value.
void setXMinimum(double x)
Set the minimum x value.
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
static Qgis::GeometryType geometryType(Qgis::WkbType type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
static bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
static Qgis::WkbType flatType(Qgis::WkbType type)
Returns the flat type for a WKB type.
Contains geos related utilities and functions.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
QVector< QgsPoint > QgsPointSequence
#define DEFAULT_QUADRANT_SEGMENTS
#define CATCH_GEOS_WITH_ERRMSG(r)
#define QgsDebugError(str)
QLineF segment(int index, QRectF rect, double radius)
Utility class for identifying a unique vertex within a geometry.
struct GEOSGeom_t GEOSGeometry