Skip to content
Permalink
Browse files

Nicer API for raster sampling

  • Loading branch information
nyalldawson committed Jul 16, 2018
1 parent 54e5119 commit 74c2ed12a5336e18b80bb82980e8ef3a21c1588e
@@ -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`

@@ -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()
@@ -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
@@ -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 );
}
@@ -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;
@@ -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 )

0 comments on commit 74c2ed1

Please sign in to comment.
You can’t perform that action at this time.