[FEATURE] Expose GEOS linear referencing function to QgsGeometry

Adds a new QgsGeometry::lineLocatePoint() function for
retrieving the distance along a linestring to the nearest
position on the linestring to a given point.

Author: nyalldawson

python/core/geometry/qgsgeometry.sip

@@ -527,8 +527,21 @@ class QgsGeometry
     /**
      * Return interpolated point on line at distance
      * @note added in 1.9
+     * @see lineLocatePoint()
+     */
+    QgsGeometry interpolate( double distance ) const;
+
+    /** Returns a distance representing the location along this linestring of the closest point
+     * on this linestring geometry to the specified point. Ie, the returned value indicates
+     * how far along this linestring you need to traverse to get to the closest location
+     * where this linestring comes to the specified point.
+     * @param point point to seek proximity to
+     * @return distance along line, or -1 on error
+     * @note only valid for linestring geometries
+     * @see interpolate()
+     * @note added in QGIS 3.0
      */
-    QgsGeometry interpolate( double distance );
+    double lineLocatePoint( const QgsGeometry& point ) const;
 
     /** Returns a geometry representing the points shared by this geometry and other. */
     QgsGeometry intersection( const QgsGeometry& geometry ) const;

src/core/geometry/qgsgeometry.cpp

@@ -1471,6 +1471,24 @@ QgsGeometry QgsGeometry::interpolate( double distance ) const
   return QgsGeometry( result );
 }
 
+double QgsGeometry::lineLocatePoint( const QgsGeometry& point ) const
+{
+  if ( type() != QgsWkbTypes::LineGeometry )
+    return -1;
+
+  if ( QgsWkbTypes::flatType( point.wkbType() ) != QgsWkbTypes::Point )
+    return -1;
+
+  QgsGeometry segmentized = *this;
+  if ( QgsWkbTypes::isCurvedType( wkbType() ) )
+  {
+    segmentized = QgsGeometry( static_cast< QgsCurve* >( d->geometry )->segmentize() );
+  }
+
+  QgsGeos geos( d->geometry );
+  return geos.lineLocatePoint( *( static_cast< QgsPointV2* >( point.d->geometry ) ) );
+}
+
 QgsGeometry QgsGeometry::intersection( const QgsGeometry& geometry ) const
 {
   if ( !d->geometry || geometry.isEmpty() )

src/core/geometry/qgsgeometry.h

@@ -564,9 +564,22 @@ class CORE_EXPORT QgsGeometry
     /**
      * Return interpolated point on line at distance
      * @note added in 1.9
+     * @see lineLocatePoint()
      */
     QgsGeometry interpolate( double distance ) const;
 
+    /** Returns a distance representing the location along this linestring of the closest point
+     * on this linestring geometry to the specified point. Ie, the returned value indicates
+     * how far along this linestring you need to traverse to get to the closest location
+     * where this linestring comes to the specified point.
+     * @param point point to seek proximity to
+     * @return distance along line, or -1 on error
+     * @note only valid for linestring geometries
+     * @see interpolate()
+     * @note added in QGIS 3.0
+     */
+    double lineLocatePoint( const QgsGeometry& point ) const;
+
     /** Returns a geometry representing the points shared by this geometry and other. */
     QgsGeometry intersection( const QgsGeometry& geometry ) const;
 

src/core/geometry/qgsgeos.cpp

@@ -1911,6 +1911,36 @@ QgsGeometry QgsGeos::shortestLine( const QgsGeometry& other, QString* errorMsg )
   return QgsGeometry( line );
 }
 
+double QgsGeos::lineLocatePoint( const QgsPointV2& point, QString* errorMsg ) const
+{
+  if ( !mGeos )
+  {
+    return -1;
+  }
+
+  GEOSGeomScopedPtr otherGeom( asGeos( &point, mPrecision ) );
+  if ( !otherGeom )
+  {
+    return -1;
+  }
+
+  double distance = -1;
+  try
+  {
+    distance = GEOSProject_r( geosinit.ctxt, mGeos, otherGeom.get() );
+  }
+  catch ( GEOSException &e )
+  {
+    if ( errorMsg )
+    {
+      *errorMsg = e.what();
+    }
+    return -1;
+  }
+
+  return distance;
+}
+
 
 /** Extract coordinates of linestring's endpoints. Returns false on error. */
 static bool _linestringEndpoints( const GEOSGeometry* linestring, double& x1, double& y1, double& x2, double& y2 )

src/core/geometry/qgsgeos.h

@@ -128,6 +128,17 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
      */
     QgsGeometry shortestLine( const QgsGeometry& other, QString* errorMsg = nullptr ) const;
 
+    /** Returns a distance representing the location along this linestring of the closest point
+     * on this linestring geometry to the specified point. Ie, the returned value indicates
+     * how far along this linestring you need to traverse to get to the closest location
+     * where this linestring comes to the specified point.
+     * @param point point to seek proximity to
+     * @param errorMsg error messages emitted, if any
+     * @note only valid for linestring geometries
+     * @return distance along line, or -1 on error
+     */
+    double lineLocatePoint( const QgsPointV2& point, QString* errorMsg = nullptr ) const;
+
     /** Create a geometry from a GEOSGeometry
      * @param geos GEOSGeometry. Ownership is NOT transferred.
      */

tests/src/python/test_qgsgeometry.py

@@ -3430,5 +3430,27 @@ def testMergeLines(self):
         exp = 'MultiLineString((0 0, 10 10),(12 2, 14 4))'
         self.assertTrue(compareWkt(result, exp, 0.00001), "Merge lines: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
 
+    def testLineLocatePoint(self):
+        """ test QgsGeometry.lineLocatePoint() """
+
+        # not a linestring
+        point = QgsGeometry.fromWkt('Point(1 2)')
+        self.assertEqual(point.lineLocatePoint(point), -1)
+
+        # not a point
+        linestring = QgsGeometry.fromWkt('LineString(0 0, 10 0)')
+        self.assertEqual(linestring.lineLocatePoint(linestring), -1)
+
+        # valid
+        self.assertEqual(linestring.lineLocatePoint(point), 1)
+        point = QgsGeometry.fromWkt('Point(9 -2)')
+        self.assertEqual(linestring.lineLocatePoint(point), 9)
+
+        # circular string
+        geom = QgsGeometry.fromWkt('CircularString (1 5, 6 2, 7 3)')
+        point = QgsGeometry.fromWkt('Point(9 -2)')
+        self.assertAlmostEqual(geom.lineLocatePoint(point), 7.372, places=3)
+
+
 if __name__ == '__main__':
     unittest.main()