1
/******************************************************************************
2
* File : PostGISRasterRasterBand.cpp
3
* Project: PostGIS Raster driver
4
* Purpose: GDAL Raster Band implementation for PostGIS Raster driver
5
* Author: Jorge Arevalo, jorgearevalo@deimos-space.com
10
******************************************************************************
11
* Copyright (c) 2009 - 2011, Jorge Arevalo, jorge.arevalo@deimos-space.com
13
* Permission is hereby granted, free of charge, to any person obtaining a
14
* copy of this software and associated documentation files (the "Software"),
15
* to deal in the Software without restriction, including without limitation
16
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
17
* and/or sell copies of the Software, and to permit persons to whom the
18
* Software is furnished to do so, subject to the following conditions:
20
* The above copyright notice and this permission notice shall be included
21
* in all copies or substantial portions of the Software.
23
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29
* DEALINGS IN THE SOFTWARE.
30
******************************************************************************/
31
#include "postgisraster.h"
33
#include "ogr_geometry.h"
34
#include "gdal_priv.h"
37
#include "cpl_string.h"
43
* - PostGISRasterDataset *: The Dataset this band belongs to
44
* - int: the band number
45
* - GDALDataType: The data type for this band
46
* - double: The nodata value. Could be any kind of data (GByte, GUInt16,
47
* GInt32...) but the variable has the bigger type.
48
* - GBool: if the data type is signed byte or not. If yes, the SIGNEDBYTE
49
* metadata value is added to IMAGE_STRUCTURE domain
50
* - int: The bit depth, to add NBITS as metadata value in IMAGE_STRUCTURE
53
PostGISRasterRasterBand::PostGISRasterRasterBand(PostGISRasterDataset *poDS,
54
int nBand, GDALDataType hDataType, double dfNodata, GBool bSignedByte,
55
int nBitDepth, int nFactor, GBool bIsOffline, char * pszSchema,
56
char * pszTable, char * pszColumn)
59
/* Basic properties */
62
this->bIsOffline = bIsOffline;
63
this->pszSchema = (pszSchema) ? pszSchema : CPLStrdup(poDS->pszSchema);
64
this->pszTable = (pszTable) ? pszTable: CPLStrdup(poDS->pszTable);
65
this->pszColumn = (pszColumn) ? pszColumn : CPLStrdup(poDS->pszColumn);
66
this->pszWhere = CPLStrdup(poDS->pszWhere);
68
eAccess = poDS->GetAccess();
69
eDataType = hDataType;
70
dfNoDataValue = dfNodata;
73
/*****************************************************
74
* Check block size issue
75
*****************************************************/
76
nBlockXSize = poDS->nBlockXSize;
77
nBlockYSize = poDS->nBlockYSize;
79
if (nBlockXSize == 0 || nBlockYSize == 0) {
80
CPLError(CE_Warning, CPLE_NotSupported,
81
"This band has irregular blocking, but is not supported yet");
85
// Add pixeltype to image structure domain
86
if (bSignedByte == true) {
87
SetMetadataItem("PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
90
// Add NBITS to metadata only for sub-byte types
92
SetMetadataItem("NBITS", CPLString().Printf( "%d", nBitDepth ),
96
nOverviewFactor = nFactor;
98
/**********************************************************
99
* Check overviews. Query RASTER_OVERVIEWS table for
100
* existing overviews, only in case we are on level 0
101
* TODO: can we do this without querying RASTER_OVERVIEWS?
102
* How do we know the number of overviews? Is an inphinite
104
**********************************************************/
105
if (nOverviewFactor == 0) {
108
PGresult * poResult = NULL;
110
int nFetchOvFactor = 0;
111
char * pszOvSchema = NULL;
112
char * pszOvTable = NULL;
113
char * pszOvColumn = NULL;
115
nRasterXSize = poDS->GetRasterXSize();
116
nRasterYSize = poDS->GetRasterYSize();
118
osCommand.Printf("select o_table_name, overview_factor, o_column, "
119
"o_table_schema from raster_overviews where r_table_schema = "
120
"'%s' and r_table_name = '%s' and r_column = '%s'",
121
poDS->pszSchema, poDS->pszTable, poDS->pszColumn);
123
poResult = PQexec(poDS->poConn, osCommand.c_str());
124
if (poResult != NULL && PQresultStatus(poResult) == PGRES_TUPLES_OK &&
125
PQntuples(poResult) > 0) {
127
/* Create overviews */
128
nOverviewCount = PQntuples(poResult);
129
papoOverviews = (PostGISRasterRasterBand **)VSICalloc(nOverviewCount,
130
sizeof(PostGISRasterRasterBand *));
131
if (papoOverviews == NULL) {
132
CPLError(CE_Warning, CPLE_OutOfMemory, "Couldn't create "
133
"overviews for band %d\n", nBand);
138
for(i = 0; i < nOverviewCount; i++) {
140
nFetchOvFactor = atoi(PQgetvalue(poResult, i, 1));
141
pszOvSchema = CPLStrdup(PQgetvalue(poResult, i, 3));
142
pszOvTable = CPLStrdup(PQgetvalue(poResult, i, 0));
143
pszOvColumn = CPLStrdup(PQgetvalue(poResult, i, 2));
146
* NOTE: Overview bands are not considered to be a part of a
147
* dataset, but we use the same dataset for all the overview
148
* bands just for simplification (we'll need to access the table
149
* and schema names). But in method GetDataset, NULL is return
150
* if we're talking about an overview band
152
papoOverviews[i] = new PostGISRasterRasterBand(poDS, nBand,
153
hDataType, dfNodata, bSignedByte, nBitDepth,
154
nFetchOvFactor, bIsOffline, pszOvSchema, pszOvTable, pszOvColumn);
165
papoOverviews = NULL;
171
/************************************
172
* We are in an overview level. Set
173
* raster size to its value
174
************************************/
178
* No overviews inside an overview (all the overviews are from original
182
papoOverviews = NULL;
184
nRasterXSize = (int) floor((double)poDS->GetRasterXSize() / nOverviewFactor);
185
nRasterYSize = (int) floor((double)poDS->GetRasterYSize() / nOverviewFactor);
188
CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand constructor: Band "
189
"created (srid = %d)", poDS->nSrid);
192
/***********************************************
193
* \brief: Band destructor
194
***********************************************/
195
PostGISRasterRasterBand::~PostGISRasterRasterBand()
209
for(i = 0; i < nOverviewCount; i++)
210
delete papoOverviews[i];
212
CPLFree(papoOverviews);
218
* \brief Set the block data to the null value if it is set, or zero if there is
219
* no null data value.
221
* - void *: the block data
224
void PostGISRasterRasterBand::NullBlock(void *pData)
226
VALIDATE_POINTER0(pData, "NullBlock");
228
int nNaturalBlockXSize = 0;
229
int nNaturalBlockYSize = 0;
230
GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
232
int nWords = nNaturalBlockXSize * nNaturalBlockYSize;
233
int nChunkSize = MAX(1, GDALGetDataTypeSize(eDataType) / 8);
236
double dfNoData = GetNoDataValue(&bNoDataSet);
239
memset(pData, 0, nWords * nChunkSize);
244
for (i = 0; i < nWords; i += nChunkSize)
245
memcpy((GByte *) pData + i, &dfNoData, nChunkSize);
250
* \brief Set the no data value for this band.
252
* - double: The nodata value
256
CPLErr PostGISRasterRasterBand::SetNoDataValue(double dfNewValue) {
257
dfNoDataValue = dfNewValue;
263
* \brief Fetch the no data value for this band.
265
* - int *: pointer to a boolean to use to indicate if a value is actually
266
* associated with this layer. May be NULL (default).
268
* - double: the nodata value for this band.
270
double PostGISRasterRasterBand::GetNoDataValue(int *pbSuccess) {
271
if (pbSuccess != NULL)
274
return dfNoDataValue;
279
* \brief Get the natural block size for this band.
281
* - int *: pointer to int to store the natural X block size
282
* - int *: pointer to int to store the natural Y block size
285
void PostGISRasterRasterBand::GetBlockSize(int * pnXSize, int *pnYSize)
287
if (nBlockXSize == 0 || nBlockYSize == 0) {
288
CPLError(CE_Failure, CPLE_AppDefined,
289
"This PostGIS Raster band has non regular blocking arrangement. \
290
This feature is under development");
299
GDALRasterBand::GetBlockSize(pnXSize, pnYSize);
304
/*****************************************************
305
* \brief Fetch the band number
306
*****************************************************/
307
int PostGISRasterRasterBand::GetBand()
309
return (nOverviewFactor) ? 0 : nBand;
312
/*****************************************************
313
* \brief Fetch the owning dataset handle
314
*****************************************************/
315
GDALDataset* PostGISRasterRasterBand::GetDataset()
317
return (nOverviewFactor) ? NULL : poDS;
320
/****************************************************
321
* \brief Check for arbitrary overviews
322
* The datastore can compute arbitrary overviews
323
* efficiently, because the overviews are tables,
324
* like the original raster. The effort is the same.
325
****************************************************/
326
int PostGISRasterRasterBand::HasArbitraryOverviews()
328
return (nOverviewFactor) ? false : true;
331
/***************************************************
332
* \brief Return the number of overview layers available
333
***************************************************/
334
int PostGISRasterRasterBand::GetOverviewCount()
336
return (nOverviewFactor) ?
341
/**********************************************************
342
* \brief Fetch overview raster band object
343
**********************************************************/
344
GDALRasterBand * PostGISRasterRasterBand::GetOverview(int i)
346
return (i >= 0 && i < GetOverviewCount()) ?
347
(GDALRasterBand *)papoOverviews[i] : GDALRasterBand::GetOverview(i);
350
/*****************************************************
351
* \brief Read a natural block of raster band data
352
*****************************************************/
353
CPLErr PostGISRasterRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void*
356
PostGISRasterDataset * poPostGISRasterDS = (PostGISRasterDataset*)poDS;
358
PGresult* poResult = NULL;
359
double adfTransform[6];
361
int nNaturalBlockXSize = 0;
362
int nNaturalBlockYSize = 0;
364
int nPixelXInitPosition;
365
int nPixelYInitPosition;
366
int nPixelXEndPosition;
367
int nPixelYEndPosition;
368
double dfProjXInit = 0.0;
369
double dfProjYInit = 0.0;
370
double dfProjXEnd = 0.0;
371
double dfProjYEnd = 0.0;
372
double dfProjLowerLeftX = 0.0;
373
double dfProjLowerLeftY = 0.0;
374
double dfProjUpperRightX = 0.0;
375
double dfProjUpperRightY = 0.0;
376
GByte* pbyData = NULL;
377
char* pbyDataToRead = NULL;
379
int nExpectedDataSize = 0;
381
/* Get pixel and block size */
382
nPixelSize = MAX(1,GDALGetDataTypeSize(eDataType)/8);
383
GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
385
/* Get pixel,line coordinates of the block */
387
* TODO: What if the georaster is rotated? Following the gdal_translate
388
* source code, you can't use -projwin option with rotated rasters
390
nPixelXInitPosition = nNaturalBlockXSize * nBlockXOff;
391
nPixelYInitPosition = nNaturalBlockYSize * nBlockYOff;
392
nPixelXEndPosition = nPixelXInitPosition + nNaturalBlockXSize;
393
nPixelYEndPosition = nPixelYInitPosition + nNaturalBlockYSize;
395
poPostGISRasterDS->GetGeoTransform(adfTransform);
397
/* Pixel size correction, in case of overviews */
398
if (nOverviewFactor) {
399
adfTransform[1] *= nOverviewFactor;
400
adfTransform[5] *= nOverviewFactor;
403
/* Calculate the "query box" */
404
dfProjXInit = adfTransform[0] +
405
nPixelXInitPosition * adfTransform[1] +
406
nPixelYInitPosition * adfTransform[2];
408
dfProjYInit = adfTransform[3] +
409
nPixelXInitPosition * adfTransform[4] +
410
nPixelYInitPosition * adfTransform[5];
412
dfProjXEnd = adfTransform[0] +
413
nPixelXEndPosition * adfTransform[1] +
414
nPixelYEndPosition * adfTransform[2];
416
dfProjYEnd = adfTransform[3] +
417
nPixelXEndPosition * adfTransform[4] +
418
nPixelYEndPosition * adfTransform[5];
420
/*************************************************************************
421
* Now we have the block coordinates transformed into coordinates of the
422
* raster reference system. This coordinates are from:
423
* - Upper left corner
424
* - Lower right corner
425
* But for ST_MakeBox2D, we'll need block's coordinates of:
426
* - Lower left corner
427
* - Upper right corner
428
*************************************************************************/
429
dfProjLowerLeftX = dfProjXInit;
430
dfProjLowerLeftY = dfProjYEnd;
432
dfProjUpperRightX = dfProjXEnd;
433
dfProjUpperRightY = dfProjYInit;
436
* Return all raster objects from our raster table that fall reside or
437
* partly reside in a coordinate bounding box.
438
* NOTE: Is it necessary to use a BB optimization like:
439
* select st_makebox2d(...) && geom and (rest of the query)...?
441
if (poPostGISRasterDS->pszWhere != NULL)
443
osCommand.Printf("select rid, %s from %s.%s where %s ~ "
444
"st_setsrid(st_makebox2d(st_point(%f, %f), st_point(%f,"
445
"%f)),%d) and %s", pszColumn, pszSchema, pszTable, pszColumn,
446
dfProjLowerLeftX, dfProjLowerLeftY, dfProjUpperRightX,
447
dfProjUpperRightY, poPostGISRasterDS->nSrid, pszWhere);
452
osCommand.Printf("select rid, %s from %s.%s where %s ~ "
453
"st_setsrid(st_makebox2d(st_point(%f, %f), st_point(%f,"
454
"%f)),%d)", pszColumn, pszSchema, pszTable, pszColumn,
455
dfProjLowerLeftX, dfProjLowerLeftY, dfProjUpperRightX,
456
dfProjUpperRightY, poPostGISRasterDS->nSrid);
459
CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IReadBlock: "
460
"The query = %s", osCommand.c_str());
462
poResult = PQexec(poPostGISRasterDS->poConn, osCommand.c_str());
463
if (poResult == NULL || PQresultStatus(poResult) != PGRES_TUPLES_OK ||
464
PQntuples(poResult) <= 0)
469
/* TODO: Raise an error and exit? */
470
CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IReadBlock: "
471
"The block (%d, %d) is empty", nBlockXOff, nBlockYOff);
476
nTuples = PQntuples(poResult);
483
* TODO: Manage this situation
487
CPLError(CE_Failure, CPLE_AppDefined, "This raster has outdb storage."
488
"This feature isn\'t still available");
492
int nRid = atoi(PQgetvalue(poResult, 0, 0));
494
/* Only data size, without payload */
495
nExpectedDataSize = nNaturalBlockXSize * nNaturalBlockYSize *
497
pbyData = CPLHexToBinary(PQgetvalue(poResult, 0, 1), &nWKBLength);
498
pbyDataToRead = (char*)GET_BAND_DATA(pbyData,nBand, nPixelSize,
501
memcpy(pImage, pbyDataToRead, nExpectedDataSize * sizeof(char));
503
CPLDebug("PostGIS_Raster", "IReadBlock: Copied %d bytes from block "
504
"(%d, %d) (rid = %d) to %p", nExpectedDataSize, nBlockXOff,
505
nBlockYOff, nRid, pImage);
513
/** Overlapping raster data.
514
* TODO: Manage this situation. Suggestion: open several datasets, because
515
* you can't manage overlapping raster data with only one dataset
519
CPLError(CE_Failure, CPLE_AppDefined,
520
"Overlapping raster data. Feature under development, not "
added
removed
frmts/postgisraster/postgisrasterrasterband.cpp
