More heavily optimised sample method for gdal provider

qgis · Jul 19, 2018 · 54e5119 · 54e5119
1 parent db8e536
commit 54e5119
Show file tree

Hide file tree

Showing 5 changed files with 34 additions and 54 deletions.
diff --git a/python/core/auto_generated/raster/qgsrasterdataprovider.sip.in b/python/core/auto_generated/raster/qgsrasterdataprovider.sip.in
@@ -286,7 +286,8 @@ 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.
+A null QVariant will be returned if the point is outside data source extent, or an invalid
+band number was specified.
 
 .. seealso:: :py:func:`identify`
 

diff --git a/src/core/raster/qgsrasterdataprovider.h b/src/core/raster/qgsrasterdataprovider.h
@@ -367,7 +367,8 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
      * 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.
+     * A null QVariant will be returned if the point is outside data source extent, or an invalid
+     * band number was specified.
      *
      * \see identify(), which is much more flexible but considerably less efficient.
      * \since QGIS 3.4

diff --git a/src/providers/gdal/qgsgdalprovider.cpp b/src/providers/gdal/qgsgdalprovider.cpp
@@ -1089,7 +1089,19 @@ 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 )
+bool QgsGdalProvider::worldToPixel( double x, double y, int &col, int &row ) const
+{
+  double div = ( mGeoTransform[2] * mGeoTransform[4] - mGeoTransform[1] * mGeoTransform[5] );
+  if ( div < 2 * std::numeric_limits<double>::epsilon() )
+    return false;
+  double doubleCol = -( mGeoTransform[2] * ( mGeoTransform[3] - y ) + mGeoTransform[5] * x - mGeoTransform[0] * mGeoTransform[5] ) / div;
+  double doubleRow = ( mGeoTransform[1] * ( mGeoTransform[3] - y ) + mGeoTransform[4] * x - mGeoTransform[0] * mGeoTransform[4] ) / div;
+  col = static_cast< int >( doubleCol );
+  row = static_cast< int >( doubleRow );
+  return true;
+};
+
+QVariant QgsGdalProvider::sample( const QgsPointXY &point, int band, const QgsRectangle &, int, int, int )
 {
   QMutexLocker locker( mpMutex );
   if ( !initIfNeeded() )
@@ -1101,60 +1113,22 @@ QVariant QgsGdalProvider::sample( const QgsPointXY &point, int band, const QgsRe
     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 ) );
+  GDALRasterBandH hBand = GDALGetRasterBand( mGdalDataset, band );
+  if ( !hBand )
+    return QVariant();
 
-  if ( !bandBlock )
-  {
-    return tr( "Cannot read data" );
-  }
+  int row;
+  int col;
+  if ( !worldToPixel( point.x(), point.y(), col, row ) )
+    return QVariant();
 
-  double value = bandBlock->value( r, c );
+  float value = 0;
+  CPLErr err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
+                             &value, 1, 1, GDT_Float32, 0, 0 );
+  if ( err != CE_None )
+    return QVariant();
 
-  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;
-  }
+  return static_cast< double >( value ) * bandScale( band ) + bandOffset( band );
 }
 
 int QgsGdalProvider::capabilities() const

diff --git a/src/providers/gdal/qgsgdalprovider.h b/src/providers/gdal/qgsgdalprovider.h
@@ -292,6 +292,7 @@ class QgsGdalProvider : public QgsRasterDataProvider, QgsGdalProviderBase
     //! Close all cached dataset for the specified provider.
     static void closeCachedGdalHandlesFor( QgsGdalProvider *provider );
 
+    bool worldToPixel( double x, double y, int &col, int &row ) const;
 };
 
 #endif

diff --git a/tests/src/core/testqgsrasterlayer.cpp b/tests/src/core/testqgsrasterlayer.cpp
@@ -731,6 +731,9 @@ void TestQgsRasterLayer::sample()
   QVERIFY( rl->isValid() );
   QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1 ).isValid() );
   QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 1 ).toInt(), 125 );
+  // bad bands
+  QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 0 ).isValid() );
+  QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 10 ).isValid() );
 
   fileName = mTestDataDir + "landsat_4326.tif";
   rasterFileInfo = QFileInfo( fileName );