Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GeoPackage Translator
4 : * Purpose: Implements GDALGeoPackageRasterBand class
5 : * Author: Even Rouault <even dot rouault at spatialys dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2014, Even Rouault <even dot rouault at spatialys dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "ogr_geopackage.h"
14 : #include "memdataset.h"
15 : #include "gdal_alg_priv.h"
16 : #include "ogrsqlitevfs.h"
17 : #include "cpl_error.h"
18 : #include "cpl_float.h"
19 :
20 : #include <algorithm>
21 : #include <cmath>
22 : #include <limits>
23 : #include <set>
24 : #include <utility>
25 :
26 : #if !defined(DEBUG_VERBOSE) && defined(DEBUG_VERBOSE_GPKG)
27 : #define DEBUG_VERBOSE
28 : #endif
29 :
30 : /************************************************************************/
31 : /* GDALGPKGMBTilesLikePseudoDataset() */
32 : /************************************************************************/
33 :
34 2876 : GDALGPKGMBTilesLikePseudoDataset::GDALGPKGMBTilesLikePseudoDataset()
35 : : m_bForceTempDBCompaction(
36 2876 : CPLTestBool(CPLGetConfigOption("GPKG_FORCE_TEMPDB_COMPACTION", "NO")))
37 : {
38 14380 : for (int i = 0; i < 4; i++)
39 : {
40 11504 : m_asCachedTilesDesc[i].nRow = -1;
41 11504 : m_asCachedTilesDesc[i].nCol = -1;
42 11504 : m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
43 11504 : m_asCachedTilesDesc[i].abBandDirty[0] = FALSE;
44 11504 : m_asCachedTilesDesc[i].abBandDirty[1] = FALSE;
45 11504 : m_asCachedTilesDesc[i].abBandDirty[2] = FALSE;
46 11504 : m_asCachedTilesDesc[i].abBandDirty[3] = FALSE;
47 : }
48 2876 : }
49 :
50 : /************************************************************************/
51 : /* ~GDALGPKGMBTilesLikePseudoDataset() */
52 : /************************************************************************/
53 :
54 2876 : GDALGPKGMBTilesLikePseudoDataset::~GDALGPKGMBTilesLikePseudoDataset()
55 : {
56 2876 : if (m_poParentDS == nullptr && m_hTempDB != nullptr)
57 : {
58 48 : sqlite3_close(m_hTempDB);
59 48 : m_hTempDB = nullptr;
60 48 : VSIUnlink(m_osTempDBFilename);
61 48 : if (m_pMyVFS)
62 : {
63 29 : sqlite3_vfs_unregister(m_pMyVFS);
64 29 : CPLFree(m_pMyVFS->pAppData);
65 29 : CPLFree(m_pMyVFS);
66 : }
67 : }
68 2876 : CPLFree(m_pabyCachedTiles);
69 2876 : delete m_poCT;
70 2876 : CPLFree(m_pabyHugeColorArray);
71 2876 : }
72 :
73 : /************************************************************************/
74 : /* SetDataType() */
75 : /************************************************************************/
76 :
77 322 : void GDALGPKGMBTilesLikePseudoDataset::SetDataType(GDALDataType eDT)
78 : {
79 322 : CPLAssert(eDT == GDT_UInt8 || eDT == GDT_Int16 || eDT == GDT_UInt16 ||
80 : eDT == GDT_Float32);
81 322 : m_eDT = eDT;
82 322 : m_nDTSize = GDALGetDataTypeSizeBytes(m_eDT);
83 322 : }
84 :
85 : /************************************************************************/
86 : /* SetGlobalOffsetScale() */
87 : /************************************************************************/
88 :
89 75 : void GDALGPKGMBTilesLikePseudoDataset::SetGlobalOffsetScale(double dfOffset,
90 : double dfScale)
91 : {
92 75 : m_dfOffset = dfOffset;
93 75 : m_dfScale = dfScale;
94 75 : }
95 :
96 : /************************************************************************/
97 : /* GDALGPKGMBTilesLikeRasterBand() */
98 : /************************************************************************/
99 :
100 : // Recent GCC versions complain about null dereference of m_poTPD in
101 : // the constructor body
102 : #ifdef __GNUC__
103 : #pragma GCC diagnostic push
104 : #pragma GCC diagnostic ignored "-Wnull-dereference"
105 : #endif
106 :
107 2559 : GDALGPKGMBTilesLikeRasterBand::GDALGPKGMBTilesLikeRasterBand(
108 2559 : GDALGPKGMBTilesLikePseudoDataset *poTPD, int nTileWidth, int nTileHeight)
109 2559 : : m_poTPD(poTPD)
110 : {
111 2559 : eDataType = m_poTPD->m_eDT;
112 2559 : m_nDTSize = m_poTPD->m_nDTSize;
113 2559 : nBlockXSize = nTileWidth;
114 2559 : nBlockYSize = nTileHeight;
115 2559 : }
116 :
117 : #ifdef __GNUC__
118 : #pragma GCC diagnostic pop
119 : #endif
120 :
121 : /************************************************************************/
122 : /* FlushCache() */
123 : /************************************************************************/
124 :
125 5518 : CPLErr GDALGPKGMBTilesLikeRasterBand::FlushCache(bool bAtClosing)
126 : {
127 5518 : m_poTPD->m_nLastSpaceCheckTimestamp = -1; // disable partial flushes
128 5518 : CPLErr eErr = GDALPamRasterBand::FlushCache(bAtClosing);
129 5518 : if (eErr == CE_None)
130 5515 : eErr = m_poTPD->IFlushCacheWithErrCode(bAtClosing);
131 5518 : m_poTPD->m_nLastSpaceCheckTimestamp = 0;
132 5518 : return eErr;
133 : }
134 :
135 : /************************************************************************/
136 : /* FlushTiles() */
137 : /************************************************************************/
138 :
139 3638 : CPLErr GDALGPKGMBTilesLikePseudoDataset::FlushTiles()
140 : {
141 3638 : CPLErr eErr = CE_None;
142 3638 : GDALGPKGMBTilesLikePseudoDataset *poMainDS =
143 3638 : m_poParentDS ? m_poParentDS : this;
144 3638 : if (poMainDS->m_nTileInsertionCount < 0)
145 7 : return CE_Failure;
146 :
147 3631 : if (IGetUpdate())
148 : {
149 2189 : if (m_nShiftXPixelsMod || m_nShiftYPixelsMod)
150 : {
151 607 : eErr = FlushRemainingShiftedTiles(false /* total flush*/);
152 : }
153 : else
154 : {
155 1582 : eErr = WriteTile();
156 : }
157 : }
158 :
159 3631 : if (poMainDS->m_nTileInsertionCount > 0)
160 : {
161 190 : if (poMainDS->ICommitTransaction() != OGRERR_NONE)
162 : {
163 0 : poMainDS->m_nTileInsertionCount = -1;
164 0 : eErr = CE_Failure;
165 : }
166 : else
167 : {
168 190 : poMainDS->m_nTileInsertionCount = 0;
169 : }
170 : }
171 3631 : return eErr;
172 : }
173 :
174 : /************************************************************************/
175 : /* GetColorTable() */
176 : /************************************************************************/
177 :
178 352 : GDALColorTable *GDALGPKGMBTilesLikeRasterBand::GetColorTable()
179 : {
180 352 : if (poDS->GetRasterCount() != 1)
181 127 : return nullptr;
182 :
183 225 : if (!m_poTPD->m_bTriedEstablishingCT)
184 : {
185 66 : m_poTPD->m_bTriedEstablishingCT = true;
186 66 : if (m_poTPD->m_poParentDS != nullptr)
187 : {
188 16 : m_poTPD->m_poCT =
189 8 : m_poTPD->m_poParentDS->IGetRasterBand(1)->GetColorTable();
190 8 : if (m_poTPD->m_poCT)
191 4 : m_poTPD->m_poCT = m_poTPD->m_poCT->Clone();
192 8 : return m_poTPD->m_poCT;
193 : }
194 :
195 70 : for (int i = 0; i < 2; i++)
196 : {
197 64 : bool bRetry = false;
198 64 : char *pszSQL = nullptr;
199 64 : if (i == 0)
200 : {
201 58 : pszSQL = sqlite3_mprintf("SELECT tile_data FROM \"%w\" "
202 : "WHERE zoom_level = %d LIMIT 1",
203 58 : m_poTPD->m_osRasterTable.c_str(),
204 58 : m_poTPD->m_nZoomLevel);
205 : }
206 : else
207 : {
208 : // Try a tile in the middle of the raster
209 6 : pszSQL = sqlite3_mprintf(
210 : "SELECT tile_data FROM \"%w\" "
211 : "WHERE zoom_level = %d AND tile_column = %d AND tile_row = "
212 : "%d",
213 6 : m_poTPD->m_osRasterTable.c_str(), m_poTPD->m_nZoomLevel,
214 6 : m_poTPD->m_nShiftXTiles + nRasterXSize / 2 / nBlockXSize,
215 6 : m_poTPD->GetRowFromIntoTopConvention(
216 6 : m_poTPD->m_nShiftYTiles +
217 6 : nRasterYSize / 2 / nBlockYSize));
218 : }
219 64 : sqlite3_stmt *hStmt = nullptr;
220 64 : int rc = SQLPrepareWithError(m_poTPD->IGetDB(), pszSQL, -1, &hStmt,
221 : nullptr);
222 64 : if (rc == SQLITE_OK)
223 : {
224 64 : rc = sqlite3_step(hStmt);
225 81 : if (rc == SQLITE_ROW &&
226 17 : sqlite3_column_type(hStmt, 0) == SQLITE_BLOB)
227 : {
228 17 : const int nBytes = sqlite3_column_bytes(hStmt, 0);
229 : GByte *pabyRawData = reinterpret_cast<GByte *>(
230 17 : const_cast<void *>(sqlite3_column_blob(hStmt, 0)));
231 : const CPLString osMemFileName(
232 34 : VSIMemGenerateHiddenFilename("gpkg_read_tile"));
233 17 : VSILFILE *fp = VSIFileFromMemBuffer(
234 : osMemFileName.c_str(), pabyRawData, nBytes, FALSE);
235 17 : VSIFCloseL(fp);
236 :
237 : // Only PNG can have color table.
238 17 : const char *const apszDrivers[] = {"PNG", nullptr};
239 : auto poDSTile = std::unique_ptr<GDALDataset>(
240 : GDALDataset::Open(osMemFileName.c_str(),
241 : GDAL_OF_RASTER | GDAL_OF_INTERNAL,
242 34 : apszDrivers, nullptr, nullptr));
243 17 : if (poDSTile != nullptr)
244 : {
245 17 : if (poDSTile->GetRasterCount() == 1)
246 : {
247 10 : m_poTPD->m_poCT =
248 5 : poDSTile->GetRasterBand(1)->GetColorTable();
249 5 : if (m_poTPD->m_poCT != nullptr)
250 1 : m_poTPD->m_poCT = m_poTPD->m_poCT->Clone();
251 : }
252 : else
253 : {
254 12 : bRetry = true;
255 : }
256 : }
257 : else
258 0 : bRetry = true;
259 :
260 17 : VSIUnlink(osMemFileName);
261 : }
262 : }
263 64 : sqlite3_free(pszSQL);
264 64 : sqlite3_finalize(hStmt);
265 64 : if (!bRetry)
266 52 : break;
267 : }
268 : }
269 :
270 217 : return m_poTPD->m_poCT;
271 : }
272 :
273 : /************************************************************************/
274 : /* SetColorTable() */
275 : /************************************************************************/
276 :
277 15 : CPLErr GDALGPKGMBTilesLikeRasterBand::SetColorTable(GDALColorTable *poCT)
278 : {
279 15 : if (m_poTPD->m_eDT != GDT_UInt8)
280 0 : return CE_Failure;
281 15 : if (poDS->GetRasterCount() != 1)
282 : {
283 1 : CPLError(CE_Failure, CPLE_NotSupported,
284 : "SetColorTable() only supported for a single band dataset");
285 1 : return CE_Failure;
286 : }
287 14 : if (!m_poTPD->m_bNew || m_poTPD->m_bTriedEstablishingCT)
288 : {
289 2 : CPLError(CE_Failure, CPLE_NotSupported,
290 : "SetColorTable() only supported on a newly created dataset");
291 2 : return CE_Failure;
292 : }
293 :
294 12 : AssignColorTable(poCT);
295 12 : return CE_None;
296 : }
297 :
298 : /************************************************************************/
299 : /* AssignColorTable() */
300 : /************************************************************************/
301 :
302 15 : void GDALGPKGMBTilesLikeRasterBand::AssignColorTable(const GDALColorTable *poCT)
303 : {
304 15 : m_poTPD->m_bTriedEstablishingCT = true;
305 15 : delete m_poTPD->m_poCT;
306 15 : if (poCT != nullptr)
307 14 : m_poTPD->m_poCT = poCT->Clone();
308 : else
309 1 : m_poTPD->m_poCT = nullptr;
310 15 : }
311 :
312 : /************************************************************************/
313 : /* GetColorInterpretation() */
314 : /************************************************************************/
315 :
316 212 : GDALColorInterp GDALGPKGMBTilesLikeRasterBand::GetColorInterpretation()
317 : {
318 212 : if (m_poTPD->m_eDT != GDT_UInt8)
319 23 : return GCI_Undefined;
320 189 : if (poDS->GetRasterCount() == 1)
321 51 : return GetColorTable() ? GCI_PaletteIndex : GCI_GrayIndex;
322 138 : else if (poDS->GetRasterCount() == 2)
323 5 : return (nBand == 1) ? GCI_GrayIndex : GCI_AlphaBand;
324 : else
325 133 : return static_cast<GDALColorInterp>(GCI_RedBand + (nBand - 1));
326 : }
327 :
328 : /************************************************************************/
329 : /* SetColorInterpretation() */
330 : /************************************************************************/
331 :
332 : CPLErr
333 31 : GDALGPKGMBTilesLikeRasterBand::SetColorInterpretation(GDALColorInterp eInterp)
334 : {
335 31 : if (eInterp == GCI_Undefined)
336 1 : return CE_None;
337 32 : if (poDS->GetRasterCount() == 1 &&
338 2 : (eInterp == GCI_GrayIndex || eInterp == GCI_PaletteIndex))
339 23 : return CE_None;
340 11 : if (poDS->GetRasterCount() == 2 &&
341 4 : ((nBand == 1 && eInterp == GCI_GrayIndex) ||
342 3 : (nBand == 2 && eInterp == GCI_AlphaBand)))
343 2 : return CE_None;
344 5 : if (poDS->GetRasterCount() >= 3 && eInterp == GCI_RedBand + nBand - 1)
345 1 : return CE_None;
346 4 : CPLError(CE_Warning, CPLE_NotSupported,
347 : "%s color interpretation not supported. Will be ignored",
348 : GDALGetColorInterpretationName(eInterp));
349 4 : return CE_Warning;
350 : }
351 :
352 : /************************************************************************/
353 : /* GPKGFindBestEntry() */
354 : /************************************************************************/
355 :
356 1 : static int GPKGFindBestEntry(GDALColorTable *poCT, GByte c1, GByte c2, GByte c3,
357 : GByte c4, int nTileBandCount)
358 : {
359 1 : const int nEntries = std::min(256, poCT->GetColorEntryCount());
360 1 : int iBestIdx = 0;
361 1 : int nBestDistance = 4 * 256 * 256;
362 257 : for (int i = 0; i < nEntries; i++)
363 : {
364 256 : const GDALColorEntry *psEntry = poCT->GetColorEntry(i);
365 256 : int nDistance = (psEntry->c1 - c1) * (psEntry->c1 - c1) +
366 256 : (psEntry->c2 - c2) * (psEntry->c2 - c2) +
367 256 : (psEntry->c3 - c3) * (psEntry->c3 - c3);
368 256 : if (nTileBandCount == 4)
369 256 : nDistance += (psEntry->c4 - c4) * (psEntry->c4 - c4);
370 256 : if (nDistance < nBestDistance)
371 : {
372 10 : iBestIdx = i;
373 10 : nBestDistance = nDistance;
374 : }
375 : }
376 1 : return iBestIdx;
377 : }
378 :
379 : /************************************************************************/
380 : /* FillBuffer() */
381 : /************************************************************************/
382 :
383 14296 : void GDALGPKGMBTilesLikePseudoDataset::FillBuffer(GByte *pabyData,
384 : size_t nPixels)
385 : {
386 14296 : int bHasNoData = FALSE;
387 14296 : const double dfNoDataValue = IGetRasterBand(1)->GetNoDataValue(&bHasNoData);
388 14296 : if (!bHasNoData || dfNoDataValue == 0.0)
389 : {
390 : // cppcheck-suppress nullPointer
391 13994 : memset(pabyData, 0, nPixels * m_nDTSize);
392 : }
393 : else
394 : {
395 302 : GDALCopyWords64(&dfNoDataValue, GDT_Float64, 0, pabyData, m_eDT,
396 : m_nDTSize, nPixels);
397 : }
398 14296 : }
399 :
400 : /************************************************************************/
401 : /* FillEmptyTile() */
402 : /************************************************************************/
403 :
404 2955 : void GDALGPKGMBTilesLikePseudoDataset::FillEmptyTile(GByte *pabyData)
405 : {
406 : int nBlockXSize, nBlockYSize;
407 2955 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
408 2955 : const int nBands = IGetRasterCount();
409 2955 : const size_t nPixels =
410 2955 : static_cast<size_t>(nBands) * nBlockXSize * nBlockYSize;
411 2955 : FillBuffer(pabyData, nPixels);
412 2955 : }
413 :
414 : /************************************************************************/
415 : /* FillEmptyTileSingleBand() */
416 : /************************************************************************/
417 :
418 1481 : void GDALGPKGMBTilesLikePseudoDataset::FillEmptyTileSingleBand(GByte *pabyData)
419 : {
420 : int nBlockXSize, nBlockYSize;
421 1481 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
422 1481 : const size_t nPixels = static_cast<size_t>(nBlockXSize) * nBlockYSize;
423 1481 : FillBuffer(pabyData, nPixels);
424 1481 : }
425 :
426 : /************************************************************************/
427 : /* ReadTile() */
428 : /************************************************************************/
429 :
430 2583 : CPLErr GDALGPKGMBTilesLikePseudoDataset::ReadTile(
431 : const CPLString &osMemFileName, GByte *pabyTileData, double dfTileOffset,
432 : double dfTileScale, bool *pbIsLossyFormat)
433 : {
434 2583 : const char *const apszDriversByte[] = {"JPEG", "PNG", "WEBP", nullptr};
435 2583 : const char *const apszDriversInt[] = {"PNG", nullptr};
436 2583 : const char *const apszDriversFloat[] = {"GTiff", nullptr};
437 : int nBlockXSize, nBlockYSize;
438 2583 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
439 2583 : const int nBands = IGetRasterCount();
440 : auto poDSTile = std::unique_ptr<GDALDataset>(GDALDataset::Open(
441 : osMemFileName.c_str(), GDAL_OF_RASTER | GDAL_OF_INTERNAL,
442 2583 : (m_eDT == GDT_UInt8) ? apszDriversByte
443 43 : : (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT) ? apszDriversFloat
444 : : apszDriversInt,
445 5166 : nullptr, nullptr));
446 2583 : if (poDSTile == nullptr)
447 : {
448 2 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot parse tile data");
449 2 : FillEmptyTile(pabyTileData);
450 2 : return CE_Failure;
451 : }
452 :
453 2581 : const int nTileBandCount = poDSTile->GetRasterCount();
454 :
455 5162 : if (!(poDSTile->GetRasterXSize() == nBlockXSize &&
456 2581 : poDSTile->GetRasterYSize() == nBlockYSize &&
457 7743 : (nTileBandCount >= 1 && nTileBandCount <= 4)) ||
458 2581 : (m_eDT != GDT_UInt8 && nTileBandCount != 1))
459 : {
460 0 : CPLError(CE_Failure, CPLE_AppDefined,
461 : "Inconsistent tiles characteristics");
462 0 : FillEmptyTile(pabyTileData);
463 0 : return CE_Failure;
464 : }
465 :
466 2581 : GDALDataType eRequestDT = GDT_UInt8;
467 2581 : if (m_eTF == GPKG_TF_PNG_16BIT)
468 : {
469 40 : CPLAssert(m_eDT == GDT_Int16 || m_eDT == GDT_UInt16 ||
470 : m_eDT == GDT_Float32);
471 40 : eRequestDT = GDT_UInt16;
472 : }
473 2541 : else if (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
474 : {
475 3 : CPLAssert(m_eDT == GDT_Float32);
476 3 : eRequestDT = GDT_Float32;
477 : }
478 :
479 2581 : if (poDSTile->RasterIO(GF_Read, 0, 0, nBlockXSize, nBlockYSize,
480 : pabyTileData, nBlockXSize, nBlockYSize, eRequestDT,
481 : poDSTile->GetRasterCount(), nullptr, 0, 0, 0,
482 2581 : nullptr) != CE_None)
483 : {
484 0 : FillEmptyTile(pabyTileData);
485 0 : return CE_Failure;
486 : }
487 :
488 2581 : if (m_eDT != GDT_UInt8)
489 : {
490 43 : int bHasNoData = FALSE;
491 : const double dfNoDataValue =
492 43 : IGetRasterBand(1)->GetNoDataValue(&bHasNoData);
493 43 : if (m_eDT == GDT_Int16)
494 : {
495 29 : CPLAssert(eRequestDT == GDT_UInt16);
496 29 : for (size_t i = 0;
497 1900570 : i < static_cast<size_t>(nBlockXSize) * nBlockYSize; i++)
498 : {
499 : GUInt16 usVal;
500 1900540 : memcpy(&usVal, pabyTileData + i * sizeof(GUInt16),
501 : sizeof(usVal));
502 : double dfVal =
503 1900540 : floor((usVal * dfTileScale + dfTileOffset) * m_dfScale +
504 1900540 : m_dfOffset + 0.5);
505 1900540 : if (bHasNoData && usVal == m_usGPKGNull)
506 195804 : dfVal = dfNoDataValue;
507 1900540 : if (dfVal > 32767)
508 0 : dfVal = 32767;
509 1900540 : else if (dfVal < -32768)
510 0 : dfVal = -32768;
511 1900540 : GInt16 sVal = static_cast<GInt16>(dfVal);
512 1900540 : memcpy(pabyTileData + i * sizeof(GUInt16), &sVal, sizeof(sVal));
513 : }
514 : }
515 14 : else if (m_eDT == GDT_UInt16 &&
516 6 : (m_dfOffset != 0.0 || m_dfScale != 1.0 ||
517 6 : dfTileOffset != 0.0 || dfTileScale != 1.0))
518 : {
519 0 : CPLAssert(eRequestDT == GDT_UInt16);
520 0 : GUInt16 *psVal = reinterpret_cast<GUInt16 *>(pabyTileData);
521 0 : for (size_t i = 0;
522 0 : i < static_cast<size_t>(nBlockXSize) * nBlockYSize; i++)
523 : {
524 0 : const GUInt16 nVal = psVal[i];
525 : double dfVal =
526 0 : floor((nVal * dfTileScale + dfTileOffset) * m_dfScale +
527 0 : m_dfOffset + 0.5);
528 0 : if (bHasNoData && nVal == m_usGPKGNull)
529 0 : dfVal = dfNoDataValue;
530 0 : if (dfVal > 65535)
531 0 : dfVal = 65535;
532 0 : else if (dfVal < 0)
533 0 : dfVal = 0;
534 0 : psVal[i] = static_cast<GUInt16>(dfVal);
535 0 : }
536 : }
537 14 : else if (m_eDT == GDT_Float32 && eRequestDT == GDT_UInt16)
538 : {
539 : // Due to non identical data type size, we need to start from the
540 : // end of the buffer.
541 327685 : for (GPtrDiff_t i =
542 5 : static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize - 1;
543 327685 : i >= 0; i--)
544 : {
545 : // Use memcpy() and not reinterpret_cast<GUInt16*> and
546 : // reinterpret_cast<float*>, otherwise compilers such as ICC
547 : // may (ab)use rules about aliasing to generate wrong code!
548 : GUInt16 usVal;
549 327680 : memcpy(&usVal, pabyTileData + i * sizeof(GUInt16),
550 : sizeof(usVal));
551 327680 : double dfVal =
552 327680 : (usVal * dfTileScale + dfTileOffset) * m_dfScale +
553 327680 : m_dfOffset;
554 327680 : if (m_dfPrecision == 1.0)
555 327680 : dfVal = floor(dfVal + 0.5);
556 327680 : if (bHasNoData && usVal == m_usGPKGNull)
557 196206 : dfVal = dfNoDataValue;
558 327680 : const float fVal = static_cast<float>(dfVal);
559 327680 : memcpy(pabyTileData + i * sizeof(float), &fVal, sizeof(fVal));
560 : }
561 : }
562 :
563 43 : return CE_None;
564 : }
565 :
566 2538 : GDALColorTable *poCT = nullptr;
567 2538 : if (nBands == 1 || nTileBandCount == 1)
568 : {
569 154 : poCT = poDSTile->GetRasterBand(1)->GetColorTable();
570 154 : IGetRasterBand(1)->GetColorTable();
571 : }
572 :
573 2538 : if (pbIsLossyFormat)
574 9 : *pbIsLossyFormat =
575 9 : !EQUAL(poDSTile->GetDriver()->GetDescription(), "PNG") ||
576 0 : (poCT != nullptr && poCT->GetColorEntryCount() == 256) /* PNG8 */;
577 :
578 : /* Map RGB(A) tile to single-band color indexed */
579 2538 : const GPtrDiff_t nBlockPixels =
580 2538 : static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize;
581 2538 : if (nBands == 1 && m_poCT != nullptr && nTileBandCount != 1)
582 : {
583 1 : std::map<GUInt32, int> oMapEntryToIndex;
584 1 : const int nEntries = std::min(256, m_poCT->GetColorEntryCount());
585 257 : for (int i = 0; i < nEntries; i++)
586 : {
587 256 : const GDALColorEntry *psEntry = m_poCT->GetColorEntry(i);
588 256 : GByte c1 = static_cast<GByte>(psEntry->c1);
589 256 : GByte c2 = static_cast<GByte>(psEntry->c2);
590 256 : GByte c3 = static_cast<GByte>(psEntry->c3);
591 256 : GUInt32 nVal = c1 + (c2 << 8) + (c3 << 16);
592 256 : if (nTileBandCount == 4)
593 256 : nVal += (static_cast<GUInt32>(psEntry->c4) << 24);
594 256 : oMapEntryToIndex[nVal] = i;
595 : }
596 : int iBestEntryFor0 =
597 1 : GPKGFindBestEntry(m_poCT, 0, 0, 0, 0, nTileBandCount);
598 10001 : for (GPtrDiff_t i = 0; i < nBlockPixels; i++)
599 : {
600 10000 : const GByte c1 = pabyTileData[i];
601 10000 : const GByte c2 = pabyTileData[i + nBlockPixels];
602 10000 : const GByte c3 = pabyTileData[i + 2 * nBlockPixels];
603 10000 : const GByte c4 = pabyTileData[i + 3 * nBlockPixels];
604 10000 : GUInt32 nVal = c1 + (c2 << 8) + (c3 << 16);
605 10000 : if (nTileBandCount == 4)
606 10000 : nVal += (c4 << 24);
607 10000 : if (nVal == 0)
608 : // In most cases we will reach that point at partial tiles.
609 5000 : pabyTileData[i] = static_cast<GByte>(iBestEntryFor0);
610 : else
611 : {
612 : std::map<GUInt32, int>::iterator oMapEntryToIndexIter =
613 5000 : oMapEntryToIndex.find(nVal);
614 5000 : if (oMapEntryToIndexIter == oMapEntryToIndex.end())
615 : /* Could happen with JPEG tiles */
616 0 : pabyTileData[i] = static_cast<GByte>(GPKGFindBestEntry(
617 : m_poCT, c1, c2, c3, c4, nTileBandCount));
618 : else
619 5000 : pabyTileData[i] =
620 5000 : static_cast<GByte>(oMapEntryToIndexIter->second);
621 : }
622 : }
623 1 : return CE_None;
624 : }
625 :
626 32 : if (nBands == 1 && nTileBandCount == 1 && poCT != nullptr &&
627 2569 : m_poCT != nullptr && !poCT->IsSame(m_poCT))
628 : {
629 0 : CPLError(CE_Warning, CPLE_NotSupported,
630 : "Different color tables. Unhandled for now");
631 : }
632 2537 : else if ((nBands == 1 && nTileBandCount >= 3) ||
633 32 : (nBands == 1 && nTileBandCount == 1 && m_poCT != nullptr &&
634 2537 : poCT == nullptr) ||
635 2537 : ((nBands == 1 || nBands == 2) && nTileBandCount == 1 &&
636 35 : m_poCT == nullptr && poCT != nullptr))
637 : {
638 0 : CPLError(CE_Failure, CPLE_AppDefined,
639 : "Inconsistent dataset and tiles band characteristics");
640 : }
641 :
642 2537 : if (nBands == 2)
643 : {
644 : // assuming that the RGB is Grey,Grey,Grey
645 307 : if (nTileBandCount == 1 || nTileBandCount == 3)
646 : {
647 : /* Create fully opaque alpha */
648 97 : memset(pabyTileData + 1 * nBlockPixels, 255, nBlockPixels);
649 : }
650 210 : else if (nTileBandCount == 4)
651 : {
652 : /* Transfer alpha band */
653 67 : memcpy(pabyTileData + 1 * nBlockPixels,
654 67 : pabyTileData + 3 * nBlockPixels, nBlockPixels);
655 : }
656 : }
657 2230 : else if (nTileBandCount == 2)
658 : {
659 : /* Do Grey+Alpha -> RGBA */
660 363 : memcpy(pabyTileData + 3 * nBlockPixels, pabyTileData + 1 * nBlockPixels,
661 : nBlockPixels);
662 363 : memcpy(pabyTileData + 1 * nBlockPixels, pabyTileData, nBlockPixels);
663 363 : memcpy(pabyTileData + 2 * nBlockPixels, pabyTileData, nBlockPixels);
664 : }
665 1867 : else if (nTileBandCount == 1 && !(nBands == 1 && m_poCT != nullptr))
666 : {
667 : /* Expand color indexed to RGB(A) */
668 122 : if (poCT != nullptr)
669 : {
670 : GByte abyCT[4 * 256];
671 15 : const int nEntries = std::min(256, poCT->GetColorEntryCount());
672 3600 : for (int i = 0; i < nEntries; i++)
673 : {
674 3585 : const GDALColorEntry *psEntry = poCT->GetColorEntry(i);
675 3585 : abyCT[4 * i] = static_cast<GByte>(psEntry->c1);
676 3585 : abyCT[4 * i + 1] = static_cast<GByte>(psEntry->c2);
677 3585 : abyCT[4 * i + 2] = static_cast<GByte>(psEntry->c3);
678 3585 : abyCT[4 * i + 3] = static_cast<GByte>(psEntry->c4);
679 : }
680 270 : for (int i = nEntries; i < 256; i++)
681 : {
682 255 : abyCT[4 * i] = 0;
683 255 : abyCT[4 * i + 1] = 0;
684 255 : abyCT[4 * i + 2] = 0;
685 255 : abyCT[4 * i + 3] = 0;
686 : }
687 967175 : for (GPtrDiff_t i = 0; i < nBlockPixels; i++)
688 : {
689 967160 : const GByte byVal = pabyTileData[i];
690 967160 : pabyTileData[i] = abyCT[4 * byVal];
691 967160 : pabyTileData[i + 1 * nBlockPixels] = abyCT[4 * byVal + 1];
692 967160 : pabyTileData[i + 2 * nBlockPixels] = abyCT[4 * byVal + 2];
693 967160 : pabyTileData[i + 3 * nBlockPixels] = abyCT[4 * byVal + 3];
694 : }
695 : }
696 : else
697 : {
698 107 : memcpy(pabyTileData + 1 * nBlockPixels, pabyTileData, nBlockPixels);
699 107 : memcpy(pabyTileData + 2 * nBlockPixels, pabyTileData, nBlockPixels);
700 107 : if (nBands == 4)
701 : {
702 96 : memset(pabyTileData + 3 * nBlockPixels, 255, nBlockPixels);
703 : }
704 122 : }
705 : }
706 1745 : else if (nTileBandCount == 3 && nBands == 4)
707 : {
708 : /* Create fully opaque alpha */
709 224 : memset(pabyTileData + 3 * nBlockPixels, 255, nBlockPixels);
710 : }
711 :
712 2537 : return CE_None;
713 : }
714 :
715 : /************************************************************************/
716 : /* ReadTile() */
717 : /************************************************************************/
718 :
719 7543 : GByte *GDALGPKGMBTilesLikePseudoDataset::ReadTile(int nRow, int nCol)
720 : {
721 : int nBlockXSize, nBlockYSize;
722 7543 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
723 7543 : const int nBands = IGetRasterCount();
724 7543 : const size_t nBandBlockSize =
725 7543 : static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
726 7543 : const int nTileBands = m_eDT == GDT_UInt8 ? 4 : 1;
727 7543 : if (m_nShiftXPixelsMod || m_nShiftYPixelsMod)
728 : {
729 4803 : GByte *pabyData = nullptr;
730 4803 : int i = 0;
731 11981 : for (; i < 4; i++)
732 : {
733 11981 : if (m_asCachedTilesDesc[i].nRow == nRow &&
734 7181 : m_asCachedTilesDesc[i].nCol == nCol)
735 : {
736 4803 : if (m_asCachedTilesDesc[i].nIdxWithinTileData >= 0)
737 : {
738 2040 : return m_pabyCachedTiles +
739 2040 : nBandBlockSize *
740 2040 : m_asCachedTilesDesc[i].nIdxWithinTileData *
741 2040 : nTileBands;
742 : }
743 : else
744 : {
745 2763 : if (i == 0)
746 194 : m_asCachedTilesDesc[i].nIdxWithinTileData =
747 194 : (m_asCachedTilesDesc[1].nIdxWithinTileData == 0)
748 194 : ? 1
749 : : 0;
750 2569 : else if (i == 1)
751 1189 : m_asCachedTilesDesc[i].nIdxWithinTileData =
752 1189 : (m_asCachedTilesDesc[0].nIdxWithinTileData == 0)
753 1189 : ? 1
754 : : 0;
755 1380 : else if (i == 2)
756 191 : m_asCachedTilesDesc[i].nIdxWithinTileData =
757 191 : (m_asCachedTilesDesc[3].nIdxWithinTileData == 2)
758 191 : ? 3
759 : : 2;
760 : else
761 1189 : m_asCachedTilesDesc[i].nIdxWithinTileData =
762 1189 : (m_asCachedTilesDesc[2].nIdxWithinTileData == 2)
763 1189 : ? 3
764 : : 2;
765 2763 : pabyData = m_pabyCachedTiles +
766 2763 : nBandBlockSize *
767 2763 : m_asCachedTilesDesc[i].nIdxWithinTileData *
768 2763 : nTileBands;
769 2763 : break;
770 : }
771 : }
772 : }
773 2763 : CPLAssert(i < 4);
774 2763 : return ReadTile(nRow, nCol, pabyData);
775 : }
776 : else
777 : {
778 2740 : GByte *pabyDest = m_pabyCachedTiles + 2 * nTileBands * nBandBlockSize;
779 2740 : bool bAllNonDirty = true;
780 11229 : for (int i = 0; i < nBands; i++)
781 : {
782 8489 : if (m_asCachedTilesDesc[0].abBandDirty[i])
783 2 : bAllNonDirty = false;
784 : }
785 2740 : if (bAllNonDirty)
786 : {
787 2738 : return ReadTile(nRow, nCol, pabyDest);
788 : }
789 :
790 : /* If some bands of the blocks are dirty/written we need to fetch */
791 : /* the tile in a temporary buffer in order not to override dirty bands*/
792 2 : GByte *pabyTemp = m_pabyCachedTiles + 3 * nTileBands * nBandBlockSize;
793 2 : if (ReadTile(nRow, nCol, pabyTemp) != nullptr)
794 : {
795 8 : for (int i = 0; i < nBands; i++)
796 : {
797 6 : if (!m_asCachedTilesDesc[0].abBandDirty[i])
798 : {
799 4 : memcpy(pabyDest + i * nBandBlockSize,
800 4 : pabyTemp + i * nBandBlockSize, nBandBlockSize);
801 : }
802 : }
803 : }
804 2 : return pabyDest;
805 : }
806 : }
807 :
808 : /************************************************************************/
809 : /* GetTileOffsetAndScale() */
810 : /************************************************************************/
811 :
812 2583 : void GDALGPKGMBTilesLikePseudoDataset::GetTileOffsetAndScale(
813 : GIntBig nTileId, double &dfTileOffset, double &dfTileScale)
814 : {
815 2583 : dfTileOffset = 0.0;
816 2583 : dfTileScale = 1.0;
817 :
818 2583 : if (m_eTF == GPKG_TF_PNG_16BIT)
819 : {
820 40 : char *pszSQL = sqlite3_mprintf(
821 : "SELECT offset, scale FROM gpkg_2d_gridded_tile_ancillary WHERE "
822 : "tpudt_name = '%q' AND tpudt_id = ?",
823 : m_osRasterTable.c_str());
824 40 : sqlite3_stmt *hStmt = nullptr;
825 40 : int rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt, nullptr);
826 40 : if (rc == SQLITE_OK)
827 : {
828 40 : sqlite3_bind_int64(hStmt, 1, nTileId);
829 40 : rc = sqlite3_step(hStmt);
830 40 : if (rc == SQLITE_ROW)
831 : {
832 40 : if (sqlite3_column_type(hStmt, 0) == SQLITE_FLOAT)
833 40 : dfTileOffset = sqlite3_column_double(hStmt, 0);
834 40 : if (sqlite3_column_type(hStmt, 1) == SQLITE_FLOAT)
835 40 : dfTileScale = sqlite3_column_double(hStmt, 1);
836 : }
837 40 : sqlite3_finalize(hStmt);
838 : }
839 40 : sqlite3_free(pszSQL);
840 : }
841 2583 : }
842 :
843 : /************************************************************************/
844 : /* ReadTile() */
845 : /************************************************************************/
846 :
847 5521 : GByte *GDALGPKGMBTilesLikePseudoDataset::ReadTile(int nRow, int nCol,
848 : GByte *pabyData,
849 : bool *pbIsLossyFormat)
850 : {
851 5521 : int nBlockXSize = 0;
852 5521 : int nBlockYSize = 0;
853 5521 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
854 5521 : const int nBands = IGetRasterCount();
855 :
856 5521 : if (pbIsLossyFormat)
857 18 : *pbIsLossyFormat = false;
858 :
859 5521 : const size_t nBandBlockSize =
860 5521 : static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
861 5521 : if (nRow < 0 || nCol < 0 || nRow >= m_nTileMatrixHeight ||
862 5464 : nCol >= m_nTileMatrixWidth)
863 : {
864 73 : FillEmptyTile(pabyData);
865 73 : return pabyData;
866 : }
867 :
868 : #ifdef DEBUG_VERBOSE
869 : CPLDebug("GPKG", "ReadTile(row=%d, col=%d)", nRow, nCol);
870 : #endif
871 :
872 16344 : char *pszSQL = sqlite3_mprintf(
873 : "SELECT tile_data%s FROM \"%w\" "
874 : "WHERE zoom_level = %d AND tile_row = %d AND tile_column = %d%s",
875 5448 : m_eDT != GDT_UInt8 ? ", id" : "", // MBTiles do not have an id
876 : m_osRasterTable.c_str(), m_nZoomLevel,
877 5448 : GetRowFromIntoTopConvention(nRow), nCol,
878 5448 : !m_osWHERE.empty() ? CPLSPrintf(" AND (%s)", m_osWHERE.c_str()) : "");
879 :
880 : #ifdef DEBUG_VERBOSE
881 : CPLDebug("GPKG", "%s", pszSQL);
882 : #endif
883 :
884 5448 : sqlite3_stmt *hStmt = nullptr;
885 5448 : int rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt, nullptr);
886 5448 : sqlite3_free(pszSQL);
887 5448 : if (rc != SQLITE_OK)
888 : {
889 0 : return nullptr;
890 : }
891 5448 : rc = sqlite3_step(hStmt);
892 :
893 5448 : if (rc == SQLITE_ROW && sqlite3_column_type(hStmt, 0) == SQLITE_BLOB)
894 : {
895 2566 : const int nBytes = sqlite3_column_bytes(hStmt, 0);
896 : GIntBig nTileId =
897 2566 : (m_eDT == GDT_UInt8) ? 0 : sqlite3_column_int64(hStmt, 1);
898 : GByte *pabyRawData = static_cast<GByte *>(
899 2566 : const_cast<void *>(sqlite3_column_blob(hStmt, 0)));
900 : const CPLString osMemFileName(
901 5132 : VSIMemGenerateHiddenFilename("gpkg_read_tile"));
902 2566 : VSILFILE *fp = VSIFileFromMemBuffer(osMemFileName.c_str(), pabyRawData,
903 : nBytes, FALSE);
904 2566 : VSIFCloseL(fp);
905 :
906 2566 : double dfTileOffset = 0.0;
907 2566 : double dfTileScale = 1.0;
908 2566 : GetTileOffsetAndScale(nTileId, dfTileOffset, dfTileScale);
909 2566 : ReadTile(osMemFileName, pabyData, dfTileOffset, dfTileScale,
910 : pbIsLossyFormat);
911 2566 : VSIUnlink(osMemFileName);
912 2566 : sqlite3_finalize(hStmt);
913 : }
914 2882 : else if (rc == SQLITE_BUSY)
915 : {
916 0 : FillEmptyTile(pabyData);
917 0 : CPLError(CE_Failure, CPLE_AppDefined,
918 : "sqlite3_step(%s) failed (SQLITE_BUSY): %s",
919 0 : sqlite3_sql(hStmt), sqlite3_errmsg(IGetDB()));
920 0 : sqlite3_finalize(hStmt);
921 0 : return pabyData;
922 : }
923 : else
924 : {
925 2882 : sqlite3_finalize(hStmt);
926 2882 : hStmt = nullptr;
927 :
928 2882 : if (m_hTempDB && (m_nShiftXPixelsMod || m_nShiftYPixelsMod))
929 : {
930 602 : const char *pszSQLNew = CPLSPrintf(
931 : "SELECT partial_flag, tile_data_band_1, tile_data_band_2, "
932 : "tile_data_band_3, tile_data_band_4 FROM partial_tiles WHERE "
933 : "zoom_level = %d AND tile_row = %d AND tile_column = %d",
934 : m_nZoomLevel, nRow, nCol);
935 :
936 : #ifdef DEBUG_VERBOSE
937 : CPLDebug("GPKG", "%s", pszSQLNew);
938 : #endif
939 :
940 602 : rc = SQLPrepareWithError(m_hTempDB, pszSQLNew, -1, &hStmt, nullptr);
941 602 : if (rc != SQLITE_OK)
942 : {
943 0 : FillEmptyTile(pabyData);
944 0 : return pabyData;
945 : }
946 :
947 602 : rc = sqlite3_step(hStmt);
948 602 : if (rc == SQLITE_ROW)
949 : {
950 2 : const int nPartialFlag = sqlite3_column_int(hStmt, 0);
951 10 : for (int iBand = 1; iBand <= nBands; iBand++)
952 : {
953 8 : GByte *pabyDestBand =
954 8 : pabyData + (iBand - 1) * nBandBlockSize;
955 8 : if (nPartialFlag & (((1 << 4) - 1) << (4 * (iBand - 1))))
956 : {
957 2 : CPLAssert(sqlite3_column_bytes(hStmt, iBand) ==
958 : static_cast<int>(nBandBlockSize));
959 2 : memcpy(pabyDestBand, sqlite3_column_blob(hStmt, iBand),
960 : nBandBlockSize);
961 : }
962 : else
963 : {
964 6 : FillEmptyTileSingleBand(pabyDestBand);
965 : }
966 : }
967 : }
968 : else
969 : {
970 600 : FillEmptyTile(pabyData);
971 : }
972 602 : sqlite3_finalize(hStmt);
973 : }
974 : else
975 : {
976 2280 : FillEmptyTile(pabyData);
977 : }
978 : }
979 :
980 5448 : return pabyData;
981 : }
982 :
983 : /************************************************************************/
984 : /* IReadBlock() */
985 : /************************************************************************/
986 :
987 3951 : CPLErr GDALGPKGMBTilesLikeRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
988 : void *pData)
989 : {
990 : #ifdef DEBUG_VERBOSE
991 : CPLDebug("GPKG",
992 : "IReadBlock(nBand=%d,nBlockXOff=%d,nBlockYOff=%d,m_nZoomLevel=%d)",
993 : nBand, nBlockXOff, nBlockYOff, m_poTPD->m_nZoomLevel);
994 : #endif
995 :
996 3951 : if (m_poTPD->m_pabyCachedTiles == nullptr)
997 0 : return CE_Failure;
998 :
999 3951 : const int nRowMin = nBlockYOff + m_poTPD->m_nShiftYTiles;
1000 3951 : int nRowMax = nRowMin;
1001 3951 : if (m_poTPD->m_nShiftYPixelsMod)
1002 1211 : nRowMax++;
1003 :
1004 3951 : const int nColMin = nBlockXOff + m_poTPD->m_nShiftXTiles;
1005 3951 : int nColMax = nColMin;
1006 3951 : if (m_poTPD->m_nShiftXPixelsMod)
1007 1189 : nColMax++;
1008 :
1009 2762 : retry:
1010 : /* Optimize for left to right reading at constant row */
1011 3954 : if (m_poTPD->m_nShiftXPixelsMod || m_poTPD->m_nShiftYPixelsMod)
1012 : {
1013 1214 : if (nRowMin == m_poTPD->m_asCachedTilesDesc[0].nRow &&
1014 1036 : nColMin == m_poTPD->m_asCachedTilesDesc[0].nCol + 1 &&
1015 1030 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData >= 0)
1016 : {
1017 1030 : CPLAssert(nRowMin == m_poTPD->m_asCachedTilesDesc[1].nRow);
1018 1030 : CPLAssert(nColMin == m_poTPD->m_asCachedTilesDesc[1].nCol);
1019 1030 : CPLAssert(m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData == 0 ||
1020 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData == 1);
1021 :
1022 : /* 0 1 --> 1 -1 */
1023 : /* 2 3 3 -1 */
1024 : /* or */
1025 : /* 1 0 --> 0 -1 */
1026 : /* 3 2 2 -1 */
1027 1030 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData =
1028 1030 : m_poTPD->m_asCachedTilesDesc[1].nIdxWithinTileData;
1029 1030 : m_poTPD->m_asCachedTilesDesc[2].nIdxWithinTileData =
1030 1030 : m_poTPD->m_asCachedTilesDesc[3].nIdxWithinTileData;
1031 : }
1032 : else
1033 : {
1034 184 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData = -1;
1035 184 : m_poTPD->m_asCachedTilesDesc[2].nIdxWithinTileData = -1;
1036 : }
1037 1214 : m_poTPD->m_asCachedTilesDesc[0].nRow = nRowMin;
1038 1214 : m_poTPD->m_asCachedTilesDesc[0].nCol = nColMin;
1039 1214 : m_poTPD->m_asCachedTilesDesc[1].nRow = nRowMin;
1040 1214 : m_poTPD->m_asCachedTilesDesc[1].nCol = nColMin + 1;
1041 1214 : m_poTPD->m_asCachedTilesDesc[2].nRow = nRowMin + 1;
1042 1214 : m_poTPD->m_asCachedTilesDesc[2].nCol = nColMin;
1043 1214 : m_poTPD->m_asCachedTilesDesc[3].nRow = nRowMin + 1;
1044 1214 : m_poTPD->m_asCachedTilesDesc[3].nCol = nColMin + 1;
1045 1214 : m_poTPD->m_asCachedTilesDesc[1].nIdxWithinTileData = -1;
1046 1214 : m_poTPD->m_asCachedTilesDesc[3].nIdxWithinTileData = -1;
1047 : }
1048 :
1049 9116 : for (int nRow = nRowMin; nRow <= nRowMax; nRow++)
1050 : {
1051 12705 : for (int nCol = nColMin; nCol <= nColMax; nCol++)
1052 : {
1053 7543 : if (m_poTPD->m_nShiftXPixelsMod == 0 &&
1054 2787 : m_poTPD->m_nShiftYPixelsMod == 0)
1055 : {
1056 2740 : if (!(nRow == m_poTPD->m_asCachedTilesDesc[0].nRow &&
1057 7 : nCol == m_poTPD->m_asCachedTilesDesc[0].nCol &&
1058 3 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData == 0))
1059 : {
1060 2738 : if (m_poTPD->WriteTile() != CE_None)
1061 0 : return CE_Failure;
1062 : }
1063 : }
1064 :
1065 7543 : GByte *pabyTileData = m_poTPD->ReadTile(nRow, nCol);
1066 7543 : if (pabyTileData == nullptr)
1067 0 : return CE_Failure;
1068 :
1069 34602 : for (int iBand = 1; iBand <= poDS->GetRasterCount(); iBand++)
1070 : {
1071 27062 : GDALRasterBlock *poBlock = nullptr;
1072 27062 : GByte *pabyDest = nullptr;
1073 27062 : if (iBand == nBand)
1074 : {
1075 7543 : pabyDest = static_cast<GByte *>(pData);
1076 : }
1077 : else
1078 : {
1079 19519 : poBlock = poDS->GetRasterBand(iBand)->GetLockedBlockRef(
1080 19519 : nBlockXOff, nBlockYOff, TRUE);
1081 19519 : if (poBlock == nullptr)
1082 0 : continue;
1083 19519 : if (poBlock->GetDirty())
1084 : {
1085 2 : poBlock->DropLock();
1086 2 : continue;
1087 : }
1088 : /* if we are short of GDAL cache max and there are dirty
1089 : * blocks */
1090 : /* of our dataset, the above GetLockedBlockRef() might have
1091 : * reset */
1092 : /* (at least part of) the 4 tiles we want to cache and have
1093 : */
1094 : /* already read */
1095 : // FIXME this is way too fragile.
1096 19517 : if ((m_poTPD->m_nShiftXPixelsMod != 0 ||
1097 5857 : m_poTPD->m_nShiftYPixelsMod != 0) &&
1098 13770 : (m_poTPD->m_asCachedTilesDesc[0].nRow != nRowMin ||
1099 13767 : m_poTPD->m_asCachedTilesDesc[0].nCol != nColMin))
1100 : {
1101 3 : poBlock->DropLock();
1102 3 : goto retry;
1103 : }
1104 19514 : pabyDest = static_cast<GByte *>(poBlock->GetDataRef());
1105 : }
1106 :
1107 : // Composite tile data into block data
1108 27057 : if (m_poTPD->m_nShiftXPixelsMod == 0 &&
1109 8641 : m_poTPD->m_nShiftYPixelsMod == 0)
1110 : {
1111 8487 : const size_t nBandBlockSize =
1112 8487 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
1113 8487 : m_nDTSize;
1114 8487 : memcpy(pabyDest,
1115 8487 : pabyTileData + (iBand - 1) * nBandBlockSize,
1116 8487 : nBandBlockSize);
1117 : #ifdef DEBUG_VERBOSE
1118 : if (eDataType == GDT_UInt8 &&
1119 : (nBlockXOff + 1) * nBlockXSize <= nRasterXSize &&
1120 : (nBlockYOff + 1) * nBlockYSize > nRasterYSize)
1121 : {
1122 : bool bFoundNonZero = false;
1123 : for (int y = nRasterYSize - nBlockYOff * nBlockYSize;
1124 : y < nBlockYSize; y++)
1125 : {
1126 : for (int x = 0; x < nBlockXSize; x++)
1127 : {
1128 : if (pabyDest[static_cast<GPtrDiff_t>(y) *
1129 : nBlockXSize +
1130 : x] != 0 &&
1131 : !bFoundNonZero)
1132 : {
1133 : CPLDebug("GPKG",
1134 : "IReadBlock(): Found non-zero "
1135 : "content in ghost part of "
1136 : "tile(nBand=%d,nBlockXOff=%d,"
1137 : "nBlockYOff=%d,m_nZoomLevel=%d)\n",
1138 : iBand, nBlockXOff, nBlockYOff,
1139 : m_poTPD->m_nZoomLevel);
1140 : bFoundNonZero = true;
1141 : }
1142 : }
1143 : }
1144 : }
1145 : #endif
1146 : }
1147 : else
1148 : {
1149 18570 : int nSrcXOffset = 0;
1150 18570 : int nSrcXSize = 0;
1151 18570 : int nDstXOffset = 0;
1152 18570 : if (nCol == nColMin)
1153 : {
1154 9362 : nSrcXOffset = m_poTPD->m_nShiftXPixelsMod;
1155 9362 : nSrcXSize = nBlockXSize - m_poTPD->m_nShiftXPixelsMod;
1156 9362 : nDstXOffset = 0;
1157 : }
1158 : else
1159 : {
1160 9208 : nSrcXOffset = 0;
1161 9208 : nSrcXSize = m_poTPD->m_nShiftXPixelsMod;
1162 9208 : nDstXOffset = nBlockXSize - m_poTPD->m_nShiftXPixelsMod;
1163 : }
1164 18570 : int nSrcYOffset = 0;
1165 18570 : int nSrcYSize = 0;
1166 18570 : int nDstYOffset = 0;
1167 18570 : if (nRow == nRowMin)
1168 : {
1169 9288 : nSrcYOffset = m_poTPD->m_nShiftYPixelsMod;
1170 9288 : nSrcYSize = nBlockYSize - m_poTPD->m_nShiftYPixelsMod;
1171 9288 : nDstYOffset = 0;
1172 : }
1173 : else
1174 : {
1175 9282 : nSrcYOffset = 0;
1176 9282 : nSrcYSize = m_poTPD->m_nShiftYPixelsMod;
1177 9282 : nDstYOffset = nBlockYSize - m_poTPD->m_nShiftYPixelsMod;
1178 : }
1179 :
1180 : #ifdef DEBUG_VERBOSE
1181 : CPLDebug("GPKG",
1182 : "Copy source tile x=%d,w=%d,y=%d,h=%d into "
1183 : "buffer at x=%d,y=%d",
1184 : nSrcXOffset, nSrcXSize, nSrcYOffset, nSrcYSize,
1185 : nDstXOffset, nDstYOffset);
1186 : #endif
1187 :
1188 2368610 : for (GPtrDiff_t y = 0; y < nSrcYSize; y++)
1189 : {
1190 2350040 : GByte *pSrc =
1191 : pabyTileData +
1192 2350040 : (static_cast<GPtrDiff_t>(iBand - 1) * nBlockXSize *
1193 2350040 : nBlockYSize +
1194 2350040 : (y + nSrcYOffset) * nBlockXSize + nSrcXOffset) *
1195 2350040 : m_nDTSize;
1196 2350040 : GByte *pDst =
1197 : pabyDest +
1198 2350040 : ((y + nDstYOffset) * nBlockXSize + nDstXOffset) *
1199 2350040 : m_nDTSize;
1200 2350040 : GDALCopyWords(pSrc, eDataType, m_nDTSize, pDst,
1201 : eDataType, m_nDTSize, nSrcXSize);
1202 : }
1203 : }
1204 :
1205 27057 : if (poBlock)
1206 19514 : poBlock->DropLock();
1207 : }
1208 : }
1209 : }
1210 :
1211 3951 : return CE_None;
1212 : }
1213 :
1214 : /************************************************************************/
1215 : /* IGetDataCoverageStatus() */
1216 : /************************************************************************/
1217 :
1218 8 : int GDALGPKGMBTilesLikeRasterBand::IGetDataCoverageStatus(int nXOff, int nYOff,
1219 : int nXSize,
1220 : int nYSize,
1221 : int nMaskFlagStop,
1222 : double *pdfDataPct)
1223 : {
1224 8 : if (eAccess == GA_Update)
1225 0 : FlushCache(false);
1226 :
1227 8 : const int iColMin = nXOff / nBlockXSize + m_poTPD->m_nShiftXTiles;
1228 16 : const int iColMax = (nXOff + nXSize - 1) / nBlockXSize +
1229 8 : m_poTPD->m_nShiftXTiles +
1230 8 : (m_poTPD->m_nShiftXPixelsMod ? 1 : 0);
1231 8 : const int iRowMin = nYOff / nBlockYSize + m_poTPD->m_nShiftYTiles;
1232 16 : const int iRowMax = (nYOff + nYSize - 1) / nBlockYSize +
1233 8 : m_poTPD->m_nShiftYTiles +
1234 8 : (m_poTPD->m_nShiftYPixelsMod ? 1 : 0);
1235 :
1236 8 : int iDBRowMin = m_poTPD->GetRowFromIntoTopConvention(iRowMin);
1237 8 : int iDBRowMax = m_poTPD->GetRowFromIntoTopConvention(iRowMax);
1238 8 : if (iDBRowMin > iDBRowMax)
1239 0 : std::swap(iDBRowMin, iDBRowMax);
1240 :
1241 16 : char *pszSQL = sqlite3_mprintf(
1242 : "SELECT tile_column, tile_row FROM \"%w\" "
1243 : "WHERE zoom_level = %d AND "
1244 : "(tile_row BETWEEN %d AND %d) AND "
1245 : "(tile_column BETWEEN %d AND %d)"
1246 : "%s",
1247 8 : m_poTPD->m_osRasterTable.c_str(), m_poTPD->m_nZoomLevel, iDBRowMin,
1248 : iDBRowMax, iColMin, iColMax,
1249 8 : !m_poTPD->m_osWHERE.empty()
1250 0 : ? CPLSPrintf(" AND (%s)", m_poTPD->m_osWHERE.c_str())
1251 : : "");
1252 :
1253 : #ifdef DEBUG_VERBOSE
1254 : CPLDebug("GPKG", "%s", pszSQL);
1255 : #endif
1256 :
1257 8 : sqlite3_stmt *hStmt = nullptr;
1258 : int rc =
1259 8 : SQLPrepareWithError(m_poTPD->IGetDB(), pszSQL, -1, &hStmt, nullptr);
1260 8 : if (rc != SQLITE_OK)
1261 : {
1262 0 : sqlite3_free(pszSQL);
1263 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
1264 0 : GDAL_DATA_COVERAGE_STATUS_DATA;
1265 : }
1266 8 : sqlite3_free(pszSQL);
1267 8 : rc = sqlite3_step(hStmt);
1268 16 : std::set<std::pair<int, int>> oSetTiles; // (col, row)
1269 12 : while (rc == SQLITE_ROW)
1270 : {
1271 4 : const int iCol = sqlite3_column_int(hStmt, 0);
1272 : const int iRow =
1273 4 : m_poTPD->GetRowFromIntoTopConvention(sqlite3_column_int(hStmt, 1));
1274 4 : oSetTiles.insert(std::pair(iCol, iRow));
1275 4 : rc = sqlite3_step(hStmt);
1276 : }
1277 8 : sqlite3_finalize(hStmt);
1278 8 : if (rc != SQLITE_DONE)
1279 : {
1280 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
1281 0 : GDAL_DATA_COVERAGE_STATUS_DATA;
1282 : }
1283 8 : if (oSetTiles.empty())
1284 : {
1285 4 : if (pdfDataPct)
1286 4 : *pdfDataPct = 0.0;
1287 4 : return GDAL_DATA_COVERAGE_STATUS_EMPTY;
1288 : }
1289 :
1290 8 : if (m_poTPD->m_nShiftXPixelsMod == 0 && m_poTPD->m_nShiftYPixelsMod == 0 &&
1291 4 : oSetTiles.size() == static_cast<size_t>(iRowMax - iRowMin + 1) *
1292 4 : (iColMax - iColMin + 1))
1293 : {
1294 2 : if (pdfDataPct)
1295 2 : *pdfDataPct = 100.0;
1296 2 : return GDAL_DATA_COVERAGE_STATUS_DATA;
1297 : }
1298 :
1299 2 : if (m_poTPD->m_nShiftXPixelsMod == 0 && m_poTPD->m_nShiftYPixelsMod == 0)
1300 : {
1301 2 : int nStatus = 0;
1302 2 : GIntBig nPixelsData = 0;
1303 7 : for (int iY = iRowMin; iY <= iRowMax; ++iY)
1304 : {
1305 21 : for (int iX = iColMin; iX <= iColMax; ++iX)
1306 : {
1307 16 : if (oSetTiles.find(std::pair(iX, iY)) == oSetTiles.end())
1308 : {
1309 14 : nStatus |= GDAL_DATA_COVERAGE_STATUS_EMPTY;
1310 : }
1311 : else
1312 : {
1313 2 : const int iXGDAL = iX - m_poTPD->m_nShiftXTiles;
1314 2 : const int iYGDAL = iY - m_poTPD->m_nShiftYTiles;
1315 2 : const int nXBlockRight =
1316 2 : (iXGDAL * nBlockXSize > INT_MAX - nBlockXSize)
1317 2 : ? INT_MAX
1318 2 : : (iXGDAL + 1) * nBlockXSize;
1319 2 : const int nYBlockBottom =
1320 2 : (iYGDAL * nBlockYSize > INT_MAX - nBlockYSize)
1321 2 : ? INT_MAX
1322 2 : : (iYGDAL + 1) * nBlockYSize;
1323 :
1324 4 : nPixelsData += (static_cast<GIntBig>(std::min(
1325 2 : nXBlockRight, nXOff + nXSize)) -
1326 2 : std::max(iXGDAL * nBlockXSize, nXOff)) *
1327 2 : (std::min(nYBlockBottom, nYOff + nYSize) -
1328 2 : std::max(iYGDAL * nBlockYSize, nYOff));
1329 2 : nStatus |= GDAL_DATA_COVERAGE_STATUS_DATA;
1330 : }
1331 16 : if (nMaskFlagStop != 0 && (nMaskFlagStop & nStatus) != 0)
1332 : {
1333 0 : if (pdfDataPct)
1334 0 : *pdfDataPct = -1.0;
1335 0 : return nStatus;
1336 : }
1337 : }
1338 : }
1339 :
1340 2 : if (pdfDataPct)
1341 : {
1342 2 : *pdfDataPct =
1343 2 : 100.0 * nPixelsData / (static_cast<GIntBig>(nXSize) * nYSize);
1344 : }
1345 2 : return nStatus;
1346 : }
1347 :
1348 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
1349 0 : GDAL_DATA_COVERAGE_STATUS_DATA;
1350 : }
1351 :
1352 : /************************************************************************/
1353 : /* WEBPSupports4Bands() */
1354 : /************************************************************************/
1355 :
1356 83 : static bool WEBPSupports4Bands()
1357 : {
1358 : static int bRes = -1;
1359 83 : if (bRes < 0)
1360 : {
1361 1 : GDALDriver *poDrv = GDALDriver::FromHandle(GDALGetDriverByName("WEBP"));
1362 2 : if (poDrv == nullptr ||
1363 1 : CPLTestBool(CPLGetConfigOption("GPKG_SIMUL_WEBP_3BAND", "FALSE")))
1364 0 : bRes = false;
1365 : else
1366 : {
1367 : // LOSSLESS and RGBA support appeared in the same version
1368 1 : bRes = strstr(poDrv->GetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST),
1369 1 : "LOSSLESS") != nullptr;
1370 : }
1371 1 : if (poDrv != nullptr && !bRes)
1372 : {
1373 0 : CPLError(
1374 : CE_Warning, CPLE_AppDefined,
1375 : "The version of WEBP available does not support 4-band RGBA");
1376 : }
1377 : }
1378 83 : return CPL_TO_BOOL(bRes);
1379 : }
1380 :
1381 : /************************************************************************/
1382 : /* GetTileId() */
1383 : /************************************************************************/
1384 :
1385 58 : GIntBig GDALGPKGMBTilesLikePseudoDataset::GetTileId(int nRow, int nCol)
1386 : {
1387 : char *pszSQL =
1388 58 : sqlite3_mprintf("SELECT id FROM \"%w\" WHERE zoom_level = %d AND "
1389 : "tile_row = %d AND tile_column = %d",
1390 : m_osRasterTable.c_str(), m_nZoomLevel,
1391 58 : GetRowFromIntoTopConvention(nRow), nCol);
1392 58 : GIntBig nRes = SQLGetInteger64(IGetDB(), pszSQL, nullptr);
1393 58 : sqlite3_free(pszSQL);
1394 58 : return nRes;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* DeleteTile() */
1399 : /************************************************************************/
1400 :
1401 460 : bool GDALGPKGMBTilesLikePseudoDataset::DeleteTile(int nRow, int nCol)
1402 : {
1403 : char *pszSQL =
1404 460 : sqlite3_mprintf("DELETE FROM \"%w\" "
1405 : "WHERE zoom_level = %d AND tile_row = %d AND "
1406 : "tile_column = %d",
1407 : m_osRasterTable.c_str(), m_nZoomLevel,
1408 460 : GetRowFromIntoTopConvention(nRow), nCol);
1409 : #ifdef DEBUG_VERBOSE
1410 : CPLDebug("GPKG", "%s", pszSQL);
1411 : #endif
1412 460 : char *pszErrMsg = nullptr;
1413 460 : int rc = sqlite3_exec(IGetDB(), pszSQL, nullptr, nullptr, &pszErrMsg);
1414 460 : if (rc != SQLITE_OK)
1415 : {
1416 0 : CPLError(CE_Failure, CPLE_AppDefined,
1417 : "Failure when deleting tile (row=%d,col=%d) "
1418 : "at zoom_level=%d : %s",
1419 0 : GetRowFromIntoTopConvention(nRow), nCol, m_nZoomLevel,
1420 0 : pszErrMsg ? pszErrMsg : "");
1421 : }
1422 460 : sqlite3_free(pszSQL);
1423 460 : sqlite3_free(pszErrMsg);
1424 460 : return rc == SQLITE_OK;
1425 : }
1426 :
1427 : /************************************************************************/
1428 : /* DeleteFromGriddedTileAncillary() */
1429 : /************************************************************************/
1430 :
1431 58 : bool GDALGPKGMBTilesLikePseudoDataset::DeleteFromGriddedTileAncillary(
1432 : GIntBig nTileId)
1433 : {
1434 : char *pszSQL =
1435 58 : sqlite3_mprintf("DELETE FROM gpkg_2d_gridded_tile_ancillary WHERE "
1436 : "tpudt_name = '%q' AND tpudt_id = ?",
1437 : m_osRasterTable.c_str());
1438 58 : sqlite3_stmt *hStmt = nullptr;
1439 58 : int rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt, nullptr);
1440 58 : if (rc == SQLITE_OK)
1441 : {
1442 58 : sqlite3_bind_int64(hStmt, 1, nTileId);
1443 58 : rc = sqlite3_step(hStmt);
1444 58 : sqlite3_finalize(hStmt);
1445 : }
1446 58 : sqlite3_free(pszSQL);
1447 58 : return rc == SQLITE_OK;
1448 : }
1449 :
1450 : /************************************************************************/
1451 : /* ProcessInt16UInt16Tile() */
1452 : /************************************************************************/
1453 :
1454 : template <class T>
1455 48 : static void ProcessInt16UInt16Tile(
1456 : const void *pabyData, GPtrDiff_t nPixels, bool bIsInt16, bool bHasNoData,
1457 : double dfNoDataValue, GUInt16 usGPKGNull, double m_dfOffset,
1458 : double m_dfScale, GUInt16 *pTempTileBuffer, double &dfTileOffset,
1459 : double &dfTileScale, double &dfTileMin, double &dfTileMax,
1460 : double &dfTileMean, double &dfTileStdDev, GPtrDiff_t &nValidPixels)
1461 : {
1462 48 : const T *pSrc = reinterpret_cast<const T *>(pabyData);
1463 48 : T nMin = 0;
1464 48 : T nMax = 0;
1465 48 : double dfM2 = 0.0;
1466 3145777 : for (int i = 0; i < nPixels; i++)
1467 : {
1468 3145726 : const T nVal = pSrc[i];
1469 3145726 : if (bHasNoData && nVal == dfNoDataValue)
1470 718161 : continue;
1471 :
1472 2427565 : if (nValidPixels == 0)
1473 : {
1474 47 : nMin = nVal;
1475 47 : nMax = nVal;
1476 : }
1477 : else
1478 : {
1479 2427515 : nMin = std::min(nMin, nVal);
1480 2427515 : nMax = std::max(nMax, nVal);
1481 : }
1482 2427565 : nValidPixels++;
1483 2427565 : const double dfDelta = nVal - dfTileMean;
1484 2427565 : dfTileMean += dfDelta / nValidPixels;
1485 2427565 : dfM2 += dfDelta * (nVal - dfTileMean);
1486 : }
1487 48 : dfTileMin = nMin;
1488 48 : dfTileMax = nMax;
1489 48 : if (nValidPixels)
1490 47 : dfTileStdDev = sqrt(dfM2 / nValidPixels);
1491 :
1492 48 : double dfGlobalMin = (nMin - m_dfOffset) / m_dfScale;
1493 48 : double dfGlobalMax = (nMax - m_dfOffset) / m_dfScale;
1494 48 : double dfRange = 65535.0;
1495 48 : if (bHasNoData && usGPKGNull == 65535 &&
1496 8 : dfGlobalMax - dfGlobalMin >= dfRange)
1497 : {
1498 0 : dfRange = 65534.0;
1499 : }
1500 :
1501 48 : if (dfGlobalMax - dfGlobalMin > dfRange)
1502 : {
1503 0 : dfTileScale = (dfGlobalMax - dfGlobalMin) / dfRange;
1504 : }
1505 48 : if (dfGlobalMin < 0.0)
1506 : {
1507 0 : dfTileOffset = -dfGlobalMin;
1508 : }
1509 48 : else if (dfGlobalMax / dfTileScale > dfRange)
1510 : {
1511 0 : dfTileOffset = dfGlobalMax - dfRange * dfTileScale;
1512 : }
1513 :
1514 48 : if (bHasNoData && std::numeric_limits<T>::min() == 0 && m_dfOffset == 0.0 &&
1515 : m_dfScale == 1.0)
1516 : {
1517 7 : dfTileOffset = 0.0;
1518 7 : dfTileScale = 1.0;
1519 : }
1520 41 : else if (bHasNoData && bIsInt16 && dfNoDataValue == -32768.0 &&
1521 1 : usGPKGNull == 65535 && m_dfOffset == -32768.0 && m_dfScale == 1.0)
1522 : {
1523 1 : dfTileOffset = 1.0;
1524 1 : dfTileScale = 1.0;
1525 : }
1526 :
1527 3145777 : for (GPtrDiff_t i = 0; i < nPixels; i++)
1528 : {
1529 3145726 : const T nVal = pSrc[i];
1530 3145726 : if (bHasNoData && nVal == dfNoDataValue)
1531 718161 : pTempTileBuffer[i] = usGPKGNull;
1532 : else
1533 : {
1534 2427565 : double dfVal =
1535 2427565 : ((nVal - m_dfOffset) / m_dfScale - dfTileOffset) / dfTileScale;
1536 2427565 : CPLAssert(dfVal >= 0.0 && dfVal < 65535.5);
1537 2427565 : pTempTileBuffer[i] = static_cast<GUInt16>(dfVal + 0.5);
1538 2427565 : if (bHasNoData && pTempTileBuffer[i] == usGPKGNull)
1539 : {
1540 0 : if (usGPKGNull > 0)
1541 0 : pTempTileBuffer[i]--;
1542 : else
1543 0 : pTempTileBuffer[i]++;
1544 : ;
1545 : }
1546 : }
1547 : }
1548 48 : }
1549 :
1550 : /************************************************************************/
1551 : /* WriteTile() */
1552 : /************************************************************************/
1553 :
1554 13633 : CPLErr GDALGPKGMBTilesLikePseudoDataset::WriteTile()
1555 : {
1556 13633 : GDALGPKGMBTilesLikePseudoDataset *poMainDS =
1557 13633 : m_poParentDS ? m_poParentDS : this;
1558 13633 : if (poMainDS->m_nTileInsertionCount < 0)
1559 501 : return CE_Failure;
1560 :
1561 13132 : if (m_bInWriteTile)
1562 : {
1563 : // Shouldn't happen in practice, but #7022 shows that the unexpected
1564 : // can happen sometimes.
1565 0 : CPLError(
1566 : CE_Failure, CPLE_AppDefined,
1567 : "Recursive call to GDALGPKGMBTilesLikePseudoDataset::WriteTile()");
1568 0 : return CE_Failure;
1569 : }
1570 13132 : GDALRasterBlock::EnterDisableDirtyBlockFlush();
1571 13132 : m_bInWriteTile = true;
1572 13132 : CPLErr eErr = WriteTileInternal();
1573 13132 : m_bInWriteTile = false;
1574 13132 : GDALRasterBlock::LeaveDisableDirtyBlockFlush();
1575 13132 : return eErr;
1576 : }
1577 :
1578 : /* should only be called by WriteTile() */
1579 13132 : CPLErr GDALGPKGMBTilesLikePseudoDataset::WriteTileInternal()
1580 : {
1581 17758 : if (!(IGetUpdate() && m_asCachedTilesDesc[0].nRow >= 0 &&
1582 4626 : m_asCachedTilesDesc[0].nCol >= 0 &&
1583 4626 : m_asCachedTilesDesc[0].nIdxWithinTileData == 0))
1584 8506 : return CE_None;
1585 :
1586 4626 : int nRow = m_asCachedTilesDesc[0].nRow;
1587 4626 : int nCol = m_asCachedTilesDesc[0].nCol;
1588 :
1589 4626 : bool bAllDirty = true;
1590 4626 : bool bAllNonDirty = true;
1591 4626 : const int nBands = IGetRasterCount();
1592 18010 : for (int i = 0; i < nBands; i++)
1593 : {
1594 13384 : if (m_asCachedTilesDesc[0].abBandDirty[i])
1595 13351 : bAllNonDirty = false;
1596 : else
1597 33 : bAllDirty = false;
1598 : }
1599 4626 : if (bAllNonDirty)
1600 0 : return CE_None;
1601 :
1602 : int nBlockXSize, nBlockYSize;
1603 4626 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
1604 :
1605 : /* If all bands for that block are not dirty/written, we need to */
1606 : /* fetch the missing ones if the tile exists */
1607 4626 : bool bIsLossyFormat = false;
1608 4626 : const size_t nBandBlockSize =
1609 4626 : static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
1610 4626 : if (!bAllDirty)
1611 : {
1612 72 : for (int i = 1; i <= 3; i++)
1613 : {
1614 54 : m_asCachedTilesDesc[i].nRow = -1;
1615 54 : m_asCachedTilesDesc[i].nCol = -1;
1616 54 : m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
1617 : }
1618 18 : const int nTileBands = m_eDT == GDT_UInt8 ? 4 : 1;
1619 18 : GByte *pabyTemp = m_pabyCachedTiles + nTileBands * nBandBlockSize;
1620 18 : ReadTile(nRow, nCol, pabyTemp, &bIsLossyFormat);
1621 82 : for (int i = 0; i < nBands; i++)
1622 : {
1623 64 : if (!m_asCachedTilesDesc[0].abBandDirty[i])
1624 : {
1625 33 : memcpy(m_pabyCachedTiles + i * nBandBlockSize,
1626 33 : pabyTemp + i * nBandBlockSize, nBandBlockSize);
1627 : }
1628 : }
1629 : }
1630 :
1631 : /* Compute origin of tile in GDAL raster space */
1632 4626 : int nXOff = (nCol - m_nShiftXTiles) * nBlockXSize - m_nShiftXPixelsMod;
1633 4626 : int nYOff = (nRow - m_nShiftYTiles) * nBlockYSize - m_nShiftYPixelsMod;
1634 :
1635 : /* Assert that the tile at least intersects some of the GDAL raster space */
1636 4626 : CPLAssert(nXOff > -nBlockXSize);
1637 4626 : CPLAssert(nYOff > -nBlockYSize);
1638 : /* Can happen if the tile of the raster is less than the block size */
1639 4626 : const int nRasterXSize = IGetRasterBand(1)->GetXSize();
1640 4626 : const int nRasterYSize = IGetRasterBand(1)->GetYSize();
1641 4626 : if (nXOff >= nRasterXSize || nYOff >= nRasterYSize)
1642 0 : return CE_None;
1643 :
1644 : #ifdef DEBUG_VERBOSE
1645 : if (m_nShiftXPixelsMod == 0 && m_nShiftYPixelsMod == 0 &&
1646 : m_eDT == GDT_UInt8)
1647 : {
1648 : int nBlockXOff = nCol;
1649 : int nBlockYOff = nRow;
1650 : if (nBlockXOff * nBlockXSize <= nRasterXSize - nBlockXSize &&
1651 : nBlockYOff * nBlockYSize > nRasterYSize - nBlockYSize)
1652 : {
1653 : for (int i = 0; i < nBands; i++)
1654 : {
1655 : bool bFoundNonZero = false;
1656 : for (GPtrDiff_t y = nRasterYSize - nBlockYOff * nBlockYSize;
1657 : y < nBlockYSize; y++)
1658 : {
1659 : for (int x = 0; x < nBlockXSize; x++)
1660 : {
1661 : if (m_pabyCachedTiles[y * nBlockXSize + x +
1662 : i * nBandBlockSize] != 0 &&
1663 : !bFoundNonZero)
1664 : {
1665 : CPLDebug("GPKG",
1666 : "WriteTileInternal(): Found non-zero "
1667 : "content in ghost part of "
1668 : "tile(band=%d,nBlockXOff=%d,nBlockYOff=%d,"
1669 : "m_nZoomLevel=%d)\n",
1670 : i + 1, nBlockXOff, nBlockYOff,
1671 : m_nZoomLevel);
1672 : bFoundNonZero = true;
1673 : }
1674 : }
1675 : }
1676 : }
1677 : }
1678 : }
1679 : #endif
1680 :
1681 : /* Validity area of tile data in intra-tile coordinate space */
1682 4626 : int iXOff = 0;
1683 4626 : GPtrDiff_t iYOff = 0;
1684 4626 : int iXCount = nBlockXSize;
1685 4626 : int iYCount = nBlockYSize;
1686 :
1687 4626 : bool bPartialTile = false;
1688 4626 : int nAlphaBand = (nBands == 2) ? 2 : (nBands == 4) ? 4 : 0;
1689 4626 : if (nAlphaBand == 0)
1690 : {
1691 3744 : if (nXOff < 0)
1692 : {
1693 4 : bPartialTile = true;
1694 4 : iXOff = -nXOff;
1695 4 : iXCount += nXOff;
1696 : }
1697 3744 : if (nXOff > nRasterXSize - nBlockXSize)
1698 : {
1699 95 : bPartialTile = true;
1700 95 : iXCount -= static_cast<int>(static_cast<GIntBig>(nXOff) +
1701 95 : nBlockXSize - nRasterXSize);
1702 : }
1703 3744 : if (nYOff < 0)
1704 : {
1705 8 : bPartialTile = true;
1706 8 : iYOff = -nYOff;
1707 8 : iYCount += nYOff;
1708 : }
1709 3744 : if (nYOff > nRasterYSize - nBlockYSize)
1710 : {
1711 109 : bPartialTile = true;
1712 109 : iYCount -= static_cast<int>(static_cast<GIntBig>(nYOff) +
1713 109 : nBlockYSize - nRasterYSize);
1714 : }
1715 3744 : CPLAssert(iXOff >= 0);
1716 3744 : CPLAssert(iYOff >= 0);
1717 3744 : CPLAssert(iXCount > 0);
1718 3744 : CPLAssert(iYCount > 0);
1719 3744 : CPLAssert(iXOff + iXCount <= nBlockXSize);
1720 3744 : CPLAssert(iYOff + iYCount <= nBlockYSize);
1721 : }
1722 :
1723 4626 : m_asCachedTilesDesc[0].nRow = -1;
1724 4626 : m_asCachedTilesDesc[0].nCol = -1;
1725 4626 : m_asCachedTilesDesc[0].nIdxWithinTileData = -1;
1726 4626 : m_asCachedTilesDesc[0].abBandDirty[0] = false;
1727 4626 : m_asCachedTilesDesc[0].abBandDirty[1] = false;
1728 4626 : m_asCachedTilesDesc[0].abBandDirty[2] = false;
1729 4626 : m_asCachedTilesDesc[0].abBandDirty[3] = false;
1730 :
1731 4626 : CPLErr eErr = CE_Failure;
1732 :
1733 4626 : int bHasNoData = FALSE;
1734 4626 : double dfNoDataValue = IGetRasterBand(1)->GetNoDataValue(&bHasNoData);
1735 4626 : const bool bHasNanNoData = bHasNoData && std::isnan(dfNoDataValue);
1736 :
1737 4626 : bool bAllOpaque = true;
1738 : // Detect fully transparent tiles, but only if all bands are dirty (that is
1739 : // the user wrote content into them) !
1740 4626 : if (bAllDirty && m_eDT == GDT_UInt8 && m_poCT == nullptr && nAlphaBand != 0)
1741 : {
1742 872 : GByte byFirstAlphaVal =
1743 872 : m_pabyCachedTiles[(nAlphaBand - 1) * nBlockXSize * nBlockYSize];
1744 872 : GPtrDiff_t i = 1;
1745 24105900 : for (; i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize; i++)
1746 : {
1747 24105400 : if (m_pabyCachedTiles[static_cast<GPtrDiff_t>(nAlphaBand - 1) *
1748 24105400 : nBlockXSize * nBlockYSize +
1749 24105400 : i] != byFirstAlphaVal)
1750 367 : break;
1751 : }
1752 872 : if (i == static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize)
1753 : {
1754 : // If tile is fully transparent, don't serialize it and remove it if
1755 : // it exists
1756 505 : if (byFirstAlphaVal == 0)
1757 : {
1758 440 : DeleteTile(nRow, nCol);
1759 :
1760 440 : return CE_None;
1761 : }
1762 65 : bAllOpaque = (byFirstAlphaVal == 255);
1763 : }
1764 : else
1765 432 : bAllOpaque = false;
1766 : }
1767 3754 : else if (bAllDirty && m_eDT == GDT_UInt8 && m_poCT == nullptr &&
1768 3656 : (!bHasNoData || dfNoDataValue == 0.0))
1769 : {
1770 3656 : bool bAllEmpty = true;
1771 3656 : const auto nPixels =
1772 3656 : static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize * nBands;
1773 1207860 : for (GPtrDiff_t i = 0; i < nPixels; i++)
1774 : {
1775 1207850 : if (m_pabyCachedTiles[i] != 0)
1776 : {
1777 3639 : bAllEmpty = false;
1778 3639 : break;
1779 : }
1780 : }
1781 3656 : if (bAllEmpty)
1782 : {
1783 : // If tile is fully transparent, don't serialize it and remove it if
1784 : // it exists
1785 17 : DeleteTile(nRow, nCol);
1786 :
1787 17 : return CE_None;
1788 3639 : }
1789 : }
1790 98 : else if (bAllDirty && m_eDT == GDT_Float32)
1791 : {
1792 12 : const float *pSrc = reinterpret_cast<float *>(m_pabyCachedTiles);
1793 : GPtrDiff_t i;
1794 12 : const float fNoDataValueOrZero =
1795 12 : bHasNoData ? static_cast<float>(dfNoDataValue) : 0.0f;
1796 131084 : for (i = 0; i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize; i++)
1797 : {
1798 131082 : const float fVal = pSrc[i];
1799 131082 : if (bHasNanNoData)
1800 : {
1801 0 : if (std::isnan(fVal))
1802 0 : continue;
1803 : }
1804 131082 : else if (fVal == fNoDataValueOrZero)
1805 : {
1806 131072 : continue;
1807 : }
1808 10 : break;
1809 : }
1810 12 : if (i == static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize)
1811 : {
1812 : // If tile is fully transparent, don't serialize it and remove it if
1813 : // it exists
1814 2 : DeleteTile(nRow, nCol);
1815 :
1816 2 : return CE_None;
1817 : }
1818 : }
1819 :
1820 4167 : if (bIsLossyFormat)
1821 : {
1822 1 : CPLDebug("GPKG",
1823 : "Had to read tile (row=%d,col=%d) at zoom_level=%d, "
1824 : "stored in a lossy format, before rewriting it, causing "
1825 : "potential extra quality loss",
1826 : nRow, nCol, m_nZoomLevel);
1827 : }
1828 :
1829 : const CPLString osMemFileName(
1830 8334 : VSIMemGenerateHiddenFilename("gpkg_write_tile"));
1831 4167 : const char *pszDriverName = "PNG";
1832 4167 : CPL_IGNORE_RET_VAL(pszDriverName); // Make CSA happy
1833 4167 : bool bTileDriverSupports1Band = false;
1834 4167 : bool bTileDriverSupports2Bands = false;
1835 4167 : bool bTileDriverSupports4Bands = false;
1836 4167 : bool bTileDriverSupportsCT = false;
1837 :
1838 4167 : if (nBands == 1 && m_eDT == GDT_UInt8)
1839 78 : IGetRasterBand(1)->GetColorTable();
1840 :
1841 4167 : GDALDataType eTileDT = GDT_UInt8;
1842 : // If not all bands are dirty, then (temporarily) use a lossless format
1843 4167 : if (m_eTF == GPKG_TF_PNG_JPEG || (!bAllDirty && m_eTF == GPKG_TF_JPEG))
1844 : {
1845 97 : bTileDriverSupports1Band = true;
1846 97 : if (bPartialTile || !bAllDirty || (nBands == 2 && !bAllOpaque) ||
1847 19 : (nBands == 4 && !bAllOpaque) || m_poCT != nullptr)
1848 : {
1849 89 : pszDriverName = "PNG";
1850 89 : bTileDriverSupports2Bands = m_bPNGSupports2Bands;
1851 89 : bTileDriverSupports4Bands = true;
1852 89 : bTileDriverSupportsCT = m_bPNGSupportsCT;
1853 : }
1854 : else
1855 8 : pszDriverName = "JPEG";
1856 : }
1857 4070 : else if (m_eTF == GPKG_TF_PNG || m_eTF == GPKG_TF_PNG8)
1858 : {
1859 3847 : pszDriverName = "PNG";
1860 3847 : bTileDriverSupports1Band = true;
1861 3847 : bTileDriverSupports2Bands = m_bPNGSupports2Bands;
1862 3847 : bTileDriverSupports4Bands = true;
1863 3847 : bTileDriverSupportsCT = m_bPNGSupportsCT;
1864 : }
1865 223 : else if (m_eTF == GPKG_TF_JPEG)
1866 : {
1867 82 : pszDriverName = "JPEG";
1868 82 : bTileDriverSupports1Band = true;
1869 : }
1870 141 : else if (m_eTF == GPKG_TF_WEBP)
1871 : {
1872 83 : pszDriverName = "WEBP";
1873 83 : bTileDriverSupports4Bands = WEBPSupports4Bands();
1874 : }
1875 58 : else if (m_eTF == GPKG_TF_PNG_16BIT)
1876 : {
1877 53 : pszDriverName = "PNG";
1878 53 : eTileDT = GDT_UInt16;
1879 53 : bTileDriverSupports1Band = true;
1880 : }
1881 5 : else if (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
1882 : {
1883 5 : pszDriverName = "GTiff";
1884 5 : eTileDT = GDT_Float32;
1885 5 : bTileDriverSupports1Band = true;
1886 : }
1887 : else
1888 : {
1889 0 : CPLAssert(false);
1890 : }
1891 :
1892 : GDALDriver *l_poDriver =
1893 4167 : GDALDriver::FromHandle(GDALGetDriverByName(pszDriverName));
1894 4167 : if (l_poDriver != nullptr)
1895 : {
1896 4166 : auto poMEMDS = MEMDataset::Create("", nBlockXSize, nBlockYSize, 0,
1897 : eTileDT, nullptr);
1898 4166 : int nTileBands = nBands;
1899 4166 : if (bPartialTile && nBands == 1 && m_poCT == nullptr &&
1900 : bTileDriverSupports2Bands)
1901 33 : nTileBands = 2;
1902 4133 : else if (bPartialTile && bTileDriverSupports4Bands)
1903 40 : nTileBands = 4;
1904 : // only use (somewhat lossy) PNG8 if all bands are dirty
1905 4093 : else if (bAllDirty && m_eTF == GPKG_TF_PNG8 && nBands >= 3 &&
1906 10 : bAllOpaque && !bPartialTile)
1907 10 : nTileBands = 1;
1908 4083 : else if (nBands == 2)
1909 : {
1910 383 : if (bAllOpaque)
1911 : {
1912 51 : if (bTileDriverSupports2Bands)
1913 30 : nTileBands = 1;
1914 : else
1915 21 : nTileBands = 3;
1916 : }
1917 332 : else if (!bTileDriverSupports2Bands)
1918 : {
1919 135 : if (bTileDriverSupports4Bands)
1920 68 : nTileBands = 4;
1921 : else
1922 67 : nTileBands = 3;
1923 : }
1924 : }
1925 3700 : else if (nBands == 4 && (bAllOpaque || !bTileDriverSupports4Bands))
1926 21 : nTileBands = 3;
1927 3679 : else if (nBands == 1 && m_poCT != nullptr && !bTileDriverSupportsCT)
1928 : {
1929 1 : nTileBands = 3;
1930 1 : if (bTileDriverSupports4Bands)
1931 : {
1932 0 : for (int i = 0; i < m_poCT->GetColorEntryCount(); i++)
1933 : {
1934 0 : const GDALColorEntry *psEntry = m_poCT->GetColorEntry(i);
1935 0 : if (psEntry->c4 == 0)
1936 : {
1937 0 : nTileBands = 4;
1938 0 : break;
1939 : }
1940 : }
1941 1 : }
1942 : }
1943 3678 : else if (nBands == 1 && m_poCT == nullptr && !bTileDriverSupports1Band)
1944 1 : nTileBands = 3;
1945 :
1946 4166 : if (bPartialTile && (nTileBands == 2 || nTileBands == 4))
1947 : {
1948 73 : int nTargetAlphaBand = nTileBands;
1949 73 : memset(m_pabyCachedTiles + (nTargetAlphaBand - 1) * nBandBlockSize,
1950 : 0, nBandBlockSize);
1951 4398 : for (GPtrDiff_t iY = iYOff; iY < iYOff + iYCount; iY++)
1952 : {
1953 4325 : memset(m_pabyCachedTiles +
1954 4325 : (static_cast<size_t>(nTargetAlphaBand - 1) *
1955 4325 : nBlockYSize +
1956 4325 : iY) *
1957 4325 : nBlockXSize +
1958 4325 : iXOff,
1959 : 255, iXCount);
1960 : }
1961 : }
1962 :
1963 4166 : GUInt16 *pTempTileBuffer = nullptr;
1964 4166 : GPtrDiff_t nValidPixels = 0;
1965 4166 : double dfTileMin = 0.0;
1966 4166 : double dfTileMax = 0.0;
1967 4166 : double dfTileMean = 0.0;
1968 4166 : double dfTileStdDev = 0.0;
1969 4166 : double dfTileOffset = 0.0;
1970 4166 : double dfTileScale = 1.0;
1971 4166 : if (m_eTF == GPKG_TF_PNG_16BIT)
1972 : {
1973 : pTempTileBuffer = static_cast<GUInt16 *>(
1974 53 : VSI_MALLOC3_VERBOSE(2, nBlockXSize, nBlockYSize));
1975 :
1976 53 : if (m_eDT == GDT_Int16)
1977 : {
1978 74 : ProcessInt16UInt16Tile<GInt16>(
1979 37 : m_pabyCachedTiles,
1980 37 : static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize, true,
1981 37 : CPL_TO_BOOL(bHasNoData), dfNoDataValue, m_usGPKGNull,
1982 : m_dfOffset, m_dfScale, pTempTileBuffer, dfTileOffset,
1983 : dfTileScale, dfTileMin, dfTileMax, dfTileMean, dfTileStdDev,
1984 : nValidPixels);
1985 : }
1986 16 : else if (m_eDT == GDT_UInt16)
1987 : {
1988 22 : ProcessInt16UInt16Tile<GUInt16>(
1989 11 : m_pabyCachedTiles,
1990 11 : static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize, false,
1991 11 : CPL_TO_BOOL(bHasNoData), dfNoDataValue, m_usGPKGNull,
1992 : m_dfOffset, m_dfScale, pTempTileBuffer, dfTileOffset,
1993 : dfTileScale, dfTileMin, dfTileMax, dfTileMean, dfTileStdDev,
1994 : nValidPixels);
1995 : }
1996 5 : else if (m_eDT == GDT_Float32)
1997 : {
1998 5 : const float *pSrc =
1999 : reinterpret_cast<float *>(m_pabyCachedTiles);
2000 5 : float fMin = 0.0f;
2001 5 : float fMax = 0.0f;
2002 5 : double dfM2 = 0.0;
2003 327685 : for (GPtrDiff_t i = 0;
2004 327685 : i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize;
2005 : i++)
2006 : {
2007 327680 : const float fVal = pSrc[i];
2008 327680 : if (bHasNanNoData)
2009 : {
2010 0 : if (std::isnan(fVal))
2011 196206 : continue;
2012 : }
2013 720494 : else if (bHasNoData &&
2014 327680 : fVal == static_cast<float>(dfNoDataValue))
2015 : {
2016 196206 : continue;
2017 : }
2018 131474 : if (std::isinf(fVal))
2019 0 : continue;
2020 :
2021 131474 : if (nValidPixels == 0)
2022 : {
2023 5 : fMin = fVal;
2024 5 : fMax = fVal;
2025 : }
2026 : else
2027 : {
2028 131469 : fMin = std::min(fMin, fVal);
2029 131469 : fMax = std::max(fMax, fVal);
2030 : }
2031 131474 : nValidPixels++;
2032 131474 : const double dfDelta = fVal - dfTileMean;
2033 131474 : dfTileMean += dfDelta / nValidPixels;
2034 131474 : dfM2 += dfDelta * (fVal - dfTileMean);
2035 : }
2036 5 : dfTileMin = fMin;
2037 5 : dfTileMax = fMax;
2038 5 : if (nValidPixels)
2039 5 : dfTileStdDev = sqrt(dfM2 / nValidPixels);
2040 :
2041 5 : double dfGlobalMin = (fMin - m_dfOffset) / m_dfScale;
2042 5 : double dfGlobalMax = (fMax - m_dfOffset) / m_dfScale;
2043 5 : if (dfGlobalMax > dfGlobalMin)
2044 : {
2045 4 : if (bHasNoData && m_usGPKGNull == 65535 &&
2046 2 : dfGlobalMax - dfGlobalMin >= 65534.0)
2047 : {
2048 1 : dfTileOffset = dfGlobalMin;
2049 1 : dfTileScale = (dfGlobalMax - dfGlobalMin) / 65534.0;
2050 : }
2051 3 : else if (bHasNoData && m_usGPKGNull == 0 &&
2052 0 : (dfNoDataValue - m_dfOffset) / m_dfScale != 0)
2053 : {
2054 0 : dfTileOffset =
2055 0 : (65535.0 * dfGlobalMin - dfGlobalMax) / 65534.0;
2056 0 : dfTileScale = dfGlobalMin - dfTileOffset;
2057 : }
2058 : else
2059 : {
2060 3 : dfTileOffset = dfGlobalMin;
2061 3 : dfTileScale = (dfGlobalMax - dfGlobalMin) / 65535.0;
2062 : }
2063 : }
2064 : else
2065 : {
2066 1 : dfTileOffset = dfGlobalMin;
2067 : }
2068 :
2069 327685 : for (GPtrDiff_t i = 0;
2070 327685 : i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize;
2071 : i++)
2072 : {
2073 327680 : const float fVal = pSrc[i];
2074 327680 : if (bHasNanNoData)
2075 : {
2076 0 : if (std::isnan(fVal))
2077 : {
2078 0 : pTempTileBuffer[i] = m_usGPKGNull;
2079 0 : continue;
2080 : }
2081 : }
2082 327680 : else if (bHasNoData)
2083 : {
2084 196608 : if (fVal == static_cast<float>(dfNoDataValue))
2085 : {
2086 196206 : pTempTileBuffer[i] = m_usGPKGNull;
2087 196206 : continue;
2088 : }
2089 : }
2090 : double dfVal =
2091 131474 : std::isfinite(fVal)
2092 131474 : ? ((fVal - m_dfOffset) / m_dfScale - dfTileOffset) /
2093 : dfTileScale
2094 : : (fVal > 0) ? 65535
2095 131474 : : 0;
2096 131474 : CPLAssert(dfVal >= 0.0 && dfVal < 65535.5);
2097 131474 : pTempTileBuffer[i] = static_cast<GUInt16>(dfVal + 0.5);
2098 131474 : if (bHasNoData && pTempTileBuffer[i] == m_usGPKGNull)
2099 : {
2100 1 : if (m_usGPKGNull > 0)
2101 1 : pTempTileBuffer[i]--;
2102 : else
2103 0 : pTempTileBuffer[i]++;
2104 : }
2105 : }
2106 : }
2107 :
2108 53 : auto hBand = MEMCreateRasterBandEx(
2109 : poMEMDS, 1, reinterpret_cast<GByte *>(pTempTileBuffer),
2110 : GDT_UInt16, 0, 0, false);
2111 53 : poMEMDS->AddMEMBand(hBand);
2112 : }
2113 4113 : else if (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
2114 : {
2115 5 : const float *pSrc = reinterpret_cast<float *>(m_pabyCachedTiles);
2116 5 : float fMin = 0.0f;
2117 5 : float fMax = 0.0f;
2118 5 : double dfM2 = 0.0;
2119 327685 : for (GPtrDiff_t i = 0;
2120 327685 : i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize; i++)
2121 : {
2122 327680 : const float fVal = pSrc[i];
2123 327680 : if (bHasNanNoData)
2124 : {
2125 0 : if (std::isnan(fVal))
2126 65136 : continue;
2127 : }
2128 458352 : else if (bHasNoData &&
2129 327680 : fVal == static_cast<float>(dfNoDataValue))
2130 : {
2131 65136 : continue;
2132 : }
2133 :
2134 262544 : if (nValidPixels == 0)
2135 : {
2136 5 : fMin = fVal;
2137 5 : fMax = fVal;
2138 : }
2139 : else
2140 : {
2141 262539 : fMin = std::min(fMin, fVal);
2142 262539 : fMax = std::max(fMax, fVal);
2143 : }
2144 262544 : nValidPixels++;
2145 262544 : const double dfDelta = fVal - dfTileMean;
2146 262544 : dfTileMean += dfDelta / nValidPixels;
2147 262544 : dfM2 += dfDelta * (fVal - dfTileMean);
2148 : }
2149 5 : dfTileMin = fMin;
2150 5 : dfTileMax = fMax;
2151 5 : if (nValidPixels)
2152 5 : dfTileStdDev = sqrt(dfM2 / nValidPixels);
2153 :
2154 5 : auto hBand = MEMCreateRasterBandEx(poMEMDS, 1, m_pabyCachedTiles,
2155 : GDT_Float32, 0, 0, false);
2156 5 : poMEMDS->AddMEMBand(hBand);
2157 : }
2158 : else
2159 : {
2160 4108 : CPLAssert(m_eDT == GDT_UInt8);
2161 16196 : for (int i = 0; i < nTileBands; i++)
2162 : {
2163 12088 : int iSrc = i;
2164 12088 : if (nBands == 1 && m_poCT == nullptr && nTileBands == 3)
2165 3 : iSrc = 0;
2166 12085 : else if (nBands == 1 && m_poCT == nullptr && bPartialTile &&
2167 : nTileBands == 4)
2168 8 : iSrc = (i < 3) ? 0 : 3;
2169 12077 : else if (nBands == 2 && nTileBands >= 3)
2170 536 : iSrc = (i < 3) ? 0 : 1;
2171 :
2172 24176 : auto hBand = MEMCreateRasterBandEx(
2173 : poMEMDS, i + 1,
2174 12088 : m_pabyCachedTiles + iSrc * nBlockXSize * nBlockYSize,
2175 : GDT_UInt8, 0, 0, false);
2176 12088 : poMEMDS->AddMEMBand(hBand);
2177 :
2178 12088 : if (i == 0 && nTileBands == 1 && m_poCT != nullptr)
2179 15 : poMEMDS->GetRasterBand(1)->SetColorTable(m_poCT);
2180 : }
2181 : }
2182 :
2183 4166 : if ((m_eTF == GPKG_TF_PNG_16BIT || m_eTF == GPKG_TF_TIFF_32BIT_FLOAT) &&
2184 58 : nValidPixels == 0)
2185 : {
2186 : // If tile is fully transparent, don't serialize it and remove
2187 : // it if it exists.
2188 1 : GIntBig nId = GetTileId(nRow, nCol);
2189 1 : if (nId > 0)
2190 : {
2191 1 : DeleteTile(nRow, nCol);
2192 :
2193 1 : DeleteFromGriddedTileAncillary(nId);
2194 : }
2195 :
2196 1 : CPLFree(pTempTileBuffer);
2197 1 : delete poMEMDS;
2198 2 : return CE_None;
2199 : }
2200 :
2201 4165 : if (m_eTF == GPKG_TF_PNG8 && nTileBands == 1 && nBands >= 3)
2202 : {
2203 10 : CPLAssert(bAllDirty);
2204 :
2205 10 : auto poMEM_RGB_DS = MEMDataset::Create("", nBlockXSize, nBlockYSize,
2206 : 0, GDT_UInt8, nullptr);
2207 40 : for (int i = 0; i < 3; i++)
2208 : {
2209 60 : auto hBand = MEMCreateRasterBandEx(
2210 30 : poMEMDS, i + 1, m_pabyCachedTiles + i * nBandBlockSize,
2211 : GDT_UInt8, 0, 0, false);
2212 30 : poMEM_RGB_DS->AddMEMBand(hBand);
2213 : }
2214 :
2215 10 : if (m_pabyHugeColorArray == nullptr)
2216 : {
2217 5 : if (nBlockXSize <= 65536 / nBlockYSize)
2218 4 : m_pabyHugeColorArray =
2219 4 : VSIMalloc(MEDIAN_CUT_AND_DITHER_BUFFER_SIZE_65536);
2220 : else
2221 1 : m_pabyHugeColorArray =
2222 1 : VSIMalloc2(256 * 256 * 256, sizeof(GUInt32));
2223 : }
2224 :
2225 10 : GDALColorTable *poCT = new GDALColorTable();
2226 20 : GDALComputeMedianCutPCTInternal(
2227 10 : poMEM_RGB_DS->GetRasterBand(1), poMEM_RGB_DS->GetRasterBand(2),
2228 10 : poMEM_RGB_DS->GetRasterBand(3),
2229 : /*NULL, NULL, NULL,*/
2230 10 : m_pabyCachedTiles, m_pabyCachedTiles + nBandBlockSize,
2231 10 : m_pabyCachedTiles + 2 * nBandBlockSize, nullptr,
2232 : 256, /* max colors */
2233 : 8, /* bit depth */
2234 : static_cast<GUInt32 *>(
2235 10 : m_pabyHugeColorArray), /* preallocated histogram */
2236 : poCT, nullptr, nullptr);
2237 :
2238 20 : GDALDitherRGB2PCTInternal(
2239 10 : poMEM_RGB_DS->GetRasterBand(1), poMEM_RGB_DS->GetRasterBand(2),
2240 10 : poMEM_RGB_DS->GetRasterBand(3), poMEMDS->GetRasterBand(1), poCT,
2241 : 8, /* bit depth */
2242 : static_cast<GInt16 *>(
2243 10 : m_pabyHugeColorArray), /* pasDynamicColorMap */
2244 10 : m_bDither, nullptr, nullptr);
2245 10 : poMEMDS->GetRasterBand(1)->SetColorTable(poCT);
2246 10 : delete poCT;
2247 10 : GDALClose(poMEM_RGB_DS);
2248 : }
2249 4155 : else if (nBands == 1 && m_poCT != nullptr && nTileBands > 1)
2250 : {
2251 : GByte abyCT[4 * 256];
2252 7 : const int nEntries = std::min(256, m_poCT->GetColorEntryCount());
2253 1799 : for (int i = 0; i < nEntries; i++)
2254 : {
2255 1792 : const GDALColorEntry *psEntry = m_poCT->GetColorEntry(i);
2256 1792 : abyCT[4 * i] = static_cast<GByte>(psEntry->c1);
2257 1792 : abyCT[4 * i + 1] = static_cast<GByte>(psEntry->c2);
2258 1792 : abyCT[4 * i + 2] = static_cast<GByte>(psEntry->c3);
2259 1792 : abyCT[4 * i + 3] = static_cast<GByte>(psEntry->c4);
2260 : }
2261 7 : for (int i = nEntries; i < 256; i++)
2262 : {
2263 0 : abyCT[4 * i] = 0;
2264 0 : abyCT[4 * i + 1] = 0;
2265 0 : abyCT[4 * i + 2] = 0;
2266 0 : abyCT[4 * i + 3] = 0;
2267 : }
2268 7 : if (iYOff > 0)
2269 : {
2270 0 : memset(m_pabyCachedTiles + 0 * nBandBlockSize, 0,
2271 0 : nBlockXSize * iYOff);
2272 0 : memset(m_pabyCachedTiles + 1 * nBandBlockSize, 0,
2273 0 : nBlockXSize * iYOff);
2274 0 : memset(m_pabyCachedTiles + 2 * nBandBlockSize, 0,
2275 0 : nBlockXSize * iYOff);
2276 0 : memset(m_pabyCachedTiles + 3 * nBandBlockSize, 0,
2277 0 : nBlockXSize * iYOff);
2278 : }
2279 1057 : for (GPtrDiff_t iY = iYOff; iY < iYOff + iYCount; iY++)
2280 : {
2281 1050 : if (iXOff > 0)
2282 : {
2283 0 : const GPtrDiff_t i = iY * nBlockXSize;
2284 0 : memset(m_pabyCachedTiles + 0 * nBandBlockSize + i, 0,
2285 : iXOff);
2286 0 : memset(m_pabyCachedTiles + 1 * nBandBlockSize + i, 0,
2287 : iXOff);
2288 0 : memset(m_pabyCachedTiles + 2 * nBandBlockSize + i, 0,
2289 : iXOff);
2290 0 : memset(m_pabyCachedTiles + 3 * nBandBlockSize + i, 0,
2291 : iXOff);
2292 : }
2293 183250 : for (int iX = iXOff; iX < iXOff + iXCount; iX++)
2294 : {
2295 182200 : const GPtrDiff_t i = iY * nBlockXSize + iX;
2296 182200 : GByte byVal = m_pabyCachedTiles[i];
2297 182200 : m_pabyCachedTiles[i] = abyCT[4 * byVal];
2298 182200 : m_pabyCachedTiles[i + 1 * nBandBlockSize] =
2299 182200 : abyCT[4 * byVal + 1];
2300 182200 : m_pabyCachedTiles[i + 2 * nBandBlockSize] =
2301 182200 : abyCT[4 * byVal + 2];
2302 182200 : m_pabyCachedTiles[i + 3 * nBandBlockSize] =
2303 182200 : abyCT[4 * byVal + 3];
2304 : }
2305 1050 : if (iXOff + iXCount < nBlockXSize)
2306 : {
2307 800 : const GPtrDiff_t i = iY * nBlockXSize + iXOff + iXCount;
2308 800 : memset(m_pabyCachedTiles + 0 * nBandBlockSize + i, 0,
2309 800 : nBlockXSize - (iXOff + iXCount));
2310 800 : memset(m_pabyCachedTiles + 1 * nBandBlockSize + i, 0,
2311 800 : nBlockXSize - (iXOff + iXCount));
2312 800 : memset(m_pabyCachedTiles + 2 * nBandBlockSize + i, 0,
2313 800 : nBlockXSize - (iXOff + iXCount));
2314 800 : memset(m_pabyCachedTiles + 3 * nBandBlockSize + i, 0,
2315 800 : nBlockXSize - (iXOff + iXCount));
2316 : }
2317 : }
2318 7 : if (iYOff + iYCount < nBlockYSize)
2319 : {
2320 7 : const GPtrDiff_t i = (iYOff + iYCount) * nBlockXSize;
2321 7 : memset(m_pabyCachedTiles + 0 * nBandBlockSize + i, 0,
2322 7 : nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
2323 7 : memset(m_pabyCachedTiles + 1 * nBandBlockSize + i, 0,
2324 7 : nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
2325 7 : memset(m_pabyCachedTiles + 2 * nBandBlockSize + i, 0,
2326 7 : nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
2327 7 : memset(m_pabyCachedTiles + 3 * nBandBlockSize + i, 0,
2328 7 : nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
2329 : }
2330 : }
2331 :
2332 : char **papszDriverOptions =
2333 4165 : CSLSetNameValue(nullptr, "_INTERNAL_DATASET", "YES");
2334 4165 : if (EQUAL(pszDriverName, "JPEG") || EQUAL(pszDriverName, "WEBP"))
2335 : {
2336 : // If not all bands are dirty, then use lossless WEBP
2337 172 : if (!bAllDirty && EQUAL(pszDriverName, "WEBP"))
2338 : {
2339 2 : papszDriverOptions =
2340 2 : CSLSetNameValue(papszDriverOptions, "LOSSLESS", "YES");
2341 : }
2342 : else
2343 : {
2344 : papszDriverOptions =
2345 170 : CSLSetNameValue(papszDriverOptions, "QUALITY",
2346 : CPLSPrintf("%d", m_nQuality));
2347 : }
2348 : }
2349 3993 : else if (EQUAL(pszDriverName, "PNG"))
2350 : {
2351 3988 : papszDriverOptions = CSLSetNameValue(papszDriverOptions, "ZLEVEL",
2352 : CPLSPrintf("%d", m_nZLevel));
2353 : }
2354 5 : else if (EQUAL(pszDriverName, "GTiff"))
2355 : {
2356 : papszDriverOptions =
2357 5 : CSLSetNameValue(papszDriverOptions, "COMPRESS", "LZW");
2358 5 : if (nBlockXSize * nBlockYSize <= 512 * 512)
2359 : {
2360 : // If tile is not too big, create it as single-strip TIFF
2361 : papszDriverOptions =
2362 5 : CSLSetNameValue(papszDriverOptions, "BLOCKYSIZE",
2363 : CPLSPrintf("%d", nBlockYSize));
2364 : }
2365 : }
2366 : #ifdef DEBUG
2367 : VSIStatBufL sStat;
2368 4165 : CPLAssert(VSIStatL(osMemFileName, &sStat) != 0);
2369 : #endif
2370 : GDALDataset *poOutDS =
2371 4165 : l_poDriver->CreateCopy(osMemFileName, poMEMDS, FALSE,
2372 : papszDriverOptions, nullptr, nullptr);
2373 4165 : CSLDestroy(papszDriverOptions);
2374 4165 : CPLFree(pTempTileBuffer);
2375 :
2376 4165 : if (poOutDS)
2377 : {
2378 4165 : GDALClose(poOutDS);
2379 4165 : vsi_l_offset nBlobSize = 0;
2380 : GByte *pabyBlob =
2381 4165 : VSIGetMemFileBuffer(osMemFileName, &nBlobSize, TRUE);
2382 :
2383 : /* Create or commit and recreate transaction */
2384 4165 : GDALGPKGMBTilesLikePseudoDataset *poMainDS =
2385 4165 : m_poParentDS ? m_poParentDS : this;
2386 4165 : if (poMainDS->m_nTileInsertionCount == 0)
2387 : {
2388 191 : poMainDS->IStartTransaction();
2389 : }
2390 3974 : else if (poMainDS->m_nTileInsertionCount == 1000)
2391 : {
2392 3 : if (poMainDS->ICommitTransaction() != OGRERR_NONE)
2393 : {
2394 1 : poMainDS->m_nTileInsertionCount = -1;
2395 1 : CPLFree(pabyBlob);
2396 1 : VSIUnlink(osMemFileName);
2397 1 : delete poMEMDS;
2398 1 : return CE_Failure;
2399 : }
2400 2 : poMainDS->IStartTransaction();
2401 2 : poMainDS->m_nTileInsertionCount = 0;
2402 : }
2403 4164 : poMainDS->m_nTileInsertionCount++;
2404 :
2405 : char *pszSQL =
2406 4164 : sqlite3_mprintf("INSERT OR REPLACE INTO \"%w\" "
2407 : "(zoom_level, tile_row, tile_column, "
2408 : "tile_data) VALUES (%d, %d, %d, ?)",
2409 : m_osRasterTable.c_str(), m_nZoomLevel,
2410 4164 : GetRowFromIntoTopConvention(nRow), nCol);
2411 : #ifdef DEBUG_VERBOSE
2412 : CPLDebug("GPKG", "%s", pszSQL);
2413 : #endif
2414 4164 : sqlite3_stmt *hStmt = nullptr;
2415 4164 : int rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt, nullptr);
2416 4164 : if (rc != SQLITE_OK)
2417 : {
2418 0 : CPLFree(pabyBlob);
2419 : }
2420 : else
2421 : {
2422 4164 : sqlite3_bind_blob(hStmt, 1, pabyBlob,
2423 : static_cast<int>(nBlobSize), CPLFree);
2424 4164 : rc = sqlite3_step(hStmt);
2425 4164 : if (rc == SQLITE_DONE)
2426 4164 : eErr = CE_None;
2427 : else
2428 : {
2429 0 : CPLError(CE_Failure, CPLE_AppDefined,
2430 : "Failure when inserting tile (row=%d,col=%d) at "
2431 : "zoom_level=%d : %s",
2432 0 : GetRowFromIntoTopConvention(nRow), nCol,
2433 0 : m_nZoomLevel, sqlite3_errmsg(IGetDB()));
2434 : }
2435 : }
2436 4164 : sqlite3_finalize(hStmt);
2437 4164 : sqlite3_free(pszSQL);
2438 :
2439 4164 : if (m_eTF == GPKG_TF_PNG_16BIT || m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
2440 : {
2441 57 : GIntBig nTileId = GetTileId(nRow, nCol);
2442 57 : if (nTileId == 0)
2443 0 : eErr = CE_Failure;
2444 : else
2445 : {
2446 57 : DeleteFromGriddedTileAncillary(nTileId);
2447 :
2448 57 : pszSQL = sqlite3_mprintf(
2449 : "INSERT INTO gpkg_2d_gridded_tile_ancillary "
2450 : "(tpudt_name, tpudt_id, scale, offset, min, max, "
2451 : "mean, std_dev) VALUES "
2452 : "('%q', ?, %.17g, %.17g, ?, ?, ?, ?)",
2453 : m_osRasterTable.c_str(), dfTileScale, dfTileOffset);
2454 : #ifdef DEBUG_VERBOSE
2455 : CPLDebug("GPKG", "%s", pszSQL);
2456 : #endif
2457 57 : hStmt = nullptr;
2458 57 : rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt,
2459 : nullptr);
2460 57 : if (rc != SQLITE_OK)
2461 : {
2462 0 : eErr = CE_Failure;
2463 : }
2464 : else
2465 : {
2466 57 : sqlite3_bind_int64(hStmt, 1, nTileId);
2467 57 : sqlite3_bind_double(hStmt, 2, dfTileMin);
2468 57 : sqlite3_bind_double(hStmt, 3, dfTileMax);
2469 57 : sqlite3_bind_double(hStmt, 4, dfTileMean);
2470 57 : sqlite3_bind_double(hStmt, 5, dfTileStdDev);
2471 57 : rc = sqlite3_step(hStmt);
2472 57 : if (rc == SQLITE_DONE)
2473 : {
2474 57 : eErr = CE_None;
2475 : }
2476 : else
2477 : {
2478 0 : CPLError(CE_Failure, CPLE_AppDefined,
2479 : "Cannot insert into "
2480 : "gpkg_2d_gridded_tile_ancillary");
2481 0 : eErr = CE_Failure;
2482 : }
2483 : }
2484 57 : sqlite3_finalize(hStmt);
2485 57 : sqlite3_free(pszSQL);
2486 : }
2487 : }
2488 : }
2489 :
2490 4164 : VSIUnlink(osMemFileName);
2491 4164 : delete poMEMDS;
2492 : }
2493 : else
2494 : {
2495 1 : CPLError(CE_Failure, CPLE_NotSupported, "Cannot find driver %s",
2496 : pszDriverName);
2497 : }
2498 :
2499 4165 : return eErr;
2500 : }
2501 :
2502 : /************************************************************************/
2503 : /* FlushRemainingShiftedTiles() */
2504 : /************************************************************************/
2505 :
2506 : CPLErr
2507 609 : GDALGPKGMBTilesLikePseudoDataset::FlushRemainingShiftedTiles(bool bPartialFlush)
2508 : {
2509 609 : if (m_hTempDB == nullptr)
2510 460 : return CE_None;
2511 :
2512 745 : for (int i = 0; i <= 3; i++)
2513 : {
2514 596 : m_asCachedTilesDesc[i].nRow = -1;
2515 596 : m_asCachedTilesDesc[i].nCol = -1;
2516 596 : m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
2517 : }
2518 :
2519 : int nBlockXSize, nBlockYSize;
2520 149 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
2521 149 : const int nBands = IGetRasterCount();
2522 149 : const int nRasterXSize = IGetRasterBand(1)->GetXSize();
2523 149 : const int nRasterYSize = IGetRasterBand(1)->GetYSize();
2524 149 : const int nXBlocks = DIV_ROUND_UP(nRasterXSize, nBlockXSize);
2525 149 : const int nYBlocks = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
2526 :
2527 149 : int nPartialActiveTiles = 0;
2528 149 : if (bPartialFlush)
2529 : {
2530 2 : sqlite3_stmt *hStmt = nullptr;
2531 4 : CPLString osSQL;
2532 : osSQL.Printf("SELECT COUNT(*) FROM partial_tiles WHERE zoom_level = %d "
2533 : "AND partial_flag != 0",
2534 2 : m_nZoomLevel);
2535 2 : if (SQLPrepareWithError(m_hTempDB, osSQL.c_str(), -1, &hStmt,
2536 2 : nullptr) == SQLITE_OK)
2537 : {
2538 2 : if (sqlite3_step(hStmt) == SQLITE_ROW)
2539 : {
2540 2 : nPartialActiveTiles = sqlite3_column_int(hStmt, 0);
2541 2 : CPLDebug("GPKG", "Active partial tiles before flush: %d",
2542 : nPartialActiveTiles);
2543 : }
2544 2 : sqlite3_finalize(hStmt);
2545 : }
2546 : }
2547 :
2548 298 : CPLString osSQL = "SELECT tile_row, tile_column, partial_flag";
2549 616 : for (int nBand = 1; nBand <= nBands; nBand++)
2550 : {
2551 467 : osSQL += CPLSPrintf(", tile_data_band_%d", nBand);
2552 : }
2553 : osSQL += CPLSPrintf(" FROM partial_tiles WHERE "
2554 : "zoom_level = %d AND partial_flag != 0",
2555 149 : m_nZoomLevel);
2556 149 : if (bPartialFlush)
2557 : {
2558 2 : osSQL += " ORDER BY age";
2559 : }
2560 149 : const char *pszSQL = osSQL.c_str();
2561 :
2562 : #ifdef DEBUG_VERBOSE
2563 : CPLDebug("GPKG", "%s", pszSQL);
2564 : #endif
2565 149 : sqlite3_stmt *hStmt = nullptr;
2566 149 : int rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
2567 149 : if (rc != SQLITE_OK)
2568 : {
2569 0 : return CE_Failure;
2570 : }
2571 :
2572 149 : CPLErr eErr = CE_None;
2573 149 : bool bGotPartialTiles = false;
2574 149 : int nCountFlushedTiles = 0;
2575 149 : const size_t nBandBlockSize =
2576 149 : static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
2577 175 : do
2578 : {
2579 324 : rc = sqlite3_step(hStmt);
2580 324 : if (rc == SQLITE_ROW)
2581 : {
2582 177 : bGotPartialTiles = true;
2583 :
2584 177 : int nRow = sqlite3_column_int(hStmt, 0);
2585 177 : int nCol = sqlite3_column_int(hStmt, 1);
2586 177 : int nPartialFlags = sqlite3_column_int(hStmt, 2);
2587 :
2588 177 : if (bPartialFlush)
2589 : {
2590 : // This method assumes that there are no dirty blocks alive
2591 : // so check this assumption.
2592 : // When called with bPartialFlush = false, FlushCache() has
2593 : // already been called, so no need to check.
2594 16 : bool bFoundDirtyBlock = false;
2595 16 : int nBlockXOff = nCol - m_nShiftXTiles;
2596 16 : int nBlockYOff = nRow - m_nShiftYTiles;
2597 96 : for (int iX = 0; !bFoundDirtyBlock &&
2598 48 : iX < ((m_nShiftXPixelsMod != 0) ? 2 : 1);
2599 : iX++)
2600 : {
2601 32 : if (nBlockXOff + iX < 0 || nBlockXOff + iX >= nXBlocks)
2602 12 : continue;
2603 120 : for (int iY = 0; !bFoundDirtyBlock &&
2604 60 : iY < ((m_nShiftYPixelsMod != 0) ? 2 : 1);
2605 : iY++)
2606 : {
2607 40 : if (nBlockYOff + iY < 0 || nBlockYOff + iY >= nYBlocks)
2608 17 : continue;
2609 69 : for (int iBand = 1;
2610 69 : !bFoundDirtyBlock && iBand <= nBands; iBand++)
2611 : {
2612 : GDALRasterBlock *poBlock =
2613 : cpl::down_cast<GDALGPKGMBTilesLikeRasterBand *>(
2614 46 : IGetRasterBand(iBand))
2615 46 : ->AccessibleTryGetLockedBlockRef(
2616 : nBlockXOff + iX, nBlockYOff + iY);
2617 46 : if (poBlock)
2618 : {
2619 0 : if (poBlock->GetDirty())
2620 0 : bFoundDirtyBlock = true;
2621 0 : poBlock->DropLock();
2622 : }
2623 : }
2624 : }
2625 : }
2626 16 : if (bFoundDirtyBlock)
2627 : {
2628 : #ifdef DEBUG_VERBOSE
2629 : CPLDebug("GPKG",
2630 : "Skipped flushing tile row = %d, column = %d "
2631 : "because it has dirty block(s) in GDAL cache",
2632 : nRow, nCol);
2633 : #endif
2634 0 : continue;
2635 : }
2636 : }
2637 :
2638 177 : nCountFlushedTiles++;
2639 177 : if (bPartialFlush && nCountFlushedTiles >= nPartialActiveTiles / 2)
2640 : {
2641 2 : CPLDebug("GPKG", "Flushed %d tiles", nCountFlushedTiles);
2642 2 : break;
2643 : }
2644 :
2645 705 : for (int nBand = 1; nBand <= nBands; nBand++)
2646 : {
2647 530 : if (nPartialFlags & (((1 << 4) - 1) << (4 * (nBand - 1))))
2648 : {
2649 498 : CPLAssert(sqlite3_column_bytes(hStmt, 2 + nBand) ==
2650 : static_cast<int>(nBandBlockSize));
2651 498 : memcpy(m_pabyCachedTiles + (nBand - 1) * nBandBlockSize,
2652 : sqlite3_column_blob(hStmt, 2 + nBand),
2653 : nBandBlockSize);
2654 : }
2655 : else
2656 : {
2657 32 : FillEmptyTileSingleBand(m_pabyCachedTiles +
2658 32 : (nBand - 1) * nBandBlockSize);
2659 : }
2660 : }
2661 :
2662 175 : int nFullFlags = (1 << (4 * nBands)) - 1;
2663 :
2664 : // In case the partial flags indicate that there's some quadrant
2665 : // missing, check in the main database if there is already a tile
2666 : // In which case, use the parts of that tile that aren't in the
2667 : // temporary database
2668 175 : if (nPartialFlags != nFullFlags)
2669 : {
2670 525 : char *pszNewSQL = sqlite3_mprintf(
2671 : "SELECT tile_data%s FROM \"%w\" "
2672 : "WHERE zoom_level = %d AND tile_row = %d AND tile_column = "
2673 : "%d%s",
2674 175 : m_eDT != GDT_UInt8 ? ", id"
2675 : : "", // MBTiles do not have an id
2676 : m_osRasterTable.c_str(), m_nZoomLevel,
2677 175 : GetRowFromIntoTopConvention(nRow), nCol,
2678 175 : !m_osWHERE.empty()
2679 0 : ? CPLSPrintf(" AND (%s)", m_osWHERE.c_str())
2680 : : "");
2681 : #ifdef DEBUG_VERBOSE
2682 : CPLDebug("GPKG", "%s", pszNewSQL);
2683 : #endif
2684 175 : sqlite3_stmt *hNewStmt = nullptr;
2685 175 : rc = SQLPrepareWithError(IGetDB(), pszNewSQL, -1, &hNewStmt,
2686 : nullptr);
2687 175 : if (rc == SQLITE_OK)
2688 : {
2689 175 : rc = sqlite3_step(hNewStmt);
2690 192 : if (rc == SQLITE_ROW &&
2691 17 : sqlite3_column_type(hNewStmt, 0) == SQLITE_BLOB)
2692 : {
2693 17 : const int nBytes = sqlite3_column_bytes(hNewStmt, 0);
2694 : GIntBig nTileId =
2695 17 : (m_eDT == GDT_UInt8)
2696 17 : ? 0
2697 0 : : sqlite3_column_int64(hNewStmt, 1);
2698 : GByte *pabyRawData =
2699 : const_cast<GByte *>(static_cast<const GByte *>(
2700 17 : sqlite3_column_blob(hNewStmt, 0)));
2701 : const CPLString osMemFileName(
2702 34 : VSIMemGenerateHiddenFilename("gpkg_read_tile"));
2703 17 : VSILFILE *fp = VSIFileFromMemBuffer(
2704 : osMemFileName.c_str(), pabyRawData, nBytes, FALSE);
2705 17 : VSIFCloseL(fp);
2706 :
2707 17 : double dfTileOffset = 0.0;
2708 17 : double dfTileScale = 1.0;
2709 17 : GetTileOffsetAndScale(nTileId, dfTileOffset,
2710 : dfTileScale);
2711 17 : const int nTileBands = m_eDT == GDT_UInt8 ? 4 : 1;
2712 17 : GByte *pabyTemp =
2713 17 : m_pabyCachedTiles + nTileBands * nBandBlockSize;
2714 17 : ReadTile(osMemFileName, pabyTemp, dfTileOffset,
2715 : dfTileScale);
2716 17 : VSIUnlink(osMemFileName);
2717 :
2718 17 : int iYQuadrantMax = (m_nShiftYPixelsMod) ? 1 : 0;
2719 17 : int iXQuadrantMax = (m_nShiftXPixelsMod) ? 1 : 0;
2720 51 : for (int iYQuadrant = 0; iYQuadrant <= iYQuadrantMax;
2721 : iYQuadrant++)
2722 : {
2723 100 : for (int iXQuadrant = 0;
2724 100 : iXQuadrant <= iXQuadrantMax; iXQuadrant++)
2725 : {
2726 226 : for (int nBand = 1; nBand <= nBands; nBand++)
2727 : {
2728 160 : int iQuadrantFlag = 0;
2729 160 : if (iXQuadrant == 0 && iYQuadrant == 0)
2730 42 : iQuadrantFlag |= (1 << 0);
2731 160 : if (iXQuadrant == iXQuadrantMax &&
2732 : iYQuadrant == 0)
2733 42 : iQuadrantFlag |= (1 << 1);
2734 160 : if (iXQuadrant == 0 &&
2735 : iYQuadrant == iYQuadrantMax)
2736 42 : iQuadrantFlag |= (1 << 2);
2737 160 : if (iXQuadrant == iXQuadrantMax &&
2738 : iYQuadrant == iYQuadrantMax)
2739 42 : iQuadrantFlag |= (1 << 3);
2740 160 : int nLocalFlag = iQuadrantFlag
2741 160 : << (4 * (nBand - 1));
2742 160 : if (!(nPartialFlags & nLocalFlag))
2743 : {
2744 : int nXOff, nYOff, nXSize, nYSize;
2745 118 : if (iXQuadrant == 0 &&
2746 60 : m_nShiftXPixelsMod != 0)
2747 : {
2748 56 : nXOff = 0;
2749 56 : nXSize = m_nShiftXPixelsMod;
2750 : }
2751 : else
2752 : {
2753 62 : nXOff = m_nShiftXPixelsMod;
2754 62 : nXSize = nBlockXSize -
2755 62 : m_nShiftXPixelsMod;
2756 : }
2757 118 : if (iYQuadrant == 0 &&
2758 63 : m_nShiftYPixelsMod != 0)
2759 : {
2760 63 : nYOff = 0;
2761 63 : nYSize = m_nShiftYPixelsMod;
2762 : }
2763 : else
2764 : {
2765 55 : nYOff = m_nShiftYPixelsMod;
2766 55 : nYSize = nBlockYSize -
2767 55 : m_nShiftYPixelsMod;
2768 : }
2769 118 : for (int iY = nYOff;
2770 13888 : iY < nYOff + nYSize; iY++)
2771 : {
2772 13770 : memcpy(m_pabyCachedTiles +
2773 13770 : ((static_cast<size_t>(
2774 13770 : nBand - 1) *
2775 13770 : nBlockYSize +
2776 13770 : iY) *
2777 13770 : nBlockXSize +
2778 13770 : nXOff) *
2779 13770 : m_nDTSize,
2780 13770 : pabyTemp +
2781 13770 : ((static_cast<size_t>(
2782 13770 : nBand - 1) *
2783 13770 : nBlockYSize +
2784 13770 : iY) *
2785 13770 : nBlockXSize +
2786 13770 : nXOff) *
2787 13770 : m_nDTSize,
2788 13770 : static_cast<size_t>(nXSize) *
2789 13770 : m_nDTSize);
2790 : }
2791 : }
2792 : }
2793 : }
2794 : }
2795 : }
2796 158 : else if (rc != SQLITE_DONE)
2797 : {
2798 0 : CPLError(CE_Failure, CPLE_AppDefined,
2799 : "sqlite3_step(%s) failed: %s", pszNewSQL,
2800 : sqlite3_errmsg(m_hTempDB));
2801 : }
2802 175 : sqlite3_finalize(hNewStmt);
2803 : }
2804 175 : sqlite3_free(pszNewSQL);
2805 : }
2806 :
2807 175 : m_asCachedTilesDesc[0].nRow = nRow;
2808 175 : m_asCachedTilesDesc[0].nCol = nCol;
2809 175 : m_asCachedTilesDesc[0].nIdxWithinTileData = 0;
2810 175 : m_asCachedTilesDesc[0].abBandDirty[0] = true;
2811 175 : m_asCachedTilesDesc[0].abBandDirty[1] = true;
2812 175 : m_asCachedTilesDesc[0].abBandDirty[2] = true;
2813 175 : m_asCachedTilesDesc[0].abBandDirty[3] = true;
2814 :
2815 175 : eErr = WriteTile();
2816 :
2817 175 : if (eErr == CE_None && bPartialFlush)
2818 : {
2819 : pszSQL =
2820 14 : CPLSPrintf("DELETE FROM partial_tiles WHERE zoom_level = "
2821 : "%d AND tile_row = %d AND tile_column = %d",
2822 : m_nZoomLevel, nRow, nCol);
2823 : #ifdef DEBUG_VERBOSE
2824 : CPLDebug("GPKG", "%s", pszSQL);
2825 : #endif
2826 14 : if (SQLCommand(m_hTempDB, pszSQL) != OGRERR_NONE)
2827 0 : eErr = CE_None;
2828 : }
2829 : }
2830 : else
2831 : {
2832 147 : if (rc != SQLITE_DONE)
2833 : {
2834 0 : CPLError(CE_Failure, CPLE_AppDefined,
2835 : "sqlite3_step(%s) failed: %s", pszSQL,
2836 : sqlite3_errmsg(m_hTempDB));
2837 : }
2838 147 : break;
2839 : }
2840 175 : } while (eErr == CE_None);
2841 :
2842 149 : sqlite3_finalize(hStmt);
2843 :
2844 149 : if (bPartialFlush && nCountFlushedTiles < nPartialActiveTiles / 2)
2845 : {
2846 0 : CPLDebug("GPKG", "Flushed %d tiles. Target was %d", nCountFlushedTiles,
2847 : nPartialActiveTiles / 2);
2848 : }
2849 :
2850 149 : if (bGotPartialTiles && !bPartialFlush)
2851 : {
2852 : #ifdef DEBUG_VERBOSE
2853 : pszSQL = CPLSPrintf("SELECT p1.id, p1.tile_row, p1.tile_column FROM "
2854 : "partial_tiles p1, partial_tiles p2 "
2855 : "WHERE p1.zoom_level = %d AND p2.zoom_level = %d "
2856 : "AND p1.tile_row = p2.tile_row AND p1.tile_column "
2857 : "= p2.tile_column AND p2.partial_flag != 0",
2858 : -1 - m_nZoomLevel, m_nZoomLevel);
2859 : rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
2860 : CPLAssert(rc == SQLITE_OK);
2861 : while ((rc = sqlite3_step(hStmt)) == SQLITE_ROW)
2862 : {
2863 : CPLDebug("GPKG",
2864 : "Conflict: existing id = %d, tile_row = %d, tile_column = "
2865 : "%d, zoom_level = %d",
2866 : sqlite3_column_int(hStmt, 0), sqlite3_column_int(hStmt, 1),
2867 : sqlite3_column_int(hStmt, 2), m_nZoomLevel);
2868 : }
2869 : sqlite3_finalize(hStmt);
2870 : #endif
2871 :
2872 110 : pszSQL = CPLSPrintf("UPDATE partial_tiles SET zoom_level = %d, "
2873 : "partial_flag = 0, age = -1 WHERE zoom_level = %d "
2874 : "AND partial_flag != 0",
2875 55 : -1 - m_nZoomLevel, m_nZoomLevel);
2876 : #ifdef DEBUG_VERBOSE
2877 : CPLDebug("GPKG", "%s", pszSQL);
2878 : #endif
2879 55 : SQLCommand(m_hTempDB, pszSQL);
2880 : }
2881 :
2882 149 : return eErr;
2883 : }
2884 :
2885 : /************************************************************************/
2886 : /* DoPartialFlushOfPartialTilesIfNecessary() */
2887 : /************************************************************************/
2888 :
2889 : CPLErr
2890 4556 : GDALGPKGMBTilesLikePseudoDataset::DoPartialFlushOfPartialTilesIfNecessary()
2891 : {
2892 4556 : time_t nCurTimeStamp = time(nullptr);
2893 4556 : if (m_nLastSpaceCheckTimestamp == 0)
2894 3 : m_nLastSpaceCheckTimestamp = nCurTimeStamp;
2895 4556 : if (m_nLastSpaceCheckTimestamp > 0 &&
2896 81 : (m_bForceTempDBCompaction ||
2897 13 : nCurTimeStamp - m_nLastSpaceCheckTimestamp > 10))
2898 : {
2899 68 : m_nLastSpaceCheckTimestamp = nCurTimeStamp;
2900 : GIntBig nFreeSpace =
2901 68 : VSIGetDiskFreeSpace(CPLGetDirnameSafe(m_osTempDBFilename).c_str());
2902 68 : bool bTryFreeing = false;
2903 68 : if (nFreeSpace >= 0 && nFreeSpace < 1024 * 1024 * 1024)
2904 : {
2905 0 : CPLDebug("GPKG",
2906 : "Free space below 1GB. Flushing part of partial tiles");
2907 0 : bTryFreeing = true;
2908 : }
2909 : else
2910 : {
2911 : VSIStatBufL sStat;
2912 68 : if (VSIStatL(m_osTempDBFilename, &sStat) == 0)
2913 : {
2914 68 : GIntBig nTempSpace = sStat.st_size;
2915 136 : if (VSIStatL((m_osTempDBFilename + "-journal").c_str(),
2916 68 : &sStat) == 0)
2917 0 : nTempSpace += sStat.st_size;
2918 136 : else if (VSIStatL((m_osTempDBFilename + "-wal").c_str(),
2919 68 : &sStat) == 0)
2920 0 : nTempSpace += sStat.st_size;
2921 :
2922 : int nBlockXSize, nBlockYSize;
2923 68 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
2924 68 : const int nBands = IGetRasterCount();
2925 :
2926 68 : if (nTempSpace >
2927 68 : 4 * static_cast<GIntBig>(IGetRasterBand(1)->GetXSize()) *
2928 68 : nBlockYSize * nBands * m_nDTSize)
2929 : {
2930 2 : CPLDebug("GPKG",
2931 : "Partial tiles DB is " CPL_FRMT_GIB
2932 : " bytes. Flushing part of partial tiles",
2933 : nTempSpace);
2934 2 : bTryFreeing = true;
2935 : }
2936 : }
2937 : }
2938 68 : if (bTryFreeing)
2939 : {
2940 2 : if (FlushRemainingShiftedTiles(true /* partial flush*/) != CE_None)
2941 : {
2942 0 : return CE_Failure;
2943 : }
2944 2 : SQLCommand(m_hTempDB,
2945 : "DELETE FROM partial_tiles WHERE zoom_level < 0");
2946 2 : SQLCommand(m_hTempDB, "VACUUM");
2947 : }
2948 : }
2949 4556 : return CE_None;
2950 : }
2951 :
2952 : /************************************************************************/
2953 : /* WriteShiftedTile() */
2954 : /************************************************************************/
2955 :
2956 4556 : CPLErr GDALGPKGMBTilesLikePseudoDataset::WriteShiftedTile(
2957 : int nRow, int nCol, int nBand, int nDstXOffset, int nDstYOffset,
2958 : int nDstXSize, int nDstYSize)
2959 : {
2960 4556 : CPLAssert(m_nShiftXPixelsMod || m_nShiftYPixelsMod);
2961 4556 : CPLAssert(nRow >= 0);
2962 4556 : CPLAssert(nCol >= 0);
2963 4556 : CPLAssert(nRow < m_nTileMatrixHeight);
2964 4556 : CPLAssert(nCol < m_nTileMatrixWidth);
2965 :
2966 4556 : if (m_hTempDB == nullptr &&
2967 55 : (m_poParentDS == nullptr || m_poParentDS->m_hTempDB == nullptr))
2968 : {
2969 : const char *pszBaseFilename =
2970 48 : m_poParentDS ? m_poParentDS->IGetFilename() : IGetFilename();
2971 : m_osTempDBFilename =
2972 48 : CPLResetExtensionSafe(pszBaseFilename, "partial_tiles.db");
2973 48 : CPLPushErrorHandler(CPLQuietErrorHandler);
2974 48 : VSIUnlink(m_osTempDBFilename);
2975 48 : CPLPopErrorHandler();
2976 48 : m_hTempDB = nullptr;
2977 48 : int rc = 0;
2978 48 : if (STARTS_WITH(m_osTempDBFilename, "/vsi"))
2979 : {
2980 33 : m_pMyVFS = OGRSQLiteCreateVFS(nullptr, nullptr);
2981 33 : sqlite3_vfs_register(m_pMyVFS, 0);
2982 33 : rc = sqlite3_open_v2(m_osTempDBFilename, &m_hTempDB,
2983 : SQLITE_OPEN_READWRITE | SQLITE_OPEN_CREATE |
2984 : SQLITE_OPEN_NOMUTEX,
2985 33 : m_pMyVFS->zName);
2986 : }
2987 : else
2988 : {
2989 15 : rc = sqlite3_open(m_osTempDBFilename, &m_hTempDB);
2990 : }
2991 48 : if (rc != SQLITE_OK || m_hTempDB == nullptr)
2992 : {
2993 0 : CPLError(CE_Failure, CPLE_AppDefined,
2994 : "Cannot create temporary database %s",
2995 : m_osTempDBFilename.c_str());
2996 0 : return CE_Failure;
2997 : }
2998 48 : SQLCommand(m_hTempDB, "PRAGMA synchronous = OFF");
2999 48 : SQLCommand(m_hTempDB, "PRAGMA journal_mode = OFF");
3000 48 : SQLCommand(m_hTempDB, "CREATE TABLE partial_tiles("
3001 : "id INTEGER PRIMARY KEY AUTOINCREMENT,"
3002 : "zoom_level INTEGER NOT NULL,"
3003 : "tile_column INTEGER NOT NULL,"
3004 : "tile_row INTEGER NOT NULL,"
3005 : "tile_data_band_1 BLOB,"
3006 : "tile_data_band_2 BLOB,"
3007 : "tile_data_band_3 BLOB,"
3008 : "tile_data_band_4 BLOB,"
3009 : "partial_flag INTEGER NOT NULL,"
3010 : "age INTEGER NOT NULL,"
3011 : "UNIQUE (zoom_level, tile_column, tile_row))");
3012 48 : SQLCommand(m_hTempDB, "CREATE INDEX partial_tiles_partial_flag_idx "
3013 : "ON partial_tiles(partial_flag)");
3014 48 : SQLCommand(m_hTempDB, "CREATE INDEX partial_tiles_age_idx "
3015 : "ON partial_tiles(age)");
3016 :
3017 48 : if (m_poParentDS != nullptr)
3018 : {
3019 4 : m_poParentDS->m_osTempDBFilename = m_osTempDBFilename;
3020 4 : m_poParentDS->m_hTempDB = m_hTempDB;
3021 : }
3022 : }
3023 :
3024 4556 : if (m_poParentDS != nullptr)
3025 4020 : m_hTempDB = m_poParentDS->m_hTempDB;
3026 :
3027 : int nBlockXSize, nBlockYSize;
3028 4556 : IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
3029 4556 : const int nBands = IGetRasterCount();
3030 4556 : const size_t nBandBlockSize =
3031 4556 : static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
3032 :
3033 4556 : int iQuadrantFlag = 0;
3034 4556 : if (nDstXOffset == 0 && nDstYOffset == 0)
3035 1031 : iQuadrantFlag |= (1 << 0);
3036 4556 : if ((nDstXOffset != 0 || nDstXOffset + nDstXSize == nBlockXSize) &&
3037 : nDstYOffset == 0)
3038 1125 : iQuadrantFlag |= (1 << 1);
3039 4556 : if (nDstXOffset == 0 &&
3040 1031 : (nDstYOffset != 0 || nDstYOffset + nDstYSize == nBlockYSize))
3041 1135 : iQuadrantFlag |= (1 << 2);
3042 4556 : if ((nDstXOffset != 0 || nDstXOffset + nDstXSize == nBlockXSize) &&
3043 1125 : (nDstYOffset != 0 || nDstYOffset + nDstYSize == nBlockYSize))
3044 1323 : iQuadrantFlag |= (1 << 3);
3045 4556 : int l_nFlags = iQuadrantFlag << (4 * (nBand - 1));
3046 4556 : int nFullFlags = (1 << (4 * nBands)) - 1;
3047 4556 : int nOldFlags = 0;
3048 :
3049 18224 : for (int i = 1; i <= 3; i++)
3050 : {
3051 13668 : m_asCachedTilesDesc[i].nRow = -1;
3052 13668 : m_asCachedTilesDesc[i].nCol = -1;
3053 13668 : m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
3054 : }
3055 :
3056 4556 : int nExistingId = 0;
3057 4556 : const char *pszSQL = CPLSPrintf(
3058 : "SELECT id, partial_flag, tile_data_band_%d FROM partial_tiles WHERE "
3059 : "zoom_level = %d AND tile_row = %d AND tile_column = %d",
3060 : nBand, m_nZoomLevel, nRow, nCol);
3061 : #ifdef DEBUG_VERBOSE
3062 : CPLDebug("GPKG", "%s", pszSQL);
3063 : #endif
3064 4556 : sqlite3_stmt *hStmt = nullptr;
3065 4556 : int rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
3066 4556 : if (rc != SQLITE_OK)
3067 : {
3068 0 : return CE_Failure;
3069 : }
3070 :
3071 4556 : rc = sqlite3_step(hStmt);
3072 4556 : const int nTileBands = m_eDT == GDT_UInt8 ? 4 : 1;
3073 4556 : GByte *pabyTemp = m_pabyCachedTiles + nTileBands * nBandBlockSize;
3074 4556 : if (rc == SQLITE_ROW)
3075 : {
3076 4134 : nExistingId = sqlite3_column_int(hStmt, 0);
3077 : #ifdef DEBUG_VERBOSE
3078 : CPLDebug("GPKG", "Using partial_tile id=%d", nExistingId);
3079 : #endif
3080 4134 : nOldFlags = sqlite3_column_int(hStmt, 1);
3081 4134 : CPLAssert(nOldFlags != 0);
3082 4134 : if ((nOldFlags & (((1 << 4) - 1) << (4 * (nBand - 1)))) == 0)
3083 : {
3084 1021 : FillEmptyTileSingleBand(pabyTemp + (nBand - 1) * nBandBlockSize);
3085 : }
3086 : else
3087 : {
3088 3113 : CPLAssert(sqlite3_column_bytes(hStmt, 2) ==
3089 : static_cast<int>(nBandBlockSize));
3090 3113 : memcpy(pabyTemp + (nBand - 1) * nBandBlockSize,
3091 : sqlite3_column_blob(hStmt, 2), nBandBlockSize);
3092 : }
3093 : }
3094 : else
3095 : {
3096 422 : FillEmptyTileSingleBand(pabyTemp + (nBand - 1) * nBandBlockSize);
3097 : }
3098 4556 : sqlite3_finalize(hStmt);
3099 4556 : hStmt = nullptr;
3100 :
3101 : /* Copy the updated rectangle into the full tile */
3102 556108 : for (int iY = nDstYOffset; iY < nDstYOffset + nDstYSize; iY++)
3103 : {
3104 551552 : memcpy(pabyTemp +
3105 551552 : (static_cast<size_t>(nBand - 1) * nBlockXSize * nBlockYSize +
3106 551552 : static_cast<size_t>(iY) * nBlockXSize + nDstXOffset) *
3107 551552 : m_nDTSize,
3108 551552 : m_pabyCachedTiles +
3109 551552 : (static_cast<size_t>(nBand - 1) * nBlockXSize * nBlockYSize +
3110 551552 : static_cast<size_t>(iY) * nBlockXSize + nDstXOffset) *
3111 551552 : m_nDTSize,
3112 551552 : static_cast<size_t>(nDstXSize) * m_nDTSize);
3113 : }
3114 :
3115 : #ifdef notdef
3116 : static int nCounter = 1;
3117 : GDALDataset *poLogDS =
3118 : ((GDALDriver *)GDALGetDriverByName("GTiff"))
3119 : ->Create(CPLSPrintf("/tmp/partial_band_%d_%d.tif", 1, nCounter++),
3120 : nBlockXSize, nBlockYSize, nBands, m_eDT, NULL);
3121 : poLogDS->RasterIO(GF_Write, 0, 0, nBlockXSize, nBlockYSize,
3122 : pabyTemp + (nBand - 1) * nBandBlockSize, nBlockXSize,
3123 : nBlockYSize, m_eDT, 1, NULL, 0, 0, 0, NULL);
3124 : GDALClose(poLogDS);
3125 : #endif
3126 :
3127 4556 : if ((nOldFlags & l_nFlags) != 0)
3128 : {
3129 0 : CPLDebug("GPKG",
3130 : "Rewriting quadrant %d of band %d of tile (row=%d,col=%d)",
3131 : iQuadrantFlag, nBand, nRow, nCol);
3132 : }
3133 :
3134 4556 : l_nFlags |= nOldFlags;
3135 4556 : if (l_nFlags == nFullFlags)
3136 : {
3137 : #ifdef DEBUG_VERBOSE
3138 : CPLDebug("GPKG", "Got all quadrants for that tile");
3139 : #endif
3140 1192 : for (int iBand = 1; iBand <= nBands; iBand++)
3141 : {
3142 945 : if (iBand != nBand)
3143 : {
3144 698 : pszSQL = CPLSPrintf(
3145 : "SELECT tile_data_band_%d FROM partial_tiles WHERE "
3146 : "id = %d",
3147 : iBand, nExistingId);
3148 : #ifdef DEBUG_VERBOSE
3149 : CPLDebug("GPKG", "%s", pszSQL);
3150 : #endif
3151 698 : hStmt = nullptr;
3152 : rc =
3153 698 : SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
3154 698 : if (rc != SQLITE_OK)
3155 : {
3156 0 : return CE_Failure;
3157 : }
3158 :
3159 698 : rc = sqlite3_step(hStmt);
3160 698 : if (rc == SQLITE_ROW)
3161 : {
3162 698 : CPLAssert(sqlite3_column_bytes(hStmt, 0) ==
3163 : static_cast<int>(nBandBlockSize));
3164 698 : memcpy(m_pabyCachedTiles + (iBand - 1) * nBandBlockSize,
3165 : sqlite3_column_blob(hStmt, 0), nBandBlockSize);
3166 : }
3167 698 : sqlite3_finalize(hStmt);
3168 698 : hStmt = nullptr;
3169 : }
3170 : else
3171 : {
3172 247 : memcpy(m_pabyCachedTiles + (iBand - 1) * nBandBlockSize,
3173 247 : pabyTemp + (iBand - 1) * nBandBlockSize, nBandBlockSize);
3174 : }
3175 : }
3176 :
3177 247 : m_asCachedTilesDesc[0].nRow = nRow;
3178 247 : m_asCachedTilesDesc[0].nCol = nCol;
3179 247 : m_asCachedTilesDesc[0].nIdxWithinTileData = 0;
3180 247 : m_asCachedTilesDesc[0].abBandDirty[0] = true;
3181 247 : m_asCachedTilesDesc[0].abBandDirty[1] = true;
3182 247 : m_asCachedTilesDesc[0].abBandDirty[2] = true;
3183 247 : m_asCachedTilesDesc[0].abBandDirty[3] = true;
3184 :
3185 494 : pszSQL = CPLSPrintf("UPDATE partial_tiles SET zoom_level = %d, "
3186 : "partial_flag = 0, age = -1 WHERE id = %d",
3187 247 : -1 - m_nZoomLevel, nExistingId);
3188 247 : SQLCommand(m_hTempDB, pszSQL);
3189 : #ifdef DEBUG_VERBOSE
3190 : CPLDebug("GPKG", "%s", pszSQL);
3191 : #endif
3192 247 : CPLErr eErr = WriteTile();
3193 :
3194 : // Call DoPartialFlushOfPartialTilesIfNecessary() after using
3195 : // m_pabyCachedTiles as it is going to mess with it.
3196 247 : if (DoPartialFlushOfPartialTilesIfNecessary() != CE_None)
3197 0 : eErr = CE_None;
3198 :
3199 247 : return eErr;
3200 : }
3201 :
3202 4309 : if (nExistingId == 0)
3203 : {
3204 : OGRErr err;
3205 844 : pszSQL = CPLSPrintf("SELECT id FROM partial_tiles WHERE "
3206 : "partial_flag = 0 AND zoom_level = %d "
3207 : "AND tile_row = %d AND tile_column = %d",
3208 422 : -1 - m_nZoomLevel, nRow, nCol);
3209 : #ifdef DEBUG_VERBOSE
3210 : CPLDebug("GPKG", "%s", pszSQL);
3211 : #endif
3212 422 : nExistingId = SQLGetInteger(m_hTempDB, pszSQL, &err);
3213 422 : if (nExistingId == 0)
3214 : {
3215 418 : pszSQL =
3216 : "SELECT id FROM partial_tiles WHERE partial_flag = 0 LIMIT 1";
3217 : #ifdef DEBUG_VERBOSE
3218 : CPLDebug("GPKG", "%s", pszSQL);
3219 : #endif
3220 418 : nExistingId = SQLGetInteger(m_hTempDB, pszSQL, &err);
3221 : }
3222 : }
3223 :
3224 4309 : const GIntBig nAge = (m_poParentDS) ? m_poParentDS->m_nAge : m_nAge;
3225 4309 : if (nExistingId == 0)
3226 : {
3227 338 : pszSQL = CPLSPrintf(
3228 : "INSERT INTO partial_tiles "
3229 : "(zoom_level, tile_row, tile_column, tile_data_band_%d, "
3230 : "partial_flag, age) VALUES (%d, %d, %d, ?, %d, " CPL_FRMT_GIB ")",
3231 : nBand, m_nZoomLevel, nRow, nCol, l_nFlags, nAge);
3232 : }
3233 : else
3234 : {
3235 3971 : pszSQL = CPLSPrintf(
3236 : "UPDATE partial_tiles SET zoom_level = %d, "
3237 : "tile_row = %d, tile_column = %d, "
3238 : "tile_data_band_%d = ?, partial_flag = %d, age = " CPL_FRMT_GIB
3239 : " WHERE id = %d",
3240 : m_nZoomLevel, nRow, nCol, nBand, l_nFlags, nAge, nExistingId);
3241 : }
3242 4309 : if (m_poParentDS)
3243 3801 : m_poParentDS->m_nAge++;
3244 : else
3245 508 : m_nAge++;
3246 :
3247 : #ifdef DEBUG_VERBOSE
3248 : CPLDebug("GPKG", "%s", pszSQL);
3249 : #endif
3250 :
3251 4309 : hStmt = nullptr;
3252 4309 : rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
3253 4309 : if (rc != SQLITE_OK)
3254 : {
3255 0 : return CE_Failure;
3256 : }
3257 :
3258 4309 : sqlite3_bind_blob(hStmt, 1, pabyTemp + (nBand - 1) * nBandBlockSize,
3259 : static_cast<int>(nBandBlockSize), SQLITE_TRANSIENT);
3260 4309 : rc = sqlite3_step(hStmt);
3261 4309 : CPLErr eErr = CE_Failure;
3262 4309 : if (rc == SQLITE_DONE)
3263 4309 : eErr = CE_None;
3264 : else
3265 : {
3266 0 : CPLError(CE_Failure, CPLE_AppDefined,
3267 : "Failure when inserting partial tile (row=%d,col=%d) at "
3268 : "zoom_level=%d : %s",
3269 : nRow, nCol, m_nZoomLevel, sqlite3_errmsg(m_hTempDB));
3270 : }
3271 :
3272 4309 : sqlite3_finalize(hStmt);
3273 :
3274 : // Call DoPartialFlushOfPartialTilesIfNecessary() after using
3275 : // m_pabyCachedTiles as it is going to mess with it.
3276 4309 : if (DoPartialFlushOfPartialTilesIfNecessary() != CE_None)
3277 0 : eErr = CE_None;
3278 :
3279 4309 : return eErr;
3280 : }
3281 :
3282 : /************************************************************************/
3283 : /* IWriteBlock() */
3284 : /************************************************************************/
3285 :
3286 5851 : CPLErr GDALGPKGMBTilesLikeRasterBand::IWriteBlock(int nBlockXOff,
3287 : int nBlockYOff, void *pData)
3288 : {
3289 : #ifdef DEBUG_VERBOSE
3290 : CPLDebug(
3291 : "GPKG",
3292 : "IWriteBlock(nBand=%d,nBlockXOff=%d,nBlockYOff=%d,m_nZoomLevel=%d)",
3293 : nBand, nBlockXOff, nBlockYOff, m_poTPD->m_nZoomLevel);
3294 : #endif
3295 5851 : if (!m_poTPD->ICanIWriteBlock())
3296 : {
3297 0 : return CE_Failure;
3298 : }
3299 5851 : if (m_poTPD->m_poParentDS)
3300 1139 : m_poTPD->m_poParentDS->m_bHasModifiedTiles = true;
3301 : else
3302 4712 : m_poTPD->m_bHasModifiedTiles = true;
3303 :
3304 5851 : int nRow = nBlockYOff + m_poTPD->m_nShiftYTiles;
3305 5851 : int nCol = nBlockXOff + m_poTPD->m_nShiftXTiles;
3306 :
3307 5851 : int nRowMin = nRow;
3308 5851 : int nRowMax = nRowMin;
3309 5851 : if (m_poTPD->m_nShiftYPixelsMod)
3310 1360 : nRowMax++;
3311 :
3312 5851 : int nColMin = nCol;
3313 5851 : int nColMax = nColMin;
3314 5851 : if (m_poTPD->m_nShiftXPixelsMod)
3315 1322 : nColMax++;
3316 :
3317 5851 : CPLErr eErr = CE_None;
3318 :
3319 13062 : for (nRow = nRowMin; eErr == CE_None && nRow <= nRowMax; nRow++)
3320 : {
3321 17066 : for (nCol = nColMin; eErr == CE_None && nCol <= nColMax; nCol++)
3322 : {
3323 9855 : if (nRow < 0 || nCol < 0 || nRow >= m_poTPD->m_nTileMatrixHeight ||
3324 9715 : nCol >= m_poTPD->m_nTileMatrixWidth)
3325 : {
3326 167 : continue;
3327 : }
3328 :
3329 9688 : if (m_poTPD->m_nShiftXPixelsMod == 0 &&
3330 4517 : m_poTPD->m_nShiftYPixelsMod == 0)
3331 : {
3332 4459 : if (!(nRow == m_poTPD->m_asCachedTilesDesc[0].nRow &&
3333 250 : nCol == m_poTPD->m_asCachedTilesDesc[0].nCol &&
3334 5 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData == 0))
3335 : {
3336 4454 : eErr = m_poTPD->WriteTile();
3337 :
3338 4454 : m_poTPD->m_asCachedTilesDesc[0].nRow = nRow;
3339 4454 : m_poTPD->m_asCachedTilesDesc[0].nCol = nCol;
3340 4454 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData = 0;
3341 : }
3342 : }
3343 :
3344 : // Composite block data into tile, and check if all bands for this
3345 : // block are dirty, and if so write the tile
3346 9688 : bool bAllDirty = true;
3347 42237 : for (int iBand = 1; iBand <= poDS->GetRasterCount(); iBand++)
3348 : {
3349 32549 : GDALRasterBlock *poBlock = nullptr;
3350 32549 : GByte *pabySrc = nullptr;
3351 32549 : if (iBand == nBand)
3352 : {
3353 9688 : pabySrc = static_cast<GByte *>(pData);
3354 : }
3355 : else
3356 : {
3357 22861 : if (!(m_poTPD->m_nShiftXPixelsMod == 0 &&
3358 8359 : m_poTPD->m_nShiftYPixelsMod == 0))
3359 14646 : continue;
3360 :
3361 : // If the block for this band is not dirty, it might be
3362 : // dirty in cache
3363 8215 : if (m_poTPD->m_asCachedTilesDesc[0].abBandDirty[iBand - 1])
3364 505 : continue;
3365 : else
3366 : {
3367 : poBlock =
3368 7710 : cpl::down_cast<GDALGPKGMBTilesLikeRasterBand *>(
3369 7710 : poDS->GetRasterBand(iBand))
3370 7710 : ->TryGetLockedBlockRef(nBlockXOff, nBlockYOff);
3371 7710 : if (poBlock && poBlock->GetDirty())
3372 : {
3373 : pabySrc =
3374 7670 : static_cast<GByte *>(poBlock->GetDataRef());
3375 7670 : poBlock->MarkClean();
3376 : }
3377 : else
3378 : {
3379 40 : if (poBlock)
3380 2 : poBlock->DropLock();
3381 40 : bAllDirty = false;
3382 40 : continue;
3383 : }
3384 : }
3385 : }
3386 :
3387 17358 : if (m_poTPD->m_nShiftXPixelsMod == 0 &&
3388 12187 : m_poTPD->m_nShiftYPixelsMod == 0)
3389 12129 : m_poTPD->m_asCachedTilesDesc[0].abBandDirty[iBand - 1] =
3390 : true;
3391 :
3392 17358 : int nDstXOffset = 0;
3393 17358 : int nDstXSize = nBlockXSize;
3394 17358 : int nDstYOffset = 0;
3395 17358 : int nDstYSize = nBlockYSize;
3396 : // Composite block data into tile data
3397 17358 : if (m_poTPD->m_nShiftXPixelsMod == 0 &&
3398 12187 : m_poTPD->m_nShiftYPixelsMod == 0)
3399 : {
3400 :
3401 : #ifdef DEBUG_VERBOSE
3402 : if (eDataType == GDT_UInt8 &&
3403 : nBlockXOff * nBlockXSize <=
3404 : nRasterXSize - nBlockXSize &&
3405 : nBlockYOff * nBlockYSize > nRasterYSize - nBlockYSize)
3406 : {
3407 : bool bFoundNonZero = false;
3408 : for (int y = nRasterYSize - nBlockYOff * nBlockYSize;
3409 : y < nBlockYSize; y++)
3410 : {
3411 : for (int x = 0; x < nBlockXSize; x++)
3412 : {
3413 : if (pabySrc[static_cast<GPtrDiff_t>(y) *
3414 : nBlockXSize +
3415 : x] != 0 &&
3416 : !bFoundNonZero)
3417 : {
3418 : CPLDebug("GPKG",
3419 : "IWriteBlock(): Found non-zero "
3420 : "content in ghost part of "
3421 : "tile(nBand=%d,nBlockXOff=%d,"
3422 : "nBlockYOff=%d,m_nZoomLevel=%d)\n",
3423 : iBand, nBlockXOff, nBlockYOff,
3424 : m_poTPD->m_nZoomLevel);
3425 : bFoundNonZero = true;
3426 : }
3427 : }
3428 : }
3429 : }
3430 : #endif
3431 :
3432 12129 : const size_t nBandBlockSize =
3433 12129 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
3434 12129 : m_nDTSize;
3435 12129 : memcpy(m_poTPD->m_pabyCachedTiles +
3436 12129 : (iBand - 1) * nBandBlockSize,
3437 : pabySrc, nBandBlockSize);
3438 :
3439 : // Make sure partial blocks are zero'ed outside of the
3440 : // validity area but do that only when know that JPEG will
3441 : // not be used so as to avoid edge effects (although we
3442 : // should probably repeat last pixels if we really want to
3443 : // do that, but that only makes sense if readers only clip
3444 : // to the gpkg_contents extent). Well, ere on the safe side
3445 : // for now
3446 12129 : if (m_poTPD->m_eTF != GPKG_TF_JPEG &&
3447 11881 : (nBlockXOff * nBlockXSize >=
3448 11881 : nRasterXSize - nBlockXSize ||
3449 11338 : nBlockYOff * nBlockYSize >=
3450 11338 : nRasterYSize - nBlockYSize))
3451 : {
3452 1039 : int nXEndValidity =
3453 1039 : nRasterXSize - nBlockXOff * nBlockXSize;
3454 1039 : if (nXEndValidity > nBlockXSize)
3455 496 : nXEndValidity = nBlockXSize;
3456 1039 : int nYEndValidity =
3457 1039 : nRasterYSize - nBlockYOff * nBlockYSize;
3458 1039 : if (nYEndValidity > nBlockYSize)
3459 294 : nYEndValidity = nBlockYSize;
3460 1039 : if (nXEndValidity < nBlockXSize)
3461 : {
3462 9820 : for (int iY = 0; iY < nYEndValidity; iY++)
3463 : {
3464 9623 : m_poTPD->FillBuffer(
3465 9623 : m_poTPD->m_pabyCachedTiles +
3466 9623 : ((static_cast<size_t>(iBand - 1) *
3467 9623 : nBlockYSize +
3468 9623 : iY) *
3469 9623 : nBlockXSize +
3470 9623 : nXEndValidity) *
3471 9623 : m_nDTSize,
3472 9623 : nBlockXSize - nXEndValidity);
3473 : }
3474 : }
3475 1039 : if (nYEndValidity < nBlockYSize)
3476 : {
3477 237 : m_poTPD->FillBuffer(
3478 237 : m_poTPD->m_pabyCachedTiles +
3479 237 : (static_cast<size_t>(iBand - 1) *
3480 237 : nBlockYSize +
3481 237 : nYEndValidity) *
3482 237 : nBlockXSize * m_nDTSize,
3483 237 : static_cast<size_t>(nBlockYSize -
3484 : nYEndValidity) *
3485 237 : nBlockXSize);
3486 : }
3487 12129 : }
3488 : }
3489 : else
3490 : {
3491 5229 : const int nXValid =
3492 5229 : (nBlockXOff * nBlockXSize > nRasterXSize - nBlockXSize)
3493 5229 : ? (nRasterXSize - nBlockXOff * nBlockXSize)
3494 : : nBlockXSize;
3495 5229 : const int nYValid =
3496 5229 : (nBlockYOff * nBlockYSize > nRasterYSize - nBlockYSize)
3497 5229 : ? (nRasterYSize - nBlockYOff * nBlockYSize)
3498 : : nBlockYSize;
3499 :
3500 5229 : int nSrcXOffset = 0;
3501 5229 : if (nCol == nColMin)
3502 : {
3503 2643 : nDstXOffset = m_poTPD->m_nShiftXPixelsMod;
3504 2643 : nDstXSize = std::min(
3505 2643 : nXValid, nBlockXSize - m_poTPD->m_nShiftXPixelsMod);
3506 : }
3507 : else
3508 : {
3509 2586 : nDstXOffset = 0;
3510 2586 : if (nXValid > nBlockXSize - m_poTPD->m_nShiftXPixelsMod)
3511 : {
3512 2212 : nDstXSize = nXValid - (nBlockXSize -
3513 2212 : m_poTPD->m_nShiftXPixelsMod);
3514 : }
3515 : else
3516 374 : nDstXSize = 0;
3517 2586 : nSrcXOffset = nBlockXSize - m_poTPD->m_nShiftXPixelsMod;
3518 : }
3519 :
3520 5229 : int nSrcYOffset = 0;
3521 5229 : if (nRow == nRowMin)
3522 : {
3523 2607 : nDstYOffset = m_poTPD->m_nShiftYPixelsMod;
3524 2607 : nDstYSize = std::min(
3525 2607 : nYValid, nBlockYSize - m_poTPD->m_nShiftYPixelsMod);
3526 : }
3527 : else
3528 : {
3529 2622 : nDstYOffset = 0;
3530 2622 : if (nYValid > nBlockYSize - m_poTPD->m_nShiftYPixelsMod)
3531 : {
3532 2232 : nDstYSize = nYValid - (nBlockYSize -
3533 2232 : m_poTPD->m_nShiftYPixelsMod);
3534 : }
3535 : else
3536 390 : nDstYSize = 0;
3537 2622 : nSrcYOffset = nBlockYSize - m_poTPD->m_nShiftYPixelsMod;
3538 : }
3539 :
3540 : #ifdef DEBUG_VERBOSE
3541 : CPLDebug("GPKG",
3542 : "Copy source tile x=%d,w=%d,y=%d,h=%d into "
3543 : "buffer at x=%d,y=%d",
3544 : nDstXOffset, nDstXSize, nDstYOffset, nDstYSize,
3545 : nSrcXOffset, nSrcYOffset);
3546 : #endif
3547 :
3548 5229 : if (nDstXSize > 0 && nDstYSize > 0)
3549 : {
3550 556108 : for (GPtrDiff_t y = 0; y < nDstYSize; y++)
3551 : {
3552 551552 : GByte *pDst = m_poTPD->m_pabyCachedTiles +
3553 551552 : (static_cast<size_t>(iBand - 1) *
3554 551552 : nBlockXSize * nBlockYSize +
3555 551552 : (y + nDstYOffset) * nBlockXSize +
3556 551552 : nDstXOffset) *
3557 551552 : m_nDTSize;
3558 551552 : GByte *pSrc =
3559 551552 : pabySrc + ((y + nSrcYOffset) * nBlockXSize +
3560 551552 : nSrcXOffset) *
3561 551552 : m_nDTSize;
3562 551552 : GDALCopyWords(pSrc, eDataType, m_nDTSize, pDst,
3563 : eDataType, m_nDTSize, nDstXSize);
3564 : }
3565 : }
3566 : }
3567 :
3568 17358 : if (poBlock)
3569 7670 : poBlock->DropLock();
3570 :
3571 17358 : if (!(m_poTPD->m_nShiftXPixelsMod == 0 &&
3572 12187 : m_poTPD->m_nShiftYPixelsMod == 0))
3573 : {
3574 5229 : m_poTPD->m_asCachedTilesDesc[0].nRow = -1;
3575 5229 : m_poTPD->m_asCachedTilesDesc[0].nCol = -1;
3576 5229 : m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData = -1;
3577 5229 : if (nDstXSize > 0 && nDstYSize > 0)
3578 : {
3579 4556 : eErr = m_poTPD->WriteShiftedTile(
3580 : nRow, nCol, iBand, nDstXOffset, nDstYOffset,
3581 : nDstXSize, nDstYSize);
3582 : }
3583 : }
3584 : }
3585 :
3586 9688 : if (m_poTPD->m_nShiftXPixelsMod == 0 &&
3587 4517 : m_poTPD->m_nShiftYPixelsMod == 0)
3588 : {
3589 4459 : if (bAllDirty)
3590 : {
3591 4437 : eErr = m_poTPD->WriteTile();
3592 : }
3593 : }
3594 : }
3595 : }
3596 :
3597 5851 : return eErr;
3598 : }
3599 :
3600 : /************************************************************************/
3601 : /* GetNoDataValue() */
3602 : /************************************************************************/
3603 :
3604 19954 : double GDALGPKGMBTilesLikeRasterBand::GetNoDataValue(int *pbSuccess)
3605 : {
3606 19954 : if (m_bHasNoData)
3607 : {
3608 402 : if (pbSuccess)
3609 401 : *pbSuccess = TRUE;
3610 402 : return m_dfNoDataValue;
3611 : }
3612 19552 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
3613 : }
3614 :
3615 : /************************************************************************/
3616 : /* SetNoDataValueInternal() */
3617 : /************************************************************************/
3618 :
3619 75 : void GDALGPKGMBTilesLikeRasterBand::SetNoDataValueInternal(double dfNoDataValue)
3620 : {
3621 75 : m_bHasNoData = true;
3622 75 : m_dfNoDataValue = dfNoDataValue;
3623 75 : }
3624 :
3625 : /************************************************************************/
3626 : /* GDALGeoPackageRasterBand() */
3627 : /************************************************************************/
3628 :
3629 1817 : GDALGeoPackageRasterBand::GDALGeoPackageRasterBand(
3630 1817 : GDALGeoPackageDataset *poDSIn, int nTileWidth, int nTileHeight)
3631 1817 : : GDALGPKGMBTilesLikeRasterBand(poDSIn, nTileWidth, nTileHeight)
3632 : {
3633 1817 : poDS = poDSIn;
3634 1817 : }
3635 :
3636 : /************************************************************************/
3637 : /* GetOverviewCount() */
3638 : /************************************************************************/
3639 :
3640 8 : int GDALGeoPackageRasterBand::GetOverviewCount()
3641 : {
3642 : GDALGeoPackageDataset *poGDS =
3643 8 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3644 8 : return static_cast<int>(poGDS->m_apoOverviewDS.size());
3645 : }
3646 :
3647 : /************************************************************************/
3648 : /* GetOverviewCount() */
3649 : /************************************************************************/
3650 :
3651 41 : GDALRasterBand *GDALGeoPackageRasterBand::GetOverview(int nIdx)
3652 : {
3653 : GDALGeoPackageDataset *poGDS =
3654 41 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3655 41 : if (nIdx < 0 || nIdx >= static_cast<int>(poGDS->m_apoOverviewDS.size()))
3656 1 : return nullptr;
3657 40 : return poGDS->m_apoOverviewDS[nIdx]->GetRasterBand(nBand);
3658 : }
3659 :
3660 : /************************************************************************/
3661 : /* SetNoDataValue() */
3662 : /************************************************************************/
3663 :
3664 23 : CPLErr GDALGeoPackageRasterBand::SetNoDataValue(double dfNoDataValue)
3665 : {
3666 : GDALGeoPackageDataset *poGDS =
3667 23 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3668 :
3669 23 : if (eDataType == GDT_UInt8)
3670 : {
3671 6 : if (!(dfNoDataValue >= 0 && dfNoDataValue <= 255 &&
3672 4 : static_cast<int>(dfNoDataValue) == dfNoDataValue))
3673 : {
3674 2 : CPLError(CE_Failure, CPLE_NotSupported,
3675 : "Invalid nodata value for a Byte band: %.17g",
3676 : dfNoDataValue);
3677 2 : return CE_Failure;
3678 : }
3679 :
3680 8 : for (int i = 1; i <= poGDS->nBands; ++i)
3681 : {
3682 5 : if (i != nBand)
3683 : {
3684 2 : int bHasNoData = FALSE;
3685 : double dfOtherNoDataValue =
3686 2 : poGDS->GetRasterBand(i)->GetNoDataValue(&bHasNoData);
3687 2 : if (bHasNoData && dfOtherNoDataValue != dfNoDataValue)
3688 : {
3689 1 : CPLError(
3690 : CE_Failure, CPLE_NotSupported,
3691 : "Only the same nodata value can be set on all bands");
3692 1 : return CE_Failure;
3693 : }
3694 : }
3695 : }
3696 :
3697 3 : SetNoDataValueInternal(dfNoDataValue);
3698 3 : poGDS->m_bMetadataDirty = true;
3699 3 : return CE_None;
3700 : }
3701 :
3702 17 : if (std::isnan(dfNoDataValue))
3703 : {
3704 0 : CPLError(CE_Warning, CPLE_NotSupported,
3705 : "A NaN nodata value cannot be recorded in "
3706 : "gpkg_2d_gridded_coverage_ancillary table");
3707 : }
3708 :
3709 17 : SetNoDataValueInternal(dfNoDataValue);
3710 :
3711 17 : char *pszSQL = sqlite3_mprintf(
3712 : "UPDATE gpkg_2d_gridded_coverage_ancillary SET data_null = ? "
3713 : "WHERE tile_matrix_set_name = '%q'",
3714 : poGDS->m_osRasterTable.c_str());
3715 17 : sqlite3_stmt *hStmt = nullptr;
3716 17 : int rc = SQLPrepareWithError(poGDS->IGetDB(), pszSQL, -1, &hStmt, nullptr);
3717 17 : if (rc == SQLITE_OK)
3718 : {
3719 17 : if (poGDS->m_eTF == GPKG_TF_PNG_16BIT)
3720 : {
3721 15 : if (eDataType == GDT_UInt16 && poGDS->m_dfOffset == 0.0 &&
3722 6 : poGDS->m_dfScale == 1.0 && dfNoDataValue >= 0 &&
3723 5 : dfNoDataValue <= 65535 &&
3724 5 : static_cast<GUInt16>(dfNoDataValue) == dfNoDataValue)
3725 : {
3726 5 : poGDS->m_usGPKGNull = static_cast<GUInt16>(dfNoDataValue);
3727 : }
3728 : else
3729 : {
3730 10 : poGDS->m_usGPKGNull = 65535;
3731 : }
3732 15 : sqlite3_bind_double(hStmt, 1, poGDS->m_usGPKGNull);
3733 : }
3734 : else
3735 : {
3736 2 : sqlite3_bind_double(hStmt, 1, static_cast<float>(dfNoDataValue));
3737 : }
3738 17 : rc = sqlite3_step(hStmt);
3739 17 : sqlite3_finalize(hStmt);
3740 : }
3741 17 : sqlite3_free(pszSQL);
3742 :
3743 17 : return (rc == SQLITE_OK) ? CE_None : CE_Failure;
3744 : }
3745 :
3746 : /************************************************************************/
3747 : /* LoadBandMetadata() */
3748 : /************************************************************************/
3749 :
3750 1245 : void GDALGeoPackageRasterBand::LoadBandMetadata()
3751 : {
3752 : GDALGeoPackageDataset *poGDS =
3753 1245 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3754 :
3755 1245 : if (m_bHasReadMetadataFromStorage)
3756 1143 : return;
3757 :
3758 477 : m_bHasReadMetadataFromStorage = true;
3759 :
3760 477 : poGDS->TryLoadXML();
3761 :
3762 477 : if (!poGDS->HasMetadataTables())
3763 375 : return;
3764 :
3765 102 : char *pszSQL = sqlite3_mprintf(
3766 : "SELECT md.metadata, md.md_standard_uri, md.mime_type "
3767 : "FROM gpkg_metadata md "
3768 : "JOIN gpkg_metadata_reference mdr ON (md.id = mdr.md_file_id ) "
3769 : "WHERE "
3770 : "mdr.reference_scope = 'table' AND lower(mdr.table_name) = "
3771 : "lower('%q') ORDER BY md.id "
3772 : "LIMIT 1000", // to avoid denial of service
3773 : poGDS->m_osRasterTable.c_str());
3774 :
3775 102 : auto oResult = SQLQuery(poGDS->hDB, pszSQL);
3776 102 : sqlite3_free(pszSQL);
3777 102 : if (!oResult)
3778 : {
3779 0 : return;
3780 : }
3781 :
3782 : /* GDAL metadata */
3783 200 : for (int i = 0; i < oResult->RowCount(); i++)
3784 : {
3785 98 : const char *pszMetadata = oResult->GetValue(0, i);
3786 98 : const char *pszMDStandardURI = oResult->GetValue(1, i);
3787 98 : const char *pszMimeType = oResult->GetValue(2, i);
3788 98 : if (pszMetadata && pszMDStandardURI && pszMimeType &&
3789 98 : EQUAL(pszMDStandardURI, "http://gdal.org") &&
3790 74 : EQUAL(pszMimeType, "text/xml"))
3791 : {
3792 74 : CPLXMLNode *psXMLNode = CPLParseXMLString(pszMetadata);
3793 74 : if (psXMLNode)
3794 : {
3795 148 : GDALMultiDomainMetadata oLocalMDMD;
3796 74 : oLocalMDMD.XMLInit(psXMLNode, FALSE);
3797 :
3798 74 : CSLConstList papszDomainList = oLocalMDMD.GetDomainList();
3799 208 : for (CSLConstList papszIter = papszDomainList;
3800 208 : papszIter && *papszIter; ++papszIter)
3801 : {
3802 134 : if (STARTS_WITH(*papszIter, "BAND_"))
3803 : {
3804 5 : int l_nBand = atoi(*papszIter + strlen("BAND_"));
3805 5 : if (l_nBand >= 1 && l_nBand <= poGDS->GetRasterCount())
3806 : {
3807 : auto l_poBand =
3808 5 : cpl::down_cast<GDALGeoPackageRasterBand *>(
3809 : poGDS->GetRasterBand(l_nBand));
3810 5 : l_poBand->m_bHasReadMetadataFromStorage = true;
3811 :
3812 10 : char **papszMD = CSLDuplicate(
3813 5 : oLocalMDMD.GetMetadata(*papszIter));
3814 10 : papszMD = CSLMerge(
3815 : papszMD,
3816 5 : GDALGPKGMBTilesLikeRasterBand::GetMetadata(""));
3817 5 : l_poBand->GDALPamRasterBand::SetMetadata(papszMD);
3818 5 : CSLDestroy(papszMD);
3819 : }
3820 : }
3821 : }
3822 :
3823 74 : CPLDestroyXMLNode(psXMLNode);
3824 : }
3825 : }
3826 : }
3827 : }
3828 :
3829 : /************************************************************************/
3830 : /* GetMetadata() */
3831 : /************************************************************************/
3832 :
3833 958 : char **GDALGeoPackageRasterBand::GetMetadata(const char *pszDomain)
3834 : {
3835 : GDALGeoPackageDataset *poGDS =
3836 958 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3837 958 : LoadBandMetadata(); /* force loading from storage if needed */
3838 :
3839 15 : if (poGDS->eAccess == GA_ReadOnly && eDataType != GDT_UInt8 &&
3840 11 : (pszDomain == nullptr || EQUAL(pszDomain, "")) &&
3841 11 : !m_bMinMaxComputedFromTileAncillary &&
3842 981 : !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MINIMUM") &&
3843 8 : !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MAXIMUM"))
3844 : {
3845 8 : m_bMinMaxComputedFromTileAncillary = true;
3846 :
3847 8 : const int nColMin = poGDS->m_nShiftXTiles;
3848 8 : const int nColMax =
3849 8 : (nRasterXSize - 1 + poGDS->m_nShiftXPixelsMod) / nBlockXSize +
3850 8 : poGDS->m_nShiftXTiles;
3851 8 : const int nRowMin = poGDS->m_nShiftYTiles;
3852 8 : const int nRowMax =
3853 8 : (nRasterYSize - 1 + poGDS->m_nShiftYPixelsMod) / nBlockYSize +
3854 8 : poGDS->m_nShiftYTiles;
3855 :
3856 8 : bool bOK = false;
3857 8 : if (poGDS->m_nShiftXPixelsMod == 0 && poGDS->m_nShiftYPixelsMod == 0 &&
3858 8 : (nRasterXSize % nBlockXSize) == 0 &&
3859 2 : (nRasterYSize % nBlockYSize) == 0)
3860 : {
3861 : // If the area of interest matches entire tiles, then we can
3862 : // use tile statistics
3863 2 : bOK = true;
3864 : }
3865 6 : else if (m_bHasNoData)
3866 : {
3867 : // Otherwise, in the case where we have nodata, we assume that
3868 : // if the area of interest is at least larger than the existing
3869 : // tiles, the tile statistics will be reliable.
3870 2 : char *pszSQL = sqlite3_mprintf(
3871 : "SELECT MIN(tile_column), MAX(tile_column), "
3872 : "MIN(tile_row), MAX(tile_row) FROM \"%w\" "
3873 : "WHERE zoom_level = %d",
3874 : poGDS->m_osRasterTable.c_str(), poGDS->m_nZoomLevel);
3875 4 : auto sResult = SQLQuery(poGDS->IGetDB(), pszSQL);
3876 2 : if (sResult && sResult->RowCount() == 1)
3877 : {
3878 2 : const char *pszMinX = sResult->GetValue(0, 0);
3879 2 : const char *pszMaxX = sResult->GetValue(1, 0);
3880 2 : const char *pszMinY = sResult->GetValue(2, 0);
3881 2 : const char *pszMaxY = sResult->GetValue(3, 0);
3882 2 : if (pszMinX && pszMaxX && pszMinY && pszMaxY)
3883 : {
3884 2 : bOK = atoi(pszMinX) >= nColMin &&
3885 2 : atoi(pszMaxX) <= nColMax &&
3886 4 : atoi(pszMinY) >= nRowMin && atoi(pszMaxY) <= nRowMax;
3887 : }
3888 : }
3889 2 : sqlite3_free(pszSQL);
3890 : }
3891 :
3892 8 : if (bOK)
3893 : {
3894 4 : char *pszSQL = sqlite3_mprintf(
3895 : "SELECT MIN(min), MAX(max) FROM "
3896 : "gpkg_2d_gridded_tile_ancillary WHERE tpudt_id "
3897 : "IN (SELECT id FROM \"%w\" WHERE "
3898 : "zoom_level = %d AND "
3899 : "tile_column >= %d AND tile_column <= %d AND "
3900 : "tile_row >= %d AND tile_row <= %d)",
3901 : poGDS->m_osRasterTable.c_str(), poGDS->m_nZoomLevel, nColMin,
3902 : nColMax, nRowMin, nRowMax);
3903 8 : auto sResult = SQLQuery(poGDS->IGetDB(), pszSQL);
3904 4 : CPLDebug("GPKG", "%s", pszSQL);
3905 4 : if (sResult && sResult->RowCount() == 1)
3906 : {
3907 4 : const char *pszMin = sResult->GetValue(0, 0);
3908 4 : const char *pszMax = sResult->GetValue(1, 0);
3909 4 : if (pszMin)
3910 : {
3911 4 : m_dfStatsMinFromTileAncillary = CPLAtof(pszMin);
3912 : }
3913 4 : if (pszMax)
3914 : {
3915 4 : m_dfStatsMaxFromTileAncillary = CPLAtof(pszMax);
3916 : }
3917 : }
3918 4 : sqlite3_free(pszSQL);
3919 : }
3920 : }
3921 :
3922 709 : if (m_bAddImplicitStatistics && m_bMinMaxComputedFromTileAncillary &&
3923 8 : (pszDomain == nullptr || EQUAL(pszDomain, "")) &&
3924 1675 : !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MINIMUM") &&
3925 8 : !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MAXIMUM"))
3926 : {
3927 : m_aosMD.Assign(CSLDuplicate(
3928 8 : GDALGPKGMBTilesLikeRasterBand::GetMetadata(pszDomain)));
3929 8 : if (!std::isnan(m_dfStatsMinFromTileAncillary))
3930 : {
3931 : m_aosMD.SetNameValue(
3932 : "STATISTICS_MINIMUM",
3933 4 : CPLSPrintf("%.14g", m_dfStatsMinFromTileAncillary));
3934 : }
3935 8 : if (!std::isnan(m_dfStatsMaxFromTileAncillary))
3936 : {
3937 : m_aosMD.SetNameValue(
3938 : "STATISTICS_MAXIMUM",
3939 4 : CPLSPrintf("%.14g", m_dfStatsMaxFromTileAncillary));
3940 : }
3941 8 : return m_aosMD.List();
3942 : }
3943 :
3944 950 : return GDALGPKGMBTilesLikeRasterBand::GetMetadata(pszDomain);
3945 : }
3946 :
3947 : /************************************************************************/
3948 : /* GetMetadataItem() */
3949 : /************************************************************************/
3950 :
3951 251 : const char *GDALGeoPackageRasterBand::GetMetadataItem(const char *pszName,
3952 : const char *pszDomain)
3953 : {
3954 251 : LoadBandMetadata(); /* force loading from storage if needed */
3955 :
3956 251 : if (m_bAddImplicitStatistics && eDataType != GDT_UInt8 &&
3957 7 : (pszDomain == nullptr || EQUAL(pszDomain, "")) &&
3958 4 : (EQUAL(pszName, "STATISTICS_MINIMUM") ||
3959 2 : EQUAL(pszName, "STATISTICS_MAXIMUM")))
3960 : {
3961 2 : return CSLFetchNameValue(GetMetadata(), pszName);
3962 : }
3963 :
3964 249 : return GDALGPKGMBTilesLikeRasterBand::GetMetadataItem(pszName, pszDomain);
3965 : }
3966 :
3967 : /************************************************************************/
3968 : /* SetMetadata() */
3969 : /************************************************************************/
3970 :
3971 16 : CPLErr GDALGeoPackageRasterBand::SetMetadata(char **papszMetadata,
3972 : const char *pszDomain)
3973 : {
3974 : GDALGeoPackageDataset *poGDS =
3975 16 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3976 16 : LoadBandMetadata(); /* force loading from storage if needed */
3977 16 : poGDS->m_bMetadataDirty = true;
3978 72 : for (CSLConstList papszIter = papszMetadata; papszIter && *papszIter;
3979 : ++papszIter)
3980 : {
3981 56 : if (STARTS_WITH(*papszIter, "STATISTICS_"))
3982 51 : m_bStatsMetadataSetInThisSession = true;
3983 : }
3984 16 : return GDALPamRasterBand::SetMetadata(papszMetadata, pszDomain);
3985 : }
3986 :
3987 : /************************************************************************/
3988 : /* SetMetadataItem() */
3989 : /************************************************************************/
3990 :
3991 20 : CPLErr GDALGeoPackageRasterBand::SetMetadataItem(const char *pszName,
3992 : const char *pszValue,
3993 : const char *pszDomain)
3994 : {
3995 : GDALGeoPackageDataset *poGDS =
3996 20 : cpl::down_cast<GDALGeoPackageDataset *>(poDS);
3997 20 : LoadBandMetadata(); /* force loading from storage if needed */
3998 20 : poGDS->m_bMetadataDirty = true;
3999 20 : if (STARTS_WITH(pszName, "STATISTICS_"))
4000 20 : m_bStatsMetadataSetInThisSession = true;
4001 20 : return GDALPamRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
4002 : }
4003 :
4004 : /************************************************************************/
4005 : /* InvalidateStatistics() */
4006 : /************************************************************************/
4007 :
4008 346 : void GDALGeoPackageRasterBand::InvalidateStatistics()
4009 : {
4010 346 : bool bModified = false;
4011 692 : CPLStringList aosMD(CSLDuplicate(GetMetadata()));
4012 367 : for (CSLConstList papszIter = GetMetadata(); papszIter && *papszIter;
4013 : ++papszIter)
4014 : {
4015 21 : if (STARTS_WITH(*papszIter, "STATISTICS_"))
4016 : {
4017 20 : char *pszKey = nullptr;
4018 20 : CPLParseNameValue(*papszIter, &pszKey);
4019 20 : CPLAssert(pszKey);
4020 20 : aosMD.SetNameValue(pszKey, nullptr);
4021 20 : CPLFree(pszKey);
4022 20 : bModified = true;
4023 : }
4024 : }
4025 346 : if (bModified)
4026 4 : SetMetadata(aosMD.List());
4027 346 : }
|
