LTP GCOV extension - code coverage report
Current view: directory - frmts/wktraster - wktrasterrasterband.cpp
Test: gdal_filtered.info
Date: 2010-07-12 Instrumented lines: 292
Code covered: 0.0 % Executed lines: 0

       1                 : 
       2                 : #include <string>
       3                 : 
       4                 : 
       5                 : #include "cpl_string.h"
       6                 : 
       7                 : /******************************************************************************
       8                 :  * File :    wktrasterrasterband.cpp
       9                 :  * Project:  WKT Raster driver
      10                 :  * Purpose:  GDAL Dataset code for WKTRaster driver 
      11                 :  * Author:   Jorge Arevalo, jorgearevalo@gis4free.org
      12                 :  * 
      13                 :  * Last changes:
      14                 :  * $Id: wktrasterrasterband.cpp 19624 2010-05-07 10:33:25Z jorgearevalo $
      15                 :  *
      16                 :  * TODO:
      17                 :  *  - Block caching, to avoid fetching the whole raster from database each time
      18                 :  *    IReadBlock is called.
      19                 :  *  - Read outdb rasters. Take into account that the outdb raster may don't have
      20                 :  *    the same block structure...
      21                 :  *  - Update raster_columns table if values read from IReadBlock are different
      22                 :  *
      23                 :  ******************************************************************************
      24                 :  * Copyright (c) 2009, Jorge Arevalo, jorgearevalo@gis4free.org
      25                 :  *
      26                 :  * Permission is hereby granted, free of charge, to any person obtaining a
      27                 :  * copy of this software and associated documentation files (the "Software"),
      28                 :  * to deal in the Software without restriction, including without limitation
      29                 :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      30                 :  * and/or sell copies of the Software, and to permit persons to whom the
      31                 :  * Software is furnished to do so, subject to the following conditions:
      32                 :  *
      33                 :  * The above copyright notice and this permission notice shall be included
      34                 :  * in all copies or substantial portions of the Software.
      35                 :  *
      36                 :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      37                 :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      38                 :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      39                 :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      40                 :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      41                 :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      42                 :  * DEALINGS IN THE SOFTWARE.
      43                 :  ******************************************************************************/
      44                 : #include "wktraster.h"
      45                 : #include "ogr_api.h"
      46                 : #include "ogr_geometry.h"
      47                 : #include <gdal_priv.h>
      48                 : 
      49                 : CPL_CVSID("$Id: wktrasterrasterband.cpp 19624 2010-05-07 10:33:25Z jorgearevalo $");
      50                 : 
      51                 : /**
      52                 :  * Constructor.
      53                 :  * Parameters:
      54                 :  *  - WKTRasterDataset *: The Dataset this band belongs to
      55                 :  *  - int: the band number
      56                 :  *  - GDALDataType: The data type for this band
      57                 :  *  - double: The nodata value.  Could be any kind of data (GByte, GUInt16,
      58                 :  *          GInt32...) but the variable has the bigger type.
      59                 :  *  - GBool: if the data type is signed byte or not. If yes, the SIGNEDBYTE
      60                 :  *          metadata value is added to IMAGE_STRUCTURE domain
      61                 :  *  - int: The bit depth, to add NBITS as metadata value in IMAGE_STRUCTURE
      62                 :  *          domain.
      63                 :  */
      64                 : WKTRasterRasterBand::WKTRasterRasterBand(WKTRasterDataset *poDS,
      65                 :         int nBand, GDALDataType hDataType, double dfNodata, GBool bSignedByte,
      66               0 :         int nBitD) {
      67                 :     
      68               0 :     VALIDATE_POINTER0(poDS, "WKTRasterRasterBand");
      69                 : 
      70               0 :     poDS = poDS;
      71               0 :     nBand = nBand;
      72                 : 
      73               0 :     nRasterXSize = poDS->GetRasterXSize();
      74               0 :     nRasterYSize = poDS->GetRasterYSize();
      75                 : 
      76               0 :     nBlockXSize = ((WKTRasterDataset *) poDS)->nBlockSizeX;
      77               0 :     nBlockYSize = ((WKTRasterDataset *) poDS)->nBlockSizeY;
      78                 : 
      79                 :      // Check Irregular blocking
      80               0 :     if (nBlockXSize == 0 || nBlockYSize == 0) {
      81                 :         CPLError(CE_Warning, CPLE_NotSupported,
      82               0 :                 "This band has irregular blocking, but is not supported yet");
      83                 :     }
      84                 : 
      85                 : 
      86               0 :     eAccess = ((WKTRasterDataset *) poDS)->GetAccess();
      87                 : 
      88                 :     // Get no data value and pixel type from dataset too
      89               0 :     eDataType = hDataType;
      90               0 :     dfNoDataValue = dfNodata;
      91                 : 
      92               0 :     nBitDepth = nBitD;
      93                 : 
      94                 :     // Add pixeltype to image structure domain
      95               0 :     if (bSignedByte == TRUE) {
      96               0 :         SetMetadataItem("PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
      97               0 :         bIsSignedByte = bSignedByte;
      98                 :     }
      99                 : 
     100                 :     // Add NBITS to metadata only for sub-byte types
     101               0 :     if (nBitDepth < 8)
     102                 :         SetMetadataItem("NBITS", CPLString().Printf( "%d", nBitDepth ),
     103               0 :             "IMAGE_STRUCTURE" );
     104               0 : }
     105                 : 
     106                 : /**
     107                 :  * Says if the datatype for this band is signedbyte.
     108                 :  * Parameters: none
     109                 :  * Returns:
     110                 :  *  - TRUE if the datatype for this band is signedbyte, FALSE otherwise
     111                 :  */
     112               0 :  GBool WKTRasterRasterBand::IsSignedByteDataType()
     113                 :  {
     114               0 :      return bIsSignedByte;
     115                 : }
     116                 : 
     117                 : /**
     118                 :  * Get the bit depth for this raster band
     119                 :  * Parameters: none.
     120                 :  * Returns: the bit depth
     121                 :  */
     122               0 :  int WKTRasterRasterBand::GetNBitDepth()
     123                 :  {
     124               0 :      return nBitDepth;
     125                 :  }
     126                 : 
     127                 : /**
     128                 :  * Write a block of data to the raster band. First establish if a corresponding
     129                 :  * row to the block already exists with a SELECT. If so, update the raster 
     130                 :  * column contents. If it does not exist create a new row for the block.
     131                 :  * Inputs:
     132                 :  *  int: horizontal block offset
     133                 :  *  int: vertical block offset
     134                 :  *  void *: The buffer wit the data to be write
     135                 :  * Output:
     136                 :  *  CE_None on success, CE_Failure on error
     137                 :  */
     138                 : CPLErr WKTRasterRasterBand::IWriteBlock(int nBlockXOff,
     139               0 :         int nBlockYOff, void * pImage) {
     140                 : 
     141                 :     // Check parameters
     142               0 :     if (pImage == NULL || nBlockXOff < 0 || nBlockYOff < 0) {
     143                 :         CPLError(CE_Failure, CPLE_NotSupported,
     144               0 :                 "Unsupported block size or NULL buffer");
     145               0 :         return CE_Failure;
     146                 :     }
     147                 : 
     148               0 :     WKTRasterDataset * poWKTRasterDS = (WKTRasterDataset *) poDS;
     149               0 :     CPLString osCommand;
     150                 :     PGresult * hPGresult;
     151                 :     int nPixelSize;
     152                 :     int nPixelXInitPosition;
     153                 :     int nPixelYInitPosition;
     154                 :     int nPixelXEndPosition;
     155                 :     int nPixelYEndPosition;
     156                 :     double adfTransform[6];
     157               0 :     double fProjXInit = 0.0;
     158               0 :     double fProjYInit = 0.0;
     159               0 :     double fProjXEnd = 0.0;
     160               0 :     double fProjYEnd = 0.0;
     161               0 :     double fProjLowerLeftX = 0.0;
     162               0 :     double fProjLowerLeftY = 0.0;
     163               0 :     double fProjUpperRightX = 0.0;
     164               0 :     double fProjUpperRightY = 0.0;
     165               0 :     GByte byMachineEndianess = NDR;
     166               0 :     int nTuples = 0;
     167               0 :     char * pszHexWkb = NULL;
     168               0 :     WKTRasterWrapper * poWKTRasterWrapper = NULL;
     169               0 :     WKTRasterBandWrapper * poWKTRasterBandWrapper = NULL;
     170               0 :     int nRid = 0;
     171               0 :     int nCurrentBandPixelSize = 0;
     172                 : 
     173                 : 
     174                 :     // Check machine endianess
     175                 : #ifdef CPL_LSB
     176               0 :     byMachineEndianess = NDR;
     177                 : #else
     178                 :     byMachineEndianess = XDR;
     179                 : #endif
     180                 : 
     181                 :     /***********************************************************************
     182                 :      * Get pixel size (divide by 8 because GDALGetDataTypeSize returns the
     183                 :      * size in bits)
     184                 :      ***********************************************************************/
     185               0 :     nPixelSize = GDALGetDataTypeSize(eDataType) / 8;
     186                 : 
     187                 :     /***********************************************************************
     188                 :      * nBlockXOff and nBlockYOff are block offsets. So, we first have to
     189                 :      * transform them in pixel/line coordinates, taking into account the
     190                 :      * size of a block.
     191                 :      ***********************************************************************/
     192               0 :     nPixelXInitPosition = nBlockXSize * nBlockXOff;
     193               0 :     nPixelYInitPosition = nBlockYSize * nBlockYOff;
     194                 : 
     195                 :     // Now, get the end of the block
     196               0 :     nPixelXEndPosition = nPixelXInitPosition + nBlockXSize;
     197               0 :     nPixelYEndPosition = nPixelYInitPosition + nBlockYSize;
     198                 : 
     199                 : 
     200                 :     /**************************************************************************
     201                 :      * Transform pixel/line coordinates into coordinates of the raster
     202                 :      * reference system.
     203                 :      * NOTE: I take the georeference information from dataset. The SQL function
     204                 :      * ST_GdalGeoTransform takes the same information from the raster row, in
     205                 :      * array format.
     206                 :      **************************************************************************/
     207               0 :     poWKTRasterDS->GetGeoTransform(adfTransform);
     208                 :     fProjXInit = adfTransform[0] +
     209                 :             nPixelXInitPosition * adfTransform[1] +
     210               0 :             nPixelYInitPosition * adfTransform[2];
     211                 : 
     212                 :     fProjYInit = adfTransform[3] +
     213                 :             nPixelXInitPosition * adfTransform[4] +
     214               0 :             nPixelYInitPosition * adfTransform[5];
     215                 : 
     216                 :     fProjXEnd = adfTransform[0] +
     217                 :             nPixelXEndPosition * adfTransform[1] +
     218               0 :             nPixelYEndPosition * adfTransform[2];
     219                 : 
     220                 :     fProjYEnd = adfTransform[3] +
     221                 :             nPixelXEndPosition * adfTransform[4] +
     222               0 :             nPixelYEndPosition * adfTransform[5];
     223                 : 
     224                 : 
     225                 :     /*************************************************************************
     226                 :      * Now we have the block coordinates transformed into coordinates of the
     227                 :      * raster reference system. This coordinates are from:
     228                 :      *  - Upper left corner
     229                 :      *  - Lower right corner
     230                 :      * But for ST_MakeBox2D, we'll need block's coordinates of:
     231                 :      *  - Lower left corner
     232                 :      *  - Upper righ corner
     233                 :      *************************************************************************/
     234               0 :     fProjLowerLeftX = fProjXInit;
     235               0 :     fProjLowerLeftY = fProjYEnd;
     236                 : 
     237               0 :     fProjUpperRightX = fProjXEnd;
     238               0 :     fProjUpperRightY = fProjYInit;
     239                 : 
     240                 : 
     241                 :     /*************************************************************************
     242                 :      * Perform a spatial query that gives the tile/block (row of raster table)
     243                 :      * or tiles/blocks (in case of non-regular blocking) that should contain
     244                 :      * this block
     245                 :      *************************************************************************/
     246                 : 
     247               0 :     if (poWKTRasterDS->pszWhereClause != NULL) {
     248                 : 
     249                 :         /**
     250                 :          * Table has GIST index-
     251                 :          * We could disable sequential scanning, but for versions of PostGIS
     252                 :          * up to 0.8, this is no necessary. PostGIS >=0.8 is correctly
     253                 :          * integrated with query planner, thus PostgreSQL will use indexes
     254                 :          * whenever appropriate.
     255                 :          * NOTE: We should add version checking here, and disable sequential
     256                 :          * scanning when needed. See OGRDataSource::Open method.
     257                 :          */
     258               0 :         if (poWKTRasterDS->bTableHasGISTIndex) {
     259                 : 
     260                 :             /**
     261                 :              * This queries for the tiles that contains the given block. When
     262                 :              * regular_blocking rasters, the blocks will have the same size of
     263                 :              * a tile, so, we can use "contain" function and "~" operator.
     264                 :              * But when non-regular_blocking, the only way this will work is
     265                 :              * setting the block size to the smallest tile size. The problem is
     266                 :              * how do we know all the tiles size?
     267                 :              */
     268                 :             osCommand.Printf(
     269                 :                     "SELECT rid, %s FROM %s.%s WHERE %s ~ ST_SetSRID(ST_MakeBox2D\
     270                 :                     (ST_Point(%f, %f),ST_Point(%f, %f)), %d) AND %s",
     271                 :                     poWKTRasterDS->pszRasterColumnName,
     272                 :                     poWKTRasterDS->pszSchemaName,
     273                 :                     poWKTRasterDS->pszTableName,
     274                 :                     poWKTRasterDS->pszRasterColumnName, fProjLowerLeftX,
     275                 :                     fProjLowerLeftY, fProjUpperRightX, fProjUpperRightY,
     276               0 :                     poWKTRasterDS->nSrid, poWKTRasterDS->pszWhereClause);
     277                 : 
     278                 :         }            /**
     279                 :              * Table hasn't a GIST index. Normal searching
     280                 :              */
     281                 :         else {
     282                 :             osCommand.Printf(
     283                 :                     "SELECT rid, %s FROM %s.%s WHERE _ST_Contains(%s, ST_SetSRID(\
     284                 :                     ST_MakeBox2D(ST_Point(%f, %f), ST_Point(%f, %f)), %d)) AND \
     285                 :                     %s",
     286                 :                     poWKTRasterDS->pszRasterColumnName,
     287                 :                     poWKTRasterDS->pszSchemaName,
     288                 :                     poWKTRasterDS->pszTableName,
     289                 :                     poWKTRasterDS->pszRasterColumnName, fProjLowerLeftX,
     290                 :                     fProjLowerLeftY, fProjUpperRightX, fProjUpperRightY,
     291               0 :                     poWKTRasterDS->nSrid, poWKTRasterDS->pszWhereClause);
     292                 :         }
     293                 : 
     294                 : 
     295                 :     }        // No where clause
     296                 :     else {
     297                 : 
     298                 : 
     299                 :         /**
     300                 :          * Table has a GIST index.
     301                 :          */
     302               0 :         if (poWKTRasterDS->bTableHasGISTIndex) {
     303                 : 
     304                 :             osCommand.Printf(
     305                 :                     "SELECT rid, %s FROM %s.%s WHERE %s ~ ST_SetSRID(ST_MakeBox2D\
     306                 :                     (ST_Point(%f, %f), ST_Point(%f, %f)), %d)",
     307                 :                     poWKTRasterDS->pszRasterColumnName,
     308                 :                     poWKTRasterDS->pszSchemaName,
     309                 :                     poWKTRasterDS->pszTableName,
     310                 :                     poWKTRasterDS->pszRasterColumnName, fProjLowerLeftX,
     311                 :                     fProjLowerLeftY, fProjUpperRightX, fProjUpperRightY,
     312               0 :                     poWKTRasterDS->nSrid);
     313                 : 
     314                 :         }            /**
     315                 :              * Table hasn't a GIST index. Normal searching
     316                 :              */
     317                 :         else {
     318                 :             osCommand.Printf(
     319                 :                     "SELECT rid, %s FROM %s.%s WHERE _ST_Contains(%s, ST_SetSRID(\
     320                 :                     ST_MakeBox2D(ST_Point(%f, %f), ST_Point(%f, %f)), %d))",
     321                 :                     poWKTRasterDS->pszRasterColumnName,
     322                 :                     poWKTRasterDS->pszSchemaName,
     323                 :                     poWKTRasterDS->pszTableName,
     324                 :                     poWKTRasterDS->pszRasterColumnName, fProjLowerLeftX, 
     325                 :                     fProjLowerLeftY, fProjUpperRightX, fProjUpperRightY,
     326               0 :                     poWKTRasterDS->nSrid);
     327                 :         }
     328                 :     }
     329                 : 
     330                 : 
     331                 :     //printf("query: %s\n", osCommand.c_str());
     332                 : 
     333               0 :     hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
     334               0 :     if (hPGresult == NULL || PQresultStatus(hPGresult) != PGRES_TUPLES_OK) {
     335               0 :         if (hPGresult)
     336               0 :             PQclear(hPGresult);
     337                 : 
     338                 :         CPLError(CE_Failure, CPLE_AppDefined,
     339                 :                 "Sorry, couldn't fetch block information from database: %s",
     340               0 :                 PQerrorMessage(poWKTRasterDS->hPGconn));
     341               0 :         return CE_Failure;
     342                 :     }
     343                 : 
     344                 : 
     345               0 :     nTuples = PQntuples(hPGresult);
     346                 : 
     347                 : 
     348                 :     /*******************************************************
     349                 :      * Block not found. We have to create a new one
     350                 :      *******************************************************/
     351               0 :     if (nTuples <= 0) {
     352                 : 
     353                 :         /**
     354                 :          * There is no block. We need to create a new one. We can use the
     355                 :          * first block of the table and modify it.
     356                 :          */
     357               0 :         PQclear(hPGresult);
     358                 :         osCommand.Printf(
     359                 :                 "SELECT %s FROM %s.%s LIMIT 1 OFFSET 0",
     360                 :                 poWKTRasterDS->pszRasterColumnName,
     361                 :                 poWKTRasterDS->pszSchemaName,
     362               0 :                 poWKTRasterDS->pszTableName);
     363                 : 
     364               0 :         hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
     365                 : 
     366                 :         /**
     367                 :          * Empty table, What should we do? Is it an error?
     368                 :          */
     369               0 :         if (hPGresult == NULL || PQresultStatus(hPGresult) != PGRES_TUPLES_OK) {
     370               0 :             if (hPGresult)
     371               0 :                 PQclear(hPGresult);
     372                 : 
     373                 :             CPLError(CE_Failure, CPLE_AppDefined,
     374                 :                     "Sorry, couldn't fetch block information from database: %s",
     375               0 :                     PQerrorMessage(poWKTRasterDS->hPGconn));
     376                 : 
     377               0 :             return CE_Failure;
     378                 :         }
     379                 :         
     380                 :         // Get HEXWKB representation of the block
     381               0 :         pszHexWkb = CPLStrdup(PQgetvalue(hPGresult, 0, 0));
     382               0 :         PQclear(hPGresult);
     383                 : 
     384                 :         // Create a wrapper object with this data
     385               0 :         poWKTRasterWrapper = new WKTRasterWrapper();
     386                 : 
     387                 :         // Should we try creating the raster block in other way?
     388               0 :         if (poWKTRasterWrapper->Initialize(pszHexWkb) == FALSE) {
     389               0 :             CPLFree(pszHexWkb);
     390               0 :             delete poWKTRasterWrapper; 
     391               0 :             return CE_Failure;
     392                 :         }
     393                 : 
     394                 :         // We won't need this
     395               0 :         CPLFree(pszHexWkb);
     396                 : 
     397                 :         // Get raster band
     398               0 :         poWKTRasterBandWrapper = poWKTRasterWrapper->GetBand((GUInt16)nBand);
     399                 : 
     400                 :         // Should we try creating the raster block in other way?
     401               0 :         if (poWKTRasterBandWrapper == NULL) {
     402               0 :             delete poWKTRasterWrapper;
     403               0 :             return CE_Failure;
     404                 :         }
     405                 : 
     406                 :         // Set raster data
     407                 :         poWKTRasterBandWrapper->SetData((GByte *)pImage,
     408               0 :                 (nBlockXSize * nBlockYSize) * nPixelSize * sizeof (GByte));
     409                 : 
     410                 :         // Insert new block into table. First, we need a new rid
     411                 :         osCommand.Printf(
     412                 :                 "SELECT rid FROM %s.%s ORDER BY rid DESC LIMIT 1 OFFSET 0",
     413               0 :                 poWKTRasterDS->pszSchemaName, poWKTRasterDS->pszTableName);
     414                 : 
     415               0 :         hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
     416               0 :         if (
     417                 :                 hPGresult == NULL ||
     418                 :                 PQresultStatus(hPGresult) != PGRES_TUPLES_OK ||
     419                 :                 PQntuples(hPGresult) <= 0) {
     420                 :             // What should we do?
     421               0 :             if (hPGresult)
     422               0 :                 PQclear(hPGresult);
     423               0 :             delete poWKTRasterWrapper;
     424               0 :             delete poWKTRasterBandWrapper;
     425                 : 
     426               0 :             return CE_Failure;
     427                 :         }
     428                 : 
     429               0 :         int nRid = atoi(PQgetvalue(hPGresult, 0, 0)) + 1;
     430               0 :         pszHexWkb = poWKTRasterWrapper->GetHexWkbRepresentation();
     431                 : 
     432                 :         // Insert the block        
     433                 :         osCommand.Printf(
     434                 :                 "INSERT INTO %s.%s (rid, %s) VALUES (%d, %s)",
     435                 :                 poWKTRasterDS->pszSchemaName,poWKTRasterDS->pszTableName,
     436               0 :                 poWKTRasterDS->pszRasterColumnName,nRid, pszHexWkb);
     437                 : 
     438               0 :         hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
     439               0 :         if (hPGresult == NULL || 
     440                 :                 PQresultStatus(hPGresult) != PGRES_COMMAND_OK) {
     441                 :             CPLError(CE_Failure, CPLE_NoWriteAccess,
     442                 :                     "Couldn't add new block to database: %s",
     443               0 :                     PQerrorMessage(poWKTRasterDS->hPGconn));
     444               0 :             if (hPGresult)
     445               0 :                 PQclear(hPGresult);
     446               0 :             delete poWKTRasterWrapper;
     447               0 :             delete poWKTRasterBandWrapper;
     448               0 :             CPLFree(pszHexWkb);
     449                 : 
     450               0 :             return CE_Failure;
     451                 :         }
     452                 : 
     453                 :         // block added
     454               0 :         CPLFree(pszHexWkb);
     455               0 :         PQclear(hPGresult);
     456                 :     }
     457                 : 
     458                 :    /****************************************************************
     459                 :     * One block found: We have to update the block data
     460                 :     ****************************************************************/
     461               0 :     else if (nTuples == 1) {
     462                 : 
     463                 :         // get data
     464               0 :         nRid = atoi(PQgetvalue(hPGresult, 0, 0));
     465               0 :         pszHexWkb = CPLStrdup(PQgetvalue(hPGresult, 0, 1));
     466               0 :         PQclear(hPGresult);
     467                 : 
     468                 :         // Create wrapper
     469                 :          // Should we try creating the raster block in other way?
     470               0 :         poWKTRasterWrapper = new WKTRasterWrapper();
     471               0 :   if (poWKTRasterWrapper->Initialize(pszHexWkb) == FALSE) {
     472               0 :             CPLFree(pszHexWkb);
     473               0 :             delete poWKTRasterWrapper; 
     474               0 :             return CE_Failure;
     475                 :         }
     476                 : 
     477                 :         // We won't need this
     478               0 :         CPLFree(pszHexWkb);
     479                 : 
     480                 :         // Get raster band
     481               0 :         poWKTRasterBandWrapper = poWKTRasterWrapper->GetBand((GUInt16)nBand);
     482                 : 
     483                 :         // Should we try creating the raster block in other way?
     484                 : 
     485               0 :         if (poWKTRasterBandWrapper == NULL) {
     486               0 :             delete poWKTRasterWrapper;
     487               0 :             return CE_Failure;
     488                 :         }
     489                 : 
     490                 :         /**
     491                 :          * Swap words if needed
     492                 :          */
     493               0 :         if (poWKTRasterWrapper->byEndianess != byMachineEndianess) {
     494                 : 
     495                 :             // Get pixel size of this band
     496               0 :             switch (poWKTRasterBandWrapper->byPixelType) {
     497                 :                 case 0: case 1: case 2: case 3: case 4:
     498               0 :                     nCurrentBandPixelSize = 1;
     499               0 :                     break;
     500                 :                 case 5: case 6: case 9:
     501               0 :                     nCurrentBandPixelSize = 2;
     502               0 :                     break;
     503                 :                 case 7: case 8: case 10:
     504               0 :                     nCurrentBandPixelSize = 4;
     505               0 :                     break;
     506                 :                 case 11:
     507               0 :                     nCurrentBandPixelSize = 8;
     508               0 :                     break;
     509                 :                 default:
     510               0 :                     nCurrentBandPixelSize = 1;
     511                 :             }
     512                 : 
     513                 :             GDALSwapWords((GByte *)pImage, nCurrentBandPixelSize,
     514                 :                     poWKTRasterBandWrapper->nDataSize / nCurrentBandPixelSize,
     515               0 :                     nCurrentBandPixelSize);
     516                 :         }
     517                 : 
     518                 :         // Set raster data
     519                 :         poWKTRasterBandWrapper->SetData((GByte *)pImage,
     520               0 :                 (nBlockXSize * nBlockYSize) * nPixelSize * sizeof (GByte));
     521                 :         
     522                 :         // Get hexwkb again, with new data
     523               0 :         pszHexWkb = poWKTRasterWrapper->GetHexWkbRepresentation();
     524                 : 
     525                 : 
     526                 :         // update register
     527                 :         osCommand.Printf(
     528                 :                 "UPDATE %s.%s SET %s = %s WHERE rid = %d",
     529                 :                 poWKTRasterDS->pszSchemaName, poWKTRasterDS->pszTableName,
     530               0 :                 poWKTRasterDS->pszRasterColumnName, pszHexWkb, nRid);
     531               0 :         hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
     532               0 :         if (hPGresult == NULL ||
     533                 :                 PQresultStatus(hPGresult) != PGRES_COMMAND_OK) {
     534               0 :             if (hPGresult)
     535               0 :                 PQclear(hPGresult);
     536               0 :             CPLFree(pszHexWkb);
     537                 :             CPLError(CE_Failure, CPLE_NoWriteAccess,
     538               0 :                     "Couldn't update the raster data");
     539               0 :             return CE_Failure;
     540                 :         }
     541                 : 
     542                 :         // ok, updated
     543               0 :         CPLFree(pszHexWkb);
     544               0 :         PQclear(hPGresult);
     545                 : 
     546                 :     }
     547                 :    
     548                 :     /*****************************************************************
     549                 :      * More than one block found. What should we do?
     550                 :      *****************************************************************/
     551                 :     else {
     552                 :         // Only regular_block supported, just now
     553               0 :         PQclear(hPGresult);
     554                 :         CPLError(CE_Failure, CPLE_NotSupported,
     555               0 :                 "Sorry, but the raster presents block overlapping. This feature\
     556                 :                 is under development");
     557                 : 
     558               0 :         return CE_Failure;
     559                 :     }
     560                 : 
     561                 : 
     562               0 :     return CE_None;
     563                 : }
     564                 : 
     565                 : 
     566                 : /**
     567                 :  * Read a block of image data
     568                 :  * Inputs:
     569                 :  *  int: horizontal block offset
     570                 :  *  int: vertical block offset
     571                 :  *  void *: The buffer into the data will be read
     572                 :  * Output:
     573                 :  *  CE_None on success, CE_Failure on error
     574                 :  */
     575                 : CPLErr WKTRasterRasterBand::IReadBlock(int nBlockXOff,
     576               0 :         int nBlockYOff, void * pImage) {
     577                 : 
     578               0 :     WKTRasterDataset * poWKTRasterDS = (WKTRasterDataset *) poDS;
     579               0 :     CPLString osCommand;
     580                 :     PGresult * hPGresult;
     581                 :     int nPixelSize;
     582                 :     int nPixelXInitPosition;
     583                 :     int nPixelYInitPosition;
     584                 :     int nPixelXEndPosition;
     585                 :     int nPixelYEndPosition;
     586                 :     double adfTransform[6];
     587               0 :     double dfProjXInit = 0.0;
     588               0 :     double dfProjYInit = 0.0;
     589               0 :     double dfProjXEnd = 0.0;
     590               0 :     double dfProjYEnd = 0.0;
     591               0 :     double dfProjLowerLeftX = 0.0;
     592               0 :     double dfProjLowerLeftY = 0.0;
     593               0 :     double dfProjUpperRightX = 0.0;
     594               0 :     double dfProjUpperRightY = 0.0;
     595               0 :     char * pszHexWkb = NULL;
     596               0 :     int nTuples = 0;
     597               0 :     GByte * pbyRasterData = NULL;
     598               0 :     WKTRasterWrapper * poWKTRasterWrapper = NULL;
     599               0 :     WKTRasterBandWrapper * poWKTRasterBandWrapper = NULL;
     600               0 :     int nNaturalXBlockSize = 0;
     601               0 :     int nNaturalYBlockSize = 0;
     602               0 :     int nPadXSize = 0;
     603               0 :     int nPadYSize = 0;
     604               0 :     int nBlockXBound = 0;
     605               0 :     int nBlockYBound = 0;
     606                 : 
     607                 :     /* Check input parameters */
     608               0 :     if (pImage == NULL || nBlockXOff < 0 || nBlockYOff < 0) {
     609                 :         CPLError(CE_Failure, CPLE_NotSupported,
     610               0 :                 "Unsupported block size or NULL buffer");
     611               0 :         return CE_Failure;
     612                 :     }
     613                 : 
     614                 :     /*************************************************************************
     615                 :      * Get pixel size (divide by 8 because GDALGetDataTypeSize returns the
     616                 :      * size in bits)
     617                 :      *************************************************************************/
     618               0 :     nPixelSize = MAX(1,GDALGetDataTypeSize(eDataType)/8);
     619                 : 
     620                 :     /*************************************************************************
     621                 :      * nBlockXOff and nBlockYOff are block offsets. So, we first have to
     622                 :      * transform them in pixel/line coordinates, taking into account the
     623                 :      * size of a block.
     624                 :      *************************************************************************/
     625               0 :     GetBlockSize(&nNaturalXBlockSize, &nNaturalYBlockSize);
     626                 : 
     627                 :     /**
     628                 :      * The end of this block is the start of the next one
     629                 :      * xxx jorgearevalo: sure?? always??
     630                 :      */
     631               0 :     nBlockXBound = (nBlockXOff * nNaturalXBlockSize) + nNaturalXBlockSize;
     632               0 :     nBlockYBound = (nBlockYOff * nNaturalYBlockSize) + nNaturalYBlockSize;
     633                 : 
     634               0 :     if (nBlockXBound > nRasterXSize)
     635               0 :         nPadXSize = nBlockXBound - nRasterXSize;            
     636               0 :     if (nBlockYBound > nRasterYSize)
     637               0 :         nPadYSize = nBlockYBound - nRasterYSize;
     638                 : 
     639                 : 
     640               0 :     nPixelXInitPosition = nBlockXOff * nNaturalXBlockSize;
     641               0 :     nPixelYInitPosition = nBlockYOff * nNaturalYBlockSize;
     642                 : 
     643               0 :     nPixelXEndPosition = nPixelXInitPosition + (nNaturalXBlockSize - nPadXSize);
     644               0 :     nPixelYEndPosition = nPixelYInitPosition + (nNaturalYBlockSize - nPadYSize);
     645                 : 
     646                 :     /**************************************************************************
     647                 :      * Transform pixel/line coordinates into coordinates of the raster
     648                 :      * reference system.
     649                 :      * NOTE: I take the georeference information from dataset. The SQL function
     650                 :      * ST_GdalGeoTransform takes the same information from the raster row, in
     651                 :      * array format.
     652                 :      **************************************************************************/
     653               0 :     poWKTRasterDS->GetGeoTransform(adfTransform);
     654                 :    
     655                 :     dfProjXInit = adfTransform[0] +
     656                 :             nPixelXInitPosition * adfTransform[1] +
     657               0 :             nPixelYInitPosition * adfTransform[2];
     658                 : 
     659                 :     dfProjYInit = adfTransform[3] +
     660                 :             nPixelXInitPosition * adfTransform[4] +
     661               0 :             nPixelYInitPosition * adfTransform[5];
     662                 : 
     663                 :     dfProjXEnd = adfTransform[0] +
     664                 :             nPixelXEndPosition * adfTransform[1] +
     665               0 :             nPixelYEndPosition * adfTransform[2];
     666                 : 
     667                 :     dfProjYEnd = adfTransform[3] +
     668                 :             nPixelXEndPosition * adfTransform[4] +
     669               0 :             nPixelYEndPosition * adfTransform[5];
     670                 : 
     671                 : 
     672                 :     /*************************************************************************
     673                 :      * Now we have the block coordinates transformed into coordinates of the
     674                 :      * raster reference system. This coordinates are from:
     675                 :      *  - Upper left corner
     676                 :      *  - Lower right corner
     677                 :      * But for ST_MakeBox2D, we'll need block's coordinates of:
     678                 :      *  - Lower left corner
     679                 :      *  - Upper right corner
     680                 :      *************************************************************************/
     681               0 :     dfProjLowerLeftX = dfProjXInit;
     682               0 :     dfProjLowerLeftY = dfProjYEnd;
     683                 : 
     684               0 :     dfProjUpperRightX = dfProjXEnd;
     685               0 :     dfProjUpperRightY = dfProjYInit;
     686                 : 
     687                 : 
     688                 :     /**************************************************************************
     689                 :      * Perform a spatial query that gives the tile/block (row of raster table)
     690                 :      * or tiles/blocks (in case of non-regular blocking) that contain this block
     691                 :      **************************************************************************/
     692               0 :     if (poWKTRasterDS->pszWhereClause != NULL)
     693                 :     {
     694                 :        osCommand.Printf(
     695                 :                "SELECT %s FROM %s.%s WHERE %s ~ ST_SetSRID(ST_MakeBox2D\
     696                 :                 (ST_Point(%f, %f), ST_Point(%f, %f)), %d) AND %s",
     697                 :                poWKTRasterDS->pszRasterColumnName, poWKTRasterDS->pszSchemaName,
     698                 :                poWKTRasterDS->pszTableName, poWKTRasterDS->pszRasterColumnName,
     699                 :                dfProjLowerLeftX, dfProjLowerLeftY, dfProjUpperRightX,
     700                 :                dfProjUpperRightY, poWKTRasterDS->nSrid,
     701               0 :                poWKTRasterDS->pszWhereClause); 
     702                 :     }
     703                 : 
     704                 :     else
     705                 :     {
     706                 :        osCommand.Printf(
     707                 :                "SELECT %s FROM %s.%s WHERE %s ~ ST_SetSRID(ST_MakeBox2D\
     708                 :                 (ST_Point(%f, %f), ST_Point(%f, %f)), %d)",
     709                 :                poWKTRasterDS->pszRasterColumnName, poWKTRasterDS->pszSchemaName,
     710                 :                poWKTRasterDS->pszTableName, poWKTRasterDS->pszRasterColumnName,
     711                 :                dfProjLowerLeftX, dfProjLowerLeftY, dfProjUpperRightX,
     712               0 :                dfProjUpperRightY, poWKTRasterDS->nSrid);
     713                 :  
     714                 :     }
     715                 : 
     716               0 :     hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
     717                 : 
     718               0 :     if (hPGresult == NULL ||
     719                 :             PQresultStatus(hPGresult) != PGRES_TUPLES_OK ||
     720                 :             PQntuples(hPGresult) < 0)
     721                 :     {
     722               0 :         if (hPGresult)
     723               0 :             PQclear(hPGresult);
     724                 :         CPLError(CE_Failure, CPLE_AppDefined,
     725                 :                 "Sorry, couldn't fetch block information from database: %s",
     726               0 :                 PQerrorMessage(poWKTRasterDS->hPGconn));
     727               0 :         return CE_Failure;
     728                 :     }
     729                 : 
     730                 : 
     731               0 :     nTuples = PQntuples(hPGresult);
     732                 : 
     733                 : 
     734                 :     /*****************************************************************
     735                 :      * No blocks found. Fill the buffer with nodata value
     736                 :      *****************************************************************/
     737               0 :     if (nTuples == 0) 
     738                 :     {
     739               0 :         NullBlock(pImage);
     740               0 :         return CE_None;
     741                 :     }
     742                 : 
     743                 :     
     744                 :     /******************************************************************
     745                 :      * One block found. Regular blocking arrangements, no overlaps
     746                 :      ******************************************************************/
     747               0 :     else if (nTuples == 1) 
     748                 :     {
     749                 :         // Get HEXWKB representation of the block
     750               0 :         pszHexWkb = CPLStrdup(PQgetvalue(hPGresult, 0, 0));
     751                 : 
     752               0 :         PQclear(hPGresult);
     753                 : 
     754                 :         // Raster hex must have an even number of characters
     755               0 :         if (pszHexWkb == NULL || strlen(pszHexWkb) % 2) 
     756                 :         {
     757                 :             CPLError(CE_Failure, CPLE_AppDefined,
     758               0 :                     "The HEXWKB data fetch from database must have an even number \
     759                 :                   of characters");
     760               0 :             return CE_Failure;
     761                 :         }
     762                 : 
     763                 : 
     764                 :         // Create a wrapper object
     765               0 :         poWKTRasterWrapper = new WKTRasterWrapper();
     766               0 :         if (poWKTRasterWrapper->Initialize(pszHexWkb) == FALSE) 
     767                 :         {
     768               0 :             CPLFree(pszHexWkb);
     769               0 :             delete poWKTRasterWrapper;
     770               0 :             return CE_Failure;
     771                 :         }
     772                 :         
     773                 : 
     774                 :         // We won't need this
     775               0 :         CPLFree(pszHexWkb);
     776                 : 
     777                 :         // Create raster band wrapper
     778               0 :         poWKTRasterBandWrapper = poWKTRasterWrapper->GetBand((GUInt16)nBand);
     779               0 :         if (poWKTRasterBandWrapper == NULL) 
     780                 :         {
     781               0 :             CPLError(CE_Failure, CPLE_ObjectNull,"Couldn't fetch band data");
     782               0 :             delete poWKTRasterWrapper;
     783               0 :             return CE_Failure;
     784                 :         }
     785                 : 
     786                 :         // Get raster data
     787               0 :         pbyRasterData = poWKTRasterBandWrapper->GetData();
     788                 : 
     789                 :        
     790                 :         
     791                 :         //printf("IReadBlock with offset: %d, %d\n", nBlockXOff, nBlockYOff);
     792                 : 
     793                 :         /**********************************************************
     794                 :          * Check if the raster is offline
     795                 :          **********************************************************/
     796               0 :         if (poWKTRasterBandWrapper->bIsOffline == TRUE) 
     797                 :         {
     798                 :             // The raster data in this case is a path to the raster file
     799               0 :             int nBandToRead = poWKTRasterBandWrapper->nOutDbBandNumber;
     800                 : 
     801                 : 
     802                 :             // Open dataset, if needed
     803               0 :             if  (poWKTRasterDS->poOutdbRasterDS == NULL) 
     804                 :             {
     805                 :                 poWKTRasterDS->poOutdbRasterDS = (GDALDataset *)
     806               0 :                         GDALOpen((char *)pbyRasterData, GA_ReadOnly);
     807                 :             }
     808                 : 
     809               0 :             if (poWKTRasterDS->poOutdbRasterDS != NULL)
     810                 :             {
     811                 : 
     812                 :                 // Read data from band
     813                 :                 /**
     814                 :                  * NOT SO SIMPLE!!!!
     815                 :                  * The outdb file may don't have the same block structure...
     816                 :                  */
     817                 : 
     818                 :                 
     819                 :                 poWKTRasterDS->poOutdbRasterDS->GetRasterBand(nBandToRead)->ReadBlock(
     820               0 :                     nBlockXOff, nBlockYOff, pImage);                             
     821                 :             }
     822                 :             else 
     823                 :             {
     824                 :                 CPLError(CE_Failure, CPLE_ObjectNull,
     825               0 :                     "Couldn't read band data from out-db raster");
     826               0 :                 delete poWKTRasterWrapper;
     827               0 :                 return CE_Failure;
     828                 :             }
     829                 :             
     830                 :         }
     831                 : 
     832                 : 
     833                 :         /****************************
     834                 :          * Indb raster
     835                 :          ****************************/
     836                 :         else 
     837                 :         {
     838                 :             /**
     839                 :              * Copy the data buffer into pImage.
     840                 :              * nBlockXSize * nBlockYSize should be equal to nDataSize/2
     841                 :              */
     842                 :             memcpy(pImage, pbyRasterData,
     843                 :                     (nNaturalXBlockSize * nNaturalYBlockSize) *
     844               0 :                     nPixelSize * sizeof (GByte));
     845                 :         }
     846                 :        
     847                 : 
     848                 : 
     849                 :         // Free resources and exit
     850               0 :         delete poWKTRasterWrapper;
     851                 : 
     852               0 :         return CE_None;
     853                 : 
     854                 :     }// end if ntuples == 1
     855                 : 
     856                 :         /*********************************************************************
     857                 :          * More than one block found. Non regular blocking arrangements
     858                 :          *********************************************************************/
     859                 :     else 
     860                 :     {
     861                 :         // Only regular_block supported, just now
     862               0 :         PQclear(hPGresult);
     863                 :         CPLError(CE_Failure, CPLE_NotSupported,
     864               0 :                 "Sorry, but the raster presents block overlapping. This feature \
     865                 :                 is under development");
     866                 : 
     867               0 :         return CE_Failure;
     868               0 :     }
     869                 : }
     870                 : 
     871                 : /**
     872                 :  * Set the block data to the null value if it is set, or zero if there is
     873                 :  * no null data value.
     874                 :  * Parameters:
     875                 :  *  - void *: the block data
     876                 :  * Returns: nothing
     877                 :  */
     878               0 : void WKTRasterRasterBand::NullBlock(void *pData) 
     879                 : {
     880               0 :     VALIDATE_POINTER0(pData, "NullBlock");
     881                 : 
     882               0 :     int nNaturalBlockXSize = 0;
     883               0 :     int nNaturalBlockYSize = 0;
     884               0 :     GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
     885                 : 
     886               0 :     int nWords = nNaturalBlockXSize * nNaturalBlockYSize;
     887               0 :     int nChunkSize = MAX(1, GDALGetDataTypeSize(eDataType) / 8);
     888                 : 
     889                 :     int bNoDataSet;
     890               0 :     double dfNoData = GetNoDataValue(&bNoDataSet);
     891               0 :     if (!bNoDataSet) 
     892                 :     {
     893               0 :         memset(pData, 0, nWords * nChunkSize);
     894                 :     } 
     895                 :     else 
     896                 :     {
     897               0 :         int i = 0;
     898               0 :         for (i = 0; i < nWords; i += nChunkSize)
     899               0 :             memcpy((GByte *) pData + i, &dfNoData, nChunkSize);
     900                 :     }
     901                 : }
     902                 : 
     903                 : /**
     904                 :  * Set the no data value for this band.
     905                 :  * Parameters:
     906                 :  *  - double: The nodata value
     907                 :  * Returns:
     908                 :  *  - CE_None.
     909                 :  */
     910               0 : CPLErr WKTRasterRasterBand::SetNoDataValue(double dfNewValue) {
     911               0 :     dfNoDataValue = dfNewValue;
     912                 : 
     913               0 :     return CE_None;
     914                 : }
     915                 : 
     916                 : /**
     917                 :  * Fetch the no data value for this band.
     918                 :  * Parameters:
     919                 :  *  - int *: pointer to a boolean to use to indicate if a value is actually
     920                 :  *          associated with this layer. May be NULL (default).
     921                 :  *  Returns:
     922                 :  *  - double: the nodata value for this band.
     923                 :  */
     924               0 : double WKTRasterRasterBand::GetNoDataValue(int *pbSuccess) {
     925               0 :     if (pbSuccess != NULL)
     926               0 :         *pbSuccess = TRUE;
     927                 : 
     928               0 :     return dfNoDataValue;
     929                 : }
     930                 : 
     931                 : /**
     932                 :  * Returns the number of overview layers available.
     933                 :  * Parameters: none
     934                 :  * Returns:
     935                 :  *  int: the number of overviews layers available
     936                 :  */
     937               0 : int WKTRasterRasterBand::GetOverviewCount() {
     938               0 :     WKTRasterDataset * poWKTRasterDS = (WKTRasterDataset *) poDS;
     939                 : 
     940                 :     return (poWKTRasterDS->nOverviews > 0) ?
     941                 :             poWKTRasterDS->nOverviews :
     942               0 :             GDALRasterBand::GetOverviewCount();
     943                 : }
     944                 : 
     945                 : /**
     946                 :  * Fetch overview raster band object.
     947                 :  * Parameters:
     948                 :  *  - int: overview index between 0 and GetOverviewCount()-1
     949                 :  * Returns:
     950                 :  *  - GDALRasterBand *: overview GDALRasterBand.
     951                 :  */
     952               0 : GDALRasterBand * WKTRasterRasterBand::GetOverview(int nOverview) {
     953               0 :     WKTRasterDataset * poWKTRasterDS = (WKTRasterDataset *) poDS;
     954                 : 
     955               0 :     if (poWKTRasterDS->nOverviews > 0) {
     956               0 :         if (nOverview < 0 || nOverview >= poWKTRasterDS->nOverviews)
     957               0 :             return NULL;
     958                 :         else
     959                 :             return
     960               0 :             poWKTRasterDS->papoWKTRasterOv[nOverview]->GetRasterBand(nBand);
     961                 :     } else
     962               0 :         return GDALRasterBand::GetOverview(nOverview);
     963                 : }
     964                 : 
     965                 : /**
     966                 :  * Get the natural block size for this band.
     967                 :  * Parameters:
     968                 :  *  - int *: pointer to int to store the natural X block size
     969                 :  *  - int *: pointer to int to store the natural Y block size
     970                 :  * Returns: nothing
     971                 :  */
     972               0 : void WKTRasterRasterBand::GetBlockSize(int * pnXSize, int *pnYSize)
     973                 :  {
     974               0 :     if (nBlockXSize == 0 || nBlockYSize == 0) {
     975                 :         CPLError(CE_Failure, CPLE_AppDefined,
     976               0 :                 "This WKT Raster band has non regular blocking arrangement. \
     977                 :                 This feature is under development");
     978                 : 
     979               0 :         if (pnXSize != NULL)
     980               0 :             *pnXSize = 0;
     981               0 :         if (pnYSize != NULL)
     982               0 :             *pnYSize = 0;
     983                 : 
     984                 :     }
     985                 :     else {
     986               0 :         GDALRasterBand::GetBlockSize(pnXSize, pnYSize);
     987                 :     }
     988               0 : }
     989                 : 
     990                 : 
     991                 : 
     992                 : 
     993                 : 

Generated by: LTP GCOV extension version 1.5