src/providers/grass/qgis.g.info.c

@@ -27,6 +27,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <float.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/display.h>
@@ -35,7 +36,7 @@
 int main( int argc, char **argv )
 {
   struct GModule *module;
-  struct Option *info_opt, *rast_opt, *vect_opt, *coor_opt;
+  struct Option *info_opt, *rast_opt, *vect_opt, *coor_opt, *north_opt, *south_opt, *east_opt, *west_opt, *rows_opt, *cols_opt;
   struct Cell_head window;
 
   /* Initialize the GIS calls */
@@ -48,7 +49,7 @@ int main( int argc, char **argv )
   info_opt->key = "info";
   info_opt->type = TYPE_STRING;
   info_opt->description = "info key";
-  info_opt->options = "proj,window,size,query,info,colors";
+  info_opt->options = "proj,window,size,query,info,colors,stats";
 
   rast_opt = G_define_standard_option( G_OPT_R_INPUT );
   rast_opt->key = "rast";
@@ -63,6 +64,30 @@ int main( int argc, char **argv )
   coor_opt->type = TYPE_DOUBLE;
   coor_opt->multiple = YES;
 
+  north_opt = G_define_option();
+  north_opt->key = "north";
+  north_opt->type = TYPE_STRING;
+
+  south_opt = G_define_option();
+  south_opt->key = "south";
+  south_opt->type = TYPE_STRING;
+
+  east_opt = G_define_option();
+  east_opt->key = "east";
+  east_opt->type = TYPE_STRING;
+
+  west_opt = G_define_option();
+  west_opt->key = "west";
+  west_opt->type = TYPE_STRING;
+
+  rows_opt = G_define_option();
+  rows_opt->key = "rows";
+  rows_opt->type = TYPE_INTEGER;
+
+  cols_opt = G_define_option();
+  cols_opt->key = "cols";
+  cols_opt->type = TYPE_INTEGER;
+
   if ( G_parser( argc, argv ) )
     exit( EXIT_FAILURE );
 
@@ -237,6 +262,111 @@ int main( int argc, char **argv )
       G_fatal_error( "Not yet supported" );
     }
   }
+  else if ( strcmp( "stats", info_opt->answer ) == 0 )
+  {
+    if ( rast_opt->answer )
+    {
+      int fd;
+      RASTER_MAP_TYPE rast_type;
+      DCELL *dcell;
+      CELL *cell;
+      int ncols, nrows;
+      int row, col;
+      void *ptr;
+      double val;
+      double min = DBL_MAX;
+      double max = -DBL_MAX;
+      double sum = 0; // sum of values
+      int count = 0; // count of non null values
+      double mean = 0;
+      double squares_sum = 0; // sum of squares
+      double stdev = 0; // standard deviation
+
+      G_get_cellhd( rast_opt->answer, "", &window );
+      window.north = atof( north_opt->answer );
+      window.south = atof( south_opt->answer );
+      window.east = atof( east_opt->answer );
+      window.west = atof( west_opt->answer );
+      window.rows = ( int ) atoi( rows_opt->answer );
+      window.cols = ( int ) atoi( cols_opt->answer );
+
+      G_set_window( &window );
+      fd = G_open_cell_old( rast_opt->answer, "" );
+
+      ncols = G_window_cols();
+      nrows = G_window_rows();
+
+#if defined(GRASS_VERSION_MAJOR) && defined(GRASS_VERSION_MINOR) && \
+    ( ( GRASS_VERSION_MAJOR == 6 && GRASS_VERSION_MINOR > 2 ) || GRASS_VERSION_MAJOR > 6 )
+      rast_type = G_get_raster_map_type( fd );
+#else
+      rast_type = G_raster_map_type( rast_opt->answer, "" );
+#endif
+      cell = G_allocate_c_raster_buf();
+      dcell = G_allocate_d_raster_buf();
+
+      // Calc stats is very slow for large rasters -> prefer optimization for speed over
+      // code length and readability (which is not currently true)
+      for ( row = 0; row < nrows; row++ )
+      {
+        if ( rast_type == CELL_TYPE )
+        {
+          if ( G_get_c_raster_row( fd, cell, row ) < 0 )
+          {
+            G_fatal_error(( "Unable to read raster map <%s> row %d" ),
+                          rast_opt->answer, row );
+          }
+        }
+        else
+        {
+          if ( G_get_d_raster_row( fd, dcell, row ) < 0 )
+          {
+            G_fatal_error(( "Unable to read raster map <%s> row %d" ),
+                          rast_opt->answer, row );
+          }
+        }
+
+        for ( col = 0; col < ncols; col++ )
+        {
+          if ( rast_type == CELL_TYPE )
+          {
+            val = cell[col];
+            ptr = &( cell[col] );
+          }
+          else
+          {
+            val = dcell[col];
+            ptr = &( dcell[col] );
+          }
+          if ( ! G_is_null_value( ptr, rast_type ) )
+          {
+            if ( val < min ) min = val;
+            if ( val > max ) max = val;
+            sum += val;
+            count++;
+            squares_sum += pow( val, 2 );
+          }
+        }
+      }
+      mean = sum / count;
+      squares_sum -= count * pow( mean, 2 );
+      stdev = sqrt( squares_sum / ( count - 1 ) );
+
+      fprintf( stdout, "MIN:%.17e\n", min );
+      fprintf( stdout, "MAX:%.17e\n", max );
+      fprintf( stdout, "SUM:%.17e\n", sum );
+      fprintf( stdout, "MEAN:%.17e\n", mean );
+      fprintf( stdout, "COUNT:%d\n", count );
+      fprintf( stdout, "STDEV:%.17e\n", stdev );
+      fprintf( stdout, "SQSUM:%.17e\n", squares_sum );
+
+      G_close_cell( fd );
+    }
+    else if ( vect_opt->answer )
+    {
+      G_fatal_error( "Not yet supported" );
+    }
+  }
 
   exit( EXIT_SUCCESS );
 }

