Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
qgsrasterdataprovider.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasterdataprovider.cpp - DataProvider Interface for raster layers
3  --------------------------------------
4  Date : Mar 11, 2005
5  Copyright : (C) 2005 by Brendan Morley
6  email : morb at ozemail dot com dot au
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgsrasterdataprovider.h"
19 #include "qgsrasterprojector.h"
20 #include "qgslogger.h"
21 
22 #include <QTime>
23 #include <QMap>
24 #include <QByteArray>
25 
26 void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, QgsCoordinateReferenceSystem theSrcCRS, QgsCoordinateReferenceSystem theDestCRS, void *data )
27 {
28  QgsDebugMsg( "Entered" );
29  QgsDebugMsg( "viewExtent = " + viewExtent.toString() );
30 
31  if ( ! theSrcCRS.isValid() || ! theDestCRS.isValid() || theSrcCRS == theDestCRS )
32  {
33  readBlock( bandNo, viewExtent, width, height, data );
34  return;
35  }
36 
37  QTime time;
38  time.start();
39 
40  double mMaxSrcXRes = 0;
41  double mMaxSrcYRes = 0;
42 
44  {
45  mMaxSrcXRes = extent().width() / xSize();
46  mMaxSrcYRes = extent().height() / ySize();
47  }
48 
49  QgsRasterProjector myProjector( theSrcCRS, theDestCRS, viewExtent, height, width, mMaxSrcXRes, mMaxSrcYRes, extent() );
50 
51  QgsDebugMsg( QString( "create projector time (ms): %1" ).arg( time.elapsed() ) );
52 
53  // TODO: init data by nulls
54 
55  // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers
56  if ( myProjector.srcRows() <= 0 || myProjector.srcCols() <= 0 )
57  return;
58 
59  // Allocate memory for not projected source data
60  int mySize = dataTypeSize( bandNo ) / 8;
61  void *mySrcData = malloc( mySize * myProjector.srcRows() * myProjector.srcCols() );
62 
63  time.restart();
64 
65  readBlock( bandNo, myProjector.srcExtent(), myProjector.srcCols(), myProjector.srcRows(), mySrcData );
66 
67  QgsDebugMsg( QString( "read not projected block time (ms): %1" ).arg( time.elapsed() ) );
68  time.restart();
69 
70  // Project data from source
71  int mySrcRow;
72  int mySrcCol;
73  int mySrcOffset;
74  int myDestOffset;
75  for ( int r = 0; r < height; r++ )
76  {
77  for ( int c = 0; c < width; c++ )
78  {
79  myProjector.srcRowCol( r, c, &mySrcRow, &mySrcCol );
80  mySrcOffset = mySize * ( mySrcRow * myProjector.srcCols() + mySrcCol );
81  myDestOffset = mySize * ( r * width + c );
82  // retype to char is just to avoid g++ warning
83  memcpy(( char* ) data + myDestOffset, ( char* )mySrcData + mySrcOffset, mySize );
84  }
85  }
86  QgsDebugMsg( QString( "reproject block time (ms): %1" ).arg( time.elapsed() ) );
87 
88  free( mySrcData );
89 }
90 
91 
93 {
94 }
95 
97  : QgsDataProvider( uri )
98  , mDpi( -1 )
99 {
100 }
101 
102 //
103 //Random Static convenience function
104 //
106 //TODO: Change these to private function or make seprate class
107 // convenience function for building metadata() HTML table cells
108 // convenience function for creating a string list from a C style string list
109 QStringList QgsRasterDataProvider::cStringList2Q_( char ** stringList )
110 {
111  QStringList strings;
112 
113  // presume null terminated string list
114  for ( size_t i = 0; stringList[i]; ++i )
115  {
116  strings.append( stringList[i] );
117  }
118 
119  return strings;
120 
121 } // cStringList2Q_
122 
123 
124 QString QgsRasterDataProvider::makeTableCell( QString const & value )
125 {
126  return "<p>\n" + value + "</p>\n";
127 } // makeTableCell_
128 
129 
130 
131 // convenience function for building metadata() HTML table cells
132 QString QgsRasterDataProvider::makeTableCells( QStringList const & values )
133 {
134  QString s( "<tr>" );
135 
136  for ( QStringList::const_iterator i = values.begin();
137  i != values.end();
138  ++i )
139  {
141  }
142 
143  s += "</tr>";
144 
145  return s;
146 } // makeTableCell_
147 
149 {
150  QString s;
151  return s;
152 }
153 
155 {
156  QStringList abilitiesList;
157 
158  int abilities = capabilities();
159 
160  if ( abilities & QgsRasterDataProvider::Identify )
161  {
162  abilitiesList += tr( "Identify" );
163  }
164 
165  if ( abilities & QgsRasterDataProvider::BuildPyramids )
166  {
167  abilitiesList += tr( "Build Pyramids" );
168  }
169 
170  QgsDebugMsg( "Capability: " + abilitiesList.join( ", " ) );
171 
172  return abilitiesList.join( ", " );
173 }
174 
175 bool QgsRasterDataProvider::identify( const QgsPoint& thePoint, QMap<QString, QString>& theResults )
176 {
177  Q_UNUSED( thePoint );
178  theResults.clear();
179  return false;
180 }
181 
183 {
184  return "text/plain";
185 }
186 
187 QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
188 {
189  int type = dataType( theBandNo );
190  int size = dataTypeSize( theBandNo ) / 8;
191  QByteArray ba;
192  ba.resize( size );
193  char * data = ba.data();
194  double noval = mNoDataValue[theBandNo-1];
195  unsigned char uc;
196  unsigned short us;
197  short s;
198  unsigned int ui;
199  int i;
200  float f;
201  double d;
202  switch ( type )
203  {
204  case Byte:
205  uc = ( unsigned char )noval;
206  memcpy( data, &uc, size );
207  break;
208  case UInt16:
209  us = ( unsigned short )noval;
210  memcpy( data, &us, size );
211  break;
212  case Int16:
213  s = ( short )noval;
214  memcpy( data, &s, size );
215  break;
216  case UInt32:
217  ui = ( unsigned int )noval;
218  memcpy( data, &ui, size );
219  break;
220  case Int32:
221  i = ( int )noval;
222  memcpy( data, &i, size );
223  break;
224  case Float32:
225  f = ( float )noval;
226  memcpy( data, &f, size );
227  break;
228  case Float64:
229  d = ( double )noval;
230  memcpy( data, &d, size );
231  break;
232  default:
233  QgsLogger::warning( "GDAL data type is not supported" );
234  }
235  return ba;
236 }
237 
238 
240 {
241  double myNoDataValue = noDataValue();
242  QgsRasterBandStats myRasterBandStats;
243  myRasterBandStats.elementCount = 0; // because we'll be counting only VALID pixels later
244  myRasterBandStats.bandName = generateBandName( theBandNo );
245  myRasterBandStats.bandNumber = theBandNo;
246 
247  int myDataType = dataType( theBandNo );
248 
249  int myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize;
250  myXBlockSize = xBlockSize();
251  myYBlockSize = yBlockSize();
252 
253  myNXBlocks = ( xSize() + myXBlockSize - 1 ) / myXBlockSize;
254  myNYBlocks = ( ySize() + myYBlockSize - 1 ) / myYBlockSize;
255 
256  void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( dataTypeSize( theBandNo ) / 8 ) );
257 
258  // unfortunately we need to make two passes through the data to calculate stddev
259  bool myFirstIterationFlag = true;
260 
261  int myBandXSize = xSize();
262  int myBandYSize = ySize();
263  for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
264  {
265  for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
266  {
267  int nXValid, nYValid;
268  readBlock( theBandNo, iXBlock, iYBlock, myData );
269 
270  // Compute the portion of the block that is valid
271  // for partial edge blocks.
272  if (( iXBlock + 1 ) * myXBlockSize > myBandXSize )
273  nXValid = myBandXSize - iXBlock * myXBlockSize;
274  else
275  nXValid = myXBlockSize;
276 
277  if (( iYBlock + 1 ) * myYBlockSize > myBandYSize )
278  nYValid = myBandYSize - iYBlock * myYBlockSize;
279  else
280  nYValid = myYBlockSize;
281 
282  // Collect the histogram counts.
283  for ( int iY = 0; iY < nYValid; iY++ )
284  {
285  for ( int iX = 0; iX < nXValid; iX++ )
286  {
287  double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
288  //QgsDebugMsg ( QString ( "%1 %2 value %3" ).arg (iX).arg(iY).arg( myValue ) );
289 
290  if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
291  {
292  continue; // NULL
293  }
294 
295  myRasterBandStats.sum += myValue;
296  ++myRasterBandStats.elementCount;
297  //only use this element if we have a non null element
298  if ( myFirstIterationFlag )
299  {
300  //this is the first iteration so initialise vars
301  myFirstIterationFlag = false;
302  myRasterBandStats.minimumValue = myValue;
303  myRasterBandStats.maximumValue = myValue;
304  } //end of true part for first iteration check
305  else
306  {
307  //this is done for all subsequent iterations
308  if ( myValue < myRasterBandStats.minimumValue )
309  {
310  myRasterBandStats.minimumValue = myValue;
311  }
312  if ( myValue > myRasterBandStats.maximumValue )
313  {
314  myRasterBandStats.maximumValue = myValue;
315  }
316  } //end of false part for first iteration check
317  }
318  }
319  } //end of column wise loop
320  } //end of row wise loop
321 
322 
323  //end of first pass through data now calculate the range
324  myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
325  //calculate the mean
326  myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount;
327 
328  //for the second pass we will get the sum of the squares / mean
329  for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
330  {
331  for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
332  {
333  int nXValid, nYValid;
334 
335  readBlock( theBandNo, iXBlock, iYBlock, myData );
336 
337  // Compute the portion of the block that is valid
338  // for partial edge blocks.
339  if (( iXBlock + 1 ) * myXBlockSize > myBandXSize )
340  nXValid = myBandXSize - iXBlock * myXBlockSize;
341  else
342  nXValid = myXBlockSize;
343 
344  if (( iYBlock + 1 ) * myYBlockSize > myBandYSize )
345  nYValid = myBandYSize - iYBlock * myYBlockSize;
346  else
347  nYValid = myYBlockSize;
348 
349  // Collect the histogram counts.
350  for ( int iY = 0; iY < nYValid; iY++ )
351  {
352  for ( int iX = 0; iX < nXValid; iX++ )
353  {
354  double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
355  //QgsDebugMsg ( "myValue = " + QString::number(myValue) );
356 
357  if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
358  {
359  continue; // NULL
360  }
361 
362  myRasterBandStats.sumOfSquares += static_cast < double >
363  ( pow( myValue - myRasterBandStats.mean, 2 ) );
364  }
365  }
366  } //end of column wise loop
367  } //end of row wise loop
368 
369  //divide result by sample size - 1 and get square root to get stdev
370  myRasterBandStats.stdDev = static_cast < double >( sqrt( myRasterBandStats.sumOfSquares /
371  ( myRasterBandStats.elementCount - 1 ) ) );
372 
373 #ifdef QGISDEBUG
374  QgsLogger::debug( "************ STATS **************", 1, __FILE__, __FUNCTION__, __LINE__ );
375  QgsLogger::debug( "VALID NODATA", mValidNoDataValue, 1, __FILE__, __FUNCTION__, __LINE__ );
376  QgsLogger::debug( "NULL", noDataValue() , 1, __FILE__, __FUNCTION__, __LINE__ );
377  QgsLogger::debug( "MIN", myRasterBandStats.minimumValue, 1, __FILE__, __FUNCTION__, __LINE__ );
378  QgsLogger::debug( "MAX", myRasterBandStats.maximumValue, 1, __FILE__, __FUNCTION__, __LINE__ );
379  QgsLogger::debug( "RANGE", myRasterBandStats.range, 1, __FILE__, __FUNCTION__, __LINE__ );
380  QgsLogger::debug( "MEAN", myRasterBandStats.mean, 1, __FILE__, __FUNCTION__, __LINE__ );
381  QgsLogger::debug( "STDDEV", myRasterBandStats.stdDev, 1, __FILE__, __FUNCTION__, __LINE__ );
382 #endif
383 
384  CPLFree( myData );
385  myRasterBandStats.statsGathered = true;
386  return myRasterBandStats;
387 }
388 
389 double QgsRasterDataProvider::readValue( void *data, int type, int index )
390 {
391  if ( !data )
392  return mValidNoDataValue ? noDataValue() : 0.0;
393 
394  switch ( type )
395  {
396  case Byte:
397  return ( double )(( GByte * )data )[index];
398  break;
399  case UInt16:
400  return ( double )(( GUInt16 * )data )[index];
401  break;
402  case Int16:
403  return ( double )(( GInt16 * )data )[index];
404  break;
405  case UInt32:
406  return ( double )(( GUInt32 * )data )[index];
407  break;
408  case Int32:
409  return ( double )(( GInt32 * )data )[index];
410  break;
411  case Float32:
412  return ( double )(( float * )data )[index];
413  break;
414  case Float64:
415  return ( double )(( double * )data )[index];
416  break;
417  default:
418  QgsLogger::warning( "GDAL data type is not supported" );
419  }
420 
421  return mValidNoDataValue ? noDataValue() : 0.0;
422 }
423 
424 // ENDS