Skip to content
Permalink
Browse files

Fix incorrect reverse coordinate transforms are created when using a

PROJ 6+ build and a project has manual coordinate operation pathways
set

Refs #33121

(cherry picked from commit 51046d9)
  • Loading branch information
nyalldawson committed Dec 13, 2019
1 parent 8a7d6d2 commit 49c7652aa97974904197cb8629b247a73f0d50d9
@@ -657,7 +657,7 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
int projResult = 0;
#if PROJ_VERSION_MAJOR>=6
proj_errno_reset( projData );
proj_trans_generic( projData, direction == ForwardTransform ? PJ_FWD : PJ_INV,
proj_trans_generic( projData, ( direction == ForwardTransform && !d->mIsReversed ) || ( direction == ReverseTransform && d->mIsReversed ) ? PJ_FWD : PJ_INV,
x, sizeof( double ), numPoints,
y, sizeof( double ), numPoints,
z, sizeof( double ), numPoints,
@@ -813,6 +813,7 @@ void QgsCoordinateTransform::setCoordinateOperation( const QString &operation )
{
d.detach();
d->mProjCoordinateOperation = operation;
d->mShouldReverseCoordinateOperation = false;
}

const char *finder( const char *name )
@@ -106,6 +106,8 @@ QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinat
, mSourceDatumTransform( other.mSourceDatumTransform )
, mDestinationDatumTransform( other.mDestinationDatumTransform )
, mProjCoordinateOperation( other.mProjCoordinateOperation )
, mShouldReverseCoordinateOperation( other.mShouldReverseCoordinateOperation )
, mIsReversed( other.mIsReversed )
{
#if PROJ_VERSION_MAJOR < 6
//must reinitialize to setup mSourceProjection and mDestinationProjection
@@ -260,6 +262,7 @@ void QgsCoordinateTransformPrivate::calculateTransforms( const QgsCoordinateTran
// recalculate datum transforms from context
#if PROJ_VERSION_MAJOR >= 6
mProjCoordinateOperation = context.calculateCoordinateOperation( mSourceCRS, mDestCRS );
mShouldReverseCoordinateOperation = context.mustReverseCoordinateOperation( mSourceCRS, mDestCRS );
#else
Q_NOWARN_DEPRECATED_PUSH
QgsDatumTransform::TransformPair transforms = context.calculateDatumTransforms( mSourceCRS, mDestCRS );
@@ -331,6 +334,8 @@ ProjData QgsCoordinateTransformPrivate::threadLocalProjData()
QStringList projErrors;
proj_log_func( context, &projErrors, proj_collecting_logger );

mIsReversed = false;

QgsProjUtils::proj_pj_unique_ptr transform;
if ( !mProjCoordinateOperation.isEmpty() )
{
@@ -355,6 +360,10 @@ ProjData QgsCoordinateTransformPrivate::threadLocalProjData()

transform.reset();
}
else
{
mIsReversed = mShouldReverseCoordinateOperation;
}
}

QString nonAvailableError;
@@ -131,6 +131,10 @@ class QgsCoordinateTransformPrivate : public QSharedData
Q_DECL_DEPRECATED int mSourceDatumTransform = -1;
Q_DECL_DEPRECATED int mDestinationDatumTransform = -1;
QString mProjCoordinateOperation;
bool mShouldReverseCoordinateOperation = false;

//! True if the proj transform corresponds to the reverse direction, and must be flipped when transforming...
bool mIsReversed = false;

#if PROJ_VERSION_MAJOR<6

@@ -45,6 +45,7 @@ class TestQgsCoordinateTransform: public QObject
void transformContextNormalize();
void transformErrorMultiplePoints();
void transformErrorOnePoint();
void testDeprecated4240to4326();
};


@@ -533,5 +534,80 @@ void TestQgsCoordinateTransform::transformErrorOnePoint()
}
}

#include <proj.h>

