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
6 : *
7 : * Last changes:
8 : * $Id: $
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2009 - 2011, Jorge Arevalo, jorge.arevalo@deimos-space.com
12 : *
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:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
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"
32 : #include "ogr_api.h"
33 : #include "ogr_geometry.h"
34 : #include "gdal_priv.h"
35 : #include "gdal.h"
36 : #include <string>
37 : #include "cpl_string.h"
38 : #include "gdal_vrt.h"
39 : #include "vrtdataset.h"
40 : #include "memdataset.h"
41 :
42 :
43 : /**
44 : * \brief Constructor.
45 : * Parameters:
46 : * - PostGISRasterDataset *: The Dataset this band belongs to
47 : * - int: the band number
48 : * - GDALDataType: The data type for this band
49 : * - double: The nodata value. Could be any kind of data (GByte, GUInt16,
50 : * GInt32...) but the variable has the bigger type.
51 : * - GBool: if the data type is signed byte or not. If yes, the SIGNEDBYTE
52 : * metadata value is added to IMAGE_STRUCTURE domain
53 : * - int: The bit depth, to add NBITS as metadata value in IMAGE_STRUCTURE
54 : * domain.
55 : * TODO: Comment the rest of parameters
56 : */
57 0 : PostGISRasterRasterBand::PostGISRasterRasterBand(PostGISRasterDataset *poDS,
58 : int nBand, GDALDataType hDataType, GBool bHasNoDataValue, double dfNodata,
59 : GBool bSignedByte,int nBitDepth, int nFactor, int nBlockXSize, int nBlockYSize,
60 0 : GBool bIsOffline, char * inPszSchema, char * inPszTable, char * inPszColumn)
61 : {
62 :
63 : /* Basic properties */
64 0 : this->poDS = poDS;
65 0 : this->nBand = nBand;
66 0 : this->bIsOffline = bIsOffline;
67 :
68 0 : eAccess = poDS->GetAccess();
69 0 : eDataType = hDataType;
70 0 : this->bHasNoDataValue = bHasNoDataValue;
71 0 : dfNoDataValue = dfNodata;
72 :
73 0 : if (poDS->nBands == 1) {
74 0 : eBandInterp = GCI_GrayIndex;
75 : }
76 :
77 0 : else if (poDS->nBands == 3) {
78 0 : if (nBand == 1)
79 0 : eBandInterp = GCI_RedBand;
80 0 : else if( nBand == 2 )
81 0 : eBandInterp = GCI_GreenBand;
82 0 : else if( nBand == 3 )
83 0 : eBandInterp = GCI_BlueBand;
84 : else
85 0 : eBandInterp = GCI_Undefined;
86 : }
87 :
88 : else {
89 0 : eBandInterp = GCI_Undefined;
90 : }
91 :
92 :
93 :
94 : /**************************************************************************
95 : * TODO: Set a non arbitrary blocksize. In case or regular blocking, this is
96 : * easy, otherwise, does it make any sense? Any single tile has its own
97 : * dimensions.
98 : *************************************************************************/
99 0 : if (poDS->bRegularBlocking) {
100 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::Constructor: "
101 0 : "Band %d has regular blocking", nBand);
102 :
103 0 : this->nBlockXSize = nBlockXSize;
104 0 : this->nBlockYSize = nBlockYSize;
105 : }
106 :
107 : else {
108 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::Constructor: "
109 0 : "Band %d does not have regular blocking", nBand);
110 :
111 : /*
112 : this->nBlockXSize = MIN(poDS->nRasterXSize, DEFAULT_BLOCK_X_SIZE);
113 : this->nBlockYSize = MIN(poDS->nRasterYSize, DEFAULT_BLOCK_Y_SIZE);
114 : */
115 :
116 0 : this->nBlockXSize = poDS->nRasterXSize;
117 0 : this->nBlockYSize = 1;
118 : }
119 :
120 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::Constructor: "
121 0 : "Block size (%dx%d)", this->nBlockXSize, this->nBlockYSize);
122 :
123 : // Add pixeltype to image structure domain
124 0 : if (bSignedByte == true) {
125 0 : SetMetadataItem("PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
126 : }
127 :
128 : // Add NBITS to metadata only for sub-byte types
129 0 : if (nBitDepth < 8)
130 : SetMetadataItem("NBITS", CPLString().Printf( "%d", nBitDepth ),
131 0 : "IMAGE_STRUCTURE" );
132 :
133 :
134 0 : nOverviewFactor = nFactor;
135 :
136 0 : this->pszSchema = (inPszSchema) ? inPszSchema : poDS->pszSchema;
137 0 : this->pszTable = (inPszTable) ? inPszTable : poDS->pszTable;
138 0 : this->pszColumn = (inPszColumn) ? inPszColumn : poDS->pszColumn;
139 :
140 : /**********************************************************
141 : * Check overviews. Query RASTER_OVERVIEWS table for
142 : * existing overviews, only in case we are on level 0
143 : * TODO: can we do this without querying RASTER_OVERVIEWS?
144 : * How do we know the number of overviews? Is an inphinite
145 : * loop...
146 : **********************************************************/
147 0 : if (nOverviewFactor == 0) {
148 :
149 0 : CPLString osCommand;
150 0 : PGresult * poResult = NULL;
151 0 : int i = 0;
152 0 : int nFetchOvFactor = 0;
153 0 : char * pszOvSchema = NULL;
154 0 : char * pszOvTable = NULL;
155 0 : char * pszOvColumn = NULL;
156 :
157 0 : nRasterXSize = poDS->GetRasterXSize();
158 0 : nRasterYSize = poDS->GetRasterYSize();
159 :
160 : osCommand.Printf("select o_table_name, overview_factor, o_raster_column, "
161 : "o_table_schema from raster_overviews where r_table_schema = "
162 : "'%s' and r_table_name = '%s' and r_raster_column = '%s'",
163 0 : poDS->pszSchema, poDS->pszTable, poDS->pszColumn);
164 :
165 0 : poResult = PQexec(poDS->poConn, osCommand.c_str());
166 0 : if (poResult != NULL && PQresultStatus(poResult) == PGRES_TUPLES_OK &&
167 : PQntuples(poResult) > 0) {
168 :
169 : /* Create overviews */
170 0 : nOverviewCount = PQntuples(poResult);
171 : papoOverviews = (PostGISRasterRasterBand **)VSICalloc(nOverviewCount,
172 0 : sizeof(PostGISRasterRasterBand *));
173 0 : if (papoOverviews == NULL) {
174 : CPLError(CE_Warning, CPLE_OutOfMemory, "Couldn't create "
175 0 : "overviews for band %d\n", nBand);
176 0 : PQclear(poResult);
177 : return;
178 : }
179 :
180 0 : for(i = 0; i < nOverviewCount; i++) {
181 :
182 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::Constructor: "
183 0 : "Creating overview for band %d", nBand);
184 :
185 0 : nFetchOvFactor = atoi(PQgetvalue(poResult, i, 1));
186 0 : pszOvSchema = CPLStrdup(PQgetvalue(poResult, i, 3));
187 0 : pszOvTable = CPLStrdup(PQgetvalue(poResult, i, 0));
188 0 : pszOvColumn = CPLStrdup(PQgetvalue(poResult, i, 2));
189 :
190 : /**
191 : * NOTE: Overview bands are not considered to be a part of a
192 : * dataset, but we use the same dataset for all the overview
193 : * bands just for simplification (we'll need to access the table
194 : * and schema names). But in method GetDataset, NULL is return
195 : * if we're talking about an overview band
196 : */
197 0 : papoOverviews[i] = new PostGISRasterRasterBand(poDS, nBand,
198 : hDataType, bHasNoDataValue, dfNodata, bSignedByte, nBitDepth,
199 : nFetchOvFactor, nBlockXSize, nBlockYSize, bIsOffline, pszOvSchema,
200 0 : pszOvTable, pszOvColumn);
201 :
202 : }
203 :
204 0 : PQclear(poResult);
205 :
206 : }
207 :
208 : else {
209 :
210 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::Constructor: "
211 0 : "Band %d does not have overviews", nBand);
212 0 : nOverviewCount = 0;
213 0 : papoOverviews = NULL;
214 0 : if (poResult)
215 0 : PQclear(poResult);
216 0 : }
217 : }
218 :
219 : /************************************
220 : * We are in an overview level. Set
221 : * raster size to its value
222 : ************************************/
223 : else {
224 :
225 : /*
226 : * No overviews inside an overview (all the overviews are from original
227 : * band
228 : */
229 0 : nOverviewCount = 0;
230 0 : papoOverviews = NULL;
231 :
232 :
233 0 : nRasterXSize = (int) floor((double)poDS->GetRasterXSize() / nOverviewFactor);
234 0 : nRasterYSize = (int) floor((double)poDS->GetRasterYSize() / nOverviewFactor);
235 : }
236 :
237 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand constructor: Band "
238 0 : "created (srid = %d)", poDS->nSrid);
239 :
240 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand constructor: Band "
241 0 : "size: (%d X %d)", nRasterXSize, nRasterYSize);
242 0 : }
243 :
244 : /***********************************************
245 : * \brief: Band destructor
246 : ***********************************************/
247 0 : PostGISRasterRasterBand::~PostGISRasterRasterBand()
248 : {
249 : int i;
250 :
251 0 : if (papoOverviews) {
252 0 : for(i = 0; i < nOverviewCount; i++)
253 0 : delete papoOverviews[i];
254 :
255 0 : CPLFree(papoOverviews);
256 : }
257 0 : }
258 :
259 :
260 : /**
261 : *
262 : **/
263 0 : GDALDataType PostGISRasterRasterBand::TranslateDataType(const char * pszDataType)
264 : {
265 0 : if (EQUALN(pszDataType, "1BB", 3 * sizeof (char)) ||
266 : EQUALN(pszDataType, "2BUI", 4 * sizeof (char)) ||
267 : EQUALN(pszDataType, "4BUI", 4 * sizeof (char)) ||
268 : EQUALN(pszDataType, "8BUI", 4 * sizeof (char)) ||
269 : EQUALN(pszDataType, "8BSI", 4 * sizeof (char)))
270 :
271 0 : return GDT_Byte;
272 :
273 0 : else if (EQUALN(pszDataType, "16BSI", 5 * sizeof (char)))
274 0 : return GDT_Int16;
275 :
276 0 : else if (EQUALN(pszDataType, "16BUI", 5 * sizeof (char)))
277 0 : return GDT_UInt16;
278 :
279 0 : else if (EQUALN(pszDataType, "32BSI", 5 * sizeof (char)))
280 0 : return GDT_Int32;
281 :
282 0 : else if (EQUALN(pszDataType, "32BUI", 5 * sizeof (char)))
283 0 : return GDT_UInt32;
284 :
285 0 : else if (EQUALN(pszDataType, "32BF", 4 * sizeof (char)))
286 0 : return GDT_Float32;
287 :
288 0 : else if (EQUALN(pszDataType, "64BF", 4 * sizeof (char)))
289 0 : return GDT_Float64;
290 :
291 : else
292 0 : return GDT_Unknown;
293 : }
294 :
295 :
296 :
297 : /**
298 : * Read/write a region of image data for this band.
299 : *
300 : * This method allows reading a region of a PostGISRasterBanda into a buffer.
301 : * The write support is still under development
302 : *
303 : * The function fetches all the raster data that intersects with the region
304 : * provided, and store the data in the GDAL cache.
305 : *
306 : * It automatically takes care of data type translation if the data type
307 : * (eBufType) of the buffer is different than that of the PostGISRasterRasterBand.
308 : *
309 : * The nPixelSpace and nLineSpace parameters allow reading into from various
310 : * organization of buffers.
311 : *
312 : * @param eRWFlag Either GF_Read to read a region of data (GF_Write, to write
313 : * a region of data, yet not supported)
314 : *
315 : * @param nXOff The pixel offset to the top left corner of the region of the
316 : * band to be accessed. This would be zero to start from the left side.
317 : *
318 : * @param nYOff The line offset to the top left corner of the region of the band
319 : * to be accessed. This would be zero to start from the top.
320 : *
321 : * @param nXSize The width of the region of the band to be accessed in pixels.
322 : *
323 : * @param nYSize The height of the region of the band to be accessed in lines.
324 : *
325 : * @param pData The buffer into which the data should be read, or from which it
326 : * should be written. This buffer must contain at least
327 : * nBufXSize * nBufYSize * nBandCount words of type eBufType. It is organized in
328 : * left to right,top to bottom pixel order. Spacing is controlled by the
329 : * nPixelSpace, and nLineSpace parameters.
330 : *
331 : * @param nBufXSize the width of the buffer image into which the desired region
332 : * is to be read, or from which it is to be written.
333 : *
334 : * @param nBufYSize the height of the buffer image into which the desired region
335 : * is to be read, or from which it is to be written.
336 : *
337 : * @param eBufType the type of the pixel values in the pData data buffer. The
338 : * pixel values will automatically be translated to/from the
339 : * PostGISRasterRasterBand data type as needed.
340 : *
341 : * @param nPixelSpace The byte offset from the start of one pixel value in pData
342 : * to the start of the next pixel value within a scanline. If defaulted (0) the
343 : * size of the datatype eBufType is used.
344 : *
345 : * @param nLineSpace The byte offset from the start of one scanline in pData to
346 : * the start of the next. If defaulted (0) the size of the datatype
347 : * eBufType * nBufXSize is used.
348 : *
349 : * @return CE_Failure if the access fails, otherwise CE_None.
350 : */
351 :
352 0 : CPLErr PostGISRasterRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
353 : int nXSize, int nYSize, void * pData, int nBufXSize, int nBufYSize,
354 : GDALDataType eBufType, int nPixelSpace, int nLineSpace)
355 : {
356 : double adfTransform[6];
357 : double adfProjWin[8];
358 : int ulx, uly, lrx, lry;
359 0 : CPLString osCommand;
360 0 : PGresult* poResult = NULL;
361 : int iTuplesIndex;
362 0 : int nTuples = 0;
363 0 : GByte* pbyData = NULL;
364 0 : GByte** ppbyBandData = NULL;
365 0 : int nWKBLength = 0;
366 : int nBandDataLength;
367 : int nBandDataSize;
368 : int nBufDataSize;
369 : int nTileWidth;
370 : int nTileHeight;
371 : double dfTileScaleX;
372 : double dfTileScaleY;
373 : double dfTileUpperLeftX;
374 : double dfTileUpperLeftY;
375 0 : char * pszDataType = NULL;
376 : GDALDataType eTileDataType;
377 : int nTileDataTypeSize;
378 : double dfTileBandNoDataValue;
379 : VRTDatasetH vrtDataset;
380 : GDALDataset ** memDatasets;
381 : GDALRasterBandH memRasterBand;
382 : GDALRasterBandH vrtRasterBand;
383 : char szMemOpenInfo[100];
384 : char szTmp[64];
385 : char szTileWidth[64];
386 : char szTileHeight[64];
387 : CPLErr err;
388 0 : PostGISRasterDataset * poPostGISRasterDS = (PostGISRasterDataset*)poDS;
389 : int nSrcXOff, nSrcYOff, nDstXOff, nDstYOff;
390 : int nDstXSize, nDstYSize;
391 : double xRes, yRes;
392 :
393 : /**
394 : * TODO: Write support not implemented yet
395 : **/
396 0 : if (eRWFlag == GF_Write) {
397 : CPLError(CE_Failure, CPLE_NotSupported,
398 0 : "Writing through PostGIS Raster band not supported yet");
399 :
400 0 : return CE_Failure;
401 : }
402 :
403 0 : nBandDataSize = GDALGetDataTypeSize(eDataType) / 8;
404 0 : nBufDataSize = GDALGetDataTypeSize( eBufType ) / 8;
405 :
406 :
407 : /**************************************************************************
408 : * Do we have overviews that would be appropriate to satisfy this request?
409 : *************************************************************************/
410 0 : if( (nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0 ) {
411 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: "
412 : "nBufXSize = %d, nBufYSize = %d, nXSize = %d, nYSize = %d "
413 0 : "- OverviewRasterIO call", nBufXSize, nBufYSize, nXSize, nYSize);
414 0 : if( OverviewRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
415 : nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace ) == CE_None )
416 :
417 0 : return CE_None;
418 : }
419 :
420 : /**************************************************************************
421 : * Get all the raster rows that are intersected by the window requested
422 : *************************************************************************/
423 : // We first construct a polygon to intersect with
424 0 : poPostGISRasterDS->GetGeoTransform(adfTransform);
425 0 : ulx = nXOff;
426 0 : uly = nYOff;
427 0 : lrx = nXOff + nXSize * nBandDataSize;
428 0 : lry = nYOff + nYSize * nBandDataSize;
429 :
430 : // Calculate right pixel resolution
431 : xRes = (nOverviewFactor == 0) ?
432 : adfTransform[GEOTRSFRM_WE_RES] :
433 0 : adfTransform[GEOTRSFRM_WE_RES] * nOverviewFactor;
434 :
435 : yRes = (nOverviewFactor == 0) ?
436 : adfTransform[GEOTRSFRM_NS_RES] :
437 0 : adfTransform[GEOTRSFRM_NS_RES] * nOverviewFactor;
438 :
439 0 : adfProjWin[0] = adfTransform[GEOTRSFRM_TOPLEFT_X] +
440 : ulx * xRes +
441 0 : uly * adfTransform[GEOTRSFRM_ROTATION_PARAM1];
442 0 : adfProjWin[1] = adfTransform[GEOTRSFRM_TOPLEFT_Y] +
443 0 : ulx * adfTransform[GEOTRSFRM_ROTATION_PARAM2] +
444 0 : uly * yRes;
445 0 : adfProjWin[2] = adfTransform[GEOTRSFRM_TOPLEFT_X] +
446 : lrx * xRes +
447 0 : uly * adfTransform[GEOTRSFRM_ROTATION_PARAM1];
448 0 : adfProjWin[3] = adfTransform[GEOTRSFRM_TOPLEFT_Y] +
449 0 : lrx * adfTransform[GEOTRSFRM_ROTATION_PARAM2] +
450 0 : uly * yRes;
451 0 : adfProjWin[4] = adfTransform[GEOTRSFRM_TOPLEFT_X] +
452 : lrx * xRes +
453 0 : lry * adfTransform[GEOTRSFRM_ROTATION_PARAM1];
454 0 : adfProjWin[5] = adfTransform[GEOTRSFRM_TOPLEFT_Y] +
455 0 : lrx * adfTransform[GEOTRSFRM_ROTATION_PARAM2] +
456 0 : lry * yRes;
457 0 : adfProjWin[6] = adfTransform[GEOTRSFRM_TOPLEFT_X] +
458 : ulx * xRes +
459 0 : lry * adfTransform[GEOTRSFRM_ROTATION_PARAM1];
460 0 : adfProjWin[7] = adfTransform[GEOTRSFRM_TOPLEFT_Y] +
461 0 : ulx * adfTransform[GEOTRSFRM_ROTATION_PARAM2] +
462 0 : lry * yRes;
463 :
464 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: "
465 : "Buffer size = (%d, %d), Region size = (%d, %d)",
466 0 : nBufXSize, nBufYSize, nXSize, nYSize);
467 :
468 0 : if (poPostGISRasterDS->pszWhere == NULL) {
469 : osCommand.Printf("SELECT st_band(%s, %d), st_width(%s), st_height(%s), st_bandpixeltype(%s, %d), "
470 : "st_bandnodatavalue(%s, %d), st_scalex(%s), st_scaley(%s), st_upperleftx(%s), st_upperlefty(%s) "
471 : "FROM %s.%s WHERE st_intersects(%s, st_polygonfromtext('POLYGON((%.17f %.17f, %.17f %.17f, "
472 : "%.17f %.17f, %.17f %.17f, %.17f %.17f))', %d))", pszColumn, nBand, pszColumn, pszColumn, pszColumn, nBand, pszColumn,
473 : nBand, pszColumn, pszColumn, pszColumn, pszColumn, pszSchema, pszTable, pszColumn,
474 : adfProjWin[0], adfProjWin[1], adfProjWin[2], adfProjWin[3], adfProjWin[4], adfProjWin[5],
475 0 : adfProjWin[6], adfProjWin[7], adfProjWin[0], adfProjWin[1], poPostGISRasterDS->nSrid);
476 : }
477 :
478 : else {
479 : osCommand.Printf("SELECT st_band(%s, %d), st_width(%s), st_height(%s), st_bandpixeltype(%s, %d), "
480 : "st_bandnodatavalue(%s, %d), st_scalex(%s), st_scaley(%s), st_upperleftx(%s), st_upperlefty(%s) "
481 : "FROM %s.%s WHERE (%s) AND st_intersects(%s, st_polygonfromtext('POLYGON((%.17f %.17f, %.17f %.17f, "
482 : "%.17f %.17f, %.17f %.17f, %.17f %.17f))', %d))", pszColumn, nBand, pszColumn, pszColumn, pszColumn, nBand, pszColumn,
483 : nBand, pszColumn, pszColumn, pszColumn, pszColumn, pszSchema, pszTable, poPostGISRasterDS->pszWhere,
484 : pszColumn, adfProjWin[0], adfProjWin[1], adfProjWin[2], adfProjWin[3], adfProjWin[4], adfProjWin[5],
485 0 : adfProjWin[6], adfProjWin[7], adfProjWin[0], adfProjWin[1], poPostGISRasterDS->nSrid);
486 : }
487 :
488 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): Query = %s", osCommand.c_str());
489 :
490 0 : poResult = PQexec(poPostGISRasterDS->poConn, osCommand.c_str());
491 0 : if (poResult == NULL || PQresultStatus(poResult) != PGRES_TUPLES_OK ||
492 : PQntuples(poResult) < 0) {
493 :
494 0 : if (poResult)
495 0 : PQclear(poResult);
496 :
497 0 : CPLError(CE_Failure, CPLE_AppDefined, "Error retrieving raster data from database");
498 :
499 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): %s",
500 0 : PQerrorMessage(poPostGISRasterDS->poConn));
501 :
502 0 : return CE_Failure;
503 : }
504 :
505 : /**
506 : * No data. Return the buffer filled with nodata values
507 : **/
508 0 : else if (PQntuples(poResult) == 0) {
509 0 : PQclear(poResult);
510 :
511 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): Null block");
512 :
513 0 : memset(pData, dfNoDataValue, nBufDataSize * nBufXSize * nBufYSize);
514 :
515 0 : return CE_None;
516 : }
517 :
518 :
519 0 : nTuples = PQntuples(poResult);
520 :
521 : /**************************************************************************
522 : * Allocate memory for MEM dataset
523 : * TODO: In case of memory error, provide a different alternative
524 : *************************************************************************/
525 0 : memDatasets = (GDALDataset **)VSICalloc(nTuples, sizeof(GDALDataset *));
526 0 : if (!memDatasets) {
527 0 : PQclear(poResult);
528 : CPLError(CE_Failure, CPLE_AppDefined, "Memory error while trying to read band data "
529 0 : "from database");
530 :
531 0 : return CE_Failure;
532 : }
533 :
534 :
535 : /**************************************************************************
536 : * Create an empty in-memory VRT dataset
537 : * TODO: In case of memory error, provide a different alternative
538 : *************************************************************************/
539 0 : vrtDataset = VRTCreate(nXSize, nYSize);
540 0 : if (!vrtDataset) {
541 0 : PQclear(poResult);
542 0 : CPLFree(memDatasets);
543 :
544 : CPLError(CE_Failure, CPLE_AppDefined, "Memory error while trying to read band data "
545 0 : "from database");
546 :
547 0 : return CE_Failure;
548 : }
549 :
550 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: VRT Dataset "
551 0 : "of (%d, %d) created", nXSize, nYSize);
552 :
553 :
554 : // NOTE: Do NOT add a Dataset description, or the VRT file will be written to disk.
555 : // This is a memory only dataset.
556 0 : GDALSetProjection(vrtDataset, GDALGetProjectionRef((GDALDatasetH)this->poDS));
557 0 : GDALSetGeoTransform(vrtDataset, adfTransform);
558 :
559 :
560 : /**
561 : * Create one VRT Raster Band. It will contain the same band of all tiles
562 : * as Simple Sources
563 : **/
564 0 : VRTAddBand(vrtDataset, eDataType, NULL);
565 0 : vrtRasterBand = GDALGetRasterBand(vrtDataset, 1);
566 :
567 :
568 : /**
569 : * Allocate memory for MEM data pointers
570 : **/
571 0 : ppbyBandData = (GByte **)VSICalloc(nTuples, sizeof(GByte *));
572 0 : if (!ppbyBandData) {
573 0 : PQclear(poResult);
574 0 : CPLFree(memDatasets);
575 0 : GDALClose(vrtDataset);
576 :
577 : CPLError(CE_Failure, CPLE_AppDefined, "Memory error while trying to read band data "
578 0 : "from database");
579 :
580 0 : return CE_Failure;
581 : }
582 :
583 : /**************************************************************************
584 : * Now, for each block, create a MEM dataset
585 : * TODO: What if whe have a really BIG amount of data fetched from db? CURSORS
586 : *************************************************************************/
587 0 : for(iTuplesIndex = 0; iTuplesIndex < nTuples; iTuplesIndex++) {
588 : /**
589 : * Fetch data from result
590 : **/
591 0 : pbyData = CPLHexToBinary(PQgetvalue(poResult, iTuplesIndex, 0), &nWKBLength);
592 0 : nTileWidth = atoi(PQgetvalue(poResult, iTuplesIndex, 1));
593 0 : nTileHeight = atoi(PQgetvalue(poResult, iTuplesIndex, 2));
594 0 : pszDataType = CPLStrdup(PQgetvalue(poResult, iTuplesIndex, 3));
595 0 : dfTileBandNoDataValue = atof(PQgetvalue(poResult, iTuplesIndex, 4));
596 0 : dfTileScaleX = atof(PQgetvalue(poResult, iTuplesIndex, 5));
597 0 : dfTileScaleY = atof(PQgetvalue(poResult, iTuplesIndex, 6));
598 0 : dfTileUpperLeftX = atof(PQgetvalue(poResult, iTuplesIndex, 7));
599 0 : dfTileUpperLeftY = atof(PQgetvalue(poResult, iTuplesIndex, 8));
600 :
601 :
602 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: Tile of "
603 : "(%d, %d) px created. With pixel size of (%f, %f) and located at "
604 : "(%f, %f)", nTileWidth, nTileHeight, dfTileScaleX, dfTileScaleY,
605 0 : dfTileUpperLeftX, dfTileUpperLeftY);
606 :
607 : /**
608 : * Calculate some useful parameters
609 : **/
610 0 : eTileDataType = TranslateDataType(pszDataType);
611 0 : nTileDataTypeSize = GDALGetDataTypeSize(eTileDataType) / 8;
612 :
613 0 : nBandDataLength = nTileWidth * nTileHeight * nTileDataTypeSize;
614 0 : ppbyBandData[iTuplesIndex] = (GByte *)
615 0 : VSIMalloc(nBandDataLength * sizeof(GByte));
616 :
617 0 : if (!ppbyBandData[iTuplesIndex]) {
618 : CPLError(CE_Warning, CPLE_AppDefined, "Could not allocate memory for "
619 0 : "MEMDataset, skipping. The result image may contain gaps");
620 0 : continue;
621 : }
622 :
623 : /**
624 : * Get the pointer to the band pixels
625 : **/
626 0 : memcpy(ppbyBandData[iTuplesIndex],
627 : GET_BAND_DATA(pbyData, 1, nTileDataTypeSize, nBandDataLength),
628 0 : nBandDataLength);
629 :
630 : /**
631 : * Create new MEM dataset, based on in-memory array, to hold the pixels.
632 : * The dataset will only have 1 band
633 : **/
634 0 : memset(szTmp, 0, sizeof(szTmp));
635 0 : CPLPrintPointer(szTmp, ppbyBandData[iTuplesIndex], sizeof(szTmp));
636 :
637 0 : memset(szTileWidth, 0, sizeof(szTileWidth));
638 0 : CPLPrintInt32(szTileWidth, (GInt32)nTileWidth, sizeof(nTileWidth));
639 0 : memset(szTileHeight, 0, sizeof(szTileHeight));
640 0 : CPLPrintInt32(szTileHeight, (GInt32)nTileHeight, sizeof(nTileHeight));
641 :
642 0 : memset(szMemOpenInfo, 0, sizeof(szMemOpenInfo));
643 : sprintf(szMemOpenInfo, "MEM:::DATAPOINTER=%s,PIXELS=%d,LINES=%d,DATATYPE=%s",
644 0 : szTmp, nTileWidth, nTileHeight, GDALGetDataTypeName(eTileDataType));
645 :
646 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: MEMDataset "
647 0 : "open info = %s", szMemOpenInfo);
648 :
649 0 : GDALOpenInfo oOpenInfo(szMemOpenInfo, GA_ReadOnly, NULL);
650 :
651 0 : memDatasets[iTuplesIndex] = MEMDataset::Open(&oOpenInfo);
652 0 : if (!memDatasets[iTuplesIndex]) {
653 : CPLError(CE_Warning, CPLE_AppDefined, "Could not create MEMDataset, "
654 0 : "skipping. The result image may contain gaps");
655 0 : continue;
656 : }
657 :
658 0 : GDALSetDescription(memDatasets[iTuplesIndex], szMemOpenInfo);
659 :
660 : /**
661 : * Get MEM raster band, to add it as simple source.
662 : **/
663 0 : memRasterBand = (GDALRasterBandH)memDatasets[iTuplesIndex]->GetRasterBand(1);
664 0 : if (!memRasterBand) {
665 : CPLError(CE_Warning, CPLE_AppDefined, "Could not get MEMRasterBand , "
666 0 : "skipping. The result image may contain gaps");
667 0 : continue;
668 : }
669 :
670 0 : ((MEMRasterBand *)memRasterBand)->SetNoDataValue(dfTileBandNoDataValue);
671 :
672 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: Adding "
673 0 : "VRT Complex Source");
674 :
675 : /**
676 : * Get source and destination windows for the simple source (first check source
677 : * and destination bounding boxes match. Otherwise, skip this data)
678 : **/
679 :
680 0 : if (dfTileUpperLeftX + nTileWidth * dfTileScaleX < adfProjWin[0]) {
681 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: "
682 : "dfTileUpperLeftX = %f, nTileWidth = %d, dfTileScaleX = %f, "
683 : "RasterDataset minx = %f", dfTileUpperLeftX, nTileWidth, dfTileScaleX,
684 0 : adfProjWin[0]);
685 0 : continue;
686 : }
687 :
688 0 : if (dfTileUpperLeftX > adfProjWin[4]) {
689 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: "
690 : "dfTileUpperLeftX = %f, RasterDataset maxx = %f", dfTileUpperLeftX,
691 0 : adfProjWin[4]);
692 0 : continue;
693 : }
694 :
695 0 : if (dfTileUpperLeftY + nTileHeight * dfTileScaleY > adfProjWin[1]) {
696 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: "
697 : "dfTileUpperLeftY = %f, nTileHeight = %d, ns res = %f, "
698 : "RasterDataset maxy = %f", dfTileUpperLeftY, nTileHeight, dfTileScaleY,
699 0 : adfProjWin[1]);
700 0 : continue;
701 : }
702 :
703 0 : if (dfTileUpperLeftY < adfProjWin[5]) {
704 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO: "
705 : "dfTileUpperLeftY = %f, RasterDataset miny = %f", dfTileUpperLeftY,
706 0 : adfProjWin[5]);
707 0 : continue;
708 : }
709 :
710 :
711 0 : if (dfTileUpperLeftX < adfProjWin[0]) {
712 0 : nSrcXOff = (int)((adfProjWin[0] - dfTileUpperLeftX) /
713 0 : dfTileScaleX + 0.5);
714 0 : nDstXOff = 0;
715 : }
716 :
717 : else {
718 0 : nSrcXOff = 0;
719 0 : nDstXOff = (int)(0.5 + (dfTileUpperLeftX - adfProjWin[0]) /
720 0 : xRes);
721 : }
722 :
723 0 : if (adfProjWin[1] < dfTileUpperLeftY) {
724 0 : nSrcYOff = (int)((dfTileUpperLeftY - adfProjWin[1]) /
725 0 : fabs(dfTileScaleY) + 0.5);
726 0 : nDstYOff = 0;
727 : }
728 :
729 : else {
730 0 : nSrcYOff = 0;
731 0 : nDstYOff = (int)(0.5 + (adfProjWin[1] - dfTileUpperLeftY) /
732 0 : fabs(yRes));
733 : }
734 :
735 0 : nDstXSize = (int)(0.5 + nTileWidth * dfTileScaleX / xRes);
736 0 : nDstYSize = (int)(0.5 + nTileHeight * fabs(dfTileScaleY) / fabs(yRes));
737 :
738 :
739 : /**
740 : * Add the mem raster band as new complex source band (so, I can specify a nodata value)
741 : **/
742 : VRTAddComplexSource(vrtRasterBand, memRasterBand, nSrcXOff, nSrcYOff, nTileWidth, nTileHeight,
743 0 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, 0, 1, dfTileBandNoDataValue);
744 :
745 0 : CPLFree(pbyData);
746 0 : CPLFree(pszDataType);
747 :
748 :
749 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): VRT complex source added");
750 : }
751 :
752 0 : PQclear(poResult);
753 :
754 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): VRT dataset created");
755 :
756 : /**
757 : * We've constructed the VRT Dataset based on the window requested. So, we always
758 : * start from 0
759 : **/
760 0 : nXOff = nYOff = 0;
761 :
762 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): The window requested is "
763 : "from (%d, %d) of size (%d, %d). Buffer of size (%d, %d)", nXOff, nYOff,
764 0 : nXSize, nYSize, nBufXSize, nBufYSize);
765 :
766 : // Execute VRT RasterIO over the band
767 : err = ((VRTRasterBand *)vrtRasterBand)->RasterIO(eRWFlag, nXOff, nYOff, nXSize,
768 0 : nYSize, pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace);
769 :
770 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): Data read");
771 :
772 0 : GDALClose(vrtDataset);
773 :
774 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): VRTDataset released");
775 :
776 : // Free resources
777 0 : for(iTuplesIndex = 0; iTuplesIndex < nTuples; iTuplesIndex++) {
778 0 : if (ppbyBandData[iTuplesIndex])
779 0 : VSIFree(ppbyBandData[iTuplesIndex]);
780 0 : delete memDatasets[iTuplesIndex];
781 : //GDALClose(memDatasets[iTuplesIndex]);
782 : }
783 0 : VSIFree(ppbyBandData);
784 0 : VSIFree(memDatasets);
785 :
786 0 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IRasterIO(): MEMDatasets were released");
787 :
788 0 : return err;
789 :
790 : }
791 :
792 : /**
793 : * \brief Set the no data value for this band.
794 : * Parameters:
795 : * - double: The nodata value
796 : * Returns:
797 : * - CE_None.
798 : */
799 0 : CPLErr PostGISRasterRasterBand::SetNoDataValue(double dfNewValue) {
800 0 : dfNoDataValue = dfNewValue;
801 :
802 0 : return CE_None;
803 : }
804 :
805 : /**
806 : * \brief Fetch the no data value for this band.
807 : * Parameters:
808 : * - int *: pointer to a boolean to use to indicate if a value is actually
809 : * associated with this layer. May be NULL (default).
810 : * Returns:
811 : * - double: the nodata value for this band.
812 : */
813 0 : double PostGISRasterRasterBand::GetNoDataValue(int *pbSuccess) {
814 0 : if (pbSuccess != NULL)
815 0 : *pbSuccess = (int) bHasNoDataValue;
816 :
817 0 : return dfNoDataValue;
818 : }
819 :
820 : /***************************************************
821 : * \brief Return the number of overview layers available
822 : ***************************************************/
823 0 : int PostGISRasterRasterBand::GetOverviewCount()
824 : {
825 : return (nOverviewFactor) ?
826 : 0 :
827 0 : nOverviewCount;
828 : }
829 :
830 : /**********************************************************
831 : * \brief Fetch overview raster band object
832 : **********************************************************/
833 0 : GDALRasterBand * PostGISRasterRasterBand::GetOverview(int i)
834 : {
835 0 : return (i >= 0 && i < GetOverviewCount()) ?
836 0 : (GDALRasterBand *)papoOverviews[i] : GDALRasterBand::GetOverview(i);
837 : }
838 :
839 : /*****************************************************
840 : * \brief Read a natural block of raster band data
841 : *****************************************************/
842 0 : CPLErr PostGISRasterRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void*
843 : pImage)
844 : {
845 0 : int nPixelSize = GDALGetDataTypeSize(eDataType)/8;
846 : int nReadXSize, nReadYSize;
847 :
848 0 : if( (nBlockXOff+1) * nBlockXSize > GetXSize() )
849 0 : nReadXSize = GetXSize() - nBlockXOff * nBlockXSize;
850 : else
851 0 : nReadXSize = nBlockXSize;
852 :
853 0 : if( (nBlockYOff+1) * nBlockYSize > GetYSize() )
854 0 : nReadYSize = GetYSize() - nBlockYOff * nBlockYSize;
855 : else
856 0 : nReadYSize = nBlockYSize;
857 :
858 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IReadBlock: "
859 0 : "Calling to IRasterIO");
860 :
861 : return IRasterIO( GF_Read,
862 : nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
863 : nReadXSize, nReadYSize,
864 : pImage, nReadXSize, nReadYSize, eDataType,
865 0 : nPixelSize, nPixelSize * nBlockXSize );
866 : }
867 :
868 : /**
869 : * \brief How should this band be interpreted as color?
870 : * GCI_Undefined is returned when the format doesn't know anything about the
871 : * color interpretation.
872 : **/
873 0 : GDALColorInterp PostGISRasterRasterBand::GetColorInterpretation()
874 : {
875 0 : return eBandInterp;
876 : }
|
