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 :
|
