[Qgis-developer] Re: Access raster data from python plugins

[Ghislain Picard]

Hi,

Following the idea in my last email. I propose a function getValue in 
QgsRatserLayer before API freezing. More efficient than identify to 
retrieve value (return double instead of string) in a single band.

The same kind of function to return a piece of raster would be 
interesting too.

Regards,
Ghislain






double QgsRasterLayer::getValue(const QgsPoint& point, int iband)
{
   double x = point.x();
   double y = point.y();

   QgsDebugMsg( QString::number( x ) + ", " + QString::number( y ) );

   if ( x < mLayerExtent.xMin() || x > mLayerExtent.xMax() || y < 
mLayerExtent.yMin() || y > mLayerExtent.yMax() )
   {
     return nan();
   }
   else
   {
     /* Calculate the row / column where the point falls */
     double xres = ( mLayerExtent.xMax() - mLayerExtent.xMin() ) / 
mRasterXDim;
     double yres = ( mLayerExtent.yMax() - mLayerExtent.yMin() ) / 
mRasterYDim;

     // Offset, not the cell index -> flor
     int col = ( int ) floor(( x - mLayerExtent.xMin() ) / xres );
     int row = ( int ) floor(( mLayerExtent.yMax() - y ) / yres );

     QgsDebugMsg( "row = " + QString::number( row ) + " col = " + 
QString::number( col ) );

     GDALRasterBandH gdalBand = GDALGetRasterBand( mGdalDataset, iband );
     GDALDataType type = GDALGetRasterDataType( gdalBand );
     int size = GDALGetDataTypeSize( type ) / 8;
     void *data = CPLMalloc( size );

     CPLErr err = GDALRasterIO( gdalBand, GF_Read, col, row, 1, 1,
			       data, 1, 1, type, 0, 0 );

     if ( err != CPLE_None )
       {
         QgsLogger::warning( "RaterIO error: " + QString( 
CPLGetLastErrorMsg() ) );
       }

     double value = readValue( data, type, 0 );
#ifdef QGISDEBUG
     QgsLogger::debug( "value", value, 1, __FILE__, __FUNCTION__, 
__LINE__ );
#endif
     if ( mValidNoDataValue && ( mNoDataValue == value || value != value ) )
       {
	return nan();
       }
     free( data );
     return value;
   }



} // void QgsRasterLayer::getValue







Ghislain Picard wrote:
> 
> Hi,
> 
> I'd like to access large part of the raster dataset to build scattergram 
> from two raster layers.
> 
> Is there any mean to access raster data from python plugins more 
> efficiently than with the QgsRasterLayer::identify method ?
> 
> A version of identify that return a single band would be interesting too.
> 
> Finally what do you thing of a public access to mGDALDataset ?
> 
> Cheers,
> Ghislain
> P.S. new on the list. Aims at enhancing qgis for the analysis of remote 
> sensing images.
> 
>