EVOLUTION-MANAGER

Edit File: postgisrasterrasterband.cpp

/***********************************************************************
 * File :    postgisrasterrasterband.cpp
 * Project:  PostGIS Raster driver
 * Purpose:  GDAL RasterBand implementation for PostGIS Raster driver
 * Author:   Jorge Arevalo, jorge.arevalo@deimos-space.com
 *                          jorgearevalo@libregis.org
 * 
 * Author:	 David Zwarg, dzwarg@azavea.com
 *
 * Last changes: $Id: $
 *
 ***********************************************************************
 * Copyright (c) 2009 - 2013, Jorge Arevalo, David Zwarg
 * Copyright (c) 2013, Even Rouault
 *
 * Permission is hereby granted, free of charge, to any person obtaining 
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be 
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 
 * NONINFRINGEMENT.IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 
 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN 
 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN 
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 
 * SOFTWARE.
 **********************************************************************/
#include "postgisraster.h"

/**
 * \brief Constructor.
 * 
 * nBand it's just necessary for overview band creation
 */
PostGISRasterRasterBand::PostGISRasterRasterBand(
    PostGISRasterDataset * poDS, int nBand,
    GDALDataType eDataType, GBool bNoDataValueSet, double dfNodata, 
    GBool bIsOffline = false) : 
    VRTSourcedRasterBand(poDS, nBand)
{
    /* Basic properties */
    this->poDS = poDS;
    this->bIsOffline = bIsOffline;
    this->nBand = nBand;
    
    this->eDataType = eDataType;
    this->bNoDataValueSet = bNoDataValueSet;
    this->dfNoDataValue = dfNodata;

this->pszSchema = poDS->pszSchema;
    this->pszTable = poDS->pszTable;
    this->pszColumn = poDS->pszColumn;

nRasterXSize = poDS->GetRasterXSize();
    nRasterYSize = poDS->GetRasterYSize();

/*******************************************************************
     * Finally, set the block size. We apply the same logic than in VRT 
     * driver.
     * 
     * We limit the size of a block with MAX_BLOCK_SIZE here to prevent
     * arrangements of just one big tile.
     * 
     * This value is just used in case whe only have 1 tile in the
     * table. Otherwise, the reading operations are performed by the
     * sources, not the PostGISRasterBand object itself.
     ******************************************************************/
    this->nBlockXSize = MIN(MAX_BLOCK_SIZE, this->nRasterXSize);
    this->nBlockYSize = MIN(MAX_BLOCK_SIZE, this->nRasterYSize);

#ifdef DEBUG_VERBOSE
    CPLDebug("PostGIS_Raster", 
        "PostGISRasterRasterBand constructor: Band created (srid = %d)", 
        poDS->nSrid);
    
    CPLDebug("PostGIS_Raster", 
        "PostGISRasterRasterBand constructor: Band size: (%d X %d)", 
        nRasterXSize, nRasterYSize);
        
    CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::Constructor: "
        "Block size (%dx%d)", this->nBlockXSize, this->nBlockYSize);
#endif
}

/***********************************************
 * \brief: Band destructor
 ***********************************************/
PostGISRasterRasterBand::~PostGISRasterRasterBand()
{
}
 
#ifdef notdef
/***********************************************************************
 * \brief Set the block data to the null value if it is set, or zero if 
 * there is no null data value.
 * Parameters:
 *  - void *: the block data
 * Returns: nothing
 **********************************************************************/
void PostGISRasterRasterBand::NullBlock(void *pData) 
{
    int nWords = nBlockXSize * nBlockYSize;
    int nDataTypeSize = GDALGetDataTypeSize(eDataType) / 8;
    
    int bNoDataSet;
    double dfNoData = GetNoDataValue(&bNoDataSet);
    if (!bNoDataSet) {
        memset(pData, 0, nWords * nDataTypeSize);
    } 
    
    else {
        GDALCopyWords(&dfNoData, GDT_Float64, 0,
                      pData, eDataType, nDataTypeSize,
                      nWords);
    }
}

/***********************************************************************
 * \brief Returns the metadata for this band
 * 
 * If the metadata is actually stored in band's properties, simply 
 * returns them. Otherwise, it raises a query to fetch metadata
 * FROM db
 **********************************************************************/ 
