mirror of
https://github.com/qgis/QGIS.git
synced 2025-04-17 00:04:02 -04:00
use QgsRasterBlock instead of GDAL in zonal statistics (addresses #8736)
This should make zonal statistics usable rasters which come from other providers, e.g. WCS.
This commit is contained in:
Alexander Bruy
parent
290758a4bc
commit
278913b7ce
@ -30,7 +30,7 @@ class QgsZonalStatistics
|
||||
|
||||
typedef QFlags<QgsZonalStatistics::Statistic> Statistics;
|
||||
|
||||
QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix = "", int rasterBand = 1,
|
||||
QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix = "", int rasterBand = 1,
|
||||
QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistics( QgsZonalStatistics::Count | QgsZonalStatistics::Sum | QgsZonalStatistics::Mean) );
|
||||
|
||||
/** Starts the calculation
|
||||
|
||||
@ -100,6 +100,7 @@ class ZonalStatisticsQgis(GeoAlgorithm):
|
||||
st = self.getParameterValue(self.STATISTICS)
|
||||
|
||||
vectorLayer = dataobjects.getObjectFromUri(vectorPath)
|
||||
rasterLayer = dataobjects.getObjectFromUri(rasterPath)
|
||||
|
||||
keys = list(self.STATS.keys())
|
||||
selectedStats = 0
|
||||
@ -107,7 +108,7 @@ class ZonalStatisticsQgis(GeoAlgorithm):
|
||||
selectedStats |= self.STATS[keys[i]]
|
||||
|
||||
zs = QgsZonalStatistics(vectorLayer,
|
||||
rasterPath,
|
||||
rasterLayer,
|
||||
columnPrefix,
|
||||
bandNumber,
|
||||
selectedStats)
|
||||
|
||||
@ -20,16 +20,17 @@
|
||||
#include "qgsgeometry.h"
|
||||
#include "qgsvectordataprovider.h"
|
||||
#include "qgsvectorlayer.h"
|
||||
#include "qgsrasterdataprovider.h"
|
||||
#include "qgsrasterlayer.h"
|
||||
#include "qgsrasterblock.h"
|
||||
#include "qmath.h"
|
||||
#include "gdal.h"
|
||||
#include "cpl_string.h"
|
||||
#include "qgslogger.h"
|
||||
|
||||
#include <QProgressDialog>
|
||||
#include <QFile>
|
||||
|
||||
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand, Statistics stats )
|
||||
: mRasterFilePath( rasterFile )
|
||||
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix, int rasterBand, Statistics stats )
|
||||
: mRasterLayer( rasterLayer )
|
||||
, mRasterBand( rasterBand )
|
||||
, mPolygonLayer( polygonLayer )
|
||||
, mAttributePrefix( attributePrefix )
|
||||
@ -40,7 +41,8 @@ QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QStr
|
||||
}
|
||||
|
||||
QgsZonalStatistics::QgsZonalStatistics()
|
||||
: mRasterBand( 0 )
|
||||
: mRasterLayer( nullptr )
|
||||
, mRasterBand( 0 )
|
||||
, mPolygonLayer( nullptr )
|
||||
, mInputNodataValue( -1 )
|
||||
, mStatistics( QgsZonalStatistics::All )
|
||||
@ -61,49 +63,33 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
|
||||
return 2;
|
||||
}
|
||||
|
||||
//open the raster layer and the raster band
|
||||
GDALAllRegister();
|
||||
GDALDatasetH inputDataset = GDALOpen( mRasterFilePath.toUtf8().constData(), GA_ReadOnly );
|
||||
if ( !inputDataset )
|
||||
if ( !mRasterLayer )
|
||||
{
|
||||
return 3;
|
||||
}
|
||||
|
||||
if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
|
||||
if ( mRasterLayer->bandCount() < mRasterBand )
|
||||
{
|
||||
GDALClose( inputDataset );
|
||||
return 4;
|
||||
}
|
||||
|
||||
GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
|
||||
if ( !rasterBand )
|
||||
{
|
||||
GDALClose( inputDataset );
|
||||
return 5;
|
||||
}
|
||||
mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
|
||||
mRasterProvider = mRasterLayer->dataProvider();
|
||||
mInputNodataValue = mRasterProvider->sourceNoDataValue( mRasterBand );
|
||||
|
||||
//get geometry info about raster layer
|
||||
int nCellsXGDAL = GDALGetRasterXSize( inputDataset );
|
||||
int nCellsYGDAL = GDALGetRasterYSize( inputDataset );
|
||||
double geoTransform[6];
|
||||
if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
|
||||
{
|
||||
GDALClose( inputDataset );
|
||||
return 6;
|
||||
}
|
||||
double cellsizeX = geoTransform[1];
|
||||
int nCellsXProvider = mRasterProvider->xSize();
|
||||
int nCellsYProvider = mRasterProvider->ySize();
|
||||
double cellsizeX = mRasterLayer->rasterUnitsPerPixelX();
|
||||
if ( cellsizeX < 0 )
|
||||
{
|
||||
cellsizeX = -cellsizeX;
|
||||
}
|
||||
double cellsizeY = geoTransform[5];
|
||||
double cellsizeY = mRasterLayer->rasterUnitsPerPixelY();
|
||||
if ( cellsizeY < 0 )
|
||||
{
|
||||
cellsizeY = -cellsizeY;
|
||||
}
|
||||
QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsYGDAL * cellsizeY ),
|
||||
geoTransform[0] + ( nCellsXGDAL * cellsizeX ), geoTransform[3] );
|
||||
QgsRectangle rasterBBox = mRasterProvider->extent();
|
||||
|
||||
//add the new fields to the provider
|
||||
QList<QgsField> newFieldList;
|
||||
@ -223,7 +209,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
|
||||
p->setMaximum( featureCount );
|
||||
}
|
||||
|
||||
|
||||
//iterate over each polygon
|
||||
QgsFeatureRequest request;
|
||||
request.setSubsetOfAttributes( QgsAttributeList() );
|
||||
@ -273,22 +258,22 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
|
||||
}
|
||||
|
||||
//avoid access to cells outside of the raster (may occur because of rounding)
|
||||
if (( offsetX + nCellsX ) > nCellsXGDAL )
|
||||
if (( offsetX + nCellsX ) > nCellsXProvider )
|
||||
{
|
||||
nCellsX = nCellsXGDAL - offsetX;
|
||||
nCellsX = nCellsXProvider - offsetX;
|
||||
}
|
||||
if (( offsetY + nCellsY ) > nCellsYGDAL )
|
||||
if (( offsetY + nCellsY ) > nCellsYProvider )
|
||||
{
|
||||
nCellsY = nCellsYGDAL - offsetY;
|
||||
nCellsY = nCellsYProvider - offsetY;
|
||||
}
|
||||
|
||||
statisticsFromMiddlePointTest( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
|
||||
statisticsFromMiddlePointTest( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
|
||||
rasterBBox, featureStats );
|
||||
|
||||
if ( featureStats.count <= 1 )
|
||||
{
|
||||
//the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
|
||||
statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
|
||||
statisticsFromPreciseIntersection( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
|
||||
rasterBBox, featureStats );
|
||||
}
|
||||
|
||||
@ -366,7 +351,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
|
||||
p->setValue( featureCount );
|
||||
}
|
||||
|
||||
GDALClose( inputDataset );
|
||||
mPolygonLayer->updateFields();
|
||||
|
||||
if ( p && p->wasCanceled() )
|
||||
@ -404,12 +388,11 @@ int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const Q
|
||||
return 0;
|
||||
}
|
||||
|
||||
void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX,
|
||||
void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry& poly, int pixelOffsetX,
|
||||
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
|
||||
{
|
||||
double cellCenterX, cellCenterY;
|
||||
|
||||
float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
|
||||
cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
|
||||
stats.reset();
|
||||
|
||||
@ -430,17 +413,16 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeo
|
||||
GEOSCoordSequence* cellCenterCoords = nullptr;
|
||||
GEOSGeometry* currentCellCenter = nullptr;
|
||||
|
||||
QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
|
||||
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );
|
||||
|
||||
QgsRasterBlock* block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
|
||||
for ( int i = 0; i < nCellsY; ++i )
|
||||
{
|
||||
if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
|
||||
!= CPLE_None )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
|
||||
for ( int j = 0; j < nCellsX; ++j )
|
||||
{
|
||||
if ( validPixel( scanLine[j] ) )
|
||||
if ( validPixel( block->value( i, j ) ) )
|
||||
{
|
||||
GEOSGeom_destroy_r( geosctxt, currentCellCenter );
|
||||
cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
|
||||
@ -449,26 +431,26 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeo
|
||||
currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
|
||||
if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
|
||||
{
|
||||
stats.addValue( scanLine[j] );
|
||||
stats.addValue( block->value( i, j ) );
|
||||
}
|
||||
}
|
||||
cellCenterX += cellSizeX;
|
||||
}
|
||||
cellCenterY -= cellSizeY;
|
||||
}
|
||||
|
||||
GEOSGeom_destroy_r( geosctxt, currentCellCenter );
|
||||
CPLFree( scanLine );
|
||||
GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
|
||||
GEOSGeom_destroy_r( geosctxt, polyGeos );
|
||||
delete block;
|
||||
}
|
||||
|
||||
void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const QgsGeometry& poly, int pixelOffsetX,
|
||||
void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry& poly, int pixelOffsetX,
|
||||
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
|
||||
{
|
||||
stats.reset();
|
||||
|
||||
double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
|
||||
float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
|
||||
QgsGeometry pixelRectGeometry;
|
||||
|
||||
double hCellSizeX = cellSizeX / 2.0;
|
||||
@ -476,18 +458,19 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
|
||||
double pixelArea = cellSizeX * cellSizeY;
|
||||
double weight = 0;
|
||||
|
||||
for ( int row = 0; row < nCellsY; ++row )
|
||||
QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
|
||||
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );
|
||||
|
||||
QgsRasterBlock* block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
|
||||
for ( int i = 0; i < nCellsY; ++i )
|
||||
{
|
||||
double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
|
||||
for ( int col = 0; col < nCellsX; ++col )
|
||||
for ( int j = 0; j < nCellsX; ++j )
|
||||
{
|
||||
if ( GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 ) != CE_None )
|
||||
if ( !validPixel( block->value( i, j ) ) )
|
||||
{
|
||||
QgsDebugMsg( "Raster IO Error" );
|
||||
}
|
||||
|
||||
if ( !validPixel( *pixelData ) )
|
||||
continue;
|
||||
}
|
||||
|
||||
pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
|
||||
if ( !pixelRectGeometry.isEmpty() )
|
||||
@ -500,7 +483,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
|
||||
if ( intersectionArea >= 0.0 )
|
||||
{
|
||||
weight = intersectionArea / pixelArea;
|
||||
stats.addValue( *pixelData, weight );
|
||||
stats.addValue( block->value( i, j ), weight );
|
||||
}
|
||||
}
|
||||
pixelRectGeometry = QgsGeometry();
|
||||
@ -509,7 +492,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
|
||||
}
|
||||
currentY -= cellSizeY;
|
||||
}
|
||||
CPLFree( pixelData );
|
||||
delete block;
|
||||
}
|
||||
|
||||
bool QgsZonalStatistics::validPixel( float value ) const
|
||||
|
||||
@ -26,6 +26,8 @@
|
||||
|
||||
class QgsGeometry;
|
||||
class QgsVectorLayer;
|
||||
class QgsRasterLayer;
|
||||
class QgsRasterDataProvider;
|
||||
class QProgressDialog;
|
||||
class QgsRectangle;
|
||||
class QgsField;
|
||||
@ -57,7 +59,7 @@ class ANALYSIS_EXPORT QgsZonalStatistics
|
||||
/**
|
||||
* Constructor for QgsZonalStatistics.
|
||||
*/
|
||||
QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix = "", int rasterBand = 1,
|
||||
QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix = "", int rasterBand = 1,
|
||||
Statistics stats = Statistics( Count | Sum | Mean ) );
|
||||
|
||||
/** Starts the calculation
|
||||
@ -114,11 +116,11 @@ class ANALYSIS_EXPORT QgsZonalStatistics
|
||||
int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const;
|
||||
|
||||
//! Returns statistics by considering the pixels where the center point is within the polygon (fast)
|
||||
void statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
|
||||
void statisticsFromMiddlePointTest( const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
|
||||
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats& stats );
|
||||
|
||||
//! Returns statistics with precise pixel - polygon intersection test (slow)
|
||||
void statisticsFromPreciseIntersection( void* band, const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
|
||||
void statisticsFromPreciseIntersection( const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
|
||||
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats& stats );
|
||||
|
||||
//! Tests whether a pixel's value should be included in the result
|
||||
@ -126,7 +128,8 @@ class ANALYSIS_EXPORT QgsZonalStatistics
|
||||
|
||||
QString getUniqueFieldName( const QString& fieldName, const QList<QgsField>& newFields );
|
||||
|
||||
QString mRasterFilePath;
|
||||
QgsRasterLayer* mRasterLayer;
|
||||
QgsRasterDataProvider* mRasterProvider;
|
||||
//! Raster band to calculate statistics from (defaults to 1)
|
||||
int mRasterBand;
|
||||
QgsVectorLayer* mPolygonLayer;
|
||||
|
||||
@ -19,6 +19,7 @@
|
||||
#include "qgsapplication.h"
|
||||
#include "qgsfeatureiterator.h"
|
||||
#include "qgsvectorlayer.h"
|
||||
#include "qgsrasterlayer.h"
|
||||
#include "qgszonalstatistics.h"
|
||||
#include "qgsproject.h"
|
||||
|
||||
@ -42,7 +43,7 @@ class TestQgsZonalStatistics : public QObject
|
||||
|
||||
private:
|
||||
QgsVectorLayer* mVectorLayer;
|
||||
QString mRasterPath;
|
||||
QgsRasterLayer* mRasterLayer;
|
||||
};
|
||||
|
||||
TestQgsZonalStatistics::TestQgsZonalStatistics()
|
||||
@ -71,10 +72,9 @@ void TestQgsZonalStatistics::initTestCase()
|
||||
}
|
||||
|
||||
mVectorLayer = new QgsVectorLayer( myTempPath + "polys.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
|
||||
mRasterLayer = new QgsRasterLayer( myTempPath + "edge_problem.asc", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
|
||||
QgsProject::instance()->addMapLayers(
|
||||
QList<QgsMapLayer *>() << mVectorLayer );
|
||||
|
||||
mRasterPath = myTempPath + "edge_problem.asc";
|
||||
QList<QgsMapLayer *>() << mVectorLayer << mRasterLayer );
|
||||
}
|
||||
|
||||
void TestQgsZonalStatistics::cleanupTestCase()
|
||||
@ -84,7 +84,7 @@ void TestQgsZonalStatistics::cleanupTestCase()
|
||||
|
||||
void TestQgsZonalStatistics::testStatistics()
|
||||
{
|
||||
QgsZonalStatistics zs( mVectorLayer, mRasterPath, QLatin1String( "" ), 1, QgsZonalStatistics::All );
|
||||
QgsZonalStatistics zs( mVectorLayer, mRasterLayer, QLatin1String( "" ), 1, QgsZonalStatistics::All );
|
||||
zs.calculateStatistics( nullptr );
|
||||
|
||||
QgsFeature f;
|
||||
@ -135,7 +135,7 @@ void TestQgsZonalStatistics::testStatistics()
|
||||
QCOMPARE( f.attribute( "variety" ).toDouble(), 2.0 );
|
||||
|
||||
// same with long prefix to ensure that field name truncation handled correctly
|
||||
QgsZonalStatistics zsl( mVectorLayer, mRasterPath, QStringLiteral( "myqgis2_" ), 1, QgsZonalStatistics::All );
|
||||
QgsZonalStatistics zsl( mVectorLayer, mRasterLayer, QStringLiteral( "myqgis2_" ), 1, QgsZonalStatistics::All );
|
||||
zsl.calculateStatistics( nullptr );
|
||||
|
||||
request.setFilterFid( 0 );
|
||||
|
||||
@ -15,7 +15,7 @@ __revision__ = '$Format:%H$'
|
||||
import qgis # NOQA
|
||||
|
||||
from qgis.PyQt.QtCore import QDir, QFile
|
||||
from qgis.core import QgsVectorLayer, QgsFeature, QgsFeatureRequest
|
||||
from qgis.core import QgsVectorLayer, QgsRasterLayer, QgsFeature, QgsFeatureRequest
|
||||
from qgis.analysis import QgsZonalStatistics
|
||||
|
||||
from qgis.testing import start_app, unittest
|
||||
@ -38,8 +38,8 @@ class TestQgsZonalStatistics(unittest.TestCase):
|
||||
QFile.copy(TEST_DATA_DIR + f, myTempPath + f)
|
||||
|
||||
myVector = QgsVectorLayer(myTempPath + "polys.shp", "poly", "ogr")
|
||||
myRasterPath = myTempPath + "edge_problem.asc"
|
||||
zs = QgsZonalStatistics(myVector, myRasterPath, "", 1, QgsZonalStatistics.All)
|
||||
myRaster = QgsRasterLayer(myTempPath + "edge_problem.asc", "raster", "gdal")
|
||||
zs = QgsZonalStatistics(myVector, myRaster, "", 1, QgsZonalStatistics.All)
|
||||
zs.calculateStatistics(None)
|
||||
|
||||
feat = QgsFeature()
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user