Add API to retrieve avaiable coordinate operations between a source

and destination QgsCoordinateReferenceSystem on proj >= 6 builds

python/core/auto_generated/qgsdatumtransform.sip.in

@@ -70,6 +70,49 @@ and ``destinationTransformId`` transforms.
 
     };
 
+
+    struct GridDetails
+    {
+      QString shortName;
+      QString fullName;
+      QString packageName;
+      QString url;
+      bool directDownload;
+      bool openLicense;
+      bool isAvailable;
+    };
+
+    struct TransformDetails
+    {
+      QString proj;
+      QString name;
+      double accuracy;
+
+      bool isAvailable;
+
+      QList< QgsDatumTransform::GridDetails > grids;
+    };
+
+    static QList< QgsDatumTransform::TransformDetails > operations( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination );
+%Docstring
+Returns a list of coordinate operations available for transforming
+coordinates from the ``source`` to ``destination`` CRS.
+
+This list is sorted in order of preference, with the most preferable
+operation listed first.
+
+Not all operations may be available for use. Check QgsDatumTransform.TransformDetails.isAvailable
+first. Operations may require grid transformation files which are not available on the local
+install.
+
+.. note::
+
+   Requires Proj 6.0 or later. Builds based on earlier Proj versions will always return an empty list,
+   and the deprecated API from QgsDatumTransform must be used instead.
+
+.. versionadded:: 3.8
+%End
+
  static QList< QgsDatumTransform::TransformPair > datumTransformations( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination ) /Deprecated/;
 %Docstring
 Returns a list of datum transformations which are available for the given ``source`` and ``destination`` CRS.

src/core/qgscoordinatereferencesystem.cpp

@@ -2971,7 +2971,7 @@ QString QgsCoordinateReferenceSystem::geographicCrsAuthId() const
 }
 
 #if PROJ_VERSION_MAJOR>=6
-PJ *QgsCoordinateReferenceSystem::projObject()
+PJ *QgsCoordinateReferenceSystem::projObject() const
 {
   return d->mPj.get();
 }

src/core/qgscoordinatereferencesystem.h

@@ -660,7 +660,7 @@ class CORE_EXPORT QgsCoordinateReferenceSystem
      * \since QGIS 3.8
      * \note Not available in Python bindings.
      */
-    PJ *projObject();
+    PJ *projObject() const;
 #endif
 #endif
 

src/core/qgsdatumtransform.cpp

@@ -20,6 +20,71 @@
 #include "qgssqliteutils.h"
 #include <sqlite3.h>
 
+#if PROJ_VERSION_MAJOR>=6
+#include "qgsprojutils.h"
+#include <proj.h>
+#endif
+
+QList<QgsDatumTransform::TransformDetails> QgsDatumTransform::operations( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination )
+{
+  QList< QgsDatumTransform::TransformDetails > res;
+#if PROJ_VERSION_MAJOR>=6
+  if ( !source.projObject() || !destination.projObject() )
+    return res;
+
+  PJ_CONTEXT *pjContext = QgsProjContext::get();
+
+  PJ_OPERATION_FACTORY_CONTEXT *operationContext = proj_create_operation_factory_context( pjContext, nullptr );
+
+  // We want to return ALL grids, not just those available for use
+  proj_operation_factory_context_set_grid_availability_use( pjContext, operationContext, PROJ_GRID_AVAILABILITY_IGNORED );
+
+  // See https://lists.osgeo.org/pipermail/proj/2019-May/008604.html
+  proj_operation_factory_context_set_spatial_criterion( pjContext,   operationContext,  PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION );
+
+  PJ_OBJ_LIST *ops = proj_create_operations( pjContext, source.projObject(), destination.projObject(), operationContext );
+  int count = proj_list_get_count( ops );
+  for ( int i = 0; i < count; ++i )
+  {
+    PJ *op = proj_list_get( pjContext, ops, i );
+    if ( !op )
+      continue;
+
+    TransformDetails details;
+    details.proj = QString( proj_as_proj_string( pjContext, op, PJ_PROJ_5, nullptr ) );
+    details.name = QString( proj_get_name( op ) );
+    details.accuracy = proj_coordoperation_get_accuracy( pjContext, op );
+    details.isAvailable = proj_coordoperation_is_instantiable( pjContext, op );
+
+    for ( int j = 0; j < proj_coordoperation_get_grid_used_count( pjContext, op ); ++j )
+    {
+      const char *shortName = nullptr;
+      const char *fullName = nullptr;
+      const char *packageName = nullptr;
+      const char *url = nullptr;
+      int directDownload = 0;
+      int openLicense = 0;
+      int isAvailable = 0;
+      proj_coordoperation_get_grid_used( pjContext, op, j, &shortName, &fullName, &packageName, &url, &directDownload, &openLicense, &isAvailable );
+      GridDetails gridDetails;
+      gridDetails.shortName = QString( shortName );
+      gridDetails.fullName = QString( fullName );
+      gridDetails.packageName = QString( packageName );
+      gridDetails.url = QString( url );
+      gridDetails.directDownload = directDownload;
+      gridDetails.openLicense = openLicense;
+      gridDetails.isAvailable = isAvailable;
+
+      details.grids.append( gridDetails );
+    }
+    res.push_back( details );
+  }
+  proj_list_destroy( ops );
+  proj_operation_factory_context_destroy( operationContext );
+#endif
+  return res;
+}
+
 QList< QgsDatumTransform::TransformPair > QgsDatumTransform::datumTransformations( const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS )
 {
   QList< QgsDatumTransform::TransformPair > transformations;

src/core/qgsdatumtransform.h

@@ -20,6 +20,7 @@
 #include "qgis_core.h"
 #include "qgis_sip.h"
 #include <QString>
+#include <QList>
 
 class QgsCoordinateReferenceSystem;
 
@@ -117,6 +118,74 @@ class CORE_EXPORT QgsDatumTransform
 
     };
 
