QGIS API Documentation  2.0.1-Dufour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgsrasterinterface.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasterface.cpp - Internal raster processing modules interface
3  --------------------------------------
4  Date : Jun 21, 2012
5  Copyright : (C) 2012 by Radim Blazek
6  email : radim dot blazek at gmail dot com
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 <limits>
19 #include <typeinfo>
20 
21 #include <QByteArray>
22 #include <QTime>
23 
24 #include <qmath.h>
25 
26 #include "qgslogger.h"
27 #include "qgsrasterbandstats.h"
28 #include "qgsrasterhistogram.h"
29 #include "qgsrasterinterface.h"
30 #include "qgsrectangle.h"
31 
33  : mInput( input )
34  , mOn( true )
35 {
36 }
37 
39 {
40 }
41 
43  int theBandNo,
44  int theStats,
45  const QgsRectangle & theExtent,
46  int theSampleSize )
47 {
48  QgsDebugMsg( QString( "theBandNo = %1 theSampleSize = %2" ).arg( theBandNo ).arg( theSampleSize ) );
49 
50  theStatistics.bandNumber = theBandNo;
51  theStatistics.statsGathered = theStats;
52 
53  QgsRectangle myExtent;
54  if ( theExtent.isEmpty() )
55  {
56  myExtent = extent();
57  }
58  else
59  {
60  myExtent = extent().intersect( &theExtent );
61  }
62  theStatistics.extent = myExtent;
63 
64  if ( theSampleSize > 0 )
65  {
66  // Calc resolution from theSampleSize
67  double xRes, yRes;
68  xRes = yRes = sqrt(( myExtent.width() * myExtent.height() ) / theSampleSize );
69 
70  // But limit by physical resolution
71  if ( capabilities() & Size )
72  {
73  double srcXRes = extent().width() / xSize();
74  double srcYRes = extent().height() / ySize();
75  if ( xRes < srcXRes ) xRes = srcXRes;
76  if ( yRes < srcYRes ) yRes = srcYRes;
77  }
78  QgsDebugMsg( QString( "xRes = %1 yRes = %2" ).arg( xRes ).arg( yRes ) );
79 
80  theStatistics.width = static_cast <int>( myExtent.width() / xRes );
81  theStatistics.height = static_cast <int>( myExtent.height() / yRes );
82  }
83  else
84  {
85  if ( capabilities() & Size )
86  {
87  theStatistics.width = xSize();
88  theStatistics.height = ySize();
89  }
90  else
91  {
92  theStatistics.width = 1000;
93  theStatistics.height = 1000;
94  }
95  }
96  QgsDebugMsg( QString( "theStatistics.width = %1 theStatistics.height = %2" ).arg( theStatistics.width ).arg( theStatistics.height ) );
97 }
98 
100  int theStats,
101  const QgsRectangle & theExtent,
102  int theSampleSize )
103 {
104  QgsDebugMsg( QString( "theBandNo = %1 theStats = %2 theSampleSize = %3" ).arg( theBandNo ).arg( theStats ).arg( theSampleSize ) );
105  if ( mStatistics.size() == 0 ) return false;
106 
107  QgsRasterBandStats myRasterBandStats;
108  initStatistics( myRasterBandStats, theBandNo, theStats, theExtent, theSampleSize );
109 
110  foreach ( QgsRasterBandStats stats, mStatistics )
111  {
112  if ( stats.contains( myRasterBandStats ) )
113  {
114  QgsDebugMsg( "Has cached statistics." );
115  return true;
116  }
117  }
118  return false;
119 }
120 
122  int theStats,
123  const QgsRectangle & theExtent,
124  int theSampleSize )
125 {
126  QgsDebugMsg( QString( "theBandNo = %1 theStats = %2 theSampleSize = %3" ).arg( theBandNo ).arg( theStats ).arg( theSampleSize ) );
127 
128  // TODO: null values set on raster layer!!!
129 
130  QgsRasterBandStats myRasterBandStats;
131  initStatistics( myRasterBandStats, theBandNo, theStats, theExtent, theSampleSize );
132 
133  foreach ( QgsRasterBandStats stats, mStatistics )
134  {
135  if ( stats.contains( myRasterBandStats ) )
136  {
137  QgsDebugMsg( "Using cached statistics." );
138  return stats;
139  }
140  }
141 
142  QgsRectangle myExtent = myRasterBandStats.extent;
143  int myWidth = myRasterBandStats.width;
144  int myHeight = myRasterBandStats.height;
145 
146  //int myDataType = dataType( theBandNo );
147 
148  int myXBlockSize = xBlockSize();
149  int myYBlockSize = yBlockSize();
150  if ( myXBlockSize == 0 ) // should not happen, but happens
151  {
152  myXBlockSize = 500;
153  }
154  if ( myYBlockSize == 0 ) // should not happen, but happens
155  {
156  myYBlockSize = 500;
157  }
158 
159  int myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
160  int myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
161 
162  double myXRes = myExtent.width() / myWidth;
163  double myYRes = myExtent.height() / myHeight;
164  // TODO: progress signals
165 
166  // used by single pass stdev
167  double myMean = 0;
168  double mySumOfSquares = 0;
169 
170  bool myFirstIterationFlag = true;
171  for ( int myYBlock = 0; myYBlock < myNYBlocks; myYBlock++ )
172  {
173  for ( int myXBlock = 0; myXBlock < myNXBlocks; myXBlock++ )
174  {
175  QgsDebugMsg( QString( "myYBlock = %1 myXBlock = %2" ).arg( myYBlock ).arg( myXBlock ) );
176  int myBlockWidth = qMin( myXBlockSize, myWidth - myXBlock * myXBlockSize );
177  int myBlockHeight = qMin( myYBlockSize, myHeight - myYBlock * myYBlockSize );
178 
179  double xmin = myExtent.xMinimum() + myXBlock * myXBlockSize * myXRes;
180  double xmax = xmin + myBlockWidth * myXRes;
181  double ymin = myExtent.yMaximum() - myYBlock * myYBlockSize * myYRes;
182  double ymax = ymin - myBlockHeight * myYRes;
183 
184  QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
185 
186  QgsRasterBlock* blk = block( theBandNo, myPartExtent, myBlockWidth, myBlockHeight );
187 
188  // Collect the histogram counts.
189  for ( size_t i = 0; i < (( size_t ) myBlockHeight ) * myBlockWidth; i++ )
190  {
191  if ( blk->isNoData( i ) ) continue; // NULL
192 
193  double myValue = blk->value( i );
194 
195  myRasterBandStats.sum += myValue;
196  myRasterBandStats.elementCount++;
197 
198  if ( myFirstIterationFlag )
199  {
200  myFirstIterationFlag = false;
201  myRasterBandStats.minimumValue = myValue;
202  myRasterBandStats.maximumValue = myValue;
203  }
204  else
205  {
206  if ( myValue < myRasterBandStats.minimumValue )
207  {
208  myRasterBandStats.minimumValue = myValue;
209  }
210  if ( myValue > myRasterBandStats.maximumValue )
211  {
212  myRasterBandStats.maximumValue = myValue;
213  }
214  }
215 
216  // Single pass stdev
217  double myDelta = myValue - myMean;
218  myMean += myDelta / myRasterBandStats.elementCount;
219  mySumOfSquares += myDelta * ( myValue - myMean );
220  }
221  delete blk;
222  }
223  }
224 
225  myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
226  myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount;
227 
228  myRasterBandStats.sumOfSquares = mySumOfSquares; // OK with single pass?
229 
230  // stdDev may differ from GDAL stats, because GDAL is using naive single pass
231  // algorithm which is more error prone (because of rounding errors)
232  // Divide result by sample size - 1 and get square root to get stdev
233  myRasterBandStats.stdDev = sqrt( mySumOfSquares / ( myRasterBandStats.elementCount - 1 ) );
234 
235  QgsDebugMsg( "************ STATS **************" );
236  QgsDebugMsg( QString( "MIN %1" ).arg( myRasterBandStats.minimumValue ) );
237  QgsDebugMsg( QString( "MAX %1" ).arg( myRasterBandStats.maximumValue ) );
238  QgsDebugMsg( QString( "RANGE %1" ).arg( myRasterBandStats.range ) );
239  QgsDebugMsg( QString( "MEAN %1" ).arg( myRasterBandStats.mean ) );
240  QgsDebugMsg( QString( "STDDEV %1" ).arg( myRasterBandStats.stdDev ) );
241 
242  myRasterBandStats.statsGathered = QgsRasterBandStats::All;
243  mStatistics.append( myRasterBandStats );
244 
245  return myRasterBandStats;
246 }
247 
249  int theBandNo,
250  int theBinCount,
251  double theMinimum, double theMaximum,
252  const QgsRectangle & theExtent,
253  int theSampleSize,
254  bool theIncludeOutOfRange )
255 {
256  theHistogram.bandNumber = theBandNo;
257  theHistogram.minimum = theMinimum;
258  theHistogram.maximum = theMaximum;
259  theHistogram.includeOutOfRange = theIncludeOutOfRange;
260 
261  int mySrcDataType = srcDataType( theBandNo );
262 
263  if ( qIsNaN( theHistogram.minimum ) )
264  {
265  // TODO: this was OK when stats/histogram were calced in provider,
266  // but what TODO in other interfaces? Check for mInput for now.
267  if ( !mInput && mySrcDataType == QGis::Byte )
268  {
269  theHistogram.minimum = 0; // see histogram() for shift for rounding
270  }
271  else
272  {
273  // We need statistics -> avoid histogramDefaults in hasHistogram if possible
274  // TODO: use approximated statistics if aproximated histogram is requested
275  // (theSampleSize > 0)
276  QgsRasterBandStats stats = bandStatistics( theBandNo, QgsRasterBandStats::Min, theExtent, theSampleSize );
277  theHistogram.minimum = stats.minimumValue;
278  }
279  }
280  if ( qIsNaN( theHistogram.maximum ) )
281  {
282  if ( !mInput && mySrcDataType == QGis::Byte )
283  {
284  theHistogram.maximum = 255;
285  }
286  else
287  {
288  QgsRasterBandStats stats = bandStatistics( theBandNo, QgsRasterBandStats::Max, theExtent, theSampleSize );
289  theHistogram.maximum = stats.maximumValue;
290  }
291  }
292 
293  QgsRectangle myExtent;
294  if ( theExtent.isEmpty() )
295  {
296  myExtent = extent();
297  }
298  else
299  {
300  myExtent = extent().intersect( &theExtent );
301  }
302  theHistogram.extent = myExtent;
303 
304  if ( theSampleSize > 0 )
305  {
306  // Calc resolution from theSampleSize
307  double xRes, yRes;
308  xRes = yRes = sqrt(( myExtent.width() * myExtent.height() ) / theSampleSize );
309 
310  // But limit by physical resolution
311  if ( capabilities() & Size )
312  {
313  double srcXRes = extent().width() / xSize();
314  double srcYRes = extent().height() / ySize();
315  if ( xRes < srcXRes ) xRes = srcXRes;
316  if ( yRes < srcYRes ) yRes = srcYRes;
317  }
318  QgsDebugMsg( QString( "xRes = %1 yRes = %2" ).arg( xRes ).arg( yRes ) );
319 
320  theHistogram.width = static_cast <int>( myExtent.width() / xRes );
321  theHistogram.height = static_cast <int>( myExtent.height() / yRes );
322  }
323  else
324  {
325  if ( capabilities() & Size )
326  {
327  theHistogram.width = xSize();
328  theHistogram.height = ySize();
329  }
330  else
331  {
332  theHistogram.width = 1000;
333  theHistogram.height = 1000;
334  }
335  }
336  QgsDebugMsg( QString( "theHistogram.width = %1 theHistogram.height = %2" ).arg( theHistogram.width ).arg( theHistogram.height ) );
337 
338  int myBinCount = theBinCount;
339  if ( myBinCount == 0 )
340  {
341  // TODO: this was OK when stats/histogram were calced in provider,
342  // but what TODO in other interfaces? Check for mInput for now.
343  if ( !mInput && mySrcDataType == QGis::Byte )
344  {
345  myBinCount = 256; // Cannot store more values in byte
346  }
347  else
348  {
349  // There is no best default value, to display something reasonable in histogram chart, binCount should be small, OTOH, to get precise data for cumulative cut, the number should be big. Because it is easier to define fixed lower value for the chart, we calc optimum binCount for higher resolution (to avoid calculating that where histogram() is used. In any any case, it does not make sense to use more than width*height;
350  myBinCount = theHistogram.width * theHistogram.height;
351  if ( myBinCount > 1000 ) myBinCount = 1000;
352  }
353  }
354  theHistogram.binCount = myBinCount;
355  QgsDebugMsg( QString( "theHistogram.binCount = %1" ).arg( theHistogram.binCount ) );
356 }
357 
358 
360  int theBinCount,
361  double theMinimum, double theMaximum,
362  const QgsRectangle & theExtent,
363  int theSampleSize,
364  bool theIncludeOutOfRange )
365 {
366  QgsDebugMsg( QString( "theBandNo = %1 theBinCount = %2 theMinimum = %3 theMaximum = %4 theSampleSize = %5" ).arg( theBandNo ).arg( theBinCount ).arg( theMinimum ).arg( theMaximum ).arg( theSampleSize ) );
367  // histogramDefaults() needs statistics if theMinimum or theMaximum is NaN ->
368  // do other checks which don't need statistics before histogramDefaults()
369  if ( mHistograms.size() == 0 ) return false;
370 
371  QgsRasterHistogram myHistogram;
372  initHistogram( myHistogram, theBandNo, theBinCount, theMinimum, theMaximum, theExtent, theSampleSize, theIncludeOutOfRange );
373 
375  {
376  if ( histogram == myHistogram )
377  {
378  QgsDebugMsg( "Has cached histogram." );
379  return true;
380  }
381  }
382  return false;
383 }
384 
386  int theBinCount,
387  double theMinimum, double theMaximum,
388  const QgsRectangle & theExtent,
389  int theSampleSize,
390  bool theIncludeOutOfRange )
391 {
392  QgsDebugMsg( QString( "theBandNo = %1 theBinCount = %2 theMinimum = %3 theMaximum = %4 theSampleSize = %5" ).arg( theBandNo ).arg( theBinCount ).arg( theMinimum ).arg( theMaximum ).arg( theSampleSize ) );
393 
394  QgsRasterHistogram myHistogram;
395  initHistogram( myHistogram, theBandNo, theBinCount, theMinimum, theMaximum, theExtent, theSampleSize, theIncludeOutOfRange );
396 
397  // Find cached
399  {
400  if ( histogram == myHistogram )
401  {
402  QgsDebugMsg( "Using cached histogram." );
403  return histogram;
404  }
405  }
406 
407  int myBinCount = myHistogram.binCount;
408  int myWidth = myHistogram.width;
409  int myHeight = myHistogram.height;
410  QgsRectangle myExtent = myHistogram.extent;
411  myHistogram.histogramVector.resize( myBinCount );
412 
413  int myXBlockSize = xBlockSize();
414  int myYBlockSize = yBlockSize();
415  if ( myXBlockSize == 0 ) // should not happen, but happens
416  {
417  myXBlockSize = 500;
418  }
419  if ( myYBlockSize == 0 ) // should not happen, but happens
420  {
421  myYBlockSize = 500;
422  }
423 
424  int myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
425  int myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
426 
427  double myXRes = myExtent.width() / myWidth;
428  double myYRes = myExtent.height() / myHeight;
429 
430  double myMinimum = myHistogram.minimum;
431  double myMaximum = myHistogram.maximum;
432 
433  // To avoid rounding errors
434  // TODO: check this
435  double myerval = ( myMaximum - myMinimum ) / myHistogram.binCount;
436  myMinimum -= 0.1 * myerval;
437  myMaximum += 0.1 * myerval;
438 
439  QgsDebugMsg( QString( "binCount = %1 myMinimum = %2 myMaximum = %3" ).arg( myHistogram.binCount ).arg( myMinimum ).arg( myMaximum ) );
440 
441  double myBinSize = ( myMaximum - myMinimum ) / myBinCount;
442 
443  // TODO: progress signals
444  for ( int myYBlock = 0; myYBlock < myNYBlocks; myYBlock++ )
445  {
446  for ( int myXBlock = 0; myXBlock < myNXBlocks; myXBlock++ )
447  {
448  int myBlockWidth = qMin( myXBlockSize, myWidth - myXBlock * myXBlockSize );
449  int myBlockHeight = qMin( myYBlockSize, myHeight - myYBlock * myYBlockSize );
450 
451  double xmin = myExtent.xMinimum() + myXBlock * myXBlockSize * myXRes;
452  double xmax = xmin + myBlockWidth * myXRes;
453  double ymin = myExtent.yMaximum() - myYBlock * myYBlockSize * myYRes;
454  double ymax = ymin - myBlockHeight * myYRes;
455 
456  QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
457 
458  QgsRasterBlock* blk = block( theBandNo, myPartExtent, myBlockWidth, myBlockHeight );
459 
460  // Collect the histogram counts.
461  for ( size_t i = 0; i < (( size_t ) myBlockHeight ) * myBlockWidth; i++ )
462  {
463  if ( blk->isNoData( i ) )
464  {
465  continue; // NULL
466  }
467  double myValue = blk->value( i );
468 
469  int myBinIndex = static_cast <int>( qFloor(( myValue - myMinimum ) / myBinSize ) ) ;
470 
471  if (( myBinIndex < 0 || myBinIndex > ( myBinCount - 1 ) ) && !theIncludeOutOfRange )
472  {
473  continue;
474  }
475  if ( myBinIndex < 0 ) myBinIndex = 0;
476  if ( myBinIndex > ( myBinCount - 1 ) ) myBinIndex = myBinCount - 1;
477 
478  myHistogram.histogramVector[myBinIndex] += 1;
479  myHistogram.nonNullCount++;
480  }
481  delete blk;
482  }
483  }
484 
485  myHistogram.valid = true;
486  mHistograms.append( myHistogram );
487 
488 #ifdef QGISDEBUG
489  QString hist;
490  for ( int i = 0; i < qMin( myHistogram.histogramVector.size(), 500 ); i++ )
491  {
492  hist += QString::number( myHistogram.histogramVector.value( i ) ) + " ";
493  }
494  QgsDebugMsg( "Histogram (max first 500 bins): " + hist );
495 #endif
496 
497  return myHistogram;
498 }
499 
501  double theLowerCount, double theUpperCount,
502  double &theLowerValue, double &theUpperValue,
503  const QgsRectangle & theExtent,
504  int theSampleSize )
505 {
506  QgsDebugMsg( QString( "theBandNo = %1 theLowerCount = %2 theUpperCount = %3 theSampleSize = %4" ).arg( theBandNo ).arg( theLowerCount ).arg( theUpperCount ).arg( theSampleSize ) );
507 
508  QgsRasterHistogram myHistogram = histogram( theBandNo, 0, std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), theExtent, theSampleSize );
509 
510  // Init to NaN is better than histogram min/max to catch errors
511  theLowerValue = std::numeric_limits<double>::quiet_NaN();
512  theUpperValue = std::numeric_limits<double>::quiet_NaN();
513 
514  double myBinXStep = ( myHistogram.maximum - myHistogram.minimum ) / myHistogram.binCount;
515  int myCount = 0;
516  int myMinCount = ( int ) qRound( theLowerCount * myHistogram.nonNullCount );
517  int myMaxCount = ( int ) qRound( theUpperCount * myHistogram.nonNullCount );
518  bool myLowerFound = false;
519  QgsDebugMsg( QString( "binCount = %1 minimum = %2 maximum = %3 myBinXStep = %4" ).arg( myHistogram.binCount ).arg( myHistogram.minimum ).arg( myHistogram.maximum ).arg( myBinXStep ) );
520  QgsDebugMsg( QString( "myMinCount = %1 myMaxCount = %2" ).arg( myMinCount ).arg( myMaxCount ) );
521 
522  for ( int myBin = 0; myBin < myHistogram.histogramVector.size(); myBin++ )
523  {
524  int myBinValue = myHistogram.histogramVector.value( myBin );
525  myCount += myBinValue;
526  if ( !myLowerFound && myCount > myMinCount )
527  {
528  theLowerValue = myHistogram.minimum + myBin * myBinXStep;
529  myLowerFound = true;
530  }
531  if ( myCount >= myMaxCount )
532  {
533  theUpperValue = myHistogram.minimum + myBin * myBinXStep;
534  break;
535  }
536  }
537 }
538 
540 {
541  QStringList abilitiesList;
542 
543  int abilities = capabilities();
544 
545  // Not all all capabilities are here (Size, IdentifyValue, IdentifyText,
546  // IdentifyHtml, IdentifyFeature) because those are quite technical and probably
547  // would be confusing for users
548 
549  if ( abilities & QgsRasterInterface::Identify )
550  {
551  abilitiesList += tr( "Identify" );
552  }
553 
554  if ( abilities & QgsRasterInterface::Create )
555  {
556  abilitiesList += tr( "Create Datasources" );
557  }
558 
559  if ( abilities & QgsRasterInterface::Remove )
560  {
561  abilitiesList += tr( "Remove Datasources" );
562  }
563 
564  if ( abilities & QgsRasterInterface::BuildPyramids )
565  {
566  abilitiesList += tr( "Build Pyramids" );
567  }
568 
569  QgsDebugMsg( "Capability: " + abilitiesList.join( ", " ) );
570 
571  return abilitiesList.join( ", " );
572 }