Add geometry methods for interpolating angle along geometry

Sponsored by Andreas Neumann

python/core/geometry/qgsgeometry.sip

@@ -545,6 +545,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;
 

src/core/geometry/qgsgeometry.cpp

@@ -1510,6 +1510,60 @@ 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( wkbType() ) )
+  {
+    segmentized = QgsGeometry( static_cast< QgsCurve* >( 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.isEmpty() )

src/core/geometry/qgsgeometry.h

@@ -582,6 +582,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;
 

src/core/geometry/qgsgeometryutils.cpp

@@ -119,6 +119,51 @@ double QgsGeometryUtils::distanceToVertex( const QgsAbstractGeometry &geom, cons
   return -1;
 }
 
+bool QgsGeometryUtils::verticesAtDistance( const QgsAbstractGeometry& geometry, double distance, QgsVertexId& previousVertex, QgsVertexId& nextVertex )
+{
+  double currentDist = 0;
+  previousVertex = QgsVertexId();
+  nextVertex = QgsVertexId();
+
+  QgsPointV2 point;
+  QgsPointV2 previousPoint;
+
+  if ( qgsDoubleNear( distance, 0.0 ) )
+  {
+    geometry.nextVertex( previousVertex, point );
+    nextVertex = previousVertex;
+    return true;
+  }
+
+  bool first = true;
+  while ( currentDist < distance && geometry.nextVertex( nextVertex, point ) )
+  {
+    if ( !first )
+    {
+      currentDist += sqrt( QgsGeometryUtils::sqrDistance2D( previousPoint, point ) );
+    }
+
+    if ( qgsDoubleNear( currentDist, distance ) )
+    {
+      // exact hit!
+      previousVertex = nextVertex;
+      return true;
+    }
+
+    if ( currentDist > distance )
+    {
+      return true;
+    }
+
+    previousVertex = nextVertex;
+    previousPoint = point;
+    first = false;
+  }
+
+  //could not find target distance
+  return false;
+}
+
 void QgsGeometryUtils::adjacentVertices( const QgsAbstractGeometry& geom, QgsVertexId atVertex, QgsVertexId& beforeVertex, QgsVertexId& afterVertex )
 {
   bool polygonType = ( geom.dimension()  == 2 );

src/core/geometry/qgsgeometryutils.h

@@ -49,6 +49,21 @@ class CORE_EXPORT QgsGeometryUtils
      */
     static double distanceToVertex( const QgsAbstractGeometry& geom, const QgsVertexId& id );
 
+    /** Retrieves the vertices which are before and after the interpolated point at a specified distance along a linestring
+     * (or polygon boundary).
+     * @param geometry line or polygon geometry
+     * @param distance distance to traverse along geometry
+     * @param previousVertex will be set to previous vertex ID
+     * @param nextVertex will be set to next vertex ID
+     * @note if the distance coincides exactly with a vertex, then both previousVertex and nextVertex will be set to this vertex
+     * @returns true if vertices were successfully retrieved
+     * @note added in QGIS 3.0
+     */
+    static bool verticesAtDistance( const QgsAbstractGeometry& geometry,
+                                    double distance,
+                                    QgsVertexId& previousVertex,
+                                    QgsVertexId& nextVertex );
+
     /** Returns vertices adjacent to a specified vertex within a geometry.
      */
     static void adjacentVertices( const QgsAbstractGeometry& geom, QgsVertexId atVertex, QgsVertexId& beforeVertex, QgsVertexId& afterVertex );

tests/src/core/testqgsgeometryutils.cpp

@@ -44,7 +44,9 @@ class TestQgsGeometryUtils: public QObject
     void testLinePerpendicularAngle();
     void testAverageAngle_data();
     void testAverageAngle();
+    void testAdjacentVertices();
     void testDistanceToVertex();
+    void testVerticesAtDistance();
     void testCircleCenterRadius_data();
     void testCircleCenterRadius();
     void testSqrDistToLine();
@@ -335,6 +337,21 @@ void TestQgsGeometryUtils::testAverageAngle()
   QVERIFY( qgsDoubleNear( averageAngle, expected, 0.0000000001 ) );
 }
 
+void TestQgsGeometryUtils::testAdjacentVertices()
+{
+  // test polygon - should wrap around!
+  QgsLineString* closedRing1 = new QgsLineString();
+  closedRing1->setPoints( QList<QgsPointV2>() << QgsPointV2( 1, 1 ) << QgsPointV2( 1, 2 ) << QgsPointV2( 2, 2 ) << QgsPointV2( 2, 1 ) << QgsPointV2( 1, 1 ) );
+  QgsPolygonV2 polygon1;
+  polygon1.setExteriorRing( closedRing1 );
+  QgsVertexId previous;
+  QgsVertexId next;
+
+  QgsGeometryUtils::adjacentVertices( polygon1, QgsVertexId( 0, 0, 0 ), previous, next );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 3 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 1 ) );
+}
+
 void TestQgsGeometryUtils::testDistanceToVertex()
 {
   //test with linestring
@@ -365,6 +382,84 @@ void TestQgsGeometryUtils::testDistanceToVertex()
   QCOMPARE( QgsGeometryUtils::distanceToVertex( point, QgsVertexId( 0, 0, 1 ) ), -1.0 );
 }
 