+
+    /**
+     * Contains information about a projection transformation grid file.
+     * \since QGIS 3.8
+     */
+    struct GridDetails
+    {
+      //! Short name of transform grid
+      QString shortName;
+      //! Full name of transform grid
+      QString fullName;
+      //! Name of package the grid is included within
+      QString packageName;
+      //! Url to download grid from
+      QString url;
+      //! TRUE if direct download of grid is possible
+      bool directDownload = false;
+      //! TRUE if grid is available under an open license
+      bool openLicense = false;
+      //! TRUE if grid is currently available for use
+      bool isAvailable = false;
+    };
+
+    /**
+     * Contains information about a coordinate transformation operation.
+     * \since QGIS 3.8
+     */
+    struct TransformDetails
+    {
+      //! Proj representation of transform operation
+      QString proj;
+      //! Display name of transform operation
+      QString name;
+      //! Transformation accuracy (in meters)
+      double accuracy = 0;
+
+      /**
+       * TRUE if operation is available.
+       *
+       * If FALSE, it likely means a transform grid is required which is not
+       * available.
+       */
+      bool isAvailable = false;
+
+      /**
+       * Contains a list of transform grids used by the operation.
+       */
+      QList< QgsDatumTransform::GridDetails > grids;
+    };
+
+    /**
+     * Returns a list of coordinate operations available for transforming
+     * coordinates from the \a source to \a destination CRS.
+     *
+     * This list is sorted in order of preference, with the most preferable
+     * operation listed first.
+     *
+     * Not all operations may be available for use. Check QgsDatumTransform::TransformDetails::isAvailable
+     * first. Operations may require grid transformation files which are not available on the local
+     * install.
+     *
+     * \note Requires Proj 6.0 or later. Builds based on earlier Proj versions will always return an empty list,
+     * and the deprecated API from QgsDatumTransform must be used instead.
+     *
+     * \since QGIS 3.8
+     */
+    static QList< QgsDatumTransform::TransformDetails > operations( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination );
+
     /**
      * Returns a list of datum transformations which are available for the given \a source and \a destination CRS.
      * \see datumTransformToProj()

tests/src/python/CMakeLists.txt

@@ -40,6 +40,7 @@ ADD_PYTHON_TEST(PyQgsDataItemGuiProviderRegistry test_qgsdataitemguiproviderregi
 ADD_PYTHON_TEST(PyQgsDataItemProviderRegistry test_qgsdataitemproviderregistry.py)
 ADD_PYTHON_TEST(PyQgsDateTimeEdit test_qgsdatetimeedit.py)
 ADD_PYTHON_TEST(PyQgsDateTimeStatisticalSummary test_qgsdatetimestatisticalsummary.py)
+ADD_PYTHON_TEST(PyQgsDatumTransform test_qgsdatumtransforms.py)
 ADD_PYTHON_TEST(PyQgsDelimitedTextProvider test_qgsdelimitedtextprovider.py)
 ADD_PYTHON_TEST(PyQgsDistanceArea test_qgsdistancearea.py)
 ADD_PYTHON_TEST(PyQgsEditFormConfig test_qgseditformconfig.py)

tests/src/python/test_qgsdatumtransforms.py

@@ -0,0 +1,148 @@
+
+# -*- coding: utf-8 -*-
+"""QGIS Unit tests for QgsDatumTransforms.
+
+.. note:: This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+"""
+__author__ = 'Nyall Dawson'
+__date__ = '2019-05-25'
+__copyright__ = 'Copyright 2019, The QGIS Project'
+
+from qgis.core import (
+    QgsProjUtils,
+    QgsCoordinateReferenceSystem,
+    QgsDatumTransform
+)
+from qgis.testing import (start_app,
+                          unittest,
+                          )
+from utilities import unitTestDataPath
+
+start_app()
+TEST_DATA_DIR = unitTestDataPath()
+
+
+class TestPyQgsDatumTransform(unittest.TestCase):
+
+    @unittest.skipIf(QgsProjUtils.projVersionMajor() < 6, 'Not a proj6 build')
+    def testOperations(self):
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem(),
+                                           QgsCoordinateReferenceSystem())
+        self.assertEqual(ops, [])
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem('EPSG:3111'),
+                                           QgsCoordinateReferenceSystem())
+        self.assertEqual(ops, [])
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem(),
+                                           QgsCoordinateReferenceSystem('EPSG:3111'))
+        self.assertEqual(ops, [])
+
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem('EPSG:3111'),
+                                           QgsCoordinateReferenceSystem('EPSG:3111'))
+        self.assertEqual(len(ops), 1)
+        self.assertEqual(ops[0].name, 'Inverse of Vicgrid + Vicgrid')
+        self.assertEqual(ops[0].proj, '+proj=noop')
+        self.assertEqual(ops[0].accuracy, 0.0)
+        self.assertTrue(ops[0].isAvailable)
+
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem('EPSG:3111'),
+                                           QgsCoordinateReferenceSystem('EPSG:4283'))
+        self.assertEqual(len(ops), 1)
+        self.assertEqual(ops[0].name, 'Inverse of Vicgrid')
+        self.assertEqual(ops[0].proj, '+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1')
+        self.assertEqual(ops[0].accuracy, -1.0)
+        self.assertTrue(ops[0].isAvailable)
+
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem('EPSG:3111'),
+                                           QgsCoordinateReferenceSystem('EPSG:28355'))
+        self.assertEqual(len(ops), 1)
+        self.assertEqual(ops[0].name, 'Inverse of Vicgrid + Map Grid of Australia zone 55')
+        self.assertEqual(ops[0].proj, '+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step +proj=utm +zone=55 +south +ellps=GRS80')
+        self.assertEqual(ops[0].accuracy, 0.0)
+        self.assertTrue(ops[0].isAvailable)
+
+        # uses a grid file
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem('EPSG:4283'),
+                                           QgsCoordinateReferenceSystem('EPSG:7844'))
+        self.assertEqual(len(ops), 5)
+        self.assertEqual(ops[0].name, 'GDA94 to GDA2020 (1)')
+        self.assertEqual(ops[0].proj, '+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1')
+        self.assertTrue(ops[0].isAvailable)
+        self.assertEqual(ops[0].accuracy, 0.01)
+        self.assertEqual(len(ops[0].grids), 0)
+        self.assertEqual(ops[1].name, 'GDA94 to GDA2020 (2)')
+        self.assertEqual(ops[1].proj, '+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GDA94_GDA2020_conformal_and_distortion.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1')
+        self.assertEqual(ops[1].accuracy, 0.05) # actually incorrect in EPSG registry, may need updating
+        self.assertEqual(len(ops[1].grids), 1)
+        self.assertEqual(ops[1].grids[0].shortName, 'GDA94_GDA2020_conformal_and_distortion.gsb')
+        self.assertEqual(ops[1].grids[0].fullName, '')
+        self.assertEqual(ops[1].grids[0].packageName, 'proj-datumgrid-oceania')
+        self.assertEqual(ops[1].grids[0].url, 'https://download.osgeo.org/proj/proj-datumgrid-oceania-latest.zip')
+        self.assertTrue(ops[1].grids[0].directDownload)
+        self.assertTrue(ops[1].grids[0].openLicense)
+        self.assertEqual(ops[2].name, 'GDA94 to GDA2020 (3)')
+        self.assertEqual(ops[2].proj, '+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GDA94_GDA2020_conformal.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1')
+        self.assertEqual(ops[2].accuracy, 0.05) # actually incorrect in EPSG registry, may need updating
+        self.assertEqual(len(ops[2].grids), 1)
+        self.assertEqual(ops[2].grids[0].shortName, 'GDA94_GDA2020_conformal.gsb')
+        self.assertEqual(ops[2].grids[0].fullName, '')
+        self.assertEqual(ops[2].grids[0].packageName, 'proj-datumgrid-oceania')
+        self.assertEqual(ops[2].grids[0].url, 'https://download.osgeo.org/proj/proj-datumgrid-oceania-latest.zip')
+        self.assertTrue(ops[2].grids[0].directDownload)
+        self.assertTrue(ops[2].grids[0].openLicense)
+        self.assertEqual(ops[3].name, 'GDA94 to GDA2020 (5)')
+        self.assertEqual(ops[3].proj, '+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GDA94_GDA2020_conformal_cocos_island.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1')
+        self.assertEqual(ops[3].accuracy, 0.05) # actually incorrect in EPSG registry, may need updating
+        self.assertEqual(len(ops[3].grids), 1)
+        self.assertEqual(ops[3].grids[0].shortName, 'GDA94_GDA2020_conformal_cocos_island.gsb')
+        self.assertEqual(ops[3].grids[0].fullName, '')
+        self.assertEqual(ops[3].grids[0].packageName, '')
+        self.assertEqual(ops[3].grids[0].url, '')
+        self.assertFalse(ops[3].grids[0].directDownload)
+        self.assertFalse(ops[3].grids[0].openLicense)
+        self.assertEqual(ops[4].name, 'GDA94 to GDA2020 (4)')
+        self.assertEqual(ops[4].proj, '+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GDA94_GDA2020_conformal_christmas_island.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1')
+        self.assertEqual(ops[4].accuracy, 0.05) # actually incorrect in EPSG registry, may need updating
+        self.assertEqual(len(ops[4].grids), 1)
+        self.assertEqual(ops[4].grids[0].shortName, 'GDA94_GDA2020_conformal_christmas_island.gsb')
+        self.assertEqual(ops[4].grids[0].fullName, '')
+        self.assertEqual(ops[4].grids[0].packageName, '')
+        self.assertEqual(ops[4].grids[0].url, '')
+        self.assertFalse(ops[4].grids[0].directDownload)
+        self.assertFalse(ops[4].grids[0].openLicense)
+
+        # uses a pivot datum
+        ops = QgsDatumTransform.operations(QgsCoordinateReferenceSystem('EPSG:3111'),
+                                           QgsCoordinateReferenceSystem('EPSG:7899'))
+        self.assertEqual(len(ops), 3)
+        self.assertEqual(ops[0].name, 'Inverse of Vicgrid + GDA94 to GDA2020 (1) + Vicgrid')
+        self.assertEqual(ops[0].proj, '+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step +proj=push +v_3 +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80')
+        self.assertTrue(ops[0].isAvailable)
+        self.assertEqual(ops[0].accuracy, 0.01)
+        self.assertEqual(len(ops[0].grids), 0)
+        self.assertEqual(ops[1].name, 'Inverse of Vicgrid + GDA94 to GDA2020 (2) + Vicgrid')
+        self.assertEqual(ops[1].proj, '+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step +proj=hgridshift +grids=GDA94_GDA2020_conformal_and_distortion.gsb +step +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80')
+        self.assertEqual(ops[1].accuracy, 0.05) # actually incorrect in EPSG registry, may need updating
+        self.assertEqual(len(ops[1].grids), 1)
+        self.assertEqual(ops[1].grids[0].shortName, 'GDA94_GDA2020_conformal_and_distortion.gsb')
+        self.assertEqual(ops[1].grids[0].fullName, '')
+        self.assertEqual(ops[1].grids[0].packageName, 'proj-datumgrid-oceania')
+        self.assertEqual(ops[1].grids[0].url, 'https://download.osgeo.org/proj/proj-datumgrid-oceania-latest.zip')
+        self.assertTrue(ops[1].grids[0].directDownload)
+        self.assertTrue(ops[1].grids[0].openLicense)
+        self.assertEqual(ops[2].name, 'Inverse of Vicgrid + GDA94 to GDA2020 (3) + Vicgrid')
+        self.assertEqual(ops[2].proj, '+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step +proj=hgridshift +grids=GDA94_GDA2020_conformal.gsb +step +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80')
+        self.assertEqual(ops[2].accuracy, 0.05) # actually incorrect in EPSG registry, may need updating
+        self.assertEqual(len(ops[2].grids), 1)
+        self.assertEqual(ops[2].grids[0].shortName, 'GDA94_GDA2020_conformal.gsb')
+        self.assertEqual(ops[2].grids[0].fullName, '')
+        self.assertEqual(ops[2].grids[0].packageName, 'proj-datumgrid-oceania')
+        self.assertEqual(ops[2].grids[0].url, 'https://download.osgeo.org/proj/proj-datumgrid-oceania-latest.zip')
+        self.assertTrue(ops[2].grids[0].directDownload)
+        self.assertTrue(ops[2].grids[0].openLicense)
+
+
+if __name__ == '__main__':
+    unittest.main()