00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "qgszonalstatistics.h"
00019 #include "qgsgeometry.h"
00020 #include "qgsvectordataprovider.h"
00021 #include "qgsvectorlayer.h"
00022 #include "gdal.h"
00023 #include "cpl_string.h"
00024 #include <QProgressDialog>
00025
00026 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
00027 #define TO8(x) (x).toUtf8().constData()
00028 #else
00029 #define TO8(x) (x).toLocal8Bit().constData()
00030 #endif
00031
00032 QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand )
00033 : mRasterFilePath( rasterFile )
00034 , mRasterBand( rasterBand )
00035 , mPolygonLayer( polygonLayer )
00036 , mAttributePrefix( attributePrefix )
00037 , mInputNodataValue( -1 )
00038 {
00039
00040 }
00041
00042 QgsZonalStatistics::QgsZonalStatistics()
00043 : mRasterBand( 0 )
00044 , mPolygonLayer( 0 )
00045 {
00046
00047 }
00048
00049 QgsZonalStatistics::~QgsZonalStatistics()
00050 {
00051
00052 }
00053
00054 int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
00055 {
00056 if ( !mPolygonLayer || mPolygonLayer->geometryType() != QGis::Polygon )
00057 {
00058 return 1;
00059 }
00060
00061 QgsVectorDataProvider* vectorProvider = mPolygonLayer->dataProvider();
00062 if ( !vectorProvider )
00063 {
00064 return 2;
00065 }
00066
00067
00068 GDALAllRegister();
00069 GDALDatasetH inputDataset = GDALOpen( TO8( mRasterFilePath ), GA_ReadOnly );
00070 if ( inputDataset == NULL )
00071 {
00072 return 3;
00073 }
00074
00075 if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
00076 {
00077 GDALClose( inputDataset );
00078 return 4;
00079 }
00080
00081 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
00082 if ( rasterBand == NULL )
00083 {
00084 GDALClose( inputDataset );
00085 return 5;
00086 }
00087 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
00088
00089
00090 int nCellsX = GDALGetRasterXSize( inputDataset );
00091 int nCellsY = GDALGetRasterYSize( inputDataset );
00092 double geoTransform[6];
00093 if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
00094 {
00095 GDALClose( inputDataset );
00096 return 6;
00097 }
00098 double cellsizeX = geoTransform[1];
00099 if ( cellsizeX < 0 )
00100 {
00101 cellsizeX = -cellsizeX;
00102 }
00103 double cellsizeY = geoTransform[5];
00104 if ( cellsizeY < 0 )
00105 {
00106 cellsizeY = -cellsizeY;
00107 }
00108 QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsY * cellsizeY ), geoTransform[0] + ( nCellsX * cellsizeX ), geoTransform[3] );
00109
00110
00111 QList<QgsField> newFieldList;
00112 QgsField countField( mAttributePrefix + "count", QVariant::Double );
00113 QgsField sumField( mAttributePrefix + "sum", QVariant::Double );
00114 QgsField meanField( mAttributePrefix + "mean", QVariant::Double );
00115 newFieldList.push_back( countField );
00116 newFieldList.push_back( sumField );
00117 newFieldList.push_back( meanField );
00118 if ( !vectorProvider->addAttributes( newFieldList ) )
00119 {
00120 return 7;
00121 }
00122
00123
00124 int countIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "count" );
00125 int sumIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "sum" );
00126 int meanIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "mean" );
00127
00128 if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )
00129 {
00130 return 8;
00131 }
00132
00133
00134 long featureCount = vectorProvider->featureCount();
00135 if ( p )
00136 {
00137 p->setMaximum( featureCount );
00138 }
00139
00140
00141
00142 vectorProvider->select( QgsAttributeList(), QgsRectangle(), true, false );
00143 vectorProvider->rewind();
00144 QgsFeature f;
00145 double count = 0;
00146 double sum = 0;
00147 double mean = 0;
00148 int featureCounter = 0;
00149
00150 while ( vectorProvider->nextFeature( f ) )
00151 {
00152 qWarning( "%d", featureCounter );
00153 if ( p )
00154 {
00155 p->setValue( featureCounter );
00156 }
00157
00158 if ( p && p->wasCanceled() )
00159 {
00160 break;
00161 }
00162
00163 QgsGeometry* featureGeometry = f.geometry();
00164 if ( !featureGeometry )
00165 {
00166 ++featureCounter;
00167 continue;
00168 }
00169
00170 int offsetX, offsetY, nCellsX, nCellsY;
00171 if ( cellInfoForBBox( rasterBBox, featureGeometry->boundingBox(), cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
00172 {
00173 ++featureCounter;
00174 continue;
00175 }
00176
00177 statisticsFromMiddlePointTest_improved( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
00178 rasterBBox, sum, count );
00179
00180 if ( count <= 1 )
00181 {
00182
00183 statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
00184 rasterBBox, sum, count );
00185 }
00186
00187
00188 if ( count == 0 )
00189 {
00190 mean = 0;
00191 }
00192 else
00193 {
00194 mean = sum / count;
00195 }
00196
00197
00198 QgsChangedAttributesMap changeMap;
00199 QgsAttributeMap changeAttributeMap;
00200 changeAttributeMap.insert( countIndex, QVariant( count ) );
00201 changeAttributeMap.insert( sumIndex, QVariant( sum ) );
00202 changeAttributeMap.insert( meanIndex, QVariant( mean ) );
00203 changeMap.insert( f.id(), changeAttributeMap );
00204 vectorProvider->changeAttributeValues( changeMap );
00205
00206 ++featureCounter;
00207 }
00208
00209 if ( p )
00210 {
00211 p->setValue( featureCount );
00212 }
00213
00214 GDALClose( inputDataset );
00215 return 0;
00216 }
00217
00218 int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY,
00219 int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const
00220 {
00221
00222 QgsRectangle intersectBox = rasterBBox.intersect( &featureBBox );
00223 if ( intersectBox.isEmpty() )
00224 {
00225 nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
00226 return 0;
00227 }
00228
00229
00230 offsetX = ( int )(( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX );
00231 offsetY = ( int )(( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY );
00232
00233 int maxColumn = ( int )(( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) + 1;
00234 int maxRow = ( int )(( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) + 1;
00235
00236 nCellsX = maxColumn - offsetX;
00237 nCellsY = maxRow - offsetY;
00238
00239 return 0;
00240 }
00241
00242 void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX,
00243 int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
00244 {
00245 double cellCenterX, cellCenterY;
00246 QgsPoint currentCellCenter;
00247
00248 float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
00249 cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
00250 count = 0;
00251 sum = 0;
00252
00253 for ( int i = 0; i < nCellsY; ++i )
00254 {
00255 GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
00256 cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
00257 for ( int j = 0; j < nCellsX; ++j )
00258 {
00259 currentCellCenter = QgsPoint( cellCenterX, cellCenterY );
00260 if ( poly->contains( ¤tCellCenter ) )
00261 {
00262 if ( scanLine[j] != mInputNodataValue )
00263 {
00264 sum += scanLine[j];
00265 ++count;
00266 }
00267 }
00268 cellCenterX += cellSizeX;
00269 }
00270 cellCenterY -= cellSizeY;
00271 }
00272 CPLFree( scanLine );
00273 }
00274
00275 void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX,
00276 int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
00277 {
00278 sum = 0;
00279 count = 0;
00280 double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
00281 float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
00282 QgsGeometry* pixelRectGeometry = 0;
00283
00284 double hCellSizeX = cellSizeX / 2.0;
00285 double hCellSizeY = cellSizeY / 2.0;
00286 double pixelArea = cellSizeX * cellSizeY;
00287 double weight = 0;
00288
00289 for ( int row = 0; row < nCellsY; ++row )
00290 {
00291 double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
00292 for ( int col = 0; col < nCellsX; ++col )
00293 {
00294 GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
00295 pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
00296 if ( pixelRectGeometry )
00297 {
00298
00299 QgsGeometry *intersectGeometry = pixelRectGeometry->intersection( poly );
00300 if ( intersectGeometry )
00301 {
00302 double intersectionArea = intersectGeometry->area();
00303 if ( intersectionArea >= 0.0 )
00304 {
00305 weight = intersectionArea / pixelArea;
00306 count += weight;
00307 sum += *pixelData * weight;
00308 }
00309 delete intersectGeometry;
00310 }
00311 }
00312 currentX += cellSizeX;
00313 }
00314 currentY -= cellSizeY;
00315 }
00316 CPLFree( pixelData );
00317 }
00318
00319 void QgsZonalStatistics::statisticsFromMiddlePointTest_improved( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
00320 double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
00321 {
00322 double cellCenterX, cellCenterY;
00323 QgsPoint currentCellCenter;
00324
00325 float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
00326 cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
00327 count = 0;
00328 sum = 0;
00329
00330 for ( int i = 0; i < nCellsY; ++i )
00331 {
00332 GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
00333 cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
00334
00335
00336 GEOSCoordSequence* scanLineSequence = GEOSCoordSeq_create( 2, 2 );
00337 GEOSCoordSeq_setX( scanLineSequence, 0, cellCenterX );
00338 GEOSCoordSeq_setY( scanLineSequence, 0, cellCenterY );
00339 GEOSCoordSeq_setX( scanLineSequence, 1, cellCenterX + nCellsX * cellSizeX );
00340 GEOSCoordSeq_setY( scanLineSequence, 1, cellCenterY );
00341 GEOSGeometry* scanLineGeos = GEOSGeom_createLineString( scanLineSequence );
00342 GEOSGeometry* polyGeos = poly->asGeos();
00343 GEOSGeometry* scanLineIntersection = GEOSIntersection( scanLineGeos, polyGeos );
00344 GEOSGeom_destroy( scanLineGeos );
00345 if ( !scanLineIntersection )
00346 {
00347 cellCenterY -= cellSizeY;
00348 continue;
00349 }
00350
00351
00352
00353
00354 int numGeoms = GEOSGetNumGeometries( scanLineIntersection );
00355 if ( numGeoms < 1 )
00356 {
00357 GEOSGeom_destroy( scanLineIntersection );
00358 cellCenterY -= cellSizeY;
00359 continue;
00360 }
00361
00362 QList<double> scanLineList;
00363 double currentValue;
00364 GEOSGeometry* currentGeom = 0;
00365 for ( int z = 0; z < numGeoms; ++z )
00366 {
00367 if ( numGeoms == 1 )
00368 {
00369 currentGeom = scanLineIntersection;
00370 }
00371 else
00372 {
00373 currentGeom = GEOSGeom_clone( GEOSGetGeometryN( scanLineIntersection, z ) );
00374 }
00375 const GEOSCoordSequence* scanLineCoordSequence = GEOSGeom_getCoordSeq( currentGeom );
00376 if ( !scanLineCoordSequence )
00377 {
00378
00379 }
00380 unsigned int scanLineIntersectionSize;
00381 GEOSCoordSeq_getSize( scanLineCoordSequence, &scanLineIntersectionSize );
00382 if ( !scanLineCoordSequence || scanLineIntersectionSize < 2 || ( scanLineIntersectionSize & 1 ) )
00383 {
00384
00385 }
00386 for ( unsigned int k = 0; k < scanLineIntersectionSize; ++k )
00387 {
00388 GEOSCoordSeq_getX( scanLineCoordSequence, k, ¤tValue );
00389 scanLineList.push_back( currentValue );
00390 }
00391
00392 if ( numGeoms != 1 )
00393 {
00394 GEOSGeom_destroy( currentGeom );
00395 }
00396 }
00397 GEOSGeom_destroy( scanLineIntersection );
00398 qSort( scanLineList );
00399
00400 if ( scanLineList.size() < 1 )
00401 {
00402 cellCenterY -= cellSizeY;
00403 continue;
00404 }
00405
00406 int listPlace = -1;
00407 for ( int j = 0; j < nCellsX; ++j )
00408 {
00409
00410
00411
00412 if ( listPlace >= scanLineList.size() - 1 )
00413 {
00414 break;
00415 }
00416 if ( cellCenterX >= scanLineList.at( listPlace + 1 ) )
00417 {
00418 ++listPlace;
00419 if ( listPlace >= scanLineList.size() )
00420 {
00421 break;
00422 }
00423 }
00424 if ( listPlace >= 0 && listPlace < ( scanLineList.size() - 1 ) && !( listPlace & 1 ) )
00425 {
00426 if ( scanLine[j] != mInputNodataValue )
00427 {
00428 sum += scanLine[j];
00429 ++count;
00430 }
00431 }
00432 cellCenterX += cellSizeX;
00433 }
00434 cellCenterY -= cellSizeY;
00435 }
00436 CPLFree( scanLine );
00437 }
00438
00439