Line data Source code
1 : /***********************************************************************
2 : * File : postgisrasterrasterband.cpp
3 : * Project: PostGIS Raster driver
4 : * Purpose: GDAL RasterBand implementation for PostGIS Raster driver
5 : * Author: Jorge Arevalo, jorge.arevalo@deimos-space.com
6 : * jorgearevalo@libregis.org
7 : *
8 : * Author: David Zwarg, dzwarg@azavea.com
9 : *
10 : *
11 : ***********************************************************************
12 : * Copyright (c) 2009 - 2013, Jorge Arevalo, David Zwarg
13 : * Copyright (c) 2013-2018, Even Rouault <even.rouault at spatialys.com>
14 : *
15 : * SPDX-License-Identifier: MIT
16 : **********************************************************************/
17 : #include "postgisraster.h"
18 :
19 : #include <algorithm>
20 : #include <limits>
21 : #include <cmath>
22 :
23 : /**
24 : * \brief Constructor.
25 : *
26 : * nBand it is just necessary for overview band creation
27 : */
28 0 : PostGISRasterRasterBand::PostGISRasterRasterBand(PostGISRasterDataset *poDSIn,
29 : int nBandIn,
30 : GDALDataType eDataTypeIn,
31 : GBool bNoDataValueSetIn,
32 0 : double dfNodata)
33 0 : : VRTSourcedRasterBand(poDSIn, nBandIn), pszSchema(poDSIn->pszSchema),
34 0 : pszTable(poDSIn->pszTable), pszColumn(poDSIn->pszColumn)
35 : {
36 : /* Basic properties */
37 0 : poDS = poDSIn;
38 0 : nBand = nBandIn;
39 :
40 0 : eDataType = eDataTypeIn;
41 0 : m_bNoDataValueSet = CPL_TO_BOOL(bNoDataValueSetIn);
42 0 : m_dfNoDataValue = dfNodata;
43 :
44 0 : nRasterXSize = poDS->GetRasterXSize();
45 0 : nRasterYSize = poDS->GetRasterYSize();
46 :
47 : /*******************************************************************
48 : * Finally, set the block size. We apply the same logic than in VRT
49 : * driver.
50 : *
51 : * We limit the size of a block with MAX_BLOCK_SIZE here to prevent
52 : * arrangements of just one big tile.
53 : *
54 : * This value is just used in case we only have 1 tile in the
55 : * table. Otherwise, the reading operations are performed by the
56 : * sources, not the PostGISRasterBand object itself.
57 : ******************************************************************/
58 0 : nBlockXSize = atoi(CPLGetConfigOption(
59 : "PR_BLOCKXSIZE",
60 0 : CPLSPrintf("%d", std::min(MAX_BLOCK_SIZE, this->nRasterXSize))));
61 0 : nBlockYSize = atoi(CPLGetConfigOption(
62 : "PR_BLOCKYSIZE",
63 0 : CPLSPrintf("%d", std::min(MAX_BLOCK_SIZE, this->nRasterYSize))));
64 :
65 : #ifdef DEBUG_VERBOSE
66 : CPLDebug("PostGIS_Raster",
67 : "PostGISRasterRasterBand constructor: Band size: (%d X %d)",
68 : nRasterXSize, nRasterYSize);
69 :
70 : CPLDebug("PostGIS_Raster",
71 : "PostGISRasterRasterBand::Constructor: "
72 : "Block size (%dx%d)",
73 : this->nBlockXSize, this->nBlockYSize);
74 : #endif
75 0 : }
76 :
77 : /***********************************************
78 : * \brief: Band destructor
79 : ***********************************************/
80 0 : PostGISRasterRasterBand::~PostGISRasterRasterBand()
81 : {
82 0 : }
83 :
84 : /********************************************************
85 : * \brief Query statistics for this band
86 : ********************************************************/
87 0 : bool PostGISRasterRasterBand::QueryStats()
88 : {
89 0 : m_dfStatsCount = std::numeric_limits<double>::quiet_NaN();
90 0 : m_dfStatsSum = std::numeric_limits<double>::quiet_NaN();
91 0 : m_dfStatsMean = std::numeric_limits<double>::quiet_NaN();
92 0 : m_dfStatsStdDev = std::numeric_limits<double>::quiet_NaN();
93 0 : m_dfStatsMin = std::numeric_limits<double>::quiet_NaN();
94 0 : m_dfStatsMax = std::numeric_limits<double>::quiet_NaN();
95 :
96 0 : m_bStatsFetched = false;
97 :
98 0 : PostGISRasterDataset *poRDS = cpl::down_cast<PostGISRasterDataset *>(poDS);
99 0 : const CPLString osSchemaI(CPLQuotedSQLIdentifier(pszSchema));
100 0 : const CPLString osTableI(CPLQuotedSQLIdentifier(pszTable));
101 0 : const CPLString osColumnI(CPLQuotedSQLIdentifier(pszColumn));
102 0 : std::string osCommand = "SELECT ST_SummaryStatsAgg(";
103 0 : osCommand += osColumnI;
104 0 : osCommand += ", ";
105 0 : osCommand += std::to_string(nBand);
106 0 : osCommand += ", TRUE) FROM ";
107 0 : osCommand += osSchemaI;
108 0 : osCommand += ".";
109 0 : osCommand += osTableI;
110 0 : PGresult *poResult = PQexec(poRDS->poConn, osCommand.c_str());
111 0 : if (PQresultStatus(poResult) != PGRES_TUPLES_OK)
112 : {
113 0 : CPLError(CE_Failure, CPLE_AppDefined,
114 : "PostGISRasterRasterBand::queryStats(): "
115 : "Error executing query: %s",
116 0 : PQerrorMessage(poRDS->poConn));
117 0 : PQclear(poResult);
118 0 : return false;
119 : }
120 : else
121 : {
122 0 : if (PQntuples(poResult) > 0 && PQgetisnull(poResult, 0, 0) == 0)
123 : {
124 0 : std::string osStats = PQgetvalue(poResult, 0, 0);
125 0 : if (osStats.size() <= 2)
126 : {
127 0 : CPLError(CE_Failure, CPLE_AppDefined,
128 : "PostGISRasterRasterBand::queryStats(): "
129 : "Unexpected result from ST_SummaryStatsAgg: %s",
130 : osStats.c_str());
131 0 : PQclear(poResult);
132 0 : return false;
133 : }
134 : // Remove trailing and ending parenthesis
135 0 : osStats.erase(0, 1);
136 0 : osStats.erase(osStats.size() - 1);
137 : const CPLStringList aosTokens(
138 0 : CSLTokenizeString2(osStats.c_str(), ",", 0), false);
139 : // count, sum, mean, stddev, min, max
140 0 : if (CSLCount(aosTokens) == 6)
141 : {
142 0 : m_dfStatsCount = CPLAtof(aosTokens[0]);
143 0 : m_dfStatsSum = CPLAtof(aosTokens[1]);
144 0 : m_dfStatsMean = CPLAtof(aosTokens[2]);
145 0 : m_dfStatsStdDev = CPLAtof(aosTokens[3]);
146 0 : m_dfStatsMin = CPLAtof(aosTokens[4]);
147 0 : m_dfStatsMax = CPLAtof(aosTokens[5]);
148 : }
149 : else
150 : {
151 0 : CPLError(CE_Failure, CPLE_AppDefined,
152 : "PostGISRasterRasterBand::queryStats(): "
153 : "Unexpected number of tokens in ST_SummaryStatsAgg "
154 : "result: %d",
155 : CSLCount(aosTokens));
156 : }
157 : }
158 0 : PQclear(poResult);
159 0 : m_bStatsFetched = !std::isnan(m_dfStatsCount);
160 0 : return m_bStatsFetched;
161 : }
162 : }
163 :
164 : /********************************************************
165 : * \brief Check if statistics have been fetched and are valid
166 : *********************************************************/
167 0 : bool PostGISRasterRasterBand::StatsFetchedAndValid() const
168 : {
169 0 : return m_bStatsFetched && !std::isnan(m_dfStatsCount) &&
170 0 : !std::isnan(m_dfStatsSum) && !std::isnan(m_dfStatsMean) &&
171 0 : !std::isnan(m_dfStatsStdDev) && !std::isnan(m_dfStatsMin) &&
172 0 : !std::isnan(m_dfStatsMax) && m_dfStatsCount > 0;
173 : }
174 :
175 : /********************************************************
176 : * \brief Set nodata value to a buffer
177 : ********************************************************/
178 0 : void PostGISRasterRasterBand::NullBuffer(void *pData, int nBufXSize,
179 : int nBufYSize, GDALDataType eBufType,
180 : int nPixelSpace, int nLineSpace)
181 : {
182 : int j;
183 0 : for (j = 0; j < nBufYSize; j++)
184 : {
185 0 : double dfVal = 0.0;
186 0 : if (m_bNoDataValueSet)
187 0 : dfVal = m_dfNoDataValue;
188 0 : GDALCopyWords(&dfVal, GDT_Float64, 0,
189 0 : static_cast<GByte *>(pData) + j * nLineSpace, eBufType,
190 : nPixelSpace, nBufXSize);
191 : }
192 0 : }
193 :
194 : /********************************************************
195 : * \brief SortTilesByPKID
196 : ********************************************************/
197 0 : static int SortTilesByPKID(const void *a, const void *b)
198 : {
199 0 : const PostGISRasterTileDataset *pa =
200 : *static_cast<const PostGISRasterTileDataset *const *>(a);
201 0 : const PostGISRasterTileDataset *pb =
202 : *static_cast<const PostGISRasterTileDataset *const *>(b);
203 0 : return strcmp(pa->GetPKID(), pb->GetPKID());
204 : }
205 :
206 : /**
207 : * Read/write a region of image data for this band.
208 : *
209 : * This method allows reading a region of a PostGISRasterBand into a buffer.
210 : * The write support is still under development
211 : *
212 : * The function fetches all the raster data that intersects with the region
213 : * provided, and store the data in the GDAL cache.
214 : *
215 : * It automatically takes care of data type translation if the data type
216 : * (eBufType) of the buffer is different than that of the
217 : * PostGISRasterRasterBand.
218 : *
219 : * The nPixelSpace and nLineSpace parameters allow reading into FROM various
220 : * organization of buffers.
221 : *
222 : * @param eRWFlag Either GF_Read to read a region of data (GF_Write, to write
223 : * a region of data, yet not supported)
224 : *
225 : * @param nXOff The pixel offset to the top left corner of the region of the
226 : * band to be accessed. This would be zero to start FROM the left side.
227 : *
228 : * @param nYOff The line offset to the top left corner of the region of the band
229 : * to be accessed. This would be zero to start FROM the top.
230 : *
231 : * @param nXSize The width of the region of the band to be accessed in pixels.
232 : *
233 : * @param nYSize The height of the region of the band to be accessed in lines.
234 : *
235 : * @param pData The buffer into which the data should be read, or FROM which it
236 : * should be written. This buffer must contain at least
237 : * nBufXSize * nBufYSize * nBandCount words of type eBufType. It is organized in
238 : * left to right,top to bottom pixel order. Spacing is controlled by the
239 : * nPixelSpace, and nLineSpace parameters.
240 : *
241 : * @param nBufXSize the width of the buffer image into which the desired region
242 : * is to be read, or FROM which it is to be written.
243 : *
244 : * @param nBufYSize the height of the buffer image into which the desired region
245 : * is to be read, or FROM which it is to be written.
246 : *
247 : * @param eBufType the type of the pixel values in the pData data buffer. The
248 : * pixel values will automatically be translated to/FROM the
249 : * PostGISRasterRasterBand data type as needed.
250 : *
251 : * @param nPixelSpace The byte offset FROM the start of one pixel value in pData
252 : * to the start of the next pixel value within a scanline. If defaulted (0) the
253 : * size of the datatype eBufType is used.
254 : *
255 : * @param nLineSpace The byte offset FROM the start of one scanline in pData to
256 : * the start of the next. If defaulted (0) the size of the datatype
257 : * eBufType * nBufXSize is used.
258 : *
259 : * @return CE_Failure if the access fails, otherwise CE_None.
260 : */
261 :
262 0 : CPLErr PostGISRasterRasterBand::IRasterIO(
263 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
264 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
265 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
266 : {
267 : /**
268 : * TODO: Write support not implemented yet
269 : **/
270 0 : if (eRWFlag == GF_Write)
271 : {
272 0 : ReportError(CE_Failure, CPLE_NotSupported,
273 : "Writing through PostGIS Raster band not supported yet");
274 :
275 0 : return CE_Failure;
276 : }
277 :
278 : /*******************************************************************
279 : * Do we have overviews that would be appropriate to satisfy this
280 : * request?
281 : ******************************************************************/
282 0 : if ((nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0)
283 : {
284 0 : if (OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
285 : nBufXSize, nBufYSize, eBufType, nPixelSpace,
286 0 : nLineSpace, psExtraArg) == CE_None)
287 :
288 0 : return CE_None;
289 : }
290 :
291 0 : PostGISRasterDataset *poRDS = cpl::down_cast<PostGISRasterDataset *>(poDS);
292 :
293 0 : int bSameWindowAsOtherBand =
294 0 : (nXOff == poRDS->nXOffPrev && nYOff == poRDS->nYOffPrev &&
295 0 : nXSize == poRDS->nXSizePrev && nYSize == poRDS->nYSizePrev);
296 0 : poRDS->nXOffPrev = nXOff;
297 0 : poRDS->nYOffPrev = nYOff;
298 0 : poRDS->nXSizePrev = nXSize;
299 0 : poRDS->nYSizePrev = nYSize;
300 :
301 : /* Logic to determine if bands are read in order 1, 2, ... N */
302 : /* If so, then use multi-band caching, otherwise do just single band caching
303 : */
304 0 : if (poRDS->bAssumeMultiBandReadPattern)
305 : {
306 0 : if (nBand != poRDS->nNextExpectedBand)
307 : {
308 0 : CPLDebug("PostGIS_Raster", "Disabling multi-band caching since "
309 : "band access pattern does not match");
310 0 : poRDS->bAssumeMultiBandReadPattern = false;
311 0 : poRDS->nNextExpectedBand = 1;
312 : }
313 : else
314 : {
315 0 : poRDS->nNextExpectedBand++;
316 0 : if (poRDS->nNextExpectedBand > poRDS->GetRasterCount())
317 0 : poRDS->nNextExpectedBand = 1;
318 : }
319 : }
320 : else
321 : {
322 0 : if (nBand == poRDS->nNextExpectedBand)
323 : {
324 0 : poRDS->nNextExpectedBand++;
325 0 : if (poRDS->nNextExpectedBand > poRDS->GetRasterCount())
326 : {
327 0 : CPLDebug("PostGIS_Raster", "Re-enabling multi-band caching");
328 0 : poRDS->bAssumeMultiBandReadPattern = true;
329 0 : poRDS->nNextExpectedBand = 1;
330 : }
331 : }
332 : }
333 :
334 : #ifdef DEBUG_VERBOSE
335 : CPLDebug("PostGIS_Raster",
336 : "PostGISRasterRasterBand::IRasterIO: "
337 : "nBand = %d, nXOff = %d, nYOff = %d, nXSize = %d, nYSize = %d, "
338 : "nBufXSize = %d, nBufYSize = %d",
339 : nBand, nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize);
340 : #endif
341 :
342 : /*******************************************************************
343 : * Several tiles: we first look in all our sources caches. Missing
344 : * blocks are queried
345 : ******************************************************************/
346 : double adfProjWin[8];
347 0 : int nFeatureCount = 0;
348 : CPLRectObj sAoi;
349 :
350 0 : poRDS->PolygonFromCoords(nXOff, nYOff, nXOff + nXSize, nYOff + nYSize,
351 : adfProjWin);
352 : // (p[6], p[7]) is the minimum (x, y), and (p[2], p[3]) the max
353 0 : sAoi.minx = adfProjWin[6];
354 0 : sAoi.maxx = adfProjWin[2];
355 0 : if (adfProjWin[7] < adfProjWin[3])
356 : {
357 0 : sAoi.miny = adfProjWin[7];
358 0 : sAoi.maxy = adfProjWin[3];
359 : }
360 : else
361 : {
362 0 : sAoi.maxy = adfProjWin[7];
363 0 : sAoi.miny = adfProjWin[3];
364 : }
365 :
366 : #ifdef DEBUG_VERBOSE
367 : CPLDebug("PostGIS_Raster",
368 : "PostGISRasterRasterBand::IRasterIO: "
369 : "Intersection box: (%f, %f) - (%f, %f)",
370 : sAoi.minx, sAoi.miny, sAoi.maxx, sAoi.maxy);
371 : #endif
372 :
373 0 : if (poRDS->hQuadTree == nullptr)
374 : {
375 0 : ReportError(CE_Failure, CPLE_AppDefined,
376 : "Could not read metadata index.");
377 0 : return CE_Failure;
378 : }
379 :
380 0 : NullBuffer(pData, nBufXSize, nBufYSize, eBufType,
381 : static_cast<int>(nPixelSpace), static_cast<int>(nLineSpace));
382 :
383 0 : if (poRDS->bBuildQuadTreeDynamically && !bSameWindowAsOtherBand)
384 : {
385 0 : if (!(poRDS->LoadSources(nXOff, nYOff, nXSize, nYSize, nBand)))
386 0 : return CE_Failure;
387 : }
388 :
389 : // Matching sources, to avoid a dumb for loop over the sources
390 : PostGISRasterTileDataset **papsMatchingTiles =
391 : reinterpret_cast<PostGISRasterTileDataset **>(
392 0 : CPLQuadTreeSearch(poRDS->hQuadTree, &sAoi, &nFeatureCount));
393 :
394 : // No blocks found. This is not an error (the raster may have holes)
395 0 : if (nFeatureCount == 0)
396 : {
397 0 : CPLFree(papsMatchingTiles);
398 :
399 0 : return CE_None;
400 : }
401 :
402 : int i;
403 :
404 : /**
405 : * We need to store the max, min coords for the missing tiles in
406 : * any place. This is as good as any other
407 : **/
408 0 : sAoi.minx = 0.0;
409 0 : sAoi.miny = 0.0;
410 0 : sAoi.maxx = 0.0;
411 0 : sAoi.maxy = 0.0;
412 :
413 0 : GIntBig nMemoryRequiredForTiles = 0;
414 0 : CPLString osIDsToFetch;
415 0 : int nTilesToFetch = 0;
416 0 : const int nBandDataTypeSize = GDALGetDataTypeSizeBytes(eDataType);
417 :
418 : // Loop just over the intersecting sources
419 0 : for (i = 0; i < nFeatureCount; i++)
420 : {
421 0 : PostGISRasterTileDataset *poTile = papsMatchingTiles[i];
422 : PostGISRasterTileRasterBand *poTileBand =
423 0 : cpl::down_cast<PostGISRasterTileRasterBand *>(
424 : poTile->GetRasterBand(nBand));
425 :
426 0 : nMemoryRequiredForTiles +=
427 0 : static_cast<GIntBig>(poTileBand->GetXSize()) *
428 0 : poTileBand->GetYSize() * nBandDataTypeSize;
429 :
430 : // Missing tile: we'll need to query for it
431 0 : if (!poTileBand->IsCached())
432 : {
433 :
434 : // If we have a PKID, add the tile PKID to the list
435 0 : if (poTile->pszPKID != nullptr)
436 : {
437 0 : if (!osIDsToFetch.empty())
438 0 : osIDsToFetch += ",";
439 0 : osIDsToFetch += "'";
440 0 : osIDsToFetch += poTile->pszPKID;
441 0 : osIDsToFetch += "'";
442 : }
443 :
444 : double dfTileMinX, dfTileMinY, dfTileMaxX, dfTileMaxY;
445 0 : poTile->GetNativeExtent(&dfTileMinX, &dfTileMinY, &dfTileMaxX,
446 : &dfTileMaxY);
447 :
448 : /**
449 : * We keep the general max and min values of all the missing
450 : * tiles, to raise a query that intersect just that area.
451 : *
452 : * TODO: In case of just a few tiles and very separated,
453 : * this strategy is clearly suboptimal. We'll get our
454 : * missing tiles, but with a lot of other not needed tiles.
455 : *
456 : * A possible optimization will be to simply rely on the
457 : * I/O method of the source (must be implemented), in case
458 : * we have minus than a reasonable amount of tiles missing.
459 : * Another criteria to decide would be how separated the
460 : * tiles are. Two queries for just two adjacent tiles is
461 : * also a dumb strategy.
462 : **/
463 0 : if (nTilesToFetch == 0)
464 : {
465 0 : sAoi.minx = dfTileMinX;
466 0 : sAoi.miny = dfTileMinY;
467 0 : sAoi.maxx = dfTileMaxX;
468 0 : sAoi.maxy = dfTileMaxY;
469 : }
470 : else
471 : {
472 0 : if (dfTileMinX < sAoi.minx)
473 0 : sAoi.minx = dfTileMinX;
474 :
475 0 : if (dfTileMinY < sAoi.miny)
476 0 : sAoi.miny = dfTileMinY;
477 :
478 0 : if (dfTileMaxX > sAoi.maxx)
479 0 : sAoi.maxx = dfTileMaxX;
480 :
481 0 : if (dfTileMaxY > sAoi.maxy)
482 0 : sAoi.maxy = dfTileMaxY;
483 : }
484 :
485 0 : nTilesToFetch++;
486 : }
487 : }
488 :
489 : /* Determine caching strategy */
490 0 : bool bAllBandCaching = false;
491 0 : if (nTilesToFetch > 0)
492 : {
493 0 : GIntBig nCacheMax = GDALGetCacheMax64();
494 0 : if (nMemoryRequiredForTiles > nCacheMax)
495 : {
496 0 : CPLDebug("PostGIS_Raster",
497 : "For best performance, the block cache should be able to "
498 : "store " CPL_FRMT_GIB
499 : " bytes for the tiles of the requested window, "
500 : "but it is only " CPL_FRMT_GIB " byte large",
501 : nMemoryRequiredForTiles, nCacheMax);
502 0 : nTilesToFetch = 0;
503 : }
504 :
505 0 : if (poRDS->GetRasterCount() > 1 && poRDS->bAssumeMultiBandReadPattern)
506 : {
507 : GIntBig nMemoryRequiredForTilesAllBands =
508 0 : nMemoryRequiredForTiles * poRDS->GetRasterCount();
509 0 : if (nMemoryRequiredForTilesAllBands <= nCacheMax)
510 : {
511 0 : bAllBandCaching = true;
512 : }
513 : else
514 : {
515 0 : CPLDebug("PostGIS_Raster",
516 : "Caching only this band, but not all bands. "
517 : "Cache should be " CPL_FRMT_GIB " byte large for that",
518 : nMemoryRequiredForTilesAllBands);
519 : }
520 : }
521 : }
522 :
523 : // Raise a query for missing tiles and cache them
524 0 : if (nTilesToFetch > 0)
525 : {
526 :
527 : /**
528 : * There are several options here, to raise the query.
529 : * - Get all the tiles which PKID is in a list of missing
530 : * PKIDs.
531 : * - Get all the tiles that intersect a polygon constructed
532 : * based on the (min - max) values calculated before.
533 : * - Get all the tiles with upper left pixel included in the
534 : * range (min - max) calculated before.
535 : *
536 : * The first option is the most efficient one when a PKID exists.
537 : * After that, the second one is the most efficient one when a
538 : * spatial index exists.
539 : * The third one is the only one available when neither a PKID or
540 : *spatial index exist.
541 : **/
542 :
543 0 : CPLString osSchemaI(CPLQuotedSQLIdentifier(pszSchema));
544 0 : CPLString osTableI(CPLQuotedSQLIdentifier(pszTable));
545 0 : CPLString osColumnI(CPLQuotedSQLIdentifier(pszColumn));
546 :
547 0 : CPLString osWHERE;
548 0 : if (!osIDsToFetch.empty() &&
549 0 : (poRDS->bIsFastPK || !(poRDS->HasSpatialIndex())))
550 : {
551 0 : if (nTilesToFetch < poRDS->m_nTiles ||
552 0 : poRDS->bBuildQuadTreeDynamically)
553 : {
554 0 : osWHERE += poRDS->pszPrimaryKeyName;
555 0 : osWHERE += " IN (";
556 0 : osWHERE += osIDsToFetch;
557 0 : osWHERE += ")";
558 : }
559 : }
560 : else
561 : {
562 0 : if (poRDS->HasSpatialIndex())
563 : {
564 : osWHERE += CPLSPrintf(
565 : "%s && "
566 : "ST_GeomFromText('POLYGON((%.18f %.18f,%.18f %.18f,%.18f "
567 : "%.18f,%.18f %.18f,%.18f %.18f))')",
568 : osColumnI.c_str(), adfProjWin[0], adfProjWin[1],
569 : adfProjWin[2], adfProjWin[3], adfProjWin[4], adfProjWin[5],
570 0 : adfProjWin[6], adfProjWin[7], adfProjWin[0], adfProjWin[1]);
571 : }
572 : else
573 : {
574 : #define EPS 1e-5
575 : osWHERE += CPLSPrintf(
576 : "ST_UpperLeftX(%s)"
577 : " BETWEEN %f AND %f AND ST_UpperLeftY(%s) BETWEEN "
578 : "%f AND %f",
579 0 : osColumnI.c_str(), sAoi.minx - EPS, sAoi.maxx + EPS,
580 0 : osColumnI.c_str(), sAoi.miny - EPS, sAoi.maxy + EPS);
581 : }
582 : }
583 :
584 0 : if (poRDS->pszWhere != nullptr)
585 : {
586 0 : if (!osWHERE.empty())
587 0 : osWHERE += " AND ";
588 0 : osWHERE += "(";
589 0 : osWHERE += poRDS->pszWhere;
590 0 : osWHERE += ")";
591 : }
592 :
593 0 : bool bCanUseClientSide = true;
594 0 : if (poRDS->eOutDBResolution == OutDBResolution::CLIENT_SIDE_IF_POSSIBLE)
595 : {
596 : bCanUseClientSide =
597 0 : poRDS->CanUseClientSideOutDB(bAllBandCaching, nBand, osWHERE);
598 : }
599 :
600 0 : CPLString osRasterToFetch;
601 0 : if (bAllBandCaching)
602 0 : osRasterToFetch = osColumnI;
603 : else
604 0 : osRasterToFetch.Printf("ST_Band(%s, %d)", osColumnI.c_str(), nBand);
605 0 : if (poRDS->eOutDBResolution == OutDBResolution::SERVER_SIDE ||
606 0 : !bCanUseClientSide)
607 : {
608 : osRasterToFetch =
609 0 : "encode(ST_AsBinary(" + osRasterToFetch + ",TRUE),'hex')";
610 : }
611 :
612 0 : CPLString osCommand;
613 : osCommand.Printf("SELECT %s, ST_Metadata(%s), %s FROM %s.%s",
614 0 : (poRDS->GetPrimaryKeyRef()) ? poRDS->GetPrimaryKeyRef()
615 : : "NULL",
616 : osColumnI.c_str(), osRasterToFetch.c_str(),
617 0 : osSchemaI.c_str(), osTableI.c_str());
618 0 : if (!osWHERE.empty())
619 : {
620 0 : osCommand += " WHERE " + osWHERE;
621 : }
622 :
623 0 : PGresult *poResult = PQexec(poRDS->poConn, osCommand.c_str());
624 :
625 : #ifdef DEBUG_QUERY
626 : CPLDebug("PostGIS_Raster",
627 : "PostGISRasterRasterBand::IRasterIO(): Query = \"%s\" --> "
628 : "number of rows = %d",
629 : osCommand.c_str(), poResult ? PQntuples(poResult) : 0);
630 : #endif
631 :
632 0 : if (poResult == nullptr ||
633 0 : PQresultStatus(poResult) != PGRES_TUPLES_OK ||
634 0 : PQntuples(poResult) < 0)
635 : {
636 :
637 0 : if (poResult)
638 0 : PQclear(poResult);
639 :
640 0 : CPLError(CE_Failure, CPLE_AppDefined,
641 : "PostGISRasterRasterBand::IRasterIO(): %s",
642 0 : PQerrorMessage(poRDS->poConn));
643 :
644 : // Free the object that holds pointers to matching tiles
645 0 : CPLFree(papsMatchingTiles);
646 0 : return CE_Failure;
647 : }
648 :
649 : /**
650 : * No data. Return the buffer filled with nodata values
651 : **/
652 0 : else if (PQntuples(poResult) == 0)
653 : {
654 0 : PQclear(poResult);
655 :
656 : // Free the object that holds pointers to matching tiles
657 0 : CPLFree(papsMatchingTiles);
658 0 : return CE_None;
659 : }
660 :
661 : /**
662 : * Ok, we loop over the results
663 : **/
664 0 : int nTuples = PQntuples(poResult);
665 0 : for (i = 0; i < nTuples; i++)
666 : {
667 0 : const char *pszPKID = PQgetvalue(poResult, i, 0);
668 0 : const char *pszMetadata = PQgetvalue(poResult, i, 1);
669 0 : const char *pszRaster = PQgetvalue(poResult, i, 2);
670 0 : poRDS->CacheTile(pszMetadata, pszRaster, pszPKID, nBand,
671 : bAllBandCaching);
672 : } // All tiles have been added to cache
673 :
674 0 : PQclear(poResult);
675 : } // End missing tiles
676 :
677 : /* -------------------------------------------------------------------- */
678 : /* Overlay each source in turn over top this. */
679 : /* -------------------------------------------------------------------- */
680 :
681 0 : CPLErr eErr = CE_None;
682 : /* Sort tiles by ascending PKID, so that the draw order is deterministic. */
683 0 : if (poRDS->GetPrimaryKeyRef() != nullptr)
684 : {
685 0 : qsort(papsMatchingTiles, nFeatureCount,
686 : sizeof(PostGISRasterTileDataset *), SortTilesByPKID);
687 : }
688 :
689 0 : VRTSource::WorkingState oWorkingState;
690 0 : for (i = 0; i < nFeatureCount && eErr == CE_None; i++)
691 : {
692 0 : PostGISRasterTileDataset *poTile = papsMatchingTiles[i];
693 : PostGISRasterTileRasterBand *poTileBand =
694 0 : cpl::down_cast<PostGISRasterTileRasterBand *>(
695 : poTile->GetRasterBand(nBand));
696 0 : eErr = poTileBand->poSource->RasterIO(
697 : eDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
698 : nBufYSize, eBufType, nPixelSpace, nLineSpace, nullptr,
699 0 : oWorkingState);
700 : }
701 :
702 : // Free the object that holds pointers to matching tiles
703 0 : CPLFree(papsMatchingTiles);
704 :
705 0 : return eErr;
706 : }
707 :
708 : /**
709 : * \brief Set the no data value for this band.
710 : * Parameters:
711 : * - double: The nodata value
712 : * Returns:
713 : * - CE_None.
714 : */
715 0 : CPLErr PostGISRasterRasterBand::SetNoDataValue(double dfNewValue)
716 : {
717 0 : m_dfNoDataValue = dfNewValue;
718 :
719 0 : return CE_None;
720 : }
721 :
722 : /**
723 : * \brief Fetch the no data value for this band.
724 : * Parameters:
725 : * - int *: pointer to a boolean to use to indicate if a value is actually
726 : * associated with this layer. May be NULL (default).
727 : * Returns:
728 : * - double: the nodata value for this band.
729 : */
730 0 : double PostGISRasterRasterBand::GetNoDataValue(int *pbSuccess)
731 : {
732 0 : if (pbSuccess != nullptr)
733 0 : *pbSuccess = m_bNoDataValueSet;
734 :
735 0 : return m_dfNoDataValue;
736 : }
737 :
738 : /***************************************************
739 : * \brief Return the number of overview layers available
740 : ***************************************************/
741 0 : int PostGISRasterRasterBand::GetOverviewCount()
742 : {
743 0 : PostGISRasterDataset *poRDS = cpl::down_cast<PostGISRasterDataset *>(poDS);
744 0 : return poRDS->GetOverviewCount();
745 : }
746 :
747 : /**********************************************************
748 : * \brief Fetch overview raster band object
749 : **********************************************************/
750 0 : GDALRasterBand *PostGISRasterRasterBand::GetOverview(int i)
751 : {
752 0 : if (i < 0 || i >= GetOverviewCount())
753 0 : return nullptr;
754 :
755 0 : PostGISRasterDataset *poRDS = cpl::down_cast<PostGISRasterDataset *>(poDS);
756 0 : PostGISRasterDataset *poOverviewDS = poRDS->GetOverviewDS(i);
757 0 : if (!poOverviewDS)
758 : {
759 0 : CPLAssert(false);
760 : return nullptr;
761 : }
762 0 : if (poOverviewDS->nBands == 0)
763 : {
764 0 : if (!poOverviewDS->SetRasterProperties(nullptr) ||
765 0 : poOverviewDS->GetRasterCount() != poRDS->GetRasterCount())
766 : {
767 0 : CPLDebug("PostGIS_Raster",
768 : "Request for overview %d of band %d failed", i, nBand);
769 0 : return nullptr;
770 : }
771 : }
772 :
773 0 : return poOverviewDS->GetRasterBand(nBand);
774 : }
775 :
776 : /**
777 : * \brief How should this band be interpreted as color?
778 : * GCI_Undefined is returned when the format doesn't know anything about the
779 : * color interpretation.
780 : **/
781 0 : GDALColorInterp PostGISRasterRasterBand::GetColorInterpretation()
782 : {
783 0 : if (poDS->GetRasterCount() == 1)
784 : {
785 0 : m_eColorInterp = GCI_GrayIndex;
786 : }
787 :
788 0 : else if (poDS->GetRasterCount() == 3)
789 : {
790 0 : if (nBand == 1)
791 0 : m_eColorInterp = GCI_RedBand;
792 0 : else if (nBand == 2)
793 0 : m_eColorInterp = GCI_GreenBand;
794 0 : else if (nBand == 3)
795 0 : m_eColorInterp = GCI_BlueBand;
796 : else
797 0 : m_eColorInterp = GCI_Undefined;
798 : }
799 :
800 : else
801 : {
802 0 : m_eColorInterp = GCI_Undefined;
803 : }
804 :
805 0 : return m_eColorInterp;
806 : }
807 :
808 : /************************************************************************/
809 : /* GetMinimum() */
810 : /************************************************************************/
811 :
812 0 : double PostGISRasterRasterBand::GetMinimum(int *pbSuccess)
813 : {
814 0 : if (!m_bStatsFetched)
815 : {
816 0 : CPL_IGNORE_RET_VAL(QueryStats());
817 : }
818 :
819 0 : if (StatsFetchedAndValid())
820 : {
821 0 : if (pbSuccess)
822 0 : *pbSuccess = TRUE;
823 0 : return m_dfStatsMin;
824 : }
825 :
826 0 : PostGISRasterDataset *poRDS = cpl::down_cast<PostGISRasterDataset *>(poDS);
827 0 : if (poRDS->bBuildQuadTreeDynamically && poRDS->m_nTiles == 0)
828 : {
829 0 : if (pbSuccess)
830 0 : *pbSuccess = FALSE;
831 0 : return 0.0;
832 : }
833 0 : return VRTSourcedRasterBand::GetMinimum(pbSuccess);
834 : }
835 :
836 : /************************************************************************/
837 : /* GetMaximum() */
838 : /************************************************************************/
839 :
840 0 : double PostGISRasterRasterBand::GetMaximum(int *pbSuccess)
841 : {
842 0 : if (!m_bStatsFetched)
843 : {
844 0 : CPL_IGNORE_RET_VAL(QueryStats());
845 : }
846 :
847 0 : if (StatsFetchedAndValid())
848 : {
849 0 : if (pbSuccess)
850 0 : *pbSuccess = TRUE;
851 0 : return m_dfStatsMax;
852 : }
853 :
854 0 : PostGISRasterDataset *poRDS = cpl::down_cast<PostGISRasterDataset *>(poDS);
855 0 : if (poRDS->bBuildQuadTreeDynamically && poRDS->m_nTiles == 0)
856 : {
857 0 : if (pbSuccess)
858 0 : *pbSuccess = FALSE;
859 0 : return 0.0;
860 : }
861 0 : return VRTSourcedRasterBand::GetMaximum(pbSuccess);
862 : }
863 :
864 : /************************************************************************/
865 : /* ComputeRasterMinMax() */
866 : /************************************************************************/
867 :
868 0 : CPLErr PostGISRasterRasterBand::ComputeRasterMinMax(int bApproxOK,
869 : double *adfMinMax)
870 : {
871 0 : if (StatsFetchedAndValid())
872 : {
873 0 : return CE_None;
874 : }
875 :
876 0 : if (nRasterXSize < 1024 && nRasterYSize < 1024)
877 0 : return VRTSourcedRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
878 :
879 0 : int nOverviewCount = GetOverviewCount();
880 0 : for (int i = 0; i < nOverviewCount; i++)
881 : {
882 0 : auto poOverview = GetOverview(i);
883 0 : if (poOverview->GetXSize() < 1024 && poOverview->GetYSize() < 1024)
884 0 : return poOverview->ComputeRasterMinMax(bApproxOK, adfMinMax);
885 : }
886 :
887 : // Try to fetch the min/max from the database
888 0 : if (QueryStats())
889 : {
890 0 : return CE_None;
891 : }
892 : else
893 : {
894 0 : return CE_Failure;
895 : }
896 : }
897 :
898 : /************************************************************************/
899 : /* ComputeStatistics() */
900 : /************************************************************************/
901 0 : CPLErr PostGISRasterRasterBand::ComputeStatistics(
902 : int bApproxOK, double *pdfMin, double *pdfMax, double *pdfMean,
903 : double *pdfStdDev, GDALProgressFunc pfnProgress, void *pProgressData,
904 : CSLConstList papszOptions)
905 : {
906 0 : if (!m_bStatsFetched)
907 : {
908 0 : CPL_IGNORE_RET_VAL(QueryStats());
909 : }
910 :
911 0 : if (StatsFetchedAndValid())
912 : {
913 0 : *pdfMin = m_dfStatsMin;
914 0 : *pdfMax = m_dfStatsMax;
915 0 : *pdfMean = m_dfStatsMean;
916 0 : *pdfStdDev = m_dfStatsStdDev;
917 0 : if (CPLFetchBool(papszOptions, "SET_STATISTICS", true))
918 : {
919 0 : SetStatistics(m_dfStatsMin, m_dfStatsMax, m_dfStatsMean,
920 0 : m_dfStatsStdDev);
921 : }
922 0 : return CE_None;
923 : }
924 : else
925 : {
926 0 : return VRTSourcedRasterBand::ComputeStatistics(
927 : bApproxOK, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
928 0 : pProgressData, papszOptions);
929 : }
930 : }
|
