Quantum GIS API Documentation  1.7.4
src/analysis/vector/qgszonalstatistics.cpp
Go to the documentation of this file.
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( &currentCellCenter ) )
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, &currentValue );
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines