Logo Search packages:      
Sourcecode: qgis version File versions

const RasterBandStats QgsRasterLayer::getRasterBandStats ( int  theBandNoInt  ) 

Get RasterBandStats for a band given its number (read only).

Private method to calculate statistics for a band. Populates rasterStatsMemArray.

Calculates:

  • myRasterBandStats.elementCountInt
  • myRasterBandStats.minValDouble
  • myRasterBandStats.maxValDouble
  • myRasterBandStats.sumDouble
  • myRasterBandStats.rangeDouble
  • myRasterBandStats.meanDouble
  • myRasterBandStats.sumSqrDevDouble
  • myRasterBandStats.stdDevDouble

RasterBandStats

Note:
That this is a cpu intensive and slow task!

Definition at line 1882 of file qgsrasterlayer.cpp.

References RasterBandStats::bandName, RasterBandStats::bandNoInt, RasterBandStats::colorTable, RasterBandStats::elementCountInt, gdalDataset, QgsMapLayer::layerName, RasterBandStats::maxValDouble, RasterBandStats::meanDouble, RasterBandStats::minValDouble, noDataValueDouble, RasterBandStats::rangeDouble, rasterLayerType, rasterStatsVector, rasterYDimInt, readValue(), QgsMapLayer::setProgress(), QgsMapLayer::setStatus(), RasterBandStats::statsGatheredFlag, RasterBandStats::stdDevDouble, RasterBandStats::sumDouble, and RasterBandStats::sumSqrDevDouble.

Referenced by drawPalettedSingleBandGray(), drawPalettedSingleBandPseudoColor(), drawSingleBandGray(), drawSingleBandPseudoColor(), getMetadata(), getRasterBandStats(), and readFile().

