Skip to content

Commit 7b05793

Browse files
blazekblazek
committed
GRASS stats upgrade
1 parent eb93568 commit 7b05793

File tree

5 files changed

+106
-21
lines changed

5 files changed

+106
-21
lines changed

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

Lines changed: 43 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include <stdlib.h>
2828
#include <stdio.h>
2929
#include <string.h>
30+
#include <float.h>
3031
#include <grass/gis.h>
3132
#include <grass/raster.h>
3233
#include <grass/display.h>
@@ -35,7 +36,7 @@
3536
int main( int argc, char **argv )
3637
{
3738
struct GModule *module;
38-
struct Option *info_opt, *rast_opt, *vect_opt, *coor_opt;
39+
struct Option *info_opt, *rast_opt, *vect_opt, *coor_opt, *north_opt, *south_opt, *east_opt, *west_opt, *rows_opt, *cols_opt;
3940
struct Cell_head window;
4041

4142
/* Initialize the GIS calls */
@@ -63,6 +64,30 @@ int main( int argc, char **argv )
6364
coor_opt->type = TYPE_DOUBLE;
6465
coor_opt->multiple = YES;
6566

67+
north_opt = G_define_option();
68+
north_opt->key = "north";
69+
north_opt->type = TYPE_STRING;
70+
71+
south_opt = G_define_option();
72+
south_opt->key = "south";
73+
south_opt->type = TYPE_STRING;
74+
75+
east_opt = G_define_option();
76+
east_opt->key = "east";
77+
east_opt->type = TYPE_STRING;
78+
79+
west_opt = G_define_option();
80+
west_opt->key = "west";
81+
west_opt->type = TYPE_STRING;
82+
83+
rows_opt = G_define_option();
84+
rows_opt->key = "rows";
85+
rows_opt->type = TYPE_INTEGER;
86+
87+
cols_opt = G_define_option();
88+
cols_opt->key = "cols";
89+
cols_opt->type = TYPE_INTEGER;
90+
6691
if ( G_parser( argc, argv ) )
6792
exit( EXIT_FAILURE );
6893

@@ -249,13 +274,22 @@ int main( int argc, char **argv )
249274
int row, col;
250275
void *ptr;
251276
double val;
277+
double min = DBL_MAX;
278+
double max = -DBL_MAX;
252279
double sum = 0; // sum of values
253280
int count = 0; // count of non null values
254281
double mean = 0;
255282
double squares_sum = 0; // sum of squares
256283
double stdev = 0; // standard deviation
257284

258285
G_get_cellhd( rast_opt->answer, "", &window );
286+
window.north = atof( north_opt->answer );
287+
window.south = atof( south_opt->answer );
288+
window.east = atof( east_opt->answer );
289+
window.west = atof( west_opt->answer );
290+
window.rows = ( int ) atoi( rows_opt->answer );
291+
window.cols = ( int ) atoi( cols_opt->answer );
292+
259293
G_set_window( &window );
260294
fd = G_open_cell_old( rast_opt->answer, "" );
261295

@@ -306,6 +340,8 @@ int main( int argc, char **argv )
306340
}
307341
if ( ! G_is_null_value( ptr, rast_type ) )
308342
{
343+
if ( val < min ) min = val;
344+
if ( val > max ) max = val;
309345
sum += val;
310346
count++;
311347
squares_sum += pow( val, 2 );
@@ -316,11 +352,13 @@ int main( int argc, char **argv )
316352
squares_sum -= count * pow( mean, 2 );
317353
stdev = sqrt( squares_sum / ( count - 1 ) );
318354

319-
fprintf( stdout, "SUM:%e\n", sum );
320-
fprintf( stdout, "MEAN:%e\n", mean );
355+
fprintf( stdout, "MIN:%.17e\n", min );
356+
fprintf( stdout, "MAX:%.17e\n", max );
357+
fprintf( stdout, "SUM:%.17e\n", sum );
358+
fprintf( stdout, "MEAN:%.17e\n", mean );
321359
fprintf( stdout, "COUNT:%d\n", count );
322-
fprintf( stdout, "STDEV:%e\n", stdev );
323-
fprintf( stdout, "SQSUM:%e\n", squares_sum );
360+
fprintf( stdout, "STDEV:%.17e\n", stdev );
361+
fprintf( stdout, "SQSUM:%.17e\n", squares_sum );
324362

325363
G_close_cell( fd );
326364
}

‎src/providers/grass/qgsgrass.cpp‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1213,7 +1213,7 @@ QByteArray GRASS_LIB_EXPORT QgsGrass::runModule( QString gisdbase, QString locat
12131213
}
12141214