GBool PostGISRasterRasterBand::GetBandMetadata(
    GDALDataType * peDataType, GBool * pbHasNoData, double * pdfNoData)
{
    // No need to raise a query
    if (eDataType != GDT_Unknown) {
        if (peDataType)
            *peDataType = eDataType;
            
        if (pbHasNoData)
            *pbHasNoData = bNoDataValueSet;
            
        if (pdfNoData)
            *pdfNoData = dfNoDataValue;
            
        return true;
    }
    
    /**
     * Queries are expensive. So, we only raise them if all parameters 
     * are not null
     **/
    if (!peDataType || !pbHasNoData || !pdfNoData) {
        return false;
    }
    
    /**
     * It is safe to assume all the tiles will have the same values for
     * metadata properties. That was checked during band's construction
     * (or we simply trusted the user, to avoid expensive checkings). 
     * So, we can limit the results to just one.
     **/
    int nTuples = 0;
    CPLString osCommand = NULL;
    PGresult * poResult = NULL;
    PostGISRasterDataset * poRDS = (PostGISRasterDataset *)poDS;
    
    osCommand.Printf("st_bandpixeltype(%s, %d), "
        "st_bandnodatavalue(%s, %d) is not null, "
        "st_bandnodatavalue(%s, %d) FROM %s.%s limit 1", pszColumn, 
        nBand, pszColumn, nBand, pszColumn, nBand, pszSchema, pszTable);
        
#ifdef DEBUG_QUERY
    CPLDebug("PostGIS_Raster", 
        "PostGISRasterRasterBand::GetBandMetadata(): Query: %s", 
        osCommand.c_str());
#endif
    
    poResult = PQexec(poRDS->poConn, osCommand.c_str());
    nTuples = PQntuples(poResult);
    
    /* Error getting info FROM database */
    if (poResult == NULL || 
        PQresultStatus(poResult) != PGRES_TUPLES_OK ||
        nTuples <= 0) {
        
        ReportError(CE_Failure, CPLE_AppDefined, 
            "Error getting band metadata while creating raster "
            "bands");
            
        CPLDebug("PostGIS_Raster", 
            "PostGISRasterDataset::GetBandMetadata(): %s", 
            PQerrorMessage(poRDS->poConn));
        
        if (poResult)
            PQclear(poResult);

return false;
    }
    
    // Fill band metadata values
    GBool bSignedByte = false;
    int nBitsDepth = 8;
    char* pszDataType = NULL;
    
    pszDataType = CPLStrdup(PQgetvalue(poResult, 0, 0));
    
    TranslateDataType(pszDataType, &eDataType, &nBitsDepth, 
            &bSignedByte);
            
    // Add pixeltype to image structure domain
    if (bSignedByte) {
        SetMetadataItem("PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
    }

// Add NBITS to metadata only for sub-byte types
    if (nBitsDepth < 8)
        SetMetadataItem("NBITS", CPLString().Printf( "%d", nBitsDepth ),
            "IMAGE_STRUCTURE" );
            
    CPLFree(pszDataType);
            
    bNoDataValueSet = 
            EQUALN(PQgetvalue(poResult, 0, 1), "t", sizeof(char));
            
    dfNoDataValue = CPLAtof(PQgetvalue(poResult, 0, 2));
    
    // Fill output arguments
    *peDataType = eDataType;
    *pbHasNoData = bNoDataValueSet;
    *pdfNoData = dfNoDataValue;
    
    return true;
}
#endif

/********************************************************
 * \brief Set nodata value to a buffer
 ********************************************************/
void PostGISRasterRasterBand::NullBuffer(void* pData,
                                         int nBufXSize,
                                         int nBufYSize,
                                         GDALDataType eBufType, 
                                         int nPixelSpace,
                                         int nLineSpace)
{
    int j;
    for(j = 0; j < nBufYSize; j++)
    {
        double dfVal = 0.0;
        if( bNoDataValueSet )
            dfVal = dfNoDataValue;
        GDALCopyWords(&dfVal, GDT_Float64, 0,
                    (GByte*)pData + j * nLineSpace, eBufType, nPixelSpace,
                    nBufXSize);
    }
}

/********************************************************
 * \brief SortTilesByPKID
 ********************************************************/
static int SortTilesByPKID(const void* a, const void* b)
{
    PostGISRasterTileDataset* pa = *(PostGISRasterTileDataset** )a;
    PostGISRasterTileDataset* pb = *(PostGISRasterTileDataset** )b;
    return strcmp(pa->GetPKID(), pb->GetPKID());
}

/**
 * Read/write a region of image data for this band.
 *
 * This method allows reading a region of a PostGISRasterBand into a buffer. 
 * The write support is still under development
 *
 * The function fetches all the raster data that intersects with the region
 * provided, and store the data in the GDAL cache.
 *
 * It automatically takes care of data type translation if the data type
 * (eBufType) of the buffer is different than that of the PostGISRasterRasterBand.
 *
 * The nPixelSpace and nLineSpace parameters allow reading into FROM various 
 * organization of buffers.
 *
 * @param eRWFlag Either GF_Read to read a region of data (GF_Write, to write
 * a region of data, yet not supported)
 *
 * @param nXOff The pixel offset to the top left corner of the region of the
 * band to be accessed. This would be zero to start FROM the left side.
 *
 * @param nYOff The line offset to the top left corner of the region of the band
 * to be accessed. This would be zero to start FROM the top.
 *
 * @param nXSize The width of the region of the band to be accessed in pixels.
 *
 * @param nYSize The height of the region of the band to be accessed in lines.
 *
 * @param pData The buffer into which the data should be read, or FROM which it
 * should be written. This buffer must contain at least
 * nBufXSize * nBufYSize * nBandCount words of type eBufType. It is organized in
 * left to right,top to bottom pixel order. Spacing is controlled by the
 * nPixelSpace, and nLineSpace parameters.
 *
 * @param nBufXSize the width of the buffer image into which the desired region
 * is to be read, or FROM which it is to be written.
 *
 * @param nBufYSize the height of the buffer image into which the desired region
 * is to be read, or FROM which it is to be written.
 *
 * @param eBufType the type of the pixel values in the pData data buffer. The
 * pixel values will automatically be translated to/FROM the
 * PostGISRasterRasterBand data type as needed.
 *
 * @param nPixelSpace The byte offset FROM the start of one pixel value in pData
 * to the start of the next pixel value within a scanline. If defaulted (0) the
 * size of the datatype eBufType is used.
 *
 * @param nLineSpace The byte offset FROM the start of one scanline in pData to
 * the start of the next. If defaulted (0) the size of the datatype
 * eBufType * nBufXSize is used.
 *
 * @return CE_Failure if the access fails, otherwise CE_None.
 */

CPLErr PostGISRasterRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, 
    int nYOff, int nXSize, int nYSize, void * pData, int nBufXSize, 
    int nBufYSize, GDALDataType eBufType, int nPixelSpace, 
    int nLineSpace)
{
    /**
     * TODO: Write support not implemented yet
     **/
    if (eRWFlag == GF_Write) {
        ReportError(CE_Failure, CPLE_NotSupported,
            "Writing through PostGIS Raster band not supported yet");
        
        return CE_Failure;
    }

/*******************************************************************
     * Do we have overviews that would be appropriate to satisfy this
     * request?                                                   
     ******************************************************************/
    if( (nBufXSize < nXSize || nBufYSize < nYSize) && 
        GetOverviewCount() > 0 )
    {
        if(OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, 
            pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, 
            nLineSpace) == CE_None)
                
        return CE_None;
    }

