Skip to content
Permalink
Browse files

[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.
  • Loading branch information
lbartoletti authored and nyalldawson committed Jan 28, 2017
1 parent 290758a commit 548461ac3cd366af1bd36a95f020fe21d7161336
@@ -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
@@ -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 );
}
@@ -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
@@ -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<const QgsPointV2*>( 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::Function*>& 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" ) )
@@ -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();
@@ -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()

0 comments on commit 548461a

Please sign in to comment.
You can’t perform that action at this time.