12151215
QString GRASS_LIB_EXPORT QgsGrass::getInfo( QString info, QString gisdbase, QString location,
1216-
QString mapset, QString map, MapType type, double x, double y, int timeOut )
1216+
QString mapset, QString map, MapType type, double x, double y, QgsRectangle extent, int sampleRows, int sampleCols, int timeOut )
12171217
{
12181218
QgsDebugMsg( QString( "gisdbase = %1 location = %2" ).arg( gisdbase ).arg( location ) );
12191219

@@ -1243,6 +1243,15 @@ QString GRASS_LIB_EXPORT QgsGrass::getInfo( QString info, QString gisdbase, QStr
12431243
{
12441244
arguments.append( QString( "coor=%1,%2" ).arg( x ).arg( y ) );
12451245
}
1246+
if ( info == "stats" )
1247+
{
1248+
arguments.append( QString( "north=%1" ).arg( extent.yMaximum() ) );
1249+
arguments.append( QString( "south=%1" ).arg( extent.yMinimum() ) );
1250+
arguments.append( QString( "east=%1" ).arg( extent.xMaximum() ) );
1251+
arguments.append( QString( "west=%1" ).arg( extent.xMinimum() ) );
1252+
arguments.append( QString( "rows=%1" ).arg( sampleRows ) );
1253+
arguments.append( QString( "cols=%1" ).arg( sampleCols ) );
1254+
}
12461255

12471256
QByteArray data = QgsGrass::runModule( gisdbase, location, cmd, arguments, timeOut );
12481257
QgsDebugMsg( data );
@@ -1358,15 +1367,15 @@ void GRASS_LIB_EXPORT QgsGrass::size( QString gisdbase, QString location, QStrin
13581367
QgsDebugMsg( QString( "raster size = %1 %2" ).arg( *cols ).arg( *rows ) );
13591368
}
13601369

1361-
QHash<QString, QString> GRASS_LIB_EXPORT QgsGrass::info( QString gisdbase, QString location, QString mapset, QString map, MapType type, QString info, int timeOut )
1370+
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 )
13621371
{
13631372
QgsDebugMsg( QString( "gisdbase = %1 location = %2" ).arg( gisdbase ).arg( location ) );
13641373
QHash<QString, QString> inf;
13651374

13661375
try
13671376
{
13681377
//QString str = QgsGrass::getInfo( "info", gisdbase, location, mapset, map, type );
1369-
QString str = QgsGrass::getInfo( info, gisdbase, location, mapset, map, type, 0, 0, timeOut );
1378+
QString str = QgsGrass::getInfo( info, gisdbase, location, mapset, map, type, 0, 0, extent, sampleRows, sampleCols, timeOut );
13701379
QgsDebugMsg( str );
13711380
QStringList list = str.split( "\n" );
13721381
for ( int i = 0; i < list.size(); i++ )

‎src/providers/grass/qgsgrass.h‎

Lines changed: 17 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ extern "C"
2525

2626
#include <stdexcept>
2727
#include "qgsexception.h"
28+
#include <qgsrectangle.h>
2829
#include <QProcess>
2930
#include <QString>
3031
#include <QMap>
@@ -192,9 +193,21 @@ class QgsGrass
192193
// ! Run a GRASS module in any gisdbase/location
193194
static GRASS_LIB_EXPORT QByteArray runModule( QString gisdbase, QString location, QString module, QStringList arguments, int timeOut = 30000 );
194195

195-
// ! Get info string from qgis.g.info module
196+
/** \brief Get info string from qgis.g.info module
197+
* @param info info type
198+
* @gisdbase GISBASE path
199+
* @location location name
200+
* @mapset mapset name
201+
* @map map name
202+
* @type map type
203+
* @x x coordinate for query
204+
* @y y coordinate for query
205+
* @extent extent for statistics
206+
* @sampleSize sample size for statistics
207+
* @timeOut timeout
208+
*/
196209
static GRASS_LIB_EXPORT QString getInfo( QString info, QString gisdbase,
197-
QString location, QString mapset = "", QString map = "", MapType type = None, double x = 0.0, double y = 0.0, int timeOut = 30000 );
210+
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 );
198211

199212
// ! Get location projection
200213
static GRASS_LIB_EXPORT QgsCoordinateReferenceSystem crs( QString gisdbase, QString location );
@@ -211,8 +224,9 @@ class QgsGrass
211224
QString mapset, QString map, int *cols, int *rows );
212225

213226
// ! Get raster info, info is either 'info' or 'stats'
227+
// extent and sampleSize are stats options
214228
static GRASS_LIB_EXPORT QHash<QString, QString> info( QString gisdbase, QString location,
215-
QString mapset, QString map, MapType type, QString info = "info", int timeOut = 30000 );
229+
QString mapset, QString map, MapType type, QString info = "info", QgsRectangle extent = QgsRectangle(), int sampleRows = 0, int sampleCols = 0, int timeOut = 30000 );
216230