PostGISRasterDataset * poRDS = (PostGISRasterDataset *)poDS;

int bSameWindowAsOtherBand =
        (nXOff == poRDS->nXOffPrev &&
         nYOff == poRDS->nYOffPrev &&
         nXSize == poRDS->nXSizePrev &&
         nYSize == poRDS->nYSizePrev);
    poRDS->nXOffPrev = nXOff;
    poRDS->nYOffPrev = nYOff;
    poRDS->nXSizePrev = nXSize;
    poRDS->nYSizePrev = nYSize;

/* Logic to determine if bands are read in order 1, 2, ... N */
    /* If so, then use multi-band caching, otherwise do just single band caching */
    if( poRDS->bAssumeMultiBandReadPattern )
    {
        if( nBand != poRDS->nNextExpectedBand )
        {
            CPLDebug("PostGIS_Raster",
                    "Disabling multi-band caching since band access pattern does not match");
            poRDS->bAssumeMultiBandReadPattern = false;
            poRDS->nNextExpectedBand = 1;
        }
        else
        {
            poRDS->nNextExpectedBand ++;
            if( poRDS->nNextExpectedBand > poRDS->GetRasterCount() )
                poRDS->nNextExpectedBand = 1;
        }
    }
    else
    {
        if( nBand == poRDS->nNextExpectedBand )
        {
            poRDS->nNextExpectedBand ++;
            if( poRDS->nNextExpectedBand > poRDS->GetRasterCount() )
            {
                CPLDebug("PostGIS_Raster", "Re-enabling multi-band caching");
                poRDS->bAssumeMultiBandReadPattern = true;
                poRDS->nNextExpectedBand = 1;
            }
        }
    }
    
