sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-03-01 00:46:20 -05:00
Code Issues Packages Projects Releases Wiki Activity

GRASS stats upgrade

Browse Source
This commit is contained in:
Radim Blazek 2012-09-05 10:05:27 +02:00
parent eb93568f76
commit 7b05793b54
5 changed files with 106 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
src/providers/grass/qgis.g.info.c
View File

@ -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 */
@ -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 );
@ -249,6 +274,8 @@ 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;
@ -256,6 +283,13 @@ int main( int argc, char **argv )
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, "" );
@ -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 );
@ -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 );
}

15
src/providers/grass/qgsgrass.cpp
View File

@ -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 ) );
@ -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 );
@ -1358,7 +1367,7 @@ 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;
@ -1366,7 +1375,7 @@ QHash<QString, QString> GRASS_LIB_EXPORT QgsGrass::info( QString gisdbase, QStri
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++ )

20
src/providers/grass/qgsgrass.h
View File

@ -25,6 +25,7 @@ extern "C"
#include <stdexcept>
#include "qgsexception.h"
#include <qgsrectangle.h>
#include <QProcess>
#include <QString>
#include <QMap>
@ -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 );
@ -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,

39
src/providers/grass/qgsgrassrasterprovider.cpp
View File

@ -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() )
{
@ -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;
}

5
src/providers/grass/qgsgrassrasterprovider.h
View File

@ -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;