Add geometry methods for interpolating angle along geometry

Sponsored by Andreas Neumann

(cherry-picked from 93c7f5f)

python/core/geometry/qgsgeometry.sip

@@ -488,6 +488,15 @@ class QgsGeometry
      */
     double lineLocatePoint( const QgsGeometry& point ) const;
 
+    /** Returns the angle parallel to the linestring or polygon boundary at the specified distance
+     * along the geometry. Angles are in radians, clockwise from north.
+     * If the distance coincides precisely at a node then the average angle from the segment either side
+     * of the node is returned.
+     * @param distance distance along geometry
+     * @note added in QGIS 3.0
+     */
+    double interpolateAngle( double distance ) const;
+
     /** Returns a geometry representing the points shared by this geometry and other. */
     QgsGeometry* intersection( const QgsGeometry* geometry ) const /Factory/;
 

src/core/geometry/qgsgeometry.cpp

@@ -1409,6 +1409,61 @@ double QgsGeometry::lineLocatePoint( const QgsGeometry& point ) const
   return geos.lineLocatePoint( *( static_cast< QgsPointV2* >( point.d->geometry ) ) );
 }
 
+double QgsGeometry::interpolateAngle( double distance ) const
+{
+  if ( !d->geometry )
+    return 0.0;
+
+  // always operate on segmentized geometries
+  QgsGeometry segmentized = *this;
+  if ( QgsWKBTypes::isCurvedType( d->geometry->wkbType() ) )
+  {
+    segmentized = QgsGeometry( static_cast< QgsCurveV2* >( d->geometry )->segmentize() );
+  }
+
+  QgsVertexId previous;
+  QgsVertexId next;
+  if ( !QgsGeometryUtils::verticesAtDistance( *segmentized.geometry(), distance, previous, next ) )
+    return 0.0;
+
+  if ( previous == next )
+  {
+    // distance coincided exactly with a vertex
+    QgsVertexId v2 = previous;
+    QgsVertexId v1;
+    QgsVertexId v3;
+    QgsGeometryUtils::adjacentVertices( *segmentized.geometry(), v2, v1, v3 );
+    if ( v1.isValid() && v3.isValid() )
+    {
+      QgsPointV2 p1 = segmentized.geometry()->vertexAt( v1 );
+      QgsPointV2 p2 = segmentized.geometry()->vertexAt( v2 );
+      QgsPointV2 p3 = segmentized.geometry()->vertexAt( v3 );
+      double angle1 = QgsGeometryUtils::lineAngle( p1.x(), p1.y(), p2.x(), p2.y() );
+      double angle2 = QgsGeometryUtils::lineAngle( p2.x(), p2.y(), p3.x(), p3.y() );
+      return QgsGeometryUtils::averageAngle( angle1, angle2 );
+    }
+    else if ( v3.isValid() )
+    {
+      QgsPointV2 p1 = segmentized.geometry()->vertexAt( v2 );
+      QgsPointV2 p2 = segmentized.geometry()->vertexAt( v3 );
+      return QgsGeometryUtils::lineAngle( p1.x(), p1.y(), p2.x(), p2.y() );
+    }
+    else
+    {
+      QgsPointV2 p1 = segmentized.geometry()->vertexAt( v1 );
+      QgsPointV2 p2 = segmentized.geometry()->vertexAt( v2 );
+      return QgsGeometryUtils::lineAngle( p1.x(), p1.y(), p2.x(), p2.y() );
+    }
+  }
+  else
+  {
+    QgsPointV2 p1 = segmentized.geometry()->vertexAt( previous );
+    QgsPointV2 p2 = segmentized.geometry()->vertexAt( next );
+    return QgsGeometryUtils::lineAngle( p1.x(), p1.y(), p2.x(), p2.y() );
+  }
+}
+
+
 QgsGeometry* QgsGeometry::intersection( const QgsGeometry* geometry ) const
 {
   if ( !d->geometry || !geometry->d->geometry )

src/core/geometry/qgsgeometry.h

@@ -531,6 +531,15 @@ class CORE_EXPORT QgsGeometry
      */
     double lineLocatePoint( const QgsGeometry& point ) const;
 
+    /** Returns the angle parallel to the linestring or polygon boundary at the specified distance
+     * along the geometry. Angles are in radians, clockwise from north.
+     * If the distance coincides precisely at a node then the average angle from the segment either side
+     * of the node is returned.
+     * @param distance distance along geometry
+     * @note added in QGIS 3.0
+     */
+    double interpolateAngle( double distance ) const;
+
     /** Returns a geometry representing the points shared by this geometry and other. */
     QgsGeometry* intersection( const QgsGeometry* geometry ) const;