From 548461ac3cd366af1bd36a95f020fe21d7161336 Mon Sep 17 00:00:00 2001 From: lbartoletti Date: Sat, 28 Jan 2017 23:29:25 +0100 Subject: [PATCH] [FEATURE] QgsPointV2 add project with 3D support Allows projection of a point with inclination to return a 3d point. Expression "project" function has been extended to support a new inclination parameter. --- python/core/geometry/qgspointv2.sip | 29 ++++++++++++++++ src/core/geometry/qgspointv2.cpp | 32 +++++++++++++++++ src/core/geometry/qgspointv2.h | 29 ++++++++++++++++ src/core/qgsexpression.cpp | 11 +++--- tests/src/core/testqgsexpression.cpp | 5 ++- tests/src/core/testqgsgeometry.cpp | 51 ++++++++++++++++++++++++++++ 6 files changed, 151 insertions(+), 6 deletions(-) diff --git a/python/core/geometry/qgspointv2.sip b/python/core/geometry/qgspointv2.sip index ec86363b4c8..f7832bc7f0b 100644 --- a/python/core/geometry/qgspointv2.sip +++ b/python/core/geometry/qgspointv2.sip @@ -179,6 +179,35 @@ class QgsPointV2: public QgsAbstractGeometry */ double azimuth( const QgsPointV2& other ) const; + /** Returns a new point which correspond to this point projected by a specified distance + * with specified angles (azimuth and inclination). + * M value is preserved. + * @param distance distance to project + * @param azimuth angle to project in X Y, clockwise in degrees starting from north + * @param inclination angle to project in Z (3D) + * @return The point projected. If a 2D point is projected a 3D point will be returned except if + * inclination is 90. A 3D point is always returned if a 3D point is projected. + * Example: + * \code{.py} + * p = QgsPointV2( 1, 2 ) # 2D point + * pr = p.project ( 1, 0 ) + * # pr is a 2D point: 'Point (1 3)' + * pr = p.project ( 1, 0, 90 ) + * # pr is a 2D point: 'Point (1 3)' + * pr = p.project (1, 0, 0 ) + * # pr is a 3D point: 'PointZ (1 2 1)' + * p = QgsPointV2( QgsWkbTypes.PointZ, 1, 2, 2 ) # 3D point + * pr = p.project ( 1, 0 ) + * # pr is a 3D point: 'PointZ (1 3 2)' + * pr = p.project ( 1, 0, 90 ) + * # pr is a 3D point: 'PointZ (1 3 2)' + * pr = p.project (1, 0, 0 ) + * # pr is a 3D point: 'PointZ (1 2 3)' + * \endcode + * @note added in QGIS 3.0 + */ + QgsPointV2 project( double distance, double azimuth, double inclination = 90.0 ) const; + /** * Calculates the vector obtained by subtracting a point from this point. * @note added in QGIS 3.0 diff --git a/src/core/geometry/qgspointv2.cpp b/src/core/geometry/qgspointv2.cpp index 15cb5194fcb..266935f13d1 100644 --- a/src/core/geometry/qgspointv2.cpp +++ b/src/core/geometry/qgspointv2.cpp @@ -462,3 +462,35 @@ double QgsPointV2::azimuth( const QgsPointV2& other ) const double dy = other.y() - mY; return ( atan2( dx, dy ) * 180.0 / M_PI ); } + + +QgsPointV2 QgsPointV2::project( double distance, double azimuth, double inclination ) const +{ + QgsWkbTypes::Type pType(QgsWkbTypes::Point); + + double rads_xy = azimuth * M_PI / 180.0; + double dx = 0.0, dy = 0.0, dz = 0.0; + + inclination = fmod( inclination, 360.0 ); + + if ( !is3D() && qgsDoubleNear(inclination, 90.0) ) + { + dx = distance * sin( rads_xy ); + dy = distance * cos( rads_xy ); + } + else + { + pType = QgsWkbTypes::addZ( pType ); + double rads_z = inclination * M_PI / 180.0; + dx = distance * sin( rads_z ) * sin( rads_xy ); + dy = distance * sin( rads_z ) * cos( rads_xy ); + dz = distance * cos( rads_z ); + } + + if ( isMeasure() ) + { + pType = QgsWkbTypes::addM( pType ); + } + + return QgsPointV2( pType, mX + dx, mY + dy, mZ + dz, mM ); +} diff --git a/src/core/geometry/qgspointv2.h b/src/core/geometry/qgspointv2.h index 34e17bd1aa8..c6635948ab3 100644 --- a/src/core/geometry/qgspointv2.h +++ b/src/core/geometry/qgspointv2.h @@ -193,6 +193,35 @@ class CORE_EXPORT QgsPointV2: public QgsAbstractGeometry */ double azimuth( const QgsPointV2& other ) const; + /** Returns a new point which correspond to this point projected by a specified distance + * with specified angles (azimuth and inclination). + * M value is preserved. + * @param distance distance to project + * @param azimuth angle to project in X Y, clockwise in degrees starting from north + * @param inclination angle to project in Z (3D) + * @return The point projected. If a 2D point is projected a 3D point will be returned except if + * inclination is 90. A 3D point is always returned if a 3D point is projected. + * Example: + * \code{.py} + * p = QgsPointV2( 1, 2 ) # 2D point + * pr = p.project ( 1, 0 ) + * # pr is a 2D point: 'Point (1 3)' + * pr = p.project ( 1, 0, 90 ) + * # pr is a 2D point: 'Point (1 3)' + * pr = p.project (1, 0, 0 ) + * # pr is a 3D point: 'PointZ (1 2 1)' + * p = QgsPointV2( QgsWkbTypes.PointZ, 1, 2, 2 ) # 3D point + * pr = p.project ( 1, 0 ) + * # pr is a 3D point: 'PointZ (1 3 2)' + * pr = p.project ( 1, 0, 90 ) + * # pr is a 3D point: 'PointZ (1 3 2)' + * pr = p.project (1, 0, 0 ) + * # pr is a 3D point: 'PointZ (1 2 3)' + * \endcode + * @note added in QGIS 3.0 + */ + QgsPointV2 project(double distance, double azimuth, double inclination = 90.0 ) const; + /** * Calculates the vector obtained by subtracting a point from this point. * @note added in QGIS 3.0 diff --git a/src/core/qgsexpression.cpp b/src/core/qgsexpression.cpp index 9e30bc07813..0fff5ea5d7f 100644 --- a/src/core/qgsexpression.cpp +++ b/src/core/qgsexpression.cpp @@ -2762,12 +2762,13 @@ static QVariant fcnProject( const QVariantList& values, const QgsExpressionConte } double distance = getDoubleValue( values.at( 1 ), parent ); - double bearing = getDoubleValue( values.at( 2 ), parent ); + double azimuth = getDoubleValue( values.at( 2 ), parent ); + double inclination = getDoubleValue( values.at( 3 ), parent ); - QgsPoint p = geom.asPoint(); - QgsPoint newPoint = p.project( distance, ( 180 * bearing ) / M_PI ); + const QgsPointV2* p = dynamic_cast( geom.geometry() ); + QgsPointV2 newPoint = p->project( distance, 180.0 * azimuth / M_PI, 180.0 * inclination / M_PI ); - return QVariant::fromValue( QgsGeometry( new QgsPointV2( newPoint.x(), newPoint.y() ) ) ); + return QVariant::fromValue( QgsGeometry( new QgsPointV2( newPoint) ) ); } static QVariant fcnExtrude( const QVariantList& values, const QgsExpressionContext*, QgsExpression* parent ) @@ -3765,7 +3766,7 @@ const QList& QgsExpression::Functions() << new StaticFunction( QStringLiteral( "radians" ), ParameterList() << Parameter( QStringLiteral( "degrees" ) ), fcnRadians, QStringLiteral( "Math" ) ) << new StaticFunction( QStringLiteral( "degrees" ), ParameterList() << Parameter( QStringLiteral( "radians" ) ), fcnDegrees, QStringLiteral( "Math" ) ) << new StaticFunction( QStringLiteral( "azimuth" ), ParameterList() << Parameter( QStringLiteral( "point_a" ) ) << Parameter( QStringLiteral( "point_b" ) ), fcnAzimuth, QStringList() << QStringLiteral( "Math" ) << QStringLiteral( "GeometryGroup" ) ) - << new StaticFunction( QStringLiteral( "project" ), ParameterList() << Parameter( QStringLiteral( "point" ) ) << Parameter( QStringLiteral( "distance" ) ) << Parameter( QStringLiteral( "bearing" ) ), fcnProject, QStringLiteral( "GeometryGroup" ) ) + << new StaticFunction( QStringLiteral( "project" ), ParameterList() << Parameter( QStringLiteral( "point" ) ) << Parameter( QStringLiteral( "distance" ) ) << Parameter( QStringLiteral( "azimuth" ) ) << Parameter( QStringLiteral( "elevation" ), true, M_PI / 2 ), fcnProject, QStringLiteral( "GeometryGroup" ) ) << new StaticFunction( QStringLiteral( "abs" ), ParameterList() << Parameter( QStringLiteral( "value" ) ), fcnAbs, QStringLiteral( "Math" ) ) << new StaticFunction( QStringLiteral( "cos" ), ParameterList() << Parameter( QStringLiteral( "angle" ) ), fcnCos, QStringLiteral( "Math" ) ) << new StaticFunction( QStringLiteral( "sin" ), ParameterList() << Parameter( QStringLiteral( "angle" ) ), fcnSin, QStringLiteral( "Math" ) ) diff --git a/tests/src/core/testqgsexpression.cpp b/tests/src/core/testqgsexpression.cpp index da7c2982cd7..2ceafbbfa1c 100644 --- a/tests/src/core/testqgsexpression.cpp +++ b/tests/src/core/testqgsexpression.cpp @@ -786,7 +786,10 @@ class TestQgsExpression: public QObject QTest::newRow( "project not geom" ) << "project( 'asd', 1, 2 )" << true << QVariant(); QTest::newRow( "project not point" ) << "project( geom_from_wkt('LINESTRING(2 0,2 2, 3 2, 3 0)'), 1, 2 )" << true << QVariant(); QTest::newRow( "project x" ) << "toint(x(project( make_point( 1, 2 ), 3, radians(270)))*1000000)" << false << QVariant( -2 * 1000000 ); - QTest::newRow( "project y" ) << "toint(y(project( point:=make_point( 1, 2 ), distance:=3, bearing:=radians(270)))*1000000)" << false << QVariant( 2 * 1000000 ); + QTest::newRow( "project y" ) << "toint(y(project( point:=make_point( 1, 2 ), distance:=3, azimuth:=radians(270)))*1000000)" << false << QVariant( 2 * 1000000 ); + QTest::newRow( "project m value preserved" ) << "geom_to_wkt(project( make_point( 1, 2, 2, 5), 1, 0.0, 0.0 ) )" << false << QVariant( "PointZM (1 2 3 5)" ); + QTest::newRow( "project 2D Point" ) << "geom_to_wkt(project( point:=make_point( 1, 2), distance:=1, azimuth:=radians(0), elevation:=0 ) )" << false << QVariant( "PointZ (1 2 1)" ); + QTest::newRow( "project 3D Point" ) << "geom_to_wkt(project( make_point( 1, 2, 2), 5, radians(450), radians (450) ) )" << false << QVariant( "PointZ (6 2 2)" ); 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(); diff --git a/tests/src/core/testqgsgeometry.cpp b/tests/src/core/testqgsgeometry.cpp index 22de32053b4..54b46290002 100644 --- a/tests/src/core/testqgsgeometry.cpp +++ b/tests/src/core/testqgsgeometry.cpp @@ -839,6 +839,57 @@ void TestQgsGeometry::point() QCOMPARE( p31 - QgsVector( 3, 5 ), QgsPointV2( 1 , 2 ) ); p31 -= QgsVector( 3, 5 ); QCOMPARE( p31, QgsPointV2( 1, 2 ) ); + + // test projecting a point + // 2D + QgsPointV2 p33 = QgsPointV2( 1, 2 ); + QCOMPARE( p33.project( 1, 0 ), QgsPointV2( 1, 3 ) ); + QCOMPARE( p33.project( 1, 0, 0 ), QgsPointV2( QgsWkbTypes::PointZ, 1, 2, 1 ) ); + QCOMPARE( p33.project( 1.5, 90 ), QgsPointV2( 2.5, 2 ) ); + QCOMPARE( p33.project( 1.5, 90, 90 ), QgsPointV2( 2.5, 2 ) ); // stay QgsWkbTypes::Point + QCOMPARE( p33.project( 2, 180 ), QgsPointV2( 1, 0 ) ); + QCOMPARE( p33.project( 2, 180, 180 ), QgsPointV2( QgsWkbTypes::PointZ, 1, 2, -2 ) ); + QCOMPARE( p33.project( 5, 270 ), QgsPointV2( -4, 2 ) ); + QCOMPARE( p33.project( 5, 270, 270 ), QgsPointV2( QgsWkbTypes::PointZ, 6, 2, 0 ) ); + QCOMPARE( p33.project( 6, 360 ), QgsPointV2( 1, 8 ) ); + QCOMPARE( p33.project( 6, 360, 360 ), QgsPointV2( QgsWkbTypes::PointZ, 1, 2, 6 ) ); + QCOMPARE( p33.project( 5, 450 ), QgsPointV2( 6, 2 ) ); + QCOMPARE( p33.project( 5, 450, 450 ), QgsPointV2( 6, 2 ) ); // stay QgsWkbTypes::Point + QCOMPARE( p33.project( -1, 0 ), QgsPointV2( 1, 1 ) ); + QCOMPARE( p33.project( -1, 0, 0 ), QgsPointV2( QgsWkbTypes::PointZ, 1, 2, -1 ) ); + QCOMPARE( p33.project( 1.5, -90 ), QgsPointV2( -0.5, 2 ) ); + QCOMPARE( p33.project( 1.5, -90, -90 ), QgsPointV2( QgsWkbTypes::PointZ, 2.5, 2, 0 ) ); + // PointM + p33.addMValue( 5.0 ); + QCOMPARE( p33.project( 1, 0 ), QgsPointV2( QgsWkbTypes::PointM, 1, 3, 0, 5 ) ); + QCOMPARE( p33.project( 1, 0, 0 ), QgsPointV2( QgsWkbTypes::PointZM, 1, 2, 1, 5 ) ); + QCOMPARE( p33.project( 5, 450, 450 ), QgsPointV2( QgsWkbTypes::PointM, 6, 2, 0, 5 ) ); + + // 3D + QgsPointV2 p34 = QgsPointV2( QgsWkbTypes::PointZ, 1, 2, 2 ); + QCOMPARE( p34.project( 1, 0 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 3, 2 ) ); + QCOMPARE( p34.project( 1, 0, 0 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 2, 3 ) ); + QCOMPARE( p34.project( 1.5, 90 ), QgsPointV2(QgsWkbTypes::PointZ, 2.5, 2, 2 ) ); + QCOMPARE( p34.project( 1.5, 90, 90 ), QgsPointV2(QgsWkbTypes::PointZ, 2.5, 2, 2 ) ); + QCOMPARE( p34.project( 2, 180 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 0, 2 ) ); + QCOMPARE( p34.project( 2, 180, 180 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 2, 0 ) ); + QCOMPARE( p34.project( 5, 270 ), QgsPointV2(QgsWkbTypes::PointZ, -4, 2, 2 ) ); + QCOMPARE( p34.project( 5, 270, 270 ), QgsPointV2(QgsWkbTypes::PointZ, 6, 2, 2 ) ); + QCOMPARE( p34.project( 6, 360 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 8, 2 ) ); + QCOMPARE( p34.project( 6, 360, 360 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 2, 8 ) ); + QCOMPARE( p34.project( 5, 450 ), QgsPointV2(QgsWkbTypes::PointZ, 6, 2, 2 ) ); + QCOMPARE( p34.project( 5, 450, 450 ), QgsPointV2(QgsWkbTypes::PointZ, 6, 2, 2 ) ); + QCOMPARE( p34.project( -1, 0 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 1, 2 ) ); + QCOMPARE( p34.project( -1, 0, 0 ), QgsPointV2(QgsWkbTypes::PointZ, 1, 2, 1 ) ); + QCOMPARE( p34.project( 1.5, -90 ), QgsPointV2(QgsWkbTypes::PointZ, -0.5, 2, 2 ) ); + QCOMPARE( p34.project( 1.5, -90, -90 ), QgsPointV2(QgsWkbTypes::PointZ, 2.5, 2, 2 ) ); + // PointM + p34.addMValue( 5.0 ); + QCOMPARE( p34.project( 1, 0 ), QgsPointV2(QgsWkbTypes::PointZM, 1, 3, 2, 5 ) ); + QCOMPARE( p34.project( 1, 0, 0 ), QgsPointV2(QgsWkbTypes::PointZM, 1, 2, 3, 5 ) ); + QCOMPARE( p34.project( 5, 450 ), QgsPointV2(QgsWkbTypes::PointZM, 6, 2, 2, 5 ) ); + QCOMPARE( p34.project( 5, 450, 450 ), QgsPointV2(QgsWkbTypes::PointZM, 6, 2, 2, 5 ) ); + } void TestQgsGeometry::lineString()