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