[FEATURE] Add closest_point and shortest_line expression functions

nyalldawson · nyalldawson · commit 275eb9463d5b · 2016-01-15T12:50:17.000+11:00
closest_point: returns closest point a geometry to a second geometry
shortest_line: returns the shortest possible line joining two
geometries
diff --git a/python/core/geometry/qgsgeometry.sip b/python/core/geometry/qgsgeometry.sip
@@ -240,6 +240,18 @@ class QgsGeometry
      */
     double sqrDistToVertexAt( QgsPoint& point /In/, int atVertex ) const;
 
+    /** Returns the nearest point on this geometry to another geometry.
+     * @note added in QGIS 2.14
+     * @see shortestLine()
+     */
+    QgsGeometry nearestPoint( const QgsGeometry& other ) const;
+
+    /** Returns the shortest line joining this geometry to another geometry.
+     * @note added in QGIS 2.14
+     * @see nearestPoint()
+     */
+    QgsGeometry shortestLine( const QgsGeometry& other ) const;
+
     /**
      * Searches for the closest vertex in this geometry to the given point.
      * @param point Specifiest the point for search
diff --git a/resources/function_help/json/closest_point b/resources/function_help/json/closest_point
@@ -0,0 +1,15 @@
+{
+  "name": "closest_point",
+  "type": "function",
+  "description": "Returns the point on geometry 1 that is closest to geometry 2.",
+  "arguments": [
+    {"arg":"geometry 1","description":"geometry to find closest point on"},
+    {"arg":"geometry 2","description":"geometry to find closest point to"}
+  ],
+  "examples": [
+    {
+      "expression":"geom_to_wkt(closest_point(geom_from_wkt('LINESTRING (20 80, 98 190, 110 180, 50 75 )'),geom_from_wkt('POINT(100 100)')))",
+      "returns":"Point(73.0769 115.384)"
+    }
+  ]
+}
diff --git a/resources/function_help/json/shortest_line b/resources/function_help/json/shortest_line
@@ -0,0 +1,15 @@
+{
+  "name": "shortest_line",
+  "type": "function",
+  "description": "Returns the shortest line joining geometry 1 to geometry 2. The resultant line will start at geometry 1 and end at geometry 2.",
+  "arguments": [
+    {"arg":"geometry 1","description":"geometry to find shortest line from"},
+    {"arg":"geometry 2","description":"geometry to find shortest line to"}
+  ],
+  "examples": [
+    {
+      "expression":"geom_to_wkt(shortest_line(geom_from_wkt('LINESTRING (20 80, 98 190, 110 180, 50 75 )'),geom_from_wkt('POINT(100 100)')))",
+      "returns":"LineString(73.0769 115.384, 100 100)"
+    }
+  ]
+}
diff --git a/src/core/geometry/qgsgeometry.cpp b/src/core/geometry/qgsgeometry.cpp
@@ -502,6 +502,18 @@ double QgsGeometry::sqrDistToVertexAt( QgsPoint& point, int atVertex ) const
   return QgsGeometryUtils::sqrDistance2D( QgsPointV2( vertexPoint.x(), vertexPoint.y() ), QgsPointV2( point.x(), point.y() ) );
 }
 
+QgsGeometry QgsGeometry::nearestPoint( const QgsGeometry& other ) const
+{
+  QgsGeos geos( d->geometry );
+  return geos.closestPoint( other );
+}
+
+QgsGeometry QgsGeometry::shortestLine( const QgsGeometry& other ) const
+{
+  QgsGeos geos( d->geometry );
+  return geos.shortestLine( other );
+}
+
 double QgsGeometry::closestVertexWithContext( const QgsPoint& point, int& atVertex ) const
 {
   if ( !d->geometry )
diff --git a/src/core/geometry/qgsgeometry.h b/src/core/geometry/qgsgeometry.h
@@ -285,6 +285,18 @@ class CORE_EXPORT QgsGeometry
      */
     double sqrDistToVertexAt( QgsPoint& point, int atVertex ) const;
 
