Fix incorrect reverse coordinate transforms are created when using a

PROJ 6+ build and a project has manual coordinate operation pathways set Refs #33121
qgis · Dec 13, 2019 · 51046d9 · 51046d9
1 parent a1ac778
commit 51046d9
Show file tree

Hide file tree

Showing 4 changed files with 91 additions and 1 deletion.
diff --git a/src/core/qgscoordinatetransform.cpp b/src/core/qgscoordinatetransform.cpp
@@ -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 )

diff --git a/src/core/qgscoordinatetransform_p.cpp b/src/core/qgscoordinatetransform_p.cpp
@@ -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;

diff --git a/src/core/qgscoordinatetransform_p.h b/src/core/qgscoordinatetransform_p.h
@@ -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
 

diff --git a/tests/src/core/testqgscoordinatetransform.cpp b/tests/src/core/testqgscoordinatetransform.cpp
@@ -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"