[FEATURE] Add QgsRasterDataProvider::sample method for efficient

sampling of rasters at a given point

This is an alternative to the ::identify method, which is less
efficient but more powerful

python/core/auto_generated/raster/qgsrasterdataprovider.sip.in

@@ -246,6 +246,50 @@ into a subset of the GUI raster properties "Metadata" tab.
 %End
 
     virtual QgsRasterIdentifyResult identify( const QgsPointXY &point, QgsRaster::IdentifyFormat format, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 );
+%Docstring
+Identify raster value(s) found on the point position. The context
+parameters extent, width and height are important to identify
+on the same zoom level as a displayed map and to do effective
+caching (WCS). If context params are not specified the highest
+resolution is used. capabilities() may be used to test if format
+is supported by provider. Values are set to 'no data' or empty string
+if point is outside data source extent.
+
+:param point: coordinates in data source CRS
+:param format: result format
+:param boundingBox: context bounding box
+:param width: context width
+:param height: context height
+:param dpi: context dpi
+
+:return: QgsRaster.IdentifyFormatValue: map of values for each band, keys are band numbers
+         (from 1).
+         QgsRaster.IdentifyFormatFeature: map of QgsRasterFeatureList for each sublayer
+         (WMS) - TODO: it is not consistent with QgsRaster.IdentifyFormatValue.
+         QgsRaster.IdentifyFormatHtml: map of HTML strings for each sublayer (WMS).
+         Empty if failed or there are no results (TODO: better error reporting).
+
+.. note::
+
+   The arbitraryness of the returned document is enforced by WMS standards
+   up to at least v1.3.0
+
+.. seealso:: :py:func:`sample`
+%End
+
+    virtual QVariant sample( const QgsPointXY &point, int band, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 );
+%Docstring
+Samples a raster value from the specified ``band`` found at the ``point`` position. The context
+parameters ``boundingBox``, ``width`` and ``height`` are important to identify
+on the same zoom level as a displayed map and to do effective
+caching (WCS). If context params are not specified the highest
+resolution is used. capabilities() may be used to test if format
+is supported by provider.
+
+A null QVariant will be returned if the point is outside data source extent.
+
+.. versionadded:: 3.4
+%End
 
     virtual QString lastErrorTitle() = 0;
 %Docstring

src/core/raster/qgsrasterdataprovider.cpp

@@ -300,6 +300,45 @@ QgsRasterIdentifyResult QgsRasterDataProvider::identify( const QgsPointXY &point
   return QgsRasterIdentifyResult( QgsRaster::IdentifyFormatValue, results );
 }
 