+    /** Returns the nearest point on this geometry to another geometry.
+     * @note added in QGIS 2.14
+     * @see shortestLine()
+     */
+    QgsGeometry nearestPoint( const QgsGeometry& other ) const;
+
+    /** Returns the shortest line joining this geometry to another geometry.
+     * @note added in QGIS 2.14
+     * @see nearestPoint()
+     */
+    QgsGeometry shortestLine( const QgsGeometry& other ) const;
+
     /**
      * Searches for the closest vertex in this geometry to the given point.
      * @param point Specifiest the point for search
diff --git a/src/core/geometry/qgsgeos.cpp b/src/core/geometry/qgsgeos.cpp
@@ -1764,6 +1764,84 @@ QgsAbstractGeometryV2* QgsGeos::reshapeGeometry( const QgsLineStringV2& reshapeW
   }
 }
 
+QgsGeometry QgsGeos::closestPoint( const QgsGeometry& other, QString* errorMsg ) const
+{
+  if ( !mGeos || other.isEmpty() )
+  {
+    return QgsGeometry();
+  }
+
+  GEOSGeomScopedPtr otherGeom( asGeos( other.geometry(), mPrecision ) );
+  if ( !otherGeom )
+  {
+    return QgsGeometry();
+  }
+
+  double nx = 0.0;
+  double ny = 0.0;
+  try
+  {
+    GEOSCoordSequence* nearestCoord = GEOSNearestPoints_r( geosinit.ctxt, mGeos, otherGeom.get() );
+
+    ( void )GEOSCoordSeq_getX_r( geosinit.ctxt, nearestCoord, 0, &nx );
+    ( void )GEOSCoordSeq_getY_r( geosinit.ctxt, nearestCoord, 0, &ny );
+    GEOSCoordSeq_destroy_r( geosinit.ctxt, nearestCoord );
+  }
+  catch ( GEOSException &e )
+  {
+    if ( errorMsg )
+    {
+      *errorMsg = e.what();
+    }
+    return QgsGeometry();
+  }
+
+  return QgsGeometry( new QgsPointV2( nx, ny ) );
+}
+
+QgsGeometry QgsGeos::shortestLine( const QgsGeometry& other, QString* errorMsg ) const
+{
+  if ( !mGeos || other.isEmpty() )
+  {
+    return QgsGeometry();
+  }
+
+  GEOSGeomScopedPtr otherGeom( asGeos( other.geometry(), mPrecision ) );
+  if ( !otherGeom )
+  {
+    return QgsGeometry();
+  }
+
+  double nx1 = 0.0;
+  double ny1 = 0.0;
+  double nx2 = 0.0;
+  double ny2 = 0.0;
+  try
+  {
+    GEOSCoordSequence* nearestCoord = GEOSNearestPoints_r( geosinit.ctxt, mGeos, otherGeom.get() );
+
+    ( void )GEOSCoordSeq_getX_r( geosinit.ctxt, nearestCoord, 0, &nx1 );
+    ( void )GEOSCoordSeq_getY_r( geosinit.ctxt, nearestCoord, 0, &ny1 );
+    ( void )GEOSCoordSeq_getX_r( geosinit.ctxt, nearestCoord, 1, &nx2 );
+    ( void )GEOSCoordSeq_getY_r( geosinit.ctxt, nearestCoord, 1, &ny2 );
+
+    GEOSCoordSeq_destroy_r( geosinit.ctxt, nearestCoord );
+  }
+  catch ( GEOSException &e )
+  {
+    if ( errorMsg )
+    {
+      *errorMsg = e.what();
+    }
+    return QgsGeometry();
+  }
+
+  QgsLineStringV2* line = new QgsLineStringV2();
+  line->addVertex( QgsPointV2( nx1, ny1 ) );
+  line->addVertex( QgsPointV2( nx2, ny2 ) );
+  return QgsGeometry( line );
+}
+
 GEOSGeometry* QgsGeos::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos , double precision )
 {
   if ( !line || !reshapeLineGeos )
diff --git a/src/core/geometry/qgsgeos.h b/src/core/geometry/qgsgeos.h
@@ -18,6 +18,7 @@ email                : marco.hugentobler at sourcepole dot com
 
 #include "qgsgeometryengine.h"
 #include "qgspointv2.h"
+#include "qgsgeometry.h"
 #include <geos_c.h>
 
 class QgsLineStringV2;
@@ -86,6 +87,18 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
     QgsAbstractGeometryV2* offsetCurve( double distance, int segments, int joinStyle, double mitreLimit, QString* errorMsg = nullptr ) const override;
     QgsAbstractGeometryV2* reshapeGeometry( const QgsLineStringV2& reshapeWithLine, int* errorCode, QString* errorMsg = nullptr ) const;
 
+    /** Returns the closest point on the geometry to the other geometry.
+     * @note added in QGIS 2.14
+     * @see shortestLine()
+     */
+    QgsGeometry closestPoint( const QgsGeometry& other, QString* errorMsg = nullptr ) const;
+
+    /** Returns the shortest line joining this geometry to the other geometry.
+     * @note added in QGIS 2.14
+     * @see closestPoint()
+     */
+    QgsGeometry shortestLine( const QgsGeometry& other, QString* errorMsg = nullptr ) const;
+
     /** Create a geometry from a GEOSGeometry
      * @param geos GEOSGeometry. Ownership is NOT transferred.
      */
