[FEATURE] Expose GEOS linear referencing function to QgsGeometry

nyalldawson · nyalldawson · commit 193e9e746f84 · 2016-08-29T13:42:39.000+10:00
Adds a new QgsGeometry::lineLocatePoint() function for retrieving the distance along a linestring to the nearest position on the linestring to a given point. (cherry-picked from 409dfdf)
diff --git a/python/core/geometry/qgsgeometry.sip b/python/core/geometry/qgsgeometry.sip
@@ -476,6 +476,18 @@ class QgsGeometry
      */
     QgsGeometry* interpolate( double distance ) /Factory/;
 
+    /** 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 /Factory/;
 
diff --git a/src/core/geometry/qgsgeometry.cpp b/src/core/geometry/qgsgeometry.cpp
@@ -1391,6 +1391,24 @@ QgsGeometry* QgsGeometry::interpolate( double distance ) const
   return new QgsGeometry( result );
 }
 
+double QgsGeometry::lineLocatePoint( const QgsGeometry& point ) const
+{
+  if ( type() != QGis::Line )
+    return -1;
+
+  if ( QgsWKBTypes::flatType( point.d->geometry->wkbType() ) != QgsWKBTypes::Point )
+    return -1;
+
+  QgsGeometry segmentized = *this;
+  if ( QgsWKBTypes::isCurvedType( d->geometry->wkbType() ) )
+  {
+    segmentized = QgsGeometry( static_cast< QgsCurveV2* >( 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->d->geometry )
diff --git a/src/core/geometry/qgsgeometry.h b/src/core/geometry/qgsgeometry.h
@@ -515,9 +515,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;
 
diff --git a/src/core/geometry/qgsgeos.cpp b/src/core/geometry/qgsgeos.cpp
@@ -1891,6 +1891,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 )
diff --git a/src/core/geometry/qgsgeos.h b/src/core/geometry/qgsgeos.h
@@ -110,6 +110,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.
      */
diff --git a/tests/src/python/test_qgsgeometry.py b/tests/src/python/test_qgsgeometry.py
@@ -3390,5 +3390,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()