#ifdef DEBUG_VERBOSE
    CPLDebug("PostGIS_Raster", 
            "PostGISRasterRasterBand::IRasterIO: "
            "nBand = %d, nXOff = %d, nYOff = %d, nXSize = %d, nYSize = %d, nBufXSize = %d, nBufYSize = %d",
             nBand, nXOff, nYOff, nXSize,  nYSize, nBufXSize, nBufYSize);
#endif

#ifdef notdef
    /*******************************************************************
     * Optimization: We just have one tile. So, we can read it with
     * IReadBlock
     * 
     * TODO: Review it. It's not working (see comment in 
     * PostGISRasterDataset::ConstructOneDatasetFromTiles)
     ******************************************************************/
    
    if (poRDS->nTiles == 1) {
        
        return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, 
            nYSize, pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, 
            nLineSpace);
    }
#endif    
    
    /*******************************************************************
     * Several tiles: we first look in all our sources caches. Missing
     * blocks are queried
     ******************************************************************/
    double adfProjWin[8];
    int nFeatureCount = 0;
    CPLRectObj sAoi;

poRDS->PolygonFromCoords(nXOff, nYOff, nXOff + nXSize, nYOff + nYSize, adfProjWin); 
    // (p[6], p[7]) is the minimum (x, y), and (p[2], p[3]) the max
    sAoi.minx = adfProjWin[6];
    sAoi.maxx = adfProjWin[2];
    if( adfProjWin[7] < adfProjWin[3] )
    {
        sAoi.miny = adfProjWin[7];
        sAoi.maxy = adfProjWin[3];
    }
    else
    {
        sAoi.maxy = adfProjWin[7];
        sAoi.miny = adfProjWin[3];
    }

#ifdef DEBUG_VERBOSE
    CPLDebug("PostGIS_Raster", 
            "PostGISRasterRasterBand::IRasterIO: "
            "Intersection box: (%f, %f) - (%f, %f)", sAoi.minx, 
            sAoi.miny, sAoi.maxx, sAoi.maxy);
#endif

if (poRDS->hQuadTree == NULL)
    {
        ReportError(CE_Failure, CPLE_AppDefined,
            "Could not read metadata index.");
        return CE_Failure; 
    }

NullBuffer(pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace);

if( poRDS->bBuildQuadTreeDynamically && !bSameWindowAsOtherBand )
    {
        if( !(poRDS->LoadSources(nXOff, nYOff, nXSize, nYSize, nBand)) )
            return CE_Failure;
    }

// Matching sources, to avoid a dumb for loop over the sources
    PostGISRasterTileDataset ** papsMatchingTiles = 
        (PostGISRasterTileDataset **)
            CPLQuadTreeSearch(poRDS->hQuadTree, &sAoi, &nFeatureCount);
         
    // No blocks found. This is not an error (the raster may have holes)
    if (nFeatureCount == 0) {
        CPLFree(papsMatchingTiles);
 
        return CE_None; 
    }

int i;

/**
     * We need to store the max, min coords for the missing tiles in
     * any place. This is as good as any other
     **/ 
    sAoi.minx = 0.0;
    sAoi.miny = 0.0;
    sAoi.maxx = 0.0;
    sAoi.maxy = 0.0;
    
    GIntBig nMemoryRequiredForTiles = 0;
    CPLString osIDsToFetch;
    int nTilesToFetch = 0;
    int nBandDataTypeSize = GDALGetDataTypeSize(eDataType) / 8;