src/providers/grass/qgsgrass.cpp

@@ -1189,14 +1189,14 @@ QProcess GRASS_LIB_EXPORT *QgsGrass::startModule( QString gisdbase, QString loca
 }
 
 QByteArray GRASS_LIB_EXPORT QgsGrass::runModule( QString gisdbase, QString location,
-    QString module, QStringList arguments )
+    QString module, QStringList arguments, int timeOut )
 {
-  QgsDebugMsg( QString( "gisdbase = %1 location = %2" ).arg( gisdbase ).arg( location ) );
+  QgsDebugMsg( QString( "gisdbase = %1 location = %2 timeOut = %3" ).arg( gisdbase ).arg( location ).arg( timeOut ) );
 
   QTemporaryFile gisrcFile;
   QProcess *process = QgsGrass::startModule( gisdbase, location, module, arguments, gisrcFile );
 
-  if ( !process->waitForFinished()
+  if ( !process->waitForFinished( timeOut )
        || ( process->exitCode() != 0 && process->exitCode() != 255 ) )
   {
     QgsDebugMsg( "process->exitCode() = " + QString::number( process->exitCode() ) );
@@ -1213,7 +1213,7 @@ QByteArray GRASS_LIB_EXPORT QgsGrass::runModule( QString gisdbase, QString locat
 }
 
 QString GRASS_LIB_EXPORT QgsGrass::getInfo( QString info, QString gisdbase, QString location,
-    QString mapset, QString map, MapType type, double x, double y )
+    QString mapset, QString map, MapType type, double x, double y, QgsRectangle extent, int sampleRows, int sampleCols, int timeOut )
 {
   QgsDebugMsg( QString( "gisdbase = %1 location = %2" ).arg( gisdbase ).arg( location ) );
 
@@ -1243,8 +1243,17 @@ QString GRASS_LIB_EXPORT QgsGrass::getInfo( QString info, QString gisdbase, QStr
   {
     arguments.append( QString( "coor=%1,%2" ).arg( x ).arg( y ) );
   }
+  if ( info == "stats" )
+  {
+    arguments.append( QString( "north=%1" ).arg( extent.yMaximum() ) );
+    arguments.append( QString( "south=%1" ).arg( extent.yMinimum() ) );
+    arguments.append( QString( "east=%1" ).arg( extent.xMaximum() ) );
+    arguments.append( QString( "west=%1" ).arg( extent.xMinimum() ) );
+    arguments.append( QString( "rows=%1" ).arg( sampleRows ) );
+    arguments.append( QString( "cols=%1" ).arg( sampleCols ) );
+  }
 
-  QByteArray data =  QgsGrass::runModule( gisdbase, location, cmd, arguments );
+  QByteArray data =  QgsGrass::runModule( gisdbase, location, cmd, arguments, timeOut );
   QgsDebugMsg( data );
   return QString( data );
 }
@@ -1358,14 +1367,15 @@ void GRASS_LIB_EXPORT QgsGrass::size( QString gisdbase, QString location, QStrin
   QgsDebugMsg( QString( "raster size = %1 %2" ).arg( *cols ).arg( *rows ) );
 }
 
-QHash<QString, QString> GRASS_LIB_EXPORT QgsGrass::info( QString gisdbase, QString location, QString mapset, QString map, MapType type )
+QHash<QString, QString> GRASS_LIB_EXPORT QgsGrass::info( QString gisdbase, QString location, QString mapset, QString map, MapType type, QString info, QgsRectangle extent, int sampleRows, int sampleCols, int timeOut )
 {
   QgsDebugMsg( QString( "gisdbase = %1 location = %2" ).arg( gisdbase ).arg( location ) );
   QHash<QString, QString> inf;
 
   try
   {
-    QString str = QgsGrass::getInfo( "info", gisdbase, location, mapset, map, type );
+    //QString str = QgsGrass::getInfo( "info", gisdbase, location, mapset, map, type );
+    QString str = QgsGrass::getInfo( info, gisdbase, location, mapset, map, type, 0, 0, extent, sampleRows, sampleCols, timeOut );
     QgsDebugMsg( str );
     QStringList list = str.split( "\n" );
     for ( int i = 0; i < list.size(); i++ )

src/providers/grass/qgsgrass.h

@@ -25,6 +25,7 @@ extern "C"
 
 #include <stdexcept>
 #include "qgsexception.h"
+#include <qgsrectangle.h>
 #include <QProcess>
 #include <QString>
 #include <QMap>
@@ -190,11 +191,23 @@ class QgsGrass
     static GRASS_LIB_EXPORT QProcess *startModule( QString gisdbase, QString location, QString module, QStringList arguments, QTemporaryFile &gisrcFile );
 
     // ! Run a GRASS module in any gisdbase/location
-    static GRASS_LIB_EXPORT QByteArray runModule( QString gisdbase, QString location, QString module, QStringList arguments );
-
-    // ! Get info string from qgis.g.info module
+    static GRASS_LIB_EXPORT QByteArray runModule( QString gisdbase, QString location, QString module, QStringList arguments, int timeOut = 30000 );
+
+    /** \brief Get info string from qgis.g.info module
+     * @param info info type
+     * @gisdbase GISBASE path
+     * @location location name
+     * @mapset mapset name
+     * @map map name
+     * @type map type
+     * @x x coordinate for query
+     * @y y coordinate for query
+     * @extent extent for statistics
+     * @sampleSize sample size for statistics
+     * @timeOut timeout
+     */
     static GRASS_LIB_EXPORT QString getInfo( QString info, QString gisdbase,
-        QString location, QString mapset = "", QString map = "", MapType type = None, double x = 0.0, double y = 0.0 );
+        QString location, QString mapset = "", QString map = "", MapType type = None, double x = 0.0, double y = 0.0, QgsRectangle extent = QgsRectangle(), int sampleRows = 0, int sampleCols = 0, int timeOut = 30000 );
 
     // ! Get location projection
     static GRASS_LIB_EXPORT QgsCoordinateReferenceSystem crs( QString gisdbase, QString location );
@@ -210,9 +223,10 @@ class QgsGrass
     static GRASS_LIB_EXPORT void size( QString gisdbase, QString location,
                                        QString mapset, QString map, int *cols, int *rows );
 
-    // ! Get raster info
+    // ! Get raster info, info is either 'info' or 'stats'
+    //   extent and sampleSize are stats options
     static GRASS_LIB_EXPORT QHash<QString, QString> info( QString gisdbase, QString location,
-        QString mapset, QString map, MapType type );
+        QString mapset, QString map, MapType type, QString info = "info", QgsRectangle extent = QgsRectangle(), int sampleRows = 0, int sampleCols = 0, int timeOut = 30000 );
 
     // ! List of Color
     static GRASS_LIB_EXPORT QList<QgsGrass::Color> colors( QString gisdbase, QString location,

src/providers/grass/qgsgrassrasterprovider.cpp

@@ -301,6 +301,60 @@ double  QgsGrassRasterProvider::maximumValue( int bandNo ) const
   return mInfo["MAX_VALUE"].toDouble();
 }
 
+QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo, int theStats, const QgsRectangle & theExtent, int theSampleSize )
+{
+  QgsDebugMsg( QString( "theBandNo = %1 theSampleSize = %2" ).arg( theBandNo ).arg( theSampleSize ) );
+  QgsRasterBandStats myRasterBandStats;
+  initStatistics( myRasterBandStats, theBandNo, theStats, theExtent, theSampleSize );
+
+  foreach ( QgsRasterBandStats stats, mStatistics )
+  {
+    if ( stats.contains( myRasterBandStats ) )
+    {
+      QgsDebugMsg( "Using cached statistics." );
+      return stats;
+    }
+  }
+
+  QgsRectangle extent = myRasterBandStats.extent;
+
+  int sampleRows = myRasterBandStats.height;
+  int sampleCols = myRasterBandStats.width;
+
+  // With stats we have to be careful about timeout, empirical value,
+  // 0.001 / cell should be sufficient using 0.005 to be sure + constant (ms)
+  int timeout = 30000 + 0.005 * xSize() * ySize();
+
+  QHash<QString, QString> info = QgsGrass::info( mGisdbase, mLocation, mMapset, mMapName, QgsGrass::Raster, "stats", extent, sampleRows, sampleCols, timeout );
+
+  if ( info.isEmpty() )
+  {
+    return myRasterBandStats;
+  }
+
+  myRasterBandStats.sum = info["SUM"].toDouble();
+  myRasterBandStats.elementCount = info["COUNT"].toInt();
+  myRasterBandStats.minimumValue = info["MIN"].toDouble();
+  myRasterBandStats.maximumValue = info["MAX"].toDouble();
+  myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
+  myRasterBandStats.sumOfSquares = info["SQSUM"].toDouble();
+  myRasterBandStats.mean = info["MEAN"].toDouble();
+  myRasterBandStats.stdDev = info["STDEV"].toDouble();
+
+  QgsDebugMsg( QString( "min = %1" ).arg( myRasterBandStats.minimumValue ) );
+  QgsDebugMsg( QString( "max = %1" ).arg( myRasterBandStats.maximumValue ) );
+  QgsDebugMsg( QString( "count = %1" ).arg( myRasterBandStats.elementCount ) );
+  QgsDebugMsg( QString( "stdev = %1" ).arg( myRasterBandStats.stdDev ) );
+
+  myRasterBandStats.statsGathered = QgsRasterBandStats::Min | QgsRasterBandStats::Max |
+                                    QgsRasterBandStats::Range | QgsRasterBandStats::Mean |
+                                    QgsRasterBandStats::Sum | QgsRasterBandStats::SumOfSquares |
+                                    QgsRasterBandStats::StdDev;
+
+  mStatistics.append( myRasterBandStats );
+  return myRasterBandStats;
+}
+
 QList<QgsColorRampShader::ColorRampItem> QgsGrassRasterProvider::colorTable( int bandNo )const
 {
   Q_UNUSED( bandNo );

src/providers/grass/qgsgrassrasterprovider.h

@@ -213,6 +213,11 @@ class QgsGrassRasterProvider : public QgsRasterDataProvider
     double minimumValue( int bandNo )const;
     double maximumValue( int bandNo )const;
 
+    QgsRasterBandStats bandStatistics( int theBandNo,
+                                       int theStats = QgsRasterBandStats::All,
+                                       const QgsRectangle & theExtent = QgsRectangle(),
+                                       int theSampleSize = 0 );
+
     QList<QgsColorRampShader::ColorRampItem> colorTable( int bandNo )const;
 
     // void buildSupportedRasterFileFilter( QString & theFileFiltersString );