23 #include "cpl_string.h"
24 #include <QProgressDialog>
27 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
28 #define TO8F(x) (x).toUtf8().constData()
30 #define TO8F(x) QFile::encodeName( x ).constData()
34 : mRasterFilePath( rasterFile )
35 , mRasterBand( rasterBand )
36 , mPolygonLayer( polygonLayer )
37 , mAttributePrefix( attributePrefix )
38 , mInputNodataValue( -1 )
63 if ( !vectorProvider )
71 if ( inputDataset == NULL )
76 if ( GDALGetRasterCount( inputDataset ) < (
mRasterBand - 1 ) )
78 GDALClose( inputDataset );
82 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset,
mRasterBand );
83 if ( rasterBand == NULL )
85 GDALClose( inputDataset );
91 int nCellsXGDAL = GDALGetRasterXSize( inputDataset );
92 int nCellsYGDAL = GDALGetRasterYSize( inputDataset );
93 double geoTransform[6];
94 if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
96 GDALClose( inputDataset );
99 double cellsizeX = geoTransform[1];
102 cellsizeX = -cellsizeX;
104 double cellsizeY = geoTransform[5];
107 cellsizeY = -cellsizeY;
109 QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsYGDAL * cellsizeY ),
110 geoTransform[0] + ( nCellsXGDAL * cellsizeX ), geoTransform[3] );
113 QList<QgsField> newFieldList;
117 QgsField countField( countFieldName, QVariant::Double,
"double precision" );
118 QgsField sumField( sumFieldName, QVariant::Double,
"double precision" );
119 QgsField meanField( meanFieldName, QVariant::Double,
"double precision" );
120 newFieldList.push_back( countField );
121 newFieldList.push_back( sumField );
122 newFieldList.push_back( meanField );
126 int countIndex = vectorProvider->
fieldNameIndex( countFieldName );
130 if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )
139 p->setMaximum( featureCount );
151 int featureCounter = 0;
157 p->setValue( featureCounter );
160 if ( p && p->wasCanceled() )
166 if ( !featureGeometry )
179 int offsetX, offsetY, nCellsX, nCellsY;
180 if (
cellInfoForBBox( rasterBBox, featureRect, cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
187 if (( offsetX + nCellsX ) > nCellsXGDAL )
189 nCellsX = nCellsXGDAL - offsetX;
191 if (( offsetY + nCellsY ) > nCellsYGDAL )
193 nCellsY = nCellsYGDAL - offsetY;
197 rasterBBox, sum, count );
203 rasterBBox, sum, count );
219 changeAttributeMap.insert( countIndex, QVariant( count ) );
220 changeAttributeMap.insert( sumIndex, QVariant( sum ) );
221 changeAttributeMap.insert( meanIndex, QVariant( mean ) );
222 changeMap.insert( f.
id(), changeAttributeMap );
230 p->setValue( featureCount );
233 GDALClose( inputDataset );
236 if ( p && p->wasCanceled() )
245 int& offsetX,
int& offsetY,
int& nCellsX,
int& nCellsY )
const
251 nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
256 offsetX = ( int )(( intersectBox.
xMinimum() - rasterBBox.
xMinimum() ) / cellSizeX );
257 offsetY = ( int )(( rasterBBox.
yMaximum() - intersectBox.
yMaximum() ) / cellSizeY );
259 int maxColumn = ( int )(( intersectBox.
xMaximum() - rasterBBox.
xMinimum() ) / cellSizeX ) + 1;
260 int maxRow = ( int )(( rasterBBox.
yMaximum() - intersectBox.
yMinimum() ) / cellSizeY ) + 1;
262 nCellsX = maxColumn - offsetX;
263 nCellsY = maxRow - offsetY;
269 int pixelOffsetY,
int nCellsX,
int nCellsY,
double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox,
double& sum,
double& count )
271 double cellCenterX, cellCenterY;
273 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
274 cellCenterY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
278 const GEOSGeometry* polyGeos = poly->
asGeos();
284 const GEOSPreparedGeometry* polyGeosPrepared = GEOSPrepare( poly->
asGeos() );
285 if ( !polyGeosPrepared )
290 GEOSCoordSequence* cellCenterCoords = 0;
291 GEOSGeometry* currentCellCenter = 0;
293 for (
int i = 0; i < nCellsY; ++i )
295 if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
300 cellCenterX = rasterBBox.
xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
301 for (
int j = 0; j < nCellsX; ++j )
303 GEOSGeom_destroy( currentCellCenter );
304 cellCenterCoords = GEOSCoordSeq_create( 1, 2 );
305 GEOSCoordSeq_setX( cellCenterCoords, 0, cellCenterX );
306 GEOSCoordSeq_setY( cellCenterCoords, 0, cellCenterY );
307 currentCellCenter = GEOSGeom_createPoint( cellCenterCoords );
311 if ( GEOSPreparedContains( polyGeosPrepared, currentCellCenter ) )
313 if ( !qIsNaN( scanLine[j] ) )
320 cellCenterX += cellSizeX;
322 cellCenterY -= cellSizeY;
325 GEOSPreparedGeom_destroy( polyGeosPrepared );
329 int pixelOffsetY,
int nCellsX,
int nCellsY,
double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox,
double& sum,
double& count )
333 double currentY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
334 float* pixelData = (
float * ) CPLMalloc(
sizeof(
float ) );
337 double hCellSizeX = cellSizeX / 2.0;
338 double hCellSizeY = cellSizeY / 2.0;
339 double pixelArea = cellSizeX * cellSizeY;
342 for (
int row = 0; row < nCellsY; ++row )
344 double currentX = rasterBBox.
xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
345 for (
int col = 0; col < nCellsX; ++col )
347 GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
349 if ( pixelRectGeometry )
353 if ( intersectGeometry )
355 double intersectionArea = intersectGeometry->
area();
356 if ( intersectionArea >= 0.0 )
358 weight = intersectionArea / pixelArea;
360 sum += *pixelData * weight;
362 delete intersectGeometry;
365 currentX += cellSizeX;
367 currentY -= cellSizeY;
369 CPLFree( pixelData );
376 if ( !dp->
storageType().contains(
"ESRI Shapefile" ) )
382 QString shortName = fieldName.mid( 0, 10 );
385 for (
int idx = 0; idx < providerFields.
count(); ++idx )
387 if ( shortName == providerFields[idx].name() )
400 shortName = QString(
"%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
405 for (
int idx = 0; idx < providerFields.
count(); ++idx )
407 if ( shortName == providerFields[idx].name() )
412 shortName = QString(
"%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
416 shortName = QString(
"%1_%2" ).arg( fieldName.mid( 0, 7 ) ).arg( n );
QgsFeatureId id() const
Get the feature id for this feature.
void updateFields()
Assembles mUpdatedFields considering provider fields, joined fields and added fields.
Wrapper for iterator of features from vector data provider or vector layer.
A rectangle specified with double values.
bool isEmpty() const
test if rectangle is empty.
QMap< int, QVariant > QgsAttributeMap
virtual bool addAttributes(const QList< QgsField > &attributes)
Adds new attributes.
double yMaximum() const
Get the y maximum value (top side of rectangle)
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
Container of fields for a vector layer.
void statisticsFromMiddlePointTest(void *band, QgsGeometry *poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, double &sum, double &count)
Returns statistics by considering the pixels where the center point is within the polygon (fast) ...
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
int fieldNameIndex(const QString &fieldName) const
Returns the index of a field name or -1 if the field does not exist.
double area()
get area of geometry using GEOS
void statisticsFromPreciseIntersection(void *band, QgsGeometry *poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, double &sum, double &count)
Returns statistics with precise pixel - polygon intersection test (slow)
QString getUniqueFieldName(QString fieldName)
const GEOSGeometry * asGeos() const
Returns a geos geomtry.
int mRasterBand
Raster band to calculate statistics from (defaults to 1)
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
double xMaximum() const
Get the x maximum value (right side of rectangle)
virtual bool changeAttributeValues(const QgsChangedAttributesMap &attr_map)
Changes attribute values of existing features.
virtual QString storageType() const
Returns the permanent storage type for this layer as a friendly name.
virtual long featureCount() const =0
Number of features in the layer.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QList< int > QgsAttributeList
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())=0
Query the provider for features specified in request.
int count() const
Return number of items.
QGis::GeometryType geometryType() const
Returns point, line or polygon.
Encapsulate a field in an attribute table or data source.
QgsGeometry * intersection(QgsGeometry *geometry)
Returns a geometry representing the points shared by this geometry and other.
int cellInfoForBBox(const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY, int &offsetX, int &offsetY, int &nCellsX, int &nCellsY) const
Analysis what cells need to be considered to cover the bounding box of a feature. ...
QgsRectangle boundingBox()
Returns the bounding box of this feature.
virtual const QgsFields & fields() const =0
Return a map of indexes with field names for this layer.
int calculateStatistics(QProgressDialog *p)
Starts the calculation.
QMap< QgsFeatureId, QgsAttributeMap > QgsChangedAttributesMap
QgsVectorLayer * mPolygonLayer
QgsRectangle intersect(const QgsRectangle *rect) const
return the intersection with the given rectangle
static QgsGeometry * fromRect(const QgsRectangle &rect)
construct geometry from a rectangle
float mInputNodataValue
The nodata value of the input layer.
QgsVectorDataProvider * dataProvider()
Returns the data provider.
bool nextFeature(QgsFeature &f)
This is the base class for vector data providers.
Represents a vector layer which manages a vector based data sets.
double xMinimum() const
Get the x minimum value (left side of rectangle)