// Loop just over the intersecting sources
    for(i = 0; i < nFeatureCount; i++) {
        PostGISRasterTileDataset *poTile = papsMatchingTiles[i];
        PostGISRasterTileRasterBand* poTileBand = 
            (PostGISRasterTileRasterBand *)poTile->GetRasterBand(nBand);

nMemoryRequiredForTiles += poTileBand->GetXSize() * poTileBand->GetYSize() *
            nBandDataTypeSize;
        
        // Missing tile: we'll need to query for it
        if (!poTileBand->IsCached()) {
            
            // If we have a PKID, add the tile PKID to the list
            if (poTile->pszPKID != NULL)
            {
                if( osIDsToFetch.size() != 0 )
                    osIDsToFetch += ",";
                osIDsToFetch += "'";
                osIDsToFetch += poTile->pszPKID;
                osIDsToFetch += "'";
            }

double dfTileMinX, dfTileMinY, dfTileMaxX, dfTileMaxY;
            poTile->GetExtent(&dfTileMinX, &dfTileMinY,
                              &dfTileMaxX, &dfTileMaxY);

/**
             * We keep the general max and min values of all the missing
             * tiles, to raise a query that intersect just that area.
             * 
             * TODO: In case of just a few tiles and very separated,
             * this strategy is clearly suboptimal. We'll get our
             * missing tiles, but with a lot of other not needed tiles.
             * 
             * A possible optimization will be to simply rely on the
             * I/O method of the source (must be implemented), in case
             * we have minus than a reasonable amount of tiles missing.
             * Another criteria to decide would be how separated the
             * tiles are. Two queries for just two adjacent tiles is 
             * also a dumb strategy.
             **/
            if( nTilesToFetch == 0 )
            {
                sAoi.minx = dfTileMinX;
                sAoi.miny = dfTileMinY;
                sAoi.maxx = dfTileMaxX;
                sAoi.maxy = dfTileMaxY;
            }
            else
            {
                if (dfTileMinX < sAoi.minx)
                    sAoi.minx = dfTileMinX;
                    
                if (dfTileMinY < sAoi.miny)
                    sAoi.miny = dfTileMinY;
                    
                if (dfTileMaxX > sAoi.maxx)
                    sAoi.maxx = dfTileMaxX;
                    
                if (dfTileMaxY > sAoi.maxy)
                    sAoi.maxy = dfTileMaxY;
            }
            
            nTilesToFetch ++;
            
        }
    }

