Quantum GIS API Documentation  1.7.4
src/core/qgsrasterdataprovider.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002     qgsrasterdataprovider.cpp - DataProvider Interface for raster layers
00003      --------------------------------------
00004     Date                 : Mar 11, 2005
00005     Copyright            : (C) 2005 by Brendan Morley
00006     email                : morb at ozemail dot com dot au
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 "qgsrasterdataprovider.h"
00019 #include "qgsrasterprojector.h"
00020 #include "qgslogger.h"
00021 
00022 #include <QTime>
00023 #include <QMap>
00024 #include <QByteArray>
00025 
00026 void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle  const & viewExtent, int width, int height, QgsCoordinateReferenceSystem theSrcCRS, QgsCoordinateReferenceSystem theDestCRS, void *data )
00027 {
00028   QgsDebugMsg( "Entered" );
00029   QgsDebugMsg( "viewExtent = " + viewExtent.toString() );
00030 
00031   if ( ! theSrcCRS.isValid() || ! theDestCRS.isValid() || theSrcCRS == theDestCRS )
00032   {
00033     readBlock( bandNo, viewExtent, width, height, data );
00034     return;
00035   }
00036 
00037   QTime time;
00038   time.start();
00039 
00040   double mMaxSrcXRes = 0;
00041   double mMaxSrcYRes = 0;
00042 
00043   if ( capabilities() & QgsRasterDataProvider::ExactResolution )
00044   {
00045     mMaxSrcXRes = extent().width() / xSize();
00046     mMaxSrcYRes = extent().height() / ySize();
00047   }
00048 
00049   QgsRasterProjector myProjector( theSrcCRS, theDestCRS, viewExtent, height, width, mMaxSrcXRes, mMaxSrcYRes, extent() );
00050 
00051   QgsDebugMsg( QString( "create projector time  (ms): %1" ).arg( time.elapsed() ) );
00052 
00053   // TODO: init data by nulls
00054 
00055   // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers
00056   if ( myProjector.srcRows() <= 0 || myProjector.srcCols() <= 0 )
00057     return;
00058 
00059   // Allocate memory for not projected source data
00060   int mySize = dataTypeSize( bandNo ) / 8;
00061   void *mySrcData = malloc( mySize * myProjector.srcRows() * myProjector.srcCols() );
00062 
00063   time.restart();
00064 
00065   readBlock( bandNo, myProjector.srcExtent(), myProjector.srcCols(), myProjector.srcRows(), mySrcData );
00066 
00067   QgsDebugMsg( QString( "read not projected block time  (ms): %1" ).arg( time.elapsed() ) );
00068   time.restart();
00069 
00070   // Project data from source
00071   int mySrcRow;
00072   int mySrcCol;
00073   int mySrcOffset;
00074   int myDestOffset;
00075   for ( int r = 0; r < height; r++ )
00076   {
00077     for ( int c = 0; c < width; c++ )
00078     {
00079       myProjector.srcRowCol( r, c, &mySrcRow, &mySrcCol );
00080       mySrcOffset = mySize * ( mySrcRow * myProjector.srcCols() + mySrcCol );
00081       myDestOffset = mySize * ( r * width + c );
00082       // retype to char is just to avoid g++ warning
00083       memcpy(( char* ) data + myDestOffset, ( char* )mySrcData + mySrcOffset, mySize );
00084     }
00085   }
00086   QgsDebugMsg( QString( "reproject block time  (ms): %1" ).arg( time.elapsed() ) );
00087 
00088   free( mySrcData );
00089 }
00090 
00091 
00092 QgsRasterDataProvider::QgsRasterDataProvider(): mDpi( -1 )
00093 {
00094 }
00095 
00096 QgsRasterDataProvider::QgsRasterDataProvider( QString const & uri )
00097     : QgsDataProvider( uri )
00098     , mDpi( -1 )
00099 {
00100 }
00101 
00102 //
00103 //Random Static convenience function
00104 //
00106 //TODO: Change these to private function or make seprate class
00107 // convenience function for building metadata() HTML table cells
00108 // convenience function for creating a string list from a C style string list
00109 QStringList QgsRasterDataProvider::cStringList2Q_( char ** stringList )
00110 {
00111   QStringList strings;
00112 
00113   // presume null terminated string list
00114   for ( size_t i = 0; stringList[i]; ++i )
00115   {
00116     strings.append( stringList[i] );
00117   }
00118 
00119   return strings;
00120 
00121 } // cStringList2Q_
00122 
00123 
00124 QString QgsRasterDataProvider::makeTableCell( QString const & value )
00125 {
00126   return "<p>\n" + value + "</p>\n";
00127 } // makeTableCell_
00128 
00129 
00130 
00131 // convenience function for building metadata() HTML table cells
00132 QString QgsRasterDataProvider::makeTableCells( QStringList const & values )
00133 {
00134   QString s( "<tr>" );
00135 
00136   for ( QStringList::const_iterator i = values.begin();
00137         i != values.end();
00138         ++i )
00139   {
00140     s += QgsRasterDataProvider::makeTableCell( *i );
00141   }
00142 
00143   s += "</tr>";
00144 
00145   return s;
00146 } // makeTableCell_
00147 
00148 QString QgsRasterDataProvider::metadata()
00149 {
00150   QString s;
00151   return s;
00152 }
00153 
00154 QString QgsRasterDataProvider::capabilitiesString() const
00155 {
00156   QStringList abilitiesList;
00157 
00158   int abilities = capabilities();
00159 
00160   if ( abilities & QgsRasterDataProvider::Identify )
00161   {
00162     abilitiesList += tr( "Identify" );
00163   }
00164 
00165   if ( abilities & QgsRasterDataProvider::BuildPyramids )
00166   {
00167     abilitiesList += tr( "Build Pyramids" );
00168   }
00169 
00170   QgsDebugMsg( "Capability: " + abilitiesList.join( ", " ) );
00171 
00172   return abilitiesList.join( ", " );
00173 }
00174 
00175 bool QgsRasterDataProvider::identify( const QgsPoint& thePoint, QMap<QString, QString>& theResults )
00176 {
00177   Q_UNUSED( thePoint );
00178   theResults.clear();
00179   return false;
00180 }
00181 
00182 QString QgsRasterDataProvider::lastErrorFormat()
00183 {
00184   return "text/plain";
00185 }
00186 
00187 QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
00188 {
00189   int type = dataType( theBandNo );
00190   int size = dataTypeSize( theBandNo ) / 8;
00191   QByteArray ba;
00192   ba.resize( size );
00193   char * data = ba.data();
00194   double noval = mNoDataValue[theBandNo-1];
00195   unsigned char uc;
00196   unsigned short us;
00197   short s;
00198   unsigned int ui;
00199   int i;
00200   float f;
00201   double d;
00202   switch ( type )
00203   {
00204     case Byte:
00205       uc = ( unsigned char )noval;
00206       memcpy( data, &uc, size );
00207       break;
00208     case UInt16:
00209       us = ( unsigned short )noval;
00210       memcpy( data, &us, size );
00211       break;
00212     case Int16:
00213       s = ( short )noval;
00214       memcpy( data, &s, size );
00215       break;
00216     case UInt32:
00217       ui = ( unsigned int )noval;
00218       memcpy( data, &ui, size );
00219       break;
00220     case Int32:
00221       i = ( int )noval;
00222       memcpy( data, &i, size );
00223       break;
00224     case Float32:
00225       f = ( float )noval;
00226       memcpy( data, &f, size );
00227       break;
00228     case Float64:
00229       d = ( double )noval;
00230       memcpy( data, &d, size );
00231       break;
00232     default:
00233       QgsLogger::warning( "GDAL data type is not supported" );
00234   }
00235   return ba;
00236 }
00237 
00238 
00239 QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo )
00240 {
00241   double myNoDataValue = noDataValue();
00242   QgsRasterBandStats myRasterBandStats;
00243   myRasterBandStats.elementCount = 0; // because we'll be counting only VALID pixels later
00244   myRasterBandStats.bandName = generateBandName( theBandNo );
00245   myRasterBandStats.bandNumber = theBandNo;
00246 
00247   int myDataType = dataType( theBandNo );
00248 
00249   int  myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize;
00250   myXBlockSize = xBlockSize();
00251   myYBlockSize = yBlockSize();
00252 
00253   myNXBlocks = ( xSize() + myXBlockSize - 1 ) / myXBlockSize;
00254   myNYBlocks = ( ySize() + myYBlockSize - 1 ) / myYBlockSize;
00255 
00256   void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( dataTypeSize( theBandNo ) / 8 ) );
00257 
00258   // unfortunately we need to make two passes through the data to calculate stddev
00259   bool myFirstIterationFlag = true;
00260 
00261   int myBandXSize = xSize();
00262   int myBandYSize = ySize();
00263   for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
00264   {
00265     for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
00266     {
00267       int  nXValid, nYValid;
00268       readBlock( theBandNo, iXBlock, iYBlock, myData );
00269 
00270       // Compute the portion of the block that is valid
00271       // for partial edge blocks.
00272       if (( iXBlock + 1 ) * myXBlockSize > myBandXSize )
00273         nXValid = myBandXSize - iXBlock * myXBlockSize;
00274       else
00275         nXValid = myXBlockSize;
00276 
00277       if (( iYBlock + 1 ) * myYBlockSize > myBandYSize )
00278         nYValid = myBandYSize - iYBlock * myYBlockSize;
00279       else
00280         nYValid = myYBlockSize;
00281 
00282       // Collect the histogram counts.
00283       for ( int iY = 0; iY < nYValid; iY++ )
00284       {
00285         for ( int iX = 0; iX < nXValid; iX++ )
00286         {
00287           double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
00288           //QgsDebugMsg ( QString ( "%1 %2 value %3" ).arg (iX).arg(iY).arg( myValue ) );
00289 
00290           if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
00291           {
00292             continue; // NULL
00293           }
00294 
00295           myRasterBandStats.sum += myValue;
00296           ++myRasterBandStats.elementCount;
00297           //only use this element if we have a non null element
00298           if ( myFirstIterationFlag )
00299           {
00300             //this is the first iteration so initialise vars
00301             myFirstIterationFlag = false;
00302             myRasterBandStats.minimumValue = myValue;
00303             myRasterBandStats.maximumValue = myValue;
00304           }               //end of true part for first iteration check
00305           else
00306           {
00307             //this is done for all subsequent iterations
00308             if ( myValue < myRasterBandStats.minimumValue )
00309             {
00310               myRasterBandStats.minimumValue = myValue;
00311             }
00312             if ( myValue > myRasterBandStats.maximumValue )
00313             {
00314               myRasterBandStats.maximumValue = myValue;
00315             }
00316           } //end of false part for first iteration check
00317         }
00318       }
00319     } //end of column wise loop
00320   } //end of row wise loop
00321 
00322 
00323   //end of first pass through data now calculate the range
00324   myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
00325   //calculate the mean
00326   myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount;
00327 
00328   //for the second pass we will get the sum of the squares / mean
00329   for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
00330   {
00331     for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
00332     {
00333       int  nXValid, nYValid;
00334 
00335       readBlock( theBandNo, iXBlock, iYBlock, myData );
00336 
00337       // Compute the portion of the block that is valid
00338       // for partial edge blocks.
00339       if (( iXBlock + 1 ) * myXBlockSize > myBandXSize )
00340         nXValid = myBandXSize - iXBlock * myXBlockSize;
00341       else
00342         nXValid = myXBlockSize;
00343 
00344       if (( iYBlock + 1 ) * myYBlockSize > myBandYSize )
00345         nYValid = myBandYSize - iYBlock * myYBlockSize;
00346       else
00347         nYValid = myYBlockSize;
00348 
00349       // Collect the histogram counts.
00350       for ( int iY = 0; iY < nYValid; iY++ )
00351       {
00352         for ( int iX = 0; iX < nXValid; iX++ )
00353         {
00354           double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
00355           //QgsDebugMsg ( "myValue = " + QString::number(myValue) );
00356 
00357           if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
00358           {
00359             continue; // NULL
00360           }
00361 
00362           myRasterBandStats.sumOfSquares += static_cast < double >
00363                                             ( pow( myValue - myRasterBandStats.mean, 2 ) );
00364         }
00365       }
00366     } //end of column wise loop
00367   } //end of row wise loop
00368 
00369   //divide result by sample size - 1 and get square root to get stdev
00370   myRasterBandStats.stdDev = static_cast < double >( sqrt( myRasterBandStats.sumOfSquares /
00371                              ( myRasterBandStats.elementCount - 1 ) ) );
00372 
00373 #ifdef QGISDEBUG
00374   QgsLogger::debug( "************ STATS **************", 1, __FILE__, __FUNCTION__, __LINE__ );
00375   QgsLogger::debug( "VALID NODATA", mValidNoDataValue, 1, __FILE__, __FUNCTION__, __LINE__ );
00376   QgsLogger::debug( "NULL", noDataValue() , 1, __FILE__, __FUNCTION__, __LINE__ );
00377   QgsLogger::debug( "MIN", myRasterBandStats.minimumValue, 1, __FILE__, __FUNCTION__, __LINE__ );
00378   QgsLogger::debug( "MAX", myRasterBandStats.maximumValue, 1, __FILE__, __FUNCTION__, __LINE__ );
00379   QgsLogger::debug( "RANGE", myRasterBandStats.range, 1, __FILE__, __FUNCTION__, __LINE__ );
00380   QgsLogger::debug( "MEAN", myRasterBandStats.mean, 1, __FILE__, __FUNCTION__, __LINE__ );
00381   QgsLogger::debug( "STDDEV", myRasterBandStats.stdDev, 1, __FILE__, __FUNCTION__, __LINE__ );
00382 #endif
00383 
00384   CPLFree( myData );
00385   myRasterBandStats.statsGathered = true;
00386   return myRasterBandStats;
00387 }
00388 
00389 double QgsRasterDataProvider::readValue( void *data, int type, int index )
00390 {
00391   if ( !data )
00392     return mValidNoDataValue ? noDataValue() : 0.0;
00393 
00394   switch ( type )
00395   {
00396     case Byte:
00397       return ( double )(( GByte * )data )[index];
00398       break;
00399     case UInt16:
00400       return ( double )(( GUInt16 * )data )[index];
00401       break;
00402     case Int16:
00403       return ( double )(( GInt16 * )data )[index];
00404       break;
00405     case UInt32:
00406       return ( double )(( GUInt32 * )data )[index];
00407       break;
00408     case Int32:
00409       return ( double )(( GInt32 * )data )[index];
00410       break;
00411     case Float32:
00412       return ( double )(( float * )data )[index];
00413       break;
00414     case Float64:
00415       return ( double )(( double * )data )[index];
00416       break;
00417     default:
00418       QgsLogger::warning( "GDAL data type is not supported" );
00419   }
00420 
00421   return mValidNoDataValue ? noDataValue() : 0.0;
00422 }
00423 
00424 // ENDS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines