Skip to content

Commit

Permalink
More heavily optimised sample method for gdal provider
Browse files Browse the repository at this point in the history
  • Loading branch information
nyalldawson committed Jul 19, 2018
1 parent db8e536 commit 54e5119
Show file tree
Hide file tree
Showing 5 changed files with 34 additions and 54 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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`

Expand Down
3 changes: 2 additions & 1 deletion src/core/raster/qgsrasterdataprovider.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
78 changes: 26 additions & 52 deletions src/providers/gdal/qgsgdalprovider.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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() )
Expand All @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/providers/gdal/qgsgdalprovider.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 3 additions & 0 deletions tests/src/core/testqgsrasterlayer.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
Expand Down

0 comments on commit 54e5119

Please sign in to comment.