+void TestQgsGeometryUtils::testVerticesAtDistance()
+{
+  //test with linestring
+  QgsLineString* outerRing1 = new QgsLineString();
+  QgsVertexId previous;
+  QgsVertexId next;
+  QVERIFY( !QgsGeometryUtils::verticesAtDistance( *outerRing1, .5, previous, next ) );
+
+  outerRing1->setPoints( QList<QgsPointV2>() << QgsPointV2( 1, 1 ) << QgsPointV2( 1, 2 ) << QgsPointV2( 2, 2 ) << QgsPointV2( 2, 1 ) << QgsPointV2( 3, 1 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, .5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 0 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 1 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 1.5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 2 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 2.5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 2 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 3 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 3.5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 3 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 4 ) );
+  QVERIFY( ! QgsGeometryUtils::verticesAtDistance( *outerRing1, 4.5, previous, next ) );
+
+  // test exact hits
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 0, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 0 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 0 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 1, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 1 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 2, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 2 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 2 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 3, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 3 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 3 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *outerRing1, 4, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 4 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 4 ) );
+
+  // test closed line
+  QgsLineString* closedRing1 = new QgsLineString();
+  closedRing1->setPoints( QList<QgsPointV2>() << QgsPointV2( 1, 1 ) << QgsPointV2( 1, 2 ) << QgsPointV2( 2, 2 ) << QgsPointV2( 2, 1 ) << QgsPointV2( 1, 1 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *closedRing1, 0, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 0 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 0 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( *closedRing1, 4, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 4 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 4 ) );
+
+  // test with polygon
+  QgsPolygonV2 polygon1;
+  polygon1.setExteriorRing( closedRing1 );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( polygon1, .5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 0 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 1 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( polygon1, 1.5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 2 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( polygon1, 2.5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 2 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 3 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( polygon1, 3.5, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 3 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 4 ) );
+  QVERIFY( ! QgsGeometryUtils::verticesAtDistance( polygon1, 4.5, previous, next ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( polygon1, 0, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 0 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 0 ) );
+  QVERIFY( QgsGeometryUtils::verticesAtDistance( polygon1, 4, previous, next ) );
+  QCOMPARE( previous, QgsVertexId( 0, 0, 4 ) );
+  QCOMPARE( next, QgsVertexId( 0, 0, 4 ) );
+
+  //test with point
+  QgsPointV2 point( 1, 2 );
+  QVERIFY( !QgsGeometryUtils::verticesAtDistance( point, .5, previous, next ) );
+}
+
 void TestQgsGeometryUtils::testCircleCenterRadius_data()
 {
   QTest::addColumn<double>( "x1" );

tests/src/python/test_qgsgeometry.py

@@ -14,6 +14,7 @@
 
 import os
 import csv
+import math
 
 from qgis.core import (
     QgsGeometry,
@@ -3490,6 +3491,49 @@ def testLineLocatePoint(self):
         point = QgsGeometry.fromWkt('Point(9 -2)')
         self.assertAlmostEqual(geom.lineLocatePoint(point), 7.372, places=3)
 
+    def testInterpolateAngle(self):
+        """ test QgsGeometry.interpolateAngle() """
+
+        empty = QgsGeometry()
+        # just test no crash
+        self.assertEqual(empty.interpolateAngle(5), 0)
+
+        # not a linestring
+        point = QgsGeometry.fromWkt('Point(1 2)')
+        # no meaning, just test no crash!
+        self.assertEqual(point.interpolateAngle(5), 0)
+
+        # linestring
+        linestring = QgsGeometry.fromWkt('LineString(0 0, 10 0, 20 10, 20 20, 10 20)')
+        self.assertAlmostEqual(linestring.interpolateAngle(0), math.radians(90), places=3)
+        self.assertAlmostEqual(linestring.interpolateAngle(5), math.radians(90), places=3)
+        self.assertAlmostEqual(linestring.interpolateAngle(10), math.radians(67.5), places=3)
+        self.assertAlmostEqual(linestring.interpolateAngle(15), math.radians(45), places=3)
+        self.assertAlmostEqual(linestring.interpolateAngle(25), math.radians(0), places=3)
+        self.assertAlmostEqual(linestring.interpolateAngle(35), math.radians(270), places=3)
+
+        # test first and last points in a linestring - angle should be angle of
+        # first/last segment
+        linestring = QgsGeometry.fromWkt('LineString(20 0, 10 0, 10 -10)')
+        self.assertAlmostEqual(linestring.interpolateAngle(0), math.radians(270), places=3)
+        self.assertAlmostEqual(linestring.interpolateAngle(20), math.radians(180), places=3)
+
+        # polygon
+        polygon = QgsGeometry.fromWkt('Polygon((0 0, 10 0, 20 10, 20 20, 10 20, 0 0))')
+        self.assertAlmostEqual(polygon.interpolateAngle(5), math.radians(90), places=3)
+        self.assertAlmostEqual(polygon.interpolateAngle(10), math.radians(67.5), places=3)
+        self.assertAlmostEqual(polygon.interpolateAngle(15), math.radians(45), places=3)
+        self.assertAlmostEqual(polygon.interpolateAngle(25), math.radians(0), places=3)
+        self.assertAlmostEqual(polygon.interpolateAngle(35), math.radians(270), places=3)
+
+        # test first/last vertex in polygon
+        polygon = QgsGeometry.fromWkt('Polygon((0 0, 10 0, 10 10, 0 10, 0 0))')
+        self.assertAlmostEqual(polygon.interpolateAngle(0), math.radians(135), places=3)
+        self.assertAlmostEqual(polygon.interpolateAngle(40), math.radians(135), places=3)
+
+        # circular string
+        geom = QgsGeometry.fromWkt('CircularString (1 5, 6 2, 7 3)')
+        self.assertAlmostEqual(geom.interpolateAngle(5), 1.69120, places=3)
 
 if __name__ == '__main__':
     unittest.main()