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 18029 2009-11-15 01:54:19Z 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 18029 2009-11-15 01:54:19Z 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 0 : 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 0 : CPLErr WKTRasterRasterBand::IWriteBlock(int nBlockXOff,
139 : 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 0 : fProjXInit = adfTransform[0] +
209 0 : nPixelXInitPosition * adfTransform[1] +
210 0 : nPixelYInitPosition * adfTransform[2];
211 :
212 0 : fProjYInit = adfTransform[3] +
213 0 : nPixelXInitPosition * adfTransform[4] +
214 0 : nPixelYInitPosition * adfTransform[5];
215 :
216 0 : fProjXEnd = adfTransform[0] +
217 0 : nPixelXEndPosition * adfTransform[1] +
218 0 : nPixelYEndPosition * adfTransform[2];
219 :
220 0 : fProjYEnd = adfTransform[3] +
221 0 : 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 : "Sorry, but the raster presents block overlapping. This feature\
556 0 : 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 0 : CPLErr WKTRasterRasterBand::IReadBlock(int nBlockXOff,
576 : 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 : GByte byMachineEndianess = NDR; // by default
597 0 : int nTuples = 0;
598 0 : GByte * pbyRasterData = NULL;
599 0 : WKTRasterWrapper * poWKTRasterWrapper = NULL;
600 0 : WKTRasterBandWrapper * poWKTRasterBandWrapper = NULL;
601 0 : int nNaturalXBlockSize = 0;
602 0 : int nNaturalYBlockSize = 0;
603 :
604 :
605 : // Check parameters
606 0 : if (pImage == NULL || nBlockXOff < 0 || nBlockYOff < 0) {
607 : CPLError(CE_Failure, CPLE_NotSupported,
608 0 : "Unsupported block size or NULL buffer");
609 0 : return CE_Failure;
610 : }
611 :
612 :
613 : // Check machine endianess
614 : #ifdef CPL_LSB
615 0 : byMachineEndianess = NDR;
616 : #else
617 : byMachineEndianess = XDR;
618 : #endif
619 :
620 : /*************************************************************************
621 : * Get pixel size (divide by 8 because GDALGetDataTypeSize returns the
622 : * size in bits)
623 : *************************************************************************/
624 0 : nPixelSize = MAX(1,GDALGetDataTypeSize(eDataType)/8);
625 :
626 : /*************************************************************************
627 : * nBlockXOff and nBlockYOff are block offsets. So, we first have to
628 : * transform them in pixel/line coordinates, taking into account the
629 : * size of a block.
630 : *************************************************************************/
631 0 : GetBlockSize(&nNaturalXBlockSize, &nNaturalYBlockSize);
632 :
633 :
634 0 : nPixelXInitPosition = nNaturalXBlockSize * nBlockXOff;
635 0 : nPixelYInitPosition = nNaturalYBlockSize * nBlockYOff;
636 :
637 : // Now, get the end of the block
638 0 : nPixelXEndPosition = nPixelXInitPosition + nNaturalXBlockSize;
639 0 : nPixelYEndPosition = nPixelYInitPosition + nNaturalYBlockSize;
640 :
641 :
642 : /**************************************************************************
643 : * Transform pixel/line coordinates into coordinates of the raster
644 : * reference system.
645 : * NOTE: I take the georeference information from dataset. The SQL function
646 : * ST_GdalGeoTransform takes the same information from the raster row, in
647 : * array format.
648 : **************************************************************************/
649 0 : poWKTRasterDS->GetGeoTransform(adfTransform);
650 0 : dfProjXInit = adfTransform[0] +
651 0 : nPixelXInitPosition * adfTransform[1] +
652 0 : nPixelYInitPosition * adfTransform[2];
653 :
654 0 : dfProjYInit = adfTransform[3] +
655 0 : nPixelXInitPosition * adfTransform[4] +
656 0 : nPixelYInitPosition * adfTransform[5];
657 :
658 0 : dfProjXEnd = adfTransform[0] +
659 0 : nPixelXEndPosition * adfTransform[1] +
660 0 : nPixelYEndPosition * adfTransform[2];
661 :
662 0 : dfProjYEnd = adfTransform[3] +
663 0 : nPixelXEndPosition * adfTransform[4] +
664 0 : nPixelYEndPosition * adfTransform[5];
665 :
666 :
667 : /*************************************************************************
668 : * Now we have the block coordinates transformed into coordinates of the
669 : * raster reference system. This coordinates are from:
670 : * - Upper left corner
671 : * - Lower right corner
672 : * But for ST_MakeBox2D, we'll need block's coordinates of:
673 : * - Lower left corner
674 : * - Upper right corner
675 : *************************************************************************/
676 0 : dfProjLowerLeftX = dfProjXInit;
677 0 : dfProjLowerLeftY = dfProjYEnd;
678 :
679 0 : dfProjUpperRightX = dfProjXEnd;
680 0 : dfProjUpperRightY = dfProjYInit;
681 :
682 :
683 : /**************************************************************************
684 : * Perform a spatial query that gives the tile/block (row of raster table)
685 : * or tiles/blocks (in case of non-regular blocking) that contain this block
686 : **************************************************************************/
687 :
688 0 : if (poWKTRasterDS->pszWhereClause != NULL) {
689 :
690 : /**
691 : * Table has GIST index-
692 : * We could disable sequential scanning, but for versions of PostGIS
693 : * up to 0.8, this is no necessary. PostGIS >=0.8 is correctly
694 : * integrated with query planner, thus PostgreSQL will use indexes
695 : * whenever appropriate.
696 : * NOTE: We should add version checking here, and disable sequential
697 : * scanning when needed. See OGRDataSource::Open method.
698 : */
699 0 : if (poWKTRasterDS->bTableHasGISTIndex) {
700 :
701 : /**
702 : * This queries for the tiles that contains the given block. When
703 : * regular_blocking rasters, the blocks will have the same size of
704 : * a tile, so, we can use "contain" function and "~" operator.
705 : * But when non-regular_blocking, the only way this will work is
706 : * setting the block size to the smallest tile size. The problem is
707 : * how do we know all the tiles size?
708 : */
709 : osCommand.Printf(
710 : "SELECT rid, %s FROM %s.%s WHERE %s ~ ST_SetSRID(ST_MakeBox2D\
711 : (ST_Point(%f, %f),ST_Point(%f, %f)), %d) AND %s",
712 : poWKTRasterDS->pszRasterColumnName,
713 : poWKTRasterDS->pszSchemaName,
714 : poWKTRasterDS->pszTableName,
715 : poWKTRasterDS->pszRasterColumnName, dfProjLowerLeftX,
716 : dfProjLowerLeftY, dfProjUpperRightX, dfProjUpperRightY,
717 0 : poWKTRasterDS->nSrid, poWKTRasterDS->pszWhereClause);
718 :
719 : } /**
720 : * Table hasn't a GIST index. Normal searching
721 : */
722 : else {
723 : osCommand.Printf(
724 : "SELECT rid, %s FROM %s.%s WHERE _ST_Contains(%s, ST_SetSRID(\
725 : ST_MakeBox2D(ST_Point(%f, %f), ST_Point(%f, %f)), %d)) AND\
726 : %s",
727 : poWKTRasterDS->pszRasterColumnName,
728 : poWKTRasterDS->pszSchemaName,
729 : poWKTRasterDS->pszTableName,
730 : poWKTRasterDS->pszRasterColumnName, dfProjLowerLeftX,
731 : dfProjLowerLeftY, dfProjUpperRightX, dfProjUpperRightY,
732 0 : poWKTRasterDS->nSrid, poWKTRasterDS->pszWhereClause);
733 : }
734 :
735 :
736 : } // No where clause
737 : else {
738 :
739 :
740 :
741 : //printf("raster band: srid from dataset: %d\n", poWKTRasterDS->nSrid);
742 :
743 : /**
744 : * Table has a GIST index.
745 : */
746 0 : if (poWKTRasterDS->bTableHasGISTIndex) {
747 :
748 : osCommand.Printf(
749 : "SELECT rid, %s FROM %s.%s WHERE %s ~ ST_SetSRID(ST_MakeBox2D(\
750 : ST_Point(%f, %f),ST_Point(%f, %f)), %d)",
751 : poWKTRasterDS->pszRasterColumnName,
752 : poWKTRasterDS->pszSchemaName,
753 : poWKTRasterDS->pszTableName,
754 : poWKTRasterDS->pszRasterColumnName, dfProjLowerLeftX,
755 : dfProjLowerLeftY,dfProjUpperRightX, dfProjUpperRightY,
756 0 : poWKTRasterDS->nSrid);
757 :
758 : } /**
759 : * Table hasn't a GIST index. Normal searching
760 : */
761 : else {
762 : osCommand.Printf(
763 : "SELECT rid, %s FROM %s.%s WHERE _ST_Contains(%s, ST_SetSRID(\
764 : ST_MakeBox2D(ST_Point(%f, %f), ST_Point(%f, %f)), %d))",
765 : poWKTRasterDS->pszRasterColumnName,
766 : poWKTRasterDS->pszSchemaName,
767 : poWKTRasterDS->pszTableName,
768 : poWKTRasterDS->pszRasterColumnName, dfProjLowerLeftX,
769 : dfProjLowerLeftY, dfProjUpperRightX, dfProjUpperRightY,
770 0 : poWKTRasterDS->nSrid);
771 : }
772 : }
773 :
774 :
775 : //printf("query: %s\n", szCommand);
776 :
777 0 : hPGresult = PQexec(poWKTRasterDS->hPGconn, osCommand.c_str());
778 :
779 0 : if (hPGresult == NULL ||
780 : PQresultStatus(hPGresult) != PGRES_TUPLES_OK ||
781 : PQntuples(hPGresult) < 0)
782 : {
783 0 : if (hPGresult)
784 0 : PQclear(hPGresult);
785 : CPLError(CE_Failure, CPLE_AppDefined,
786 : "Sorry, couldn't fetch block information from database: %s",
787 0 : PQerrorMessage(poWKTRasterDS->hPGconn));
788 0 : return CE_Failure;
789 : }
790 :
791 :
792 0 : nTuples = PQntuples(hPGresult);
793 :
794 : /*****************************************************************
795 : * No blocks found. Fill the buffer with nodata value
796 : *****************************************************************/
797 0 : if (nTuples == 0) {
798 0 : NullBlock(pImage);
799 0 : return CE_None;
800 : }
801 :
802 :
803 : /******************************************************************
804 : * One block found. Regular blocking arrangements, no overlaps
805 : ******************************************************************/
806 0 : else if (nTuples == 1) {
807 0 : int nRid = atoi(PQgetvalue(hPGresult, 0, 0));
808 :
809 : // Get HEXWKB representation of the block
810 0 : pszHexWkb = CPLStrdup(PQgetvalue(hPGresult, 0, 1));
811 :
812 0 : PQclear(hPGresult);
813 :
814 : // Raster hex must have an even number of characters
815 0 : if (pszHexWkb == NULL || strlen(pszHexWkb) % 2) {
816 : CPLError(CE_Failure, CPLE_AppDefined,
817 : "The HEXWKB data fetch from database must have an even number \
818 0 : of characters");
819 0 : return CE_Failure;
820 : }
821 :
822 :
823 : // Create a wrapper object
824 0 : poWKTRasterWrapper = new WKTRasterWrapper();
825 0 : if (poWKTRasterWrapper->Initialize(pszHexWkb) == FALSE) {
826 0 : CPLFree(pszHexWkb);
827 0 : delete poWKTRasterWrapper;
828 0 : return CE_Failure;
829 : }
830 :
831 :
832 : // We won't need this
833 0 : CPLFree(pszHexWkb);
834 :
835 : // Create raster band wrapper
836 0 : poWKTRasterBandWrapper = poWKTRasterWrapper->GetBand((GUInt16)nBand);
837 0 : if (poWKTRasterBandWrapper == NULL) {
838 0 : CPLError(CE_Failure, CPLE_ObjectNull,"Couldn't fetch band data");
839 0 : delete poWKTRasterWrapper;
840 0 : return CE_Failure;
841 : }
842 :
843 : // Get raster data
844 0 : pbyRasterData = poWKTRasterBandWrapper->GetData();
845 :
846 :
847 :
848 : //printf("IReadBlock with offset: %d, %d\n", nBlockXOff, nBlockYOff);
849 :
850 : /**********************************************************
851 : * Check if the raster is offline
852 : **********************************************************/
853 0 : if (poWKTRasterBandWrapper->bIsOffline == TRUE) {
854 : // The raster data in this case is a path to the raster file
855 0 : int nBandToRead = poWKTRasterBandWrapper->nOutDbBandNumber;
856 :
857 :
858 : // Open dataset, if needed
859 0 : if (poWKTRasterDS->poOutdbRasterDS == NULL) {
860 : poWKTRasterDS->poOutdbRasterDS = (GDALDataset *)
861 0 : GDALOpen((char *)pbyRasterData, GA_ReadOnly);
862 : }
863 :
864 0 : if (poWKTRasterDS->poOutdbRasterDS != NULL)
865 :
866 : // Read data from band
867 : /**
868 : * NOT SO SIMPLE!!!!
869 : * The outdb file may don't have the same block structure...
870 : */
871 :
872 : /*
873 : poWKTRasterDS->poOutdbRasterDS->GetRasterBand(nBandToRead)->ReadBlock(
874 : nBlockXOff, nBlockYOff, pImage);
875 : */
876 : ;
877 : else {
878 : CPLError(CE_Failure, CPLE_ObjectNull,
879 0 : "Couldn't read band data from out-db raster");
880 0 : delete poWKTRasterWrapper;
881 0 : return CE_Failure;
882 : }
883 :
884 : }
885 :
886 :
887 : else {
888 : /**
889 : * Copy the data buffer into pImage.
890 : * nBlockXSize * nBlockYSize should be equal to nDataSize/2
891 : */
892 : memcpy(pImage, pbyRasterData,
893 : (nNaturalXBlockSize * nNaturalYBlockSize) *
894 0 : nPixelSize * sizeof (GByte));
895 : }
896 :
897 :
898 :
899 : // Free resources and exit
900 0 : delete poWKTRasterWrapper;
901 :
902 0 : return CE_None;
903 :
904 : }// end if ntuples == 1
905 :
906 : /*********************************************************************
907 : * More than one block found. Non regular blocking arrangements
908 : *********************************************************************/
909 : else {
910 : // Only regular_block supported, just now
911 0 : PQclear(hPGresult);
912 : CPLError(CE_Failure, CPLE_NotSupported,
913 : "Sorry, but the raster presents block overlapping. This feature \
914 0 : is under development");
915 :
916 0 : return CE_Failure;
917 0 : }
918 : }
919 :
920 : /**
921 : * Set the block data to the null value if it is set, or zero if there is
922 : * no null data value.
923 : * Parameters:
924 : * - void *: the block data
925 : * Returns: nothing
926 : */
927 0 : void WKTRasterRasterBand::NullBlock(void *pData) {
928 0 : VALIDATE_POINTER0(pData, "NullBlock");
929 :
930 0 : int nNaturalBlockXSize = 0;
931 0 : int nNaturalBlockYSize = 0;
932 0 : GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
933 :
934 0 : int nWords = nNaturalBlockXSize * nNaturalBlockYSize;
935 0 : int nChunkSize = MAX(1, GDALGetDataTypeSize(eDataType) / 8);
936 :
937 : int bNoDataSet;
938 0 : double dfNoData = GetNoDataValue(&bNoDataSet);
939 0 : if (!bNoDataSet) {
940 0 : memset(pData, 0, nWords * nChunkSize);
941 : } else {
942 0 : int i = 0;
943 0 : for (i = 0; i < nWords; i += nChunkSize)
944 0 : memcpy((GByte *) pData + i, &dfNoData, nChunkSize);
945 : }
946 : }
947 :
948 : /**
949 : * Set the no data value for this band.
950 : * Parameters:
951 : * - double: The nodata value
952 : * Returns:
953 : * - CE_None.
954 : */
955 0 : CPLErr WKTRasterRasterBand::SetNoDataValue(double dfNewValue) {
956 0 : dfNoDataValue = dfNewValue;
957 :
958 0 : return CE_None;
959 : }
960 :
961 : /**
962 : * Fetch the no data value for this band.
963 : * Parameters:
964 : * - int *: pointer to a boolean to use to indicate if a value is actually
965 : * associated with this layer. May be NULL (default).
966 : * Returns:
967 : * - double: the nodata value for this band.
968 : */
969 0 : double WKTRasterRasterBand::GetNoDataValue(int *pbSuccess) {
970 0 : if (pbSuccess != NULL)
971 0 : *pbSuccess = TRUE;
972 :
973 0 : return dfNoDataValue;
974 : }
975 :
976 : /**
977 : * Returns the number of overview layers available.
978 : * Parameters: none
979 : * Returns:
980 : * int: the number of overviews layers available
981 : */
982 0 : int WKTRasterRasterBand::GetOverviewCount() {
983 0 : WKTRasterDataset * poWKTRasterDS = (WKTRasterDataset *) poDS;
984 :
985 : return (poWKTRasterDS->nOverviews > 0) ?
986 : poWKTRasterDS->nOverviews :
987 0 : GDALRasterBand::GetOverviewCount();
988 : }
989 :
990 : /**
991 : * Fetch overview raster band object.
992 : * Parameters:
993 : * - int: overview index between 0 and GetOverviewCount()-1
994 : * Returns:
995 : * - GDALRasterBand *: overview GDALRasterBand.
996 : */
997 0 : GDALRasterBand * WKTRasterRasterBand::GetOverview(int nOverview) {
998 0 : WKTRasterDataset * poWKTRasterDS = (WKTRasterDataset *) poDS;
999 :
1000 0 : if (poWKTRasterDS->nOverviews > 0) {
1001 0 : if (nOverview < 0 || nOverview >= poWKTRasterDS->nOverviews)
1002 0 : return NULL;
1003 : else
1004 : return
1005 0 : poWKTRasterDS->papoWKTRasterOv[nOverview]->GetRasterBand(nBand);
1006 : } else
1007 0 : return GDALRasterBand::GetOverview(nOverview);
1008 : }
1009 :
1010 : /**
1011 : * Get the natural block size for this band.
1012 : * Parameters:
1013 : * - int *: pointer to int to store the natural X block size
1014 : * - int *: pointer to int to store the natural Y block size
1015 : * Returns: nothing
1016 : */
1017 0 : void WKTRasterRasterBand::GetBlockSize(int * pnXSize, int *pnYSize)
1018 : {
1019 0 : if (nBlockXSize == 0 || nBlockYSize == 0) {
1020 : CPLError(CE_Failure, CPLE_AppDefined,
1021 : "This WKT Raster band has non regular blocking arrangement. \
1022 0 : This feature is under development");
1023 :
1024 0 : if (pnXSize != NULL)
1025 0 : *pnXSize = 0;
1026 0 : if (pnYSize != NULL)
1027 0 : *pnYSize = 0;
1028 :
1029 : }
1030 : else {
1031 0 : GDALRasterBand::GetBlockSize(pnXSize, pnYSize);
1032 : }
1033 0 : }
1034 :
1035 :
1036 :
1037 :
1038 :
|
