Quantum GIS API Documentation
1.7.4
|
00001 /*************************************************************************** 00002 qgszonalstatistics.cpp - description 00003 ---------------------------- 00004 begin : August 29th, 2009 00005 copyright : (C) 2009 by Marco Hugentobler 00006 email : marco at hugis dot net 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 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 //open the raster layer and the raster band 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 //get geometry info about raster layer 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 //add the new count, sum, mean fields to the provider 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 vectorProvider->addAttributes( newFieldList ); 00119 00120 //index of the new fields 00121 int countIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "count" ); 00122 int sumIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "sum" ); 00123 int meanIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "mean" ); 00124 00125 if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 ) 00126 { 00127 return 8; 00128 } 00129 00130 //progress dialog 00131 long featureCount = vectorProvider->featureCount(); 00132 if ( p ) 00133 { 00134 p->setMaximum( featureCount ); 00135 } 00136 00137 00138 //iterate over each polygon 00139 vectorProvider->select( QgsAttributeList(), QgsRectangle(), true, false ); 00140 QgsFeature f; 00141 double count = 0; 00142 double sum = 0; 00143 double mean = 0; 00144 int featureCounter = 0; 00145 00146 while ( vectorProvider->nextFeature( f ) ) 00147 { 00148 qWarning( "%d", featureCounter ); 00149 if ( p ) 00150 { 00151 p->setValue( featureCounter ); 00152 } 00153 00154 if ( p && p->wasCanceled() ) 00155 { 00156 break; 00157 } 00158 00159 QgsGeometry* featureGeometry = f.geometry(); 00160 if ( !featureGeometry ) 00161 { 00162 ++featureCounter; 00163 continue; 00164 } 00165 00166 int offsetX, offsetY, nCellsX, nCellsY; 00167 if ( cellInfoForBBox( rasterBBox, featureGeometry->boundingBox(), cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 ) 00168 { 00169 ++featureCounter; 00170 continue; 00171 } 00172 00173 statisticsFromMiddlePointTest_improved( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY, 00174 rasterBBox, sum, count ); 00175 00176 if ( count <= 1 ) 00177 { 00178 //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case 00179 statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY, 00180 rasterBBox, sum, count ); 00181 } 00182 00183 00184 if ( count == 0 ) 00185 { 00186 mean = 0; 00187 } 00188 else 00189 { 00190 mean = sum / count; 00191 } 00192 00193 //write the statistics value to the vector data provider 00194 QgsChangedAttributesMap changeMap; 00195 QgsAttributeMap changeAttributeMap; 00196 changeAttributeMap.insert( countIndex, QVariant( count ) ); 00197 changeAttributeMap.insert( sumIndex, QVariant( sum ) ); 00198 changeAttributeMap.insert( meanIndex, QVariant( mean ) ); 00199 changeMap.insert( f.id(), changeAttributeMap ); 00200 vectorProvider->changeAttributeValues( changeMap ); 00201 00202 ++featureCounter; 00203 } 00204 00205 if ( p ) 00206 { 00207 p->setValue( featureCount ); 00208 } 00209 00210 GDALClose( inputDataset ); 00211 mPolygonLayer->updateFieldMap(); 00212 00213 if ( p && p->wasCanceled() ) 00214 { 00215 return 9; 00216 } 00217 00218 return 0; 00219 } 00220 00221 int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY, 00222 int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const 00223 { 00224 //get intersecting bbox 00225 QgsRectangle intersectBox = rasterBBox.intersect( &featureBBox ); 00226 if ( intersectBox.isEmpty() ) 00227 { 00228 nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0; 00229 return 0; 00230 } 00231 00232 //get offset in pixels in x- and y- direction 00233 offsetX = ( int )(( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX ); 00234 offsetY = ( int )(( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY ); 00235 00236 int maxColumn = ( int )(( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) + 1; 00237 int maxRow = ( int )(( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) + 1; 00238 00239 nCellsX = maxColumn - offsetX; 00240 nCellsY = maxRow - offsetY; 00241 00242 return 0; 00243 } 00244 00245 void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX, 00246 int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count ) 00247 { 00248 double cellCenterX, cellCenterY; 00249 QgsPoint currentCellCenter; 00250 00251 float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX ); 00252 cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2; 00253 count = 0; 00254 sum = 0; 00255 00256 for ( int i = 0; i < nCellsY; ++i ) 00257 { 00258 GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 ); 00259 cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2; 00260 for ( int j = 0; j < nCellsX; ++j ) 00261 { 00262 currentCellCenter = QgsPoint( cellCenterX, cellCenterY ); 00263 if ( poly->contains( ¤tCellCenter ) ) 00264 { 00265 if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values 00266 { 00267 sum += scanLine[j]; 00268 ++count; 00269 } 00270 } 00271 cellCenterX += cellSizeX; 00272 } 00273 cellCenterY -= cellSizeY; 00274 } 00275 CPLFree( scanLine ); 00276 } 00277 00278 void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX, 00279 int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count ) 00280 { 00281 sum = 0; 00282 count = 0; 00283 double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2; 00284 float* pixelData = ( float * ) CPLMalloc( sizeof( float ) ); 00285 QgsGeometry* pixelRectGeometry = 0; 00286 00287 double hCellSizeX = cellSizeX / 2.0; 00288 double hCellSizeY = cellSizeY / 2.0; 00289 double pixelArea = cellSizeX * cellSizeY; 00290 double weight = 0; 00291 00292 for ( int row = 0; row < nCellsY; ++row ) 00293 { 00294 double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX; 00295 for ( int col = 0; col < nCellsX; ++col ) 00296 { 00297 GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 ); 00298 pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) ); 00299 if ( pixelRectGeometry ) 00300 { 00301 //intersection 00302 QgsGeometry *intersectGeometry = pixelRectGeometry->intersection( poly ); 00303 if ( intersectGeometry ) 00304 { 00305 double intersectionArea = intersectGeometry->area(); 00306 if ( intersectionArea >= 0.0 ) 00307 { 00308 weight = intersectionArea / pixelArea; 00309 count += weight; 00310 sum += *pixelData * weight; 00311 } 00312 delete intersectGeometry; 00313 } 00314 } 00315 currentX += cellSizeX; 00316 } 00317 currentY -= cellSizeY; 00318 } 00319 CPLFree( pixelData ); 00320 } 00321 00322 void QgsZonalStatistics::statisticsFromMiddlePointTest_improved( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, 00323 double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count ) 00324 { 00325 double cellCenterX, cellCenterY; 00326 QgsPoint currentCellCenter; 00327 00328 float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX ); 00329 cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2; 00330 count = 0; 00331 sum = 0; 00332 00333 for ( int i = 0; i < nCellsY; ++i ) 00334 { 00335 GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 ); 00336 cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2; 00337 00338 //do intersection of scanline with geometry 00339 GEOSCoordSequence* scanLineSequence = GEOSCoordSeq_create( 2, 2 ); 00340 GEOSCoordSeq_setX( scanLineSequence, 0, cellCenterX ); 00341 GEOSCoordSeq_setY( scanLineSequence, 0, cellCenterY ); 00342 GEOSCoordSeq_setX( scanLineSequence, 1, cellCenterX + nCellsX * cellSizeX ); 00343 GEOSCoordSeq_setY( scanLineSequence, 1, cellCenterY ); 00344 GEOSGeometry* scanLineGeos = GEOSGeom_createLineString( scanLineSequence ); //todo: delete 00345 GEOSGeometry* polyGeos = poly->asGeos(); 00346 GEOSGeometry* scanLineIntersection = GEOSIntersection( scanLineGeos, polyGeos ); 00347 GEOSGeom_destroy( scanLineGeos ); 00348 if ( !scanLineIntersection ) 00349 { 00350 cellCenterY -= cellSizeY; 00351 continue; 00352 } 00353 00354 //debug 00355 //char* scanLineIntersectionType = GEOSGeomType( scanLineIntersection ); 00356 00357 int numGeoms = GEOSGetNumGeometries( scanLineIntersection ); 00358 if ( numGeoms < 1 ) 00359 { 00360 GEOSGeom_destroy( scanLineIntersection ); 00361 cellCenterY -= cellSizeY; 00362 continue; 00363 } 00364 00365 QList<double> scanLineList; 00366 double currentValue; 00367 GEOSGeometry* currentGeom = 0; 00368 for ( int z = 0; z < numGeoms; ++z ) 00369 { 00370 if ( numGeoms == 1 ) 00371 { 00372 currentGeom = scanLineIntersection; 00373 } 00374 else 00375 { 00376 currentGeom = GEOSGeom_clone( GEOSGetGeometryN( scanLineIntersection, z ) ); 00377 } 00378 const GEOSCoordSequence* scanLineCoordSequence = GEOSGeom_getCoordSeq( currentGeom ); 00379 if ( !scanLineCoordSequence ) 00380 { 00381 //error 00382 } 00383 unsigned int scanLineIntersectionSize; 00384 GEOSCoordSeq_getSize( scanLineCoordSequence, &scanLineIntersectionSize ); 00385 if ( !scanLineCoordSequence || scanLineIntersectionSize < 2 || ( scanLineIntersectionSize & 1 ) ) 00386 { 00387 //error 00388 } 00389 for ( unsigned int k = 0; k < scanLineIntersectionSize; ++k ) 00390 { 00391 GEOSCoordSeq_getX( scanLineCoordSequence, k, ¤tValue ); 00392 scanLineList.push_back( currentValue ); 00393 } 00394 00395 if ( numGeoms != 1 ) 00396 { 00397 GEOSGeom_destroy( currentGeom ); 00398 } 00399 } 00400 GEOSGeom_destroy( scanLineIntersection ); 00401 qSort( scanLineList ); 00402 00403 if ( scanLineList.size() < 1 ) 00404 { 00405 cellCenterY -= cellSizeY; 00406 continue; 00407 } 00408 00409 int listPlace = -1; 00410 for ( int j = 0; j < nCellsX; ++j ) 00411 { 00412 //currentCellCenter = QgsPoint( cellCenterX, cellCenterY ); 00413 00414 //instead of doing a contained test every time, find the place of scanLineList and check if even / odd 00415 if ( listPlace >= scanLineList.size() - 1 ) 00416 { 00417 break; 00418 } 00419 if ( cellCenterX >= scanLineList.at( listPlace + 1 ) ) 00420 { 00421 ++listPlace; 00422 if ( listPlace >= scanLineList.size() ) 00423 { 00424 break; 00425 } 00426 } 00427 if ( listPlace >= 0 && listPlace < ( scanLineList.size() - 1 ) && !( listPlace & 1 ) ) 00428 { 00429 if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values 00430 { 00431 sum += scanLine[j]; 00432 ++count; 00433 } 00434 } 00435 cellCenterX += cellSizeX; 00436 } 00437 cellCenterY -= cellSizeY; 00438 } 00439 CPLFree( scanLine ); 00440 } 00441 00442