calc raster stats in grass module, time reduced to 20%, solves part of …

…#4862
qgis · Sep 5, 2012 · eb93568 · eb93568
1 parent 6662601
commit eb93568
Show file tree

Hide file tree

Showing 5 changed files with 140 additions and 12 deletions.
diff --git a/src/providers/grass/qgis.g.info.c b/src/providers/grass/qgis.g.info.c
@@ -48,7 +48,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";
@@ -237,6 +237,98 @@ 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 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 );
+      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 ) )
+          {
+            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, "SUM:%e\n", sum );
+      fprintf( stdout, "MEAN:%e\n", mean );
+      fprintf( stdout, "COUNT:%d\n", count );
+      fprintf( stdout, "STDEV:%e\n", stdev );
+      fprintf( stdout, "SQSUM:%e\n", squares_sum );
+
+      G_close_cell( fd );
+    }
+    else if ( vect_opt->answer )
+    {
+      G_fatal_error( "Not yet supported" );
+    }
+  }
 
   exit( EXIT_SUCCESS );
 }

diff --git a/src/providers/grass/qgsgrass.cpp b/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, int timeOut )
 {
   QgsDebugMsg( QString( "gisdbase = %1 location = %2" ).arg( gisdbase ).arg( location ) );
 
@@ -1244,7 +1244,7 @@ QString GRASS_LIB_EXPORT QgsGrass::getInfo( QString info, QString gisdbase, QStr
     arguments.append( QString( "coor=%1,%2" ).arg( x ).arg( y ) );
   }
 
-  QByteArray data =  QgsGrass::runModule( gisdbase, location, cmd, arguments );
+  QByteArray data =  QgsGrass::runModule( gisdbase, location, cmd, arguments, timeOut );
   QgsDebugMsg( data );
   return QString( data );
 }
@@ -1358,14 +1358,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, 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, timeOut );
     QgsDebugMsg( str );
     QStringList list = str.split( "\n" );
     for ( int i = 0; i < list.size(); i++ )

diff --git a/src/providers/grass/qgsgrass.h b/src/providers/grass/qgsgrass.h
@@ -190,11 +190,11 @@ 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 );
+    static GRASS_LIB_EXPORT QByteArray runModule( QString gisdbase, QString location, QString module, QStringList arguments, int timeOut = 30000 );
 
     // ! Get info string from qgis.g.info module
     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, int timeOut = 30000 );
 
     // ! Get location projection
     static GRASS_LIB_EXPORT QgsCoordinateReferenceSystem crs( QString gisdbase, QString location );
@@ -210,9 +210,9 @@ 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'
     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", int timeOut = 30000 );
 
     // ! List of Color
     static GRASS_LIB_EXPORT QList<QgsGrass::Color> colors( QString gisdbase, QString location,

diff --git a/src/providers/grass/qgsgrassrasterprovider.cpp b/src/providers/grass/qgsgrassrasterprovider.cpp
@@ -301,6 +301,39 @@ double  QgsGrassRasterProvider::maximumValue( int bandNo ) const
   return mInfo["MAX_VALUE"].toDouble();
 }
 
+QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo )
+{
+  QgsRasterBandStats myRasterBandStats;
+  myRasterBandStats.bandName = generateBandName( theBandNo );
+  myRasterBandStats.bandNumber = theBandNo;
+
+  // 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", timeout );
+
+  if ( info.isEmpty() )
+  {
+    return myRasterBandStats;
+  }
+
+  myRasterBandStats.sum = info["SUM"].toDouble();
+  myRasterBandStats.elementCount = info["COUNT"].toInt();
+  myRasterBandStats.minimumValue = minimumValue( theBandNo );
+  myRasterBandStats.maximumValue = maximumValue( theBandNo );
+  myRasterBandStats.range = maximumValue( theBandNo ) - minimumValue( theBandNo );
+  myRasterBandStats.sumOfSquares = info["SQSUM"].toDouble();
+  myRasterBandStats.mean = info["MEAN"].toDouble();
+  myRasterBandStats.stdDev = info["STDEV"].toDouble();
+
+  QgsDebugMsg( QString( "sum = %1" ).arg( myRasterBandStats.sum ) );
+  QgsDebugMsg( QString( "count = %1" ).arg( myRasterBandStats.elementCount ) );
+  QgsDebugMsg( QString( "stdev = %1" ).arg( myRasterBandStats.stdDev ) );
+
+  myRasterBandStats.statsGathered = true;
+  return myRasterBandStats;
+}
+
 QList<QgsColorRampShader::ColorRampItem> QgsGrassRasterProvider::colorTable( int bandNo )const
 {
   Q_UNUSED( bandNo );

diff --git a/src/providers/grass/qgsgrassrasterprovider.h b/src/providers/grass/qgsgrassrasterprovider.h
@@ -213,6 +213,8 @@ class QgsGrassRasterProvider : public QgsRasterDataProvider
     double minimumValue( int bandNo )const;
     double maximumValue( int bandNo )const;
 
+    QgsRasterBandStats bandStatistics( int theBandNo );
+
     QList<QgsColorRampShader::ColorRampItem> colorTable( int bandNo )const;
 
     // void buildSupportedRasterFileFilter( QString & theFileFiltersString );