Nicer API for raster sampling

This commit is contained in:
Nyall Dawson 2018-07-16 19:28:08 +10:00
parent 54e511960c
commit 74c2ed12a5
6 changed files with 55 additions and 28 deletions

View File

@ -277,17 +277,19 @@ if point is outside data source extent.
.. 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 );
virtual double sample( const QgsPointXY &point, int band,
bool *ok /Out/ = 0,
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.
resolution is used.
A null QVariant will be returned if the point is outside data source extent, or an invalid
band number was specified.
If ``ok`` is specified and the point is outside data source extent, or an invalid
band number was specified, then ``ok`` will be set to false. In this case the function will return
a NaN value.
.. seealso:: :py:func:`identify`

View File

@ -300,12 +300,16 @@ 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 )
double QgsRasterDataProvider::sample( const QgsPointXY &point, int band,
bool *ok, const QgsRectangle &boundingBox, int width, int height, int )
{
if ( ok )
*ok = false;
if ( !extent().contains( point ) )
{
// Outside the raster
return QVariant();
return std::numeric_limits<double>::quiet_NaN();
}
QgsRectangle finalExtent = boundingBox;
@ -335,8 +339,10 @@ QVariant QgsRasterDataProvider::sample( const QgsPointXY &point, int band, const
QgsRectangle pixelExtent( xMin, yMin, xMax, yMax );
std::unique_ptr< QgsRasterBlock > bandBlock( block( band, pixelExtent, 1, 1 ) );
if ( bandBlock && ok )
*ok = true;
return bandBlock ? bandBlock->value( 0 ) : QVariant();
return bandBlock ? bandBlock->value( 0 ) : std::numeric_limits<double>::quiet_NaN();
}
QString QgsRasterDataProvider::lastErrorFormat()

View File

@ -364,16 +364,18 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
* 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.
* resolution is used.
*
* A null QVariant will be returned if the point is outside data source extent, or an invalid
* band number was specified.
* If \a ok is specified and the point is outside data source extent, or an invalid
* band number was specified, then \a ok will be set to false. In this case the function will return
* a NaN value.
*
* \see identify(), which is much more flexible but considerably less efficient.
* \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 );
virtual double sample( const QgsPointXY &point, int band,
bool *ok SIP_OUT = nullptr,
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

View File

@ -1101,32 +1101,38 @@ bool QgsGdalProvider::worldToPixel( double x, double y, int &col, int &row ) con
return true;
};
QVariant QgsGdalProvider::sample( const QgsPointXY &point, int band, const QgsRectangle &, int, int, int )
double QgsGdalProvider::sample( const QgsPointXY &point, int band, bool *ok, const QgsRectangle &, int, int, int )
{
if ( ok )
*ok = false;
QMutexLocker locker( mpMutex );
if ( !initIfNeeded() )
return tr( "Cannot read data" );
return std::numeric_limits<double>::quiet_NaN();
if ( !extent().contains( point ) )
{
// Outside the raster
return QVariant();
return std::numeric_limits<double>::quiet_NaN();
}
GDALRasterBandH hBand = GDALGetRasterBand( mGdalDataset, band );
if ( !hBand )
return QVariant();
return std::numeric_limits<double>::quiet_NaN();
int row;
int col;
if ( !worldToPixel( point.x(), point.y(), col, row ) )
return QVariant();
return std::numeric_limits<double>::quiet_NaN();
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();
return std::numeric_limits<double>::quiet_NaN();
if ( ok )
*ok = true;
return static_cast< double >( value ) * bandScale( band ) + bandOffset( band );
}

View File

@ -103,7 +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 );
double sample( const QgsPointXY &point, int band, bool *ok = nullptr, const QgsRectangle &boundingBox = QgsRectangle(), int width = 0, int height = 0, int dpi = 96 );
QString lastErrorTitle() override;
QString lastError() override;
int capabilities() const override;

View File

@ -729,21 +729,32 @@ void TestQgsRasterLayer::sample()
std::unique_ptr< QgsRasterLayer > rl = qgis::make_unique< QgsRasterLayer> ( rasterFileInfo.filePath(),
rasterFileInfo.completeBaseName() );
QVERIFY( rl->isValid() );
QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1 ).isValid() );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 1 ).toInt(), 125 );
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1 ) ) );
bool ok = false;
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1, &ok ) ) );
QVERIFY( !ok );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 1 ), 125.0 );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 1, &ok ), 125.0 );
QVERIFY( ok );
// bad bands
QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 0 ).isValid() );
QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 10 ).isValid() );
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 0 ) ) );
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 0, &ok ) ) );
QVERIFY( !ok );
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 788461, 3344957 ), 10, &ok ) ) );
QVERIFY( !ok );
fileName = mTestDataDir + "landsat_4326.tif";
rasterFileInfo = QFileInfo( fileName );
rl = qgis::make_unique< QgsRasterLayer> ( rasterFileInfo.filePath(),
rasterFileInfo.completeBaseName() );
QVERIFY( rl->isValid() );
QVERIFY( !rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1 ).isValid() );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 17.943731, 30.230791 ), 1 ).toInt(), 125 );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 17.943731, 30.230791 ), 2 ).toInt(), 139 );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 17.943731, 30.230791 ), 3 ).toInt(), 111 );
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1 ) ) );
QVERIFY( std::isnan( rl->dataProvider()->sample( QgsPointXY( 0, 0 ), 1, &ok ) ) );
QVERIFY( !ok );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 17.943731, 30.230791 ), 1, &ok ), 125.0 );
QVERIFY( ok );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 17.943731, 30.230791 ), 2 ), 139.0 );
QCOMPARE( rl->dataProvider()->sample( QgsPointXY( 17.943731, 30.230791 ), 3 ), 111.0 );
}
QGSTEST_MAIN( TestQgsRasterLayer )