Add QgsRasterDataProvider::transformCoordinates()

wonder-sk · nyalldawson · commit 68bbf46ebb5f · 2020-05-18T09:20:32.000+10:00
This is useful when client needs to find out image space coordinates of a point
in map layer coordinates or vice versa. For warped VRT rasters this can't be
simply done by using geotransform matrix because the transform may be more complex.

This may be also useful functionality for identify tool to show source raster
image coordinates.
diff --git a/python/core/auto_generated/raster/qgsrasterdataprovider.sip.in b/python/core/auto_generated/raster/qgsrasterdataprovider.sip.in
@@ -495,6 +495,26 @@ Returns true if the extents reported by the data provider are not reliable
 and it's possible that there is renderable content outside of these extents.
 
 .. versionadded:: 3.10.0
+%End
+
+    enum TransformType
+    {
+      TransformImageToLayer,
+      TransformLayerToImage,
+    };
+
+    virtual QgsPoint transformCoordinates( const QgsPoint &point, TransformType type );
+%Docstring
+Transforms coordinates between source image coordinate space [0..width]x[0..height] and
+layer coordinate space (georeferenced coordinates). Often this transformation is a simple
+2D affine transformation (offset and scaling), but rasters with different georeferencing
+methods like GCPs (ground control points) or RPCs (rational polynomial coefficients) may
+require a more complex transform.
+
+If the transform fails (input coordinates are outside of the valid range or data provider
+does not support this functionality), an empty point is returned.
+
+.. versionadded:: 3.14
 %End
 
   signals:
diff --git a/src/core/providers/gdal/qgsgdalprovider.cpp b/src/core/providers/gdal/qgsgdalprovider.cpp
@@ -457,6 +457,9 @@ QgsGdalProvider::~QgsGdalProvider()
 {
   QMutexLocker locker( sGdalProviderMutex() );
 
+  if ( mGdalTransformerArg )
+    GDALDestroyTransformer( mGdalTransformerArg );
+
   int lightRefCounter = -- ( *mpLightRefCounter );
   int refCounter = -- ( *mpRefCounter );
   if ( refCounter == 0 )
@@ -2686,6 +2689,8 @@ void QgsGdalProvider::initBaseDataset()
     mGdalDataset = mGdalBaseDataset;
   }
 
+  mGdalTransformerArg = GDALCreateGenImgProjTransformer( mGdalBaseDataset, nullptr, nullptr, nullptr, TRUE, 1.0, 0 );
+
   if ( !hasGeoTransform )
   {
     // Initialize the affine transform matrix
@@ -3156,6 +3161,20 @@ QString QgsGdalProvider::validatePyramidsConfigOptions( QgsRaster::RasterPyramid
   return QString();
 }
 
+QgsPoint QgsGdalProvider::transformCoordinates( const QgsPoint &point, QgsRasterDataProvider::TransformType type )
+{
+  if ( !mGdalTransformerArg )
+    return QgsPoint();
+
+  int success;
+  double x = point.x(), y = point.y(), z = point.is3D() ? point.z() : 0;
+  GDALUseTransformer( mGdalTransformerArg, type == TransformLayerToImage, 1, &x, &y, &z, &success );
+  if ( !success )
+    return QgsPoint();
+
+  return QgsPoint( x, y, z );
+}
+
 bool QgsGdalProvider::isEditable() const
 {
   return mUpdate;
diff --git a/src/core/providers/gdal/qgsgdalprovider.h b/src/core/providers/gdal/qgsgdalprovider.h
@@ -203,6 +203,8 @@ class QgsGdalProvider final: public QgsRasterDataProvider, QgsGdalProviderBase
     QString validatePyramidsConfigOptions( QgsRaster::RasterPyramidsFormat pyramidsFormat,
                                            const QStringList &configOptions, const QString &fileFormat ) override;
 
+    QgsPoint transformCoordinates( const QgsPoint &point, TransformType type ) override;
+
   private:
     QgsGdalProvider( const QgsGdalProvider &other );
 
@@ -338,6 +340,9 @@ class QgsGdalProvider final: public QgsRasterDataProvider, QgsGdalProviderBase
      * Closes and reinits dataset
     */
     void reloadProviderData() override;
+
+    //! Instance of GDAL transformer function used in transformCoordinates() for conversion between image and layer coordinates
+    void *mGdalTransformerArg = nullptr;
 };
 
 /**
diff --git a/src/core/raster/qgsrasterdataprovider.cpp b/src/core/raster/qgsrasterdataprovider.cpp
@@ -512,6 +512,13 @@ bool QgsRasterDataProvider::ignoreExtents() const
   return false;
 }
 
+QgsPoint QgsRasterDataProvider::transformCoordinates( const QgsPoint &point, QgsRasterDataProvider::TransformType type )
+{
+  Q_UNUSED( point )
+  Q_UNUSED( type )
+  return QgsPoint();
+}
+
 bool QgsRasterDataProvider::userNoDataValuesContains( int bandNo, double value ) const
 {
   QgsRasterRangeList rangeList = mUserNoDataValue.value( bandNo - 1 );
diff --git a/src/core/raster/qgsrasterdataprovider.h b/src/core/raster/qgsrasterdataprovider.h
@@ -557,6 +557,30 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
      */
     virtual bool ignoreExtents() const;
 
+    /**
+     * Types of transformation in transformCoordinates() function.
+     * \since QGIS 3.14
+     */
+    enum TransformType
+    {
+      TransformImageToLayer,  //!< Transforms image coordinates to layer (georeferenced) coordinates
+      TransformLayerToImage,  //!< Transforms layer (georeferenced) coordinates to image coordinates
+    };
+
+    /**
+     * Transforms coordinates between source image coordinate space [0..width]x[0..height] and
+     * layer coordinate space (georeferenced coordinates). Often this transformation is a simple
+     * 2D affine transformation (offset and scaling), but rasters with different georeferencing
+     * methods like GCPs (ground control points) or RPCs (rational polynomial coefficients) may
+     * require a more complex transform.
+     *
+     * If the transform fails (input coordinates are outside of the valid range or data provider
+     * does not support this functionality), an empty point is returned.
+     *
+     * \since QGIS 3.14
+     */
+    virtual QgsPoint transformCoordinates( const QgsPoint &point, TransformType type );
+
   signals:
 
     /**

Original file line number	Diff line number	Diff line change
`@@ -512,6 +512,13 @@ bool QgsRasterDataProvider::ignoreExtents() const`
`512`	`512`	`return false;`
`513`	`513`	`}`
`514`	`514`
	`515`	`+QgsPoint QgsRasterDataProvider::transformCoordinates( const QgsPoint &point, QgsRasterDataProvider::TransformType type )`
	`516`	`+{`
	`517`	`+ Q_UNUSED( point )`
	`518`	`+ Q_UNUSED( type )`
	`519`	`+ return QgsPoint();`
	`520`	`+}`
	`521`	`+`
`515`	`522`	`bool QgsRasterDataProvider::userNoDataValuesContains( int bandNo, double value ) const`
`516`	`523`	`{`
`517`	`524`	`QgsRasterRangeList rangeList = mUserNoDataValue.value( bandNo - 1 );`