/* Determine caching strategy */
    int bAllBandCaching = FALSE;
    if (nTilesToFetch > 0)
    {
        GIntBig nCacheMax = (GIntBig) GDALGetCacheMax64();
        if( nMemoryRequiredForTiles > nCacheMax )
        {
            CPLDebug("PostGIS_Raster",
                    "For best performance, the block cache should be able to store " CPL_FRMT_GIB
                    " bytes for the tiles of the requested window, "
                    "but it is only " CPL_FRMT_GIB " byte large",
                    nMemoryRequiredForTiles, nCacheMax );
            nTilesToFetch = 0;
        }
        
        if( poRDS->GetRasterCount() > 1 && poRDS->bAssumeMultiBandReadPattern )
        {
            GIntBig nMemoryRequiredForTilesAllBands =
                nMemoryRequiredForTiles * poRDS->GetRasterCount();
            if( nMemoryRequiredForTilesAllBands <= nCacheMax )
            {
                bAllBandCaching = TRUE;
            }
            else
            {
                CPLDebug("PostGIS_Raster", "Caching only this band, but not all bands. "
                         "Cache should be " CPL_FRMT_GIB " byte large for that",
                         nMemoryRequiredForTilesAllBands);
            }
        }
    }
    
    // Raise a query for missing tiles and cache them
    if (nTilesToFetch > 0) {
        
        /**
         * There are several options here, to raise the query.
         * - Get all the tiles which PKID is in a list of missing 
         *   PKIDs.
         * - Get all the tiles that intersect a polygon constructed
         *   based on the (min - max) values calculated before.
         * - Get all the tiles with upper left pixel included in the
         *   range (min - max) calculated before.
         * 
         * The first option is the most efficient one when a PKID exists.
         * After that, the second one is the most efficient one when a
         * spatial index exists.
         * The third one is the only one available when neither a PKID or spatial
         * index exist.
         **/ 
        CPLString osCommand;
        PGresult * poResult;

CPLString osRasterToFetch;
        if (bAllBandCaching)
            osRasterToFetch = pszColumn;
        else
            osRasterToFetch.Printf("ST_Band(%s, %d)", pszColumn, nBand);

int bHasWhere = FALSE;
        if (osIDsToFetch.size() && (poRDS->bIsFastPK || !(poRDS->HasSpatialIndex())) ) {
            osCommand.Printf("SELECT %s, "
                "ST_Metadata(%s), %s FROM %s.%s",
                osRasterToFetch.c_str(), pszColumn,
                poRDS->GetPrimaryKeyRef(), pszSchema, pszTable);
            if( nTilesToFetch < poRDS->nTiles || poRDS->bBuildQuadTreeDynamically )
            {
                bHasWhere = TRUE;
                osCommand += " WHERE ";
                osCommand += poRDS->pszPrimaryKeyName;
                osCommand += " IN (";
                osCommand += osIDsToFetch;
                osCommand += ")";
            }
        } 
        
        else {
            CPLLocaleC oCLocale; // Force C locale to avoid commas instead of decimal points (for QGIS e.g.)
            bHasWhere = TRUE;
            osCommand.Printf("SELECT %s, ST_Metadata(%s), %s FROM %s.%s WHERE ",
                             osRasterToFetch.c_str(), pszColumn,
                             (poRDS->GetPrimaryKeyRef()) ? poRDS->GetPrimaryKeyRef() : "'foo'",
                             pszSchema, pszTable);
            if( poRDS->HasSpatialIndex() )
            {
                osCommand += CPLSPrintf("%s && "
                        "ST_GeomFromText('POLYGON((%.18f %.18f,%.18f %.18f,%.18f %.18f,%.18f %.18f,%.18f %.18f))')",
                        pszColumn,
                        adfProjWin[0], adfProjWin[1],
                        adfProjWin[2], adfProjWin[3],
                        adfProjWin[4], adfProjWin[5],
                        adfProjWin[6], adfProjWin[7],
                        adfProjWin[0], adfProjWin[1]);
            }
            else
            {
                #define EPS 1e-5
                osCommand += CPLSPrintf("ST_UpperLeftX(%s)"
                    " BETWEEN %f AND %f AND ST_UpperLeftY(%s) BETWEEN "
                    "%f AND %f", pszColumn, sAoi.minx-EPS, sAoi.maxx+EPS, 
                    pszColumn, sAoi.miny-EPS, sAoi.maxy+EPS);
            }
        }
    
        if( poRDS->pszWhere != NULL )
        {
            if( bHasWhere )
                osCommand += " AND (";
            else
                osCommand += " WHERE (";
            osCommand += poRDS->pszWhere;
            osCommand += ")";
        }
        
        poResult = PQexec(poRDS->poConn, osCommand.c_str());
         
#ifdef DEBUG_QUERY
        CPLDebug("PostGIS_Raster", 
            "PostGISRasterRasterBand::IRasterIO(): Query = \"%s\" --> number of rows = %d", 
            osCommand.c_str(), poResult ? PQntuples(poResult) : 0 );
#endif
        
        if (poResult == NULL || 
            PQresultStatus(poResult) != PGRES_TUPLES_OK || 
            PQntuples(poResult) < 0) {
            
            if (poResult)
                PQclear(poResult);
     
            CPLError(CE_Failure, CPLE_AppDefined, 
                "PostGISRasterRasterBand::IRasterIO(): %s", 
                PQerrorMessage(poRDS->poConn));
            
            // Free the object that holds pointers to matching tiles
            CPLFree(papsMatchingTiles);
            return CE_Failure;  
        }

/**
         * No data. Return the buffer filled with nodata values
         **/
        else if (PQntuples(poResult) == 0) {
            PQclear(poResult);

// Free the object that holds pointers to matching tiles
            CPLFree(papsMatchingTiles);
            return CE_None; 
        }
    
        /**
         * Ok, we loop over the results
         **/
        int nTuples = PQntuples(poResult);
        for(i = 0; i < nTuples; i++)
        {
            const char* pszMetadata = PQgetvalue(poResult, i, 1);
            const char* pszRaster = PQgetvalue(poResult, i, 0);
            const char *pszPKID = (poRDS->GetPrimaryKeyRef() != NULL) ?  PQgetvalue(poResult, i, 2) : NULL;
            poRDS->CacheTile(pszMetadata, pszRaster, pszPKID, nBand, bAllBandCaching);
        } // All tiles have been added to cache

PQclear(poResult);
    } // End missing tiles

/* -------------------------------------------------------------------- */
/*      Overlay each source in turn over top this.                      */
/* -------------------------------------------------------------------- */

CPLErr eErr = CE_None;
    /* Sort tiles by ascending PKID, so that the draw order is determinist */
    if( poRDS->GetPrimaryKeyRef() != NULL )
    {
        qsort(papsMatchingTiles, nFeatureCount, sizeof(PostGISRasterTileDataset*),
              SortTilesByPKID);
    }

for(i = 0; i < nFeatureCount && eErr == CE_None; i++)
    {
        PostGISRasterTileDataset *poTile = papsMatchingTiles[i];
        PostGISRasterTileRasterBand* poTileBand = 
            (PostGISRasterTileRasterBand *)poTile->GetRasterBand(nBand);
        eErr = 
            poTileBand->poSource->RasterIO( nXOff, nYOff, nXSize, nYSize, 
                                            pData, nBufXSize, nBufYSize, 
                                            eBufType, nPixelSpace, nLineSpace);
    }
    
    // Free the object that holds pointers to matching tiles
    CPLFree(papsMatchingTiles);

return eErr;
}