+QVariant QgsRasterDataProvider::sample( const QgsPointXY &point, int band, const QgsRectangle &boundingBox, int width, int height, int )
+{
+  if ( !extent().contains( point ) )
+  {
+    // Outside the raster
+    return QVariant();
+  }
+
+  QgsRectangle finalExtent = boundingBox;
+  if ( finalExtent.isEmpty() )
+    finalExtent = extent();
+
+  if ( width == 0 )
+  {
+    width = capabilities() & Size ? xSize() : 1000;
+  }
+  if ( height == 0 )
+  {
+    height = capabilities() & Size ? ySize() : 1000;
+  }
+
+  // Calculate the row / column where the point falls
+  double xres = ( finalExtent.width() ) / width;
+  double yres = ( finalExtent.height() ) / height;
+
+  int col = static_cast< int >( std::floor( ( point.x() - finalExtent.xMinimum() ) / xres ) );
+  int row = static_cast< int >( std::floor( ( finalExtent.yMaximum() - point.y() ) / yres ) );
+
+  double xMin = finalExtent.xMinimum() + col * xres;
+  double xMax = xMin + xres;
+  double yMax = finalExtent.yMaximum() - row * yres;
+  double yMin = yMax - yres;
+  QgsRectangle pixelExtent( xMin, yMin, xMax, yMax );
+
+  std::unique_ptr< QgsRasterBlock > bandBlock( block( band, pixelExtent, 1, 1 ) );
+
+  return bandBlock ? bandBlock->value( 0 ) : QVariant();
+}
+
 QString QgsRasterDataProvider::lastErrorFormat()
 {
   return QStringLiteral( "text/plain" );

src/core/raster/qgsrasterdataprovider.h

@@ -355,10 +355,24 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
      *         Empty if failed or there are no results (TODO: better error reporting).
      * \note The arbitraryness of the returned document is enforced by WMS standards
      *       up to at least v1.3.0
+     * \see sample()
      */
-    //virtual QMap<int, QVariant> identify( const QgsPointXY & point, QgsRaster::IdentifyFormat format, const QgsRectangle &extent = QgsRectangle(), int width = 0, int height = 0 );
     virtual QgsRasterIdentifyResult identify( const QgsPointXY &point, QgsRaster::IdentifyFormat format, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 );
 
+    /**
+     * Samples a raster value from the specified \a band found at the \a point position. The context
+     * parameters \a boundingBox, \a width and \a height are important to identify
+     * on the same zoom level as a displayed map and to do effective
+     * caching (WCS). If context params are not specified the highest
+     * resolution is used. capabilities() may be used to test if format
+     * is supported by provider.
+     *
+     * A null QVariant will be returned if the point is outside data source extent.
+     *
+     * \since QGIS 3.4
+     */
+    virtual QVariant sample( const QgsPointXY &point, int band, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 );
+
     /**
      * \brief Returns the caption error text for the last error in this provider
      *

src/providers/gdal/qgsgdalprovider.cpp

@@ -1089,6 +1089,74 @@ QgsRasterIdentifyResult QgsGdalProvider::identify( const QgsPointXY &point, QgsR
   return QgsRasterIdentifyResult( QgsRaster::IdentifyFormatValue, results );
 }
 
+QVariant QgsGdalProvider::sample( const QgsPointXY &point, int band, const QgsRectangle &boundingBox, int width, int height, int )
+{
+  QMutexLocker locker( mpMutex );
+  if ( !initIfNeeded() )
+    return tr( "Cannot read data" );
+
+  if ( !extent().contains( point ) )
+  {
+    // Outside the raster
+    return QVariant();
+  }
+
+  QgsRectangle finalExtent = boundingBox;
+  if ( finalExtent.isEmpty() )
+    finalExtent = extent();
+
+  if ( width == 0 )
+    width = xSize();
+  if ( height == 0 )
+    height = ySize();
+
+  // Calculate the row / column where the point falls
+  double xres = ( finalExtent.width() ) / width;
+  double yres = ( finalExtent.height() ) / height;
+
+  // Offset, not the cell index -> std::floor
+  int col = static_cast< int >( std::floor( ( point.x() - finalExtent.xMinimum() ) / xres ) );
+  int row = static_cast< int >( std::floor( ( finalExtent.yMaximum() - point.y() ) / yres ) );
+
+  int r = 0;
+  int c = 0;
+  int w = 1;
+  int h = 1;
+
+  double xMin = finalExtent.xMinimum() + col * xres;
+  double xMax = xMin + xres * w;
+  double yMax = finalExtent.yMaximum() - row * yres;
+  double yMin = yMax - yres * h;
+  QgsRectangle pixelExtent( xMin, yMin, xMax, yMax );
+
+  std::unique_ptr< QgsRasterBlock > bandBlock( block( band, pixelExtent, w, h ) );
+
+  if ( !bandBlock )
+  {
+    return tr( "Cannot read data" );
+  }
+
+  double value = bandBlock->value( r, c );
+
+  if ( ( sourceHasNoDataValue( band ) && useSourceNoDataValue( band ) &&
+         ( std::isnan( value ) || qgsDoubleNear( value, sourceNoDataValue( band ) ) ) ) ||
+       ( QgsRasterRange::contains( value, userNoDataValues( band ) ) ) )
+  {
+    return QVariant(); // null QVariant represents no data
+  }
+  else
+  {
+    if ( sourceDataType( band ) == Qgis::Float32 )
+    {
+      // Insert a float QVariant so that QgsMapToolIdentify::identifyRasterLayer()
+      // can print a string without an excessive precision
+      return static_cast<float>( value );
+    }
+    else
+      return value;
+  }
+}
+
 int QgsGdalProvider::capabilities() const
 {
   QMutexLocker locker( mpMutex );

src/providers/gdal/qgsgdalprovider.h

@@ -103,6 +103,7 @@ class QgsGdalProvider : public QgsRasterDataProvider, QgsGdalProviderBase
     QgsRectangle extent() const override;
     bool isValid() const override;
     QgsRasterIdentifyResult identify( const QgsPointXY &point, QgsRaster::IdentifyFormat format, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 ) override;
+    QVariant sample( const QgsPointXY &point, int band, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 );
     QString lastErrorTitle() override;
     QString lastError() override;
     int capabilities() const override;