diff --git a/src/core/qgsexpression.cpp b/src/core/qgsexpression.cpp
@@ -2168,6 +2168,28 @@ static QVariant fcnExtrude( const QVariantList& values, const QgsExpressionConte
   return result;
 }
 
+static QVariant fcnClosestPoint( const QVariantList& values, const QgsExpressionContext*, QgsExpression* parent )
+{
+  QgsGeometry fromGeom = getGeometry( values.at( 0 ), parent );
+  QgsGeometry toGeom = getGeometry( values.at( 1 ), parent );
+
+  QgsGeometry geom = fromGeom.nearestPoint( toGeom );
+
+  QVariant result = !geom.isEmpty() ? QVariant::fromValue( geom ) : QVariant();
+  return result;
+}
+
+static QVariant fcnShortestLine( const QVariantList& values, const QgsExpressionContext*, QgsExpression* parent )
+{
+  QgsGeometry fromGeom = getGeometry( values.at( 0 ), parent );
+  QgsGeometry toGeom = getGeometry( values.at( 1 ), parent );
+
+  QgsGeometry geom = fromGeom.shortestLine( toGeom );
+
+  QVariant result = !geom.isEmpty() ? QVariant::fromValue( geom ) : QVariant();
+  return result;
+}
+
 static QVariant fcnRound( const QVariantList& values, const QgsExpressionContext *, QgsExpression* parent )
 {
   if ( values.length() == 2 )
@@ -2765,6 +2787,7 @@ const QStringList& QgsExpression::BuiltinFunctions()
     << "overlaps" << "within" << "buffer" << "centroid" << "bounds" << "reverse" << "exterior_ring"
     << "bounds_width" << "bounds_height" << "is_closed" << "convex_hull" << "difference"
     << "distance" << "intersection" << "sym_difference" << "combine"
+    << "extrude" << "azimuth" << "closest_point" << "shortest_line"
     << "union" << "geom_to_wkt" << "geomToWKT" << "geometry"
     << "transform" << "get_feature" << "getFeature"
     << "levenshtein" << "longest_common_substring" << "hamming_distance"
@@ -2931,6 +2954,8 @@ const QList<QgsExpression::Function*>& QgsExpression::Functions()
     << new StaticFunction( "geometry", 1, fcnGetGeometry, "GeometryGroup", QString(), true )
     << new StaticFunction( "transform", 3, fcnTransformGeometry, "GeometryGroup" )
     << new StaticFunction( "extrude", 3, fcnExtrude, "GeometryGroup", QString() )
+    << new StaticFunction( "closest_point", 2, fcnClosestPoint, "GeometryGroup" )
+    << new StaticFunction( "shortest_line", 2, fcnShortestLine, "GeometryGroup" )
     << new StaticFunction( "$rownum", 0, fcnRowNumber, "deprecated" )
     << new StaticFunction( "$id", 0, fcnFeatureId, "Record" )
     << new StaticFunction( "$currentfeature", 0, fcnFeature, "Record" )
diff --git a/tests/src/core/testqgsexpression.cpp b/tests/src/core/testqgsexpression.cpp
@@ -583,6 +583,15 @@ class TestQgsExpression: public QObject
       QTest::newRow( "relate pattern false" ) << "relate( geom_from_wkt( 'LINESTRING(40 40,120 120)' ), geom_from_wkt( 'LINESTRING(40 40,60 120)' ), '**1F002**' )" << false << QVariant( false );
       QTest::newRow( "azimuth" ) << "toint(degrees(azimuth( make_point(25, 45), make_point(75, 100)))*1000000)" << false << QVariant( 42273689 );
       QTest::newRow( "azimuth" ) << "toint(degrees( azimuth( make_point(75, 100), make_point(25,45) ) )*1000000)" << false << QVariant( 222273689 );
+      QTest::newRow( "extrude geom" ) << "geom_to_wkt(extrude( geom_from_wkt('LineString( 1 2, 3 2, 4 3)'),1,2))" << false << QVariant( "Polygon ((1 2, 3 2, 4 3, 5 5, 4 4, 2 4, 1 2))" );
+      QTest::newRow( "extrude not geom" ) << "extrude('g',5,6)" << true << QVariant();
+      QTest::newRow( "extrude null" ) << "extrude(NULL,5,6)" << false << QVariant();
+      QTest::newRow( "closest_point geom" ) << "geom_to_wkt(closest_point( geom_from_wkt('LineString( 1 1, 5 1, 5 5 )'),geom_from_wkt('Point( 6 3 )')))" << false << QVariant( "Point (5 3)" );
+      QTest::newRow( "closest_point not geom" ) << "closest_point('g','b')" << true << QVariant();
+      QTest::newRow( "closest_point null" ) << "closest_point(NULL,NULL)" << false << QVariant();
+      QTest::newRow( "shortest_line geom" ) << "geom_to_wkt(shortest_line( geom_from_wkt('LineString( 1 1, 5 1, 5 5 )'),geom_from_wkt('Point( 6 3 )')))" << false << QVariant( "LineString (5 3, 6 3)" );
+      QTest::newRow( "shortest_line not geom" ) << "shortest_line('g','a')" << true << QVariant();
+      QTest::newRow( "shortest_line null" ) << "shortest_line(NULL,NULL)" << false << QVariant();
 
       // string functions
       QTest::newRow( "lower" ) << "lower('HeLLo')" << false << QVariant( "hello" );
diff --git a/tests/src/python/test_qgsgeometry.py b/tests/src/python/test_qgsgeometry.py
@@ -1192,6 +1192,10 @@ def testTranslate(self):
         wkt = polygon.exportToWkt()
 
     def testExtrude(self):
+        # test with empty geometry
+        g = QgsGeometry()
+        self.assertTrue(g.extrude(1, 2).isEmpty())
+
         points = [QgsPoint(1, 2), QgsPoint(3, 2), QgsPoint(4, 3)]
         line = QgsGeometry.fromPolyline(points)
         expected = QgsGeometry.fromWkt('Polygon ((1 2, 3 2, 4 3, 5 5, 4 4, 2 4, 1 2))')
@@ -1202,6 +1206,68 @@ def testExtrude(self):
         expected = QgsGeometry.fromWkt('MultiPolygon (((1 2, 3 2, 4 4, 2 4, 1 2)),((4 3, 8 3, 9 5, 5 5, 4 3)))')
         self.assertEqual(multiline.extrude(1, 2).exportToWkt(), expected.exportToWkt())
 
+    def testNearestPoint(self):
+        # test with empty geometries
+        g1 = QgsGeometry()
+        g2 = QgsGeometry()
+        self.assertTrue(g1.nearestPoint(g2).isEmpty())
+        g1 = QgsGeometry.fromWkt('LineString( 1 1, 5 1, 5 5 )')
+        self.assertTrue(g1.nearestPoint(g2).isEmpty())
+        self.assertTrue(g2.nearestPoint(g1).isEmpty())
+
+        g2 = QgsGeometry.fromWkt('Point( 6 3 )')
+        expWkt = 'Point( 5 3 )'
+        wkt = g1.nearestPoint(g2).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+        expWkt = 'Point( 6 3 )'
+        wkt = g2.nearestPoint(g1).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+
+        g1 = QgsGeometry.fromWkt('Polygon ((1 1, 5 1, 5 5, 1 5, 1 1))')
+        g2 = QgsGeometry.fromWkt('Point( 6 3 )')
+        expWkt = 'Point( 5 3 )'
+        wkt = g1.nearestPoint(g2).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+
+        expWkt = 'Point( 6 3 )'
+        wkt = g2.nearestPoint(g1).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+        g2 = QgsGeometry.fromWkt('Point( 2 3 )')
+        expWkt = 'Point( 2 3 )'
+        wkt = g1.nearestPoint(g2).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+
+    def testShortestLine(self):
+        # test with empty geometries
+        g1 = QgsGeometry()
+        g2 = QgsGeometry()
+        self.assertTrue(g1.shortestLine(g2).isEmpty())
+        g1 = QgsGeometry.fromWkt('LineString( 1 1, 5 1, 5 5 )')
+        self.assertTrue(g1.shortestLine(g2).isEmpty())
+        self.assertTrue(g2.shortestLine(g1).isEmpty())
+
+        g2 = QgsGeometry.fromWkt('Point( 6 3 )')
+        expWkt = 'LineString( 5 3, 6 3 )'
+        wkt = g1.shortestLine(g2).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+        expWkt = 'LineString( 6 3, 5 3 )'
+        wkt = g2.shortestLine(g1).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+
+        g1 = QgsGeometry.fromWkt('Polygon ((1 1, 5 1, 5 5, 1 5, 1 1))')
+        g2 = QgsGeometry.fromWkt('Point( 6 3 )')
+        expWkt = 'LineString( 5 3, 6 3 )'
+        wkt = g1.shortestLine(g2).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+
+        expWkt = 'LineString( 6 3, 5 3 )'
+        wkt = g2.shortestLine(g1).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+        g2 = QgsGeometry.fromWkt('Point( 2 3 )')
+        expWkt = 'LineString( 2 3, 2 3 )'
+        wkt = g1.shortestLine(g2).exportToWkt()
+        self.assertTrue(compareWkt(expWkt, wkt), "Expected:\n%s\nGot:\n%s\n" % (expWkt, wkt))
+
     def testBoundingBox(self):
         # 2-+-+-+-+-3
         # |         |