{
  //check if we have received a valid band number
  if ((gdalDataset->GetRasterCount() < theBandNoInt) && rasterLayerType != PALETTE)
  {
    //invalid band id, return nothing
    RasterBandStats myNullReturnStats;
    return myNullReturnStats;
  }
  if (rasterLayerType == PALETTE && (theBandNoInt > 3))
  {
    //invalid band id, return nothing
    RasterBandStats myNullReturnStats;
    return myNullReturnStats;
  }
  //check if we have previously gathered stats for this band...

  RasterBandStats myRasterBandStats = rasterStatsVector[theBandNoInt - 1];
  myRasterBandStats.bandNoInt = theBandNoInt;

  //dont bother with this if we already have stats
  if (myRasterBandStats.statsGatheredFlag)
  {
    return myRasterBandStats;
  }
  //onyl print message if we are actually gathering the stats
  emit setStatus(QString("Retrieving stats for ")+layerName);
  qApp->processEvents();
#ifdef QGISDEBUG
  std::cout << "QgsRasterLayer::retrieve stats for band " << theBandNoInt << std::endl;
#endif


  GDALRasterBand *myGdalBand = gdalDataset->GetRasterBand(theBandNoInt);
  QString myColorInterpretation = GDALGetColorInterpretationName(myGdalBand->GetColorInterpretation());

  //declare a colorTable to hold a palette - will only be used if the layer color interp is palette ???
  //get the palette colour table
  QgsColorTable *myColorTable = &(myRasterBandStats.colorTable);
  if (rasterLayerType == PALETTE)
  {

    //override the band name - palette images are really only one band
    //so we are faking three band behaviour
    switch (theBandNoInt)
    {
      // a "Red" layer
        case 1:
            myRasterBandStats.bandName = redTranslatedQString;
            break;
        case 2:
            myRasterBandStats.bandName = blueTranslatedQString;
            break;
        case 3:
            myRasterBandStats.bandName = greenTranslatedQString;
            break;
        default:
            //invalid band id so return
            RasterBandStats myNullReturnStats;
            return myNullReturnStats;
            break;
    }
  }
  else if (rasterLayerType==GRAY_OR_UNDEFINED)
  {
    myRasterBandStats.bandName = myColorInterpretation;
  }
  else //rasterLayerType is MULTIBAND
  {
    //do nothing
  }

  // XXX this sets the element count to a sensible value; but then you ADD to
  // XXX it later while iterating through all the pixels?
  //myRasterBandStats.elementCountInt = rasterXDimInt * rasterYDimInt;

  myRasterBandStats.elementCountInt = 0; // because we'll be counting only VALID pixels later

  emit setStatus(QString("Calculating stats for ")+layerName);
  //reset the main app progress bar
  emit setProgress(0,0);

  // let the user know we're going to possibly be taking a while
  QApplication::setOverrideCursor(QCursor(Qt::WaitCursor));

  GDALDataType myDataType = myGdalBand->GetRasterDataType();

  int  myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize;
  myGdalBand->GetBlockSize( &myXBlockSize, &myYBlockSize );

  myNXBlocks = (myGdalBand->GetXSize() + myXBlockSize - 1) / myXBlockSize;
  myNYBlocks = (myGdalBand->GetYSize() + myYBlockSize - 1) / myYBlockSize;

  void *myData = CPLMalloc ( myXBlockSize * myYBlockSize * GDALGetDataTypeSize(myDataType)/8 );

  // unfortunately we need to make two passes through the data to calculate stddev
  bool myFirstIterationFlag = true;

  // Comparison value for equality; i.e., we shouldn't directly compare two
  // floats so it's better to take their difference and see if they're within
  // a certain range -- in this case twenty times the smallest value that
  // doubles can take for the current system.  (Yes, 20 was arbitrary.)
  double myPrecision = std::numeric_limits<double>::epsilon() * 20;
  int success;
  double GDALminimum = myGdalBand->GetMinimum( &success );

  if ( ! success )
  {
      std::cerr << __FILE__ << " : " << __LINE__
                << " myGdalBand->GetMinimum() failed\n";
  }

  double GDALmaximum = myGdalBand->GetMaximum( &success );

  if ( ! success )
  {
      std::cerr << __FILE__ << " : " << __LINE__
                << " myGdalBand->GetMaximum() failed\n";
  }

  double GDALnodata = myGdalBand->GetNoDataValue( &success );

  if ( ! success )
  {
      std::cerr << __FILE__ << " : " << __LINE__
                << " myGdalBand->GetNoDataValue() failed\n";
  }

  std::cerr << "GDALminium:\t" << GDALminimum << "\n";
  std::cerr << "GDALmaxium:\t" << GDALmaximum << "\n";
  std::cerr << "GDALnodata:\t" << GDALnodata  << "\n";

  double GDALrange[2];          // calculated min/max, as opposed to the
                                // dataset provided

  GDALComputeRasterMinMax( myGdalBand, 0, GDALrange );

  std::cerr << "approximate computed GDALminium:\t" << GDALrange[0] << "\n";
  std::cerr << "approximate computed GDALmaxium:\t" << GDALrange[1] << "\n";

  GDALComputeRasterMinMax( myGdalBand, 1, GDALrange );

  std::cerr << "exactly computed GDALminium:\t" << GDALrange[0] << "\n";
  std::cerr << "exactly computed GDALmaxium:\t" << GDALrange[1] << "\n";


  for( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
  {
    emit setProgress ( iYBlock, myNYBlocks * 2 );

    for( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
    {
      int  nXValid, nYValid;
      myGdalBand->ReadBlock( iXBlock, iYBlock, myData );

      // Compute the portion of the block that is valid
      // for partial edge blocks.
      if( (iXBlock+1) * myXBlockSize > myGdalBand->GetXSize() )
          nXValid = myGdalBand->GetXSize() - iXBlock * myXBlockSize;
      else
          nXValid = myXBlockSize;

      if( (iYBlock+1) * myYBlockSize > myGdalBand->GetYSize() )
          nYValid = myGdalBand->GetYSize() - iYBlock * myYBlockSize;
      else
          nYValid = myYBlockSize;

      // Collect the histogram counts.
      for( int iY = 0; iY < nYValid; iY++ )
      {
        for( int iX = 0; iX < nXValid; iX++ )
        {
           double myDouble = readValue ( myData, myDataType, iX + iY * myXBlockSize );

           if ( fabs(myDouble - noDataValueDouble) < myPrecision ||
                myDouble < GDALminimum )
           {
               continue; // NULL
           }

           //get the nth element from the current row
           if (myColorInterpretation == "Palette") // dont translate this its a gdal string
           {
              //this is a palette layer so red / green / blue 'layers are 'virtual'
              //in that we need to obtain the palette entry and then get the r,g or g
              //component from that palette entry

              int c1, c2, c3;
              bool found = myColorTable->color ( myDouble, &c1, &c2, &c3 );
              if ( !found ) continue;

              //check for alternate color mappings
              switch (theBandNoInt)
              {
                  case 1:
                      myDouble = c1;
                      break;
                  case 2:
                      myDouble = c2;
                      break;
                  case 3:
                      myDouble = c3;
                      break;
              }
           }

           //only use this element if we have a non null element
           if (myFirstIterationFlag)
           {
             //this is the first iteration so initialise vars
             myFirstIterationFlag = false;
             myRasterBandStats.minValDouble = myDouble;
             myRasterBandStats.maxValDouble = myDouble;
           }               //end of true part for first iteration check
           else
           {
             //this is done for all subsequent iterations
             if (myDouble < myRasterBandStats.minValDouble)
             {
               myRasterBandStats.minValDouble = myDouble;
             }
             if (myDouble > myRasterBandStats.maxValDouble)
             {
               myRasterBandStats.maxValDouble = myDouble;
             }
             //only increment the running total if it is not a nodata value
             if (myDouble != noDataValueDouble)
             {
               myRasterBandStats.sumDouble += myDouble;
               ++myRasterBandStats.elementCountInt;
             }
           }               //end of false part for first iteration check
        }
      }
    }                       //end of column wise loop
  }                           //end of row wise loop


  //end of first pass through data now calculate the range
  myRasterBandStats.rangeDouble = myRasterBandStats.maxValDouble - myRasterBandStats.minValDouble;
  //calculate the mean
  myRasterBandStats.meanDouble = myRasterBandStats.sumDouble / myRasterBandStats.elementCountInt;

  //for the second pass we will get the sum of the squares / mean
  for( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
  {
    emit setProgress ( iYBlock+myNYBlocks, myNYBlocks * 2 );

    for( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
    {
      int  nXValid, nYValid;

      myGdalBand->ReadBlock( iXBlock, iYBlock, myData );

      // Compute the portion of the block that is valid
      // for partial edge blocks.
      if( (iXBlock+1) * myXBlockSize > myGdalBand->GetXSize() )
          nXValid = myGdalBand->GetXSize() - iXBlock * myXBlockSize;
      else
          nXValid = myXBlockSize;

      if( (iYBlock+1) * myYBlockSize > myGdalBand->GetYSize() )
          nYValid = myGdalBand->GetYSize() - iYBlock * myYBlockSize;
      else
          nYValid = myYBlockSize;

      // Collect the histogram counts.
      for( int iY = 0; iY < nYValid; iY++ )
      {
        for( int iX = 0; iX < nXValid; iX++ )
        {
           double myDouble = readValue ( myData, myDataType, iX + iY * myXBlockSize );

           if ( fabs(myDouble - noDataValueDouble) < myPrecision ||
                myDouble < GDALminimum )
           {
               continue; // NULL
           }

           //get the nth element from the current row
           if (myColorInterpretation == "Palette") // dont translate this its a gdal string
           {
              //this is a palette layer so red / green / blue 'layers are 'virtual'
              //in that we need to obtain the palette entry and then get the r,g or g
              //component from that palette entry

              int c1, c2, c3;
              bool found = myColorTable->color ( myDouble, &c1, &c2, &c3 );
              if ( !found ) continue;

              //check for alternate color mappings
              switch (theBandNoInt)
              {
                  case 1:
                      myDouble = c1;
                      break;
                  case 2:
                      myDouble = c2;
                      break;
                  case 3:
                      myDouble = c3;
                      break;
              }
           }

           myRasterBandStats.sumSqrDevDouble += static_cast < double >
                                                (pow(myDouble - myRasterBandStats.meanDouble, 2));
        }
      }
    }                       //end of column wise loop
  }                           //end of row wise loop

  //divide result by sample size - 1 and get square root to get stdev
  myRasterBandStats.stdDevDouble = static_cast < double >(sqrt(myRasterBandStats.sumSqrDevDouble /
              (myRasterBandStats.elementCountInt - 1)));

#ifdef QGISDEBUG
  std::cout << "************ STATS **************" << std::endl;
  std::cout << "NULL   = " << noDataValueDouble << std::endl;
  std::cout << "MIN    = " << myRasterBandStats.minValDouble << std::endl;
  std::cout << "MAX    = " << myRasterBandStats.maxValDouble << std::endl;
  std::cout << "RANGE  = " << myRasterBandStats.rangeDouble << std::endl;
  std::cout << "MEAN   = " << myRasterBandStats.meanDouble << std::endl;
  std::cout << "STDDEV = " << myRasterBandStats.stdDevDouble << std::endl;
#endif

  CPLFree(myData);
  myRasterBandStats.statsGatheredFlag = true;

  //add this band to the class stats map
  rasterStatsVector[theBandNoInt - 1] = myRasterBandStats;
  emit setProgress(rasterYDimInt, rasterYDimInt); //reset progress
  QApplication::restoreOverrideCursor(); //restore the cursor

  return myRasterBandStats;

} // QgsRasterLayer::getRasterBandStats


Generated by  Doxygen 1.6.0   Back to index