LCOV - code coverage report
Current view: top level - ogr/ogrsf_frmts/gpkg - gdalgeopackagerasterband.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1870 2011 93.0 %
Date: 2025-05-11 14:50:45 Functions: 46 47 97.9 %

          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        2649 : GDALGPKGMBTilesLikePseudoDataset::GDALGPKGMBTilesLikePseudoDataset()
      35             :     : m_bForceTempDBCompaction(
      36        2649 :           CPLTestBool(CPLGetConfigOption("GPKG_FORCE_TEMPDB_COMPACTION", "NO")))
      37             : {
      38       13245 :     for (int i = 0; i < 4; i++)
      39             :     {
      40       10596 :         m_asCachedTilesDesc[i].nRow = -1;
      41       10596 :         m_asCachedTilesDesc[i].nCol = -1;
      42       10596 :         m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
      43       10596 :         m_asCachedTilesDesc[i].abBandDirty[0] = FALSE;
      44       10596 :         m_asCachedTilesDesc[i].abBandDirty[1] = FALSE;
      45       10596 :         m_asCachedTilesDesc[i].abBandDirty[2] = FALSE;
      46       10596 :         m_asCachedTilesDesc[i].abBandDirty[3] = FALSE;
      47             :     }
      48        2649 : }
      49             : 
      50             : /************************************************************************/
      51             : /*                 ~GDALGPKGMBTilesLikePseudoDataset()                  */
      52             : /************************************************************************/
      53             : 
      54        2649 : GDALGPKGMBTilesLikePseudoDataset::~GDALGPKGMBTilesLikePseudoDataset()
      55             : {
      56        2649 :     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        2649 :     CPLFree(m_pabyCachedTiles);
      69        2649 :     delete m_poCT;
      70        2649 :     CPLFree(m_pabyHugeColorArray);
      71        2649 : }
      72             : 
      73             : /************************************************************************/
      74             : /*                            SetDataType()                             */
      75             : /************************************************************************/
      76             : 
      77         319 : void GDALGPKGMBTilesLikePseudoDataset::SetDataType(GDALDataType eDT)
      78             : {
      79         319 :     CPLAssert(eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16 ||
      80             :               eDT == GDT_Float32);
      81         319 :     m_eDT = eDT;
      82         319 :     m_nDTSize = GDALGetDataTypeSizeBytes(m_eDT);
      83         319 : }
      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        2551 : GDALGPKGMBTilesLikeRasterBand::GDALGPKGMBTilesLikeRasterBand(
     108        2551 :     GDALGPKGMBTilesLikePseudoDataset *poTPD, int nTileWidth, int nTileHeight)
     109        2551 :     : m_poTPD(poTPD)
     110             : {
     111        2551 :     eDataType = m_poTPD->m_eDT;
     112        2551 :     m_nDTSize = m_poTPD->m_nDTSize;
     113        2551 :     nBlockXSize = nTileWidth;
     114        2551 :     nBlockYSize = nTileHeight;
     115        2551 : }
     116             : 
     117             : #ifdef __GNUC__
     118             : #pragma GCC diagnostic pop
     119             : #endif
     120             : 
     121             : /************************************************************************/
     122             : /*                              FlushCache()                            */
     123             : /************************************************************************/
     124             : 
     125        5477 : CPLErr GDALGPKGMBTilesLikeRasterBand::FlushCache(bool bAtClosing)
     126             : {
     127        5477 :     m_poTPD->m_nLastSpaceCheckTimestamp = -1;  // disable partial flushes
     128        5477 :     CPLErr eErr = GDALPamRasterBand::FlushCache(bAtClosing);
     129        5477 :     if (eErr == CE_None)
     130        5474 :         eErr = m_poTPD->IFlushCacheWithErrCode(bAtClosing);
     131        5477 :     m_poTPD->m_nLastSpaceCheckTimestamp = 0;
     132        5477 :     return eErr;
     133             : }
     134             : 
     135             : /************************************************************************/
     136             : /*                              FlushTiles()                            */
     137             : /************************************************************************/
     138             : 
     139        3383 : CPLErr GDALGPKGMBTilesLikePseudoDataset::FlushTiles()
     140             : {
     141        3383 :     CPLErr eErr = CE_None;
     142        3383 :     GDALGPKGMBTilesLikePseudoDataset *poMainDS =
     143        3383 :         m_poParentDS ? m_poParentDS : this;
     144        3383 :     if (poMainDS->m_nTileInsertionCount < 0)
     145           7 :         return CE_Failure;
     146             : 
     147        3376 :     if (IGetUpdate())
     148             :     {
     149        2017 :         if (m_nShiftXPixelsMod || m_nShiftYPixelsMod)
     150             :         {
     151         607 :             eErr = FlushRemainingShiftedTiles(false /* total flush*/);
     152             :         }
     153             :         else
     154             :         {
     155        1410 :             eErr = WriteTile();
     156             :         }
     157             :     }
     158             : 
     159        3376 :     if (poMainDS->m_nTileInsertionCount > 0)
     160             :     {
     161         188 :         if (poMainDS->ICommitTransaction() != OGRERR_NONE)
     162             :         {
     163           0 :             poMainDS->m_nTileInsertionCount = -1;
     164           0 :             eErr = CE_Failure;
     165             :         }
     166             :         else
     167             :         {
     168         188 :             poMainDS->m_nTileInsertionCount = 0;
     169             :         }
     170             :     }
     171        3376 :     return eErr;
     172             : }
     173             : 
     174             : /************************************************************************/
     175             : /*                             GetColorTable()                          */
     176             : /************************************************************************/
     177             : 
     178         345 : GDALColorTable *GDALGPKGMBTilesLikeRasterBand::GetColorTable()
     179             : {
     180         345 :     if (poDS->GetRasterCount() != 1)
     181         127 :         return nullptr;
     182             : 
     183         218 :     if (!m_poTPD->m_bTriedEstablishingCT)
     184             :     {
     185          64 :         m_poTPD->m_bTriedEstablishingCT = true;
     186          64 :         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          68 :         for (int i = 0; i < 2; i++)
     196             :         {
     197          62 :             bool bRetry = false;
     198          62 :             char *pszSQL = nullptr;
     199          62 :             if (i == 0)
     200             :             {
     201          56 :                 pszSQL = sqlite3_mprintf("SELECT tile_data FROM \"%w\" "
     202             :                                          "WHERE zoom_level = %d LIMIT 1",
     203          56 :                                          m_poTPD->m_osRasterTable.c_str(),
     204          56 :                                          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          62 :             sqlite3_stmt *hStmt = nullptr;
     220          62 :             int rc = SQLPrepareWithError(m_poTPD->IGetDB(), pszSQL, -1, &hStmt,
     221             :                                          nullptr);
     222          62 :             if (rc == SQLITE_OK)
     223             :             {
     224          62 :                 rc = sqlite3_step(hStmt);
     225          79 :                 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          62 :             sqlite3_free(pszSQL);
     264          62 :             sqlite3_finalize(hStmt);
     265          62 :             if (!bRetry)
     266          50 :                 break;
     267             :         }
     268             :     }
     269             : 
     270         210 :     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_Byte)
     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         209 : GDALColorInterp GDALGPKGMBTilesLikeRasterBand::GetColorInterpretation()
     317             : {
     318         209 :     if (m_poTPD->m_eDT != GDT_Byte)
     319          22 :         return GCI_Undefined;
     320         187 :     if (poDS->GetRasterCount() == 1)
     321          49 :         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       14254 : void GDALGPKGMBTilesLikePseudoDataset::FillBuffer(GByte *pabyData,
     384             :                                                   size_t nPixels)
     385             : {
     386       14254 :     int bHasNoData = FALSE;
     387       14254 :     const double dfNoDataValue = IGetRasterBand(1)->GetNoDataValue(&bHasNoData);
     388       14254 :     if (!bHasNoData || dfNoDataValue == 0.0)
     389             :     {
     390             :         // cppcheck-suppress nullPointer
     391       13952 :         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       14254 : }
     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        2582 : CPLErr GDALGPKGMBTilesLikePseudoDataset::ReadTile(
     431             :     const CPLString &osMemFileName, GByte *pabyTileData, double dfTileOffset,
     432             :     double dfTileScale, bool *pbIsLossyFormat)
     433             : {
     434        2582 :     const char *const apszDriversByte[] = {"JPEG", "PNG", "WEBP", nullptr};
     435        2582 :     const char *const apszDriversInt[] = {"PNG", nullptr};
     436        2582 :     const char *const apszDriversFloat[] = {"GTiff", nullptr};
     437             :     int nBlockXSize, nBlockYSize;
     438        2582 :     IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
     439        2582 :     const int nBands = IGetRasterCount();
     440             :     auto poDSTile = std::unique_ptr<GDALDataset>(GDALDataset::Open(
     441             :         osMemFileName.c_str(), GDAL_OF_RASTER | GDAL_OF_INTERNAL,
     442        2582 :         (m_eDT == GDT_Byte)                   ? apszDriversByte
     443          43 :         : (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT) ? apszDriversFloat
     444             :                                               : apszDriversInt,
     445        5164 :         nullptr, nullptr));
     446        2582 :     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        2580 :     const int nTileBandCount = poDSTile->GetRasterCount();
     454             : 
     455        5160 :     if (!(poDSTile->GetRasterXSize() == nBlockXSize &&
     456        2580 :           poDSTile->GetRasterYSize() == nBlockYSize &&
     457        7740 :           (nTileBandCount >= 1 && nTileBandCount <= 4)) ||
     458        2580 :         (m_eDT != GDT_Byte && 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        2580 :     GDALDataType eRequestDT = GDT_Byte;
     467        2580 :     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        2540 :     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        2580 :     if (poDSTile->RasterIO(GF_Read, 0, 0, nBlockXSize, nBlockYSize,
     480             :                            pabyTileData, nBlockXSize, nBlockYSize, eRequestDT,
     481             :                            poDSTile->GetRasterCount(), nullptr, 0, 0, 0,
     482        2580 :                            nullptr) != CE_None)
     483             :     {
     484           0 :         FillEmptyTile(pabyTileData);
     485           0 :         return CE_Failure;
     486             :     }
     487             : 
     488        2580 :     if (m_eDT != GDT_Byte)
     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        2537 :     GDALColorTable *poCT = nullptr;
     567        2537 :     if (nBands == 1 || nTileBandCount == 1)
     568             :     {
     569         153 :         poCT = poDSTile->GetRasterBand(1)->GetColorTable();
     570         153 :         IGetRasterBand(1)->GetColorTable();
     571             :     }
     572             : 
     573        2537 :     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        2537 :     const GPtrDiff_t nBlockPixels =
     580        2537 :         static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize;
     581        2537 :     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          31 :     if (nBands == 1 && nTileBandCount == 1 && poCT != nullptr &&
     627        2567 :         m_poCT != nullptr && !poCT->IsSame(m_poCT))
     628             :     {
     629           0 :         CPLError(CE_Warning, CPLE_NotSupported,
     630             :                  "Different color tables. Unhandled for now");
     631             :     }
     632        2536 :     else if ((nBands == 1 && nTileBandCount >= 3) ||
     633          31 :              (nBands == 1 && nTileBandCount == 1 && m_poCT != nullptr &&
     634        2536 :               poCT == nullptr) ||
     635        2536 :              ((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        2536 :     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        2229 :     else if (nTileBandCount == 2)
     658             :     {
     659             :         /* Do Grey+Alpha -> RGBA */
     660         362 :         memcpy(pabyTileData + 3 * nBlockPixels, pabyTileData + 1 * nBlockPixels,
     661             :                nBlockPixels);
     662         362 :         memcpy(pabyTileData + 1 * nBlockPixels, pabyTileData, nBlockPixels);
     663         362 :         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        2536 :     return CE_None;
     713             : }
     714             : 
     715             : /************************************************************************/
     716             : /*                           ReadTile()                                 */
     717             : /************************************************************************/
     718             : 
     719        7542 : GByte *GDALGPKGMBTilesLikePseudoDataset::ReadTile(int nRow, int nCol)
     720             : {
     721             :     int nBlockXSize, nBlockYSize;
     722        7542 :     IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
     723        7542 :     const int nBands = IGetRasterCount();
     724        7542 :     const size_t nBandBlockSize =
     725        7542 :         static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
     726        7542 :     const int nTileBands = m_eDT == GDT_Byte ? 4 : 1;
     727        7542 :     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        2739 :         GByte *pabyDest = m_pabyCachedTiles + 2 * nTileBands * nBandBlockSize;
     779        2739 :         bool bAllNonDirty = true;
     780       11227 :         for (int i = 0; i < nBands; i++)
     781             :         {
     782        8488 :             if (m_asCachedTilesDesc[0].abBandDirty[i])
     783           2 :                 bAllNonDirty = false;
     784             :         }
     785        2739 :         if (bAllNonDirty)
     786             :         {
     787        2737 :             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        2582 : void GDALGPKGMBTilesLikePseudoDataset::GetTileOffsetAndScale(
     813             :     GIntBig nTileId, double &dfTileOffset, double &dfTileScale)
     814             : {
     815        2582 :     dfTileOffset = 0.0;
     816        2582 :     dfTileScale = 1.0;
     817             : 
     818        2582 :     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        2582 : }
     842             : 
     843             : /************************************************************************/
     844             : /*                           ReadTile()                                 */
     845             : /************************************************************************/
     846             : 
     847        5520 : GByte *GDALGPKGMBTilesLikePseudoDataset::ReadTile(int nRow, int nCol,
     848             :                                                   GByte *pabyData,
     849             :                                                   bool *pbIsLossyFormat)
     850             : {
     851        5520 :     int nBlockXSize = 0;
     852        5520 :     int nBlockYSize = 0;
     853        5520 :     IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
     854        5520 :     const int nBands = IGetRasterCount();
     855             : 
     856        5520 :     if (pbIsLossyFormat)
     857          18 :         *pbIsLossyFormat = false;
     858             : 
     859        5520 :     const size_t nBandBlockSize =
     860        5520 :         static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
     861        5520 :     if (nRow < 0 || nCol < 0 || nRow >= m_nTileMatrixHeight ||
     862        5463 :         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       16341 :     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        5447 :         m_eDT != GDT_Byte ? ", id" : "",  // MBTiles do not have an id
     876             :         m_osRasterTable.c_str(), m_nZoomLevel,
     877        5447 :         GetRowFromIntoTopConvention(nRow), nCol,
     878        5447 :         !m_osWHERE.empty() ? CPLSPrintf(" AND (%s)", m_osWHERE.c_str()) : "");
     879             : 
     880             : #ifdef DEBUG_VERBOSE
     881             :     CPLDebug("GPKG", "%s", pszSQL);
     882             : #endif
     883             : 
     884        5447 :     sqlite3_stmt *hStmt = nullptr;
     885        5447 :     int rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt, nullptr);
     886        5447 :     sqlite3_free(pszSQL);
     887        5447 :     if (rc != SQLITE_OK)
     888             :     {
     889           0 :         return nullptr;
     890             :     }
     891        5447 :     rc = sqlite3_step(hStmt);
     892             : 
     893        5447 :     if (rc == SQLITE_ROW && sqlite3_column_type(hStmt, 0) == SQLITE_BLOB)
     894             :     {
     895        2565 :         const int nBytes = sqlite3_column_bytes(hStmt, 0);
     896             :         GIntBig nTileId =
     897        2565 :             (m_eDT == GDT_Byte) ? 0 : sqlite3_column_int64(hStmt, 1);
     898             :         GByte *pabyRawData = static_cast<GByte *>(
     899        2565 :             const_cast<void *>(sqlite3_column_blob(hStmt, 0)));
     900             :         const CPLString osMemFileName(
     901        5130 :             VSIMemGenerateHiddenFilename("gpkg_read_tile"));
     902        2565 :         VSILFILE *fp = VSIFileFromMemBuffer(osMemFileName.c_str(), pabyRawData,
     903             :                                             nBytes, FALSE);
     904        2565 :         VSIFCloseL(fp);
     905             : 
     906        2565 :         double dfTileOffset = 0.0;
     907        2565 :         double dfTileScale = 1.0;
     908        2565 :         GetTileOffsetAndScale(nTileId, dfTileOffset, dfTileScale);
     909        2565 :         ReadTile(osMemFileName, pabyData, dfTileOffset, dfTileScale,
     910             :                  pbIsLossyFormat);
     911        2565 :         VSIUnlink(osMemFileName);
     912        2565 :         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        5447 :     return pabyData;
     981             : }
     982             : 
     983             : /************************************************************************/
     984             : /*                         IReadBlock()                                 */
     985             : /************************************************************************/
     986             : 
     987        3950 : 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        3950 :     if (m_poTPD->m_pabyCachedTiles == nullptr)
     997           0 :         return CE_Failure;
     998             : 
     999        3950 :     const int nRowMin = nBlockYOff + m_poTPD->m_nShiftYTiles;
    1000        3950 :     int nRowMax = nRowMin;
    1001        3950 :     if (m_poTPD->m_nShiftYPixelsMod)
    1002        1211 :         nRowMax++;
    1003             : 
    1004        3950 :     const int nColMin = nBlockXOff + m_poTPD->m_nShiftXTiles;
    1005        3950 :     int nColMax = nColMin;
    1006        3950 :     if (m_poTPD->m_nShiftXPixelsMod)
    1007        1189 :         nColMax++;
    1008             : 
    1009        2761 : retry:
    1010             :     /* Optimize for left to right reading at constant row */
    1011        3953 :     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        9114 :     for (int nRow = nRowMin; nRow <= nRowMax; nRow++)
    1050             :     {
    1051       12703 :         for (int nCol = nColMin; nCol <= nColMax; nCol++)
    1052             :         {
    1053        7542 :             if (m_poTPD->m_nShiftXPixelsMod == 0 &&
    1054        2786 :                 m_poTPD->m_nShiftYPixelsMod == 0)
    1055             :             {
    1056        2739 :                 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        2737 :                     if (m_poTPD->WriteTile() != CE_None)
    1061           0 :                         return CE_Failure;
    1062             :                 }
    1063             :             }
    1064             : 
    1065        7542 :             GByte *pabyTileData = m_poTPD->ReadTile(nRow, nCol);
    1066        7542 :             if (pabyTileData == nullptr)
    1067           0 :                 return CE_Failure;
    1068             : 
    1069       34600 :             for (int iBand = 1; iBand <= poDS->GetRasterCount(); iBand++)
    1070             :             {
    1071       27061 :                 GDALRasterBlock *poBlock = nullptr;
    1072       27061 :                 GByte *pabyDest = nullptr;
    1073       27061 :                 if (iBand == nBand)
    1074             :                 {
    1075        7542 :                     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       27056 :                 if (m_poTPD->m_nShiftXPixelsMod == 0 &&
    1109        8640 :                     m_poTPD->m_nShiftYPixelsMod == 0)
    1110             :                 {
    1111        8486 :                     const size_t nBandBlockSize =
    1112        8486 :                         static_cast<size_t>(nBlockXSize) * nBlockYSize *
    1113        8486 :                         m_nDTSize;
    1114        8486 :                     memcpy(pabyDest,
    1115        8486 :                            pabyTileData + (iBand - 1) * nBandBlockSize,
    1116        8486 :                            nBandBlockSize);
    1117             : #ifdef DEBUG_VERBOSE
    1118             :                     if (eDataType == GDT_Byte &&
    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       27056 :                 if (poBlock)
    1206       19514 :                     poBlock->DropLock();
    1207             :             }
    1208             :         }
    1209             :     }
    1210             : 
    1211        3950 :     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_row, tile_column 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;  // (row, col)
    1269          12 :     while (rc == SQLITE_ROW)
    1270             :     {
    1271           0 :         oSetTiles.insert(std::pair(
    1272           4 :             m_poTPD->GetRowFromIntoTopConvention(sqlite3_column_int(hStmt, 0)),
    1273           8 :             sqlite3_column_int(hStmt, 1)));
    1274           4 :         rc = sqlite3_step(hStmt);
    1275             :     }
    1276           8 :     sqlite3_finalize(hStmt);
    1277           8 :     if (rc != SQLITE_DONE)
    1278             :     {
    1279             :         return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
    1280           0 :                GDAL_DATA_COVERAGE_STATUS_DATA;
    1281             :     }
    1282           8 :     if (oSetTiles.empty())
    1283             :     {
    1284           4 :         if (pdfDataPct)
    1285           4 :             *pdfDataPct = 0.0;
    1286           4 :         return GDAL_DATA_COVERAGE_STATUS_EMPTY;
    1287             :     }
    1288             : 
    1289           8 :     if (m_poTPD->m_nShiftXPixelsMod == 0 && m_poTPD->m_nShiftYPixelsMod == 0 &&
    1290           4 :         oSetTiles.size() == static_cast<size_t>(iRowMax - iRowMin + 1) *
    1291           4 :                                 (iColMax - iColMin + 1))
    1292             :     {
    1293           2 :         if (pdfDataPct)
    1294           2 :             *pdfDataPct = 100.0;
    1295           2 :         return GDAL_DATA_COVERAGE_STATUS_DATA;
    1296             :     }
    1297             : 
    1298           2 :     if (m_poTPD->m_nShiftXPixelsMod == 0 && m_poTPD->m_nShiftYPixelsMod == 0)
    1299             :     {
    1300           2 :         int nStatus = 0;
    1301           2 :         GIntBig nPixelsData = 0;
    1302           7 :         for (int iY = iRowMin; iY <= iRowMax; ++iY)
    1303             :         {
    1304          21 :             for (int iX = iColMin; iX <= iColMax; ++iX)
    1305             :             {
    1306             :                 // coverity[swapped_arguments]
    1307          16 :                 if (oSetTiles.find(std::pair(iY, iX)) == 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       13456 : CPLErr GDALGPKGMBTilesLikePseudoDataset::WriteTile()
    1555             : {
    1556       13456 :     GDALGPKGMBTilesLikePseudoDataset *poMainDS =
    1557       13456 :         m_poParentDS ? m_poParentDS : this;
    1558       13456 :     if (poMainDS->m_nTileInsertionCount < 0)
    1559         501 :         return CE_Failure;
    1560             : 
    1561       12955 :     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       12955 :     GDALRasterBlock::EnterDisableDirtyBlockFlush();
    1571       12955 :     m_bInWriteTile = true;
    1572       12955 :     CPLErr eErr = WriteTileInternal();
    1573       12955 :     m_bInWriteTile = false;
    1574       12955 :     GDALRasterBlock::LeaveDisableDirtyBlockFlush();
    1575       12955 :     return eErr;
    1576             : }
    1577             : 
    1578             : /* should only be called by WriteTile() */
    1579       12955 : CPLErr GDALGPKGMBTilesLikePseudoDataset::WriteTileInternal()
    1580             : {
    1581       17579 :     if (!(IGetUpdate() && m_asCachedTilesDesc[0].nRow >= 0 &&
    1582        4624 :           m_asCachedTilesDesc[0].nCol >= 0 &&
    1583        4624 :           m_asCachedTilesDesc[0].nIdxWithinTileData == 0))
    1584        8331 :         return CE_None;
    1585             : 
    1586        4624 :     int nRow = m_asCachedTilesDesc[0].nRow;
    1587        4624 :     int nCol = m_asCachedTilesDesc[0].nCol;
    1588             : 
    1589        4624 :     bool bAllDirty = true;
    1590        4624 :     bool bAllNonDirty = true;
    1591        4624 :     const int nBands = IGetRasterCount();
    1592       18006 :     for (int i = 0; i < nBands; i++)
    1593             :     {
    1594       13382 :         if (m_asCachedTilesDesc[0].abBandDirty[i])
    1595       13349 :             bAllNonDirty = false;
    1596             :         else
    1597          33 :             bAllDirty = false;
    1598             :     }
    1599        4624 :     if (bAllNonDirty)
    1600           0 :         return CE_None;
    1601             : 
    1602             :     int nBlockXSize, nBlockYSize;
    1603        4624 :     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        4624 :     bool bIsLossyFormat = false;
    1608        4624 :     const size_t nBandBlockSize =
    1609        4624 :         static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
    1610        4624 :     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_Byte ? 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        4624 :     int nXOff = (nCol - m_nShiftXTiles) * nBlockXSize - m_nShiftXPixelsMod;
    1633        4624 :     int nYOff = (nRow - m_nShiftYTiles) * nBlockYSize - m_nShiftYPixelsMod;
    1634             : 
    1635             :     /* Assert that the tile at least intersects some of the GDAL raster space */
    1636        4624 :     CPLAssert(nXOff > -nBlockXSize);
    1637        4624 :     CPLAssert(nYOff > -nBlockYSize);
    1638             :     /* Can happen if the tile of the raster is less than the block size */
    1639        4624 :     const int nRasterXSize = IGetRasterBand(1)->GetXSize();
    1640        4624 :     const int nRasterYSize = IGetRasterBand(1)->GetYSize();
    1641        4624 :     if (nXOff >= nRasterXSize || nYOff >= nRasterYSize)
    1642           0 :         return CE_None;
    1643             : 
    1644             : #ifdef DEBUG_VERBOSE
    1645             :     if (m_nShiftXPixelsMod == 0 && m_nShiftYPixelsMod == 0 && m_eDT == GDT_Byte)
    1646             :     {
    1647             :         int nBlockXOff = nCol;
    1648             :         int nBlockYOff = nRow;
    1649             :         if (nBlockXOff * nBlockXSize <= nRasterXSize - nBlockXSize &&
    1650             :             nBlockYOff * nBlockYSize > nRasterYSize - nBlockYSize)
    1651             :         {
    1652             :             for (int i = 0; i < nBands; i++)
    1653             :             {
    1654             :                 bool bFoundNonZero = false;
    1655             :                 for (GPtrDiff_t y = nRasterYSize - nBlockYOff * nBlockYSize;
    1656             :                      y < nBlockYSize; y++)
    1657             :                 {
    1658             :                     for (int x = 0; x < nBlockXSize; x++)
    1659             :                     {
    1660             :                         if (m_pabyCachedTiles[y * nBlockXSize + x +
    1661             :                                               i * nBandBlockSize] != 0 &&
    1662             :                             !bFoundNonZero)
    1663             :                         {
    1664             :                             CPLDebug("GPKG",
    1665             :                                      "WriteTileInternal(): Found non-zero "
    1666             :                                      "content in ghost part of "
    1667             :                                      "tile(band=%d,nBlockXOff=%d,nBlockYOff=%d,"
    1668             :                                      "m_nZoomLevel=%d)\n",
    1669             :                                      i + 1, nBlockXOff, nBlockYOff,
    1670             :                                      m_nZoomLevel);
    1671             :                             bFoundNonZero = true;
    1672             :                         }
    1673             :                     }
    1674             :                 }
    1675             :             }
    1676             :         }
    1677             :     }
    1678             : #endif
    1679             : 
    1680             :     /* Validity area of tile data in intra-tile coordinate space */
    1681        4624 :     int iXOff = 0;
    1682        4624 :     GPtrDiff_t iYOff = 0;
    1683        4624 :     int iXCount = nBlockXSize;
    1684        4624 :     int iYCount = nBlockYSize;
    1685             : 
    1686        4624 :     bool bPartialTile = false;
    1687        4624 :     int nAlphaBand = (nBands == 2) ? 2 : (nBands == 4) ? 4 : 0;
    1688        4624 :     if (nAlphaBand == 0)
    1689             :     {
    1690        3742 :         if (nXOff < 0)
    1691             :         {
    1692           4 :             bPartialTile = true;
    1693           4 :             iXOff = -nXOff;
    1694           4 :             iXCount += nXOff;
    1695             :         }
    1696        3742 :         if (nXOff > nRasterXSize - nBlockXSize)
    1697             :         {
    1698          93 :             bPartialTile = true;
    1699          93 :             iXCount -= static_cast<int>(static_cast<GIntBig>(nXOff) +
    1700          93 :                                         nBlockXSize - nRasterXSize);
    1701             :         }
    1702        3742 :         if (nYOff < 0)
    1703             :         {
    1704           8 :             bPartialTile = true;
    1705           8 :             iYOff = -nYOff;
    1706           8 :             iYCount += nYOff;
    1707             :         }
    1708        3742 :         if (nYOff > nRasterYSize - nBlockYSize)
    1709             :         {
    1710         107 :             bPartialTile = true;
    1711         107 :             iYCount -= static_cast<int>(static_cast<GIntBig>(nYOff) +
    1712         107 :                                         nBlockYSize - nRasterYSize);
    1713             :         }
    1714        3742 :         CPLAssert(iXOff >= 0);
    1715        3742 :         CPLAssert(iYOff >= 0);
    1716        3742 :         CPLAssert(iXCount > 0);
    1717        3742 :         CPLAssert(iYCount > 0);
    1718        3742 :         CPLAssert(iXOff + iXCount <= nBlockXSize);
    1719        3742 :         CPLAssert(iYOff + iYCount <= nBlockYSize);
    1720             :     }
    1721             : 
    1722        4624 :     m_asCachedTilesDesc[0].nRow = -1;
    1723        4624 :     m_asCachedTilesDesc[0].nCol = -1;
    1724        4624 :     m_asCachedTilesDesc[0].nIdxWithinTileData = -1;
    1725        4624 :     m_asCachedTilesDesc[0].abBandDirty[0] = false;
    1726        4624 :     m_asCachedTilesDesc[0].abBandDirty[1] = false;
    1727        4624 :     m_asCachedTilesDesc[0].abBandDirty[2] = false;
    1728        4624 :     m_asCachedTilesDesc[0].abBandDirty[3] = false;
    1729             : 
    1730        4624 :     CPLErr eErr = CE_Failure;
    1731             : 
    1732        4624 :     int bHasNoData = FALSE;
    1733        4624 :     double dfNoDataValue = IGetRasterBand(1)->GetNoDataValue(&bHasNoData);
    1734        4624 :     const bool bHasNanNoData = bHasNoData && std::isnan(dfNoDataValue);
    1735             : 
    1736        4624 :     bool bAllOpaque = true;
    1737             :     // Detect fully transparent tiles, but only if all bands are dirty (that is
    1738             :     // the user wrote content into them) !
    1739        4624 :     if (bAllDirty && m_eDT == GDT_Byte && m_poCT == nullptr && nAlphaBand != 0)
    1740             :     {
    1741         872 :         GByte byFirstAlphaVal =
    1742         872 :             m_pabyCachedTiles[(nAlphaBand - 1) * nBlockXSize * nBlockYSize];
    1743         872 :         GPtrDiff_t i = 1;
    1744    24105900 :         for (; i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize; i++)
    1745             :         {
    1746    24105400 :             if (m_pabyCachedTiles[static_cast<GPtrDiff_t>(nAlphaBand - 1) *
    1747    24105400 :                                       nBlockXSize * nBlockYSize +
    1748    24105400 :                                   i] != byFirstAlphaVal)
    1749         367 :                 break;
    1750             :         }
    1751         872 :         if (i == static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize)
    1752             :         {
    1753             :             // If tile is fully transparent, don't serialize it and remove it if
    1754             :             // it exists
    1755         505 :             if (byFirstAlphaVal == 0)
    1756             :             {
    1757         440 :                 DeleteTile(nRow, nCol);
    1758             : 
    1759         440 :                 return CE_None;
    1760             :             }
    1761          65 :             bAllOpaque = (byFirstAlphaVal == 255);
    1762             :         }
    1763             :         else
    1764         432 :             bAllOpaque = false;
    1765             :     }
    1766        3752 :     else if (bAllDirty && m_eDT == GDT_Byte && m_poCT == nullptr &&
    1767        3654 :              (!bHasNoData || dfNoDataValue == 0.0))
    1768             :     {
    1769        3654 :         bool bAllEmpty = true;
    1770        3654 :         const auto nPixels =
    1771        3654 :             static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize * nBands;
    1772     1207860 :         for (GPtrDiff_t i = 0; i < nPixels; i++)
    1773             :         {
    1774     1207840 :             if (m_pabyCachedTiles[i] != 0)
    1775             :             {
    1776        3637 :                 bAllEmpty = false;
    1777        3637 :                 break;
    1778             :             }
    1779             :         }
    1780        3654 :         if (bAllEmpty)
    1781             :         {
    1782             :             // If tile is fully transparent, don't serialize it and remove it if
    1783             :             // it exists
    1784          17 :             DeleteTile(nRow, nCol);
    1785             : 
    1786          17 :             return CE_None;
    1787        3637 :         }
    1788             :     }
    1789          98 :     else if (bAllDirty && m_eDT == GDT_Float32)
    1790             :     {
    1791          12 :         const float *pSrc = reinterpret_cast<float *>(m_pabyCachedTiles);
    1792             :         GPtrDiff_t i;
    1793          12 :         const float fNoDataValueOrZero =
    1794          12 :             bHasNoData ? static_cast<float>(dfNoDataValue) : 0.0f;
    1795      131084 :         for (i = 0; i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize; i++)
    1796             :         {
    1797      131082 :             const float fVal = pSrc[i];
    1798      131082 :             if (bHasNanNoData)
    1799             :             {
    1800           0 :                 if (std::isnan(fVal))
    1801           0 :                     continue;
    1802             :             }
    1803      131082 :             else if (fVal == fNoDataValueOrZero)
    1804             :             {
    1805      131072 :                 continue;
    1806             :             }
    1807          10 :             break;
    1808             :         }
    1809          12 :         if (i == static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize)
    1810             :         {
    1811             :             // If tile is fully transparent, don't serialize it and remove it if
    1812             :             // it exists
    1813           2 :             DeleteTile(nRow, nCol);
    1814             : 
    1815           2 :             return CE_None;
    1816             :         }
    1817             :     }
    1818             : 
    1819        4165 :     if (bIsLossyFormat)
    1820             :     {
    1821           1 :         CPLDebug("GPKG",
    1822             :                  "Had to read tile (row=%d,col=%d) at zoom_level=%d, "
    1823             :                  "stored in a lossy format, before rewriting it, causing "
    1824             :                  "potential extra quality loss",
    1825             :                  nRow, nCol, m_nZoomLevel);
    1826             :     }
    1827             : 
    1828             :     const CPLString osMemFileName(
    1829        8330 :         VSIMemGenerateHiddenFilename("gpkg_write_tile"));
    1830        4165 :     const char *pszDriverName = "PNG";
    1831        4165 :     CPL_IGNORE_RET_VAL(pszDriverName);  // Make CSA happy
    1832        4165 :     bool bTileDriverSupports1Band = false;
    1833        4165 :     bool bTileDriverSupports2Bands = false;
    1834        4165 :     bool bTileDriverSupports4Bands = false;
    1835        4165 :     bool bTileDriverSupportsCT = false;
    1836             : 
    1837        4165 :     if (nBands == 1 && m_eDT == GDT_Byte)
    1838          76 :         IGetRasterBand(1)->GetColorTable();
    1839             : 
    1840        4165 :     GDALDataType eTileDT = GDT_Byte;
    1841             :     // If not all bands are dirty, then (temporarily) use a lossless format
    1842        4165 :     if (m_eTF == GPKG_TF_PNG_JPEG || (!bAllDirty && m_eTF == GPKG_TF_JPEG))
    1843             :     {
    1844          97 :         bTileDriverSupports1Band = true;
    1845          97 :         if (bPartialTile || !bAllDirty || (nBands == 2 && !bAllOpaque) ||
    1846          19 :             (nBands == 4 && !bAllOpaque) || m_poCT != nullptr)
    1847             :         {
    1848          89 :             pszDriverName = "PNG";
    1849          89 :             bTileDriverSupports2Bands = m_bPNGSupports2Bands;
    1850          89 :             bTileDriverSupports4Bands = true;
    1851          89 :             bTileDriverSupportsCT = m_bPNGSupportsCT;
    1852             :         }
    1853             :         else
    1854           8 :             pszDriverName = "JPEG";
    1855             :     }
    1856        4068 :     else if (m_eTF == GPKG_TF_PNG || m_eTF == GPKG_TF_PNG8)
    1857             :     {
    1858        3845 :         pszDriverName = "PNG";
    1859        3845 :         bTileDriverSupports1Band = true;
    1860        3845 :         bTileDriverSupports2Bands = m_bPNGSupports2Bands;
    1861        3845 :         bTileDriverSupports4Bands = true;
    1862        3845 :         bTileDriverSupportsCT = m_bPNGSupportsCT;
    1863             :     }
    1864         223 :     else if (m_eTF == GPKG_TF_JPEG)
    1865             :     {
    1866          82 :         pszDriverName = "JPEG";
    1867          82 :         bTileDriverSupports1Band = true;
    1868             :     }
    1869         141 :     else if (m_eTF == GPKG_TF_WEBP)
    1870             :     {
    1871          83 :         pszDriverName = "WEBP";
    1872          83 :         bTileDriverSupports4Bands = WEBPSupports4Bands();
    1873             :     }
    1874          58 :     else if (m_eTF == GPKG_TF_PNG_16BIT)
    1875             :     {
    1876          53 :         pszDriverName = "PNG";
    1877          53 :         eTileDT = GDT_UInt16;
    1878          53 :         bTileDriverSupports1Band = true;
    1879             :     }
    1880           5 :     else if (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
    1881             :     {
    1882           5 :         pszDriverName = "GTiff";
    1883           5 :         eTileDT = GDT_Float32;
    1884           5 :         bTileDriverSupports1Band = true;
    1885             :     }
    1886             :     else
    1887             :     {
    1888           0 :         CPLAssert(false);
    1889             :     }
    1890             : 
    1891             :     GDALDriver *l_poDriver =
    1892        4165 :         GDALDriver::FromHandle(GDALGetDriverByName(pszDriverName));
    1893        4165 :     if (l_poDriver != nullptr)
    1894             :     {
    1895        4164 :         auto poMEMDS = MEMDataset::Create("", nBlockXSize, nBlockYSize, 0,
    1896             :                                           eTileDT, nullptr);
    1897        4164 :         int nTileBands = nBands;
    1898        4164 :         if (bPartialTile && nBands == 1 && m_poCT == nullptr &&
    1899             :             bTileDriverSupports2Bands)
    1900          31 :             nTileBands = 2;
    1901        4133 :         else if (bPartialTile && bTileDriverSupports4Bands)
    1902          40 :             nTileBands = 4;
    1903             :         // only use (somewhat lossy) PNG8 if all bands are dirty
    1904        4093 :         else if (bAllDirty && m_eTF == GPKG_TF_PNG8 && nBands >= 3 &&
    1905          10 :                  bAllOpaque && !bPartialTile)
    1906          10 :             nTileBands = 1;
    1907        4083 :         else if (nBands == 2)
    1908             :         {
    1909         383 :             if (bAllOpaque)
    1910             :             {
    1911          51 :                 if (bTileDriverSupports2Bands)
    1912          30 :                     nTileBands = 1;
    1913             :                 else
    1914          21 :                     nTileBands = 3;
    1915             :             }
    1916         332 :             else if (!bTileDriverSupports2Bands)
    1917             :             {
    1918         135 :                 if (bTileDriverSupports4Bands)
    1919          68 :                     nTileBands = 4;
    1920             :                 else
    1921          67 :                     nTileBands = 3;
    1922             :             }
    1923             :         }
    1924        3700 :         else if (nBands == 4 && (bAllOpaque || !bTileDriverSupports4Bands))
    1925          21 :             nTileBands = 3;
    1926        3679 :         else if (nBands == 1 && m_poCT != nullptr && !bTileDriverSupportsCT)
    1927             :         {
    1928           1 :             nTileBands = 3;
    1929           1 :             if (bTileDriverSupports4Bands)
    1930             :             {
    1931           0 :                 for (int i = 0; i < m_poCT->GetColorEntryCount(); i++)
    1932             :                 {
    1933           0 :                     const GDALColorEntry *psEntry = m_poCT->GetColorEntry(i);
    1934           0 :                     if (psEntry->c4 == 0)
    1935             :                     {
    1936           0 :                         nTileBands = 4;
    1937           0 :                         break;
    1938             :                     }
    1939             :                 }
    1940           1 :             }
    1941             :         }
    1942        3678 :         else if (nBands == 1 && m_poCT == nullptr && !bTileDriverSupports1Band)
    1943           1 :             nTileBands = 3;
    1944             : 
    1945        4164 :         if (bPartialTile && (nTileBands == 2 || nTileBands == 4))
    1946             :         {
    1947          71 :             int nTargetAlphaBand = nTileBands;
    1948          71 :             memset(m_pabyCachedTiles + (nTargetAlphaBand - 1) * nBandBlockSize,
    1949             :                    0, nBandBlockSize);
    1950        4356 :             for (GPtrDiff_t iY = iYOff; iY < iYOff + iYCount; iY++)
    1951             :             {
    1952        4285 :                 memset(m_pabyCachedTiles +
    1953        4285 :                            (static_cast<size_t>(nTargetAlphaBand - 1) *
    1954        4285 :                                 nBlockYSize +
    1955        4285 :                             iY) *
    1956        4285 :                                nBlockXSize +
    1957        4285 :                            iXOff,
    1958             :                        255, iXCount);
    1959             :             }
    1960             :         }
    1961             : 
    1962        4164 :         GUInt16 *pTempTileBuffer = nullptr;
    1963        4164 :         GPtrDiff_t nValidPixels = 0;
    1964        4164 :         double dfTileMin = 0.0;
    1965        4164 :         double dfTileMax = 0.0;
    1966        4164 :         double dfTileMean = 0.0;
    1967        4164 :         double dfTileStdDev = 0.0;
    1968        4164 :         double dfTileOffset = 0.0;
    1969        4164 :         double dfTileScale = 1.0;
    1970        4164 :         if (m_eTF == GPKG_TF_PNG_16BIT)
    1971             :         {
    1972             :             pTempTileBuffer = static_cast<GUInt16 *>(
    1973          53 :                 VSI_MALLOC3_VERBOSE(2, nBlockXSize, nBlockYSize));
    1974             : 
    1975          53 :             if (m_eDT == GDT_Int16)
    1976             :             {
    1977          74 :                 ProcessInt16UInt16Tile<GInt16>(
    1978          37 :                     m_pabyCachedTiles,
    1979          37 :                     static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize, true,
    1980          37 :                     CPL_TO_BOOL(bHasNoData), dfNoDataValue, m_usGPKGNull,
    1981             :                     m_dfOffset, m_dfScale, pTempTileBuffer, dfTileOffset,
    1982             :                     dfTileScale, dfTileMin, dfTileMax, dfTileMean, dfTileStdDev,
    1983             :                     nValidPixels);
    1984             :             }
    1985          16 :             else if (m_eDT == GDT_UInt16)
    1986             :             {
    1987          22 :                 ProcessInt16UInt16Tile<GUInt16>(
    1988          11 :                     m_pabyCachedTiles,
    1989          11 :                     static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize, false,
    1990          11 :                     CPL_TO_BOOL(bHasNoData), dfNoDataValue, m_usGPKGNull,
    1991             :                     m_dfOffset, m_dfScale, pTempTileBuffer, dfTileOffset,
    1992             :                     dfTileScale, dfTileMin, dfTileMax, dfTileMean, dfTileStdDev,
    1993             :                     nValidPixels);
    1994             :             }
    1995           5 :             else if (m_eDT == GDT_Float32)
    1996             :             {
    1997           5 :                 const float *pSrc =
    1998             :                     reinterpret_cast<float *>(m_pabyCachedTiles);
    1999           5 :                 float fMin = 0.0f;
    2000           5 :                 float fMax = 0.0f;
    2001           5 :                 double dfM2 = 0.0;
    2002      327685 :                 for (GPtrDiff_t i = 0;
    2003      327685 :                      i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize;
    2004             :                      i++)
    2005             :                 {
    2006      327680 :                     const float fVal = pSrc[i];
    2007      327680 :                     if (bHasNanNoData)
    2008             :                     {
    2009           0 :                         if (std::isnan(fVal))
    2010      196206 :                             continue;
    2011             :                     }
    2012      720494 :                     else if (bHasNoData &&
    2013      327680 :                              fVal == static_cast<float>(dfNoDataValue))
    2014             :                     {
    2015      196206 :                         continue;
    2016             :                     }
    2017      131474 :                     if (std::isinf(fVal))
    2018           0 :                         continue;
    2019             : 
    2020      131474 :                     if (nValidPixels == 0)
    2021             :                     {
    2022           5 :                         fMin = fVal;
    2023           5 :                         fMax = fVal;
    2024             :                     }
    2025             :                     else
    2026             :                     {
    2027      131469 :                         fMin = std::min(fMin, fVal);
    2028      131469 :                         fMax = std::max(fMax, fVal);
    2029             :                     }
    2030      131474 :                     nValidPixels++;
    2031      131474 :                     const double dfDelta = fVal - dfTileMean;
    2032      131474 :                     dfTileMean += dfDelta / nValidPixels;
    2033      131474 :                     dfM2 += dfDelta * (fVal - dfTileMean);
    2034             :                 }
    2035           5 :                 dfTileMin = fMin;
    2036           5 :                 dfTileMax = fMax;
    2037           5 :                 if (nValidPixels)
    2038           5 :                     dfTileStdDev = sqrt(dfM2 / nValidPixels);
    2039             : 
    2040           5 :                 double dfGlobalMin = (fMin - m_dfOffset) / m_dfScale;
    2041           5 :                 double dfGlobalMax = (fMax - m_dfOffset) / m_dfScale;
    2042           5 :                 if (dfGlobalMax > dfGlobalMin)
    2043             :                 {
    2044           4 :                     if (bHasNoData && m_usGPKGNull == 65535 &&
    2045           2 :                         dfGlobalMax - dfGlobalMin >= 65534.0)
    2046             :                     {
    2047           1 :                         dfTileOffset = dfGlobalMin;
    2048           1 :                         dfTileScale = (dfGlobalMax - dfGlobalMin) / 65534.0;
    2049             :                     }
    2050           3 :                     else if (bHasNoData && m_usGPKGNull == 0 &&
    2051           0 :                              (dfNoDataValue - m_dfOffset) / m_dfScale != 0)
    2052             :                     {
    2053           0 :                         dfTileOffset =
    2054           0 :                             (65535.0 * dfGlobalMin - dfGlobalMax) / 65534.0;
    2055           0 :                         dfTileScale = dfGlobalMin - dfTileOffset;
    2056             :                     }
    2057             :                     else
    2058             :                     {
    2059           3 :                         dfTileOffset = dfGlobalMin;
    2060           3 :                         dfTileScale = (dfGlobalMax - dfGlobalMin) / 65535.0;
    2061             :                     }
    2062             :                 }
    2063             :                 else
    2064             :                 {
    2065           1 :                     dfTileOffset = dfGlobalMin;
    2066             :                 }
    2067             : 
    2068      327685 :                 for (GPtrDiff_t i = 0;
    2069      327685 :                      i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize;
    2070             :                      i++)
    2071             :                 {
    2072      327680 :                     const float fVal = pSrc[i];
    2073      327680 :                     if (bHasNanNoData)
    2074             :                     {
    2075           0 :                         if (std::isnan(fVal))
    2076             :                         {
    2077           0 :                             pTempTileBuffer[i] = m_usGPKGNull;
    2078           0 :                             continue;
    2079             :                         }
    2080             :                     }
    2081      327680 :                     else if (bHasNoData)
    2082             :                     {
    2083      196608 :                         if (fVal == static_cast<float>(dfNoDataValue))
    2084             :                         {
    2085      196206 :                             pTempTileBuffer[i] = m_usGPKGNull;
    2086      196206 :                             continue;
    2087             :                         }
    2088             :                     }
    2089             :                     double dfVal =
    2090      131474 :                         std::isfinite(fVal)
    2091      131474 :                             ? ((fVal - m_dfOffset) / m_dfScale - dfTileOffset) /
    2092             :                                   dfTileScale
    2093             :                         : (fVal > 0) ? 65535
    2094      131474 :                                      : 0;
    2095      131474 :                     CPLAssert(dfVal >= 0.0 && dfVal < 65535.5);
    2096      131474 :                     pTempTileBuffer[i] = static_cast<GUInt16>(dfVal + 0.5);
    2097      131474 :                     if (bHasNoData && pTempTileBuffer[i] == m_usGPKGNull)
    2098             :                     {
    2099           1 :                         if (m_usGPKGNull > 0)
    2100           1 :                             pTempTileBuffer[i]--;
    2101             :                         else
    2102           0 :                             pTempTileBuffer[i]++;
    2103             :                     }
    2104             :                 }
    2105             :             }
    2106             : 
    2107          53 :             auto hBand = MEMCreateRasterBandEx(
    2108             :                 poMEMDS, 1, reinterpret_cast<GByte *>(pTempTileBuffer),
    2109             :                 GDT_UInt16, 0, 0, false);
    2110          53 :             poMEMDS->AddMEMBand(hBand);
    2111             :         }
    2112        4111 :         else if (m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
    2113             :         {
    2114           5 :             const float *pSrc = reinterpret_cast<float *>(m_pabyCachedTiles);
    2115           5 :             float fMin = 0.0f;
    2116           5 :             float fMax = 0.0f;
    2117           5 :             double dfM2 = 0.0;
    2118      327685 :             for (GPtrDiff_t i = 0;
    2119      327685 :                  i < static_cast<GPtrDiff_t>(nBlockXSize) * nBlockYSize; i++)
    2120             :             {
    2121      327680 :                 const float fVal = pSrc[i];
    2122      327680 :                 if (bHasNanNoData)
    2123             :                 {
    2124           0 :                     if (std::isnan(fVal))
    2125       65136 :                         continue;
    2126             :                 }
    2127      458352 :                 else if (bHasNoData &&
    2128      327680 :                          fVal == static_cast<float>(dfNoDataValue))
    2129             :                 {
    2130       65136 :                     continue;
    2131             :                 }
    2132             : 
    2133      262544 :                 if (nValidPixels == 0)
    2134             :                 {
    2135           5 :                     fMin = fVal;
    2136           5 :                     fMax = fVal;
    2137             :                 }
    2138             :                 else
    2139             :                 {
    2140      262539 :                     fMin = std::min(fMin, fVal);
    2141      262539 :                     fMax = std::max(fMax, fVal);
    2142             :                 }
    2143      262544 :                 nValidPixels++;
    2144      262544 :                 const double dfDelta = fVal - dfTileMean;
    2145      262544 :                 dfTileMean += dfDelta / nValidPixels;
    2146      262544 :                 dfM2 += dfDelta * (fVal - dfTileMean);
    2147             :             }
    2148           5 :             dfTileMin = fMin;
    2149           5 :             dfTileMax = fMax;
    2150           5 :             if (nValidPixels)
    2151           5 :                 dfTileStdDev = sqrt(dfM2 / nValidPixels);
    2152             : 
    2153           5 :             auto hBand = MEMCreateRasterBandEx(poMEMDS, 1, m_pabyCachedTiles,
    2154             :                                                GDT_Float32, 0, 0, false);
    2155           5 :             poMEMDS->AddMEMBand(hBand);
    2156             :         }
    2157             :         else
    2158             :         {
    2159        4106 :             CPLAssert(m_eDT == GDT_Byte);
    2160       16190 :             for (int i = 0; i < nTileBands; i++)
    2161             :             {
    2162       12084 :                 int iSrc = i;
    2163       12084 :                 if (nBands == 1 && m_poCT == nullptr && nTileBands == 3)
    2164           3 :                     iSrc = 0;
    2165       12081 :                 else if (nBands == 1 && m_poCT == nullptr && bPartialTile &&
    2166             :                          nTileBands == 4)
    2167           8 :                     iSrc = (i < 3) ? 0 : 3;
    2168       12073 :                 else if (nBands == 2 && nTileBands >= 3)
    2169         536 :                     iSrc = (i < 3) ? 0 : 1;
    2170             : 
    2171       24168 :                 auto hBand = MEMCreateRasterBandEx(
    2172             :                     poMEMDS, i + 1,
    2173       12084 :                     m_pabyCachedTiles + iSrc * nBlockXSize * nBlockYSize,
    2174             :                     GDT_Byte, 0, 0, false);
    2175       12084 :                 poMEMDS->AddMEMBand(hBand);
    2176             : 
    2177       12084 :                 if (i == 0 && nTileBands == 1 && m_poCT != nullptr)
    2178          15 :                     poMEMDS->GetRasterBand(1)->SetColorTable(m_poCT);
    2179             :             }
    2180             :         }
    2181             : 
    2182        4164 :         if ((m_eTF == GPKG_TF_PNG_16BIT || m_eTF == GPKG_TF_TIFF_32BIT_FLOAT) &&
    2183          58 :             nValidPixels == 0)
    2184             :         {
    2185             :             // If tile is fully transparent, don't serialize it and remove
    2186             :             // it if it exists.
    2187           1 :             GIntBig nId = GetTileId(nRow, nCol);
    2188           1 :             if (nId > 0)
    2189             :             {
    2190           1 :                 DeleteTile(nRow, nCol);
    2191             : 
    2192           1 :                 DeleteFromGriddedTileAncillary(nId);
    2193             :             }
    2194             : 
    2195           1 :             CPLFree(pTempTileBuffer);
    2196           1 :             delete poMEMDS;
    2197           2 :             return CE_None;
    2198             :         }
    2199             : 
    2200        4163 :         if (m_eTF == GPKG_TF_PNG8 && nTileBands == 1 && nBands >= 3)
    2201             :         {
    2202          10 :             CPLAssert(bAllDirty);
    2203             : 
    2204          10 :             auto poMEM_RGB_DS = MEMDataset::Create("", nBlockXSize, nBlockYSize,
    2205             :                                                    0, GDT_Byte, nullptr);
    2206          40 :             for (int i = 0; i < 3; i++)
    2207             :             {
    2208          60 :                 auto hBand = MEMCreateRasterBandEx(
    2209          30 :                     poMEMDS, i + 1, m_pabyCachedTiles + i * nBandBlockSize,
    2210             :                     GDT_Byte, 0, 0, false);
    2211          30 :                 poMEM_RGB_DS->AddMEMBand(hBand);
    2212             :             }
    2213             : 
    2214          10 :             if (m_pabyHugeColorArray == nullptr)
    2215             :             {
    2216           5 :                 if (nBlockXSize <= 65536 / nBlockYSize)
    2217           4 :                     m_pabyHugeColorArray =
    2218           4 :                         VSIMalloc(MEDIAN_CUT_AND_DITHER_BUFFER_SIZE_65536);
    2219             :                 else
    2220           1 :                     m_pabyHugeColorArray =
    2221           1 :                         VSIMalloc2(256 * 256 * 256, sizeof(GUInt32));
    2222             :             }
    2223             : 
    2224          10 :             GDALColorTable *poCT = new GDALColorTable();
    2225          20 :             GDALComputeMedianCutPCTInternal(
    2226          10 :                 poMEM_RGB_DS->GetRasterBand(1), poMEM_RGB_DS->GetRasterBand(2),
    2227          10 :                 poMEM_RGB_DS->GetRasterBand(3),
    2228             :                 /*NULL, NULL, NULL,*/
    2229          10 :                 m_pabyCachedTiles, m_pabyCachedTiles + nBandBlockSize,
    2230          10 :                 m_pabyCachedTiles + 2 * nBandBlockSize, nullptr,
    2231             :                 256, /* max colors */
    2232             :                 8,   /* bit depth */
    2233             :                 static_cast<GUInt32 *>(
    2234          10 :                     m_pabyHugeColorArray), /* preallocated histogram */
    2235             :                 poCT, nullptr, nullptr);
    2236             : 
    2237          20 :             GDALDitherRGB2PCTInternal(
    2238          10 :                 poMEM_RGB_DS->GetRasterBand(1), poMEM_RGB_DS->GetRasterBand(2),
    2239          10 :                 poMEM_RGB_DS->GetRasterBand(3), poMEMDS->GetRasterBand(1), poCT,
    2240             :                 8, /* bit depth */
    2241             :                 static_cast<GInt16 *>(
    2242          10 :                     m_pabyHugeColorArray), /* pasDynamicColorMap */
    2243          10 :                 m_bDither, nullptr, nullptr);
    2244          10 :             poMEMDS->GetRasterBand(1)->SetColorTable(poCT);
    2245          10 :             delete poCT;
    2246          10 :             GDALClose(poMEM_RGB_DS);
    2247             :         }
    2248        4153 :         else if (nBands == 1 && m_poCT != nullptr && nTileBands > 1)
    2249             :         {
    2250             :             GByte abyCT[4 * 256];
    2251           7 :             const int nEntries = std::min(256, m_poCT->GetColorEntryCount());
    2252        1799 :             for (int i = 0; i < nEntries; i++)
    2253             :             {
    2254        1792 :                 const GDALColorEntry *psEntry = m_poCT->GetColorEntry(i);
    2255        1792 :                 abyCT[4 * i] = static_cast<GByte>(psEntry->c1);
    2256        1792 :                 abyCT[4 * i + 1] = static_cast<GByte>(psEntry->c2);
    2257        1792 :                 abyCT[4 * i + 2] = static_cast<GByte>(psEntry->c3);
    2258        1792 :                 abyCT[4 * i + 3] = static_cast<GByte>(psEntry->c4);
    2259             :             }
    2260           7 :             for (int i = nEntries; i < 256; i++)
    2261             :             {
    2262           0 :                 abyCT[4 * i] = 0;
    2263           0 :                 abyCT[4 * i + 1] = 0;
    2264           0 :                 abyCT[4 * i + 2] = 0;
    2265           0 :                 abyCT[4 * i + 3] = 0;
    2266             :             }
    2267           7 :             if (iYOff > 0)
    2268             :             {
    2269           0 :                 memset(m_pabyCachedTiles + 0 * nBandBlockSize, 0,
    2270           0 :                        nBlockXSize * iYOff);
    2271           0 :                 memset(m_pabyCachedTiles + 1 * nBandBlockSize, 0,
    2272           0 :                        nBlockXSize * iYOff);
    2273           0 :                 memset(m_pabyCachedTiles + 2 * nBandBlockSize, 0,
    2274           0 :                        nBlockXSize * iYOff);
    2275           0 :                 memset(m_pabyCachedTiles + 3 * nBandBlockSize, 0,
    2276           0 :                        nBlockXSize * iYOff);
    2277             :             }
    2278        1057 :             for (GPtrDiff_t iY = iYOff; iY < iYOff + iYCount; iY++)
    2279             :             {
    2280        1050 :                 if (iXOff > 0)
    2281             :                 {
    2282           0 :                     const GPtrDiff_t i = iY * nBlockXSize;
    2283           0 :                     memset(m_pabyCachedTiles + 0 * nBandBlockSize + i, 0,
    2284             :                            iXOff);
    2285           0 :                     memset(m_pabyCachedTiles + 1 * nBandBlockSize + i, 0,
    2286             :                            iXOff);
    2287           0 :                     memset(m_pabyCachedTiles + 2 * nBandBlockSize + i, 0,
    2288             :                            iXOff);
    2289           0 :                     memset(m_pabyCachedTiles + 3 * nBandBlockSize + i, 0,
    2290             :                            iXOff);
    2291             :                 }
    2292      183250 :                 for (int iX = iXOff; iX < iXOff + iXCount; iX++)
    2293             :                 {
    2294      182200 :                     const GPtrDiff_t i = iY * nBlockXSize + iX;
    2295      182200 :                     GByte byVal = m_pabyCachedTiles[i];
    2296      182200 :                     m_pabyCachedTiles[i] = abyCT[4 * byVal];
    2297      182200 :                     m_pabyCachedTiles[i + 1 * nBandBlockSize] =
    2298      182200 :                         abyCT[4 * byVal + 1];
    2299      182200 :                     m_pabyCachedTiles[i + 2 * nBandBlockSize] =
    2300      182200 :                         abyCT[4 * byVal + 2];
    2301      182200 :                     m_pabyCachedTiles[i + 3 * nBandBlockSize] =
    2302      182200 :                         abyCT[4 * byVal + 3];
    2303             :                 }
    2304        1050 :                 if (iXOff + iXCount < nBlockXSize)
    2305             :                 {
    2306         800 :                     const GPtrDiff_t i = iY * nBlockXSize + iXOff + iXCount;
    2307         800 :                     memset(m_pabyCachedTiles + 0 * nBandBlockSize + i, 0,
    2308         800 :                            nBlockXSize - (iXOff + iXCount));
    2309         800 :                     memset(m_pabyCachedTiles + 1 * nBandBlockSize + i, 0,
    2310         800 :                            nBlockXSize - (iXOff + iXCount));
    2311         800 :                     memset(m_pabyCachedTiles + 2 * nBandBlockSize + i, 0,
    2312         800 :                            nBlockXSize - (iXOff + iXCount));
    2313         800 :                     memset(m_pabyCachedTiles + 3 * nBandBlockSize + i, 0,
    2314         800 :                            nBlockXSize - (iXOff + iXCount));
    2315             :                 }
    2316             :             }
    2317           7 :             if (iYOff + iYCount < nBlockYSize)
    2318             :             {
    2319           7 :                 const GPtrDiff_t i = (iYOff + iYCount) * nBlockXSize;
    2320           7 :                 memset(m_pabyCachedTiles + 0 * nBandBlockSize + i, 0,
    2321           7 :                        nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
    2322           7 :                 memset(m_pabyCachedTiles + 1 * nBandBlockSize + i, 0,
    2323           7 :                        nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
    2324           7 :                 memset(m_pabyCachedTiles + 2 * nBandBlockSize + i, 0,
    2325           7 :                        nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
    2326           7 :                 memset(m_pabyCachedTiles + 3 * nBandBlockSize + i, 0,
    2327           7 :                        nBlockXSize * (nBlockYSize - (iYOff + iYCount)));
    2328             :             }
    2329             :         }
    2330             : 
    2331             :         char **papszDriverOptions =
    2332        4163 :             CSLSetNameValue(nullptr, "_INTERNAL_DATASET", "YES");
    2333        4163 :         if (EQUAL(pszDriverName, "JPEG") || EQUAL(pszDriverName, "WEBP"))
    2334             :         {
    2335             :             // If not all bands are dirty, then use lossless WEBP
    2336         172 :             if (!bAllDirty && EQUAL(pszDriverName, "WEBP"))
    2337             :             {
    2338           2 :                 papszDriverOptions =
    2339           2 :                     CSLSetNameValue(papszDriverOptions, "LOSSLESS", "YES");
    2340             :             }
    2341             :             else
    2342             :             {
    2343             :                 papszDriverOptions =
    2344         170 :                     CSLSetNameValue(papszDriverOptions, "QUALITY",
    2345             :                                     CPLSPrintf("%d", m_nQuality));
    2346             :             }
    2347             :         }
    2348        3991 :         else if (EQUAL(pszDriverName, "PNG"))
    2349             :         {
    2350        3986 :             papszDriverOptions = CSLSetNameValue(papszDriverOptions, "ZLEVEL",
    2351             :                                                  CPLSPrintf("%d", m_nZLevel));
    2352             :         }
    2353           5 :         else if (EQUAL(pszDriverName, "GTiff"))
    2354             :         {
    2355             :             papszDriverOptions =
    2356           5 :                 CSLSetNameValue(papszDriverOptions, "COMPRESS", "LZW");
    2357           5 :             if (nBlockXSize * nBlockYSize <= 512 * 512)
    2358             :             {
    2359             :                 // If tile is not too big, create it as single-strip TIFF
    2360             :                 papszDriverOptions =
    2361           5 :                     CSLSetNameValue(papszDriverOptions, "BLOCKYSIZE",
    2362             :                                     CPLSPrintf("%d", nBlockYSize));
    2363             :             }
    2364             :         }
    2365             : #ifdef DEBUG
    2366             :         VSIStatBufL sStat;
    2367        4163 :         CPLAssert(VSIStatL(osMemFileName, &sStat) != 0);
    2368             : #endif
    2369             :         GDALDataset *poOutDS =
    2370        4163 :             l_poDriver->CreateCopy(osMemFileName, poMEMDS, FALSE,
    2371             :                                    papszDriverOptions, nullptr, nullptr);
    2372        4163 :         CSLDestroy(papszDriverOptions);
    2373        4163 :         CPLFree(pTempTileBuffer);
    2374             : 
    2375        4163 :         if (poOutDS)
    2376             :         {
    2377        4163 :             GDALClose(poOutDS);
    2378        4163 :             vsi_l_offset nBlobSize = 0;
    2379             :             GByte *pabyBlob =
    2380        4163 :                 VSIGetMemFileBuffer(osMemFileName, &nBlobSize, TRUE);
    2381             : 
    2382             :             /* Create or commit and recreate transaction */
    2383        4163 :             GDALGPKGMBTilesLikePseudoDataset *poMainDS =
    2384        4163 :                 m_poParentDS ? m_poParentDS : this;
    2385        4163 :             if (poMainDS->m_nTileInsertionCount == 0)
    2386             :             {
    2387         189 :                 poMainDS->IStartTransaction();
    2388             :             }
    2389        3974 :             else if (poMainDS->m_nTileInsertionCount == 1000)
    2390             :             {
    2391           3 :                 if (poMainDS->ICommitTransaction() != OGRERR_NONE)
    2392             :                 {
    2393           1 :                     poMainDS->m_nTileInsertionCount = -1;
    2394           1 :                     CPLFree(pabyBlob);
    2395           1 :                     VSIUnlink(osMemFileName);
    2396           1 :                     delete poMEMDS;
    2397           1 :                     return CE_Failure;
    2398             :                 }
    2399           2 :                 poMainDS->IStartTransaction();
    2400           2 :                 poMainDS->m_nTileInsertionCount = 0;
    2401             :             }
    2402        4162 :             poMainDS->m_nTileInsertionCount++;
    2403             : 
    2404             :             char *pszSQL =
    2405        4162 :                 sqlite3_mprintf("INSERT OR REPLACE INTO \"%w\" "
    2406             :                                 "(zoom_level, tile_row, tile_column, "
    2407             :                                 "tile_data) VALUES (%d, %d, %d, ?)",
    2408             :                                 m_osRasterTable.c_str(), m_nZoomLevel,
    2409        4162 :                                 GetRowFromIntoTopConvention(nRow), nCol);
    2410             : #ifdef DEBUG_VERBOSE
    2411             :             CPLDebug("GPKG", "%s", pszSQL);
    2412             : #endif
    2413        4162 :             sqlite3_stmt *hStmt = nullptr;
    2414        4162 :             int rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt, nullptr);
    2415        4162 :             if (rc != SQLITE_OK)
    2416             :             {
    2417           0 :                 CPLFree(pabyBlob);
    2418             :             }
    2419             :             else
    2420             :             {
    2421        4162 :                 sqlite3_bind_blob(hStmt, 1, pabyBlob,
    2422             :                                   static_cast<int>(nBlobSize), CPLFree);
    2423        4162 :                 rc = sqlite3_step(hStmt);
    2424        4162 :                 if (rc == SQLITE_DONE)
    2425        4162 :                     eErr = CE_None;
    2426             :                 else
    2427             :                 {
    2428           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    2429             :                              "Failure when inserting tile (row=%d,col=%d) at "
    2430             :                              "zoom_level=%d : %s",
    2431           0 :                              GetRowFromIntoTopConvention(nRow), nCol,
    2432           0 :                              m_nZoomLevel, sqlite3_errmsg(IGetDB()));
    2433             :                 }
    2434             :             }
    2435        4162 :             sqlite3_finalize(hStmt);
    2436        4162 :             sqlite3_free(pszSQL);
    2437             : 
    2438        4162 :             if (m_eTF == GPKG_TF_PNG_16BIT || m_eTF == GPKG_TF_TIFF_32BIT_FLOAT)
    2439             :             {
    2440          57 :                 GIntBig nTileId = GetTileId(nRow, nCol);
    2441          57 :                 if (nTileId == 0)
    2442           0 :                     eErr = CE_Failure;
    2443             :                 else
    2444             :                 {
    2445          57 :                     DeleteFromGriddedTileAncillary(nTileId);
    2446             : 
    2447          57 :                     pszSQL = sqlite3_mprintf(
    2448             :                         "INSERT INTO gpkg_2d_gridded_tile_ancillary "
    2449             :                         "(tpudt_name, tpudt_id, scale, offset, min, max, "
    2450             :                         "mean, std_dev) VALUES "
    2451             :                         "('%q', ?, %.17g, %.17g, ?, ?, ?, ?)",
    2452             :                         m_osRasterTable.c_str(), dfTileScale, dfTileOffset);
    2453             : #ifdef DEBUG_VERBOSE
    2454             :                     CPLDebug("GPKG", "%s", pszSQL);
    2455             : #endif
    2456          57 :                     hStmt = nullptr;
    2457          57 :                     rc = SQLPrepareWithError(IGetDB(), pszSQL, -1, &hStmt,
    2458             :                                              nullptr);
    2459          57 :                     if (rc != SQLITE_OK)
    2460             :                     {
    2461           0 :                         eErr = CE_Failure;
    2462             :                     }
    2463             :                     else
    2464             :                     {
    2465          57 :                         sqlite3_bind_int64(hStmt, 1, nTileId);
    2466          57 :                         sqlite3_bind_double(hStmt, 2, dfTileMin);
    2467          57 :                         sqlite3_bind_double(hStmt, 3, dfTileMax);
    2468          57 :                         sqlite3_bind_double(hStmt, 4, dfTileMean);
    2469          57 :                         sqlite3_bind_double(hStmt, 5, dfTileStdDev);
    2470          57 :                         rc = sqlite3_step(hStmt);
    2471          57 :                         if (rc == SQLITE_DONE)
    2472             :                         {
    2473          57 :                             eErr = CE_None;
    2474             :                         }
    2475             :                         else
    2476             :                         {
    2477           0 :                             CPLError(CE_Failure, CPLE_AppDefined,
    2478             :                                      "Cannot insert into "
    2479             :                                      "gpkg_2d_gridded_tile_ancillary");
    2480           0 :                             eErr = CE_Failure;
    2481             :                         }
    2482             :                     }
    2483          57 :                     sqlite3_finalize(hStmt);
    2484          57 :                     sqlite3_free(pszSQL);
    2485             :                 }
    2486             :             }
    2487             :         }
    2488             : 
    2489        4162 :         VSIUnlink(osMemFileName);
    2490        4162 :         delete poMEMDS;
    2491             :     }
    2492             :     else
    2493             :     {
    2494           1 :         CPLError(CE_Failure, CPLE_NotSupported, "Cannot find driver %s",
    2495             :                  pszDriverName);
    2496             :     }
    2497             : 
    2498        4163 :     return eErr;
    2499             : }
    2500             : 
    2501             : /************************************************************************/
    2502             : /*                     FlushRemainingShiftedTiles()                     */
    2503             : /************************************************************************/
    2504             : 
    2505             : CPLErr
    2506         609 : GDALGPKGMBTilesLikePseudoDataset::FlushRemainingShiftedTiles(bool bPartialFlush)
    2507             : {
    2508         609 :     if (m_hTempDB == nullptr)
    2509         460 :         return CE_None;
    2510             : 
    2511         745 :     for (int i = 0; i <= 3; i++)
    2512             :     {
    2513         596 :         m_asCachedTilesDesc[i].nRow = -1;
    2514         596 :         m_asCachedTilesDesc[i].nCol = -1;
    2515         596 :         m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
    2516             :     }
    2517             : 
    2518             :     int nBlockXSize, nBlockYSize;
    2519         149 :     IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
    2520         149 :     const int nBands = IGetRasterCount();
    2521         149 :     const int nRasterXSize = IGetRasterBand(1)->GetXSize();
    2522         149 :     const int nRasterYSize = IGetRasterBand(1)->GetYSize();
    2523         149 :     const int nXBlocks = DIV_ROUND_UP(nRasterXSize, nBlockXSize);
    2524         149 :     const int nYBlocks = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
    2525             : 
    2526         149 :     int nPartialActiveTiles = 0;
    2527         149 :     if (bPartialFlush)
    2528             :     {
    2529           2 :         sqlite3_stmt *hStmt = nullptr;
    2530           4 :         CPLString osSQL;
    2531             :         osSQL.Printf("SELECT COUNT(*) FROM partial_tiles WHERE zoom_level = %d "
    2532             :                      "AND partial_flag != 0",
    2533           2 :                      m_nZoomLevel);
    2534           2 :         if (SQLPrepareWithError(m_hTempDB, osSQL.c_str(), -1, &hStmt,
    2535           2 :                                 nullptr) == SQLITE_OK)
    2536             :         {
    2537           2 :             if (sqlite3_step(hStmt) == SQLITE_ROW)
    2538             :             {
    2539           2 :                 nPartialActiveTiles = sqlite3_column_int(hStmt, 0);
    2540           2 :                 CPLDebug("GPKG", "Active partial tiles before flush: %d",
    2541             :                          nPartialActiveTiles);
    2542             :             }
    2543           2 :             sqlite3_finalize(hStmt);
    2544             :         }
    2545             :     }
    2546             : 
    2547         298 :     CPLString osSQL = "SELECT tile_row, tile_column, partial_flag";
    2548         616 :     for (int nBand = 1; nBand <= nBands; nBand++)
    2549             :     {
    2550         467 :         osSQL += CPLSPrintf(", tile_data_band_%d", nBand);
    2551             :     }
    2552             :     osSQL += CPLSPrintf(" FROM partial_tiles WHERE "
    2553             :                         "zoom_level = %d AND partial_flag != 0",
    2554         149 :                         m_nZoomLevel);
    2555         149 :     if (bPartialFlush)
    2556             :     {
    2557           2 :         osSQL += " ORDER BY age";
    2558             :     }
    2559         149 :     const char *pszSQL = osSQL.c_str();
    2560             : 
    2561             : #ifdef DEBUG_VERBOSE
    2562             :     CPLDebug("GPKG", "%s", pszSQL);
    2563             : #endif
    2564         149 :     sqlite3_stmt *hStmt = nullptr;
    2565         149 :     int rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
    2566         149 :     if (rc != SQLITE_OK)
    2567             :     {
    2568           0 :         return CE_Failure;
    2569             :     }
    2570             : 
    2571         149 :     CPLErr eErr = CE_None;
    2572         149 :     bool bGotPartialTiles = false;
    2573         149 :     int nCountFlushedTiles = 0;
    2574         149 :     const size_t nBandBlockSize =
    2575         149 :         static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
    2576         175 :     do
    2577             :     {
    2578         324 :         rc = sqlite3_step(hStmt);
    2579         324 :         if (rc == SQLITE_ROW)
    2580             :         {
    2581         177 :             bGotPartialTiles = true;
    2582             : 
    2583         177 :             int nRow = sqlite3_column_int(hStmt, 0);
    2584         177 :             int nCol = sqlite3_column_int(hStmt, 1);
    2585         177 :             int nPartialFlags = sqlite3_column_int(hStmt, 2);
    2586             : 
    2587         177 :             if (bPartialFlush)
    2588             :             {
    2589             :                 // This method assumes that there are no dirty blocks alive
    2590             :                 // so check this assumption.
    2591             :                 // When called with bPartialFlush = false, FlushCache() has
    2592             :                 // already been called, so no need to check.
    2593          16 :                 bool bFoundDirtyBlock = false;
    2594          16 :                 int nBlockXOff = nCol - m_nShiftXTiles;
    2595          16 :                 int nBlockYOff = nRow - m_nShiftYTiles;
    2596          96 :                 for (int iX = 0; !bFoundDirtyBlock &&
    2597          48 :                                  iX < ((m_nShiftXPixelsMod != 0) ? 2 : 1);
    2598             :                      iX++)
    2599             :                 {
    2600          32 :                     if (nBlockXOff + iX < 0 || nBlockXOff + iX >= nXBlocks)
    2601          12 :                         continue;
    2602         120 :                     for (int iY = 0; !bFoundDirtyBlock &&
    2603          60 :                                      iY < ((m_nShiftYPixelsMod != 0) ? 2 : 1);
    2604             :                          iY++)
    2605             :                     {
    2606          40 :                         if (nBlockYOff + iY < 0 || nBlockYOff + iY >= nYBlocks)
    2607          17 :                             continue;
    2608          69 :                         for (int iBand = 1;
    2609          69 :                              !bFoundDirtyBlock && iBand <= nBands; iBand++)
    2610             :                         {
    2611             :                             GDALRasterBlock *poBlock =
    2612             :                                 cpl::down_cast<GDALGPKGMBTilesLikeRasterBand *>(
    2613          46 :                                     IGetRasterBand(iBand))
    2614          46 :                                     ->AccessibleTryGetLockedBlockRef(
    2615             :                                         nBlockXOff + iX, nBlockYOff + iY);
    2616          46 :                             if (poBlock)
    2617             :                             {
    2618           0 :                                 if (poBlock->GetDirty())
    2619           0 :                                     bFoundDirtyBlock = true;
    2620           0 :                                 poBlock->DropLock();
    2621             :                             }
    2622             :                         }
    2623             :                     }
    2624             :                 }
    2625          16 :                 if (bFoundDirtyBlock)
    2626             :                 {
    2627             : #ifdef DEBUG_VERBOSE
    2628             :                     CPLDebug("GPKG",
    2629             :                              "Skipped flushing tile row = %d, column = %d "
    2630             :                              "because it has dirty block(s) in GDAL cache",
    2631             :                              nRow, nCol);
    2632             : #endif
    2633           0 :                     continue;
    2634             :                 }
    2635             :             }
    2636             : 
    2637         177 :             nCountFlushedTiles++;
    2638         177 :             if (bPartialFlush && nCountFlushedTiles >= nPartialActiveTiles / 2)
    2639             :             {
    2640           2 :                 CPLDebug("GPKG", "Flushed %d tiles", nCountFlushedTiles);
    2641           2 :                 break;
    2642             :             }
    2643             : 
    2644         705 :             for (int nBand = 1; nBand <= nBands; nBand++)
    2645             :             {
    2646         530 :                 if (nPartialFlags & (((1 << 4) - 1) << (4 * (nBand - 1))))
    2647             :                 {
    2648         498 :                     CPLAssert(sqlite3_column_bytes(hStmt, 2 + nBand) ==
    2649             :                               static_cast<int>(nBandBlockSize));
    2650         498 :                     memcpy(m_pabyCachedTiles + (nBand - 1) * nBandBlockSize,
    2651             :                            sqlite3_column_blob(hStmt, 2 + nBand),
    2652             :                            nBandBlockSize);
    2653             :                 }
    2654             :                 else
    2655             :                 {
    2656          32 :                     FillEmptyTileSingleBand(m_pabyCachedTiles +
    2657          32 :                                             (nBand - 1) * nBandBlockSize);
    2658             :                 }
    2659             :             }
    2660             : 
    2661         175 :             int nFullFlags = (1 << (4 * nBands)) - 1;
    2662             : 
    2663             :             // In case the partial flags indicate that there's some quadrant
    2664             :             // missing, check in the main database if there is already a tile
    2665             :             // In which case, use the parts of that tile that aren't in the
    2666             :             // temporary database
    2667         175 :             if (nPartialFlags != nFullFlags)
    2668             :             {
    2669         525 :                 char *pszNewSQL = sqlite3_mprintf(
    2670             :                     "SELECT tile_data%s FROM \"%w\" "
    2671             :                     "WHERE zoom_level = %d AND tile_row = %d AND tile_column = "
    2672             :                     "%d%s",
    2673         175 :                     m_eDT != GDT_Byte ? ", id"
    2674             :                                       : "",  // MBTiles do not have an id
    2675             :                     m_osRasterTable.c_str(), m_nZoomLevel,
    2676         175 :                     GetRowFromIntoTopConvention(nRow), nCol,
    2677         175 :                     !m_osWHERE.empty()
    2678           0 :                         ? CPLSPrintf(" AND (%s)", m_osWHERE.c_str())
    2679             :                         : "");
    2680             : #ifdef DEBUG_VERBOSE
    2681             :                 CPLDebug("GPKG", "%s", pszNewSQL);
    2682             : #endif
    2683         175 :                 sqlite3_stmt *hNewStmt = nullptr;
    2684         175 :                 rc = SQLPrepareWithError(IGetDB(), pszNewSQL, -1, &hNewStmt,
    2685             :                                          nullptr);
    2686         175 :                 if (rc == SQLITE_OK)
    2687             :                 {
    2688         175 :                     rc = sqlite3_step(hNewStmt);
    2689         192 :                     if (rc == SQLITE_ROW &&
    2690          17 :                         sqlite3_column_type(hNewStmt, 0) == SQLITE_BLOB)
    2691             :                     {
    2692          17 :                         const int nBytes = sqlite3_column_bytes(hNewStmt, 0);
    2693             :                         GIntBig nTileId =
    2694          17 :                             (m_eDT == GDT_Byte)
    2695          17 :                                 ? 0
    2696           0 :                                 : sqlite3_column_int64(hNewStmt, 1);
    2697             :                         GByte *pabyRawData =
    2698             :                             const_cast<GByte *>(static_cast<const GByte *>(
    2699          17 :                                 sqlite3_column_blob(hNewStmt, 0)));
    2700             :                         const CPLString osMemFileName(
    2701          34 :                             VSIMemGenerateHiddenFilename("gpkg_read_tile"));
    2702          17 :                         VSILFILE *fp = VSIFileFromMemBuffer(
    2703             :                             osMemFileName.c_str(), pabyRawData, nBytes, FALSE);
    2704          17 :                         VSIFCloseL(fp);
    2705             : 
    2706          17 :                         double dfTileOffset = 0.0;
    2707          17 :                         double dfTileScale = 1.0;
    2708          17 :                         GetTileOffsetAndScale(nTileId, dfTileOffset,
    2709             :                                               dfTileScale);
    2710          17 :                         const int nTileBands = m_eDT == GDT_Byte ? 4 : 1;
    2711          17 :                         GByte *pabyTemp =
    2712          17 :                             m_pabyCachedTiles + nTileBands * nBandBlockSize;
    2713          17 :                         ReadTile(osMemFileName, pabyTemp, dfTileOffset,
    2714             :                                  dfTileScale);
    2715          17 :                         VSIUnlink(osMemFileName);
    2716             : 
    2717          17 :                         int iYQuadrantMax = (m_nShiftYPixelsMod) ? 1 : 0;
    2718          17 :                         int iXQuadrantMax = (m_nShiftXPixelsMod) ? 1 : 0;
    2719          51 :                         for (int iYQuadrant = 0; iYQuadrant <= iYQuadrantMax;
    2720             :                              iYQuadrant++)
    2721             :                         {
    2722         100 :                             for (int iXQuadrant = 0;
    2723         100 :                                  iXQuadrant <= iXQuadrantMax; iXQuadrant++)
    2724             :                             {
    2725         226 :                                 for (int nBand = 1; nBand <= nBands; nBand++)
    2726             :                                 {
    2727         160 :                                     int iQuadrantFlag = 0;
    2728         160 :                                     if (iXQuadrant == 0 && iYQuadrant == 0)
    2729          42 :                                         iQuadrantFlag |= (1 << 0);
    2730         160 :                                     if (iXQuadrant == iXQuadrantMax &&
    2731             :                                         iYQuadrant == 0)
    2732          42 :                                         iQuadrantFlag |= (1 << 1);
    2733         160 :                                     if (iXQuadrant == 0 &&
    2734             :                                         iYQuadrant == iYQuadrantMax)
    2735          42 :                                         iQuadrantFlag |= (1 << 2);
    2736         160 :                                     if (iXQuadrant == iXQuadrantMax &&
    2737             :                                         iYQuadrant == iYQuadrantMax)
    2738          42 :                                         iQuadrantFlag |= (1 << 3);
    2739         160 :                                     int nLocalFlag = iQuadrantFlag
    2740         160 :                                                      << (4 * (nBand - 1));
    2741         160 :                                     if (!(nPartialFlags & nLocalFlag))
    2742             :                                     {
    2743             :                                         int nXOff, nYOff, nXSize, nYSize;
    2744         118 :                                         if (iXQuadrant == 0 &&
    2745          60 :                                             m_nShiftXPixelsMod != 0)
    2746             :                                         {
    2747          56 :                                             nXOff = 0;
    2748          56 :                                             nXSize = m_nShiftXPixelsMod;
    2749             :                                         }
    2750             :                                         else
    2751             :                                         {
    2752          62 :                                             nXOff = m_nShiftXPixelsMod;
    2753          62 :                                             nXSize = nBlockXSize -
    2754          62 :                                                      m_nShiftXPixelsMod;
    2755             :                                         }
    2756         118 :                                         if (iYQuadrant == 0 &&
    2757          63 :                                             m_nShiftYPixelsMod != 0)
    2758             :                                         {
    2759          63 :                                             nYOff = 0;
    2760          63 :                                             nYSize = m_nShiftYPixelsMod;
    2761             :                                         }
    2762             :                                         else
    2763             :                                         {
    2764          55 :                                             nYOff = m_nShiftYPixelsMod;
    2765          55 :                                             nYSize = nBlockYSize -
    2766          55 :                                                      m_nShiftYPixelsMod;
    2767             :                                         }
    2768         118 :                                         for (int iY = nYOff;
    2769       13888 :                                              iY < nYOff + nYSize; iY++)
    2770             :                                         {
    2771       13770 :                                             memcpy(m_pabyCachedTiles +
    2772       13770 :                                                        ((static_cast<size_t>(
    2773       13770 :                                                              nBand - 1) *
    2774       13770 :                                                              nBlockYSize +
    2775       13770 :                                                          iY) *
    2776       13770 :                                                             nBlockXSize +
    2777       13770 :                                                         nXOff) *
    2778       13770 :                                                            m_nDTSize,
    2779       13770 :                                                    pabyTemp +
    2780       13770 :                                                        ((static_cast<size_t>(
    2781       13770 :                                                              nBand - 1) *
    2782       13770 :                                                              nBlockYSize +
    2783       13770 :                                                          iY) *
    2784       13770 :                                                             nBlockXSize +
    2785       13770 :                                                         nXOff) *
    2786       13770 :                                                            m_nDTSize,
    2787       13770 :                                                    static_cast<size_t>(nXSize) *
    2788       13770 :                                                        m_nDTSize);
    2789             :                                         }
    2790             :                                     }
    2791             :                                 }
    2792             :                             }
    2793             :                         }
    2794             :                     }
    2795         158 :                     else if (rc != SQLITE_DONE)
    2796             :                     {
    2797           0 :                         CPLError(CE_Failure, CPLE_AppDefined,
    2798             :                                  "sqlite3_step(%s) failed: %s", pszNewSQL,
    2799             :                                  sqlite3_errmsg(m_hTempDB));
    2800             :                     }
    2801         175 :                     sqlite3_finalize(hNewStmt);
    2802             :                 }
    2803         175 :                 sqlite3_free(pszNewSQL);
    2804             :             }
    2805             : 
    2806         175 :             m_asCachedTilesDesc[0].nRow = nRow;
    2807         175 :             m_asCachedTilesDesc[0].nCol = nCol;
    2808         175 :             m_asCachedTilesDesc[0].nIdxWithinTileData = 0;
    2809         175 :             m_asCachedTilesDesc[0].abBandDirty[0] = true;
    2810         175 :             m_asCachedTilesDesc[0].abBandDirty[1] = true;
    2811         175 :             m_asCachedTilesDesc[0].abBandDirty[2] = true;
    2812         175 :             m_asCachedTilesDesc[0].abBandDirty[3] = true;
    2813             : 
    2814         175 :             eErr = WriteTile();
    2815             : 
    2816         175 :             if (eErr == CE_None && bPartialFlush)
    2817             :             {
    2818             :                 pszSQL =
    2819          14 :                     CPLSPrintf("DELETE FROM partial_tiles WHERE zoom_level = "
    2820             :                                "%d AND tile_row = %d AND tile_column = %d",
    2821             :                                m_nZoomLevel, nRow, nCol);
    2822             : #ifdef DEBUG_VERBOSE
    2823             :                 CPLDebug("GPKG", "%s", pszSQL);
    2824             : #endif
    2825          14 :                 if (SQLCommand(m_hTempDB, pszSQL) != OGRERR_NONE)
    2826           0 :                     eErr = CE_None;
    2827             :             }
    2828             :         }
    2829             :         else
    2830             :         {
    2831         147 :             if (rc != SQLITE_DONE)
    2832             :             {
    2833           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2834             :                          "sqlite3_step(%s) failed: %s", pszSQL,
    2835             :                          sqlite3_errmsg(m_hTempDB));
    2836             :             }
    2837         147 :             break;
    2838             :         }
    2839         175 :     } while (eErr == CE_None);
    2840             : 
    2841         149 :     sqlite3_finalize(hStmt);
    2842             : 
    2843         149 :     if (bPartialFlush && nCountFlushedTiles < nPartialActiveTiles / 2)
    2844             :     {
    2845           0 :         CPLDebug("GPKG", "Flushed %d tiles. Target was %d", nCountFlushedTiles,
    2846             :                  nPartialActiveTiles / 2);
    2847             :     }
    2848             : 
    2849         149 :     if (bGotPartialTiles && !bPartialFlush)
    2850             :     {
    2851             : #ifdef DEBUG_VERBOSE
    2852             :         pszSQL = CPLSPrintf("SELECT p1.id, p1.tile_row, p1.tile_column FROM "
    2853             :                             "partial_tiles p1, partial_tiles p2 "
    2854             :                             "WHERE p1.zoom_level = %d AND p2.zoom_level = %d "
    2855             :                             "AND p1.tile_row = p2.tile_row AND p1.tile_column "
    2856             :                             "= p2.tile_column AND p2.partial_flag != 0",
    2857             :                             -1 - m_nZoomLevel, m_nZoomLevel);
    2858             :         rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
    2859             :         CPLAssert(rc == SQLITE_OK);
    2860             :         while ((rc = sqlite3_step(hStmt)) == SQLITE_ROW)
    2861             :         {
    2862             :             CPLDebug("GPKG",
    2863             :                      "Conflict: existing id = %d, tile_row = %d, tile_column = "
    2864             :                      "%d, zoom_level = %d",
    2865             :                      sqlite3_column_int(hStmt, 0), sqlite3_column_int(hStmt, 1),
    2866             :                      sqlite3_column_int(hStmt, 2), m_nZoomLevel);
    2867             :         }
    2868             :         sqlite3_finalize(hStmt);
    2869             : #endif
    2870             : 
    2871         110 :         pszSQL = CPLSPrintf("UPDATE partial_tiles SET zoom_level = %d, "
    2872             :                             "partial_flag = 0, age = -1 WHERE zoom_level = %d "
    2873             :                             "AND partial_flag != 0",
    2874          55 :                             -1 - m_nZoomLevel, m_nZoomLevel);
    2875             : #ifdef DEBUG_VERBOSE
    2876             :         CPLDebug("GPKG", "%s", pszSQL);
    2877             : #endif
    2878          55 :         SQLCommand(m_hTempDB, pszSQL);
    2879             :     }
    2880             : 
    2881         149 :     return eErr;
    2882             : }
    2883             : 
    2884             : /************************************************************************/
    2885             : /*                DoPartialFlushOfPartialTilesIfNecessary()             */
    2886             : /************************************************************************/
    2887             : 
    2888             : CPLErr
    2889        4556 : GDALGPKGMBTilesLikePseudoDataset::DoPartialFlushOfPartialTilesIfNecessary()
    2890             : {
    2891        4556 :     time_t nCurTimeStamp = time(nullptr);
    2892        4556 :     if (m_nLastSpaceCheckTimestamp == 0)
    2893           3 :         m_nLastSpaceCheckTimestamp = nCurTimeStamp;
    2894        4556 :     if (m_nLastSpaceCheckTimestamp > 0 &&
    2895          81 :         (m_bForceTempDBCompaction ||
    2896          13 :          nCurTimeStamp - m_nLastSpaceCheckTimestamp > 10))
    2897             :     {
    2898          68 :         m_nLastSpaceCheckTimestamp = nCurTimeStamp;
    2899             :         GIntBig nFreeSpace =
    2900          68 :             VSIGetDiskFreeSpace(CPLGetDirnameSafe(m_osTempDBFilename).c_str());
    2901          68 :         bool bTryFreeing = false;
    2902          68 :         if (nFreeSpace >= 0 && nFreeSpace < 1024 * 1024 * 1024)
    2903             :         {
    2904           0 :             CPLDebug("GPKG",
    2905             :                      "Free space below 1GB. Flushing part of partial tiles");
    2906           0 :             bTryFreeing = true;
    2907             :         }
    2908             :         else
    2909             :         {
    2910             :             VSIStatBufL sStat;
    2911          68 :             if (VSIStatL(m_osTempDBFilename, &sStat) == 0)
    2912             :             {
    2913          68 :                 GIntBig nTempSpace = sStat.st_size;
    2914         136 :                 if (VSIStatL((m_osTempDBFilename + "-journal").c_str(),
    2915          68 :                              &sStat) == 0)
    2916           0 :                     nTempSpace += sStat.st_size;
    2917         136 :                 else if (VSIStatL((m_osTempDBFilename + "-wal").c_str(),
    2918          68 :                                   &sStat) == 0)
    2919           0 :                     nTempSpace += sStat.st_size;
    2920             : 
    2921             :                 int nBlockXSize, nBlockYSize;
    2922          68 :                 IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
    2923          68 :                 const int nBands = IGetRasterCount();
    2924             : 
    2925          68 :                 if (nTempSpace >
    2926          68 :                     4 * static_cast<GIntBig>(IGetRasterBand(1)->GetXSize()) *
    2927          68 :                         nBlockYSize * nBands * m_nDTSize)
    2928             :                 {
    2929           2 :                     CPLDebug("GPKG",
    2930             :                              "Partial tiles DB is " CPL_FRMT_GIB
    2931             :                              " bytes. Flushing part of partial tiles",
    2932             :                              nTempSpace);
    2933           2 :                     bTryFreeing = true;
    2934             :                 }
    2935             :             }
    2936             :         }
    2937          68 :         if (bTryFreeing)
    2938             :         {
    2939           2 :             if (FlushRemainingShiftedTiles(true /* partial flush*/) != CE_None)
    2940             :             {
    2941           0 :                 return CE_Failure;
    2942             :             }
    2943           2 :             SQLCommand(m_hTempDB,
    2944             :                        "DELETE FROM partial_tiles WHERE zoom_level < 0");
    2945           2 :             SQLCommand(m_hTempDB, "VACUUM");
    2946             :         }
    2947             :     }
    2948        4556 :     return CE_None;
    2949             : }
    2950             : 
    2951             : /************************************************************************/
    2952             : /*                         WriteShiftedTile()                           */
    2953             : /************************************************************************/
    2954             : 
    2955        4556 : CPLErr GDALGPKGMBTilesLikePseudoDataset::WriteShiftedTile(
    2956             :     int nRow, int nCol, int nBand, int nDstXOffset, int nDstYOffset,
    2957             :     int nDstXSize, int nDstYSize)
    2958             : {
    2959        4556 :     CPLAssert(m_nShiftXPixelsMod || m_nShiftYPixelsMod);
    2960        4556 :     CPLAssert(nRow >= 0);
    2961        4556 :     CPLAssert(nCol >= 0);
    2962        4556 :     CPLAssert(nRow < m_nTileMatrixHeight);
    2963        4556 :     CPLAssert(nCol < m_nTileMatrixWidth);
    2964             : 
    2965        4556 :     if (m_hTempDB == nullptr &&
    2966          55 :         (m_poParentDS == nullptr || m_poParentDS->m_hTempDB == nullptr))
    2967             :     {
    2968             :         const char *pszBaseFilename =
    2969          48 :             m_poParentDS ? m_poParentDS->IGetFilename() : IGetFilename();
    2970             :         m_osTempDBFilename =
    2971          48 :             CPLResetExtensionSafe(pszBaseFilename, "partial_tiles.db");
    2972          48 :         CPLPushErrorHandler(CPLQuietErrorHandler);
    2973          48 :         VSIUnlink(m_osTempDBFilename);
    2974          48 :         CPLPopErrorHandler();
    2975          48 :         m_hTempDB = nullptr;
    2976          48 :         int rc = 0;
    2977          48 :         if (STARTS_WITH(m_osTempDBFilename, "/vsi"))
    2978             :         {
    2979          33 :             m_pMyVFS = OGRSQLiteCreateVFS(nullptr, nullptr);
    2980          33 :             sqlite3_vfs_register(m_pMyVFS, 0);
    2981          33 :             rc = sqlite3_open_v2(m_osTempDBFilename, &m_hTempDB,
    2982             :                                  SQLITE_OPEN_READWRITE | SQLITE_OPEN_CREATE |
    2983             :                                      SQLITE_OPEN_NOMUTEX,
    2984          33 :                                  m_pMyVFS->zName);
    2985             :         }
    2986             :         else
    2987             :         {
    2988          15 :             rc = sqlite3_open(m_osTempDBFilename, &m_hTempDB);
    2989             :         }
    2990          48 :         if (rc != SQLITE_OK || m_hTempDB == nullptr)
    2991             :         {
    2992           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    2993             :                      "Cannot create temporary database %s",
    2994             :                      m_osTempDBFilename.c_str());
    2995           0 :             return CE_Failure;
    2996             :         }
    2997          48 :         SQLCommand(m_hTempDB, "PRAGMA synchronous = OFF");
    2998          48 :         SQLCommand(m_hTempDB, "PRAGMA journal_mode = OFF");
    2999          48 :         SQLCommand(m_hTempDB, "CREATE TABLE partial_tiles("
    3000             :                               "id INTEGER PRIMARY KEY AUTOINCREMENT,"
    3001             :                               "zoom_level INTEGER NOT NULL,"
    3002             :                               "tile_column INTEGER NOT NULL,"
    3003             :                               "tile_row INTEGER NOT NULL,"
    3004             :                               "tile_data_band_1 BLOB,"
    3005             :                               "tile_data_band_2 BLOB,"
    3006             :                               "tile_data_band_3 BLOB,"
    3007             :                               "tile_data_band_4 BLOB,"
    3008             :                               "partial_flag INTEGER NOT NULL,"
    3009             :                               "age INTEGER NOT NULL,"
    3010             :                               "UNIQUE (zoom_level, tile_column, tile_row))");
    3011          48 :         SQLCommand(m_hTempDB, "CREATE INDEX partial_tiles_partial_flag_idx "
    3012             :                               "ON partial_tiles(partial_flag)");
    3013          48 :         SQLCommand(m_hTempDB, "CREATE INDEX partial_tiles_age_idx "
    3014             :                               "ON partial_tiles(age)");
    3015             : 
    3016          48 :         if (m_poParentDS != nullptr)
    3017             :         {
    3018           4 :             m_poParentDS->m_osTempDBFilename = m_osTempDBFilename;
    3019           4 :             m_poParentDS->m_hTempDB = m_hTempDB;
    3020             :         }
    3021             :     }
    3022             : 
    3023        4556 :     if (m_poParentDS != nullptr)
    3024        4020 :         m_hTempDB = m_poParentDS->m_hTempDB;
    3025             : 
    3026             :     int nBlockXSize, nBlockYSize;
    3027        4556 :     IGetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
    3028        4556 :     const int nBands = IGetRasterCount();
    3029        4556 :     const size_t nBandBlockSize =
    3030        4556 :         static_cast<size_t>(nBlockXSize) * nBlockYSize * m_nDTSize;
    3031             : 
    3032        4556 :     int iQuadrantFlag = 0;
    3033        4556 :     if (nDstXOffset == 0 && nDstYOffset == 0)
    3034        1031 :         iQuadrantFlag |= (1 << 0);
    3035        4556 :     if ((nDstXOffset != 0 || nDstXOffset + nDstXSize == nBlockXSize) &&
    3036             :         nDstYOffset == 0)
    3037        1125 :         iQuadrantFlag |= (1 << 1);
    3038        4556 :     if (nDstXOffset == 0 &&
    3039        1031 :         (nDstYOffset != 0 || nDstYOffset + nDstYSize == nBlockYSize))
    3040        1135 :         iQuadrantFlag |= (1 << 2);
    3041        4556 :     if ((nDstXOffset != 0 || nDstXOffset + nDstXSize == nBlockXSize) &&
    3042        1125 :         (nDstYOffset != 0 || nDstYOffset + nDstYSize == nBlockYSize))
    3043        1323 :         iQuadrantFlag |= (1 << 3);
    3044        4556 :     int l_nFlags = iQuadrantFlag << (4 * (nBand - 1));
    3045        4556 :     int nFullFlags = (1 << (4 * nBands)) - 1;
    3046        4556 :     int nOldFlags = 0;
    3047             : 
    3048       18224 :     for (int i = 1; i <= 3; i++)
    3049             :     {
    3050       13668 :         m_asCachedTilesDesc[i].nRow = -1;
    3051       13668 :         m_asCachedTilesDesc[i].nCol = -1;
    3052       13668 :         m_asCachedTilesDesc[i].nIdxWithinTileData = -1;
    3053             :     }
    3054             : 
    3055        4556 :     int nExistingId = 0;
    3056        4556 :     const char *pszSQL = CPLSPrintf(
    3057             :         "SELECT id, partial_flag, tile_data_band_%d FROM partial_tiles WHERE "
    3058             :         "zoom_level = %d AND tile_row = %d AND tile_column = %d",
    3059             :         nBand, m_nZoomLevel, nRow, nCol);
    3060             : #ifdef DEBUG_VERBOSE
    3061             :     CPLDebug("GPKG", "%s", pszSQL);
    3062             : #endif
    3063        4556 :     sqlite3_stmt *hStmt = nullptr;
    3064        4556 :     int rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
    3065        4556 :     if (rc != SQLITE_OK)
    3066             :     {
    3067           0 :         return CE_Failure;
    3068             :     }
    3069             : 
    3070        4556 :     rc = sqlite3_step(hStmt);
    3071        4556 :     const int nTileBands = m_eDT == GDT_Byte ? 4 : 1;
    3072        4556 :     GByte *pabyTemp = m_pabyCachedTiles + nTileBands * nBandBlockSize;
    3073        4556 :     if (rc == SQLITE_ROW)
    3074             :     {
    3075        4134 :         nExistingId = sqlite3_column_int(hStmt, 0);
    3076             : #ifdef DEBUG_VERBOSE
    3077             :         CPLDebug("GPKG", "Using partial_tile id=%d", nExistingId);
    3078             : #endif
    3079        4134 :         nOldFlags = sqlite3_column_int(hStmt, 1);
    3080        4134 :         CPLAssert(nOldFlags != 0);
    3081        4134 :         if ((nOldFlags & (((1 << 4) - 1) << (4 * (nBand - 1)))) == 0)
    3082             :         {
    3083        1021 :             FillEmptyTileSingleBand(pabyTemp + (nBand - 1) * nBandBlockSize);
    3084             :         }
    3085             :         else
    3086             :         {
    3087        3113 :             CPLAssert(sqlite3_column_bytes(hStmt, 2) ==
    3088             :                       static_cast<int>(nBandBlockSize));
    3089        3113 :             memcpy(pabyTemp + (nBand - 1) * nBandBlockSize,
    3090             :                    sqlite3_column_blob(hStmt, 2), nBandBlockSize);
    3091             :         }
    3092             :     }
    3093             :     else
    3094             :     {
    3095         422 :         FillEmptyTileSingleBand(pabyTemp + (nBand - 1) * nBandBlockSize);
    3096             :     }
    3097        4556 :     sqlite3_finalize(hStmt);
    3098        4556 :     hStmt = nullptr;
    3099             : 
    3100             :     /* Copy the updated rectangle into the full tile */
    3101      556108 :     for (int iY = nDstYOffset; iY < nDstYOffset + nDstYSize; iY++)
    3102             :     {
    3103      551552 :         memcpy(pabyTemp +
    3104      551552 :                    (static_cast<size_t>(nBand - 1) * nBlockXSize * nBlockYSize +
    3105      551552 :                     static_cast<size_t>(iY) * nBlockXSize + nDstXOffset) *
    3106      551552 :                        m_nDTSize,
    3107      551552 :                m_pabyCachedTiles +
    3108      551552 :                    (static_cast<size_t>(nBand - 1) * nBlockXSize * nBlockYSize +
    3109      551552 :                     static_cast<size_t>(iY) * nBlockXSize + nDstXOffset) *
    3110      551552 :                        m_nDTSize,
    3111      551552 :                static_cast<size_t>(nDstXSize) * m_nDTSize);
    3112             :     }
    3113             : 
    3114             : #ifdef notdef
    3115             :     static int nCounter = 1;
    3116             :     GDALDataset *poLogDS =
    3117             :         ((GDALDriver *)GDALGetDriverByName("GTiff"))
    3118             :             ->Create(CPLSPrintf("/tmp/partial_band_%d_%d.tif", 1, nCounter++),
    3119             :                      nBlockXSize, nBlockYSize, nBands, m_eDT, NULL);
    3120             :     poLogDS->RasterIO(GF_Write, 0, 0, nBlockXSize, nBlockYSize,
    3121             :                       pabyTemp + (nBand - 1) * nBandBlockSize, nBlockXSize,
    3122             :                       nBlockYSize, m_eDT, 1, NULL, 0, 0, 0, NULL);
    3123             :     GDALClose(poLogDS);
    3124             : #endif
    3125             : 
    3126        4556 :     if ((nOldFlags & l_nFlags) != 0)
    3127             :     {
    3128           0 :         CPLDebug("GPKG",
    3129             :                  "Rewriting quadrant %d of band %d of tile (row=%d,col=%d)",
    3130             :                  iQuadrantFlag, nBand, nRow, nCol);
    3131             :     }
    3132             : 
    3133        4556 :     l_nFlags |= nOldFlags;
    3134        4556 :     if (l_nFlags == nFullFlags)
    3135             :     {
    3136             : #ifdef DEBUG_VERBOSE
    3137             :         CPLDebug("GPKG", "Got all quadrants for that tile");
    3138             : #endif
    3139        1192 :         for (int iBand = 1; iBand <= nBands; iBand++)
    3140             :         {
    3141         945 :             if (iBand != nBand)
    3142             :             {
    3143         698 :                 pszSQL = CPLSPrintf(
    3144             :                     "SELECT tile_data_band_%d FROM partial_tiles WHERE "
    3145             :                     "id = %d",
    3146             :                     iBand, nExistingId);
    3147             : #ifdef DEBUG_VERBOSE
    3148             :                 CPLDebug("GPKG", "%s", pszSQL);
    3149             : #endif
    3150         698 :                 hStmt = nullptr;
    3151             :                 rc =
    3152         698 :                     SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
    3153         698 :                 if (rc != SQLITE_OK)
    3154             :                 {
    3155           0 :                     return CE_Failure;
    3156             :                 }
    3157             : 
    3158         698 :                 rc = sqlite3_step(hStmt);
    3159         698 :                 if (rc == SQLITE_ROW)
    3160             :                 {
    3161         698 :                     CPLAssert(sqlite3_column_bytes(hStmt, 0) ==
    3162             :                               static_cast<int>(nBandBlockSize));
    3163         698 :                     memcpy(m_pabyCachedTiles + (iBand - 1) * nBandBlockSize,
    3164             :                            sqlite3_column_blob(hStmt, 0), nBandBlockSize);
    3165             :                 }
    3166         698 :                 sqlite3_finalize(hStmt);
    3167         698 :                 hStmt = nullptr;
    3168             :             }
    3169             :             else
    3170             :             {
    3171         247 :                 memcpy(m_pabyCachedTiles + (iBand - 1) * nBandBlockSize,
    3172         247 :                        pabyTemp + (iBand - 1) * nBandBlockSize, nBandBlockSize);
    3173             :             }
    3174             :         }
    3175             : 
    3176         247 :         m_asCachedTilesDesc[0].nRow = nRow;
    3177         247 :         m_asCachedTilesDesc[0].nCol = nCol;
    3178         247 :         m_asCachedTilesDesc[0].nIdxWithinTileData = 0;
    3179         247 :         m_asCachedTilesDesc[0].abBandDirty[0] = true;
    3180         247 :         m_asCachedTilesDesc[0].abBandDirty[1] = true;
    3181         247 :         m_asCachedTilesDesc[0].abBandDirty[2] = true;
    3182         247 :         m_asCachedTilesDesc[0].abBandDirty[3] = true;
    3183             : 
    3184         494 :         pszSQL = CPLSPrintf("UPDATE partial_tiles SET zoom_level = %d, "
    3185             :                             "partial_flag = 0, age = -1 WHERE id = %d",
    3186         247 :                             -1 - m_nZoomLevel, nExistingId);
    3187         247 :         SQLCommand(m_hTempDB, pszSQL);
    3188             : #ifdef DEBUG_VERBOSE
    3189             :         CPLDebug("GPKG", "%s", pszSQL);
    3190             : #endif
    3191         247 :         CPLErr eErr = WriteTile();
    3192             : 
    3193             :         // Call DoPartialFlushOfPartialTilesIfNecessary() after using
    3194             :         // m_pabyCachedTiles as it is going to mess with it.
    3195         247 :         if (DoPartialFlushOfPartialTilesIfNecessary() != CE_None)
    3196           0 :             eErr = CE_None;
    3197             : 
    3198         247 :         return eErr;
    3199             :     }
    3200             : 
    3201        4309 :     if (nExistingId == 0)
    3202             :     {
    3203             :         OGRErr err;
    3204         844 :         pszSQL = CPLSPrintf("SELECT id FROM partial_tiles WHERE "
    3205             :                             "partial_flag = 0 AND zoom_level = %d "
    3206             :                             "AND tile_row = %d AND tile_column = %d",
    3207         422 :                             -1 - m_nZoomLevel, nRow, nCol);
    3208             : #ifdef DEBUG_VERBOSE
    3209             :         CPLDebug("GPKG", "%s", pszSQL);
    3210             : #endif
    3211         422 :         nExistingId = SQLGetInteger(m_hTempDB, pszSQL, &err);
    3212         422 :         if (nExistingId == 0)
    3213             :         {
    3214         418 :             pszSQL =
    3215             :                 "SELECT id FROM partial_tiles WHERE partial_flag = 0 LIMIT 1";
    3216             : #ifdef DEBUG_VERBOSE
    3217             :             CPLDebug("GPKG", "%s", pszSQL);
    3218             : #endif
    3219         418 :             nExistingId = SQLGetInteger(m_hTempDB, pszSQL, &err);
    3220             :         }
    3221             :     }
    3222             : 
    3223        4309 :     const GIntBig nAge = (m_poParentDS) ? m_poParentDS->m_nAge : m_nAge;
    3224        4309 :     if (nExistingId == 0)
    3225             :     {
    3226         338 :         pszSQL = CPLSPrintf(
    3227             :             "INSERT INTO partial_tiles "
    3228             :             "(zoom_level, tile_row, tile_column, tile_data_band_%d, "
    3229             :             "partial_flag, age) VALUES (%d, %d, %d, ?, %d, " CPL_FRMT_GIB ")",
    3230             :             nBand, m_nZoomLevel, nRow, nCol, l_nFlags, nAge);
    3231             :     }
    3232             :     else
    3233             :     {
    3234        3971 :         pszSQL = CPLSPrintf(
    3235             :             "UPDATE partial_tiles SET zoom_level = %d, "
    3236             :             "tile_row = %d, tile_column = %d, "
    3237             :             "tile_data_band_%d = ?, partial_flag = %d, age = " CPL_FRMT_GIB
    3238             :             " WHERE id = %d",
    3239             :             m_nZoomLevel, nRow, nCol, nBand, l_nFlags, nAge, nExistingId);
    3240             :     }
    3241        4309 :     if (m_poParentDS)
    3242        3801 :         m_poParentDS->m_nAge++;
    3243             :     else
    3244         508 :         m_nAge++;
    3245             : 
    3246             : #ifdef DEBUG_VERBOSE
    3247             :     CPLDebug("GPKG", "%s", pszSQL);
    3248             : #endif
    3249             : 
    3250        4309 :     hStmt = nullptr;
    3251        4309 :     rc = SQLPrepareWithError(m_hTempDB, pszSQL, -1, &hStmt, nullptr);
    3252        4309 :     if (rc != SQLITE_OK)
    3253             :     {
    3254           0 :         return CE_Failure;
    3255             :     }
    3256             : 
    3257        4309 :     sqlite3_bind_blob(hStmt, 1, pabyTemp + (nBand - 1) * nBandBlockSize,
    3258             :                       static_cast<int>(nBandBlockSize), SQLITE_TRANSIENT);
    3259        4309 :     rc = sqlite3_step(hStmt);
    3260        4309 :     CPLErr eErr = CE_Failure;
    3261        4309 :     if (rc == SQLITE_DONE)
    3262        4309 :         eErr = CE_None;
    3263             :     else
    3264             :     {
    3265           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    3266             :                  "Failure when inserting partial tile (row=%d,col=%d) at "
    3267             :                  "zoom_level=%d : %s",
    3268             :                  nRow, nCol, m_nZoomLevel, sqlite3_errmsg(m_hTempDB));
    3269             :     }
    3270             : 
    3271        4309 :     sqlite3_finalize(hStmt);
    3272             : 
    3273             :     // Call DoPartialFlushOfPartialTilesIfNecessary() after using
    3274             :     // m_pabyCachedTiles as it is going to mess with it.
    3275        4309 :     if (DoPartialFlushOfPartialTilesIfNecessary() != CE_None)
    3276           0 :         eErr = CE_None;
    3277             : 
    3278        4309 :     return eErr;
    3279             : }
    3280             : 
    3281             : /************************************************************************/
    3282             : /*                         IWriteBlock()                                */
    3283             : /************************************************************************/
    3284             : 
    3285        5849 : CPLErr GDALGPKGMBTilesLikeRasterBand::IWriteBlock(int nBlockXOff,
    3286             :                                                   int nBlockYOff, void *pData)
    3287             : {
    3288             : #ifdef DEBUG_VERBOSE
    3289             :     CPLDebug(
    3290             :         "GPKG",
    3291             :         "IWriteBlock(nBand=%d,nBlockXOff=%d,nBlockYOff=%d,m_nZoomLevel=%d)",
    3292             :         nBand, nBlockXOff, nBlockYOff, m_poTPD->m_nZoomLevel);
    3293             : #endif
    3294        5849 :     if (!m_poTPD->ICanIWriteBlock())
    3295             :     {
    3296           0 :         return CE_Failure;
    3297             :     }
    3298        5849 :     if (m_poTPD->m_poParentDS)
    3299        1139 :         m_poTPD->m_poParentDS->m_bHasModifiedTiles = true;
    3300             :     else
    3301        4710 :         m_poTPD->m_bHasModifiedTiles = true;
    3302             : 
    3303        5849 :     int nRow = nBlockYOff + m_poTPD->m_nShiftYTiles;
    3304        5849 :     int nCol = nBlockXOff + m_poTPD->m_nShiftXTiles;
    3305             : 
    3306        5849 :     int nRowMin = nRow;
    3307        5849 :     int nRowMax = nRowMin;
    3308        5849 :     if (m_poTPD->m_nShiftYPixelsMod)
    3309        1360 :         nRowMax++;
    3310             : 
    3311        5849 :     int nColMin = nCol;
    3312        5849 :     int nColMax = nColMin;
    3313        5849 :     if (m_poTPD->m_nShiftXPixelsMod)
    3314        1322 :         nColMax++;
    3315             : 
    3316        5849 :     CPLErr eErr = CE_None;
    3317             : 
    3318       13058 :     for (nRow = nRowMin; eErr == CE_None && nRow <= nRowMax; nRow++)
    3319             :     {
    3320       17062 :         for (nCol = nColMin; eErr == CE_None && nCol <= nColMax; nCol++)
    3321             :         {
    3322        9853 :             if (nRow < 0 || nCol < 0 || nRow >= m_poTPD->m_nTileMatrixHeight ||
    3323        9713 :                 nCol >= m_poTPD->m_nTileMatrixWidth)
    3324             :             {
    3325         167 :                 continue;
    3326             :             }
    3327             : 
    3328        9686 :             if (m_poTPD->m_nShiftXPixelsMod == 0 &&
    3329        4515 :                 m_poTPD->m_nShiftYPixelsMod == 0)
    3330             :             {
    3331        4457 :                 if (!(nRow == m_poTPD->m_asCachedTilesDesc[0].nRow &&
    3332         250 :                       nCol == m_poTPD->m_asCachedTilesDesc[0].nCol &&
    3333           5 :                       m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData == 0))
    3334             :                 {
    3335        4452 :                     eErr = m_poTPD->WriteTile();
    3336             : 
    3337        4452 :                     m_poTPD->m_asCachedTilesDesc[0].nRow = nRow;
    3338        4452 :                     m_poTPD->m_asCachedTilesDesc[0].nCol = nCol;
    3339        4452 :                     m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData = 0;
    3340             :                 }
    3341             :             }
    3342             : 
    3343             :             // Composite block data into tile, and check if all bands for this
    3344             :             // block are dirty, and if so write the tile
    3345        9686 :             bool bAllDirty = true;
    3346       42233 :             for (int iBand = 1; iBand <= poDS->GetRasterCount(); iBand++)
    3347             :             {
    3348       32547 :                 GDALRasterBlock *poBlock = nullptr;
    3349       32547 :                 GByte *pabySrc = nullptr;
    3350       32547 :                 if (iBand == nBand)
    3351             :                 {
    3352        9686 :                     pabySrc = static_cast<GByte *>(pData);
    3353             :                 }
    3354             :                 else
    3355             :                 {
    3356       22861 :                     if (!(m_poTPD->m_nShiftXPixelsMod == 0 &&
    3357        8359 :                           m_poTPD->m_nShiftYPixelsMod == 0))
    3358       14646 :                         continue;
    3359             : 
    3360             :                     // If the block for this band is not dirty, it might be
    3361             :                     // dirty in cache
    3362        8215 :                     if (m_poTPD->m_asCachedTilesDesc[0].abBandDirty[iBand - 1])
    3363         505 :                         continue;
    3364             :                     else
    3365             :                     {
    3366             :                         poBlock =
    3367        7710 :                             cpl::down_cast<GDALGPKGMBTilesLikeRasterBand *>(
    3368        7710 :                                 poDS->GetRasterBand(iBand))
    3369        7710 :                                 ->TryGetLockedBlockRef(nBlockXOff, nBlockYOff);
    3370        7710 :                         if (poBlock && poBlock->GetDirty())
    3371             :                         {
    3372             :                             pabySrc =
    3373        7670 :                                 static_cast<GByte *>(poBlock->GetDataRef());
    3374        7670 :                             poBlock->MarkClean();
    3375             :                         }
    3376             :                         else
    3377             :                         {
    3378          40 :                             if (poBlock)
    3379           2 :                                 poBlock->DropLock();
    3380          40 :                             bAllDirty = false;
    3381          40 :                             continue;
    3382             :                         }
    3383             :                     }
    3384             :                 }
    3385             : 
    3386       17356 :                 if (m_poTPD->m_nShiftXPixelsMod == 0 &&
    3387       12185 :                     m_poTPD->m_nShiftYPixelsMod == 0)
    3388       12127 :                     m_poTPD->m_asCachedTilesDesc[0].abBandDirty[iBand - 1] =
    3389             :                         true;
    3390             : 
    3391       17356 :                 int nDstXOffset = 0;
    3392       17356 :                 int nDstXSize = nBlockXSize;
    3393       17356 :                 int nDstYOffset = 0;
    3394       17356 :                 int nDstYSize = nBlockYSize;
    3395             :                 // Composite block data into tile data
    3396       17356 :                 if (m_poTPD->m_nShiftXPixelsMod == 0 &&
    3397       12185 :                     m_poTPD->m_nShiftYPixelsMod == 0)
    3398             :                 {
    3399             : 
    3400             : #ifdef DEBUG_VERBOSE
    3401             :                     if (eDataType == GDT_Byte &&
    3402             :                         nBlockXOff * nBlockXSize <=
    3403             :                             nRasterXSize - nBlockXSize &&
    3404             :                         nBlockYOff * nBlockYSize > nRasterYSize - nBlockYSize)
    3405             :                     {
    3406             :                         bool bFoundNonZero = false;
    3407             :                         for (int y = nRasterYSize - nBlockYOff * nBlockYSize;
    3408             :                              y < nBlockYSize; y++)
    3409             :                         {
    3410             :                             for (int x = 0; x < nBlockXSize; x++)
    3411             :                             {
    3412             :                                 if (pabySrc[static_cast<GPtrDiff_t>(y) *
    3413             :                                                 nBlockXSize +
    3414             :                                             x] != 0 &&
    3415             :                                     !bFoundNonZero)
    3416             :                                 {
    3417             :                                     CPLDebug("GPKG",
    3418             :                                              "IWriteBlock(): Found non-zero "
    3419             :                                              "content in ghost part of "
    3420             :                                              "tile(nBand=%d,nBlockXOff=%d,"
    3421             :                                              "nBlockYOff=%d,m_nZoomLevel=%d)\n",
    3422             :                                              iBand, nBlockXOff, nBlockYOff,
    3423             :                                              m_poTPD->m_nZoomLevel);
    3424             :                                     bFoundNonZero = true;
    3425             :                                 }
    3426             :                             }
    3427             :                         }
    3428             :                     }
    3429             : #endif
    3430             : 
    3431       12127 :                     const size_t nBandBlockSize =
    3432       12127 :                         static_cast<size_t>(nBlockXSize) * nBlockYSize *
    3433       12127 :                         m_nDTSize;
    3434       12127 :                     memcpy(m_poTPD->m_pabyCachedTiles +
    3435       12127 :                                (iBand - 1) * nBandBlockSize,
    3436             :                            pabySrc, nBandBlockSize);
    3437             : 
    3438             :                     // Make sure partial blocks are zero'ed outside of the
    3439             :                     // validity area but do that only when know that JPEG will
    3440             :                     // not be used so as to avoid edge effects (although we
    3441             :                     // should probably repeat last pixels if we really want to
    3442             :                     // do that, but that only makes sense if readers only clip
    3443             :                     // to the gpkg_contents extent). Well, ere on the safe side
    3444             :                     // for now
    3445       12127 :                     if (m_poTPD->m_eTF != GPKG_TF_JPEG &&
    3446       11879 :                         (nBlockXOff * nBlockXSize >=
    3447       11879 :                              nRasterXSize - nBlockXSize ||
    3448       11338 :                          nBlockYOff * nBlockYSize >=
    3449       11338 :                              nRasterYSize - nBlockYSize))
    3450             :                     {
    3451        1037 :                         int nXEndValidity =
    3452        1037 :                             nRasterXSize - nBlockXOff * nBlockXSize;
    3453        1037 :                         if (nXEndValidity > nBlockXSize)
    3454         496 :                             nXEndValidity = nBlockXSize;
    3455        1037 :                         int nYEndValidity =
    3456        1037 :                             nRasterYSize - nBlockYOff * nBlockYSize;
    3457        1037 :                         if (nYEndValidity > nBlockYSize)
    3458         294 :                             nYEndValidity = nBlockYSize;
    3459        1037 :                         if (nXEndValidity < nBlockXSize)
    3460             :                         {
    3461        9778 :                             for (int iY = 0; iY < nYEndValidity; iY++)
    3462             :                             {
    3463        9583 :                                 m_poTPD->FillBuffer(
    3464        9583 :                                     m_poTPD->m_pabyCachedTiles +
    3465        9583 :                                         ((static_cast<size_t>(iBand - 1) *
    3466        9583 :                                               nBlockYSize +
    3467        9583 :                                           iY) *
    3468        9583 :                                              nBlockXSize +
    3469        9583 :                                          nXEndValidity) *
    3470        9583 :                                             m_nDTSize,
    3471        9583 :                                     nBlockXSize - nXEndValidity);
    3472             :                             }
    3473             :                         }
    3474        1037 :                         if (nYEndValidity < nBlockYSize)
    3475             :                         {
    3476         235 :                             m_poTPD->FillBuffer(
    3477         235 :                                 m_poTPD->m_pabyCachedTiles +
    3478         235 :                                     (static_cast<size_t>(iBand - 1) *
    3479         235 :                                          nBlockYSize +
    3480         235 :                                      nYEndValidity) *
    3481         235 :                                         nBlockXSize * m_nDTSize,
    3482         235 :                                 static_cast<size_t>(nBlockYSize -
    3483             :                                                     nYEndValidity) *
    3484         235 :                                     nBlockXSize);
    3485             :                         }
    3486       12127 :                     }
    3487             :                 }
    3488             :                 else
    3489             :                 {
    3490        5229 :                     const int nXValid =
    3491        5229 :                         (nBlockXOff * nBlockXSize > nRasterXSize - nBlockXSize)
    3492        5229 :                             ? (nRasterXSize - nBlockXOff * nBlockXSize)
    3493             :                             : nBlockXSize;
    3494        5229 :                     const int nYValid =
    3495        5229 :                         (nBlockYOff * nBlockYSize > nRasterYSize - nBlockYSize)
    3496        5229 :                             ? (nRasterYSize - nBlockYOff * nBlockYSize)
    3497             :                             : nBlockYSize;
    3498             : 
    3499        5229 :                     int nSrcXOffset = 0;
    3500        5229 :                     if (nCol == nColMin)
    3501             :                     {
    3502        2643 :                         nDstXOffset = m_poTPD->m_nShiftXPixelsMod;
    3503        2643 :                         nDstXSize = std::min(
    3504        2643 :                             nXValid, nBlockXSize - m_poTPD->m_nShiftXPixelsMod);
    3505             :                     }
    3506             :                     else
    3507             :                     {
    3508        2586 :                         nDstXOffset = 0;
    3509        2586 :                         if (nXValid > nBlockXSize - m_poTPD->m_nShiftXPixelsMod)
    3510             :                         {
    3511        2212 :                             nDstXSize = nXValid - (nBlockXSize -
    3512        2212 :                                                    m_poTPD->m_nShiftXPixelsMod);
    3513             :                         }
    3514             :                         else
    3515         374 :                             nDstXSize = 0;
    3516        2586 :                         nSrcXOffset = nBlockXSize - m_poTPD->m_nShiftXPixelsMod;
    3517             :                     }
    3518             : 
    3519        5229 :                     int nSrcYOffset = 0;
    3520        5229 :                     if (nRow == nRowMin)
    3521             :                     {
    3522        2607 :                         nDstYOffset = m_poTPD->m_nShiftYPixelsMod;
    3523        2607 :                         nDstYSize = std::min(
    3524        2607 :                             nYValid, nBlockYSize - m_poTPD->m_nShiftYPixelsMod);
    3525             :                     }
    3526             :                     else
    3527             :                     {
    3528        2622 :                         nDstYOffset = 0;
    3529        2622 :                         if (nYValid > nBlockYSize - m_poTPD->m_nShiftYPixelsMod)
    3530             :                         {
    3531        2232 :                             nDstYSize = nYValid - (nBlockYSize -
    3532        2232 :                                                    m_poTPD->m_nShiftYPixelsMod);
    3533             :                         }
    3534             :                         else
    3535         390 :                             nDstYSize = 0;
    3536        2622 :                         nSrcYOffset = nBlockYSize - m_poTPD->m_nShiftYPixelsMod;
    3537             :                     }
    3538             : 
    3539             : #ifdef DEBUG_VERBOSE
    3540             :                     CPLDebug("GPKG",
    3541             :                              "Copy source tile x=%d,w=%d,y=%d,h=%d into "
    3542             :                              "buffer at x=%d,y=%d",
    3543             :                              nDstXOffset, nDstXSize, nDstYOffset, nDstYSize,
    3544             :                              nSrcXOffset, nSrcYOffset);
    3545             : #endif
    3546             : 
    3547        5229 :                     if (nDstXSize > 0 && nDstYSize > 0)
    3548             :                     {
    3549      556108 :                         for (GPtrDiff_t y = 0; y < nDstYSize; y++)
    3550             :                         {
    3551      551552 :                             GByte *pDst = m_poTPD->m_pabyCachedTiles +
    3552      551552 :                                           (static_cast<size_t>(iBand - 1) *
    3553      551552 :                                                nBlockXSize * nBlockYSize +
    3554      551552 :                                            (y + nDstYOffset) * nBlockXSize +
    3555      551552 :                                            nDstXOffset) *
    3556      551552 :                                               m_nDTSize;
    3557      551552 :                             GByte *pSrc =
    3558      551552 :                                 pabySrc + ((y + nSrcYOffset) * nBlockXSize +
    3559      551552 :                                            nSrcXOffset) *
    3560      551552 :                                               m_nDTSize;
    3561      551552 :                             GDALCopyWords(pSrc, eDataType, m_nDTSize, pDst,
    3562             :                                           eDataType, m_nDTSize, nDstXSize);
    3563             :                         }
    3564             :                     }
    3565             :                 }
    3566             : 
    3567       17356 :                 if (poBlock)
    3568        7670 :                     poBlock->DropLock();
    3569             : 
    3570       17356 :                 if (!(m_poTPD->m_nShiftXPixelsMod == 0 &&
    3571       12185 :                       m_poTPD->m_nShiftYPixelsMod == 0))
    3572             :                 {
    3573        5229 :                     m_poTPD->m_asCachedTilesDesc[0].nRow = -1;
    3574        5229 :                     m_poTPD->m_asCachedTilesDesc[0].nCol = -1;
    3575        5229 :                     m_poTPD->m_asCachedTilesDesc[0].nIdxWithinTileData = -1;
    3576        5229 :                     if (nDstXSize > 0 && nDstYSize > 0)
    3577             :                     {
    3578        4556 :                         eErr = m_poTPD->WriteShiftedTile(
    3579             :                             nRow, nCol, iBand, nDstXOffset, nDstYOffset,
    3580             :                             nDstXSize, nDstYSize);
    3581             :                     }
    3582             :                 }
    3583             :             }
    3584             : 
    3585        9686 :             if (m_poTPD->m_nShiftXPixelsMod == 0 &&
    3586        4515 :                 m_poTPD->m_nShiftYPixelsMod == 0)
    3587             :             {
    3588        4457 :                 if (bAllDirty)
    3589             :                 {
    3590        4435 :                     eErr = m_poTPD->WriteTile();
    3591             :                 }
    3592             :             }
    3593             :         }
    3594             :     }
    3595             : 
    3596        5849 :     return eErr;
    3597             : }
    3598             : 
    3599             : /************************************************************************/
    3600             : /*                           GetNoDataValue()                           */
    3601             : /************************************************************************/
    3602             : 
    3603       19907 : double GDALGPKGMBTilesLikeRasterBand::GetNoDataValue(int *pbSuccess)
    3604             : {
    3605       19907 :     if (m_bHasNoData)
    3606             :     {
    3607         401 :         if (pbSuccess)
    3608         400 :             *pbSuccess = TRUE;
    3609         401 :         return m_dfNoDataValue;
    3610             :     }
    3611       19506 :     return GDALPamRasterBand::GetNoDataValue(pbSuccess);
    3612             : }
    3613             : 
    3614             : /************************************************************************/
    3615             : /*                        SetNoDataValueInternal()                      */
    3616             : /************************************************************************/
    3617             : 
    3618          75 : void GDALGPKGMBTilesLikeRasterBand::SetNoDataValueInternal(double dfNoDataValue)
    3619             : {
    3620          75 :     m_bHasNoData = true;
    3621          75 :     m_dfNoDataValue = dfNoDataValue;
    3622          75 : }
    3623             : 
    3624             : /************************************************************************/
    3625             : /*                      GDALGeoPackageRasterBand()                      */
    3626             : /************************************************************************/
    3627             : 
    3628        1817 : GDALGeoPackageRasterBand::GDALGeoPackageRasterBand(
    3629        1817 :     GDALGeoPackageDataset *poDSIn, int nTileWidth, int nTileHeight)
    3630        1817 :     : GDALGPKGMBTilesLikeRasterBand(poDSIn, nTileWidth, nTileHeight)
    3631             : {
    3632        1817 :     poDS = poDSIn;
    3633        1817 : }
    3634             : 
    3635             : /************************************************************************/
    3636             : /*                         GetOverviewCount()                           */
    3637             : /************************************************************************/
    3638             : 
    3639           8 : int GDALGeoPackageRasterBand::GetOverviewCount()
    3640             : {
    3641             :     GDALGeoPackageDataset *poGDS =
    3642           8 :         cpl::down_cast<GDALGeoPackageDataset *>(poDS);
    3643           8 :     return static_cast<int>(poGDS->m_apoOverviewDS.size());
    3644             : }
    3645             : 
    3646             : /************************************************************************/
    3647             : /*                         GetOverviewCount()                           */
    3648             : /************************************************************************/
    3649             : 
    3650          41 : GDALRasterBand *GDALGeoPackageRasterBand::GetOverview(int nIdx)
    3651             : {
    3652          41 :     GDALGeoPackageDataset *poGDS =
    3653             :         reinterpret_cast<GDALGeoPackageDataset *>(poDS);
    3654          41 :     if (nIdx < 0 || nIdx >= static_cast<int>(poGDS->m_apoOverviewDS.size()))
    3655           1 :         return nullptr;
    3656          40 :     return poGDS->m_apoOverviewDS[nIdx]->GetRasterBand(nBand);
    3657             : }
    3658             : 
    3659             : /************************************************************************/
    3660             : /*                           SetNoDataValue()                           */
    3661             : /************************************************************************/
    3662             : 
    3663          23 : CPLErr GDALGeoPackageRasterBand::SetNoDataValue(double dfNoDataValue)
    3664             : {
    3665             :     GDALGeoPackageDataset *poGDS =
    3666          23 :         cpl::down_cast<GDALGeoPackageDataset *>(poDS);
    3667             : 
    3668          23 :     if (eDataType == GDT_Byte)
    3669             :     {
    3670           6 :         if (!(dfNoDataValue >= 0 && dfNoDataValue <= 255 &&
    3671           4 :               static_cast<int>(dfNoDataValue) == dfNoDataValue))
    3672             :         {
    3673           2 :             CPLError(CE_Failure, CPLE_NotSupported,
    3674             :                      "Invalid nodata value for a Byte band: %.17g",
    3675             :                      dfNoDataValue);
    3676           2 :             return CE_Failure;
    3677             :         }
    3678             : 
    3679           8 :         for (int i = 1; i <= poGDS->nBands; ++i)
    3680             :         {
    3681           5 :             if (i != nBand)
    3682             :             {
    3683           2 :                 int bHasNoData = FALSE;
    3684             :                 double dfOtherNoDataValue =
    3685           2 :                     poGDS->GetRasterBand(i)->GetNoDataValue(&bHasNoData);
    3686           2 :                 if (bHasNoData && dfOtherNoDataValue != dfNoDataValue)
    3687             :                 {
    3688           1 :                     CPLError(
    3689             :                         CE_Failure, CPLE_NotSupported,
    3690             :                         "Only the same nodata value can be set on all bands");
    3691           1 :                     return CE_Failure;
    3692             :                 }
    3693             :             }
    3694             :         }
    3695             : 
    3696           3 :         SetNoDataValueInternal(dfNoDataValue);
    3697           3 :         poGDS->m_bMetadataDirty = true;
    3698           3 :         return CE_None;
    3699             :     }
    3700             : 
    3701          17 :     if (std::isnan(dfNoDataValue))
    3702             :     {
    3703           0 :         CPLError(CE_Warning, CPLE_NotSupported,
    3704             :                  "A NaN nodata value cannot be recorded in "
    3705             :                  "gpkg_2d_gridded_coverage_ancillary table");
    3706             :     }
    3707             : 
    3708          17 :     SetNoDataValueInternal(dfNoDataValue);
    3709             : 
    3710          17 :     char *pszSQL = sqlite3_mprintf(
    3711             :         "UPDATE gpkg_2d_gridded_coverage_ancillary SET data_null = ? "
    3712             :         "WHERE tile_matrix_set_name = '%q'",
    3713             :         poGDS->m_osRasterTable.c_str());
    3714          17 :     sqlite3_stmt *hStmt = nullptr;
    3715          17 :     int rc = SQLPrepareWithError(poGDS->IGetDB(), pszSQL, -1, &hStmt, nullptr);
    3716          17 :     if (rc == SQLITE_OK)
    3717             :     {
    3718          17 :         if (poGDS->m_eTF == GPKG_TF_PNG_16BIT)
    3719             :         {
    3720          15 :             if (eDataType == GDT_UInt16 && poGDS->m_dfOffset == 0.0 &&
    3721           6 :                 poGDS->m_dfScale == 1.0 && dfNoDataValue >= 0 &&
    3722           5 :                 dfNoDataValue <= 65535 &&
    3723           5 :                 static_cast<GUInt16>(dfNoDataValue) == dfNoDataValue)
    3724             :             {
    3725           5 :                 poGDS->m_usGPKGNull = static_cast<GUInt16>(dfNoDataValue);
    3726             :             }
    3727             :             else
    3728             :             {
    3729          10 :                 poGDS->m_usGPKGNull = 65535;
    3730             :             }
    3731          15 :             sqlite3_bind_double(hStmt, 1, poGDS->m_usGPKGNull);
    3732             :         }
    3733             :         else
    3734             :         {
    3735           2 :             sqlite3_bind_double(hStmt, 1, static_cast<float>(dfNoDataValue));
    3736             :         }
    3737          17 :         rc = sqlite3_step(hStmt);
    3738          17 :         sqlite3_finalize(hStmt);
    3739             :     }
    3740          17 :     sqlite3_free(pszSQL);
    3741             : 
    3742          17 :     return (rc == SQLITE_OK) ? CE_None : CE_Failure;
    3743             : }
    3744             : 
    3745             : /************************************************************************/
    3746             : /*                         LoadBandMetadata()                           */
    3747             : /************************************************************************/
    3748             : 
    3749        1239 : void GDALGeoPackageRasterBand::LoadBandMetadata()
    3750             : {
    3751             :     GDALGeoPackageDataset *poGDS =
    3752        1239 :         cpl::down_cast<GDALGeoPackageDataset *>(poDS);
    3753             : 
    3754        1239 :     if (m_bHasReadMetadataFromStorage)
    3755        1137 :         return;
    3756             : 
    3757         475 :     m_bHasReadMetadataFromStorage = true;
    3758             : 
    3759         475 :     poGDS->TryLoadXML();
    3760             : 
    3761         475 :     if (!poGDS->HasMetadataTables())
    3762         373 :         return;
    3763             : 
    3764         102 :     char *pszSQL = sqlite3_mprintf(
    3765             :         "SELECT md.metadata, md.md_standard_uri, md.mime_type "
    3766             :         "FROM gpkg_metadata md "
    3767             :         "JOIN gpkg_metadata_reference mdr ON (md.id = mdr.md_file_id ) "
    3768             :         "WHERE "
    3769             :         "mdr.reference_scope = 'table' AND lower(mdr.table_name) = "
    3770             :         "lower('%q') ORDER BY md.id "
    3771             :         "LIMIT 1000",  // to avoid denial of service
    3772             :         poGDS->m_osRasterTable.c_str());
    3773             : 
    3774         102 :     auto oResult = SQLQuery(poGDS->hDB, pszSQL);
    3775         102 :     sqlite3_free(pszSQL);
    3776         102 :     if (!oResult)
    3777             :     {
    3778           0 :         return;
    3779             :     }
    3780             : 
    3781             :     /* GDAL metadata */
    3782         200 :     for (int i = 0; i < oResult->RowCount(); i++)
    3783             :     {
    3784          98 :         const char *pszMetadata = oResult->GetValue(0, i);
    3785          98 :         const char *pszMDStandardURI = oResult->GetValue(1, i);
    3786          98 :         const char *pszMimeType = oResult->GetValue(2, i);
    3787          98 :         if (pszMetadata && pszMDStandardURI && pszMimeType &&
    3788          98 :             EQUAL(pszMDStandardURI, "http://gdal.org") &&
    3789          74 :             EQUAL(pszMimeType, "text/xml"))
    3790             :         {
    3791          74 :             CPLXMLNode *psXMLNode = CPLParseXMLString(pszMetadata);
    3792          74 :             if (psXMLNode)
    3793             :             {
    3794         148 :                 GDALMultiDomainMetadata oLocalMDMD;
    3795          74 :                 oLocalMDMD.XMLInit(psXMLNode, FALSE);
    3796             : 
    3797          74 :                 CSLConstList papszDomainList = oLocalMDMD.GetDomainList();
    3798         208 :                 for (CSLConstList papszIter = papszDomainList;
    3799         208 :                      papszIter && *papszIter; ++papszIter)
    3800             :                 {
    3801         134 :                     if (STARTS_WITH(*papszIter, "BAND_"))
    3802             :                     {
    3803           5 :                         int l_nBand = atoi(*papszIter + strlen("BAND_"));
    3804           5 :                         if (l_nBand >= 1 && l_nBand <= poGDS->GetRasterCount())
    3805             :                         {
    3806             :                             auto l_poBand =
    3807           5 :                                 cpl::down_cast<GDALGeoPackageRasterBand *>(
    3808             :                                     poGDS->GetRasterBand(l_nBand));
    3809           5 :                             l_poBand->m_bHasReadMetadataFromStorage = true;
    3810             : 
    3811          10 :                             char **papszMD = CSLDuplicate(
    3812           5 :                                 oLocalMDMD.GetMetadata(*papszIter));
    3813          10 :                             papszMD = CSLMerge(
    3814             :                                 papszMD,
    3815           5 :                                 GDALGPKGMBTilesLikeRasterBand::GetMetadata(""));
    3816           5 :                             l_poBand->GDALPamRasterBand::SetMetadata(papszMD);
    3817           5 :                             CSLDestroy(papszMD);
    3818             :                         }
    3819             :                     }
    3820             :                 }
    3821             : 
    3822          74 :                 CPLDestroyXMLNode(psXMLNode);
    3823             :             }
    3824             :         }
    3825             :     }
    3826             : }
    3827             : 
    3828             : /************************************************************************/
    3829             : /*                            GetMetadata()                             */
    3830             : /************************************************************************/
    3831             : 
    3832         952 : char **GDALGeoPackageRasterBand::GetMetadata(const char *pszDomain)
    3833             : {
    3834         952 :     GDALGeoPackageDataset *poGDS =
    3835             :         reinterpret_cast<GDALGeoPackageDataset *>(poDS);
    3836         952 :     LoadBandMetadata(); /* force loading from storage if needed */
    3837             : 
    3838          15 :     if (poGDS->eAccess == GA_ReadOnly && eDataType != GDT_Byte &&
    3839          11 :         (pszDomain == nullptr || EQUAL(pszDomain, "")) &&
    3840          11 :         !m_bMinMaxComputedFromTileAncillary &&
    3841         975 :         !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MINIMUM") &&
    3842           8 :         !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MAXIMUM"))
    3843             :     {
    3844           8 :         m_bMinMaxComputedFromTileAncillary = true;
    3845             : 
    3846           8 :         const int nColMin = poGDS->m_nShiftXTiles;
    3847           8 :         const int nColMax =
    3848           8 :             (nRasterXSize - 1 + poGDS->m_nShiftXPixelsMod) / nBlockXSize +
    3849           8 :             poGDS->m_nShiftXTiles;
    3850           8 :         const int nRowMin = poGDS->m_nShiftYTiles;
    3851           8 :         const int nRowMax =
    3852           8 :             (nRasterYSize - 1 + poGDS->m_nShiftYPixelsMod) / nBlockYSize +
    3853           8 :             poGDS->m_nShiftYTiles;
    3854             : 
    3855           8 :         bool bOK = false;
    3856           8 :         if (poGDS->m_nShiftXPixelsMod == 0 && poGDS->m_nShiftYPixelsMod == 0 &&
    3857           8 :             (nRasterXSize % nBlockXSize) == 0 &&
    3858           2 :             (nRasterYSize % nBlockYSize) == 0)
    3859             :         {
    3860             :             // If the area of interest matches entire tiles, then we can
    3861             :             // use tile statistics
    3862           2 :             bOK = true;
    3863             :         }
    3864           6 :         else if (m_bHasNoData)
    3865             :         {
    3866             :             // Otherwise, in the case where we have nodata, we assume that
    3867             :             // if the area of interest is at least larger than the existing
    3868             :             // tiles, the tile statistics will be reliable.
    3869           2 :             char *pszSQL = sqlite3_mprintf(
    3870             :                 "SELECT MIN(tile_column), MAX(tile_column), "
    3871             :                 "MIN(tile_row), MAX(tile_row) FROM \"%w\" "
    3872             :                 "WHERE zoom_level = %d",
    3873             :                 poGDS->m_osRasterTable.c_str(), poGDS->m_nZoomLevel);
    3874           4 :             auto sResult = SQLQuery(poGDS->IGetDB(), pszSQL);
    3875           2 :             if (sResult && sResult->RowCount() == 1)
    3876             :             {
    3877           2 :                 const char *pszMinX = sResult->GetValue(0, 0);
    3878           2 :                 const char *pszMaxX = sResult->GetValue(1, 0);
    3879           2 :                 const char *pszMinY = sResult->GetValue(2, 0);
    3880           2 :                 const char *pszMaxY = sResult->GetValue(3, 0);
    3881           2 :                 if (pszMinX && pszMaxX && pszMinY && pszMaxY)
    3882             :                 {
    3883           2 :                     bOK = atoi(pszMinX) >= nColMin &&
    3884           2 :                           atoi(pszMaxX) <= nColMax &&
    3885           4 :                           atoi(pszMinY) >= nRowMin && atoi(pszMaxY) <= nRowMax;
    3886             :                 }
    3887             :             }
    3888           2 :             sqlite3_free(pszSQL);
    3889             :         }
    3890             : 
    3891           8 :         if (bOK)
    3892             :         {
    3893           4 :             char *pszSQL = sqlite3_mprintf(
    3894             :                 "SELECT MIN(min), MAX(max) FROM "
    3895             :                 "gpkg_2d_gridded_tile_ancillary WHERE tpudt_id "
    3896             :                 "IN (SELECT id FROM \"%w\" WHERE "
    3897             :                 "zoom_level = %d AND "
    3898             :                 "tile_column >= %d AND tile_column <= %d AND "
    3899             :                 "tile_row >= %d AND tile_row <= %d)",
    3900             :                 poGDS->m_osRasterTable.c_str(), poGDS->m_nZoomLevel, nColMin,
    3901             :                 nColMax, nRowMin, nRowMax);
    3902           8 :             auto sResult = SQLQuery(poGDS->IGetDB(), pszSQL);
    3903           4 :             CPLDebug("GPKG", "%s", pszSQL);
    3904           4 :             if (sResult && sResult->RowCount() == 1)
    3905             :             {
    3906           4 :                 const char *pszMin = sResult->GetValue(0, 0);
    3907           4 :                 const char *pszMax = sResult->GetValue(1, 0);
    3908           4 :                 if (pszMin)
    3909             :                 {
    3910           4 :                     m_dfStatsMinFromTileAncillary = CPLAtof(pszMin);
    3911             :                 }
    3912           4 :                 if (pszMax)
    3913             :                 {
    3914           4 :                     m_dfStatsMaxFromTileAncillary = CPLAtof(pszMax);
    3915             :                 }
    3916             :             }
    3917           4 :             sqlite3_free(pszSQL);
    3918             :         }
    3919             :     }
    3920             : 
    3921         705 :     if (m_bAddImplicitStatistics && m_bMinMaxComputedFromTileAncillary &&
    3922           8 :         (pszDomain == nullptr || EQUAL(pszDomain, "")) &&
    3923        1665 :         !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MINIMUM") &&
    3924           8 :         !GDALGPKGMBTilesLikeRasterBand::GetMetadataItem("STATISTICS_MAXIMUM"))
    3925             :     {
    3926             :         m_aosMD.Assign(CSLDuplicate(
    3927           8 :             GDALGPKGMBTilesLikeRasterBand::GetMetadata(pszDomain)));
    3928           8 :         if (!std::isnan(m_dfStatsMinFromTileAncillary))
    3929             :         {
    3930             :             m_aosMD.SetNameValue(
    3931             :                 "STATISTICS_MINIMUM",
    3932           4 :                 CPLSPrintf("%.14g", m_dfStatsMinFromTileAncillary));
    3933             :         }
    3934           8 :         if (!std::isnan(m_dfStatsMaxFromTileAncillary))
    3935             :         {
    3936             :             m_aosMD.SetNameValue(
    3937             :                 "STATISTICS_MAXIMUM",
    3938           4 :                 CPLSPrintf("%.14g", m_dfStatsMaxFromTileAncillary));
    3939             :         }
    3940           8 :         return m_aosMD.List();
    3941             :     }
    3942             : 
    3943         944 :     return GDALGPKGMBTilesLikeRasterBand::GetMetadata(pszDomain);
    3944             : }
    3945             : 
    3946             : /************************************************************************/
    3947             : /*                          GetMetadataItem()                           */
    3948             : /************************************************************************/
    3949             : 
    3950         251 : const char *GDALGeoPackageRasterBand::GetMetadataItem(const char *pszName,
    3951             :                                                       const char *pszDomain)
    3952             : {
    3953         251 :     LoadBandMetadata(); /* force loading from storage if needed */
    3954             : 
    3955         251 :     if (m_bAddImplicitStatistics && eDataType != GDT_Byte &&
    3956           7 :         (pszDomain == nullptr || EQUAL(pszDomain, "")) &&
    3957           4 :         (EQUAL(pszName, "STATISTICS_MINIMUM") ||
    3958           2 :          EQUAL(pszName, "STATISTICS_MAXIMUM")))
    3959             :     {
    3960           2 :         return CSLFetchNameValue(GetMetadata(), pszName);
    3961             :     }
    3962             : 
    3963         249 :     return GDALGPKGMBTilesLikeRasterBand::GetMetadataItem(pszName, pszDomain);
    3964             : }
    3965             : 
    3966             : /************************************************************************/
    3967             : /*                            SetMetadata()                             */
    3968             : /************************************************************************/
    3969             : 
    3970          16 : CPLErr GDALGeoPackageRasterBand::SetMetadata(char **papszMetadata,
    3971             :                                              const char *pszDomain)
    3972             : {
    3973             :     GDALGeoPackageDataset *poGDS =
    3974          16 :         cpl::down_cast<GDALGeoPackageDataset *>(poDS);
    3975          16 :     LoadBandMetadata(); /* force loading from storage if needed */
    3976          16 :     poGDS->m_bMetadataDirty = true;
    3977          72 :     for (CSLConstList papszIter = papszMetadata; papszIter && *papszIter;
    3978             :          ++papszIter)
    3979             :     {
    3980          56 :         if (STARTS_WITH(*papszIter, "STATISTICS_"))
    3981          51 :             m_bStatsMetadataSetInThisSession = true;
    3982             :     }
    3983          16 :     return GDALPamRasterBand::SetMetadata(papszMetadata, pszDomain);
    3984             : }
    3985             : 
    3986             : /************************************************************************/
    3987             : /*                          SetMetadataItem()                           */
    3988             : /************************************************************************/
    3989             : 
    3990          20 : CPLErr GDALGeoPackageRasterBand::SetMetadataItem(const char *pszName,
    3991             :                                                  const char *pszValue,
    3992             :                                                  const char *pszDomain)
    3993             : {
    3994             :     GDALGeoPackageDataset *poGDS =
    3995          20 :         cpl::down_cast<GDALGeoPackageDataset *>(poDS);
    3996          20 :     LoadBandMetadata(); /* force loading from storage if needed */
    3997          20 :     poGDS->m_bMetadataDirty = true;
    3998          20 :     if (STARTS_WITH(pszName, "STATISTICS_"))
    3999          20 :         m_bStatsMetadataSetInThisSession = true;
    4000          20 :     return GDALPamRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
    4001             : }
    4002             : 
    4003             : /************************************************************************/
    4004             : /*                         InvalidateStatistics()                       */
    4005             : /************************************************************************/
    4006             : 
    4007         344 : void GDALGeoPackageRasterBand::InvalidateStatistics()
    4008             : {
    4009         344 :     bool bModified = false;
    4010         688 :     CPLStringList aosMD(CSLDuplicate(GetMetadata()));
    4011         365 :     for (CSLConstList papszIter = GetMetadata(); papszIter && *papszIter;
    4012             :          ++papszIter)
    4013             :     {
    4014          21 :         if (STARTS_WITH(*papszIter, "STATISTICS_"))
    4015             :         {
    4016          20 :             char *pszKey = nullptr;
    4017          20 :             CPLParseNameValue(*papszIter, &pszKey);
    4018          20 :             CPLAssert(pszKey);
    4019          20 :             aosMD.SetNameValue(pszKey, nullptr);
    4020          20 :             CPLFree(pszKey);
    4021          20 :             bModified = true;
    4022             :         }
    4023             :     }
    4024         344 :     if (bModified)
    4025           4 :         SetMetadata(aosMD.List());
    4026         344 : }

Generated by: LCOV version 1.14