23 #include "cpl_string.h"
24 #include <QProgressDialog>
26 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
27 #define TO8(x) (x).toUtf8().constData()
29 #define TO8(x) (x).toLocal8Bit().constData()
33 : mRasterFilePath( rasterFile )
34 , mRasterBand( rasterBand )
35 , mPolygonLayer( polygonLayer )
36 , mAttributePrefix( attributePrefix )
37 , mInputNodataValue( -1 )
62 if ( !vectorProvider )
70 if ( inputDataset == NULL )
75 if ( GDALGetRasterCount( inputDataset ) < (
mRasterBand - 1 ) )
77 GDALClose( inputDataset );
81 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset,
mRasterBand );
82 if ( rasterBand == NULL )
84 GDALClose( inputDataset );
90 int nCellsX = GDALGetRasterXSize( inputDataset );
91 int nCellsY = GDALGetRasterYSize( inputDataset );
92 double geoTransform[6];
93 if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
95 GDALClose( inputDataset );
98 double cellsizeX = geoTransform[1];
101 cellsizeX = -cellsizeX;
103 double cellsizeY = geoTransform[5];
106 cellsizeY = -cellsizeY;
108 QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsY * cellsizeY ), geoTransform[0] + ( nCellsX * cellsizeX ), geoTransform[3] );
111 QList<QgsField> newFieldList;
115 newFieldList.push_back( countField );
116 newFieldList.push_back( sumField );
117 newFieldList.push_back( meanField );
125 if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )
134 p->setMaximum( featureCount );
144 int featureCounter = 0;
148 qWarning(
"%d", featureCounter );
151 p->setValue( featureCounter );
154 if ( p && p->wasCanceled() )
160 if ( !featureGeometry )
166 int offsetX, offsetY, nCellsX, nCellsY;
167 if (
cellInfoForBBox( rasterBBox, featureGeometry->
boundingBox(), cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
174 rasterBBox, sum, count );
180 rasterBBox, sum, count );
196 changeAttributeMap.insert( countIndex, QVariant( count ) );
197 changeAttributeMap.insert( sumIndex, QVariant( sum ) );
198 changeAttributeMap.insert( meanIndex, QVariant( mean ) );
199 changeMap.insert( f.
id(), changeAttributeMap );
207 p->setValue( featureCount );
210 GDALClose( inputDataset );
213 if ( p && p->wasCanceled() )
222 int& offsetX,
int& offsetY,
int& nCellsX,
int& nCellsY )
const
228 nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
233 offsetX = ( int )(( intersectBox.
xMinimum() - rasterBBox.
xMinimum() ) / cellSizeX );
234 offsetY = ( int )(( rasterBBox.
yMaximum() - intersectBox.
yMaximum() ) / cellSizeY );
236 int maxColumn = ( int )(( intersectBox.
xMaximum() - rasterBBox.
xMinimum() ) / cellSizeX ) + 1;
237 int maxRow = ( int )(( rasterBBox.
yMaximum() - intersectBox.
yMinimum() ) / cellSizeY ) + 1;
239 nCellsX = maxColumn - offsetX;
240 nCellsY = maxRow - offsetY;
246 int pixelOffsetY,
int nCellsX,
int nCellsY,
double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox,
double& sum,
double& count )
248 double cellCenterX, cellCenterY;
251 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
252 cellCenterY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
256 for (
int i = 0; i < nCellsY; ++i )
258 GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
259 cellCenterX = rasterBBox.
xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
260 for (
int j = 0; j < nCellsX; ++j )
262 currentCellCenter =
QgsPoint( cellCenterX, cellCenterY );
263 if ( poly->
contains( ¤tCellCenter ) )
271 cellCenterX += cellSizeX;
273 cellCenterY -= cellSizeY;
279 int pixelOffsetY,
int nCellsX,
int nCellsY,
double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox,
double& sum,
double& count )
283 double currentY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
284 float* pixelData = (
float * ) CPLMalloc(
sizeof(
float ) );
287 double hCellSizeX = cellSizeX / 2.0;
288 double hCellSizeY = cellSizeY / 2.0;
289 double pixelArea = cellSizeX * cellSizeY;
292 for (
int row = 0; row < nCellsY; ++row )
294 double currentX = rasterBBox.
xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
295 for (
int col = 0; col < nCellsX; ++col )
297 GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
299 if ( pixelRectGeometry )
303 if ( intersectGeometry )
305 double intersectionArea = intersectGeometry->
area();
306 if ( intersectionArea >= 0.0 )
308 weight = intersectionArea / pixelArea;
310 sum += *pixelData * weight;
312 delete intersectGeometry;
315 currentX += cellSizeX;
317 currentY -= cellSizeY;
319 CPLFree( pixelData );
323 double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox,
double& sum,
double& count )
325 double cellCenterX, cellCenterY;
328 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
329 cellCenterY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
333 for (
int i = 0; i < nCellsY; ++i )
335 GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
336 cellCenterX = rasterBBox.
xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
339 GEOSCoordSequence* scanLineSequence = GEOSCoordSeq_create( 2, 2 );
340 GEOSCoordSeq_setX( scanLineSequence, 0, cellCenterX );
341 GEOSCoordSeq_setY( scanLineSequence, 0, cellCenterY );
342 GEOSCoordSeq_setX( scanLineSequence, 1, cellCenterX + nCellsX * cellSizeX );
343 GEOSCoordSeq_setY( scanLineSequence, 1, cellCenterY );
344 GEOSGeometry* scanLineGeos = GEOSGeom_createLineString( scanLineSequence );
345 GEOSGeometry* polyGeos = poly->
asGeos();
346 GEOSGeometry* scanLineIntersection = GEOSIntersection( scanLineGeos, polyGeos );
347 GEOSGeom_destroy( scanLineGeos );
348 if ( !scanLineIntersection )
350 cellCenterY -= cellSizeY;
357 int numGeoms = GEOSGetNumGeometries( scanLineIntersection );
360 GEOSGeom_destroy( scanLineIntersection );
361 cellCenterY -= cellSizeY;
365 QList<double> scanLineList;
367 GEOSGeometry* currentGeom = 0;
368 for (
int z = 0; z < numGeoms; ++z )
372 currentGeom = scanLineIntersection;
376 currentGeom = GEOSGeom_clone( GEOSGetGeometryN( scanLineIntersection, z ) );
378 const GEOSCoordSequence* scanLineCoordSequence = GEOSGeom_getCoordSeq( currentGeom );
379 if ( !scanLineCoordSequence )
383 unsigned int scanLineIntersectionSize;
384 GEOSCoordSeq_getSize( scanLineCoordSequence, &scanLineIntersectionSize );
385 if ( !scanLineCoordSequence || scanLineIntersectionSize < 2 || ( scanLineIntersectionSize & 1 ) )
389 for (
unsigned int k = 0; k < scanLineIntersectionSize; ++k )
391 GEOSCoordSeq_getX( scanLineCoordSequence, k, ¤tValue );
392 scanLineList.push_back( currentValue );
397 GEOSGeom_destroy( currentGeom );
400 GEOSGeom_destroy( scanLineIntersection );
401 qSort( scanLineList );
403 if ( scanLineList.size() < 1 )
405 cellCenterY -= cellSizeY;
410 for (
int j = 0; j < nCellsX; ++j )
415 if ( listPlace >= scanLineList.size() - 1 )
419 if ( cellCenterX >= scanLineList.at( listPlace + 1 ) )
422 if ( listPlace >= scanLineList.size() )
427 if ( listPlace >= 0 && listPlace < ( scanLineList.size() - 1 ) && !( listPlace & 1 ) )
435 cellCenterX += cellSizeX;
437 cellCenterY -= cellSizeY;