Skip to content

Commit

Permalink
GRASS stats upgrade
Browse files Browse the repository at this point in the history
  • Loading branch information
@blazek
blazek committed Sep 5, 2012
1 parent eb93568 commit 7b05793
Show file tree
Hide file tree
Showing 5 changed files with 106 additions and 21 deletions.
48 changes: 43 additions & 5 deletions src/providers/grass/qgis.g.info.c
View file
Expand Up @@ -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>
Expand All @@ -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 */
Expand Down Expand Up @@ -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 );

Expand Down Expand Up @@ -249,13 +274,22 @@ int main( int argc, char **argv )
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, "" );

Expand Down Expand Up @@ -306,6 +340,8 @@ int main( int argc, char **argv )
}
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 );
Expand All @@ -316,11 +352,13 @@ int main( int argc, char **argv )
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, "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:%e\n", stdev );
fprintf( stdout, "SQSUM:%e\n", squares_sum );
fprintf( stdout, "STDEV:%.17e\n", stdev );
fprintf( stdout, "SQSUM:%.17e\n", squares_sum );

G_close_cell( fd );
}
Expand Down
15 changes: 12 additions & 3 deletions src/providers/grass/qgsgrass.cpp
View file
Expand Up @@ -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, int timeOut )
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 ) );

Expand Down Expand Up @@ -1243,6 +1243,15 @@ 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, timeOut );
QgsDebugMsg( data );
Expand Down Expand Up @@ -1358,15 +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, QString info, int timeOut )
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, 0, 0, timeOut );
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++ )
Expand Down
20 changes: 17 additions & 3 deletions src/providers/grass/qgsgrass.h
View file
Expand Up @@ -25,6 +25,7 @@ extern "C"

#include <stdexcept>
#include "qgsexception.h"
#include <qgsrectangle.h>
#include <QProcess>
#include <QString>
#include <QMap>
Expand Down Expand Up @@ -192,9 +193,21 @@ class QgsGrass
// ! Run a GRASS module in any gisdbase/location
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
/** \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, int timeOut = 30000 );
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 );
Expand All @@ -211,8 +224,9 @@ class QgsGrass
QString mapset, QString map, int *cols, int *rows );

// ! 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 info = "info", int timeOut = 30000 );
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,
Expand Down
39 changes: 30 additions & 9 deletions src/providers/grass/qgsgrassrasterprovider.cpp
View file
Expand Up @@ -301,16 +301,31 @@ double QgsGrassRasterProvider::maximumValue( int bandNo ) const
return mInfo["MAX_VALUE"].toDouble();
}

QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo )
QgsRasterBandStats QgsGrassRasterProvider::bandStatistics( int theBandNo, int theStats, const QgsRectangle & theExtent, int theSampleSize )
{
QgsDebugMsg( QString( "theBandNo = %1 theSampleSize = %2" ).arg( theBandNo ).arg( theSampleSize ) );
QgsRasterBandStats myRasterBandStats;
myRasterBandStats.bandName = generateBandName( theBandNo );
myRasterBandStats.bandNumber = theBandNo;
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", timeout );

QHash<QString, QString> info = QgsGrass::info( mGisdbase, mLocation, mMapset, mMapName, QgsGrass::Raster, "stats", extent, sampleRows, sampleCols, timeout );

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

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.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( "sum = %1" ).arg( myRasterBandStats.sum ) );
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 = true;
myRasterBandStats.statsGathered = QgsRasterBandStats::Min | QgsRasterBandStats::Max |
QgsRasterBandStats::Range | QgsRasterBandStats::Mean |
QgsRasterBandStats::Sum | QgsRasterBandStats::SumOfSquares |
QgsRasterBandStats::StdDev;

mStatistics.append( myRasterBandStats );
return myRasterBandStats;
}

Expand Down
5 changes: 4 additions & 1 deletion src/providers/grass/qgsgrassrasterprovider.h
View file
Expand Up @@ -213,7 +213,10 @@ class QgsGrassRasterProvider : public QgsRasterDataProvider
double minimumValue( int bandNo )const;
double maximumValue( int bandNo )const;

QgsRasterBandStats bandStatistics( int theBandNo );
QgsRasterBandStats bandStatistics( int theBandNo,
int theStats = QgsRasterBandStats::All,
const QgsRectangle & theExtent = QgsRectangle(),
int theSampleSize = 0 );

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

Expand Down

0 comments on commit 7b05793

Please sign in to comment.