void TestQgsCoordinateTransform::testDeprecated4240to4326()
{
#if PROJ_VERSION_MAJOR >= 6
// test creating a coordinate transform between EPSG 4240 and EPSG 4326 using a deprecated coordinate operation
// see https://github.com/qgis/QGIS/issues/33121

QgsCoordinateTransformContext context;
QgsCoordinateReferenceSystem src( QStringLiteral( "EPSG:4240" ) );
QgsCoordinateReferenceSystem dest( QStringLiteral( "EPSG:4326" ) );

// first use default transform
QgsCoordinateTransform defaultTransform( src, dest, context );
QCOMPARE( defaultTransform.coordinateOperation(), QString() );
QCOMPARE( defaultTransform.instantiatedCoordinateOperationDetails().proj, QStringLiteral( "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=evrst30 +step +proj=helmert +x=293 +y=836 +z=318 +rx=0.5 +ry=1.6 +rz=-2.8 +s=2.1 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg" ) );
QVERIFY( defaultTransform.isValid() );

const QgsPointXY p( 102.5, 7.5 );
QgsPointXY p2 = defaultTransform.transform( p );
QGSCOMPARENEAR( p2.x(), 102.494938, 0.000001 );
QGSCOMPARENEAR( p2.y(), 7.502624, 0.000001 );

QgsPointXY p3 = defaultTransform.transform( p2, QgsCoordinateTransform::ReverseTransform );
QGSCOMPARENEAR( p3.x(), 102.5, 0.000001 );
QGSCOMPARENEAR( p3.y(), 7.5, 0.000001 );

// and in reverse
QgsCoordinateTransform defaultTransformRev( dest, src, context );
QVERIFY( defaultTransformRev.isValid() );
QCOMPARE( defaultTransformRev.coordinateOperation(), QString() );
QgsDebugMsg( defaultTransformRev.instantiatedCoordinateOperationDetails().proj );
QCOMPARE( defaultTransformRev.instantiatedCoordinateOperationDetails().proj, QStringLiteral( "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=293 +y=836 +z=318 +rx=0.5 +ry=1.6 +rz=-2.8 +s=2.1 +convention=position_vector +step +inv +proj=cart +ellps=evrst30 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg" ) );

p2 = defaultTransformRev.transform( QgsPointXY( 102.494938, 7.502624 ) );
QGSCOMPARENEAR( p2.x(), 102.5, 0.000001 );
QGSCOMPARENEAR( p2.y(), 7.5, 0.000001 );

p3 = defaultTransformRev.transform( p2, QgsCoordinateTransform::ReverseTransform );
QGSCOMPARENEAR( p3.x(), 102.494938, 0.000001 );
QGSCOMPARENEAR( p3.y(), 7.502624, 0.000001 );

// now force use of deprecated transform
const QString deprecatedProj = QStringLiteral( "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=evrst30 +step +proj=helmert +x=209 +y=818 +z=290 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg" );
context.addCoordinateOperation( src, dest, deprecatedProj );

QgsCoordinateTransform deprecatedTransform( src, dest, context );
QVERIFY( deprecatedTransform.isValid() );
QCOMPARE( deprecatedTransform.coordinateOperation(), deprecatedProj );
QCOMPARE( deprecatedTransform.instantiatedCoordinateOperationDetails().proj, deprecatedProj );

p2 = deprecatedTransform.transform( p );
QGSCOMPARENEAR( p2.x(), 102.496547, 0.000001 );
QGSCOMPARENEAR( p2.y(), 7.502139, 0.000001 );

p3 = deprecatedTransform.transform( p2, QgsCoordinateTransform::ReverseTransform );
QGSCOMPARENEAR( p3.x(), 102.5, 0.000001 );
QGSCOMPARENEAR( p3.y(), 7.5, 0.000001 );

// and in reverse
QgsCoordinateTransform deprecatedTransformRev( dest, src, context );
QVERIFY( deprecatedTransformRev.isValid() );
QCOMPARE( deprecatedTransformRev.coordinateOperation(), deprecatedProj );
QCOMPARE( deprecatedTransformRev.instantiatedCoordinateOperationDetails().proj, deprecatedProj );

p2 = deprecatedTransformRev.transform( QgsPointXY( 102.496547, 7.502139 ) );
QGSCOMPARENEAR( p2.x(), 102.5, 0.000001 );
QGSCOMPARENEAR( p2.y(), 7.5, 0.000001 );

p3 = deprecatedTransformRev.transform( p2, QgsCoordinateTransform::ReverseTransform );
QGSCOMPARENEAR( p3.x(), 102.496547, 0.000001 );
QGSCOMPARENEAR( p3.y(), 7.502139, 0.000001 );
#endif
}

QGSTEST_MAIN( TestQgsCoordinateTransform )
#include "testqgscoordinatetransform.moc"

0 comments on commit 49c7652

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