/**
 * \brief Set the no data value for this band.
 * Parameters:
 *  - double: The nodata value
 * Returns:
 *  - CE_None.
 */
CPLErr PostGISRasterRasterBand::SetNoDataValue(double dfNewValue) {
    dfNoDataValue = dfNewValue;

return CE_None;
}

/**
 * \brief Fetch the no data value for this band.
 * Parameters:
 *  - int *: pointer to a boolean to use to indicate if a value is actually
 *          associated with this layer. May be NULL (default).
 *  Returns:
 *  - double: the nodata value for this band.
 */
double PostGISRasterRasterBand::GetNoDataValue(int *pbSuccess) {
    if (pbSuccess != NULL)
        *pbSuccess = (int) bNoDataValueSet;

return dfNoDataValue;
}

/***************************************************
 * \brief Return the number of overview layers available
 ***************************************************/
int PostGISRasterRasterBand::GetOverviewCount()
{
    PostGISRasterDataset * poRDS = (PostGISRasterDataset *)poDS;
    return poRDS->GetOverviewCount();
}

/**********************************************************
 * \brief Fetch overview raster band object
 **********************************************************/
GDALRasterBand * PostGISRasterRasterBand::GetOverview(int i)
{
    if (i < 0 || i >= GetOverviewCount())
        return NULL;
    
    PostGISRasterDataset* poRDS = (PostGISRasterDataset *)poDS;
    PostGISRasterDataset* poOverviewDS = poRDS->GetOverviewDS(i);
    if( poOverviewDS->nBands == 0 )
    {
        if (!poOverviewDS->SetRasterProperties(NULL) ||
             poOverviewDS->GetRasterCount() != poRDS->GetRasterCount())
        {
            CPLDebug("PostGIS_Raster",
                     "Request for overview %d of band %d failed", i, nBand);
            return NULL;
        }
    }

return poOverviewDS->GetRasterBand(nBand);
}

#ifdef notdef
/*****************************************************
 * \brief Read a natural block of raster band data
 *****************************************************/
CPLErr PostGISRasterRasterBand::IReadBlock(int nBlockXOff, 
        int nBlockYOff, void * pImage)
{
    PGresult * poResult = NULL;
    CPLString osCommand;
    int nXOff = 0;
    int nYOff = 0;
    int nNaturalBlockXSize = 0;
    int nNaturalBlockYSize = 0;
    double adfProjWin[8];
    PostGISRasterDataset * poRDS = (PostGISRasterDataset *)poDS;
    
    int nPixelSize = GDALGetDataTypeSize(eDataType) / 8;
    
    // Construct a polygon to intersect with
    GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
    
    nXOff = nBlockXOff * nNaturalBlockXSize;
    nYOff = nBlockYOff * nNaturalBlockYSize;
    
    poRDS->PolygonFromCoords(nXOff, nYOff, nXOff + nNaturalBlockXSize, nYOff + nNaturalBlockYSize, adfProjWin);
    
    // Raise the query
    if (poRDS->pszWhere == NULL) {
        osCommand.Printf("SELECT st_band(%s, %d) FROM %s.%s "
            "WHERE st_intersects(%s, ST_PolygonFromText"
            "('POLYGON((%.17f %.17f, %.17f %.17f, %.17f %.17f, %.17f "
            "%.17f, %.17f %.17f))', %d))", pszColumn, nBand, pszSchema, 
            pszTable, pszColumn, adfProjWin[0], adfProjWin[1], 
            adfProjWin[2], adfProjWin[3],  adfProjWin[4], adfProjWin[5], 
            adfProjWin[6], adfProjWin[7], adfProjWin[0], adfProjWin[1], 
            poRDS->nSrid);
    }

else {
        osCommand.Printf("SELECT st_band(%s, %d) FROM %s.%s WHERE (%s) "
            "AND st_intersects(%s, ST_PolygonFromText"
            "('POLYGON((%.17f %.17f, %.17f %.17f, %.17f %.17f, %.17f "
            "%.17f, %.17f %.17f))', %d))", pszColumn, nBand, pszSchema, 
            pszTable, poRDS->pszWhere, pszColumn, adfProjWin[0], 
            adfProjWin[1], adfProjWin[2], adfProjWin[3], adfProjWin[4], 
            adfProjWin[5], adfProjWin[6], adfProjWin[7], adfProjWin[0], 
            adfProjWin[1], poRDS->nSrid);
    }

#ifdef DEBUG_QUERY
    CPLDebug("PostGIS_Raster", 
        "PostGISRasterRasterBand::IReadBlock(): Query = %s", 
            osCommand.c_str());
#endif

poResult = PQexec(poRDS->poConn, osCommand.c_str());
    if (poResult == NULL || 
        PQresultStatus(poResult) != PGRES_TUPLES_OK || 
        PQntuples(poResult) < 0) {
        
        if (poResult)
            PQclear(poResult);
 
        ReportError(CE_Failure, CPLE_AppDefined, 
            "Error retrieving raster data FROM database");

CPLDebug("PostGIS_Raster", 
            "PostGISRasterRasterBand::IRasterIO(): %s", 
                PQerrorMessage(poRDS->poConn));
        
        return CE_Failure;  
    }

/**
     * No data. Return the buffer filled with nodata values
     **/
    else if (PQntuples(poResult) == 0) {
        PQclear(poResult);
        
        CPLDebug("PostGIS_Raster", 
            "PostGISRasterRasterBand::IRasterIO(): Null block");

NullBlock(pImage);

return CE_None; 
    }
    
    
    /**
     * Ok, we get the data. Only data size, without payload 
     * 
     * TODO: Check byte order
     **/
    int nExpectedDataSize = nNaturalBlockXSize * nNaturalBlockYSize *
        nPixelSize;
        
    int nWKBLength = 0;
    
    GByte * pbyData = CPLHexToBinary(PQgetvalue(poResult, 0, 0), 
        &nWKBLength);
        
    char * pbyDataToRead = (char*)GET_BAND_DATA(pbyData,nBand, 
            nPixelSize, nExpectedDataSize);

memcpy(pImage, pbyDataToRead, nExpectedDataSize * sizeof(char));
    
    CPLDebug("PostGIS_Raster", "IReadBlock: Copied %d bytes FROM block "
            "(%d, %d) to %p", nExpectedDataSize, nBlockXOff, 
            nBlockYOff, pImage);

CPLFree(pbyData);
    PQclear(poResult);

return CE_None;
}
#endif

