Skip to content
Permalink
Browse files

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

  • Loading branch information
blazek committed Sep 5, 2012
1 parent 6662601 commit eb93568f760e14797dfb0bc46a468224b0332192
@@ -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 );
}
@@ -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++ )
@@ -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,
@@ -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 );
@@ -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 );

0 comments on commit eb93568

Please sign in to comment.
You can’t perform that action at this time.