217231
// ! List of Color
218232
static GRASS_LIB_EXPORT QList<QgsGrass::Color> colors( QString gisdbase, QString location,

‎src/providers/grass/qgsgrassrasterprovider.cpp‎

Lines changed: 30 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -301,16 +301,31 @@ double QgsGrassRasterProvider::maximumValue( int bandNo ) const
301301
return mInfo["MAX_VALUE"].toDouble();
302302
}
303303

304-
QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo )
304+
QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo, int theStats, const QgsRectangle & theExtent, int theSampleSize )
305305
{
306+
QgsDebugMsg( QString( "theBandNo = %1 theSampleSize = %2" ).arg( theBandNo ).arg( theSampleSize ) );
306307
QgsRasterBandStats myRasterBandStats;
307-
myRasterBandStats.bandName = generateBandName( theBandNo );
308-
myRasterBandStats.bandNumber = theBandNo;
308+
initStatistics( myRasterBandStats, theBandNo, theStats, theExtent, theSampleSize );
309+
310+
foreach ( QgsRasterBandStats stats, mStatistics )
311+
{
312+
if ( stats.contains( myRasterBandStats ) )
313+
{
314+
QgsDebugMsg( "Using cached statistics." );
315+
return stats;
316+
}
317+
}
318+
319+
QgsRectangle extent = myRasterBandStats.extent;
320+
321+
int sampleRows = myRasterBandStats.height;
322+
int sampleCols = myRasterBandStats.width;
309323

310324
// With stats we have to be careful about timeout, empirical value,
311325
// 0.001 / cell should be sufficient using 0.005 to be sure + constant (ms)
312326
int timeout = 30000 + 0.005 * xSize() * ySize();
313-
QHash<QString, QString> info = QgsGrass::info( mGisdbase, mLocation, mMapset, mMapName, QgsGrass::Raster, "stats", timeout );
327+
328+
QHash<QString, QString> info = QgsGrass::info( mGisdbase, mLocation, mMapset, mMapName, QgsGrass::Raster, "stats", extent, sampleRows, sampleCols, timeout );
314329

315330
if ( info.isEmpty() )
316331
{
@@ -319,18 +334,24 @@ QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo )
319334

320335
myRasterBandStats.sum = info["SUM"].toDouble();
321336
myRasterBandStats.elementCount = info["COUNT"].toInt();
322-
myRasterBandStats.minimumValue = minimumValue( theBandNo );
323-
myRasterBandStats.maximumValue = maximumValue( theBandNo );
324-
myRasterBandStats.range = maximumValue( theBandNo ) - minimumValue( theBandNo );
337+
myRasterBandStats.minimumValue = info["MIN"].toDouble();
338+
myRasterBandStats.maximumValue = info["MAX"].toDouble();
339+
myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
325340
myRasterBandStats.sumOfSquares = info["SQSUM"].toDouble();
326341
myRasterBandStats.mean = info["MEAN"].toDouble();
327342
myRasterBandStats.stdDev = info["STDEV"].toDouble();
328343

329-
QgsDebugMsg( QString( "sum = %1" ).arg( myRasterBandStats.sum ) );
344+
QgsDebugMsg( QString( "min = %1" ).arg( myRasterBandStats.minimumValue ) );
345+
QgsDebugMsg( QString( "max = %1" ).arg( myRasterBandStats.maximumValue ) );
330346
QgsDebugMsg( QString( "count = %1" ).arg( myRasterBandStats.elementCount ) );
331347
QgsDebugMsg( QString( "stdev = %1" ).arg( myRasterBandStats.stdDev ) );
332348

333-
myRasterBandStats.statsGathered = true;
349+
myRasterBandStats.statsGathered = QgsRasterBandStats::Min | QgsRasterBandStats::Max |
350+
QgsRasterBandStats::Range | QgsRasterBandStats::Mean |
351+
QgsRasterBandStats::Sum | QgsRasterBandStats::SumOfSquares |
352+
QgsRasterBandStats::StdDev;
353+
354+
mStatistics.append( myRasterBandStats );
334355
return myRasterBandStats;
335356
}
336357

‎src/providers/grass/qgsgrassrasterprovider.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -213,7 +213,10 @@ class QgsGrassRasterProvider : public QgsRasterDataProvider
213213
double minimumValue( int bandNo )const;
214214
double maximumValue( int bandNo )const;
215215

216-
QgsRasterBandStats bandStatistics( int theBandNo );
216+
QgsRasterBandStats bandStatistics( int theBandNo,
217+
int theStats = QgsRasterBandStats::All,
218+
const QgsRectangle & theExtent = QgsRectangle(),
219+
int theSampleSize = 0 );
217220

218221
QList<QgsColorRampShader::ColorRampItem> colorTable( int bandNo )const;
219222

0 commit comments

Comments
 (0)