Correctly handle swapped axis order for proj 6.0

qgis · Apr 5, 2019 · 83dc159 · 83dc159
1 parent 905ccc0
commit 83dc159
Show file tree

Hide file tree

Showing 7 changed files with 125 additions and 2 deletions.
diff --git a/src/core/qgscoordinatetransform.cpp b/src/core/qgscoordinatetransform.cpp
@@ -648,12 +648,24 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
 
   int projResult = 0;
 #if PROJ_VERSION_MAJOR>=6
+  const bool sourceAxisOrderSwapped = d->mSourceAxisOrderSwapped;
+  const bool destinationAxisOrderSwapped = d->mDestAxisOrderSwapped;
+
   proj_trans_generic( projData, direction == ForwardTransform ? PJ_FWD : PJ_INV,
-                      x, sizeof( double ), numPoints,
-                      y, sizeof( double ), numPoints,
+                      !sourceAxisOrderSwapped ? x : y, sizeof( double ), numPoints,
+                      !sourceAxisOrderSwapped ? y : x, sizeof( double ), numPoints,
                       z, sizeof( double ), numPoints,
                       nullptr, sizeof( double ), 0 );
   projResult = proj_errno( projData );
+  if ( projResult == 0 && destinationAxisOrderSwapped )
+  {
+    size_t size = sizeof( double ) * numPoints;
+    void *tmp = malloc( size );
+    memcpy( tmp, x, size );
+    memcpy( x, y, size );
+    memcpy( y, tmp, size );
+    free( tmp );
+  }
 #else
   if ( direction == ReverseTransform )
   {

diff --git a/src/core/qgscoordinatetransform_p.cpp b/src/core/qgscoordinatetransform_p.cpp
@@ -84,6 +84,8 @@ QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinat
   , mDestCRS( other.mDestCRS )
   , mSourceDatumTransform( other.mSourceDatumTransform )
   , mDestinationDatumTransform( other.mDestinationDatumTransform )
+  , mSourceAxisOrderSwapped( other.mSourceAxisOrderSwapped )
+  , mDestAxisOrderSwapped( other.mDestAxisOrderSwapped )
 {
   //must reinitialize to setup mSourceProjection and mDestinationProjection
   initialize();
@@ -168,6 +170,16 @@ bool QgsCoordinateTransformPrivate::initialize()
   // create proj projections for current thread
   ProjData res = threadLocalProjData();
 
+#if PROJ_VERSION_MAJOR>=6
+  PJ_CONTEXT *context = QgsProjContext::get();
+  QgsProjUtils::proj_pj_unique_ptr sourceCrs( proj_get_source_crs( context, res ) );
+  if ( sourceCrs )
+    mSourceAxisOrderSwapped = QgsProjUtils::axisOrderIsSwapped( sourceCrs.get() );
+  QgsProjUtils::proj_pj_unique_ptr destCrs( proj_get_target_crs( context, res ) );
+  if ( destCrs )
+    mDestAxisOrderSwapped = QgsProjUtils::axisOrderIsSwapped( destCrs.get() );
+#endif
+
 #ifdef COORDINATE_TRANSFORM_VERBOSE
   QgsDebugMsg( "From proj : " + mSourceCRS.toProj4() );
   QgsDebugMsg( "To proj   : " + mDestCRS.toProj4() );

diff --git a/src/core/qgscoordinatetransform_p.h b/src/core/qgscoordinatetransform_p.h
@@ -126,6 +126,9 @@ class QgsCoordinateTransformPrivate : public QSharedData
     int mSourceDatumTransform = -1;
     int mDestinationDatumTransform = -1;
 
+    bool mSourceAxisOrderSwapped = false;
+    bool mDestAxisOrderSwapped = false;
+
 #if PROJ_VERSION_MAJOR<6
 
     /**

diff --git a/src/core/qgsprojutils.cpp b/src/core/qgsprojutils.cpp
@@ -112,4 +112,36 @@ bool QgsProjUtils::usesAngularUnit( const QString &projDef )
   return false;
 }
 
+bool QgsProjUtils::axisOrderIsSwapped( const PJ *crs )
+{
+  //ported from https://github.com/pramsey/postgis/blob/7ecf6839c57a838e2c8540001a3cd35b78a730db/liblwgeom/lwgeom_transform.c#L299
+  if ( !crs )
+    return false;
+
+  PJ_CONTEXT *context = QgsProjContext::get();
+  QgsProjUtils::proj_pj_unique_ptr pjCrs( proj_crs_get_coordinate_system( context, crs ) );
+  if ( !pjCrs )
+    return false;
+
+  const int axisCount = proj_cs_get_axis_count( context, pjCrs.get() );
+  if ( axisCount > 0 )
+  {
+    const char *outDirection = nullptr;
+    // Read only first axis, see if it is degrees / north
+
+    proj_cs_get_axis_info( context, pjCrs.get(), 0,
+                           nullptr,
+                           nullptr,
+                           &outDirection,
+                           nullptr,
+                           nullptr,
+                           nullptr,
+                           nullptr
+                         );
+    return QString( outDirection ).compare( QLatin1String( "north" ), Qt::CaseInsensitive ) == 0;
+  }
+  return false;
+}
+
 #endif
+
diff --git a/src/core/qgsprojutils.h b/src/core/qgsprojutils.h
@@ -79,7 +79,13 @@ class CORE_EXPORT QgsProjUtils
      */
     static bool usesAngularUnit( const QString &projDef );
 
+    //TODO - remove when proj 6.1 is minimum supported version, and replace with proj_normalize_for_visualization
 
+    /**
+     * Returns true if the given proj coordinate system uses requires y/x coordinate
+     * order instead of x/y.
+     */
+    static bool axisOrderIsSwapped( const PJ *crs );
 
 #endif
 #endif

diff --git a/tests/src/core/testqgscoordinatetransform.cpp b/tests/src/core/testqgscoordinatetransform.cpp
@@ -38,6 +38,8 @@ class TestQgsCoordinateTransform: public QObject
     void contextShared();
     void scaleFactor();
     void scaleFactor_data();
+    void transform_data();
+    void transform();
 
   private:
 
@@ -252,6 +254,44 @@ void TestQgsCoordinateTransform::scaleFactor_data()
       << 1.0;
 }
 
+void TestQgsCoordinateTransform::transform_data()
+{
+  QTest::addColumn<QgsCoordinateReferenceSystem>( "sourceCrs" );
+  QTest::addColumn<QgsCoordinateReferenceSystem>( "destCrs" );
+  QTest::addColumn<double>( "x" );
+  QTest::addColumn<double>( "y" );
+  QTest::addColumn<double>( "outX" );
+  QTest::addColumn<double>( "outY" );
+  QTest::addColumn<double>( "precision" );
+
+  QTest::newRow( "To geographic" )
+      << QgsCoordinateReferenceSystem::fromEpsgId( 3111 )
+      << QgsCoordinateReferenceSystem::fromEpsgId( 4326 )
+      << 2545059.0 << 2393190.0 << 145.512750 << -37.961375 << 0.000001;
+  QTest::newRow( "From geographic" )
+      << QgsCoordinateReferenceSystem::fromEpsgId( 4326 )
+      << QgsCoordinateReferenceSystem::fromEpsgId( 3111 )
+      << 145.512750 <<  -37.961375 << 2545059.0 << 2393190.0 << 0.1;
+}
+
+void TestQgsCoordinateTransform::transform()
+{
+  QFETCH( QgsCoordinateReferenceSystem, sourceCrs );
+  QFETCH( QgsCoordinateReferenceSystem, destCrs );
+  QFETCH( double, x );
+  QFETCH( double, y );
+  QFETCH( double, outX );
+  QFETCH( double, outY );
+  QFETCH( double, precision );
+
+  double z = 0;
+  QgsCoordinateTransform ct( sourceCrs, destCrs, QgsProject::instance() );
+
+  ct.transformInPlace( x, y, z );
+  QGSCOMPARENEAR( x, outX, precision );
+  QGSCOMPARENEAR( y, outY, precision );
+}
+
 void TestQgsCoordinateTransform::transformBoundingBox()
 {
   //test transforming a bounding box which crosses the 180 degree longitude line

diff --git a/tests/src/core/testqgsprojutils.cpp b/tests/src/core/testqgsprojutils.cpp
@@ -18,6 +18,10 @@
 #include "qgsapplication.h"
 #include "qgslogger.h"
 
+#if PROJ_VERSION_MAJOR>=6
+#include <proj.h>
+#endif
+
 //header for class being tested
 #include "qgsprojutils.h"
 #include <QtConcurrent>
@@ -30,6 +34,7 @@ class TestQgsProjUtils: public QObject
     void cleanupTestCase();
     void threadSafeContext();
     void usesAngularUnits();
+    void axisOrderIsSwapped();
 
 };
 
@@ -79,5 +84,18 @@ void TestQgsProjUtils::usesAngularUnits()
 #endif
 }
 
+void TestQgsProjUtils::axisOrderIsSwapped()
+{
+#if PROJ_VERSION_MAJOR>=6
+  PJ_CONTEXT *context = QgsProjContext::get();
+  QVERIFY( !QgsProjUtils::axisOrderIsSwapped( nullptr ) );
+
+  QgsProjUtils::proj_pj_unique_ptr crs( proj_create( context, "urn:ogc:def:crs:EPSG::3111" ) );
+  QVERIFY( !QgsProjUtils::axisOrderIsSwapped( crs.get() ) );
+  crs.reset( proj_create( context, "urn:ogc:def:crs:EPSG::4326" ) );
+  QVERIFY( QgsProjUtils::axisOrderIsSwapped( crs.get() ) );
+#endif
+}
+
 QGSTEST_MAIN( TestQgsProjUtils )
 #include "testqgsprojutils.moc"