LCOV - code coverage report
Current view: top level - frmts/gtiff - cogdriver.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 757 805 94.0 %
Date: 2025-07-09 17:50:03 Functions: 21 21 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  COG Driver
       4             :  * Purpose:  Cloud optimized GeoTIFF write support.
       5             :  * Author:   Even Rouault <even dot rouault at spatialys dot com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2019, Even Rouault <even dot rouault at spatialys dot com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_port.h"
      14             : 
      15             : #include "gdal_priv.h"
      16             : #include "gtiff.h"
      17             : #include "gt_overview.h"
      18             : #include "gdal_utils.h"
      19             : #include "gdalwarper.h"
      20             : #include "cogdriver.h"
      21             : #include "geotiff.h"
      22             : 
      23             : #include "tilematrixset.hpp"
      24             : 
      25             : #include <algorithm>
      26             : #include <memory>
      27             : #include <mutex>
      28             : #include <vector>
      29             : 
      30             : static bool gbHasLZW = false;
      31             : 
      32             : extern "C" CPL_DLL void GDALRegister_COG();
      33             : 
      34             : /************************************************************************/
      35             : /*                        HasZSTDCompression()                          */
      36             : /************************************************************************/
      37             : 
      38         174 : static bool HasZSTDCompression()
      39             : {
      40         174 :     TIFFCodec *codecs = TIFFGetConfiguredCODECs();
      41         174 :     bool bHasZSTD = false;
      42        3478 :     for (TIFFCodec *c = codecs; c->name; ++c)
      43             :     {
      44        3478 :         if (c->scheme == COMPRESSION_ZSTD)
      45             :         {
      46         174 :             bHasZSTD = true;
      47         174 :             break;
      48             :         }
      49             :     }
      50         174 :     _TIFFfree(codecs);
      51         174 :     return bHasZSTD;
      52             : }
      53             : 
      54             : /************************************************************************/
      55             : /*                           GetTmpFilename()                           */
      56             : /************************************************************************/
      57             : 
      58          74 : static CPLString GetTmpFilename(const char *pszFilename, const char *pszExt)
      59             : {
      60             :     const bool bSupportsRandomWrite =
      61          74 :         VSISupportsRandomWrite(pszFilename, false);
      62          74 :     CPLString osTmpFilename;
      63         146 :     if (!bSupportsRandomWrite ||
      64          72 :         CPLGetConfigOption("CPL_TMPDIR", nullptr) != nullptr)
      65             :     {
      66           2 :         osTmpFilename = CPLGenerateTempFilenameSafe(
      67           4 :             CPLGetBasenameSafe(pszFilename).c_str());
      68             :     }
      69             :     else
      70          72 :         osTmpFilename = pszFilename;
      71          74 :     osTmpFilename += '.';
      72          74 :     osTmpFilename += pszExt;
      73          74 :     VSIUnlink(osTmpFilename);
      74          74 :     return osTmpFilename;
      75             : }
      76             : 
      77             : /************************************************************************/
      78             : /*                             GetResampling()                          */
      79             : /************************************************************************/
      80             : 
      81          83 : static const char *GetResampling(GDALDataset *poSrcDS)
      82             : {
      83          83 :     return poSrcDS->GetRasterBand(1)->GetColorTable() ? "NEAREST" : "CUBIC";
      84             : }
      85             : 
      86             : /************************************************************************/
      87             : /*                             GetPredictor()                          */
      88             : /************************************************************************/
      89         294 : static const char *GetPredictor(GDALDataset *poSrcDS, const char *pszPredictor)
      90             : {
      91         294 :     if (pszPredictor == nullptr)
      92           0 :         return nullptr;
      93             : 
      94         294 :     if (EQUAL(pszPredictor, "YES") || EQUAL(pszPredictor, "ON") ||
      95         293 :         EQUAL(pszPredictor, "TRUE"))
      96             :     {
      97           1 :         if (GDALDataTypeIsFloating(
      98           1 :                 poSrcDS->GetRasterBand(1)->GetRasterDataType()))
      99           0 :             return "3";
     100             :         else
     101           1 :             return "2";
     102             :     }
     103         293 :     else if (EQUAL(pszPredictor, "STANDARD") || EQUAL(pszPredictor, "2"))
     104             :     {
     105           1 :         return "2";
     106             :     }
     107         292 :     else if (EQUAL(pszPredictor, "FLOATING_POINT") || EQUAL(pszPredictor, "3"))
     108             :     {
     109           0 :         return "3";
     110             :     }
     111         292 :     return nullptr;
     112             : }
     113             : 
     114             : /************************************************************************/
     115             : /*                            COGGetTargetSRS()                         */
     116             : /************************************************************************/
     117             : 
     118          45 : static bool COGGetTargetSRS(const char *const *papszOptions,
     119             :                             CPLString &osTargetSRS,
     120             :                             std::unique_ptr<gdal::TileMatrixSet> &poTM)
     121             : {
     122          45 :     osTargetSRS = CSLFetchNameValueDef(papszOptions, "TARGET_SRS", "");
     123             :     CPLString osTilingScheme(
     124          90 :         CSLFetchNameValueDef(papszOptions, "TILING_SCHEME", "CUSTOM"));
     125          45 :     if (EQUAL(osTargetSRS, "") && EQUAL(osTilingScheme, "CUSTOM"))
     126           5 :         return false;
     127             : 
     128          40 :     if (!EQUAL(osTilingScheme, "CUSTOM"))
     129             :     {
     130          33 :         poTM = gdal::TileMatrixSet::parse(osTilingScheme);
     131          33 :         if (poTM == nullptr)
     132           0 :             return false;
     133          33 :         if (!poTM->haveAllLevelsSameTopLeft())
     134             :         {
     135           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     136             :                      "Unsupported tiling scheme: not all zoom levels have same "
     137             :                      "top left corner");
     138           0 :             return false;
     139             :         }
     140          33 :         if (!poTM->haveAllLevelsSameTileSize())
     141             :         {
     142           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     143             :                      "Unsupported tiling scheme: not all zoom levels have same "
     144             :                      "tile size");
     145           0 :             return false;
     146             :         }
     147          33 :         if (poTM->hasVariableMatrixWidth())
     148             :         {
     149           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     150             :                      "Unsupported tiling scheme: some levels have variable "
     151             :                      "matrix width");
     152           0 :             return false;
     153             :         }
     154          33 :         if (!osTargetSRS.empty())
     155             :         {
     156           0 :             CPLError(CE_Warning, CPLE_AppDefined, "Ignoring TARGET_SRS option");
     157             :         }
     158          33 :         osTargetSRS = poTM->crs();
     159             : 
     160             :         // "Normalize" SRS as AUTH:CODE
     161          66 :         OGRSpatialReference oTargetSRS;
     162          33 :         oTargetSRS.SetFromUserInput(
     163             :             osTargetSRS,
     164             :             OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
     165          33 :         const char *pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
     166          33 :         const char *pszAuthName = oTargetSRS.GetAuthorityName(nullptr);
     167          33 :         if (pszAuthName && pszAuthCode)
     168             :         {
     169          33 :             osTargetSRS = pszAuthName;
     170          33 :             osTargetSRS += ':';
     171          33 :             osTargetSRS += pszAuthCode;
     172             :         }
     173             :     }
     174             : 
     175          40 :     return true;
     176             : }
     177             : 
     178             : // Used by gdalwarp
     179           8 : bool COGGetTargetSRS(const char *const *papszOptions, CPLString &osTargetSRS)
     180             : {
     181           8 :     std::unique_ptr<gdal::TileMatrixSet> poTM;
     182          16 :     return COGGetTargetSRS(papszOptions, osTargetSRS, poTM);
     183             : }
     184             : 
     185             : // Used by gdalwarp
     186          34 : std::string COGGetResampling(GDALDataset *poSrcDS,
     187             :                              const char *const *papszOptions)
     188             : {
     189             :     return CSLFetchNameValueDef(papszOptions, "WARP_RESAMPLING",
     190             :                                 CSLFetchNameValueDef(papszOptions, "RESAMPLING",
     191          34 :                                                      GetResampling(poSrcDS)));
     192             : }
     193             : 
     194             : /************************************************************************/
     195             : /*                     COGGetWarpingCharacteristics()                   */
     196             : /************************************************************************/
     197             : 
     198          37 : static bool COGGetWarpingCharacteristics(
     199             :     GDALDataset *poSrcDS, const char *const *papszOptions,
     200             :     CPLString &osResampling, CPLString &osTargetSRS, int &nXSize, int &nYSize,
     201             :     double &dfMinX, double &dfMinY, double &dfMaxX, double &dfMaxY,
     202             :     double &dfRes, std::unique_ptr<gdal::TileMatrixSet> &poTM, int &nZoomLevel,
     203             :     int &nAlignedLevels)
     204             : {
     205          37 :     if (!COGGetTargetSRS(papszOptions, osTargetSRS, poTM))
     206           2 :         return false;
     207             : 
     208          70 :     CPLStringList aosTO;
     209          35 :     aosTO.SetNameValue("DST_SRS", osTargetSRS);
     210          35 :     void *hTransformArg = nullptr;
     211             : 
     212          70 :     OGRSpatialReference oTargetSRS;
     213          35 :     oTargetSRS.SetFromUserInput(
     214             :         osTargetSRS,
     215             :         OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
     216          35 :     const char *pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
     217          35 :     const int nEPSGCode = pszAuthCode ? atoi(pszAuthCode) : 0;
     218             : 
     219             :     // Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
     220             :     // ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
     221             :     // EPSG:3857.
     222          35 :     GDALGeoTransform srcGT;
     223          35 :     std::unique_ptr<GDALDataset> poTmpDS;
     224          27 :     if (nEPSGCode == 3857 && poSrcDS->GetGeoTransform(srcGT) == CE_None &&
     225          62 :         srcGT[2] == 0 && srcGT[4] == 0 && srcGT[5] < 0)
     226             :     {
     227          27 :         const auto poSrcSRS = poSrcDS->GetSpatialRef();
     228          31 :         if (poSrcSRS && poSrcSRS->IsGeographic() &&
     229           4 :             !poSrcSRS->IsDerivedGeographic())
     230             :         {
     231           4 :             double maxLat = srcGT[3];
     232           4 :             double minLat = srcGT[3] + poSrcDS->GetRasterYSize() * srcGT[5];
     233             :             // Corresponds to the latitude of below MAX_GM
     234           4 :             constexpr double MAX_LAT = 85.0511287798066;
     235           4 :             bool bModified = false;
     236           4 :             if (maxLat > MAX_LAT)
     237             :             {
     238           4 :                 maxLat = MAX_LAT;
     239           4 :                 bModified = true;
     240             :             }
     241           4 :             if (minLat < -MAX_LAT)
     242             :             {
     243           4 :                 minLat = -MAX_LAT;
     244           4 :                 bModified = true;
     245             :             }
     246           4 :             if (bModified)
     247             :             {
     248           4 :                 CPLStringList aosOptions;
     249           4 :                 aosOptions.AddString("-of");
     250           4 :                 aosOptions.AddString("VRT");
     251           4 :                 aosOptions.AddString("-projwin");
     252           4 :                 aosOptions.AddString(CPLSPrintf("%.17g", srcGT[0]));
     253           4 :                 aosOptions.AddString(CPLSPrintf("%.17g", maxLat));
     254             :                 aosOptions.AddString(CPLSPrintf(
     255           4 :                     "%.17g", srcGT[0] + poSrcDS->GetRasterXSize() * srcGT[1]));
     256           4 :                 aosOptions.AddString(CPLSPrintf("%.17g", minLat));
     257             :                 auto psOptions =
     258           4 :                     GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     259           4 :                 poTmpDS.reset(GDALDataset::FromHandle(GDALTranslate(
     260             :                     "", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
     261           4 :                 GDALTranslateOptionsFree(psOptions);
     262           4 :                 if (poTmpDS)
     263             :                 {
     264           4 :                     hTransformArg = GDALCreateGenImgProjTransformer2(
     265           4 :                         GDALDataset::FromHandle(poTmpDS.get()), nullptr,
     266           4 :                         aosTO.List());
     267           4 :                     if (hTransformArg == nullptr)
     268             :                     {
     269           0 :                         return false;
     270             :                     }
     271             :                 }
     272             :             }
     273             :         }
     274             :     }
     275          35 :     if (hTransformArg == nullptr)
     276             :     {
     277             :         hTransformArg =
     278          31 :             GDALCreateGenImgProjTransformer2(poSrcDS, nullptr, aosTO.List());
     279          31 :         if (hTransformArg == nullptr)
     280             :         {
     281           0 :             return false;
     282             :         }
     283             :     }
     284             : 
     285          35 :     GDALTransformerInfo *psInfo =
     286             :         static_cast<GDALTransformerInfo *>(hTransformArg);
     287             :     double adfGeoTransform[6];
     288             :     double adfExtent[4];
     289             : 
     290          35 :     if (GDALSuggestedWarpOutput2(poTmpDS ? poTmpDS.get() : poSrcDS,
     291             :                                  psInfo->pfnTransform, hTransformArg,
     292             :                                  adfGeoTransform, &nXSize, &nYSize, adfExtent,
     293          35 :                                  0) != CE_None)
     294             :     {
     295           0 :         GDALDestroyGenImgProjTransformer(hTransformArg);
     296           0 :         return false;
     297             :     }
     298             : 
     299          35 :     GDALDestroyGenImgProjTransformer(hTransformArg);
     300          35 :     hTransformArg = nullptr;
     301          35 :     poTmpDS.reset();
     302             : 
     303          35 :     dfMinX = adfExtent[0];
     304          35 :     dfMinY = adfExtent[1];
     305          35 :     dfMaxX = adfExtent[2];
     306          35 :     dfMaxY = adfExtent[3];
     307          35 :     dfRes = adfGeoTransform[1];
     308             : 
     309          70 :     const CPLString osExtent(CSLFetchNameValueDef(papszOptions, "EXTENT", ""));
     310          70 :     const CPLString osRes(CSLFetchNameValueDef(papszOptions, "RES", ""));
     311          35 :     if (poTM)
     312             :     {
     313          28 :         if (!osExtent.empty())
     314             :         {
     315           0 :             CPLError(CE_Warning, CPLE_AppDefined, "Ignoring EXTENT option");
     316             :         }
     317          28 :         if (!osRes.empty())
     318             :         {
     319           0 :             CPLError(CE_Warning, CPLE_AppDefined, "Ignoring RES option");
     320             :         }
     321             :         const bool bInvertAxis =
     322          56 :             oTargetSRS.EPSGTreatsAsLatLong() != FALSE ||
     323          28 :             oTargetSRS.EPSGTreatsAsNorthingEasting() != FALSE;
     324             : 
     325          28 :         const auto &bbox = poTM->bbox();
     326          28 :         if (bbox.mCrs == poTM->crs())
     327             :         {
     328          28 :             if (dfMaxX <
     329          56 :                     (bInvertAxis ? bbox.mLowerCornerY : bbox.mLowerCornerX) ||
     330          28 :                 dfMinX >
     331          56 :                     (bInvertAxis ? bbox.mUpperCornerY : bbox.mUpperCornerX) ||
     332          28 :                 dfMaxY <
     333          56 :                     (bInvertAxis ? bbox.mLowerCornerX : bbox.mLowerCornerY) ||
     334          28 :                 dfMinY >
     335          28 :                     (bInvertAxis ? bbox.mUpperCornerX : bbox.mUpperCornerY))
     336             :             {
     337           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     338             :                          "Raster extent completely outside of tile matrix set "
     339             :                          "bounding box");
     340           2 :                 return false;
     341             :             }
     342             :         }
     343             : 
     344          28 :         const auto &tmList = poTM->tileMatrixList();
     345          28 :         const int nBlockSize = atoi(CSLFetchNameValueDef(
     346          28 :             papszOptions, "BLOCKSIZE", CPLSPrintf("%d", tmList[0].mTileWidth)));
     347          28 :         dfRes = 0.0;
     348             : 
     349             :         const char *pszZoomLevel =
     350          28 :             CSLFetchNameValue(papszOptions, "ZOOM_LEVEL");
     351          28 :         if (pszZoomLevel)
     352             :         {
     353           3 :             nZoomLevel = atoi(pszZoomLevel);
     354           3 :             if (nZoomLevel < 0 || nZoomLevel >= static_cast<int>(tmList.size()))
     355             :             {
     356           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
     357             :                          "Invalid zoom level: should be in [0,%d]",
     358           2 :                          static_cast<int>(tmList.size()) - 1);
     359           2 :                 return false;
     360             :             }
     361             :         }
     362             :         else
     363             :         {
     364          25 :             double dfComputedRes = adfGeoTransform[1];
     365          25 :             double dfPrevRes = 0.0;
     366         280 :             for (; nZoomLevel < static_cast<int>(tmList.size()); nZoomLevel++)
     367             :             {
     368         280 :                 dfRes = tmList[nZoomLevel].mResX * tmList[0].mTileWidth /
     369             :                         nBlockSize;
     370         280 :                 if (dfComputedRes > dfRes ||
     371         260 :                     fabs(dfComputedRes - dfRes) / dfRes <= 1e-8)
     372             :                     break;
     373         255 :                 dfPrevRes = dfRes;
     374             :             }
     375          25 :             if (nZoomLevel == static_cast<int>(tmList.size()))
     376             :             {
     377           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     378             :                          "Could not find an appropriate zoom level");
     379           0 :                 return false;
     380             :             }
     381             : 
     382          25 :             if (nZoomLevel > 0 && fabs(dfComputedRes - dfRes) / dfRes > 1e-8)
     383             :             {
     384          19 :                 const char *pszZoomLevelStrategy = CSLFetchNameValueDef(
     385             :                     papszOptions, "ZOOM_LEVEL_STRATEGY", "AUTO");
     386          19 :                 if (EQUAL(pszZoomLevelStrategy, "LOWER"))
     387             :                 {
     388           1 :                     nZoomLevel--;
     389             :                 }
     390          18 :                 else if (EQUAL(pszZoomLevelStrategy, "UPPER"))
     391             :                 {
     392             :                     /* do nothing */
     393             :                 }
     394             :                 else
     395             :                 {
     396          17 :                     if (dfPrevRes / dfComputedRes < dfComputedRes / dfRes)
     397          15 :                         nZoomLevel--;
     398             :                 }
     399             :             }
     400             :         }
     401          26 :         CPLDebug("COG", "Using ZOOM_LEVEL %d", nZoomLevel);
     402          26 :         dfRes = tmList[nZoomLevel].mResX * tmList[0].mTileWidth / nBlockSize;
     403             : 
     404             :         const double dfOriX =
     405          26 :             bInvertAxis ? tmList[0].mTopLeftY : tmList[0].mTopLeftX;
     406             :         const double dfOriY =
     407          26 :             bInvertAxis ? tmList[0].mTopLeftX : tmList[0].mTopLeftY;
     408          26 :         const double dfTileExtent = dfRes * nBlockSize;
     409          26 :         constexpr double TOLERANCE_IN_PIXEL = 0.499;
     410          26 :         const double dfEps = TOLERANCE_IN_PIXEL * dfRes;
     411          26 :         int nTLTileX = static_cast<int>(
     412          26 :             std::floor((dfMinX - dfOriX + dfEps) / dfTileExtent));
     413          26 :         int nTLTileY = static_cast<int>(
     414          26 :             std::floor((dfOriY - dfMaxY + dfEps) / dfTileExtent));
     415          26 :         int nBRTileX = static_cast<int>(
     416          26 :             std::ceil((dfMaxX - dfOriX - dfEps) / dfTileExtent));
     417          26 :         int nBRTileY = static_cast<int>(
     418          26 :             std::ceil((dfOriY - dfMinY - dfEps) / dfTileExtent));
     419             : 
     420          26 :         nAlignedLevels =
     421          26 :             std::min(std::min(10, atoi(CSLFetchNameValueDef(
     422             :                                       papszOptions, "ALIGNED_LEVELS", "0"))),
     423          26 :                      nZoomLevel);
     424          26 :         int nAccDivisor = 1;
     425          34 :         for (int i = 0; i < nAlignedLevels - 1; i++)
     426             :         {
     427           8 :             const int nCurLevel = nZoomLevel - i;
     428             :             const double dfResRatio =
     429           8 :                 tmList[nCurLevel - 1].mResX / tmList[nCurLevel].mResX;
     430             :             // Magical number that has a great number of divisors
     431             :             // For example if previous scale denom was 50K and current one
     432             :             // is 20K, then dfResRatio = 2.5 and dfScaledInvResRatio = 24
     433             :             // We must then simplify 60 / 24 as 5 / 2, and make sure to
     434             :             // align tile coordinates on multiple of the 5 numerator
     435           8 :             constexpr int MAGICAL = 60;
     436           8 :             const double dfScaledInvResRatio = MAGICAL / dfResRatio;
     437          16 :             if (dfScaledInvResRatio < 1 || dfScaledInvResRatio > 60 ||
     438           8 :                 std::abs(std::round(dfScaledInvResRatio) -
     439             :                          dfScaledInvResRatio) > 1e-10)
     440             :             {
     441           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     442             :                          "Unsupported ratio of resolution for "
     443             :                          "ALIGNED_LEVELS between zoom level %d and %d = %g",
     444             :                          nCurLevel - 1, nCurLevel, dfResRatio);
     445           0 :                 return false;
     446             :             }
     447           8 :             const int nScaledInvResRatio =
     448           8 :                 static_cast<int>(std::round(dfScaledInvResRatio));
     449           8 :             int nNumerator = 0;
     450          32 :             for (int nDivisor = nScaledInvResRatio; nDivisor >= 2; --nDivisor)
     451             :             {
     452          32 :                 if ((MAGICAL % nDivisor) == 0 &&
     453          12 :                     (nScaledInvResRatio % nDivisor) == 0)
     454             :                 {
     455           8 :                     nNumerator = MAGICAL / nDivisor;
     456           8 :                     break;
     457             :                 }
     458             :             }
     459           8 :             if (nNumerator == 0)
     460             :             {
     461           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     462             :                          "Unsupported ratio of resolution for "
     463             :                          "ALIGNED_LEVELS between zoom level %d and %d = %g",
     464             :                          nCurLevel - 1, nCurLevel, dfResRatio);
     465           0 :                 return false;
     466             :             }
     467           8 :             nAccDivisor *= nNumerator;
     468             :         }
     469          26 :         if (nAccDivisor > 1)
     470             :         {
     471           4 :             nTLTileX = (nTLTileX / nAccDivisor) * nAccDivisor;
     472           4 :             nTLTileY = (nTLTileY / nAccDivisor) * nAccDivisor;
     473           4 :             nBRTileY = DIV_ROUND_UP(nBRTileY, nAccDivisor) * nAccDivisor;
     474           4 :             nBRTileX = DIV_ROUND_UP(nBRTileX, nAccDivisor) * nAccDivisor;
     475             :         }
     476             : 
     477          26 :         if (nTLTileX < 0 || nTLTileY < 0 ||
     478          74 :             nBRTileX > tmList[nZoomLevel].mMatrixWidth ||
     479          22 :             nBRTileY > tmList[nZoomLevel].mMatrixHeight)
     480             :         {
     481           4 :             CPLError(CE_Warning, CPLE_AppDefined,
     482             :                      "Raster extent partially outside of tile matrix "
     483             :                      "bounding box. Clamping it to it");
     484             :         }
     485          26 :         nTLTileX = std::max(0, nTLTileX);
     486          26 :         nTLTileY = std::max(0, nTLTileY);
     487          26 :         nBRTileX = std::min(tmList[nZoomLevel].mMatrixWidth, nBRTileX);
     488          26 :         nBRTileY = std::min(tmList[nZoomLevel].mMatrixHeight, nBRTileY);
     489             : 
     490          26 :         dfMinX = dfOriX + nTLTileX * dfTileExtent;
     491          26 :         dfMinY = dfOriY - nBRTileY * dfTileExtent;
     492          26 :         dfMaxX = dfOriX + nBRTileX * dfTileExtent;
     493          26 :         dfMaxY = dfOriY - nTLTileY * dfTileExtent;
     494             :     }
     495           7 :     else if (!osExtent.empty() || !osRes.empty())
     496             :     {
     497           1 :         CPLStringList aosTokens(CSLTokenizeString2(osExtent, ",", 0));
     498           1 :         if (aosTokens.size() != 4)
     499             :         {
     500           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for EXTENT");
     501           0 :             return false;
     502             :         }
     503           1 :         dfMinX = CPLAtof(aosTokens[0]);
     504           1 :         dfMinY = CPLAtof(aosTokens[1]);
     505           1 :         dfMaxX = CPLAtof(aosTokens[2]);
     506           1 :         dfMaxY = CPLAtof(aosTokens[3]);
     507           1 :         if (!osRes.empty())
     508           1 :             dfRes = CPLAtof(osRes);
     509             :     }
     510             : 
     511          33 :     nXSize = static_cast<int>(std::round((dfMaxX - dfMinX) / dfRes));
     512          33 :     nYSize = static_cast<int>(std::round((dfMaxY - dfMinY) / dfRes));
     513             : 
     514          33 :     osResampling = COGGetResampling(poSrcDS, papszOptions);
     515             : 
     516          33 :     return true;
     517             : }
     518             : 
     519             : // Used by gdalwarp
     520           4 : bool COGGetWarpingCharacteristics(GDALDataset *poSrcDS,
     521             :                                   const char *const *papszOptions,
     522             :                                   CPLString &osResampling,
     523             :                                   CPLString &osTargetSRS, int &nXSize,
     524             :                                   int &nYSize, double &dfMinX, double &dfMinY,
     525             :                                   double &dfMaxX, double &dfMaxY)
     526             : {
     527           4 :     std::unique_ptr<gdal::TileMatrixSet> poTM;
     528           4 :     int nZoomLevel = 0;
     529           4 :     int nAlignedLevels = 0;
     530             :     double dfRes;
     531           4 :     return COGGetWarpingCharacteristics(poSrcDS, papszOptions, osResampling,
     532             :                                         osTargetSRS, nXSize, nYSize, dfMinX,
     533             :                                         dfMinY, dfMaxX, dfMaxY, dfRes, poTM,
     534           8 :                                         nZoomLevel, nAlignedLevels);
     535             : }
     536             : 
     537             : /************************************************************************/
     538             : /*                        COGHasWarpingOptions()                        */
     539             : /************************************************************************/
     540             : 
     541         154 : bool COGHasWarpingOptions(CSLConstList papszOptions)
     542             : {
     543         301 :     return CSLFetchNameValue(papszOptions, "TARGET_SRS") != nullptr ||
     544         301 :            !EQUAL(CSLFetchNameValueDef(papszOptions, "TILING_SCHEME", "CUSTOM"),
     545             :                   "CUSTOM");
     546             : }
     547             : 
     548             : /************************************************************************/
     549             : /*                      COGRemoveWarpingOptions()                       */
     550             : /************************************************************************/
     551             : 
     552           2 : void COGRemoveWarpingOptions(CPLStringList &aosOptions)
     553             : {
     554           2 :     aosOptions.SetNameValue("TARGET_SRS", nullptr);
     555           2 :     aosOptions.SetNameValue("TILING_SCHEME", nullptr);
     556           2 :     aosOptions.SetNameValue("EXTENT", nullptr);
     557           2 :     aosOptions.SetNameValue("RES", nullptr);
     558           2 :     aosOptions.SetNameValue("ALIGNED_LEVELS", nullptr);
     559           2 :     aosOptions.SetNameValue("ZOOM_LEVEL_STRATEGY", nullptr);
     560           2 : }
     561             : 
     562             : /************************************************************************/
     563             : /*                        CreateReprojectedDS()                         */
     564             : /************************************************************************/
     565             : 
     566          26 : static std::unique_ptr<GDALDataset> CreateReprojectedDS(
     567             :     const char *pszDstFilename, GDALDataset *poSrcDS,
     568             :     const char *const *papszOptions, const CPLString &osResampling,
     569             :     const CPLString &osTargetSRS, const int nXSize, const int nYSize,
     570             :     const double dfMinX, const double dfMinY, const double dfMaxX,
     571             :     const double dfMaxY, const double dfRes, GDALProgressFunc pfnProgress,
     572             :     void *pProgressData, double &dfCurPixels, double &dfTotalPixelsToProcess)
     573             : {
     574          26 :     char **papszArg = nullptr;
     575             :     // We could have done a warped VRT, but overview building on it might be
     576             :     // slow, so materialize as GTiff
     577          26 :     papszArg = CSLAddString(papszArg, "-of");
     578          26 :     papszArg = CSLAddString(papszArg, "GTiff");
     579          26 :     papszArg = CSLAddString(papszArg, "-co");
     580          26 :     papszArg = CSLAddString(papszArg, "TILED=YES");
     581          26 :     papszArg = CSLAddString(papszArg, "-co");
     582          26 :     papszArg = CSLAddString(papszArg, "SPARSE_OK=YES");
     583          26 :     const char *pszBIGTIFF = CSLFetchNameValue(papszOptions, "BIGTIFF");
     584          26 :     if (pszBIGTIFF)
     585             :     {
     586           0 :         papszArg = CSLAddString(papszArg, "-co");
     587           0 :         papszArg = CSLAddString(papszArg,
     588           0 :                                 (CPLString("BIGTIFF=") + pszBIGTIFF).c_str());
     589             :     }
     590          26 :     papszArg = CSLAddString(papszArg, "-co");
     591          26 :     papszArg = CSLAddString(papszArg, HasZSTDCompression() ? "COMPRESS=ZSTD"
     592             :                                                            : "COMPRESS=LZW");
     593          26 :     papszArg = CSLAddString(papszArg, "-t_srs");
     594          26 :     papszArg = CSLAddString(papszArg, osTargetSRS);
     595          26 :     papszArg = CSLAddString(papszArg, "-te");
     596          26 :     papszArg = CSLAddString(papszArg, CPLSPrintf("%.17g", dfMinX));
     597          26 :     papszArg = CSLAddString(papszArg, CPLSPrintf("%.17g", dfMinY));
     598          26 :     papszArg = CSLAddString(papszArg, CPLSPrintf("%.17g", dfMaxX));
     599          26 :     papszArg = CSLAddString(papszArg, CPLSPrintf("%.17g", dfMaxY));
     600          26 :     papszArg = CSLAddString(papszArg, "-ts");
     601          26 :     papszArg = CSLAddString(papszArg, CPLSPrintf("%d", nXSize));
     602          26 :     papszArg = CSLAddString(papszArg, CPLSPrintf("%d", nYSize));
     603             : 
     604             :     // to be kept in sync with gdalwarp_lib.cpp
     605          26 :     constexpr double RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP = 1e-8;
     606          26 :     if (fabs((dfMaxX - dfMinX) / dfRes - nXSize) <=
     607          26 :             RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP &&
     608          26 :         fabs((dfMaxY - dfMinY) / dfRes - nYSize) <=
     609             :             RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP)
     610             :     {
     611             :         // Try to produce exactly square pixels
     612          26 :         papszArg = CSLAddString(papszArg, "-tr");
     613          26 :         papszArg = CSLAddString(papszArg, CPLSPrintf("%.17g", dfRes));
     614          26 :         papszArg = CSLAddString(papszArg, CPLSPrintf("%.17g", dfRes));
     615             :     }
     616             :     else
     617             :     {
     618           0 :         CPLDebug("COG", "Cannot pass -tr option to GDALWarp() due to extent, "
     619             :                         "size and resolution not consistent enough");
     620             :     }
     621             : 
     622          26 :     int bHasNoData = FALSE;
     623          26 :     poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoData);
     624          52 :     if (!bHasNoData &&
     625          26 :         CPLTestBool(CSLFetchNameValueDef(papszOptions, "ADD_ALPHA", "YES")))
     626             :     {
     627          26 :         papszArg = CSLAddString(papszArg, "-dstalpha");
     628             :     }
     629          26 :     papszArg = CSLAddString(papszArg, "-r");
     630          26 :     papszArg = CSLAddString(papszArg, osResampling);
     631          26 :     papszArg = CSLAddString(papszArg, "-wo");
     632          26 :     papszArg = CSLAddString(papszArg, "SAMPLE_GRID=YES");
     633          26 :     const char *pszNumThreads = CSLFetchNameValue(papszOptions, "NUM_THREADS");
     634          26 :     if (pszNumThreads)
     635             :     {
     636           0 :         papszArg = CSLAddString(papszArg, "-wo");
     637           0 :         papszArg = CSLAddString(
     638           0 :             papszArg, (CPLString("NUM_THREADS=") + pszNumThreads).c_str());
     639             :     }
     640             : 
     641          26 :     const auto poFirstBand = poSrcDS->GetRasterBand(1);
     642          26 :     const bool bHasMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
     643             : 
     644          26 :     const int nBands = poSrcDS->GetRasterCount();
     645             :     const char *pszOverviews =
     646          26 :         CSLFetchNameValueDef(papszOptions, "OVERVIEWS", "AUTO");
     647          52 :     const bool bUseExistingOrNone = EQUAL(pszOverviews, "FORCE_USE_EXISTING") ||
     648          26 :                                     EQUAL(pszOverviews, "NONE");
     649          26 :     dfTotalPixelsToProcess =
     650          26 :         double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) +
     651          26 :         ((bHasMask && !bUseExistingOrNone) ? double(nXSize) * nYSize / 3 : 0) +
     652          26 :         (!bUseExistingOrNone ? double(nXSize) * nYSize * nBands / 3 : 0) +
     653          26 :         double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) * 4. / 3;
     654             : 
     655          26 :     auto psOptions = GDALWarpAppOptionsNew(papszArg, nullptr);
     656          26 :     CSLDestroy(papszArg);
     657          26 :     if (psOptions == nullptr)
     658           1 :         return nullptr;
     659             : 
     660          25 :     const double dfNextPixels =
     661          25 :         double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0));
     662          50 :     void *pScaledProgress = GDALCreateScaledProgress(
     663          25 :         dfCurPixels / dfTotalPixelsToProcess,
     664          25 :         dfNextPixels / dfTotalPixelsToProcess, pfnProgress, pProgressData);
     665          25 :     dfCurPixels = dfNextPixels;
     666             : 
     667          25 :     CPLDebug("COG", "Reprojecting source dataset: start");
     668          25 :     GDALWarpAppOptionsSetProgress(psOptions, GDALScaledProgress,
     669             :                                   pScaledProgress);
     670          50 :     CPLString osTmpFile(GetTmpFilename(pszDstFilename, "warped.tif.tmp"));
     671          25 :     auto hSrcDS = GDALDataset::ToHandle(poSrcDS);
     672             : 
     673          25 :     std::unique_ptr<CPLConfigOptionSetter> poWarpThreadSetter;
     674          25 :     if (pszNumThreads)
     675             :     {
     676           0 :         poWarpThreadSetter.reset(new CPLConfigOptionSetter(
     677           0 :             "GDAL_NUM_THREADS", pszNumThreads, false));
     678             :     }
     679             : 
     680          25 :     auto hRet = GDALWarp(osTmpFile, nullptr, 1, &hSrcDS, psOptions, nullptr);
     681          25 :     GDALWarpAppOptionsFree(psOptions);
     682          25 :     CPLDebug("COG", "Reprojecting source dataset: end");
     683             : 
     684          25 :     GDALDestroyScaledProgress(pScaledProgress);
     685             : 
     686          25 :     return std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(hRet));
     687             : }
     688             : 
     689             : /************************************************************************/
     690             : /*                            GDALCOGCreator                            */
     691             : /************************************************************************/
     692             : 
     693             : struct GDALCOGCreator final
     694             : {
     695             :     std::unique_ptr<GDALDataset> m_poReprojectedDS{};
     696             :     std::unique_ptr<GDALDataset> m_poRGBMaskDS{};
     697             :     std::unique_ptr<GDALDataset> m_poVRTWithOrWithoutStats{};
     698             :     CPLString m_osTmpOverviewFilename{};
     699             :     CPLString m_osTmpMskOverviewFilename{};
     700             : 
     701             :     ~GDALCOGCreator();
     702             : 
     703             :     GDALDataset *Create(const char *pszFilename, GDALDataset *const poSrcDS,
     704             :                         char **papszOptions, GDALProgressFunc pfnProgress,
     705             :                         void *pProgressData);
     706             : };
     707             : 
     708             : /************************************************************************/
     709             : /*                    GDALCOGCreator::~GDALCOGCreator()                 */
     710             : /************************************************************************/
     711             : 
     712         153 : GDALCOGCreator::~GDALCOGCreator()
     713             : {
     714             :     // Destroy m_poRGBMaskDS before m_poReprojectedDS since the former
     715             :     // may reference the later
     716         153 :     m_poRGBMaskDS.reset();
     717             : 
     718             :     // Config option just for testing purposes
     719             :     const bool bDeleteTempFiles =
     720         153 :         CPLTestBool(CPLGetConfigOption("COG_DELETE_TEMP_FILES", "YES"));
     721         153 :     if (bDeleteTempFiles)
     722             :     {
     723         152 :         if (m_poReprojectedDS)
     724             :         {
     725          50 :             CPLString osProjectedDSName(m_poReprojectedDS->GetDescription());
     726          25 :             m_poReprojectedDS.reset();
     727          25 :             VSIUnlink(osProjectedDSName);
     728             :         }
     729         152 :         if (!m_osTmpOverviewFilename.empty())
     730             :         {
     731          38 :             VSIUnlink(m_osTmpOverviewFilename);
     732             :         }
     733         152 :         if (!m_osTmpMskOverviewFilename.empty())
     734             :         {
     735           9 :             VSIUnlink(m_osTmpMskOverviewFilename);
     736             :         }
     737             :     }
     738         153 : }
     739             : 
     740             : /************************************************************************/
     741             : /*                    GDALCOGCreator::Create()                          */
     742             : /************************************************************************/
     743             : 
     744         153 : GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
     745             :                                     GDALDataset *const poSrcDS,
     746             :                                     char **papszOptions,
     747             :                                     GDALProgressFunc pfnProgress,
     748             :                                     void *pProgressData)
     749             : {
     750         153 :     if (pfnProgress == nullptr)
     751           0 :         pfnProgress = GDALDummyProgress;
     752             : 
     753         153 :     if (poSrcDS->GetRasterCount() == 0)
     754             :     {
     755           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     756             :                  "COG driver does not support 0-band source raster");
     757           1 :         return nullptr;
     758             :     }
     759             : 
     760             :     const CPLString osCompress = CSLFetchNameValueDef(
     761         304 :         papszOptions, "COMPRESS", gbHasLZW ? "LZW" : "NONE");
     762             : 
     763             :     const char *pszInterleave =
     764         152 :         CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "PIXEL");
     765         152 :     if (EQUAL(osCompress, "WEBP"))
     766             :     {
     767           4 :         if (!EQUAL(pszInterleave, "PIXEL"))
     768             :         {
     769           1 :             CPLError(CE_Failure, CPLE_NotSupported,
     770             :                      "COMPRESS=WEBP only supported for INTERLEAVE=PIXEL");
     771           1 :             return nullptr;
     772             :         }
     773             :     }
     774             : 
     775             :     CPLConfigOptionSetter oSetterReportDirtyBlockFlushing(
     776         302 :         "GDAL_REPORT_DIRTY_BLOCK_FLUSHING", "NO", true);
     777             : 
     778             :     const char *pszStatistics =
     779         151 :         CSLFetchNameValueDef(papszOptions, "STATISTICS", "AUTO");
     780         151 :     auto poSrcFirstBand = poSrcDS->GetRasterBand(1);
     781             :     const bool bSrcHasStatistics =
     782         151 :         poSrcFirstBand->GetMetadataItem("STATISTICS_MINIMUM") &&
     783           8 :         poSrcFirstBand->GetMetadataItem("STATISTICS_MAXIMUM") &&
     784         167 :         poSrcFirstBand->GetMetadataItem("STATISTICS_MEAN") &&
     785           8 :         poSrcFirstBand->GetMetadataItem("STATISTICS_STDDEV");
     786         151 :     bool bNeedStats = false;
     787         151 :     bool bRemoveStats = false;
     788         151 :     bool bWrkHasStatistics = bSrcHasStatistics;
     789         151 :     if (EQUAL(pszStatistics, "AUTO"))
     790             :     {
     791             :         // nothing
     792             :     }
     793          10 :     else if (CPLTestBool(pszStatistics))
     794             :     {
     795           5 :         bNeedStats = true;
     796             :     }
     797             :     else
     798             :     {
     799           5 :         bRemoveStats = true;
     800             :     }
     801             : 
     802         151 :     double dfCurPixels = 0;
     803         151 :     double dfTotalPixelsToProcess = 0;
     804         151 :     GDALDataset *poCurDS = poSrcDS;
     805             : 
     806         151 :     std::unique_ptr<gdal::TileMatrixSet> poTM;
     807         151 :     int nZoomLevel = 0;
     808         151 :     int nAlignedLevels = 0;
     809         151 :     if (COGHasWarpingOptions(papszOptions))
     810             :     {
     811          33 :         CPLString osTargetResampling;
     812          33 :         CPLString osTargetSRS;
     813          33 :         int nTargetXSize = 0;
     814          33 :         int nTargetYSize = 0;
     815          33 :         double dfTargetMinX = 0;
     816          33 :         double dfTargetMinY = 0;
     817          33 :         double dfTargetMaxX = 0;
     818          33 :         double dfTargetMaxY = 0;
     819          33 :         double dfRes = 0;
     820          33 :         if (!COGGetWarpingCharacteristics(
     821             :                 poCurDS, papszOptions, osTargetResampling, osTargetSRS,
     822             :                 nTargetXSize, nTargetYSize, dfTargetMinX, dfTargetMinY,
     823             :                 dfTargetMaxX, dfTargetMaxY, dfRes, poTM, nZoomLevel,
     824             :                 nAlignedLevels))
     825             :         {
     826           2 :             return nullptr;
     827             :         }
     828             : 
     829             :         // Collect information on source dataset to see if it already
     830             :         // matches the warping specifications
     831          31 :         CPLString osSrcSRS;
     832          31 :         const auto poSrcSRS = poCurDS->GetSpatialRef();
     833          31 :         if (poSrcSRS)
     834             :         {
     835          31 :             const char *pszAuthName = poSrcSRS->GetAuthorityName(nullptr);
     836          31 :             const char *pszAuthCode = poSrcSRS->GetAuthorityCode(nullptr);
     837          31 :             if (pszAuthName && pszAuthCode)
     838             :             {
     839          31 :                 osSrcSRS = pszAuthName;
     840          31 :                 osSrcSRS += ':';
     841          31 :                 osSrcSRS += pszAuthCode;
     842             :             }
     843             :         }
     844          31 :         double dfSrcMinX = 0;
     845          31 :         double dfSrcMinY = 0;
     846          31 :         double dfSrcMaxX = 0;
     847          31 :         double dfSrcMaxY = 0;
     848          31 :         GDALGeoTransform srcGT;
     849          31 :         const int nSrcXSize = poCurDS->GetRasterXSize();
     850          31 :         const int nSrcYSize = poCurDS->GetRasterYSize();
     851          31 :         if (poCurDS->GetGeoTransform(srcGT) == CE_None)
     852             :         {
     853          31 :             dfSrcMinX = srcGT[0];
     854          31 :             dfSrcMaxY = srcGT[3];
     855          31 :             dfSrcMaxX = srcGT[0] + nSrcXSize * srcGT[1];
     856          31 :             dfSrcMinY = srcGT[3] + nSrcYSize * srcGT[5];
     857             :         }
     858             : 
     859          24 :         if (nTargetXSize == nSrcXSize && nTargetYSize == nSrcYSize &&
     860          12 :             osTargetSRS == osSrcSRS &&
     861           6 :             fabs(dfSrcMinX - dfTargetMinX) < 1e-10 * fabs(dfSrcMinX) &&
     862           5 :             fabs(dfSrcMinY - dfTargetMinY) < 1e-10 * fabs(dfSrcMinY) &&
     863          48 :             fabs(dfSrcMaxX - dfTargetMaxX) < 1e-10 * fabs(dfSrcMaxX) &&
     864           5 :             fabs(dfSrcMaxY - dfTargetMaxY) < 1e-10 * fabs(dfSrcMaxY))
     865             :         {
     866           5 :             CPLDebug("COG",
     867             :                      "Skipping reprojection step: "
     868             :                      "source dataset matches reprojection specifications");
     869             :         }
     870             :         else
     871             :         {
     872          52 :             m_poReprojectedDS = CreateReprojectedDS(
     873             :                 pszFilename, poCurDS, papszOptions, osTargetResampling,
     874             :                 osTargetSRS, nTargetXSize, nTargetYSize, dfTargetMinX,
     875             :                 dfTargetMinY, dfTargetMaxX, dfTargetMaxY, dfRes, pfnProgress,
     876          26 :                 pProgressData, dfCurPixels, dfTotalPixelsToProcess);
     877          26 :             if (!m_poReprojectedDS)
     878           1 :                 return nullptr;
     879          25 :             poCurDS = m_poReprojectedDS.get();
     880             : 
     881          25 :             if (bSrcHasStatistics && !bNeedStats && !bRemoveStats)
     882             :             {
     883           1 :                 bNeedStats = true;
     884             :             }
     885          25 :             bWrkHasStatistics = false;
     886             :         }
     887             :     }
     888             : 
     889         160 :     if (EQUAL(osCompress, "JPEG") && EQUAL(pszInterleave, "PIXEL") &&
     890         170 :         (poCurDS->GetRasterCount() == 2 || poCurDS->GetRasterCount() == 4) &&
     891          10 :         poCurDS->GetRasterBand(poCurDS->GetRasterCount())
     892          10 :                 ->GetColorInterpretation() == GCI_AlphaBand)
     893             :     {
     894          10 :         char **papszArg = nullptr;
     895          10 :         papszArg = CSLAddString(papszArg, "-of");
     896          10 :         papszArg = CSLAddString(papszArg, "VRT");
     897          10 :         papszArg = CSLAddString(papszArg, "-b");
     898          10 :         papszArg = CSLAddString(papszArg, "1");
     899          10 :         if (poCurDS->GetRasterCount() == 2)
     900             :         {
     901           1 :             papszArg = CSLAddString(papszArg, "-mask");
     902           1 :             papszArg = CSLAddString(papszArg, "2");
     903             :         }
     904             :         else
     905             :         {
     906           9 :             CPLAssert(poCurDS->GetRasterCount() == 4);
     907           9 :             papszArg = CSLAddString(papszArg, "-b");
     908           9 :             papszArg = CSLAddString(papszArg, "2");
     909           9 :             papszArg = CSLAddString(papszArg, "-b");
     910           9 :             papszArg = CSLAddString(papszArg, "3");
     911           9 :             papszArg = CSLAddString(papszArg, "-mask");
     912           9 :             papszArg = CSLAddString(papszArg, "4");
     913             :         }
     914             :         GDALTranslateOptions *psOptions =
     915          10 :             GDALTranslateOptionsNew(papszArg, nullptr);
     916          10 :         CSLDestroy(papszArg);
     917          10 :         GDALDatasetH hRGBMaskDS = GDALTranslate(
     918             :             "", GDALDataset::ToHandle(poCurDS), psOptions, nullptr);
     919          10 :         GDALTranslateOptionsFree(psOptions);
     920          10 :         if (!hRGBMaskDS)
     921             :         {
     922           0 :             return nullptr;
     923             :         }
     924          10 :         m_poRGBMaskDS.reset(GDALDataset::FromHandle(hRGBMaskDS));
     925          10 :         poCurDS = m_poRGBMaskDS.get();
     926             : 
     927          10 :         if (bSrcHasStatistics && !bNeedStats && !bRemoveStats)
     928             :         {
     929           1 :             bNeedStats = true;
     930             :         }
     931           9 :         else if (bRemoveStats && bWrkHasStatistics)
     932             :         {
     933           1 :             poCurDS->ClearStatistics();
     934           1 :             bRemoveStats = false;
     935             :         }
     936             :     }
     937             : 
     938         148 :     const int nBands = poCurDS->GetRasterCount();
     939         148 :     const int nXSize = poCurDS->GetRasterXSize();
     940         148 :     const int nYSize = poCurDS->GetRasterYSize();
     941             : 
     942          12 :     const auto CreateVRTWithOrWithoutStats = [this, &poCurDS]()
     943             :     {
     944           3 :         const char *const apszOptions[] = {"-of", "VRT", nullptr};
     945             :         GDALTranslateOptions *psOptions =
     946           3 :             GDALTranslateOptionsNew(const_cast<char **>(apszOptions), nullptr);
     947           3 :         GDALDatasetH hVRTDS = GDALTranslate("", GDALDataset::ToHandle(poCurDS),
     948             :                                             psOptions, nullptr);
     949           3 :         GDALTranslateOptionsFree(psOptions);
     950           3 :         if (!hVRTDS)
     951           0 :             return false;
     952           3 :         m_poVRTWithOrWithoutStats.reset(GDALDataset::FromHandle(hVRTDS));
     953           3 :         poCurDS = m_poVRTWithOrWithoutStats.get();
     954           3 :         return true;
     955         148 :     };
     956             : 
     957         148 :     if (bNeedStats)
     958             :     {
     959           7 :         if (poSrcDS == poCurDS && !CreateVRTWithOrWithoutStats())
     960             :         {
     961           0 :             return nullptr;
     962             :         }
     963           7 :         poCurDS->ClearStatistics();
     964             : 
     965             :         // Avoid source files to be modified
     966             :         CPLConfigOptionSetter enablePamDirtyDisabler(
     967          14 :             "GDAL_PAM_ENABLE_MARK_DIRTY", "NO", true);
     968             : 
     969          21 :         for (int i = 1; i <= nBands; ++i)
     970             :         {
     971          14 :             poCurDS->GetRasterBand(i)->ComputeStatistics(
     972             :                 /*bApproxOK=*/FALSE, nullptr, nullptr, nullptr, nullptr,
     973          14 :                 nullptr, nullptr);
     974             :         }
     975             :     }
     976         141 :     else if (bRemoveStats && bWrkHasStatistics)
     977             :     {
     978           1 :         if (!CreateVRTWithOrWithoutStats())
     979           0 :             return nullptr;
     980             : 
     981           1 :         m_poVRTWithOrWithoutStats->ClearStatistics();
     982             :     }
     983             : 
     984         296 :     CPLString osBlockSize(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", ""));
     985         148 :     if (osBlockSize.empty())
     986             :     {
     987         122 :         if (poTM)
     988             :         {
     989          20 :             osBlockSize.Printf("%d", poTM->tileMatrixList()[0].mTileWidth);
     990             :         }
     991             :         else
     992             :         {
     993         102 :             osBlockSize = "512";
     994             :         }
     995             :     }
     996             : 
     997         148 :     const int nOvrThresholdSize = atoi(osBlockSize);
     998             : 
     999         148 :     const auto poFirstBand = poCurDS->GetRasterBand(1);
    1000         148 :     const bool bHasMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
    1001             : 
    1002             :     CPLString osOverviews =
    1003         296 :         CSLFetchNameValueDef(papszOptions, "OVERVIEWS", "AUTO");
    1004             :     const bool bUseExistingOrNone =
    1005         148 :         EQUAL(osOverviews, "FORCE_USE_EXISTING") || EQUAL(osOverviews, "NONE");
    1006             : 
    1007             :     const int nOverviewCount =
    1008         148 :         atoi(CSLFetchNameValueDef(papszOptions, "OVERVIEW_COUNT", "-1"));
    1009             : 
    1010             :     const bool bGenerateMskOvr =
    1011         144 :         !bUseExistingOrNone && bHasMask &&
    1012           7 :         (nXSize > nOvrThresholdSize || nYSize > nOvrThresholdSize ||
    1013         292 :          nOverviewCount > 0) &&
    1014          13 :         (EQUAL(osOverviews, "IGNORE_EXISTING") ||
    1015          13 :          poFirstBand->GetMaskBand()->GetOverviewCount() == 0);
    1016             :     const bool bGenerateOvr =
    1017         144 :         !bUseExistingOrNone &&
    1018         101 :         (nXSize > nOvrThresholdSize || nYSize > nOvrThresholdSize ||
    1019         292 :          nOverviewCount > 0) &&
    1020          47 :         (EQUAL(osOverviews, "IGNORE_EXISTING") ||
    1021          45 :          poFirstBand->GetOverviewCount() == 0);
    1022             : 
    1023         296 :     std::vector<std::pair<int, int>> asOverviewDims;
    1024         148 :     int nTmpXSize = nXSize;
    1025         148 :     int nTmpYSize = nYSize;
    1026         148 :     if (poTM)
    1027             :     {
    1028          23 :         const auto &tmList = poTM->tileMatrixList();
    1029          23 :         int nCurLevel = nZoomLevel;
    1030             :         while (true)
    1031             :         {
    1032          47 :             if (nOverviewCount < 0)
    1033             :             {
    1034          33 :                 if (nTmpXSize <= nOvrThresholdSize &&
    1035          21 :                     nTmpYSize <= nOvrThresholdSize)
    1036          20 :                     break;
    1037             :             }
    1038          14 :             else if (static_cast<int>(asOverviewDims.size()) ==
    1039          26 :                          nOverviewCount ||
    1040          12 :                      (nTmpXSize == 1 && nTmpYSize == 1))
    1041             :             {
    1042           3 :                 break;
    1043             :             }
    1044             :             const double dfResRatio =
    1045             :                 (nCurLevel >= 1)
    1046          24 :                     ? tmList[nCurLevel - 1].mResX / tmList[nCurLevel].mResX
    1047          24 :                     : 2;
    1048          24 :             nTmpXSize = static_cast<int>(nTmpXSize / dfResRatio + 0.5);
    1049          24 :             nTmpYSize = static_cast<int>(nTmpYSize / dfResRatio + 0.5);
    1050          24 :             if (nTmpXSize == 0)
    1051           0 :                 nTmpXSize = 1;
    1052          24 :             if (nTmpYSize == 0)
    1053           0 :                 nTmpYSize = 1;
    1054          24 :             asOverviewDims.emplace_back(std::pair(nTmpXSize, nTmpYSize));
    1055          24 :             nCurLevel--;
    1056          24 :         }
    1057             :     }
    1058         125 :     else if (bGenerateMskOvr || bGenerateOvr)
    1059             :     {
    1060          33 :         if (!bGenerateOvr)
    1061             :         {
    1062             :             // If generating only .msk.ovr, use the exact overview size as
    1063             :             // the overviews of the imagery.
    1064           1 :             int nIters = poFirstBand->GetOverviewCount();
    1065           1 :             if (nOverviewCount >= 0 && nOverviewCount < nIters)
    1066           0 :                 nIters = nOverviewCount;
    1067           2 :             for (int i = 0; i < nIters; i++)
    1068             :             {
    1069           1 :                 auto poOvrBand = poFirstBand->GetOverview(i);
    1070             :                 asOverviewDims.emplace_back(
    1071           1 :                     std::pair(poOvrBand->GetXSize(), poOvrBand->GetYSize()));
    1072             :             }
    1073             :         }
    1074             :         else
    1075             :         {
    1076             :             while (true)
    1077             :             {
    1078          92 :                 if (nOverviewCount < 0)
    1079             :                 {
    1080          74 :                     if (nTmpXSize <= nOvrThresholdSize &&
    1081          28 :                         nTmpYSize <= nOvrThresholdSize)
    1082          28 :                         break;
    1083             :                 }
    1084          18 :                 else if (static_cast<int>(asOverviewDims.size()) ==
    1085          33 :                              nOverviewCount ||
    1086          15 :                          (nTmpXSize == 1 && nTmpYSize == 1))
    1087             :                 {
    1088           4 :                     break;
    1089             :                 }
    1090          60 :                 nTmpXSize /= 2;
    1091          60 :                 nTmpYSize /= 2;
    1092          60 :                 if (nTmpXSize == 0)
    1093           0 :                     nTmpXSize = 1;
    1094          60 :                 if (nTmpYSize == 0)
    1095           1 :                     nTmpYSize = 1;
    1096          60 :                 asOverviewDims.emplace_back(std::pair(nTmpXSize, nTmpYSize));
    1097             :             }
    1098             :         }
    1099             :     }
    1100             : 
    1101         148 :     if (dfTotalPixelsToProcess == 0.0)
    1102             :     {
    1103         123 :         dfTotalPixelsToProcess =
    1104         123 :             (bGenerateMskOvr ? double(nXSize) * nYSize / 3 : 0) +
    1105         123 :             (bGenerateOvr ? double(nXSize) * nYSize * nBands / 3 : 0) +
    1106         123 :             double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) * 4. / 3;
    1107             :     }
    1108             : 
    1109         296 :     CPLStringList aosOverviewOptions;
    1110             :     aosOverviewOptions.SetNameValue(
    1111             :         "COMPRESS",
    1112             :         CPLGetConfigOption("COG_TMP_COMPRESSION",  // only for debug purposes
    1113         148 :                            HasZSTDCompression() ? "ZSTD" : "LZW"));
    1114             :     aosOverviewOptions.SetNameValue(
    1115         148 :         "NUM_THREADS", CSLFetchNameValue(papszOptions, "NUM_THREADS"));
    1116         148 :     aosOverviewOptions.SetNameValue("BIGTIFF", "YES");
    1117         148 :     aosOverviewOptions.SetNameValue("SPARSE_OK", "YES");
    1118             : 
    1119         148 :     if (bGenerateMskOvr)
    1120             :     {
    1121          10 :         CPLDebug("COG", "Generating overviews of the mask: start");
    1122          10 :         m_osTmpMskOverviewFilename = GetTmpFilename(pszFilename, "msk.ovr.tmp");
    1123          10 :         GDALRasterBand *poSrcMask = poFirstBand->GetMaskBand();
    1124          10 :         const char *pszResampling = CSLFetchNameValueDef(
    1125             :             papszOptions, "OVERVIEW_RESAMPLING",
    1126             :             CSLFetchNameValueDef(papszOptions, "RESAMPLING",
    1127             :                                  GetResampling(poSrcDS)));
    1128             : 
    1129          10 :         double dfNextPixels = dfCurPixels + double(nXSize) * nYSize / 3;
    1130          10 :         void *pScaledProgress = GDALCreateScaledProgress(
    1131             :             dfCurPixels / dfTotalPixelsToProcess,
    1132             :             dfNextPixels / dfTotalPixelsToProcess, pfnProgress, pProgressData);
    1133          10 :         dfCurPixels = dfNextPixels;
    1134             : 
    1135          10 :         CPLErr eErr = GTIFFBuildOverviewsEx(
    1136             :             m_osTmpMskOverviewFilename, 1, &poSrcMask,
    1137          10 :             static_cast<int>(asOverviewDims.size()), nullptr,
    1138          10 :             asOverviewDims.data(), pszResampling, aosOverviewOptions.List(),
    1139             :             GDALScaledProgress, pScaledProgress);
    1140          10 :         CPLDebug("COG", "Generating overviews of the mask: end");
    1141             : 
    1142          10 :         GDALDestroyScaledProgress(pScaledProgress);
    1143          10 :         if (eErr != CE_None)
    1144             :         {
    1145           0 :             return nullptr;
    1146             :         }
    1147             :     }
    1148             : 
    1149         148 :     if (bGenerateOvr)
    1150             :     {
    1151          39 :         CPLDebug("COG", "Generating overviews of the imagery: start");
    1152          39 :         m_osTmpOverviewFilename = GetTmpFilename(pszFilename, "ovr.tmp");
    1153          39 :         std::vector<GDALRasterBand *> apoSrcBands;
    1154         121 :         for (int i = 0; i < nBands; i++)
    1155          82 :             apoSrcBands.push_back(poCurDS->GetRasterBand(i + 1));
    1156          39 :         const char *pszResampling = CSLFetchNameValueDef(
    1157             :             papszOptions, "OVERVIEW_RESAMPLING",
    1158             :             CSLFetchNameValueDef(papszOptions, "RESAMPLING",
    1159             :                                  GetResampling(poSrcDS)));
    1160             : 
    1161          39 :         double dfNextPixels =
    1162          39 :             dfCurPixels + double(nXSize) * nYSize * nBands / 3;
    1163          39 :         void *pScaledProgress = GDALCreateScaledProgress(
    1164             :             dfCurPixels / dfTotalPixelsToProcess,
    1165             :             dfNextPixels / dfTotalPixelsToProcess, pfnProgress, pProgressData);
    1166          39 :         dfCurPixels = dfNextPixels;
    1167             : 
    1168          39 :         if (nBands > 1)
    1169             :         {
    1170          22 :             aosOverviewOptions.SetNameValue("INTERLEAVE", "PIXEL");
    1171             :         }
    1172          39 :         if (!m_osTmpMskOverviewFilename.empty())
    1173             :         {
    1174             :             aosOverviewOptions.SetNameValue("MASK_OVERVIEW_DATASET",
    1175           9 :                                             m_osTmpMskOverviewFilename);
    1176             :         }
    1177          39 :         CPLErr eErr = GTIFFBuildOverviewsEx(
    1178          39 :             m_osTmpOverviewFilename, nBands, &apoSrcBands[0],
    1179          39 :             static_cast<int>(asOverviewDims.size()), nullptr,
    1180          39 :             asOverviewDims.data(), pszResampling, aosOverviewOptions.List(),
    1181             :             GDALScaledProgress, pScaledProgress);
    1182          39 :         CPLDebug("COG", "Generating overviews of the imagery: end");
    1183             : 
    1184          39 :         GDALDestroyScaledProgress(pScaledProgress);
    1185          39 :         if (eErr != CE_None)
    1186             :         {
    1187           1 :             return nullptr;
    1188             :         }
    1189             :     }
    1190             : 
    1191         294 :     CPLStringList aosOptions;
    1192         147 :     aosOptions.SetNameValue("COPY_SRC_OVERVIEWS", "YES");
    1193         147 :     aosOptions.SetNameValue("COMPRESS", osCompress);
    1194         147 :     aosOptions.SetNameValue("TILED", "YES");
    1195         147 :     aosOptions.SetNameValue("BLOCKXSIZE", osBlockSize);
    1196         147 :     aosOptions.SetNameValue("BLOCKYSIZE", osBlockSize);
    1197             :     const char *pszPredictor =
    1198         147 :         CSLFetchNameValueDef(papszOptions, "PREDICTOR", "FALSE");
    1199         147 :     const char *pszPredictorValue = GetPredictor(poSrcDS, pszPredictor);
    1200         147 :     if (pszPredictorValue != nullptr)
    1201             :     {
    1202           2 :         aosOptions.SetNameValue("PREDICTOR", pszPredictorValue);
    1203             :     }
    1204             : 
    1205         147 :     const char *pszQuality = CSLFetchNameValue(papszOptions, "QUALITY");
    1206         147 :     if (EQUAL(osCompress, "JPEG"))
    1207             :     {
    1208          12 :         aosOptions.SetNameValue("JPEG_QUALITY", pszQuality);
    1209          12 :         if (nBands == 3 && EQUAL(pszInterleave, "PIXEL"))
    1210          10 :             aosOptions.SetNameValue("PHOTOMETRIC", "YCBCR");
    1211             :     }
    1212         135 :     else if (EQUAL(osCompress, "WEBP"))
    1213             :     {
    1214           3 :         if (pszQuality && atoi(pszQuality) == 100)
    1215           2 :             aosOptions.SetNameValue("WEBP_LOSSLESS", "YES");
    1216           3 :         aosOptions.SetNameValue("WEBP_LEVEL", pszQuality);
    1217             :     }
    1218         132 :     else if (EQUAL(osCompress, "DEFLATE") || EQUAL(osCompress, "LERC_DEFLATE"))
    1219             :     {
    1220             :         aosOptions.SetNameValue("ZLEVEL",
    1221           8 :                                 CSLFetchNameValue(papszOptions, "LEVEL"));
    1222             :     }
    1223         124 :     else if (EQUAL(osCompress, "ZSTD") || EQUAL(osCompress, "LERC_ZSTD"))
    1224             :     {
    1225             :         aosOptions.SetNameValue("ZSTD_LEVEL",
    1226           3 :                                 CSLFetchNameValue(papszOptions, "LEVEL"));
    1227             :     }
    1228         121 :     else if (EQUAL(osCompress, "LZMA"))
    1229             :     {
    1230             :         aosOptions.SetNameValue("LZMA_PRESET",
    1231           1 :                                 CSLFetchNameValue(papszOptions, "LEVEL"));
    1232             :     }
    1233             : 
    1234         147 :     if (STARTS_WITH_CI(osCompress, "LERC"))
    1235             :     {
    1236             :         aosOptions.SetNameValue("MAX_Z_ERROR",
    1237           8 :                                 CSLFetchNameValue(papszOptions, "MAX_Z_ERROR"));
    1238             :         aosOptions.SetNameValue(
    1239             :             "MAX_Z_ERROR_OVERVIEW",
    1240           8 :             CSLFetchNameValue(papszOptions, "MAX_Z_ERROR_OVERVIEW"));
    1241             :     }
    1242             : 
    1243         147 :     if (STARTS_WITH_CI(osCompress, "JXL"))
    1244             :     {
    1245           8 :         for (const char *pszKey : {"JXL_LOSSLESS", "JXL_EFFORT", "JXL_DISTANCE",
    1246          10 :                                    "JXL_ALPHA_DISTANCE"})
    1247             :         {
    1248           8 :             const char *pszValue = CSLFetchNameValue(papszOptions, pszKey);
    1249           8 :             if (pszValue)
    1250           3 :                 aosOptions.SetNameValue(pszKey, pszValue);
    1251             :         }
    1252             :     }
    1253             : 
    1254             :     aosOptions.SetNameValue("BIGTIFF",
    1255         147 :                             CSLFetchNameValue(papszOptions, "BIGTIFF"));
    1256             :     aosOptions.SetNameValue("NUM_THREADS",
    1257         147 :                             CSLFetchNameValue(papszOptions, "NUM_THREADS"));
    1258             :     aosOptions.SetNameValue("GEOTIFF_VERSION",
    1259         147 :                             CSLFetchNameValue(papszOptions, "GEOTIFF_VERSION"));
    1260             :     aosOptions.SetNameValue("SPARSE_OK",
    1261         147 :                             CSLFetchNameValue(papszOptions, "SPARSE_OK"));
    1262         147 :     aosOptions.SetNameValue("NBITS", CSLFetchNameValue(papszOptions, "NBITS"));
    1263             : 
    1264         147 :     if (EQUAL(osOverviews, "NONE"))
    1265             :     {
    1266           2 :         aosOptions.SetNameValue("@OVERVIEW_DATASET", "");
    1267             :     }
    1268             :     else
    1269             :     {
    1270         145 :         if (!m_osTmpOverviewFilename.empty())
    1271             :         {
    1272             :             aosOptions.SetNameValue("@OVERVIEW_DATASET",
    1273          38 :                                     m_osTmpOverviewFilename);
    1274             :         }
    1275         145 :         if (!m_osTmpMskOverviewFilename.empty())
    1276             :         {
    1277             :             aosOptions.SetNameValue("@MASK_OVERVIEW_DATASET",
    1278          10 :                                     m_osTmpMskOverviewFilename);
    1279             :         }
    1280             :         aosOptions.SetNameValue(
    1281             :             "@OVERVIEW_COUNT",
    1282         145 :             CSLFetchNameValue(papszOptions, "OVERVIEW_COUNT"));
    1283             :     }
    1284             : 
    1285             :     const CPLString osTilingScheme(
    1286         294 :         CSLFetchNameValueDef(papszOptions, "TILING_SCHEME", "CUSTOM"));
    1287         147 :     if (osTilingScheme != "CUSTOM")
    1288             :     {
    1289          23 :         aosOptions.SetNameValue("@TILING_SCHEME_NAME", osTilingScheme);
    1290             :         aosOptions.SetNameValue("@TILING_SCHEME_ZOOM_LEVEL",
    1291          23 :                                 CPLSPrintf("%d", nZoomLevel));
    1292          23 :         if (nAlignedLevels > 0)
    1293             :         {
    1294             :             aosOptions.SetNameValue("@TILING_SCHEME_ALIGNED_LEVELS",
    1295           4 :                                     CPLSPrintf("%d", nAlignedLevels));
    1296             :         }
    1297             :     }
    1298         147 :     const char *pszOverviewCompress = CSLFetchNameValueDef(
    1299             :         papszOptions, "OVERVIEW_COMPRESS", osCompress.c_str());
    1300             : 
    1301             :     CPLConfigOptionSetter ovrCompressSetter("COMPRESS_OVERVIEW",
    1302         294 :                                             pszOverviewCompress, true);
    1303             :     const char *pszOverviewQuality =
    1304         147 :         CSLFetchNameValue(papszOptions, "OVERVIEW_QUALITY");
    1305             :     CPLConfigOptionSetter ovrQualityJpegSetter("JPEG_QUALITY_OVERVIEW",
    1306         294 :                                                pszOverviewQuality, true);
    1307             : 
    1308         147 :     std::unique_ptr<CPLConfigOptionSetter> poWebpLosslessSetter;
    1309         147 :     std::unique_ptr<CPLConfigOptionSetter> poWebpLevelSetter;
    1310         147 :     if (EQUAL(pszOverviewCompress, "WEBP"))
    1311             :     {
    1312           3 :         if (pszOverviewQuality && CPLAtof(pszOverviewQuality) == 100.0)
    1313             :         {
    1314           0 :             poWebpLosslessSetter.reset(new CPLConfigOptionSetter(
    1315           0 :                 "WEBP_LOSSLESS_OVERVIEW", "TRUE", true));
    1316             :         }
    1317             :         else
    1318             :         {
    1319           3 :             poWebpLosslessSetter.reset(new CPLConfigOptionSetter(
    1320           3 :                 "WEBP_LOSSLESS_OVERVIEW", "FALSE", true));
    1321           3 :             poWebpLevelSetter.reset(new CPLConfigOptionSetter(
    1322           3 :                 "WEBP_LEVEL_OVERVIEW", pszOverviewQuality, true));
    1323             :         }
    1324             :     }
    1325             : 
    1326         147 :     std::unique_ptr<CPLConfigOptionSetter> poPhotometricSetter;
    1327         147 :     if (nBands == 3 && EQUAL(pszOverviewCompress, "JPEG") &&
    1328          11 :         EQUAL(pszInterleave, "PIXEL"))
    1329             :     {
    1330          10 :         poPhotometricSetter.reset(
    1331          10 :             new CPLConfigOptionSetter("PHOTOMETRIC_OVERVIEW", "YCBCR", true));
    1332             :     }
    1333             : 
    1334             :     const char *osOvrPredictor =
    1335         147 :         CSLFetchNameValueDef(papszOptions, "OVERVIEW_PREDICTOR", "FALSE");
    1336         147 :     const char *pszOvrPredictorValue = GetPredictor(poSrcDS, osOvrPredictor);
    1337             :     CPLConfigOptionSetter ovrPredictorSetter("PREDICTOR_OVERVIEW",
    1338         294 :                                              pszOvrPredictorValue, true);
    1339             : 
    1340             :     GDALDriver *poGTiffDrv =
    1341         147 :         GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
    1342         147 :     if (!poGTiffDrv)
    1343           0 :         return nullptr;
    1344         147 :     void *pScaledProgress = GDALCreateScaledProgress(
    1345             :         dfCurPixels / dfTotalPixelsToProcess, 1.0, pfnProgress, pProgressData);
    1346             : 
    1347             :     CPLConfigOptionSetter oSetterInternalMask("GDAL_TIFF_INTERNAL_MASK", "YES",
    1348         147 :                                               false);
    1349             : 
    1350         147 :     const char *pszCopySrcMDD = CSLFetchNameValue(papszOptions, "COPY_SRC_MDD");
    1351         147 :     if (pszCopySrcMDD)
    1352           2 :         aosOptions.SetNameValue("COPY_SRC_MDD", pszCopySrcMDD);
    1353         147 :     char **papszSrcMDD = CSLFetchNameValueMultiple(papszOptions, "SRC_MDD");
    1354         150 :     for (CSLConstList papszSrcMDDIter = papszSrcMDD;
    1355         150 :          papszSrcMDDIter && *papszSrcMDDIter; ++papszSrcMDDIter)
    1356           3 :         aosOptions.AddNameValue("SRC_MDD", *papszSrcMDDIter);
    1357         147 :     CSLDestroy(papszSrcMDD);
    1358             : 
    1359         147 :     if (EQUAL(pszInterleave, "TILE"))
    1360             :     {
    1361           5 :         aosOptions.SetNameValue("INTERLEAVE", "BAND");
    1362           5 :         aosOptions.SetNameValue("@TILE_INTERLEAVE", "YES");
    1363             :     }
    1364             :     else
    1365             :     {
    1366         142 :         aosOptions.SetNameValue("INTERLEAVE", pszInterleave);
    1367             :     }
    1368             : 
    1369         147 :     CPLDebug("COG", "Generating final product: start");
    1370             :     auto poRet =
    1371         147 :         poGTiffDrv->CreateCopy(pszFilename, poCurDS, false, aosOptions.List(),
    1372             :                                GDALScaledProgress, pScaledProgress);
    1373             : 
    1374         147 :     GDALDestroyScaledProgress(pScaledProgress);
    1375             : 
    1376         147 :     if (poRet)
    1377         135 :         poRet->FlushCache(false);
    1378             : 
    1379         147 :     CPLDebug("COG", "Generating final product: end");
    1380         147 :     return poRet;
    1381             : }
    1382             : 
    1383             : /************************************************************************/
    1384             : /*                            COGCreateCopy()                           */
    1385             : /************************************************************************/
    1386             : 
    1387         153 : static GDALDataset *COGCreateCopy(const char *pszFilename, GDALDataset *poSrcDS,
    1388             :                                   int /*bStrict*/, char **papszOptions,
    1389             :                                   GDALProgressFunc pfnProgress,
    1390             :                                   void *pProgressData)
    1391             : {
    1392         306 :     return GDALCOGCreator().Create(pszFilename, poSrcDS, papszOptions,
    1393         306 :                                    pfnProgress, pProgressData);
    1394             : }
    1395             : 
    1396             : /************************************************************************/
    1397             : /*                          GDALRegister_COG()                          */
    1398             : /************************************************************************/
    1399             : 
    1400             : class GDALCOGDriver final : public GDALDriver
    1401             : {
    1402             :     std::mutex m_oMutex{};
    1403             :     bool m_bInitialized = false;
    1404             : 
    1405             :     bool bHasLZW = false;
    1406             :     bool bHasDEFLATE = false;
    1407             :     bool bHasLZMA = false;
    1408             :     bool bHasZSTD = false;
    1409             :     bool bHasJPEG = false;
    1410             :     bool bHasWebP = false;
    1411             :     bool bHasLERC = false;
    1412             :     std::string osCompressValues{};
    1413             : 
    1414             :     void InitializeCreationOptionList();
    1415             : 
    1416             :   public:
    1417             :     GDALCOGDriver();
    1418             : 
    1419             :     const char *GetMetadataItem(const char *pszName,
    1420             :                                 const char *pszDomain) override;
    1421             : 
    1422         467 :     char **GetMetadata(const char *pszDomain) override
    1423             :     {
    1424         934 :         std::lock_guard oLock(m_oMutex);
    1425         467 :         InitializeCreationOptionList();
    1426         934 :         return GDALDriver::GetMetadata(pszDomain);
    1427             :     }
    1428             : };
    1429             : 
    1430        1653 : GDALCOGDriver::GDALCOGDriver()
    1431             : {
    1432             :     // We could defer this in InitializeCreationOptionList() but with currently
    1433             :     // released libtiff versions where there was a bug (now fixed) in
    1434             :     // TIFFGetConfiguredCODECs(), this wouldn't work properly if the LERC codec
    1435             :     // had been registered in between
    1436        1653 :     osCompressValues = GTiffGetCompressValues(bHasLZW, bHasDEFLATE, bHasLZMA,
    1437        1653 :                                               bHasZSTD, bHasJPEG, bHasWebP,
    1438        1653 :                                               bHasLERC, true /* bForCOG */);
    1439        1653 :     gbHasLZW = bHasLZW;
    1440        1653 : }
    1441             : 
    1442       43993 : const char *GDALCOGDriver::GetMetadataItem(const char *pszName,
    1443             :                                            const char *pszDomain)
    1444             : {
    1445       87986 :     std::lock_guard oLock(m_oMutex);
    1446       43993 :     if (EQUAL(pszName, GDAL_DMD_CREATIONOPTIONLIST))
    1447             :     {
    1448         384 :         InitializeCreationOptionList();
    1449             :     }
    1450       87986 :     return GDALDriver::GetMetadataItem(pszName, pszDomain);
    1451             : }
    1452             : 
    1453         851 : void GDALCOGDriver::InitializeCreationOptionList()
    1454             : {
    1455         851 :     if (m_bInitialized)
    1456         639 :         return;
    1457         212 :     m_bInitialized = true;
    1458             : 
    1459         424 :     CPLString osOptions;
    1460             :     osOptions = "<CreationOptionList>"
    1461         212 :                 "   <Option name='COMPRESS' type='string-select' default='";
    1462         212 :     osOptions += bHasLZW ? "LZW" : "NONE";
    1463         212 :     osOptions += "'>";
    1464         212 :     osOptions += osCompressValues;
    1465         212 :     osOptions += "   </Option>";
    1466             : 
    1467             :     osOptions +=
    1468         212 :         "   <Option name='OVERVIEW_COMPRESS' type='string-select' default='";
    1469         212 :     osOptions += bHasLZW ? "LZW" : "NONE";
    1470         212 :     osOptions += "'>";
    1471         212 :     osOptions += osCompressValues;
    1472         212 :     osOptions += "   </Option>";
    1473             : 
    1474         212 :     if (bHasLZW || bHasDEFLATE || bHasZSTD || bHasLZMA)
    1475             :     {
    1476         212 :         const char *osPredictorOptions =
    1477             :             "     <Value>YES</Value>"
    1478             :             "     <Value>NO</Value>"
    1479             :             "     <Value alias='2'>STANDARD</Value>"
    1480             :             "     <Value alias='3'>FLOATING_POINT</Value>";
    1481             : 
    1482             :         osOptions +=
    1483             :             "   <Option name='LEVEL' type='int' "
    1484         212 :             "description='DEFLATE/ZSTD/LZMA compression level: 1 (fastest)'/>";
    1485             : 
    1486             :         osOptions +=
    1487         212 :             "   <Option name='PREDICTOR' type='string-select' default='FALSE'>";
    1488         212 :         osOptions += osPredictorOptions;
    1489             :         osOptions += "   </Option>"
    1490             :                      "   <Option name='OVERVIEW_PREDICTOR' "
    1491         212 :                      "type='string-select' default='FALSE'>";
    1492         212 :         osOptions += osPredictorOptions;
    1493         212 :         osOptions += "   </Option>";
    1494             :     }
    1495         212 :     if (bHasJPEG || bHasWebP)
    1496             :     {
    1497         212 :         std::string osJPEG_WEBP;
    1498         212 :         if (bHasJPEG)
    1499         212 :             osJPEG_WEBP = "JPEG";
    1500         212 :         if (bHasWebP)
    1501             :         {
    1502         212 :             if (!osJPEG_WEBP.empty())
    1503         212 :                 osJPEG_WEBP += '/';
    1504         212 :             osJPEG_WEBP += "WEBP";
    1505             :         }
    1506             :         osOptions += "   <Option name='QUALITY' type='int' "
    1507         424 :                      "description='" +
    1508         424 :                      osJPEG_WEBP +
    1509             :                      " quality 1-100' min='1' max='100' default='75'/>"
    1510             :                      "   <Option name='OVERVIEW_QUALITY' type='int' "
    1511         424 :                      "description='Overview " +
    1512         424 :                      osJPEG_WEBP +
    1513             :                      " quality 1-100' min='1' max='100' "
    1514         212 :                      "default='75'/>";
    1515             :     }
    1516         212 :     if (bHasLERC)
    1517             :     {
    1518             :         osOptions +=
    1519             :             ""
    1520             :             "   <Option name='MAX_Z_ERROR' type='float' description='Maximum "
    1521             :             "error for LERC compression' default='0'/>"
    1522             :             "   <Option name='MAX_Z_ERROR_OVERVIEW' type='float' "
    1523             :             "description='Maximum error for LERC compression in overviews' "
    1524         212 :             "default='0'/>";
    1525             :     }
    1526             : #ifdef HAVE_JXL
    1527             :     osOptions +=
    1528             :         ""
    1529             :         "   <Option name='JXL_LOSSLESS' type='boolean' description='Whether "
    1530             :         "JPEGXL compression should be lossless' default='YES'/>"
    1531             :         "   <Option name='JXL_EFFORT' type='int' description='Level of effort "
    1532             :         "1(fast)-9(slow)' min='1' max='9' default='5'/>"
    1533             :         "   <Option name='JXL_DISTANCE' type='float' description='Distance "
    1534             :         "level for lossy compression (0=mathematically lossless, 1.0=visually "
    1535         212 :         "lossless, usual range [0.5,3])' default='1.0' min='0.01' max='25.0'/>";
    1536             : #ifdef HAVE_JxlEncoderSetExtraChannelDistance
    1537             :     osOptions += "   <Option name='JXL_ALPHA_DISTANCE' type='float' "
    1538             :                  "description='Distance level for alpha channel "
    1539             :                  "(-1=same as non-alpha channels, "
    1540             :                  "0=mathematically lossless, 1.0=visually lossless, "
    1541         212 :                  "usual range [0.5,3])' default='-1' min='-1' max='25.0'/>";
    1542             : #endif
    1543             : #endif
    1544             :     osOptions +=
    1545             :         "   <Option name='NUM_THREADS' type='string' "
    1546             :         "description='Number of worker threads for compression. "
    1547             :         "Can be set to ALL_CPUS' default='1'/>"
    1548             :         "   <Option name='NBITS' type='int' description='BITS for sub-byte "
    1549             :         "files (1-7), sub-uint16_t (9-15), sub-uint32_t (17-31), or float32 "
    1550             :         "(16)'/>"
    1551             :         "   <Option name='BLOCKSIZE' type='int' "
    1552             :         "description='Tile size in pixels' min='128' default='512'/>"
    1553             :         "   <Option name='INTERLEAVE' type='string-select' default='PIXEL'>"
    1554             :         "       <Value>BAND</Value>"
    1555             :         "       <Value>PIXEL</Value>"
    1556             :         "       <Value>TILE</Value>"
    1557             :         "   </Option>"
    1558             :         "   <Option name='BIGTIFF' type='string-select' description='"
    1559             :         "Force creation of BigTIFF file'>"
    1560             :         "     <Value>YES</Value>"
    1561             :         "     <Value>NO</Value>"
    1562             :         "     <Value>IF_NEEDED</Value>"
    1563             :         "     <Value>IF_SAFER</Value>"
    1564             :         "   </Option>"
    1565             :         "   <Option name='RESAMPLING' type='string' "
    1566             :         "description='Resampling method for overviews or warping'/>"
    1567             :         "   <Option name='OVERVIEW_RESAMPLING' type='string' "
    1568             :         "description='Resampling method for overviews'/>"
    1569             :         "   <Option name='WARP_RESAMPLING' type='string' "
    1570             :         "description='Resampling method for warping'/>"
    1571             :         "   <Option name='OVERVIEWS' type='string-select' description='"
    1572             :         "Behavior regarding overviews'>"
    1573             :         "     <Value>AUTO</Value>"
    1574             :         "     <Value>IGNORE_EXISTING</Value>"
    1575             :         "     <Value>FORCE_USE_EXISTING</Value>"
    1576             :         "     <Value>NONE</Value>"
    1577             :         "   </Option>"
    1578             :         "  <Option name='OVERVIEW_COUNT' type='int' min='0' "
    1579             :         "description='Number of overviews'/>"
    1580             :         "  <Option name='TILING_SCHEME' type='string-select' description='"
    1581             :         "Which tiling scheme to use pre-defined value or custom inline/outline "
    1582             :         "JSON definition' default='CUSTOM'>"
    1583         212 :         "    <Value>CUSTOM</Value>";
    1584             : 
    1585         424 :     const auto tmsList = gdal::TileMatrixSet::listPredefinedTileMatrixSets();
    1586        2120 :     for (const auto &tmsName : tmsList)
    1587             :     {
    1588        3816 :         const auto poTM = gdal::TileMatrixSet::parse(tmsName.c_str());
    1589        3816 :         if (poTM && poTM->haveAllLevelsSameTopLeft() &&
    1590        5724 :             poTM->haveAllLevelsSameTileSize() &&
    1591        1908 :             !poTM->hasVariableMatrixWidth())
    1592             :         {
    1593        1908 :             osOptions += "    <Value>";
    1594        1908 :             osOptions += tmsName;
    1595        1908 :             osOptions += "</Value>";
    1596             :         }
    1597             :     }
    1598             : 
    1599             :     osOptions +=
    1600             :         "  </Option>"
    1601             :         "  <Option name='ZOOM_LEVEL' type='int' description='Target zoom "
    1602             :         "level. "
    1603             :         "Only used for TILING_SCHEME != CUSTOM'/>"
    1604             :         "  <Option name='ZOOM_LEVEL_STRATEGY' type='string-select' "
    1605             :         "description='Strategy to determine zoom level. "
    1606             :         "Only used for TILING_SCHEME != CUSTOM' default='AUTO'>"
    1607             :         "    <Value>AUTO</Value>"
    1608             :         "    <Value>LOWER</Value>"
    1609             :         "    <Value>UPPER</Value>"
    1610             :         "  </Option>"
    1611             :         "   <Option name='TARGET_SRS' type='string' "
    1612             :         "description='Target SRS as EPSG:XXXX, WKT or PROJ string for "
    1613             :         "reprojection'/>"
    1614             :         "  <Option name='RES' type='float' description='"
    1615             :         "Target resolution for reprojection'/>"
    1616             :         "  <Option name='EXTENT' type='string' description='"
    1617             :         "Target extent as minx,miny,maxx,maxy for reprojection'/>"
    1618             :         "  <Option name='ALIGNED_LEVELS' type='int' description='"
    1619             :         "Number of resolution levels for which the tiles from GeoTIFF and the "
    1620             :         "specified tiling scheme match'/>"
    1621             :         "  <Option name='ADD_ALPHA' type='boolean' description='Can be set to "
    1622             :         "NO to "
    1623             :         "disable the addition of an alpha band in case of reprojection' "
    1624             :         "default='YES'/>"
    1625             : #if LIBGEOTIFF_VERSION >= 1600
    1626             :         "   <Option name='GEOTIFF_VERSION' type='string-select' default='AUTO' "
    1627             :         "description='Which version of GeoTIFF must be used'>"
    1628             :         "       <Value>AUTO</Value>"
    1629             :         "       <Value>1.0</Value>"
    1630             :         "       <Value>1.1</Value>"
    1631             :         "   </Option>"
    1632             : #endif
    1633             :         "   <Option name='SPARSE_OK' type='boolean' description='Should empty "
    1634             :         "blocks be omitted on disk?' default='FALSE'/>"
    1635             :         "   <Option name='STATISTICS' type='string-select' default='AUTO' "
    1636             :         "description='Which to add statistics to the output file'>"
    1637             :         "       <Value>AUTO</Value>"
    1638             :         "       <Value>YES</Value>"
    1639             :         "       <Value>NO</Value>"
    1640             :         "   </Option>"
    1641         212 :         "</CreationOptionList>";
    1642             : 
    1643         212 :     SetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST, osOptions.c_str());
    1644             : }
    1645             : 
    1646        1935 : void GDALRegister_COG()
    1647             : 
    1648             : {
    1649        1935 :     if (GDALGetDriverByName("COG") != nullptr)
    1650         282 :         return;
    1651             : 
    1652        1653 :     auto poDriver = new GDALCOGDriver();
    1653        1653 :     poDriver->SetDescription("COG");
    1654        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1655        1653 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
    1656        1653 :                               "Cloud optimized GeoTIFF generator");
    1657        1653 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/cog.html");
    1658        1653 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "tif tiff");
    1659             : 
    1660        1653 :     poDriver->SetMetadataItem(
    1661             :         GDAL_DMD_CREATIONDATATYPES,
    1662             :         "Byte Int8 UInt16 Int16 UInt32 Int32 UInt64 Int64 Float32 "
    1663        1653 :         "Float64 CInt16 CInt32 CFloat32 CFloat64");
    1664             : 
    1665        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1666             : 
    1667        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_COORDINATE_EPOCH, "YES");
    1668             : 
    1669        1653 :     poDriver->pfnCreateCopy = COGCreateCopy;
    1670             : 
    1671        1653 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1672             : }

Generated by: LCOV version 1.14