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 :
39 :
40 : /**
41 : * \brief Constructor.
42 : * Parameters:
43 : * - PostGISRasterDataset *: The Dataset this band belongs to
44 : * - int: the band number
45 : * - GDALDataType: The data type for this band
46 : * - double: The nodata value. Could be any kind of data (GByte, GUInt16,
47 : * GInt32...) but the variable has the bigger type.
48 : * - GBool: if the data type is signed byte or not. If yes, the SIGNEDBYTE
49 : * metadata value is added to IMAGE_STRUCTURE domain
50 : * - int: The bit depth, to add NBITS as metadata value in IMAGE_STRUCTURE
51 : * domain.
52 : */
53 0 : PostGISRasterRasterBand::PostGISRasterRasterBand(PostGISRasterDataset *poDS,
54 : int nBand, GDALDataType hDataType, GBool bHasNoDataValue, double dfNodata, GBool bSignedByte,
55 : int nBitDepth, int nFactor, GBool bIsOffline, char * pszSchema,
56 0 : char * pszTable, char * pszColumn)
57 : {
58 :
59 : /* Basic properties */
60 0 : this->poDS = poDS;
61 0 : this->nBand = nBand;
62 0 : this->bIsOffline = bIsOffline;
63 0 : this->pszSchema = (pszSchema) ? pszSchema : CPLStrdup(poDS->pszSchema);
64 0 : this->pszTable = (pszTable) ? pszTable: CPLStrdup(poDS->pszTable);
65 0 : this->pszColumn = (pszColumn) ? pszColumn : CPLStrdup(poDS->pszColumn);
66 0 : this->pszWhere = CPLStrdup(poDS->pszWhere);
67 :
68 0 : eAccess = poDS->GetAccess();
69 0 : eDataType = hDataType;
70 0 : this->bHasNoDataValue = bHasNoDataValue;
71 0 : dfNoDataValue = dfNodata;
72 :
73 :
74 : /*****************************************************
75 : * Check block size issue
76 : *****************************************************/
77 0 : nBlockXSize = poDS->nBlockXSize;
78 0 : nBlockYSize = poDS->nBlockYSize;
79 :
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 : // Add pixeltype to image structure domain
87 0 : if (bSignedByte == true) {
88 0 : SetMetadataItem("PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
89 : }
90 :
91 : // Add NBITS to metadata only for sub-byte types
92 0 : if (nBitDepth < 8)
93 : SetMetadataItem("NBITS", CPLString().Printf( "%d", nBitDepth ),
94 0 : "IMAGE_STRUCTURE" );
95 :
96 :
97 0 : nOverviewFactor = nFactor;
98 :
99 : /**********************************************************
100 : * Check overviews. Query RASTER_OVERVIEWS table for
101 : * existing overviews, only in case we are on level 0
102 : * TODO: can we do this without querying RASTER_OVERVIEWS?
103 : * How do we know the number of overviews? Is an inphinite
104 : * loop...
105 : **********************************************************/
106 0 : if (nOverviewFactor == 0) {
107 :
108 0 : CPLString osCommand;
109 0 : PGresult * poResult = NULL;
110 0 : int i = 0;
111 0 : int nFetchOvFactor = 0;
112 0 : char * pszOvSchema = NULL;
113 0 : char * pszOvTable = NULL;
114 0 : char * pszOvColumn = NULL;
115 :
116 0 : nRasterXSize = poDS->GetRasterXSize();
117 0 : nRasterYSize = poDS->GetRasterYSize();
118 :
119 : osCommand.Printf("select o_table_name, overview_factor, o_raster_column, "
120 : "o_table_schema from raster_overviews where r_table_schema = "
121 : "'%s' and r_table_name = '%s' and r_raster_column = '%s'",
122 0 : poDS->pszSchema, poDS->pszTable, poDS->pszColumn);
123 :
124 0 : poResult = PQexec(poDS->poConn, osCommand.c_str());
125 0 : if (poResult != NULL && PQresultStatus(poResult) == PGRES_TUPLES_OK &&
126 : PQntuples(poResult) > 0) {
127 :
128 : /* Create overviews */
129 0 : nOverviewCount = PQntuples(poResult);
130 : papoOverviews = (PostGISRasterRasterBand **)VSICalloc(nOverviewCount,
131 0 : sizeof(PostGISRasterRasterBand *));
132 0 : if (papoOverviews == NULL) {
133 : CPLError(CE_Warning, CPLE_OutOfMemory, "Couldn't create "
134 0 : "overviews for band %d\n", nBand);
135 0 : PQclear(poResult);
136 : return;
137 : }
138 :
139 0 : for(i = 0; i < nOverviewCount; i++) {
140 :
141 0 : nFetchOvFactor = atoi(PQgetvalue(poResult, i, 1));
142 0 : pszOvSchema = CPLStrdup(PQgetvalue(poResult, i, 3));
143 0 : pszOvTable = CPLStrdup(PQgetvalue(poResult, i, 0));
144 0 : pszOvColumn = CPLStrdup(PQgetvalue(poResult, i, 2));
145 :
146 : /**
147 : * NOTE: Overview bands are not considered to be a part of a
148 : * dataset, but we use the same dataset for all the overview
149 : * bands just for simplification (we'll need to access the table
150 : * and schema names). But in method GetDataset, NULL is return
151 : * if we're talking about an overview band
152 : */
153 0 : papoOverviews[i] = new PostGISRasterRasterBand(poDS, nBand,
154 : hDataType, bHasNoDataValue, dfNodata, bSignedByte, nBitDepth,
155 0 : nFetchOvFactor, bIsOffline, pszOvSchema, pszOvTable, pszOvColumn);
156 :
157 : }
158 :
159 0 : PQclear(poResult);
160 :
161 : }
162 :
163 : else {
164 :
165 0 : nOverviewCount = 0;
166 0 : papoOverviews = NULL;
167 0 : if (poResult)
168 0 : PQclear(poResult);
169 0 : }
170 : }
171 :
172 : /************************************
173 : * We are in an overview level. Set
174 : * raster size to its value
175 : ************************************/
176 : else {
177 :
178 : /*
179 : * No overviews inside an overview (all the overviews are from original
180 : * band
181 : */
182 0 : nOverviewCount = 0;
183 0 : papoOverviews = NULL;
184 :
185 0 : nRasterXSize = (int) floor((double)poDS->GetRasterXSize() / nOverviewFactor);
186 0 : nRasterYSize = (int) floor((double)poDS->GetRasterYSize() / nOverviewFactor);
187 : }
188 :
189 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand constructor: Band "
190 0 : "created (srid = %d)", poDS->nSrid);
191 0 : }
192 :
193 : /***********************************************
194 : * \brief: Band destructor
195 : ***********************************************/
196 0 : PostGISRasterRasterBand::~PostGISRasterRasterBand()
197 : {
198 : int i;
199 :
200 0 : if (pszSchema)
201 0 : CPLFree(pszSchema);
202 0 : if (pszTable)
203 0 : CPLFree(pszTable);
204 0 : if (pszColumn)
205 0 : CPLFree(pszColumn);
206 0 : if (pszWhere)
207 0 : CPLFree(pszWhere);
208 :
209 0 : if (papoOverviews) {
210 0 : for(i = 0; i < nOverviewCount; i++)
211 0 : delete papoOverviews[i];
212 :
213 0 : CPLFree(papoOverviews);
214 : }
215 0 : }
216 :
217 :
218 : /**
219 : * \brief Set the block data to the null value if it is set, or zero if there is
220 : * no null data value.
221 : * Parameters:
222 : * - void *: the block data
223 : * Returns: nothing
224 : */
225 0 : void PostGISRasterRasterBand::NullBlock(void *pData)
226 : {
227 0 : VALIDATE_POINTER0(pData, "NullBlock");
228 :
229 0 : int nNaturalBlockXSize = 0;
230 0 : int nNaturalBlockYSize = 0;
231 0 : GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
232 :
233 0 : int nWords = nNaturalBlockXSize * nNaturalBlockYSize;
234 0 : int nChunkSize = MAX(1, GDALGetDataTypeSize(eDataType) / 8);
235 :
236 : int bNoDataSet;
237 0 : double dfNoData = GetNoDataValue(&bNoDataSet);
238 0 : if (!bNoDataSet)
239 : {
240 0 : memset(pData, 0, nWords * nChunkSize);
241 : }
242 : else
243 : {
244 0 : int i = 0;
245 0 : for (i = 0; i < nWords; i += nChunkSize)
246 0 : memcpy((GByte *) pData + i, &dfNoData, nChunkSize);
247 : }
248 : }
249 :
250 : /**
251 : * \brief Set the no data value for this band.
252 : * Parameters:
253 : * - double: The nodata value
254 : * Returns:
255 : * - CE_None.
256 : */
257 0 : CPLErr PostGISRasterRasterBand::SetNoDataValue(double dfNewValue) {
258 0 : dfNoDataValue = dfNewValue;
259 :
260 0 : return CE_None;
261 : }
262 :
263 : /**
264 : * \brief Fetch the no data value for this band.
265 : * Parameters:
266 : * - int *: pointer to a boolean to use to indicate if a value is actually
267 : * associated with this layer. May be NULL (default).
268 : * Returns:
269 : * - double: the nodata value for this band.
270 : */
271 0 : double PostGISRasterRasterBand::GetNoDataValue(int *pbSuccess) {
272 0 : if (pbSuccess != NULL)
273 0 : *pbSuccess = (int) bHasNoDataValue;
274 :
275 0 : return dfNoDataValue;
276 : }
277 :
278 :
279 : /**
280 : * \brief Get the natural block size for this band.
281 : * Parameters:
282 : * - int *: pointer to int to store the natural X block size
283 : * - int *: pointer to int to store the natural Y block size
284 : * Returns: nothing
285 : */
286 0 : void PostGISRasterRasterBand::GetBlockSize(int * pnXSize, int *pnYSize)
287 : {
288 0 : if (nBlockXSize == 0 || nBlockYSize == 0) {
289 : CPLError(CE_Failure, CPLE_AppDefined,
290 : "This PostGIS Raster band has non regular blocking arrangement. \
291 0 : This feature is under development");
292 :
293 0 : if (pnXSize != NULL)
294 0 : *pnXSize = 0;
295 0 : if (pnYSize != NULL)
296 0 : *pnYSize = 0;
297 :
298 : }
299 : else {
300 0 : GDALRasterBand::GetBlockSize(pnXSize, pnYSize);
301 : }
302 0 : }
303 :
304 :
305 : /*****************************************************
306 : * \brief Fetch the band number
307 : *****************************************************/
308 0 : int PostGISRasterRasterBand::GetBand()
309 : {
310 0 : return (nOverviewFactor) ? 0 : nBand;
311 : }
312 :
313 : /*****************************************************
314 : * \brief Fetch the owning dataset handle
315 : *****************************************************/
316 0 : GDALDataset* PostGISRasterRasterBand::GetDataset()
317 : {
318 0 : return (nOverviewFactor) ? NULL : poDS;
319 : }
320 :
321 : /****************************************************
322 : * \brief Check for arbitrary overviews
323 : * The datastore can compute arbitrary overviews
324 : * efficiently, because the overviews are tables,
325 : * like the original raster. The effort is the same.
326 : ****************************************************/
327 0 : int PostGISRasterRasterBand::HasArbitraryOverviews()
328 : {
329 0 : return (nOverviewFactor) ? false : true;
330 : }
331 :
332 : /***************************************************
333 : * \brief Return the number of overview layers available
334 : ***************************************************/
335 0 : int PostGISRasterRasterBand::GetOverviewCount()
336 : {
337 : return (nOverviewFactor) ?
338 : 0 :
339 0 : nOverviewCount;
340 : }
341 :
342 : /**********************************************************
343 : * \brief Fetch overview raster band object
344 : **********************************************************/
345 0 : GDALRasterBand * PostGISRasterRasterBand::GetOverview(int i)
346 : {
347 0 : return (i >= 0 && i < GetOverviewCount()) ?
348 0 : (GDALRasterBand *)papoOverviews[i] : GDALRasterBand::GetOverview(i);
349 : }
350 :
351 : /*****************************************************
352 : * \brief Read a natural block of raster band data
353 : *****************************************************/
354 0 : CPLErr PostGISRasterRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void*
355 : pImage)
356 : {
357 0 : PostGISRasterDataset * poPostGISRasterDS = (PostGISRasterDataset*)poDS;
358 0 : CPLString osCommand;
359 0 : PGresult* poResult = NULL;
360 : double adfTransform[6];
361 0 : int nTuples = 0;
362 0 : int nNaturalBlockXSize = 0;
363 0 : int nNaturalBlockYSize = 0;
364 0 : int nPixelSize = 0;
365 : int nPixelXInitPosition;
366 : int nPixelYInitPosition;
367 : int nPixelXEndPosition;
368 : int nPixelYEndPosition;
369 0 : double dfProjXInit = 0.0;
370 0 : double dfProjYInit = 0.0;
371 0 : double dfProjXEnd = 0.0;
372 0 : double dfProjYEnd = 0.0;
373 0 : double dfProjLowerLeftX = 0.0;
374 0 : double dfProjLowerLeftY = 0.0;
375 0 : double dfProjUpperRightX = 0.0;
376 0 : double dfProjUpperRightY = 0.0;
377 0 : GByte* pbyData = NULL;
378 0 : char* pbyDataToRead = NULL;
379 0 : int nWKBLength = 0;
380 0 : int nExpectedDataSize = 0;
381 :
382 : /* Get pixel and block size */
383 0 : nPixelSize = MAX(1,GDALGetDataTypeSize(eDataType)/8);
384 0 : GetBlockSize(&nNaturalBlockXSize, &nNaturalBlockYSize);
385 :
386 : /* Get pixel,line coordinates of the block */
387 : /**
388 : * TODO: What if the georaster is rotated? Following the gdal_translate
389 : * source code, you can't use -projwin option with rotated rasters
390 : **/
391 0 : nPixelXInitPosition = nNaturalBlockXSize * nBlockXOff;
392 0 : nPixelYInitPosition = nNaturalBlockYSize * nBlockYOff;
393 0 : nPixelXEndPosition = nPixelXInitPosition + nNaturalBlockXSize;
394 0 : nPixelYEndPosition = nPixelYInitPosition + nNaturalBlockYSize;
395 :
396 0 : poPostGISRasterDS->GetGeoTransform(adfTransform);
397 :
398 : /* Pixel size correction, in case of overviews */
399 0 : if (nOverviewFactor) {
400 0 : adfTransform[1] *= nOverviewFactor;
401 0 : adfTransform[5] *= nOverviewFactor;
402 : }
403 :
404 : /* Calculate the "query box" */
405 0 : dfProjXInit = adfTransform[0] +
406 0 : nPixelXInitPosition * adfTransform[1] +
407 0 : nPixelYInitPosition * adfTransform[2];
408 :
409 0 : dfProjYInit = adfTransform[3] +
410 0 : nPixelXInitPosition * adfTransform[4] +
411 0 : nPixelYInitPosition * adfTransform[5];
412 :
413 0 : dfProjXEnd = adfTransform[0] +
414 0 : nPixelXEndPosition * adfTransform[1] +
415 0 : nPixelYEndPosition * adfTransform[2];
416 :
417 0 : dfProjYEnd = adfTransform[3] +
418 0 : nPixelXEndPosition * adfTransform[4] +
419 0 : nPixelYEndPosition * adfTransform[5];
420 :
421 : /*************************************************************************
422 : * Now we have the block coordinates transformed into coordinates of the
423 : * raster reference system. This coordinates are from:
424 : * - Upper left corner
425 : * - Lower right corner
426 : * But for ST_MakeBox2D, we'll need block's coordinates of:
427 : * - Lower left corner
428 : * - Upper right corner
429 : *************************************************************************/
430 0 : dfProjLowerLeftX = dfProjXInit;
431 0 : dfProjLowerLeftY = dfProjYEnd;
432 :
433 0 : dfProjUpperRightX = dfProjXEnd;
434 0 : dfProjUpperRightY = dfProjYInit;
435 :
436 : /**
437 : * Return all raster objects from our raster table that fall reside or
438 : * partly reside in a coordinate bounding box.
439 : * NOTE: Is it necessary to use a BB optimization like:
440 : * select st_makebox2d(...) && geom and (rest of the query)...?
441 : **/
442 0 : if (poPostGISRasterDS->pszWhere != NULL)
443 : {
444 : osCommand.Printf("select rid, %s from %s.%s where %s ~ "
445 : "st_setsrid(st_makebox2d(st_point(%f, %f), st_point(%f,"
446 : "%f)),%d) and %s", pszColumn, pszSchema, pszTable, pszColumn,
447 : dfProjLowerLeftX, dfProjLowerLeftY, dfProjUpperRightX,
448 0 : dfProjUpperRightY, poPostGISRasterDS->nSrid, pszWhere);
449 : }
450 :
451 : else
452 : {
453 : osCommand.Printf("select rid, %s from %s.%s where %s ~ "
454 : "st_setsrid(st_makebox2d(st_point(%f, %f), st_point(%f,"
455 : "%f)),%d)", pszColumn, pszSchema, pszTable, pszColumn,
456 : dfProjLowerLeftX, dfProjLowerLeftY, dfProjUpperRightX,
457 0 : dfProjUpperRightY, poPostGISRasterDS->nSrid);
458 : }
459 :
460 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IReadBlock: "
461 0 : "The query = %s", osCommand.c_str());
462 :
463 0 : poResult = PQexec(poPostGISRasterDS->poConn, osCommand.c_str());
464 0 : if (poResult == NULL || PQresultStatus(poResult) != PGRES_TUPLES_OK ||
465 : PQntuples(poResult) <= 0)
466 : {
467 0 : if (poResult)
468 0 : PQclear(poResult);
469 :
470 : /* TODO: Raise an error and exit? */
471 : CPLDebug("PostGIS_Raster", "PostGISRasterRasterBand::IReadBlock: "
472 0 : "The block (%d, %d) is empty", nBlockXOff, nBlockYOff);
473 0 : NullBlock(pImage);
474 0 : return CE_None;
475 : }
476 :
477 0 : nTuples = PQntuples(poResult);
478 :
479 : /* No overlapping */
480 0 : if (nTuples == 1)
481 : {
482 : /**
483 : * Out db band.
484 : * TODO: Manage this situation
485 : **/
486 0 : if (bIsOffline)
487 : {
488 : CPLError(CE_Failure, CPLE_AppDefined, "This raster has outdb storage."
489 0 : "This feature isn\'t still available");
490 0 : return CE_Failure;
491 : }
492 :
493 0 : int nRid = atoi(PQgetvalue(poResult, 0, 0));
494 :
495 : /* Only data size, without payload */
496 : nExpectedDataSize = nNaturalBlockXSize * nNaturalBlockYSize *
497 0 : nPixelSize;
498 0 : pbyData = CPLHexToBinary(PQgetvalue(poResult, 0, 1), &nWKBLength);
499 : pbyDataToRead = (char*)GET_BAND_DATA(pbyData,nBand, nPixelSize,
500 0 : nExpectedDataSize);
501 :
502 0 : memcpy(pImage, pbyDataToRead, nExpectedDataSize * sizeof(char));
503 :
504 : CPLDebug("PostGIS_Raster", "IReadBlock: Copied %d bytes from block "
505 : "(%d, %d) (rid = %d) to %p", nExpectedDataSize, nBlockXOff,
506 0 : nBlockYOff, nRid, pImage);
507 :
508 0 : CPLFree(pbyData);
509 0 : PQclear(poResult);
510 :
511 0 : return CE_None;
512 : }
513 :
514 : /** Overlapping raster data.
515 : * TODO: Manage this situation. Suggestion: open several datasets, because
516 : * you can't manage overlapping raster data with only one dataset
517 : **/
518 : else
519 : {
520 : CPLError(CE_Failure, CPLE_AppDefined,
521 : "Overlapping raster data. Feature under development, not "
522 0 : "available yet");
523 0 : if (poResult)
524 0 : PQclear(poResult);
525 :
526 0 : return CE_Failure;
527 :
528 0 : }
529 : }
530 :
531 :
|
