31 #define DEFAULT_QUADRANT_SEGMENTS 8 33 #define CATCH_GEOS(r) \ 34 catch (GEOSException &e) \ 36 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr("GEOS") ); \ 40 #define CATCH_GEOS_WITH_ERRMSG(r) \ 41 catch (GEOSException &e) \ 43 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr("GEOS") ); \ 46 *errorMsg = e.what(); \ 53 static void throwGEOSException(
const char *fmt, ... )
59 vsnprintf( buffer,
sizeof buffer, fmt, ap );
62 qWarning(
"GEOS exception: %s", buffer );
64 throw GEOSException( QString::fromUtf8( buffer ) );
67 static void printGEOSNotice(
const char *fmt, ... )
69 #if defined(QGISDEBUG) 74 vsnprintf( buffer,
sizeof buffer, fmt, ap );
77 QgsDebugMsg( QString(
"GEOS notice: %1" ).arg( QString::fromUtf8( buffer ) ) );
86 GEOSContextHandle_t ctxt;
90 ctxt = initGEOS_r( printGEOSNotice, throwGEOSException );
98 GEOSInit(
const GEOSInit &rh ) =
delete;
99 GEOSInit &operator=(
const GEOSInit &rh ) =
delete;
102 static GEOSInit geosinit;
106 GEOSGeom_destroy_r( geosinit.ctxt, geom );
111 GEOSPreparedGeom_destroy_r( geosinit.ctxt, geom );
116 GEOSBufferParams_destroy_r( geosinit.ctxt, params );
121 GEOSCoordSeq_destroy_r( geosinit.ctxt, sequence );
131 , mPrecision( precision )
139 mGeosPrepared.reset();
145 mGeosPrepared.reset();
148 mGeosPrepared.reset( GEOSPrepare_r( geosinit.ctxt, mGeos.get() ) );
152 void QgsGeos::cacheGeos()
const 164 return overlay( geom, OverlayIntersection, errorMsg ).release();
169 return overlay( geom, OverlayDifference, errorMsg ).release();
184 catch ( GEOSException &e )
188 *errorMsg = e.what();
199 int partType = GEOSGeomTypeId_r( geosinit.ctxt, currentPart );
202 if ( partType == GEOS_POINT )
213 if ( partType == GEOS_MULTILINESTRING || partType == GEOS_MULTIPOLYGON || partType == GEOS_GEOMETRYCOLLECTION )
215 int partCount = GEOSGetNumGeometries_r( geosinit.ctxt, currentPart );
216 for (
int i = 0; i < partCount; ++i )
218 subdivideRecursive( GEOSGetGeometryN_r( geosinit.ctxt, currentPart, i ), maxNodes, depth, parts, clipRect );
229 int vertexCount = GEOSGetNumCoordinates_r( geosinit.ctxt, currentPart );
230 if ( vertexCount == 0 )
234 else if ( vertexCount < maxNodes )
241 double width = clipRect.
width();
242 double height = clipRect.
height();
245 if ( width > height )
278 subdivideRecursive( clipPart1.get(), maxNodes, depth, parts, halfClipRect1 );
282 subdivideRecursive( clipPart2.get(), maxNodes, depth, parts, halfClipRect2 );
294 maxNodes = std::max( maxNodes, 8 );
303 return std::move( parts );
308 return overlay( geom, OverlayUnion, errorMsg ).release();
313 QVector< GEOSGeometry * > geosGeometries;
314 geosGeometries.reserve( geomList.size() );
320 geosGeometries <<
asGeos( g, mPrecision ).release();
326 geos::unique_ptr geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geosGeometries );
327 geomUnion.reset( GEOSUnaryUnion_r( geosinit.ctxt, geomCollection.get() ) );
331 std::unique_ptr< QgsAbstractGeometry > result =
fromGeos( geomUnion.get() );
332 return result.release();
337 QVector< GEOSGeometry * > geosGeometries;
338 geosGeometries.reserve( geomList.size() );
344 geosGeometries <<
asGeos( g.constGet(), mPrecision ).release();
350 geos::unique_ptr geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geosGeometries );
351 geomUnion.reset( GEOSUnaryUnion_r( geosinit.ctxt, geomCollection.get() ) );
355 std::unique_ptr< QgsAbstractGeometry > result =
fromGeos( geomUnion.get() );
356 return result.release();
361 return overlay( geom, OverlaySymDifference, errorMsg ).release();
373 if ( !otherGeosGeom )
380 GEOSDistance_r( geosinit.ctxt, mGeos.get(), otherGeosGeom.get(), &
distance );
396 if ( !otherGeosGeom )
403 GEOSHausdorffDistance_r( geosinit.ctxt, mGeos.get(), otherGeosGeom.get(), &
distance );
419 if ( !otherGeosGeom )
426 GEOSHausdorffDistanceDensify_r( geosinit.ctxt, mGeos.get(), otherGeosGeom.get(), densifyFraction, &
distance );
435 return relation( geom, RelationIntersects, errorMsg );
440 return relation( geom, RelationTouches, errorMsg );
445 return relation( geom, RelationCrosses, errorMsg );
450 return relation( geom, RelationWithin, errorMsg );
455 return relation( geom, RelationOverlaps, errorMsg );
460 return relation( geom, RelationContains, errorMsg );
465 return relation( geom, RelationDisjoint, errorMsg );
484 char *r = GEOSRelate_r( geosinit.ctxt, mGeos.get(), geosGeom.get() );
487 result = QString( r );
488 GEOSFree_r( geosinit.ctxt, r );
491 catch ( GEOSException &e )
495 *errorMsg = e.what();
504 if ( !mGeos || !geom )
518 result = ( GEOSRelatePattern_r( geosinit.ctxt, mGeos.get(), geosGeom.get(), pattern.toLocal8Bit().constData() ) == 1 );
520 catch ( GEOSException &e )
524 *errorMsg = e.what();
541 if ( GEOSArea_r( geosinit.ctxt, mGeos.get(), &
area ) != 1 )
557 if ( GEOSLength_r( geosinit.ctxt, mGeos.get(), &
length ) != 1 )
565 QVector<QgsGeometry> &newGeometries,
568 QString *errorMsg )
const 583 if ( !GEOSisValid_r( geosinit.ctxt, mGeos.get() ) )
591 newGeometries.clear();
598 splitLineGeos = createGeosLinestring( &splitLine, mPrecision );
602 splitLineGeos = createGeosPointXY( splitLine.
xAt( 0 ), splitLine.
yAt( 0 ),
false, 0,
false, 0, 2, mPrecision );
609 if ( !GEOSisValid_r( geosinit.ctxt, splitLineGeos.get() ) || !GEOSisSimple_r( geosinit.ctxt, splitLineGeos.get() ) )
617 if ( !topologicalTestPointsSplit( splitLineGeos.get(), topologyTestPoints ) )
626 returnCode = splitLinearGeometry( splitLineGeos.get(), newGeometries );
630 returnCode = splitPolygonGeometry( splitLineGeos.get(), newGeometries );
644 bool QgsGeos::topologicalTestPointsSplit(
const GEOSGeometry *splitLine,
QgsPointSequence &testPoints, QString *errorMsg )
const 658 geos::unique_ptr intersectionGeom( GEOSIntersection_r( geosinit.ctxt, mGeos.get(), splitLine ) );
659 if ( !intersectionGeom )
663 int nIntersectGeoms = 1;
664 if ( GEOSGeomTypeId_r( geosinit.ctxt, intersectionGeom.get() ) == GEOS_LINESTRING
665 || GEOSGeomTypeId_r( geosinit.ctxt, intersectionGeom.get() ) == GEOS_POINT )
669 nIntersectGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, intersectionGeom.get() );
671 for (
int i = 0; i < nIntersectGeoms; ++i )
673 const GEOSGeometry *currentIntersectGeom =
nullptr;
675 currentIntersectGeom = intersectionGeom.get();
677 currentIntersectGeom = GEOSGetGeometryN_r( geosinit.ctxt, intersectionGeom.get(), i );
679 const GEOSCoordSequence *lineSequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, currentIntersectGeom );
680 unsigned int sequenceSize = 0;
682 if ( GEOSCoordSeq_getSize_r( geosinit.ctxt, lineSequence, &sequenceSize ) != 0 )
684 for (
unsigned int i = 0; i < sequenceSize; ++i )
686 if ( GEOSCoordSeq_getX_r( geosinit.ctxt, lineSequence, i, &x ) != 0 )
688 if ( GEOSCoordSeq_getY_r( geosinit.ctxt, lineSequence, i, &y ) != 0 )
690 testPoints.push_back(
QgsPoint( x, y ) );
702 geos::unique_ptr QgsGeos::linePointDifference( GEOSGeometry *GEOSsplitPoint )
const 704 int type = GEOSGeomTypeId_r( geosinit.ctxt, mGeos.get() );
706 std::unique_ptr< QgsMultiCurve > multiCurve;
707 if ( type == GEOS_MULTILINESTRING )
709 multiCurve.reset( qgsgeometry_cast<QgsMultiCurve *>(
mGeometry->
clone() ) );
711 else if ( type == GEOS_LINESTRING )
727 std::unique_ptr< QgsAbstractGeometry > splitGeom(
fromGeos( GEOSsplitPoint ) );
737 for (
int i = 0; i < multiCurve->numGeometries(); ++i )
746 for (
int j = 1; j < ( nVertices - 1 ); ++j )
750 if ( currentPoint == *splitPoint )
762 return asGeos( &lines, mPrecision );
774 if ( !GEOSIntersects_r( geosinit.ctxt, splitLine, mGeos.get() ) )
778 int linearIntersect = GEOSRelatePattern_r( geosinit.ctxt, mGeos.get(), splitLine,
"1********" );
779 if ( linearIntersect > 0 )
782 int splitGeomType = GEOSGeomTypeId_r( geosinit.ctxt, splitLine );
785 if ( splitGeomType == GEOS_POINT )
787 splitGeom = linePointDifference( splitLine );
791 splitGeom.reset( GEOSDifference_r( geosinit.ctxt, mGeos.get(), splitLine ) );
793 QVector<GEOSGeometry *> lineGeoms;
795 int splitType = GEOSGeomTypeId_r( geosinit.ctxt, splitGeom.get() );
796 if ( splitType == GEOS_MULTILINESTRING )
798 int nGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, splitGeom.get() );
799 lineGeoms.reserve( nGeoms );
800 for (
int i = 0; i < nGeoms; ++i )
801 lineGeoms << GEOSGeom_clone_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, splitGeom.get(), i ) );
806 lineGeoms << GEOSGeom_clone_r( geosinit.ctxt, splitGeom.get() );
809 mergeGeometriesMultiTypeSplit( lineGeoms );
811 for (
int i = 0; i < lineGeoms.size(); ++i )
814 GEOSGeom_destroy_r( geosinit.ctxt, lineGeoms[i] );
829 if ( !GEOSIntersects_r( geosinit.ctxt, splitLine, mGeos.get() ) )
834 if ( !nodedGeometry )
837 const GEOSGeometry *noded = nodedGeometry.get();
839 if ( !polygons || numberOfGeometries( polygons.get() ) == 0 )
846 QVector<GEOSGeometry *> testedGeometries;
852 for (
int i = 0; i < numberOfGeometries( polygons.get() ); i++ )
854 const GEOSGeometry *polygon = GEOSGetGeometryN_r( geosinit.ctxt, polygons.get(), i );
855 intersectGeometry.reset( GEOSIntersection_r( geosinit.ctxt, mGeos.get(), polygon ) );
856 if ( !intersectGeometry )
862 double intersectionArea;
863 GEOSArea_r( geosinit.ctxt, intersectGeometry.get(), &intersectionArea );
866 GEOSArea_r( geosinit.ctxt, polygon, &polygonArea );
868 const double areaRatio = intersectionArea / polygonArea;
869 if ( areaRatio > 0.99 && areaRatio < 1.01 )
870 testedGeometries << GEOSGeom_clone_r( geosinit.ctxt, polygon );
873 int nGeometriesThis = numberOfGeometries( mGeos.get() );
874 if ( testedGeometries.empty() || testedGeometries.size() == nGeometriesThis )
877 for (
int i = 0; i < testedGeometries.size(); ++i )
879 GEOSGeom_destroy_r( geosinit.ctxt, testedGeometries[i] );
884 mergeGeometriesMultiTypeSplit( testedGeometries );
887 for ( i = 0; i < testedGeometries.size() && GEOSisValid_r( geosinit.ctxt, testedGeometries[i] ); ++i )
890 if ( i < testedGeometries.size() )
892 for ( i = 0; i < testedGeometries.size(); ++i )
893 GEOSGeom_destroy_r( geosinit.ctxt, testedGeometries[i] );
898 for ( i = 0; i < testedGeometries.size(); ++i )
901 GEOSGeom_destroy_r( geosinit.ctxt, testedGeometries[i] );
907 geos::unique_ptr QgsGeos::nodeGeometries(
const GEOSGeometry *splitLine,
const GEOSGeometry *geom )
909 if ( !splitLine || !geom )
913 if ( GEOSGeomTypeId_r( geosinit.ctxt, geom ) == GEOS_POLYGON || GEOSGeomTypeId_r( geosinit.ctxt, geom ) == GEOS_MULTIPOLYGON )
914 geometryBoundary.reset( GEOSBoundary_r( geosinit.ctxt, geom ) );
916 geometryBoundary.reset( GEOSGeom_clone_r( geosinit.ctxt, geom ) );
918 geos::unique_ptr splitLineClone( GEOSGeom_clone_r( geosinit.ctxt, splitLine ) );
919 geos::unique_ptr unionGeometry( GEOSUnion_r( geosinit.ctxt, splitLineClone.get(), geometryBoundary.get() ) );
921 return unionGeometry;
924 int QgsGeos::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry *> &splitResult )
const 930 int type = GEOSGeomTypeId_r( geosinit.ctxt, mGeos.get() );
931 if ( type != GEOS_GEOMETRYCOLLECTION &&
932 type != GEOS_MULTILINESTRING &&
933 type != GEOS_MULTIPOLYGON &&
934 type != GEOS_MULTIPOINT )
937 QVector<GEOSGeometry *> copyList = splitResult;
941 QVector<GEOSGeometry *> unionGeom;
943 for (
int i = 0; i < copyList.size(); ++i )
947 for (
int j = 0; j < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos.get() ); j++ )
949 if ( GEOSEquals_r( geosinit.ctxt, copyList[i], GEOSGetGeometryN_r( geosinit.ctxt, mGeos.get(), j ) ) )
958 unionGeom << copyList[i];
962 QVector<GEOSGeometry *> geomVector;
963 geomVector << copyList[i];
965 if ( type == GEOS_MULTILINESTRING )
966 splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector ).release();
967 else if ( type == GEOS_MULTIPOLYGON )
968 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector ).release();
970 GEOSGeom_destroy_r( geosinit.ctxt, copyList[i] );
975 if ( !unionGeom.isEmpty() )
977 if ( type == GEOS_MULTILINESTRING )
978 splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom ).release();
979 else if ( type == GEOS_MULTIPOLYGON )
980 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom ).release();
990 geos::unique_ptr QgsGeos::createGeosCollection(
int typeId,
const QVector<GEOSGeometry *> &geoms )
992 int nNullGeoms = geoms.count(
nullptr );
993 int nNotNullGeoms = geoms.size() - nNullGeoms;
995 GEOSGeometry **geomarr =
new GEOSGeometry*[ nNotNullGeoms ];
1002 QVector<GEOSGeometry *>::const_iterator geomIt = geoms.constBegin();
1003 for ( ; geomIt != geoms.constEnd(); ++geomIt )
1007 geomarr[i] = *geomIt;
1015 geom.reset( GEOSGeom_createCollection_r( geosinit.ctxt, typeId, geomarr, nNotNullGeoms ) );
1017 catch ( GEOSException &e )
1034 int nCoordDims = GEOSGeom_getCoordinateDimension_r( geosinit.ctxt, geos );
1035 int nDims = GEOSGeom_getDimensions_r( geosinit.ctxt, geos );
1036 bool hasZ = ( nCoordDims == 3 );
1037 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1039 switch ( GEOSGeomTypeId_r( geosinit.ctxt, geos ) )
1043 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, geos );
1044 return std::unique_ptr<QgsAbstractGeometry>(
coordSeqPoint( cs, 0, hasZ, hasM ).
clone() );
1046 case GEOS_LINESTRING:
1048 return sequenceToLinestring( geos, hasZ, hasM );
1054 case GEOS_MULTIPOINT:
1056 std::unique_ptr< QgsMultiPoint > multiPoint(
new QgsMultiPoint() );
1057 int nParts = GEOSGetNumGeometries_r( geosinit.ctxt, geos );
1058 for (
int i = 0; i < nParts; ++i )
1060 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, geos, i ) );
1063 multiPoint->addGeometry(
coordSeqPoint( cs, 0, hasZ, hasM ).clone() );
1066 return std::move( multiPoint );
1068 case GEOS_MULTILINESTRING:
1071 int nParts = GEOSGetNumGeometries_r( geosinit.ctxt, geos );
1072 for (
int i = 0; i < nParts; ++i )
1074 std::unique_ptr< QgsLineString >line( sequenceToLinestring( GEOSGetGeometryN_r( geosinit.ctxt, geos, i ), hasZ, hasM ) );
1077 multiLineString->addGeometry( line.release() );
1080 return std::move( multiLineString );
1082 case GEOS_MULTIPOLYGON:
1084 std::unique_ptr< QgsMultiPolygon > multiPolygon(
new QgsMultiPolygon() );
1086 int nParts = GEOSGetNumGeometries_r( geosinit.ctxt, geos );
1087 for (
int i = 0; i < nParts; ++i )
1089 std::unique_ptr< QgsPolygon > poly =
fromGeosPolygon( GEOSGetGeometryN_r( geosinit.ctxt, geos, i ) );
1092 multiPolygon->addGeometry( poly.release() );
1095 return std::move( multiPolygon );
1097 case GEOS_GEOMETRYCOLLECTION:
1100 int nParts = GEOSGetNumGeometries_r( geosinit.ctxt, geos );
1101 for (
int i = 0; i < nParts; ++i )
1103 std::unique_ptr< QgsAbstractGeometry > geom(
fromGeos( GEOSGetGeometryN_r( geosinit.ctxt, geos, i ) ) );
1106 geomCollection->addGeometry( geom.release() );
1109 return std::move( geomCollection );
1117 if ( GEOSGeomTypeId_r( geosinit.ctxt, geos ) != GEOS_POLYGON )
1122 int nCoordDims = GEOSGeom_getCoordinateDimension_r( geosinit.ctxt, geos );
1123 int nDims = GEOSGeom_getDimensions_r( geosinit.ctxt, geos );
1124 bool hasZ = ( nCoordDims == 3 );
1125 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1127 std::unique_ptr< QgsPolygon > polygon(
new QgsPolygon() );
1129 const GEOSGeometry *ring = GEOSGetExteriorRing_r( geosinit.ctxt, geos );
1132 polygon->setExteriorRing( sequenceToLinestring( ring, hasZ, hasM ).release() );
1135 QVector<QgsCurve *> interiorRings;
1136 for (
int i = 0; i < GEOSGetNumInteriorRings_r( geosinit.ctxt, geos ); ++i )
1138 ring = GEOSGetInteriorRingN_r( geosinit.ctxt, geos, i );
1141 interiorRings.push_back( sequenceToLinestring( ring, hasZ, hasM ).release() );
1144 polygon->setInteriorRings( interiorRings );
1149 std::unique_ptr<QgsLineString> QgsGeos::sequenceToLinestring(
const GEOSGeometry *
geos,
bool hasZ,
bool hasM )
1151 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt,
geos );
1152 unsigned int nPoints;
1153 GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &nPoints );
1154 QVector< double > xOut;
1155 xOut.reserve( nPoints );
1156 QVector< double > yOut;
1157 yOut.reserve( nPoints );
1158 QVector< double > zOut;
1160 zOut.reserve( nPoints );
1161 QVector< double > mOut;
1163 mOut.reserve( nPoints );
1168 for (
unsigned int i = 0; i < nPoints; ++i )
1170 GEOSCoordSeq_getX_r( geosinit.ctxt, cs, i, &x );
1172 GEOSCoordSeq_getY_r( geosinit.ctxt, cs, i, &y );
1176 GEOSCoordSeq_getZ_r( geosinit.ctxt, cs, i, &z );
1181 GEOSCoordSeq_getOrdinate_r( geosinit.ctxt, cs, i, 3, &m );
1185 std::unique_ptr< QgsLineString > line(
new QgsLineString( xOut, yOut, zOut, mOut ) );
1189 int QgsGeos::numberOfGeometries( GEOSGeometry *g )
1194 int geometryType = GEOSGeomTypeId_r( geosinit.ctxt, g );
1195 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
1196 || geometryType == GEOS_POLYGON )
1200 return GEOSGetNumGeometries_r( geosinit.ctxt, g );
1213 GEOSCoordSeq_getX_r( geosinit.ctxt, cs, i, &x );
1214 GEOSCoordSeq_getY_r( geosinit.ctxt, cs, i, &y );
1217 GEOSCoordSeq_getZ_r( geosinit.ctxt, cs, i, &z );
1221 GEOSCoordSeq_getOrdinate_r( geosinit.ctxt, cs, i, 3, &m );
1257 int geosType = GEOS_GEOMETRYCOLLECTION;
1264 geosType = GEOS_MULTIPOINT;
1268 geosType = GEOS_MULTILINESTRING;
1272 geosType = GEOS_MULTIPOLYGON;
1293 return createGeosCollection( geosType, geomVector );
1300 return createGeosPoint( static_cast<const QgsPoint *>( geom ), coordDims, precision );
1304 return createGeosLinestring( static_cast<const QgsLineString *>( geom ), precision );
1308 return createGeosPolygon( static_cast<const QgsPolygon *>( geom ), precision );
1320 std::unique_ptr<QgsAbstractGeometry> QgsGeos::overlay(
const QgsAbstractGeometry *geom, Overlay op, QString *errorMsg )
const 1322 if ( !mGeos || !geom )
1338 case OverlayIntersection:
1339 opGeom.reset( GEOSIntersection_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) );
1341 case OverlayDifference:
1342 opGeom.reset( GEOSDifference_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) );
1346 geos::unique_ptr unionGeometry( GEOSUnion_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) );
1348 if ( unionGeometry && GEOSGeomTypeId_r( geosinit.ctxt, unionGeometry.get() ) == GEOS_MULTILINESTRING )
1350 geos::unique_ptr mergedLines( GEOSLineMerge_r( geosinit.ctxt, unionGeometry.get() ) );
1353 unionGeometry = std::move( mergedLines );
1357 opGeom = std::move( unionGeometry );
1360 case OverlaySymDifference:
1361 opGeom.reset( GEOSSymDifference_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) );
1368 catch ( GEOSException &e )
1372 *errorMsg = e.what();
1378 bool QgsGeos::relation(
const QgsAbstractGeometry *geom, Relation r, QString *errorMsg )
const 1380 if ( !mGeos || !geom )
1391 bool result =
false;
1394 if ( mGeosPrepared )
1398 case RelationIntersects:
1399 result = ( GEOSPreparedIntersects_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1401 case RelationTouches:
1402 result = ( GEOSPreparedTouches_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1404 case RelationCrosses:
1405 result = ( GEOSPreparedCrosses_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1407 case RelationWithin:
1408 result = ( GEOSPreparedWithin_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1410 case RelationContains:
1411 result = ( GEOSPreparedContains_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1413 case RelationDisjoint:
1414 result = ( GEOSPreparedDisjoint_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1416 case RelationOverlaps:
1417 result = ( GEOSPreparedOverlaps_r( geosinit.ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1427 case RelationIntersects:
1428 result = ( GEOSIntersects_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1430 case RelationTouches:
1431 result = ( GEOSTouches_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1433 case RelationCrosses:
1434 result = ( GEOSCrosses_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1436 case RelationWithin:
1437 result = ( GEOSWithin_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1439 case RelationContains:
1440 result = ( GEOSContains_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1442 case RelationDisjoint:
1443 result = ( GEOSDisjoint_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1445 case RelationOverlaps:
1446 result = ( GEOSOverlaps_r( geosinit.ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1452 catch ( GEOSException &e )
1456 *errorMsg = e.what();
1474 geos.reset( GEOSBuffer_r( geosinit.ctxt, mGeos.get(),
distance, segments ) );
1477 return fromGeos( geos.get() ).release();
1490 geos.reset( GEOSBufferWithStyle_r( geosinit.ctxt, mGeos.get(),
distance, segments, endCapStyle, joinStyle, miterLimit ) );
1493 return fromGeos( geos.get() ).release();
1505 geos.reset( GEOSTopologyPreserveSimplify_r( geosinit.ctxt, mGeos.get(), tolerance ) );
1508 return fromGeos( geos.get() ).release();
1520 geos.reset( GEOSInterpolate_r( geosinit.ctxt, mGeos.get(),
distance ) );
1523 return fromGeos( geos.get() ).release();
1539 geos.reset( GEOSGetCentroid_r( geosinit.ctxt, mGeos.get() ) );
1544 GEOSGeomGetX_r( geosinit.ctxt, geos.get(), &x );
1545 GEOSGeomGetY_r( geosinit.ctxt, geos.get(), &y );
1561 geos.reset( GEOSEnvelope_r( geosinit.ctxt, mGeos.get() ) );
1564 return fromGeos( geos.get() ).release();
1580 geos.reset( GEOSPointOnSurface_r( geosinit.ctxt, mGeos.get() ) );
1582 if ( !geos || GEOSisEmpty_r( geosinit.ctxt, geos.get() ) != 0 )
1587 GEOSGeomGetX_r( geosinit.ctxt, geos.get(), &x );
1588 GEOSGeomGetY_r( geosinit.ctxt, geos.get(), &y );
1605 std::unique_ptr< QgsAbstractGeometry > cHullGeom =
fromGeos( cHull.get() );
1606 return cHullGeom.release();
1620 return GEOSisValid_r( geosinit.ctxt, mGeos.get() );
1627 if ( !mGeos || !geom )
1639 bool equal = GEOSEquals_r( geosinit.ctxt, mGeos.get(), geosGeom.get() );
1654 return GEOSisEmpty_r( geosinit.ctxt, mGeos.get() );
1668 return GEOSisSimple_r( geosinit.ctxt, mGeos.get() );
1673 GEOSCoordSequence *QgsGeos::createCoordinateSequence(
const QgsCurve *curve,
double precision,
bool forceClose )
1675 std::unique_ptr< QgsLineString > segmentized;
1681 line = segmentized.get();
1689 bool hasZ = line->
is3D();
1703 int numOutPoints = numPoints;
1704 if ( forceClose && ( line->
pointN( 0 ) != line->
pointN( numPoints - 1 ) ) )
1709 GEOSCoordSequence *coordSeq =
nullptr;
1712 coordSeq = GEOSCoordSeq_create_r( geosinit.ctxt, numOutPoints, coordDims );
1715 QgsMessageLog::logMessage( QObject::tr(
"Could not create coordinate sequence for %1 points in %2 dimensions" ).arg( numPoints ).arg( coordDims ), QObject::tr(
"GEOS" ) );
1718 if ( precision > 0. )
1720 for (
int i = 0; i < numOutPoints; ++i )
1722 GEOSCoordSeq_setX_r( geosinit.ctxt, coordSeq, i, std::round( line->
xAt( i % numPoints ) / precision ) * precision );
1723 GEOSCoordSeq_setY_r( geosinit.ctxt, coordSeq, i, std::round( line->
yAt( i % numPoints ) / precision ) * precision );
1726 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, i, 2, std::round( line->
zAt( i % numPoints ) / precision ) * precision );
1730 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, i, 3, line->
mAt( i % numPoints ) );
1736 for (
int i = 0; i < numOutPoints; ++i )
1738 GEOSCoordSeq_setX_r( geosinit.ctxt, coordSeq, i, line->
xAt( i % numPoints ) );
1739 GEOSCoordSeq_setY_r( geosinit.ctxt, coordSeq, i, line->
yAt( i % numPoints ) );
1742 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, i, 2, line->
zAt( i % numPoints ) );
1746 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, i, 3, line->
mAt( i % numPoints ) );
1762 return createGeosPointXY( pt->
x(), pt->
y(), pt->
is3D(), pt->
z(), pt->
isMeasure(), pt->
m(), coordDims, precision );
1765 geos::unique_ptr QgsGeos::createGeosPointXY(
double x,
double y,
bool hasZ,
double z,
bool hasM,
double m,
int coordDims,
double precision )
1774 GEOSCoordSequence *coordSeq = GEOSCoordSeq_create_r( geosinit.ctxt, 1, coordDims );
1777 QgsMessageLog::logMessage( QObject::tr(
"Could not create coordinate sequence for point with %1 dimensions" ).arg( coordDims ), QObject::tr(
"GEOS" ) );
1780 if ( precision > 0. )
1782 GEOSCoordSeq_setX_r( geosinit.ctxt, coordSeq, 0, std::round( x / precision ) * precision );
1783 GEOSCoordSeq_setY_r( geosinit.ctxt, coordSeq, 0, std::round( y / precision ) * precision );
1786 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, 0, 2, std::round( z / precision ) * precision );
1791 GEOSCoordSeq_setX_r( geosinit.ctxt, coordSeq, 0, x );
1792 GEOSCoordSeq_setY_r( geosinit.ctxt, coordSeq, 0, y );
1795 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, 0, 2, z );
1798 #if 0 //disabled until geos supports m-coordinates 1801 GEOSCoordSeq_setOrdinate_r( geosinit.ctxt, coordSeq, 0, 3, m );
1804 geosPoint.reset( GEOSGeom_createPoint_r( geosinit.ctxt, coordSeq ) );
1816 GEOSCoordSequence *coordSeq = createCoordinateSequence( c, precision );
1823 geosGeom.reset( GEOSGeom_createLineString_r( geosinit.ctxt, coordSeq ) );
1836 if ( !exteriorRing )
1844 geos::unique_ptr exteriorRingGeos( GEOSGeom_createLinearRing_r( geosinit.ctxt, createCoordinateSequence( exteriorRing, precision,
true ) ) );
1847 GEOSGeometry **holes =
nullptr;
1850 holes =
new GEOSGeometry*[ nHoles ];
1853 for (
int i = 0; i < nHoles; ++i )
1856 holes[i] = GEOSGeom_createLinearRing_r( geosinit.ctxt, createCoordinateSequence( interiorRing, precision,
true ) );
1858 geosPolygon.reset( GEOSGeom_createPolygon_r( geosinit.ctxt, exteriorRingGeos.release(), holes, nHoles ) );
1874 offset.reset( GEOSOffsetCurve_r( geosinit.ctxt, mGeos.get(),
distance, segments, joinStyle, miterLimit ) );
1877 std::unique_ptr< QgsAbstractGeometry > offsetGeom =
fromGeos( offset.get() );
1878 return offsetGeom.release();
1892 GEOSBufferParams_setSingleSided_r( geosinit.ctxt, bp.get(), 1 );
1893 GEOSBufferParams_setQuadrantSegments_r( geosinit.ctxt, bp.get(), segments );
1894 GEOSBufferParams_setJoinStyle_r( geosinit.ctxt, bp.get(), joinStyle );
1895 GEOSBufferParams_setMitreLimit_r( geosinit.ctxt, bp.get(), miterLimit );
1901 geos.reset( GEOSBufferWithParams_r( geosinit.ctxt, mGeos.get(), bp.get(),
distance ) );
1921 geos::unique_ptr reshapeLineGeos = createGeosLinestring( &reshapeWithLine, mPrecision );
1924 int numGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos.get() );
1925 if ( numGeoms == -1 )
1934 bool isMultiGeom =
false;
1935 int geosTypeId = GEOSGeomTypeId_r( geosinit.ctxt, mGeos.get() );
1936 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
1946 reshapedGeometry = reshapeLine( mGeos.get(), reshapeLineGeos.get(), mPrecision );
1950 reshapedGeometry = reshapePolygon( mGeos.get(), reshapeLineGeos.get(), mPrecision );
1955 std::unique_ptr< QgsAbstractGeometry > reshapeResult =
fromGeos( reshapedGeometry.get() );
1956 return reshapeResult;
1963 bool reshapeTookPlace =
false;
1966 GEOSGeometry **newGeoms =
new GEOSGeometry*[numGeoms];
1968 for (
int i = 0; i < numGeoms; ++i )
1971 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( geosinit.ctxt, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
1973 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( geosinit.ctxt, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
1975 if ( currentReshapeGeometry )
1977 newGeoms[i] = currentReshapeGeometry.release();
1978 reshapeTookPlace =
true;
1982 newGeoms[i] = GEOSGeom_clone_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, mGeos.get(), i ) );
1989 newMultiGeom.reset( GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTILINESTRING, newGeoms, numGeoms ) );
1993 newMultiGeom.reset( GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTIPOLYGON, newGeoms, numGeoms ) );
1997 if ( !newMultiGeom )
2003 if ( reshapeTookPlace )
2007 std::unique_ptr< QgsAbstractGeometry > reshapedMultiGeom =
fromGeos( newMultiGeom.get() );
2008 return reshapedMultiGeom;
2030 if ( GEOSGeomTypeId_r( geosinit.ctxt, mGeos.get() ) != GEOS_MULTILINESTRING )
2036 geos.reset( GEOSLineMerge_r( geosinit.ctxt, mGeos.get() ) );
2044 if ( !mGeos || other.
isNull() )
2061 ( void )GEOSCoordSeq_getX_r( geosinit.ctxt, nearestCoord.get(), 0, &nx );
2062 ( void )GEOSCoordSeq_getY_r( geosinit.ctxt, nearestCoord.get(), 0, &ny );
2064 catch ( GEOSException &e )
2068 *errorMsg = e.what();
2078 if ( !mGeos || other.
isNull() )
2097 ( void )GEOSCoordSeq_getX_r( geosinit.ctxt, nearestCoord.get(), 0, &nx1 );
2098 ( void )GEOSCoordSeq_getY_r( geosinit.ctxt, nearestCoord.get(), 0, &ny1 );
2099 ( void )GEOSCoordSeq_getX_r( geosinit.ctxt, nearestCoord.get(), 1, &nx2 );
2100 ( void )GEOSCoordSeq_getY_r( geosinit.ctxt, nearestCoord.get(), 1, &ny2 );
2102 catch ( GEOSException &e )
2106 *errorMsg = e.what();
2133 distance = GEOSProject_r( geosinit.ctxt, mGeos.get(), otherGeom.get() );
2135 catch ( GEOSException &e )
2139 *errorMsg = e.what();
2149 GEOSGeometry **
const lineGeosGeometries =
new GEOSGeometry*[ geometries.size()];
2156 lineGeosGeometries[validLines] = l.release();
2163 geos::unique_ptr result( GEOSPolygonize_r( geosinit.ctxt, lineGeosGeometries, validLines ) );
2164 for (
int i = 0; i < validLines; ++i )
2166 GEOSGeom_destroy_r( geosinit.ctxt, lineGeosGeometries[i] );
2168 delete[] lineGeosGeometries;
2171 catch ( GEOSException &e )
2175 *errorMsg = e.what();
2177 for (
int i = 0; i < validLines; ++i )
2179 GEOSGeom_destroy_r( geosinit.ctxt, lineGeosGeometries[i] );
2181 delete[] lineGeosGeometries;
2196 extentGeosGeom =
asGeos( extent, mPrecision );
2197 if ( !extentGeosGeom )
2206 geos.reset( GEOSVoronoiDiagram_r( geosinit.ctxt, mGeos.get(), extentGeosGeom.get(), tolerance, edgesOnly ) );
2208 if ( !geos || GEOSisEmpty_r( geosinit.ctxt, geos.get() ) != 0 )
2228 geos.reset( GEOSDelaunayTriangulation_r( geosinit.ctxt, mGeos.get(), tolerance, edgesOnly ) );
2230 if ( !geos || GEOSisEmpty_r( geosinit.ctxt, geos.get() ) != 0 )
2242 static bool _linestringEndpoints(
const GEOSGeometry *linestring,
double &x1,
double &y1,
double &x2,
double &y2 )
2244 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, linestring );
2248 unsigned int coordSeqSize;
2249 if ( GEOSCoordSeq_getSize_r( geosinit.ctxt, coordSeq, &coordSeqSize ) == 0 )
2252 if ( coordSeqSize < 2 )
2255 GEOSCoordSeq_getX_r( geosinit.ctxt, coordSeq, 0, &x1 );
2256 GEOSCoordSeq_getY_r( geosinit.ctxt, coordSeq, 0, &y1 );
2257 GEOSCoordSeq_getX_r( geosinit.ctxt, coordSeq, coordSeqSize - 1, &x2 );
2258 GEOSCoordSeq_getY_r( geosinit.ctxt, coordSeq, coordSeqSize - 1, &y2 );
2264 static geos::unique_ptr _mergeLinestrings(
const GEOSGeometry *line1,
const GEOSGeometry *line2,
const QgsPointXY &intersectionPoint )
2266 double x1, y1, x2, y2;
2267 if ( !_linestringEndpoints( line1, x1, y1, x2, y2 ) )
2270 double rx1, ry1, rx2, ry2;
2271 if ( !_linestringEndpoints( line2, rx1, ry1, rx2, ry2 ) )
2274 bool intersectionAtOrigLineEndpoint =
2275 ( intersectionPoint.
x() == x1 && intersectionPoint.
y() == y1 ) ||
2276 ( intersectionPoint.
x() == x2 && intersectionPoint.
y() == y2 );
2277 bool intersectionAtReshapeLineEndpoint =
2278 ( intersectionPoint.
x() == rx1 && intersectionPoint.
y() == ry1 ) ||
2279 ( intersectionPoint.
x() == rx2 && intersectionPoint.
y() == ry2 );
2282 if ( intersectionAtOrigLineEndpoint && intersectionAtReshapeLineEndpoint )
2286 GEOSGeometry *geoms[2] = { g1.release(), g2.release() };
2287 geos::unique_ptr multiGeom( GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTILINESTRING, geoms, 2 ) );
2296 geos::unique_ptr QgsGeos::reshapeLine(
const GEOSGeometry *line,
const GEOSGeometry *reshapeLineGeos,
double precision )
2298 if ( !line || !reshapeLineGeos )
2301 bool atLeastTwoIntersections =
false;
2302 bool oneIntersection =
false;
2308 geos::unique_ptr intersectGeom( GEOSIntersection_r( geosinit.ctxt, line, reshapeLineGeos ) );
2309 if ( intersectGeom )
2311 atLeastTwoIntersections = ( GEOSGeomTypeId_r( geosinit.ctxt, intersectGeom.get() ) == GEOS_MULTIPOINT
2312 && GEOSGetNumGeometries_r( geosinit.ctxt, intersectGeom.get() ) > 1 );
2314 if ( GEOSGeomTypeId_r( geosinit.ctxt, intersectGeom.get() ) == GEOS_POINT )
2316 const GEOSCoordSequence *intersectionCoordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, intersectGeom.get() );
2318 GEOSCoordSeq_getX_r( geosinit.ctxt, intersectionCoordSeq, 0, &xi );
2319 GEOSCoordSeq_getY_r( geosinit.ctxt, intersectionCoordSeq, 0, &yi );
2320 oneIntersection =
true;
2325 catch ( GEOSException &e )
2328 atLeastTwoIntersections =
false;
2332 if ( oneIntersection )
2333 return _mergeLinestrings( line, reshapeLineGeos, oneIntersectionPoint );
2335 if ( !atLeastTwoIntersections )
2339 double x1, y1, x2, y2;
2340 if ( !_linestringEndpoints( line, x1, y1, x2, y2 ) )
2343 geos::unique_ptr beginLineVertex = createGeosPointXY( x1, y1,
false, 0,
false, 0, 2, precision );
2344 geos::unique_ptr endLineVertex = createGeosPointXY( x2, y2,
false, 0,
false, 0, 2, precision );
2346 bool isRing =
false;
2347 if ( GEOSGeomTypeId_r( geosinit.ctxt, line ) == GEOS_LINEARRING
2348 || GEOSEquals_r( geosinit.ctxt, beginLineVertex.get(), endLineVertex.get() ) == 1 )
2353 if ( !nodedGeometry )
2359 geos::unique_ptr mergedLines( GEOSLineMerge_r( geosinit.ctxt, nodedGeometry.get() ) );
2365 int numMergedLines = GEOSGetNumGeometries_r( geosinit.ctxt, mergedLines.get() );
2366 if ( numMergedLines < 2 )
2368 if ( numMergedLines == 1 )
2370 geos::unique_ptr result( GEOSGeom_clone_r( geosinit.ctxt, reshapeLineGeos ) );
2377 QVector<GEOSGeometry *> resultLineParts;
2378 QVector<GEOSGeometry *> probableParts;
2380 for (
int i = 0; i < numMergedLines; ++i )
2382 const GEOSGeometry *currentGeom =
nullptr;
2384 currentGeom = GEOSGetGeometryN_r( geosinit.ctxt, mergedLines.get(), i );
2385 const GEOSCoordSequence *currentCoordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, currentGeom );
2386 unsigned int currentCoordSeqSize;
2387 GEOSCoordSeq_getSize_r( geosinit.ctxt, currentCoordSeq, ¤tCoordSeqSize );
2388 if ( currentCoordSeqSize < 2 )
2392 double xBegin, xEnd, yBegin, yEnd;
2393 GEOSCoordSeq_getX_r( geosinit.ctxt, currentCoordSeq, 0, &xBegin );
2394 GEOSCoordSeq_getY_r( geosinit.ctxt, currentCoordSeq, 0, &yBegin );
2395 GEOSCoordSeq_getX_r( geosinit.ctxt, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
2396 GEOSCoordSeq_getY_r( geosinit.ctxt, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
2397 geos::unique_ptr beginCurrentGeomVertex = createGeosPointXY( xBegin, yBegin,
false, 0,
false, 0, 2, precision );
2398 geos::unique_ptr endCurrentGeomVertex = createGeosPointXY( xEnd, yEnd,
false, 0,
false, 0, 2, precision );
2401 int nEndpointsOnOriginalLine = 0;
2402 if ( pointContainedInLine( beginCurrentGeomVertex.get(), line ) == 1 )
2403 nEndpointsOnOriginalLine += 1;
2405 if ( pointContainedInLine( endCurrentGeomVertex.get(), line ) == 1 )
2406 nEndpointsOnOriginalLine += 1;
2409 int nEndpointsSameAsOriginalLine = 0;
2410 if ( GEOSEquals_r( geosinit.ctxt, beginCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
2411 || GEOSEquals_r( geosinit.ctxt, beginCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
2412 nEndpointsSameAsOriginalLine += 1;
2414 if ( GEOSEquals_r( geosinit.ctxt, endCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
2415 || GEOSEquals_r( geosinit.ctxt, endCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
2416 nEndpointsSameAsOriginalLine += 1;
2419 bool currentGeomOverlapsOriginalGeom =
false;
2420 bool currentGeomOverlapsReshapeLine =
false;
2421 if ( lineContainedInLine( currentGeom, line ) == 1 )
2422 currentGeomOverlapsOriginalGeom =
true;
2424 if ( lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
2425 currentGeomOverlapsReshapeLine =
true;
2428 if ( !isRing && nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
2430 resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
2433 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
2435 probableParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
2437 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
2439 resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
2441 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
2443 resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
2445 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
2447 resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
2452 if ( isRing && !probableParts.isEmpty() )
2455 GEOSGeometry *currentGeom =
nullptr;
2456 double maxLength = -std::numeric_limits<double>::max();
2457 double currentLength = 0;
2458 for (
int i = 0; i < probableParts.size(); ++i )
2460 currentGeom = probableParts.at( i );
2461 GEOSLength_r( geosinit.ctxt, currentGeom, ¤tLength );
2462 if ( currentLength > maxLength )
2464 maxLength = currentLength;
2465 maxGeom.reset( currentGeom );
2469 GEOSGeom_destroy_r( geosinit.ctxt, currentGeom );
2472 resultLineParts.push_back( maxGeom.release() );
2476 if ( resultLineParts.empty() )
2479 if ( resultLineParts.size() == 1 )
2481 result.reset( resultLineParts[0] );
2485 GEOSGeometry **lineArray =
new GEOSGeometry*[resultLineParts.size()];
2486 for (
int i = 0; i < resultLineParts.size(); ++i )
2488 lineArray[i] = resultLineParts[i];
2492 geos::unique_ptr multiLineGeom( GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ) );
2493 delete [] lineArray;
2496 result.reset( GEOSLineMerge_r( geosinit.ctxt, multiLineGeom.get() ) );
2500 if ( GEOSGeomTypeId_r( geosinit.ctxt, result.get() ) != GEOS_LINESTRING )
2508 geos::unique_ptr QgsGeos::reshapePolygon(
const GEOSGeometry *polygon,
const GEOSGeometry *reshapeLineGeos,
double precision )
2511 int nIntersections = 0;
2512 int lastIntersectingRing = -2;
2513 const GEOSGeometry *lastIntersectingGeom =
nullptr;
2515 int nRings = GEOSGetNumInteriorRings_r( geosinit.ctxt, polygon );
2520 const GEOSGeometry *outerRing = GEOSGetExteriorRing_r( geosinit.ctxt, polygon );
2521 if ( GEOSIntersects_r( geosinit.ctxt, outerRing, reshapeLineGeos ) == 1 )
2524 lastIntersectingRing = -1;
2525 lastIntersectingGeom = outerRing;
2529 const GEOSGeometry **innerRings =
new const GEOSGeometry*[nRings];
2533 for (
int i = 0; i < nRings; ++i )
2535 innerRings[i] = GEOSGetInteriorRingN_r( geosinit.ctxt, polygon, i );
2536 if ( GEOSIntersects_r( geosinit.ctxt, innerRings[i], reshapeLineGeos ) == 1 )
2539 lastIntersectingRing = i;
2540 lastIntersectingGeom = innerRings[i];
2544 catch ( GEOSException &e )
2550 if ( nIntersections != 1 )
2552 delete [] innerRings;
2557 geos::unique_ptr reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos, precision );
2558 if ( !reshapeResult )
2560 delete [] innerRings;
2565 GEOSGeometry *newRing =
nullptr;
2566 const GEOSCoordSequence *reshapeSequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, reshapeResult.get() );
2567 GEOSCoordSequence *newCoordSequence = GEOSCoordSeq_clone_r( geosinit.ctxt, reshapeSequence );
2569 reshapeResult.reset();
2571 newRing = GEOSGeom_createLinearRing_r( geosinit.ctxt, newCoordSequence );
2574 delete [] innerRings;
2578 GEOSGeometry *newOuterRing =
nullptr;
2579 if ( lastIntersectingRing == -1 )
2580 newOuterRing = newRing;
2582 newOuterRing = GEOSGeom_clone_r( geosinit.ctxt, outerRing );
2585 QVector<GEOSGeometry *> ringList;
2588 GEOSGeometry *outerRingPoly = GEOSGeom_createPolygon_r( geosinit.ctxt, GEOSGeom_clone_r( geosinit.ctxt, newOuterRing ),
nullptr, 0 );
2589 if ( outerRingPoly )
2591 GEOSGeometry *currentRing =
nullptr;
2592 for (
int i = 0; i < nRings; ++i )
2594 if ( lastIntersectingRing == i )
2595 currentRing = newRing;
2597 currentRing = GEOSGeom_clone_r( geosinit.ctxt, innerRings[i] );
2600 if ( GEOSContains_r( geosinit.ctxt, outerRingPoly, currentRing ) == 1 )
2601 ringList.push_back( currentRing );
2603 GEOSGeom_destroy_r( geosinit.ctxt, currentRing );
2606 GEOSGeom_destroy_r( geosinit.ctxt, outerRingPoly );
2609 GEOSGeometry **newInnerRings =
new GEOSGeometry*[ringList.size()];
2610 for (
int i = 0; i < ringList.size(); ++i )
2611 newInnerRings[i] = ringList.at( i );
2613 delete [] innerRings;
2615 geos::unique_ptr reshapedPolygon( GEOSGeom_createPolygon_r( geosinit.ctxt, newOuterRing, newInnerRings, ringList.size() ) );
2616 delete[] newInnerRings;
2618 return reshapedPolygon;
2621 int QgsGeos::lineContainedInLine(
const GEOSGeometry *line1,
const GEOSGeometry *line2 )
2623 if ( !line1 || !line2 )
2628 double bufferDistance = std::pow( 10.0L, geomDigits( line2 ) - 11 );
2634 geos::unique_ptr intersectionGeom( GEOSIntersection_r( geosinit.ctxt, bufferGeom.get(), line1 ) );
2637 double intersectGeomLength;
2640 GEOSLength_r( geosinit.ctxt, intersectionGeom.get(), &intersectGeomLength );
2641 GEOSLength_r( geosinit.ctxt, line1, &line1Length );
2643 double intersectRatio = line1Length / intersectGeomLength;
2644 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
2650 int QgsGeos::pointContainedInLine(
const GEOSGeometry *point,
const GEOSGeometry *line )
2652 if ( !point || !line )
2655 double bufferDistance = std::pow( 10.0L, geomDigits( line ) - 11 );
2657 geos::unique_ptr lineBuffer( GEOSBuffer_r( geosinit.ctxt, line, bufferDistance, 8 ) );
2661 bool contained =
false;
2662 if ( GEOSContains_r( geosinit.ctxt, lineBuffer.get(), point ) == 1 )
2668 int QgsGeos::geomDigits(
const GEOSGeometry *geom )
2674 const GEOSGeometry *bBoxRing = GEOSGetExteriorRing_r( geosinit.ctxt, bbox.get() );
2678 const GEOSCoordSequence *bBoxCoordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, bBoxRing );
2680 if ( !bBoxCoordSeq )
2683 unsigned int nCoords = 0;
2684 if ( !GEOSCoordSeq_getSize_r( geosinit.ctxt, bBoxCoordSeq, &nCoords ) )
2688 for (
unsigned int i = 0; i < nCoords - 1; ++i )
2691 GEOSCoordSeq_getX_r( geosinit.ctxt, bBoxCoordSeq, i, &t );
2694 digits = std::ceil( std::log10( std::fabs( t ) ) );
2695 if ( digits > maxDigits )
2698 GEOSCoordSeq_getY_r( geosinit.ctxt, bBoxCoordSeq, i, &t );
2699 digits = std::ceil( std::log10( std::fabs( t ) ) );
2700 if ( digits > maxDigits )
2709 return geosinit.ctxt;
bool isMeasure() const
Returns true if the geometry contains m values.
const QgsCurve * interiorRing(int i) const
A rectangle specified with double values.
bool disjoint(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom is disjoint from this.
QgsAbstractGeometry * intersection(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculate the intersection of this and geom.
QgsGeometry mergeLines(QString *errorMsg=nullptr) const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
double zAt(int index) const
Returns the z-coordinate of the specified node in the line string.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
QgsGeos(const QgsAbstractGeometry *geometry, double precision=0)
GEOS geometry engine constructor.
void setXMaximum(double x)
Set the maximum x value.
Multi point geometry collection.
static bool isMultiType(Type type)
Returns true if the WKB type is a multi type.
bool isEmpty(QString *errorMsg=nullptr) const override
QgsPoint * pointOnSurface(QString *errorMsg=nullptr) const override
Calculate a point that is guaranteed to be on the surface of this.
QgsGeometry shortestLine(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the shortest line joining this geometry to the other geometry.
double area(QString *errorMsg=nullptr) const override
A class to represent a 2D point.
const QgsAbstractGeometry * mGeometry
QgsAbstractGeometry * buffer(double distance, int segments, QString *errorMsg=nullptr) const override
Multi line string geometry collection.
Curve polygon geometry type.
static std::unique_ptr< QgsGeometryCollection > createCollectionOfType(QgsWkbTypes::Type type)
Returns a new geometry collection matching a specified WKB type.
double length(QString *errorMsg=nullptr) const override
A geometry is the spatial representation of a feature.
bool isEqual(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if this is equal to geom.
bool isValid(QString *errorMsg=nullptr) const override
QgsGeometry closestPoint(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the closest point on the geometry to the other geometry.
virtual QgsRectangle boundingBox() const =0
Returns the minimal bounding box for the geometry.
void CORE_EXPORT operator()(GEOSGeometry *geom)
Destroys the GEOS geometry geom, using the static QGIS geos context.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning)
add a message to the instance (and create it if necessary)
std::unique_ptr< QgsAbstractGeometry > singleSidedBuffer(double distance, int segments, int side, int joinStyle, double miterLimit, QString *errorMsg=nullptr) const
Returns a single sided buffer for a geometry.
static GEOSContextHandle_t getGEOSHandler()
int numPoints() const override
Returns the number of points in the curve.
QgsAbstractGeometry * interpolate(double distance, QString *errorMsg=nullptr) const override
bool within(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom is within this.
#define CATCH_GEOS_WITH_ERRMSG(r)
int numInteriorRings() const
Nothing happened, without any error.
Type
The WKB type describes the number of dimensions a geometry has.
QgsAbstractGeometry * envelope(QString *errorMsg=nullptr) const override
std::unique_ptr< QgsAbstractGeometry > reshapeGeometry(const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg=nullptr) const
Reshapes the geometry using a line.
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.
bool touches(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom touches this.
bool isEmpty() const
Returns true if the rectangle is empty.
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...
void prepareGeometry() override
Prepares the geometry, so that subsequent calls to spatial relation methods are much faster...
double width() const
Returns the width of the rectangle.
double mAt(int index) const
Returns the m value of the specified node in the line string.
void setYMinimum(double y)
Set the minimum y value.
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
static GeometryType geometryType(Type type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
double hausdorffDistanceDensify(const QgsAbstractGeometry *geom, double densifyFraction, QString *errorMsg=nullptr) const
Returns the Hausdorff distance between this geometry and geom.
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.
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
static geos::unique_ptr asGeos(const QgsAbstractGeometry *geom, double precision=0)
Multi curve geometry collection.
double distance(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculates the distance between this and geom.
static QgsGeometry polygonize(const QVector< const QgsAbstractGeometry *> &geometries, QString *errorMsg=nullptr)
Creates a GeometryCollection geometry containing possible polygons formed from the constituent linewo...
QgsPoint * clone() const override
Clones the geometry by performing a deep copy.
Abstract base class for curved geometry type.
Abstract base class for all geometries.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
const QgsCurve * exteriorRing() const
std::unique_ptr< GEOSBufferParams, GeosDeleter > buffer_params_unique_ptr
Scoped GEOS buffer params pointer.
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
QgsGeometry delaunayTriangulation(double tolerance=0.0, bool edgesOnly=false, QString *errorMsg=nullptr) const
Returns the Delaunay triangulation for the vertices of the geometry.
Error occurred while creating a noded geometry.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
Error occurred in the geometry engine.
int numGeometries() const
Returns the number of geometries within the collection.
QgsAbstractGeometry * symDifference(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculate the symmetric difference of this and geom.
Contains geos related utilities and functions.
std::unique_ptr< QgsAbstractGeometry > clip(const QgsRectangle &rectangle, QString *errorMsg=nullptr) const
Performs a fast, non-robust intersection between the geometry and a rectangle.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
QgsAbstractGeometry * combine(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculate the combination of this and geom.
double xMaximum() const
Returns the x maximum value (right side of rectangle).
void geometryChanged() override
Should be called whenever the geometry associated with the engine has been modified and the engine mu...
QVector< QgsPoint > QgsPointSequence
static QgsPoint coordSeqPoint(const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM)
QgsAbstractGeometry * offsetCurve(double distance, int segments, int joinStyle, double miterLimit, QString *errorMsg=nullptr) const override
QgsAbstractGeometry * convexHull(QString *errorMsg=nullptr) const override
Calculate the convex hull of this.
Multi polygon geometry collection.
bool addGeometry(QgsAbstractGeometry *g) override
Adds a geometry and takes ownership. Returns true in case of success.
bool contains(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom contains this.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
void setYMaximum(double y)
Set the maximum y value.
std::unique_ptr< QgsAbstractGeometry > subdivide(int maxNodes, QString *errorMsg=nullptr) const
Subdivides the geometry.
Line string geometry type, with support for z-dimension and m-values.
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...
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 isNull() const
Test if the rectangle is null (all coordinates zero or after call to setMinimal()).
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
EngineOperationResult splitGeometry(const QgsLineString &splitLine, QVector< QgsGeometry > &newGeometries, bool topological, QgsPointSequence &topologyTestPoints, QString *errorMsg=nullptr) const override
Splits this geometry according to a given line.
bool intersects(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom intersects this.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double hausdorffDistance(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const
Returns the Hausdorff distance between this geometry and geom.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
QgsAbstractGeometry * difference(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculate the difference of this and geom.
bool overlaps(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom overlaps this.
Contains geometry relation and modification algorithms.
double yMaximum() const
Returns the y maximum value (top side of rectangle).
static std::unique_ptr< QgsPolygon > fromGeosPolygon(const GEOSGeometry *geos)
std::unique_ptr< GEOSCoordSequence, GeosDeleter > coord_sequence_unique_ptr
Scoped GEOS coordinate sequence pointer.
EngineOperationResult
Success or failure of a geometry operation.
QgsGeometry 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.
static Type flatType(Type type)
Returns the flat type for a WKB type.
The geometry on which the operation occurs is not valid.
#define DEFAULT_QUADRANT_SEGMENTS
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
bool isSimple(QString *errorMsg=nullptr) const override
Determines whether the geometry is simple (according to OGC definition).
bool crosses(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom crosses this.
void setXMinimum(double x)
Set the minimum x value.
QgsAbstractGeometry * simplify(double tolerance, QString *errorMsg=nullptr) const override
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
double height() const
Returns the height of the rectangle.