Add QgsGeometry method to get bisector of angle at vertex

Sponsored by Andreas Neumann

python/core/geometry/qgsgeometry.sip

@@ -194,6 +194,15 @@ class QgsGeometry
      */
     double distanceToVertex( int vertex ) const;
 
+    /**
+     * Returns the bisector angle for this geometry at the specified vertex.
+     * @param vertex vertex index to calculate bisector angle at
+     * @returns bisector angle, in radians clockwise from north
+     * @note added in QGIS 3.0
+     * @see interpolateAngle()
+     */
+    double angleAtVertex( int vertex ) const;
+
     /**
      * Returns the indexes of the vertices before and after the given vertex index.
      *
@@ -551,6 +560,7 @@ class QgsGeometry
      * of the node is returned.
      * @param distance distance along geometry
      * @note added in QGIS 3.0
+     * @see angleAtVertex()
      */
     double interpolateAngle( double distance ) const;
 

src/core/geometry/qgsgeometry.cpp

@@ -404,6 +404,46 @@ double QgsGeometry::distanceToVertex( int vertex ) const
   return QgsGeometryUtils::distanceToVertex( *( d->geometry ), id );
 }
 
+double QgsGeometry::angleAtVertex( int vertex ) const
+{
+  if ( !d->geometry )
+  {
+    return 0;
+  }
+
+  QgsVertexId v2;
+  if ( !vertexIdFromVertexNr( vertex, v2 ) )
+  {
+    return 0;
+  }
+
+  QgsVertexId v1;
+  QgsVertexId v3;
+  QgsGeometryUtils::adjacentVertices( *d->geometry, v2, v1, v3 );
+  if ( v1.isValid() && v3.isValid() )
+  {
+    QgsPointV2 p1 = d->geometry->vertexAt( v1 );
+    QgsPointV2 p2 = d->geometry->vertexAt( v2 );
+    QgsPointV2 p3 = d->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 = d->geometry->vertexAt( v2 );
+    QgsPointV2 p2 = d->geometry->vertexAt( v3 );
+    return QgsGeometryUtils::lineAngle( p1.x(), p1.y(), p2.x(), p2.y() );
+  }
+  else if ( v1.isValid() )
+  {
+    QgsPointV2 p1 = d->geometry->vertexAt( v1 );
+    QgsPointV2 p2 = d->geometry->vertexAt( v2 );
+    return QgsGeometryUtils::lineAngle( p1.x(), p1.y(), p2.x(), p2.y() );
+  }
+  return 0.0;
+}
+
 void QgsGeometry::adjacentVertices( int atVertex, int& beforeVertex, int& afterVertex ) const
 {
   if ( !d->geometry )

src/core/geometry/qgsgeometry.h

@@ -237,6 +237,15 @@ class CORE_EXPORT QgsGeometry
      */
     double distanceToVertex( int vertex ) const;
 
+    /**
+     * Returns the bisector angle for this geometry at the specified vertex.
+     * @param vertex vertex index to calculate bisector angle at
+     * @returns bisector angle, in radians clockwise from north
+     * @note added in QGIS 3.0
+     * @see interpolateAngle()
+     */
+    double angleAtVertex( int vertex ) const;
+
     /**
      * Returns the indexes of the vertices before and after the given vertex index.
      *
@@ -588,6 +597,7 @@ class CORE_EXPORT QgsGeometry
      * of the node is returned.
      * @param distance distance along geometry
      * @note added in QGIS 3.0
+     * @see angleAtVertex()
      */
     double interpolateAngle( double distance ) const;
 

src/core/geometry/qgsgeometryutils.cpp

@@ -201,9 +201,16 @@ void QgsGeometryUtils::adjacentVertices( const QgsAbstractGeometry& geom, QgsVer
   }
   else if ( atVertex.vertex == 0 )
   {
-    afterVertex.part = atVertex.part;
-    afterVertex.ring = atVertex.ring;
-    afterVertex.vertex = atVertex.vertex + 1;
+    if ( ring.size() > 1 )
+    {
+      afterVertex.part = atVertex.part;
+      afterVertex.ring = atVertex.ring;
+      afterVertex.vertex = atVertex.vertex + 1;
+    }
+    else
+    {
+      afterVertex = QgsVertexId(); //after vertex invalid
+    }
     if ( polygonType && ring.size() > 3 )
     {
       beforeVertex.part = atVertex.part;

tests/src/core/testqgsgeometryutils.cpp

@@ -350,6 +350,12 @@ void TestQgsGeometryUtils::testAdjacentVertices()
   QgsGeometryUtils::adjacentVertices( polygon1, QgsVertexId( 0, 0, 0 ), previous, next );
   QCOMPARE( previous, QgsVertexId( 0, 0, 3 ) );
   QCOMPARE( next, QgsVertexId( 0, 0, 1 ) );
+
+  // test point - both vertices should be invalid
+  QgsPointV2 point( 1, 2 );
+  QgsGeometryUtils::adjacentVertices( point, QgsVertexId( 0, 0, 0 ), previous, next );
+  QVERIFY( !previous.isValid() );
+  QVERIFY( !next.isValid() );
 }
 
 void TestQgsGeometryUtils::testDistanceToVertex()

tests/src/python/test_qgsgeometry.py

@@ -3560,6 +3560,39 @@ def testInterpolate(self):
         self.assertTrue(compareWkt(result, exp, 0.00001),
                         "Interpolate: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
 
+    def testAngleAtVertex(self):
+        """ test QgsGeometry.angleAtVertex """
+
+        empty = QgsGeometry()
+        # just test no crash
+        self.assertEqual(empty.angleAtVertex(0), 0)
+
+        # not a linestring
+        point = QgsGeometry.fromWkt('Point(1 2)')
+        # no meaning, just test no crash!
+        self.assertEqual(point.angleAtVertex(0), 0)
+
+        # linestring
+        linestring = QgsGeometry.fromWkt('LineString(0 0, 10 0, 20 10, 20 20, 10 20)')
+        self.assertAlmostEqual(linestring.angleAtVertex(1), math.radians(67.5), places=3)
+        self.assertAlmostEqual(linestring.angleAtVertex(2), math.radians(22.5), places=3)
+        self.assertAlmostEqual(linestring.angleAtVertex(3), math.radians(315.0), places=3)
+        self.assertAlmostEqual(linestring.angleAtVertex(5), 0, places=3)
+        self.assertAlmostEqual(linestring.angleAtVertex(-1), 0, 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.angleAtVertex(0), math.radians(270), places=3)
+        self.assertAlmostEqual(linestring.angleAtVertex(2), math.radians(180), places=3)
+
+        # polygon
+        polygon = QgsGeometry.fromWkt('Polygon((0 0, 10 0, 10 10, 0 10, 0 0))')
+        self.assertAlmostEqual(polygon.angleAtVertex(0), math.radians(135.0), places=3)
+        self.assertAlmostEqual(polygon.angleAtVertex(1), math.radians(45.0), places=3)
+        self.assertAlmostEqual(polygon.angleAtVertex(2), math.radians(315.0), places=3)
+        self.assertAlmostEqual(polygon.angleAtVertex(3), math.radians(225.0), places=3)
+        self.assertAlmostEqual(polygon.angleAtVertex(4), math.radians(135.0), places=3)
 
 if __name__ == '__main__':
     unittest.main()