Quantum GIS API Documentation
1.7.4
|
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