Fix elevation profile on intersecting polyline

python/PyQt6/core/auto_generated/geometry/qgslinestring.sip.in

@@ -643,6 +643,17 @@ If ``useZValues`` is ``True`` then z values will also be considered when testing
 
 
 
+    QVector<QgsLineString *> splitToDisjointXYParts() const /Factory/;
+%Docstring
+Divides the linestring into parts that don't share any points or lines.
+
+This method throws away Z and M coordinates.
+
+The ownership of returned pointers is transferred to the caller.
+
+.. versionadded:: 3.40
+%End
+
     double length3D() const /HoldGIL/;
 %Docstring
 Returns the length in 3D world of the line string.

python/core/auto_generated/geometry/qgslinestring.sip.in

@@ -643,6 +643,17 @@ If ``useZValues`` is ``True`` then z values will also be considered when testing
 
 
 
+    QVector<QgsLineString *> splitToDisjointXYParts() const /Factory/;
+%Docstring
+Divides the linestring into parts that don't share any points or lines.
+
+This method throws away Z and M coordinates.
+
+The ownership of returned pointers is transferred to the caller.
+
+.. versionadded:: 3.40
+%End
+
     double length3D() const /HoldGIL/;
 %Docstring
 Returns the length in 3D world of the line string.

src/core/geometry/qgslinestring.cpp

@@ -1086,6 +1086,40 @@ std::tuple<std::unique_ptr<QgsCurve>, std::unique_ptr<QgsCurve> > QgsLineString:
     return std::make_tuple( std::make_unique< QgsLineString >( x1, y1, z1, m1 ), std::make_unique< QgsLineString >( x2, y2, z2, m2 ) );
 }
 
+QVector<QgsLineString *> QgsLineString::splitToDisjointXYParts() const
+{
+  const double *allPointsX = xData();
+  const double *allPointsY = yData();
+  size_t allPointsCount = numPoints();
+  QVector<double> partX;
+  QVector<double> partY;
+  QSet<QgsPointXY> partPointSet;
+
+  QVector<QgsLineString *> disjointParts;
+  for ( size_t i = 0; i < allPointsCount; i++ )
+  {
+    const QgsPointXY point( *allPointsX++, *allPointsY++ );
+    if ( partPointSet.contains( point ) )
+    {
+      // This point is used multiple times, cut the curve and add the
+      // current part
+      disjointParts.push_back( new QgsLineString( partX, partY ) );
+      // Now start a new part containing the last line
+      partX = { partX.last() };
+      partY = { partY.last() };
+      partPointSet = { QgsPointXY( partX[0], partY[0] ) };
+    }
+    partX.push_back( point.x() );
+    partY.push_back( point.y() );
+    partPointSet.insert( point );
+  }
+  // Add the last part (if we didn't stop by closing the loop)
+  if ( partX.size() > 1 || disjointParts.size() == 0 )
+    disjointParts.push_back( new QgsLineString( partX, partY ) );
+
+  return disjointParts;
+}
+
 double QgsLineString::length3D() const
 {
   if ( is3D() )

src/core/geometry/qgslinestring.h

@@ -991,6 +991,16 @@ class CORE_EXPORT QgsLineString: public QgsCurve
     std::tuple< std::unique_ptr< QgsCurve >, std::unique_ptr< QgsCurve > > splitCurveAtVertex( int index ) const final;
 #endif
 
+    /**
+     * Divides the linestring into parts that don't share any points or lines.
+     *
+     * This method throws away Z and M coordinates.
+     *
+     * The ownership of returned pointers is transferred to the caller.
+     * \since QGIS 3.40
+     */
+    QVector<QgsLineString *> splitToDisjointXYParts() const SIP_FACTORY;
+
     /**
      * Returns the length in 3D world of the line string.
      * If it is not a 3D line string, return its 2D length.

src/core/vector/qgsvectorlayerprofilegenerator.cpp

@@ -15,6 +15,7 @@
  *                                                                         *
  ***************************************************************************/
 #include "qgsvectorlayerprofilegenerator.h"
+#include "qgsabstractgeometry.h"
 #include "qgspolyhedralsurface.h"
 #include "qgsprofilerequest.h"
 #include "qgscurve.h"
@@ -721,6 +722,72 @@ bool QgsVectorLayerProfileGenerator::generateProfile( const QgsProfileGeneration
   if ( !mProfileCurve || mFeedback->isCanceled() )
     return false;
 
+  if ( QgsLineString *profileLine =
+         qgsgeometry_cast<QgsLineString *>( mProfileCurve.get() ) )
+  {
+    // The profile generation code can't deal with curves that enter a single
+    // point multiple times. We handle this for line strings by splitting them
+    // into multiple parts, each with no repeated points, and computing the
+    // profile for each by itself.
+    std::unique_ptr< QgsCurve > origCurve = std::move( mProfileCurve );
+    std::unique_ptr< QgsVectorLayerProfileResults > totalResults;
+    double distanceProcessed = 0;
+
+    QVector<QgsLineString *> disjointParts = profileLine->splitToDisjointXYParts();
+    for ( int i = 0; i < disjointParts.size(); i++ )
+    {
+      mProfileCurve.reset( disjointParts[i] );
+      if ( !generateProfileInner() )
+      {
+        mProfileCurve = std::move( origCurve );
+
+        // Free the rest of the parts
+        for ( int j = i + 1; j < disjointParts.size(); j++ )
+          delete disjointParts[j];
+
+        return false;
+      }
+
+      if ( !totalResults )
+        // Use the first result set as a base
+        totalResults.reset( mResults.release() );
+      else
+      {
+        // Merge the results, shifting them by distanceProcessed
+        totalResults->mRawPoints.append( mResults->mRawPoints );
+        totalResults->minZ = std::min( totalResults->minZ, mResults->minZ );
+        totalResults->maxZ = std::max( totalResults->maxZ, mResults->maxZ );
+        for ( auto it = mResults->mDistanceToHeightMap.constKeyValueBegin();
+              it != mResults->mDistanceToHeightMap.constKeyValueEnd();
+              ++it )
+        {
+          totalResults->mDistanceToHeightMap[it->first + distanceProcessed] = it->second;
+        }
+        for ( auto it = mResults->features.constKeyValueBegin();
+              it != mResults->features.constKeyValueEnd();
+              ++it )
+        {
+          for ( QgsVectorLayerProfileResults::Feature feature : it->second )
+          {
+            feature.crossSectionGeometry.translate( distanceProcessed, 0 );
+            totalResults->features[it->first].push_back( feature );
+          }
+        }
+      }
+
+      distanceProcessed += mProfileCurve->length();
+    }
+
+    mProfileCurve = std::move( origCurve );
+    mResults.reset( totalResults.release() );
+    return true;
+  }
+
+  return generateProfileInner();
+}
+
+bool QgsVectorLayerProfileGenerator::generateProfileInner( const QgsProfileGenerationContext & )
+{
   // we need to transform the profile curve to the vector's CRS
   mTransformedCurve.reset( mProfileCurve->clone() );
   mLayerToTargetTransform = QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext );

src/core/vector/qgsvectorlayerprofilegenerator.h

@@ -120,6 +120,10 @@ class CORE_EXPORT QgsVectorLayerProfileGenerator : public QgsAbstractProfileSurf
 
   private:
 
+    // We may need to split mProfileCurve into multiple parts, this will be
+    // called for each part.
+    bool generateProfileInner( const QgsProfileGenerationContext &context = QgsProfileGenerationContext() );
+
     bool generateProfileForPoints();
     bool generateProfileForLines();
     bool generateProfileForPolygons();

tests/src/core/geometry/testqgsgeometry.cpp

@@ -89,6 +89,7 @@ class TestQgsGeometry : public QgsTest
     void curveIndexOf();
     void splitCurve_data();
     void splitCurve();
+    void splitToDisjointXYParts();
 
     void fromBox3d();
     void fromPoint();
@@ -751,6 +752,38 @@ void TestQgsGeometry::splitCurve()
   QCOMPARE( p2->asWkt(), curve2 );
 }
 
+void TestQgsGeometry::splitToDisjointXYParts()
+{
+  QgsLineString onePartLine( QgsPoint( 1.0, 1.0 ), QgsPoint( 2.0, 2.0 ) );
+  QVector<QgsLineString *> onePartParts = onePartLine.splitToDisjointXYParts();
+  QCOMPARE( onePartParts.size(), 1 );
+  QCOMPARE( onePartParts[0]->asWkt(), onePartLine.asWkt() );
+  for ( QgsLineString *part : onePartParts )
+    delete part;
+
+  QgsLineString onePointLine( QVector<QgsPoint> {QgsPoint( 1.0, 1.0 )} );
+  QVector<QgsLineString *> onePointParts = onePointLine.splitToDisjointXYParts();
+  QCOMPARE( onePointParts.size(), 1 );
+  QCOMPARE( onePointParts[0]->asWkt(), onePointLine.asWkt() );
+  for ( QgsLineString *part : onePointParts )
+    delete part;
+
+  QgsLineString emptyLine( QVector<QgsPoint> { } );
+  QVector<QgsLineString *> emptyParts = emptyLine.splitToDisjointXYParts();
+  QCOMPARE( emptyParts.size(), 1 );
+  QCOMPARE( emptyParts[0]->asWkt(), emptyLine.asWkt() );
+  for ( QgsLineString *part : emptyParts )
+    delete part;
+
+  QgsLineString triangle( QVector<QgsPoint> {QgsPoint( 0.0, 0.0 ), QgsPoint( 1.0, 0.0 ), QgsPoint( 1.0, 1.0 ), QgsPoint( 0.0, 0.0 )} );
+  QVector<QgsLineString *> triangleParts = triangle.splitToDisjointXYParts();
+  QCOMPARE( triangleParts.size(), 2 );
+  QCOMPARE( triangleParts[0]->asWkt(), "LineString (0 0, 1 0, 1 1)" );
+  QCOMPARE( triangleParts[1]->asWkt(), "LineString (1 1, 0 0)" );
+  for ( QgsLineString *part : triangleParts )
+    delete part;
+}
+
 void TestQgsGeometry::fromBox3d()
 {
   QgsBox3D box3d( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 );