/**
 * \brief How should this band be interpreted as color?
 * GCI_Undefined is returned when the format doesn't know anything about the 
 * color interpretation. 
 **/
GDALColorInterp PostGISRasterRasterBand::GetColorInterpretation()
{
    if (poDS->GetRasterCount() == 1) {
        eColorInterp = GCI_GrayIndex;
    }

else if (poDS->GetRasterCount() == 3) {
        if (nBand == 1)
            eColorInterp = GCI_RedBand;
        else if( nBand == 2 )
            eColorInterp = GCI_GreenBand;
        else if( nBand == 3 )
            eColorInterp = GCI_BlueBand;
        else
            eColorInterp = GCI_Undefined;
    }

else {
        eColorInterp = GCI_Undefined;
    }
    
    return eColorInterp; 
}

/************************************************************************/
/*                             GetMinimum()                             */
/************************************************************************/

double PostGISRasterRasterBand::GetMinimum( int *pbSuccess )
{
    PostGISRasterDataset * poRDS = (PostGISRasterDataset *)poDS;
    if( poRDS->bBuildQuadTreeDynamically && poRDS->nTiles == 0 )
    {
        if( pbSuccess )
            *pbSuccess = FALSE;
        return 0.0;
    }
    return VRTSourcedRasterBand::GetMaximum(pbSuccess);
}

/************************************************************************/
/*                             GetMaximum()                             */
/************************************************************************/

double PostGISRasterRasterBand::GetMaximum( int *pbSuccess )
{
    PostGISRasterDataset * poRDS = (PostGISRasterDataset *)poDS;
    if( poRDS->bBuildQuadTreeDynamically && poRDS->nTiles == 0 )
    {
        if( pbSuccess )
            *pbSuccess = FALSE;
        return 0.0;
    }
    return VRTSourcedRasterBand::GetMaximum(pbSuccess);
}

/************************************************************************/
/*                       ComputeRasterMinMax()                          */
/************************************************************************/

CPLErr PostGISRasterRasterBand::ComputeRasterMinMax( int bApproxOK, double* adfMinMax )
{
    if( nRasterXSize < 1024 && nRasterYSize < 1024 )
        return VRTSourcedRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);

int nOverviewCount = GetOverviewCount();
    for(int i = 0; i < nOverviewCount; i++)
    {
        if( GetOverview(i)->GetXSize() < 1024 && GetOverview(i)->GetYSize() < 1024 )
            return GetOverview(i)->ComputeRasterMinMax(bApproxOK, adfMinMax);
    }
    
    return CE_Failure;
}