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  * Permission is hereby granted, free of charge, to any person obtaining a
11  * copy of this software and associated documentation files (the "Software"),
12  * to deal in the Software without restriction, including without limitation
13  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  * and/or sell copies of the Software, and to permit persons to whom the
15  * Software is furnished to do so, subject to the following conditions:
16  *
17  * The above copyright notice and this permission notice shall be included
18  * in all copies or substantial portions of the Software.
19  *
20  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  * DEALINGS IN THE SOFTWARE.
27  ****************************************************************************/
28 
29 #include "gdal_priv.h"
30 #include "gtiff.h"
31 #include "gt_overview.h"
32 #include "gdal_utils.h"
33 #include "gdalwarper.h"
34 #include "cogdriver.h"
35 #include "geotiff.h"
36 
37 #include "tilematrixset.hpp"
38 
39 #include <algorithm>
40 #include <memory>
41 #include <vector>
42 
43 extern "C" CPL_DLL void GDALRegister_COG();
44 
45 /************************************************************************/
46 /*                        HasZSTDCompression()                          */
47 /************************************************************************/
48 
HasZSTDCompression()49 static bool HasZSTDCompression()
50 {
51     TIFFCodec *codecs = TIFFGetConfiguredCODECs();
52     bool bHasZSTD = false;
53     for( TIFFCodec *c = codecs; c->name; ++c )
54     {
55         if( c->scheme == COMPRESSION_ZSTD )
56         {
57             bHasZSTD = true;
58             break;
59         }
60     }
61     _TIFFfree( codecs );
62     return bHasZSTD;
63 }
64 
65 /************************************************************************/
66 /*                           GetTmpFilename()                           */
67 /************************************************************************/
68 
GetTmpFilename(const char * pszFilename,const char * pszExt)69 static CPLString GetTmpFilename(const char* pszFilename,
70                                 const char* pszExt)
71 {
72     CPLString osTmpFilename;
73     osTmpFilename.Printf("%s.%s", pszFilename, pszExt);
74     VSIUnlink(osTmpFilename);
75     return osTmpFilename;
76 }
77 
78 /************************************************************************/
79 /*                             GetResampling()                          */
80 /************************************************************************/
81 
GetResampling(GDALDataset * poSrcDS)82 static const char* GetResampling(GDALDataset* poSrcDS)
83 {
84     return poSrcDS->GetRasterBand(1)->GetColorTable() ? "NEAREST" : "CUBIC";
85 }
86 
87 /************************************************************************/
88 /*                             GetPredictor()                          */
89 /************************************************************************/
GetPredictor(GDALDataset * poSrcDS,const char * pszPredictor)90 static const char* GetPredictor(GDALDataset* poSrcDS,
91                                 const char* pszPredictor)
92 {
93     if (pszPredictor == nullptr) return nullptr;
94 
95     if( EQUAL(pszPredictor, "YES") || EQUAL(pszPredictor, "ON") || EQUAL(pszPredictor, "TRUE") )
96     {
97         if( GDALDataTypeIsFloating(poSrcDS->GetRasterBand(1)->GetRasterDataType()) )
98             return "3";
99         else
100             return "2";
101     }
102     else if( EQUAL(pszPredictor, "STANDARD") || EQUAL(pszPredictor, "2") )
103     {
104         return "2";
105     }
106     else if( EQUAL(pszPredictor, "FLOATING_POINT") || EQUAL(pszPredictor, "3") )
107     {
108         return "3";
109     }
110     return nullptr;
111 }
112 
113 /************************************************************************/
114 /*                     COGGetWarpingCharacteristics()                   */
115 /************************************************************************/
116 
117 static
COGGetWarpingCharacteristics(GDALDataset * poSrcDS,const char * const * papszOptions,CPLString & osResampling,CPLString & osTargetSRS,int & nXSize,int & nYSize,double & dfMinX,double & dfMinY,double & dfMaxX,double & dfMaxY,std::unique_ptr<gdal::TileMatrixSet> & poTM,int & nZoomLevel,int & nAlignedLevels)118 bool COGGetWarpingCharacteristics(GDALDataset* poSrcDS,
119                                   const char * const* papszOptions,
120                                   CPLString& osResampling,
121                                   CPLString& osTargetSRS,
122                                   int& nXSize,
123                                   int& nYSize,
124                                   double& dfMinX,
125                                   double& dfMinY,
126                                   double& dfMaxX,
127                                   double& dfMaxY,
128                                   std::unique_ptr<gdal::TileMatrixSet>& poTM,
129                                   int& nZoomLevel,
130                                   int& nAlignedLevels)
131 {
132     osTargetSRS = CSLFetchNameValueDef(papszOptions, "TARGET_SRS", "");
133     CPLString osTilingScheme(CSLFetchNameValueDef(papszOptions,
134                                                   "TILING_SCHEME", "CUSTOM"));
135     if( EQUAL(osTargetSRS, "") && EQUAL(osTilingScheme, "CUSTOM") )
136         return false;
137 
138     CPLString osExtent(CSLFetchNameValueDef(papszOptions, "EXTENT", ""));
139     CPLString osRes(CSLFetchNameValueDef(papszOptions, "RES", ""));
140     if( !EQUAL(osTilingScheme, "CUSTOM") )
141     {
142         poTM = gdal::TileMatrixSet::parse(osTilingScheme);
143         if( poTM == nullptr )
144             return false;
145         if( !poTM->haveAllLevelsSameTopLeft() )
146         {
147             CPLError(CE_Failure, CPLE_NotSupported,
148                     "Unsupported tiling scheme: not all zoom levels have same top left corner");
149             return false;
150         }
151         if( !poTM->haveAllLevelsSameTileSize() )
152         {
153             CPLError(CE_Failure, CPLE_NotSupported,
154                     "Unsupported tiling scheme: not all zoom levels have same tile size");
155             return false;
156         }
157         if( poTM->hasVariableMatrixWidth() )
158         {
159             CPLError(CE_Failure, CPLE_NotSupported,
160                     "Unsupported tiling scheme: some levels have variable matrix width");
161             return false;
162         }
163         osTargetSRS = poTM->crs();
164 
165         // "Normalize" SRS as AUTH:CODE
166         OGRSpatialReference oTargetSRS;
167         oTargetSRS.SetFromUserInput(osTargetSRS);
168         const char* pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
169         const char* pszAuthName = oTargetSRS.GetAuthorityName(nullptr);
170         if( pszAuthName && pszAuthCode )
171         {
172             osTargetSRS = pszAuthName;
173             osTargetSRS += ':';
174             osTargetSRS += pszAuthCode;
175         }
176     }
177 
178     CPLStringList aosTO;
179     aosTO.SetNameValue( "DST_SRS", osTargetSRS );
180     void* hTransformArg = nullptr;
181 
182     OGRSpatialReference oTargetSRS;
183     oTargetSRS.SetFromUserInput(osTargetSRS);
184     const char* pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
185     const int nEPSGCode = pszAuthCode ? atoi(pszAuthCode) : 0;
186 
187     // Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
188     // ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
189     // EPSG:3857.
190     double adfSrcGeoTransform[6];
191     std::unique_ptr<GDALDataset> poTmpDS;
192     if( nEPSGCode == 3857 &&
193         poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None &&
194         adfSrcGeoTransform[2] == 0 &&
195         adfSrcGeoTransform[4] == 0 &&
196         adfSrcGeoTransform[5] < 0 )
197     {
198         const auto poSrcSRS = poSrcDS->GetSpatialRef();
199         if( poSrcSRS && poSrcSRS->IsGeographic() )
200         {
201             double maxLat = adfSrcGeoTransform[3];
202             double minLat = adfSrcGeoTransform[3] +
203                                     poSrcDS->GetRasterYSize() *
204                                     adfSrcGeoTransform[5];
205             // Corresponds to the latitude of below MAX_GM
206             constexpr double MAX_LAT = 85.0511287798066;
207             bool bModified = false;
208             if( maxLat > MAX_LAT )
209             {
210                 maxLat = MAX_LAT;
211                 bModified = true;
212             }
213             if( minLat < -MAX_LAT )
214             {
215                 minLat = -MAX_LAT;
216                 bModified = true;
217             }
218             if( bModified )
219             {
220                 CPLStringList aosOptions;
221                 aosOptions.AddString("-of");
222                 aosOptions.AddString("VRT");
223                 aosOptions.AddString("-projwin");
224                 aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0]));
225                 aosOptions.AddString(CPLSPrintf("%.18g", maxLat));
226                 aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0] + poSrcDS->GetRasterXSize() * adfSrcGeoTransform[1]));
227                 aosOptions.AddString(CPLSPrintf("%.18g", minLat));
228                 auto psOptions = GDALTranslateOptionsNew(aosOptions.List(), nullptr);
229                 poTmpDS.reset(GDALDataset::FromHandle(
230                     GDALTranslate("", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
231                 GDALTranslateOptionsFree(psOptions);
232                 if( poTmpDS )
233                 {
234                     hTransformArg = GDALCreateGenImgProjTransformer2(
235                         GDALDataset::FromHandle(poTmpDS.get()), nullptr, aosTO.List() );
236                     if( hTransformArg == nullptr )
237                     {
238                         return false;
239                     }
240                 }
241             }
242         }
243     }
244     if( hTransformArg == nullptr )
245     {
246         hTransformArg =
247             GDALCreateGenImgProjTransformer2( poSrcDS, nullptr, aosTO.List() );
248         if( hTransformArg == nullptr )
249         {
250             return false;
251         }
252     }
253 
254     GDALTransformerInfo* psInfo = static_cast<GDALTransformerInfo*>(hTransformArg);
255     double adfGeoTransform[6];
256     double adfExtent[4];
257 
258     if ( GDALSuggestedWarpOutput2( poSrcDS,
259                                   psInfo->pfnTransform, hTransformArg,
260                                   adfGeoTransform,
261                                   &nXSize, &nYSize,
262                                   adfExtent, 0 ) != CE_None )
263     {
264         GDALDestroyGenImgProjTransformer( hTransformArg );
265         return false;
266     }
267 
268     GDALDestroyGenImgProjTransformer( hTransformArg );
269     hTransformArg = nullptr;
270     poTmpDS.reset();
271 
272     dfMinX = adfExtent[0];
273     dfMinY = adfExtent[1];
274     dfMaxX = adfExtent[2];
275     dfMaxY = adfExtent[3];
276     double dfRes = adfGeoTransform[1];
277 
278     if( poTM )
279     {
280         const bool bInvertAxis =
281             oTargetSRS.EPSGTreatsAsLatLong() != FALSE ||
282             oTargetSRS.EPSGTreatsAsNorthingEasting() != FALSE;
283 
284         const auto& bbox = poTM->bbox();
285         if( bbox.mCrs == poTM->crs() )
286         {
287             if( dfMaxX < (bInvertAxis ? bbox.mLowerCornerY : bbox.mLowerCornerX) ||
288                 dfMinX > (bInvertAxis ? bbox.mUpperCornerY : bbox.mUpperCornerX) ||
289                 dfMaxY < (bInvertAxis ? bbox.mLowerCornerX : bbox.mLowerCornerY) ||
290                 dfMinY > (bInvertAxis ? bbox.mUpperCornerX : bbox.mUpperCornerY) )
291             {
292                 CPLError(CE_Failure, CPLE_AppDefined,
293                          "Raster extent completely outside of tile matrix set bounding box");
294                 return false;
295             }
296         }
297 
298         const auto& tmList = poTM->tileMatrixList();
299         const int nBlockSize = atoi(CSLFetchNameValueDef(
300             papszOptions, "BLOCKSIZE", CPLSPrintf("%d", tmList[0].mTileWidth)));
301         const double dfOriX = bInvertAxis ? tmList[0].mTopLeftY : tmList[0].mTopLeftX;
302         const double dfOriY = bInvertAxis ? tmList[0].mTopLeftX : tmList[0].mTopLeftY;
303         double dfComputedRes = adfGeoTransform[1];
304         double dfPrevRes = 0.0;
305         dfRes = 0.0;
306         for( ; nZoomLevel < static_cast<int>(tmList.size()); nZoomLevel++ )
307         {
308             dfRes = tmList[nZoomLevel].mResX * tmList[0].mTileWidth / nBlockSize;
309             if( dfComputedRes > dfRes || fabs( dfComputedRes - dfRes ) / dfRes <= 1e-8 )
310                 break;
311             dfPrevRes = dfRes;
312         }
313         if( nZoomLevel == static_cast<int>(tmList.size()) )
314         {
315             CPLError(CE_Failure, CPLE_AppDefined,
316                     "Could not find an appropriate zoom level");
317             return false;
318         }
319 
320         if( nZoomLevel > 0 && fabs( dfComputedRes - dfRes ) / dfRes > 1e-8 )
321         {
322             const char* pszZoomLevelStrategy = CSLFetchNameValueDef(papszOptions,
323                                                                     "ZOOM_LEVEL_STRATEGY",
324                                                                     "AUTO");
325             if( EQUAL(pszZoomLevelStrategy, "LOWER") )
326             {
327                 nZoomLevel --;
328             }
329             else if( EQUAL(pszZoomLevelStrategy, "UPPER") )
330             {
331                 /* do nothing */
332             }
333             else
334             {
335                 if( dfPrevRes / dfComputedRes < dfComputedRes / dfRes )
336                     nZoomLevel --;
337             }
338             dfRes = tmList[nZoomLevel].mResX * tmList[0].mTileWidth / nBlockSize;
339         }
340 
341         CPLDebug("COG", "Using ZOOM_LEVEL %d", nZoomLevel);
342 
343         const double dfTileExtent = dfRes * nBlockSize;
344         int nTLTileX = static_cast<int>(std::floor((dfMinX - dfOriX) / dfTileExtent + 1e-10));
345         int nTLTileY = static_cast<int>(std::floor((dfOriY - dfMaxY) / dfTileExtent + 1e-10));
346         int nBRTileX = static_cast<int>(std::ceil((dfMaxX - dfOriX) / dfTileExtent - 1e-10));
347         int nBRTileY = static_cast<int>(std::ceil((dfOriY - dfMinY) / dfTileExtent - 1e-10));
348 
349         nAlignedLevels = std::min(std::min(10, atoi(
350             CSLFetchNameValueDef(papszOptions, "ALIGNED_LEVELS", "0"))), nZoomLevel);
351         int nAccDivisor = 1;
352         for( int i = 0; i < nAlignedLevels - 1; i++ )
353         {
354             const int nCurLevel = nZoomLevel - i;
355             const double dfResRatio =
356                 tmList[nCurLevel-1].mResX / tmList[nCurLevel].mResX;
357             // Magical number that has a great number of divisors
358             // For example if previous scale denom was 50K and current one
359             // is 20K, then dfResRatio = 2.5 and dfScaledInvResRatio = 24
360             // We must then simplify 60 / 24 as 5 / 2, and make sure to
361             // align tile coordinates on multiple of the 5 numerator
362             constexpr int MAGICAL = 60;
363             const double dfScaledInvResRatio = MAGICAL / dfResRatio;
364             if( dfScaledInvResRatio < 1 || dfScaledInvResRatio > 60 ||
365                 std::abs(std::round(dfScaledInvResRatio) - dfScaledInvResRatio) > 1e-10 )
366             {
367                 CPLError(CE_Failure, CPLE_AppDefined,
368                             "Unsupported ratio of resolution for "
369                             "ALIGNED_LEVELS between zoom level %d and %d = %g",
370                             nCurLevel-1, nCurLevel, dfResRatio);
371                 return false;
372             }
373             const int nScaledInvResRatio = static_cast<int>(
374                 std::round(dfScaledInvResRatio));
375             int nNumerator = 0;
376             for( int nDivisor = nScaledInvResRatio; nDivisor >= 2; --nDivisor )
377             {
378                 if( (MAGICAL % nDivisor) == 0 && (nScaledInvResRatio % nDivisor) == 0 )
379                 {
380                     nNumerator = MAGICAL / nDivisor;
381                     break;
382                 }
383             }
384             if( nNumerator == 0 )
385             {
386                 CPLError(CE_Failure, CPLE_AppDefined,
387                             "Unsupported ratio of resolution for "
388                             "ALIGNED_LEVELS between zoom level %d and %d = %g",
389                             nCurLevel-1, nCurLevel, dfResRatio);
390                 return false;
391             }
392             nAccDivisor *= nNumerator;
393         }
394         if( nAccDivisor > 1 )
395         {
396             nTLTileX = (nTLTileX / nAccDivisor) * nAccDivisor;
397             nTLTileY = (nTLTileY / nAccDivisor) * nAccDivisor;
398             nBRTileY = ((nBRTileY + nAccDivisor - 1) / nAccDivisor) * nAccDivisor;
399             nBRTileX = ((nBRTileX + nAccDivisor - 1) / nAccDivisor) * nAccDivisor;
400         }
401 
402         if( nTLTileX < 0 || nTLTileY < 0 ||
403             nBRTileX > tmList[nZoomLevel].mMatrixWidth ||
404             nBRTileY > tmList[nZoomLevel].mMatrixHeight )
405         {
406             CPLError(CE_Warning, CPLE_AppDefined,
407                      "Raster extent partially outside of tile matrix "
408                      "bounding box. Clamping it to it");
409         }
410         nTLTileX = std::max(0, nTLTileX);
411         nTLTileY = std::max(0, nTLTileY);
412         nBRTileX = std::min(tmList[nZoomLevel].mMatrixWidth, nBRTileX);
413         nBRTileY = std::min(tmList[nZoomLevel].mMatrixHeight, nBRTileY);
414 
415         dfMinX = dfOriX + nTLTileX * dfTileExtent;
416         dfMinY = dfOriY - nBRTileY * dfTileExtent;
417         dfMaxX = dfOriX + nBRTileX * dfTileExtent;
418         dfMaxY = dfOriY - nTLTileY * dfTileExtent;
419         nXSize = static_cast<int>(std::round((dfMaxX - dfMinX) / dfRes));
420         nYSize = static_cast<int>(std::round((dfMaxY - dfMinY) / dfRes));
421     }
422     else if( !osExtent.empty() || !osRes.empty() )
423     {
424         CPLStringList aosTokens(CSLTokenizeString2(osExtent, ",", 0));
425         if( aosTokens.size() != 4 )
426         {
427             CPLError(CE_Failure, CPLE_AppDefined,
428                      "Invalid value for EXTENT");
429             return false;
430         }
431         dfMinX = CPLAtof(aosTokens[0]);
432         dfMinY = CPLAtof(aosTokens[1]);
433         dfMaxX = CPLAtof(aosTokens[2]);
434         dfMaxY = CPLAtof(aosTokens[3]);
435         if( !osRes.empty() )
436             dfRes = CPLAtof(osRes);
437     }
438 
439     nXSize = static_cast<int>(std::round((dfMaxX - dfMinX) / dfRes));
440     nYSize = static_cast<int>(std::round((dfMaxY - dfMinY) / dfRes));
441 
442     osResampling = CSLFetchNameValueDef(papszOptions,
443         "WARP_RESAMPLING",
444         CSLFetchNameValueDef(papszOptions,
445             "RESAMPLING", GetResampling(poSrcDS)));
446 
447     return true;
448 }
449 
450 // Used by gdalwarp
COGGetWarpingCharacteristics(GDALDataset * poSrcDS,const char * const * papszOptions,CPLString & osResampling,CPLString & osTargetSRS,int & nXSize,int & nYSize,double & dfMinX,double & dfMinY,double & dfMaxX,double & dfMaxY)451 bool COGGetWarpingCharacteristics(GDALDataset* poSrcDS,
452                                   const char * const* papszOptions,
453                                   CPLString& osResampling,
454                                   CPLString& osTargetSRS,
455                                   int& nXSize,
456                                   int& nYSize,
457                                   double& dfMinX,
458                                   double& dfMinY,
459                                   double& dfMaxX,
460                                   double& dfMaxY)
461 {
462     std::unique_ptr<gdal::TileMatrixSet> poTM;
463     int nZoomLevel = 0;
464     int nAlignedLevels = 0;
465     return COGGetWarpingCharacteristics(poSrcDS,
466                                         papszOptions,
467                                         osResampling,
468                                         osTargetSRS,
469                                         nXSize, nYSize,
470                                         dfMinX, dfMinY, dfMaxX, dfMaxY,
471                                         poTM, nZoomLevel, nAlignedLevels);
472 }
473 
474 /************************************************************************/
475 /*                        COGHasWarpingOptions()                        */
476 /************************************************************************/
477 
COGHasWarpingOptions(CSLConstList papszOptions)478 bool COGHasWarpingOptions(CSLConstList papszOptions)
479 {
480     return CSLFetchNameValue(papszOptions, "TARGET_SRS") != nullptr ||
481            !EQUAL(CSLFetchNameValueDef(papszOptions, "TILING_SCHEME", "CUSTOM"),
482                   "CUSTOM");
483 }
484 
485 /************************************************************************/
486 /*                      COGRemoveWarpingOptions()                       */
487 /************************************************************************/
488 
COGRemoveWarpingOptions(CPLStringList & aosOptions)489 void COGRemoveWarpingOptions(CPLStringList& aosOptions)
490 {
491     aosOptions.SetNameValue("TARGET_SRS", nullptr);
492     aosOptions.SetNameValue("TILING_SCHEME", nullptr);
493     aosOptions.SetNameValue("EXTENT", nullptr);
494     aosOptions.SetNameValue("RES", nullptr);
495     aosOptions.SetNameValue("ALIGNED_LEVELS", nullptr);
496     aosOptions.SetNameValue("ZOOM_LEVEL_STRATEGY", nullptr);
497 }
498 
499 /************************************************************************/
500 /*                        CreateReprojectedDS()                         */
501 /************************************************************************/
502 
CreateReprojectedDS(const char * pszDstFilename,GDALDataset * poSrcDS,const char * const * papszOptions,const CPLString & osResampling,const CPLString & osTargetSRS,const int nXSize,const int nYSize,const double dfMinX,const double dfMinY,const double dfMaxX,const double dfMaxY,GDALProgressFunc pfnProgress,void * pProgressData,double & dfCurPixels,double & dfTotalPixelsToProcess)503 static std::unique_ptr<GDALDataset> CreateReprojectedDS(
504                                 const char* pszDstFilename,
505                                 GDALDataset *poSrcDS,
506                                 const char * const* papszOptions,
507                                 const CPLString& osResampling,
508                                 const CPLString& osTargetSRS,
509                                 const int nXSize,
510                                 const int nYSize,
511                                 const double dfMinX,
512                                 const double dfMinY,
513                                 const double dfMaxX,
514                                 const double dfMaxY,
515                                 GDALProgressFunc pfnProgress,
516                                 void * pProgressData,
517                                 double& dfCurPixels,
518                                 double& dfTotalPixelsToProcess)
519 {
520     char** papszArg = nullptr;
521     // We could have done a warped VRT, but overview building on it might be
522     // slow, so materialize as GTiff
523     papszArg = CSLAddString(papszArg, "-of");
524     papszArg = CSLAddString(papszArg, "GTiff");
525     papszArg = CSLAddString(papszArg, "-co");
526     papszArg = CSLAddString(papszArg, "TILED=YES");
527     papszArg = CSLAddString(papszArg, "-co");
528     papszArg = CSLAddString(papszArg, "SPARSE_OK=YES");
529     const char* pszBIGTIFF = CSLFetchNameValue(papszOptions, "BIGTIFF");
530     if( pszBIGTIFF )
531     {
532         papszArg = CSLAddString(papszArg, "-co");
533         papszArg = CSLAddString(papszArg, (CPLString("BIGTIFF=") + pszBIGTIFF).c_str());
534     }
535     papszArg = CSLAddString(papszArg, "-co");
536     papszArg = CSLAddString(papszArg,
537                     HasZSTDCompression() ? "COMPRESS=ZSTD" : "COMPRESS=LZW");
538     papszArg = CSLAddString(papszArg, "-t_srs");
539     papszArg = CSLAddString(papszArg, osTargetSRS);
540     papszArg = CSLAddString(papszArg, "-te");
541     papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMinX));
542     papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMinY));
543     papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMaxX));
544     papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMaxY));
545     papszArg = CSLAddString(papszArg, "-ts");
546     papszArg = CSLAddString(papszArg, CPLSPrintf("%d", nXSize));
547     papszArg = CSLAddString(papszArg, CPLSPrintf("%d", nYSize));
548     int bHasNoData = FALSE;
549     poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoData);
550     if( !bHasNoData && CPLTestBool(CSLFetchNameValueDef(
551             papszOptions, "ADD_ALPHA", "YES")) )
552     {
553         papszArg = CSLAddString(papszArg, "-dstalpha");
554     }
555     papszArg = CSLAddString(papszArg, "-r");
556     papszArg = CSLAddString(papszArg, osResampling);
557     papszArg = CSLAddString(papszArg, "-wo");
558     papszArg = CSLAddString(papszArg, "SAMPLE_GRID=YES");
559     const char* pszNumThreads = CSLFetchNameValue(papszOptions, "NUM_THREADS");
560     if( pszNumThreads )
561     {
562         papszArg = CSLAddString(papszArg, "-wo");
563         papszArg = CSLAddString(papszArg, (CPLString("NUM_THREADS=") + pszNumThreads).c_str());
564     }
565 
566     const auto poFirstBand = poSrcDS->GetRasterBand(1);
567     const bool bHasMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
568 
569     const int nBands = poSrcDS->GetRasterCount();
570     const char* pszOverviews = CSLFetchNameValueDef(
571         papszOptions, "OVERVIEWS", "AUTO");
572     const bool bRecreateOvr = EQUAL(pszOverviews, "FORCE_USE_EXISTING") ||
573                               EQUAL(pszOverviews, "NONE");
574     dfTotalPixelsToProcess =
575         double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) +
576         ((bHasMask && !bRecreateOvr) ? double(nXSize) * nYSize / 3 : 0) +
577         (!bRecreateOvr ? double(nXSize) * nYSize * nBands / 3: 0) +
578         double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) * 4. / 3;
579 
580     auto psOptions = GDALWarpAppOptionsNew(papszArg, nullptr);
581     CSLDestroy(papszArg);
582     if( psOptions == nullptr )
583         return nullptr;
584 
585     const double dfNextPixels =
586         double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0));
587     void* pScaledProgress = GDALCreateScaledProgress(
588                 dfCurPixels / dfTotalPixelsToProcess,
589                 dfNextPixels / dfTotalPixelsToProcess,
590                 pfnProgress, pProgressData );
591     dfCurPixels = dfNextPixels;
592 
593     CPLDebug("COG", "Reprojecting source dataset: start");
594     GDALWarpAppOptionsSetProgress(psOptions, GDALScaledProgress, pScaledProgress );
595     CPLString osTmpFile(GetTmpFilename(pszDstFilename, "warped.tif.tmp"));
596     auto hSrcDS = GDALDataset::ToHandle(poSrcDS);
597     auto hRet = GDALWarp( osTmpFile, nullptr,
598                           1, &hSrcDS,
599                           psOptions, nullptr);
600     GDALWarpAppOptionsFree(psOptions);
601     CPLDebug("COG", "Reprojecting source dataset: end");
602 
603     GDALDestroyScaledProgress(pScaledProgress);
604 
605     return std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(hRet));
606 }
607 
608 /************************************************************************/
609 /*                            GDALCOGCreator                            */
610 /************************************************************************/
611 
612 struct GDALCOGCreator final
613 {
614     std::unique_ptr<GDALDataset> m_poReprojectedDS{};
615     std::unique_ptr<GDALDataset> m_poRGBMaskDS{};
616     CPLString                    m_osTmpOverviewFilename{};
617     CPLString                    m_osTmpMskOverviewFilename{};
618 
619     ~GDALCOGCreator();
620 
621     GDALDataset* Create(const char * pszFilename,
622                         GDALDataset * const poSrcDS,
623                         char ** papszOptions,
624                         GDALProgressFunc pfnProgress,
625                         void * pProgressData );
626 };
627 
628 /************************************************************************/
629 /*                    GDALCOGCreator::~GDALCOGCreator()                 */
630 /************************************************************************/
631 
~GDALCOGCreator()632 GDALCOGCreator::~GDALCOGCreator()
633 {
634     if( m_poReprojectedDS )
635     {
636         CPLString osProjectedDSName(m_poReprojectedDS->GetDescription());
637         // Destroy m_poRGBMaskDS before m_poReprojectedDS since the former
638         // references the later
639         m_poRGBMaskDS.reset();
640         m_poReprojectedDS.reset();
641         VSIUnlink(osProjectedDSName);
642     }
643     if( !m_osTmpOverviewFilename.empty() )
644     {
645         VSIUnlink(m_osTmpOverviewFilename);
646     }
647     if( !m_osTmpMskOverviewFilename.empty() )
648     {
649         VSIUnlink(m_osTmpMskOverviewFilename);
650     }
651 }
652 
653 /************************************************************************/
654 /*                    GDALCOGCreator::Create()                          */
655 /************************************************************************/
656 
Create(const char * pszFilename,GDALDataset * const poSrcDS,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)657 GDALDataset* GDALCOGCreator::Create(const char * pszFilename,
658                                     GDALDataset * const poSrcDS,
659                                     char ** papszOptions,
660                                     GDALProgressFunc pfnProgress,
661                                     void * pProgressData )
662 {
663     if( pfnProgress == nullptr )
664         pfnProgress = GDALDummyProgress;
665 
666     if( poSrcDS->GetRasterCount() == 0 )
667     {
668         CPLError(CE_Failure, CPLE_NotSupported,
669                  "COG driver does not support 0-band source raster");
670         return nullptr;
671     }
672 
673     CPLConfigOptionSetter oSetterReportDirtyBlockFlushing(
674         "GDAL_REPORT_DIRTY_BLOCK_FLUSHING", "NO", true);
675 
676     double dfCurPixels = 0;
677     double dfTotalPixelsToProcess = 0;
678     GDALDataset* poCurDS = poSrcDS;
679 
680     std::unique_ptr<gdal::TileMatrixSet> poTM;
681     int nZoomLevel = 0;
682     int nAlignedLevels = 0;
683     if( COGHasWarpingOptions(papszOptions) )
684     {
685         CPLString osTargetResampling;
686         CPLString osTargetSRS;
687         int nTargetXSize = 0;
688         int nTargetYSize = 0;
689         double dfTargetMinX = 0;
690         double dfTargetMinY = 0;
691         double dfTargetMaxX = 0;
692         double dfTargetMaxY = 0;
693         if( !COGGetWarpingCharacteristics(poCurDS, papszOptions,
694                                           osTargetResampling,
695                                           osTargetSRS,
696                                           nTargetXSize, nTargetYSize,
697                                           dfTargetMinX, dfTargetMinY,
698                                           dfTargetMaxX, dfTargetMaxY,
699                                           poTM, nZoomLevel, nAlignedLevels) )
700         {
701             return nullptr;
702         }
703 
704         // Collect information on source dataset to see if it already
705         // matches the warping specifications
706         CPLString osSrcSRS;
707         const auto poSrcSRS = poCurDS->GetSpatialRef();
708         if( poSrcSRS )
709         {
710             const char* pszAuthName = poSrcSRS->GetAuthorityName(nullptr);
711             const char* pszAuthCode = poSrcSRS->GetAuthorityCode(nullptr);
712             if( pszAuthName && pszAuthCode )
713             {
714                 osSrcSRS = pszAuthName;
715                 osSrcSRS += ':';
716                 osSrcSRS += pszAuthCode;
717             }
718         }
719         double dfSrcMinX = 0;
720         double dfSrcMinY = 0;
721         double dfSrcMaxX = 0;
722         double dfSrcMaxY = 0;
723         double adfSrcGT[6];
724         const int nSrcXSize = poCurDS->GetRasterXSize();
725         const int nSrcYSize = poCurDS->GetRasterYSize();
726         if( poCurDS->GetGeoTransform(adfSrcGT) == CE_None )
727         {
728             dfSrcMinX = adfSrcGT[0];
729             dfSrcMaxY = adfSrcGT[3];
730             dfSrcMaxX = adfSrcGT[0] + nSrcXSize * adfSrcGT[1];
731             dfSrcMinY = adfSrcGT[3] + nSrcYSize * adfSrcGT[5];
732         }
733 
734         if( nTargetXSize == nSrcXSize &&
735             nTargetYSize == nSrcYSize &&
736             osTargetSRS == osSrcSRS &&
737             fabs(dfSrcMinX - dfTargetMinX) < 1e-10 * fabs(dfSrcMinX) &&
738             fabs(dfSrcMinY - dfTargetMinY) < 1e-10 * fabs(dfSrcMinY) &&
739             fabs(dfSrcMaxX - dfTargetMaxX) < 1e-10 * fabs(dfSrcMaxX) &&
740             fabs(dfSrcMaxY - dfTargetMaxY) < 1e-10 * fabs(dfSrcMaxY) )
741         {
742             CPLDebug("COG", "Skipping reprojection step: "
743                      "source dataset matches reprojection specifications");
744         }
745         else
746         {
747             m_poReprojectedDS =
748                 CreateReprojectedDS(pszFilename, poCurDS,
749                                     papszOptions,
750                                     osTargetResampling,
751                                     osTargetSRS,
752                                     nTargetXSize, nTargetYSize,
753                                     dfTargetMinX, dfTargetMinY,
754                                     dfTargetMaxX, dfTargetMaxY,
755                                     pfnProgress, pProgressData,
756                                     dfCurPixels, dfTotalPixelsToProcess);
757             if( !m_poReprojectedDS )
758                 return nullptr;
759             poCurDS = m_poReprojectedDS.get();
760         }
761     }
762 
763     CPLString osCompress = CSLFetchNameValueDef(papszOptions, "COMPRESS", "NONE");
764     if( EQUAL(osCompress, "JPEG") &&
765         poCurDS->GetRasterCount() == 4 )
766     {
767         char** papszArg = nullptr;
768         papszArg = CSLAddString(papszArg, "-of");
769         papszArg = CSLAddString(papszArg, "VRT");
770         papszArg = CSLAddString(papszArg, "-b");
771         papszArg = CSLAddString(papszArg, "1");
772         papszArg = CSLAddString(papszArg, "-b");
773         papszArg = CSLAddString(papszArg, "2");
774         papszArg = CSLAddString(papszArg, "-b");
775         papszArg = CSLAddString(papszArg, "3");
776         papszArg = CSLAddString(papszArg, "-mask");
777         papszArg = CSLAddString(papszArg, "4");
778         GDALTranslateOptions* psOptions = GDALTranslateOptionsNew(papszArg, nullptr);
779         CSLDestroy(papszArg);
780         GDALDatasetH hRGBMaskDS = GDALTranslate("",
781                                                 GDALDataset::ToHandle(poCurDS),
782                                                 psOptions,
783                                                 nullptr);
784         GDALTranslateOptionsFree(psOptions);
785         if( !hRGBMaskDS )
786         {
787             return nullptr;
788         }
789         m_poRGBMaskDS.reset( GDALDataset::FromHandle(hRGBMaskDS) );
790         poCurDS = m_poRGBMaskDS.get();
791     }
792 
793     const int nBands = poCurDS->GetRasterCount();
794     const int nXSize = poCurDS->GetRasterXSize();
795     const int nYSize = poCurDS->GetRasterYSize();
796 
797     CPLString osBlockSize(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", ""));
798     if( osBlockSize.empty() )
799     {
800         if( poTM )
801         {
802             osBlockSize.Printf("%d", poTM->tileMatrixList()[0].mTileWidth);
803         }
804         else
805         {
806             osBlockSize = "512";
807         }
808     }
809 
810     const int nOvrThresholdSize = atoi(osBlockSize);
811 
812     const auto poFirstBand = poCurDS->GetRasterBand(1);
813     const bool bHasMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
814 
815     CPLString osOverviews = CSLFetchNameValueDef(
816         papszOptions, "OVERVIEWS", "AUTO");
817     const bool bRecreateOvr = EQUAL(osOverviews, "FORCE_USE_EXISTING") ||
818                               EQUAL(osOverviews, "NONE");
819     const bool bGenerateMskOvr =
820         !bRecreateOvr &&
821         bHasMask &&
822         (nXSize > nOvrThresholdSize || nYSize > nOvrThresholdSize) &&
823         (EQUAL(osOverviews, "IGNORE_EXISTING") ||
824          poFirstBand->GetMaskBand()->GetOverviewCount() == 0);
825     const bool bGenerateOvr =
826         !bRecreateOvr &&
827         (nXSize > nOvrThresholdSize || nYSize > nOvrThresholdSize) &&
828         (EQUAL(osOverviews, "IGNORE_EXISTING") ||
829          poFirstBand->GetOverviewCount() == 0);
830 
831     std::vector<std::pair<int, int>> asOverviewDims;
832     int nTmpXSize = nXSize;
833     int nTmpYSize = nYSize;
834     if( poTM )
835     {
836         const auto& tmList = poTM->tileMatrixList();
837         int nCurLevel = nZoomLevel;
838         while( nTmpXSize > nOvrThresholdSize || nTmpYSize > nOvrThresholdSize )
839         {
840             const double dfResRatio = (nCurLevel >= 1) ?
841                 tmList[nCurLevel-1].mResX / tmList[nCurLevel].mResX : 2;
842             nTmpXSize = static_cast<int>(nTmpXSize / dfResRatio + 0.5);
843             nTmpYSize = static_cast<int>(nTmpYSize / dfResRatio + 0.5);
844             asOverviewDims.push_back(std::pair<int,int>(nTmpXSize, nTmpYSize));
845             nCurLevel --;
846         }
847     }
848     else
849     {
850         while( nTmpXSize > nOvrThresholdSize || nTmpYSize > nOvrThresholdSize )
851         {
852             nTmpXSize /= 2;
853             nTmpYSize /= 2;
854             asOverviewDims.push_back(std::pair<int,int>(nTmpXSize, nTmpYSize));
855         }
856     }
857 
858     if( dfTotalPixelsToProcess == 0.0 )
859     {
860         dfTotalPixelsToProcess =
861             (bGenerateMskOvr ? double(nXSize) * nYSize / 3 : 0) +
862             (bGenerateOvr ? double(nXSize) * nYSize * nBands / 3: 0) +
863             double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) * 4. / 3;
864     }
865 
866     CPLStringList aosOverviewOptions;
867     aosOverviewOptions.SetNameValue("COMPRESS",
868         CPLGetConfigOption("COG_TMP_COMPRESSION", // only for debug purposes
869                         HasZSTDCompression() ? "ZSTD" : "LZW"));
870     aosOverviewOptions.SetNameValue("NUM_THREADS",
871                         CSLFetchNameValue(papszOptions, "NUM_THREADS"));
872     aosOverviewOptions.SetNameValue("BIGTIFF", "YES");
873     aosOverviewOptions.SetNameValue("SPARSE_OK", "YES");
874 
875     if( bGenerateMskOvr )
876     {
877         CPLDebug("COG", "Generating overviews of the mask: start");
878         m_osTmpMskOverviewFilename = GetTmpFilename(pszFilename, "msk.ovr.tmp");
879         GDALRasterBand* poSrcMask = poFirstBand->GetMaskBand();
880         const char* pszResampling = CSLFetchNameValueDef(papszOptions,
881             "OVERVIEW_RESAMPLING",
882                 CSLFetchNameValueDef(papszOptions,
883                     "RESAMPLING", GetResampling(poSrcDS)));
884 
885         double dfNextPixels = dfCurPixels + double(nXSize) * nYSize / 3;
886         void* pScaledProgress = GDALCreateScaledProgress(
887                 dfCurPixels / dfTotalPixelsToProcess,
888                 dfNextPixels / dfTotalPixelsToProcess,
889                 pfnProgress, pProgressData );
890         dfCurPixels = dfNextPixels;
891 
892         // Used by GDALRegenerateOverviews() and GDALRegenerateOverviewsMultiBand()
893         CPLConfigOptionSetter oSetterRegeneratedBandIsMask(
894             "GDAL_REGENERATED_BAND_IS_MASK", "YES", true);
895 
896         CPLErr eErr = GTIFFBuildOverviewsEx(
897             m_osTmpMskOverviewFilename,
898             1, &poSrcMask,
899             static_cast<int>(asOverviewDims.size()),
900             nullptr, asOverviewDims.data(),
901             pszResampling,
902             aosOverviewOptions.List(),
903             GDALScaledProgress, pScaledProgress );
904         CPLDebug("COG", "Generating overviews of the mask: end");
905 
906         GDALDestroyScaledProgress(pScaledProgress);
907         if( eErr != CE_None )
908         {
909             return nullptr;
910         }
911     }
912 
913     if( bGenerateOvr )
914     {
915         CPLDebug("COG", "Generating overviews of the imagery: start");
916         m_osTmpOverviewFilename = GetTmpFilename(pszFilename, "ovr.tmp");
917         std::vector<GDALRasterBand*> apoSrcBands;
918         for( int i = 0; i < nBands; i++ )
919             apoSrcBands.push_back( poCurDS->GetRasterBand(i+1) );
920         const char* pszResampling = CSLFetchNameValueDef(papszOptions,
921             "OVERVIEW_RESAMPLING",
922                 CSLFetchNameValueDef(papszOptions,
923                     "RESAMPLING", GetResampling(poSrcDS)));
924 
925         double dfNextPixels = dfCurPixels + double(nXSize) * nYSize * nBands / 3;
926         void* pScaledProgress = GDALCreateScaledProgress(
927                 dfCurPixels / dfTotalPixelsToProcess,
928                 dfNextPixels / dfTotalPixelsToProcess,
929                 pfnProgress, pProgressData );
930         dfCurPixels = dfNextPixels;
931 
932         if( nBands > 1 )
933         {
934             aosOverviewOptions.SetNameValue("INTERLEAVE", "PIXEL");
935         }
936         if( !m_osTmpMskOverviewFilename.empty() )
937         {
938             aosOverviewOptions.SetNameValue("MASK_OVERVIEW_DATASET",
939                                             m_osTmpMskOverviewFilename);
940         }
941         CPLErr eErr = GTIFFBuildOverviewsEx(
942             m_osTmpOverviewFilename,
943             nBands, &apoSrcBands[0],
944             static_cast<int>(asOverviewDims.size()),
945             nullptr, asOverviewDims.data(),
946             pszResampling,
947             aosOverviewOptions.List(),
948             GDALScaledProgress, pScaledProgress );
949         CPLDebug("COG", "Generating overviews of the imagery: end");
950 
951         GDALDestroyScaledProgress(pScaledProgress);
952         if( eErr != CE_None )
953         {
954             return nullptr;
955         }
956     }
957 
958     CPLStringList aosOptions;
959     aosOptions.SetNameValue("COPY_SRC_OVERVIEWS", "YES");
960     aosOptions.SetNameValue("COMPRESS", osCompress);
961     aosOptions.SetNameValue("TILED", "YES");
962     aosOptions.SetNameValue("BLOCKXSIZE", osBlockSize);
963     aosOptions.SetNameValue("BLOCKYSIZE", osBlockSize);
964     const char* pszPredictor = CSLFetchNameValueDef(papszOptions, "PREDICTOR", "FALSE");
965     const char* pszPredictorValue = GetPredictor(poSrcDS, pszPredictor);
966     if (pszPredictorValue != nullptr)
967     {
968         aosOptions.SetNameValue("PREDICTOR", pszPredictorValue);
969     }
970 
971     const char* pszQuality = CSLFetchNameValue(papszOptions, "QUALITY");
972     if( EQUAL(osCompress, "JPEG") )
973     {
974         aosOptions.SetNameValue("JPEG_QUALITY", pszQuality);
975         if( nBands == 3 )
976             aosOptions.SetNameValue("PHOTOMETRIC", "YCBCR");
977     }
978     else if( EQUAL(osCompress, "WEBP") )
979     {
980         if( pszQuality && atoi(pszQuality) == 100 )
981             aosOptions.SetNameValue("WEBP_LOSSLESS", "YES");
982         aosOptions.SetNameValue("WEBP_LEVEL", pszQuality);
983     }
984     else if( EQUAL(osCompress, "DEFLATE") || EQUAL(osCompress, "LERC_DEFLATE") )
985     {
986         aosOptions.SetNameValue("ZLEVEL",
987                                 CSLFetchNameValue(papszOptions, "LEVEL"));
988     }
989     else if( EQUAL(osCompress, "ZSTD") || EQUAL(osCompress, "LERC_ZSTD")  )
990     {
991         aosOptions.SetNameValue("ZSTD_LEVEL",
992                                 CSLFetchNameValue(papszOptions, "LEVEL"));
993     }
994 
995     if( STARTS_WITH_CI(osCompress, "LERC") )
996     {
997         aosOptions.SetNameValue("MAX_Z_ERROR",
998                                 CSLFetchNameValue(papszOptions, "MAX_Z_ERROR"));
999     }
1000     aosOptions.SetNameValue("BIGTIFF",
1001                                 CSLFetchNameValue(papszOptions, "BIGTIFF"));
1002     aosOptions.SetNameValue("NUM_THREADS",
1003                                 CSLFetchNameValue(papszOptions, "NUM_THREADS"));
1004     aosOptions.SetNameValue("GEOTIFF_VERSION",
1005                             CSLFetchNameValue(papszOptions, "GEOTIFF_VERSION"));
1006     aosOptions.SetNameValue("SPARSE_OK",
1007                             CSLFetchNameValue(papszOptions, "SPARSE_OK"));
1008 
1009     if( EQUAL( osOverviews, "NONE") )
1010     {
1011         aosOptions.SetNameValue("@OVERVIEW_DATASET", "");
1012     }
1013     else
1014     {
1015         if( !m_osTmpOverviewFilename.empty() )
1016         {
1017             aosOptions.SetNameValue("@OVERVIEW_DATASET", m_osTmpOverviewFilename);
1018         }
1019         if( !m_osTmpMskOverviewFilename.empty() )
1020         {
1021             aosOptions.SetNameValue("@MASK_OVERVIEW_DATASET", m_osTmpMskOverviewFilename);
1022         }
1023     }
1024 
1025     const CPLString osTilingScheme(CSLFetchNameValueDef(papszOptions,
1026                                                 "TILING_SCHEME", "CUSTOM"));
1027     if( osTilingScheme != "CUSTOM" )
1028     {
1029          aosOptions.SetNameValue("@TILING_SCHEME_NAME", osTilingScheme);
1030          aosOptions.SetNameValue("@TILING_SCHEME_ZOOM_LEVEL",
1031                                  CPLSPrintf("%d", nZoomLevel));
1032          if( nAlignedLevels > 0 )
1033          {
1034              aosOptions.SetNameValue("@TILING_SCHEME_ALIGNED_LEVELS",
1035                                      CPLSPrintf("%d", nAlignedLevels));
1036          }
1037     }
1038     const char* pszOverviewCompress = CSLFetchNameValue(papszOptions, "OVERVIEW_COMPRESS");
1039 
1040     CPLConfigOptionSetter ovrCompressSetter("COMPRESS_OVERVIEW", pszOverviewCompress, true);
1041     CPLConfigOptionSetter ovrQualityJpegSetter("JPEG_QUALITY_OVERVIEW", CSLFetchNameValue(papszOptions, "OVERVIEW_QUALITY"), true);
1042     CPLConfigOptionSetter ovrQualityWebpSetter("WEBP_LEVEL_OVERVIEW", CSLFetchNameValue(papszOptions, "OVERVIEW_QUALITY"), true);
1043 
1044     std::unique_ptr<CPLConfigOptionSetter> poPhotometricSetter;
1045     if (pszOverviewCompress != nullptr && nBands == 3 && EQUAL(pszOverviewCompress, "JPEG") )
1046     {
1047         poPhotometricSetter.reset(new CPLConfigOptionSetter("PHOTOMETRIC_OVERVIEW", "YCBCR", true));
1048     }
1049 
1050     const char* osOvrPredictor = CSLFetchNameValueDef(papszOptions, "OVERVIEW_PREDICTOR", "FALSE");
1051     const char* pszOvrPredictorValue = GetPredictor(poSrcDS, osOvrPredictor);
1052     CPLConfigOptionSetter ovrPredictorSetter("PREDICTOR_OVERVIEW", pszOvrPredictorValue, true);
1053 
1054     GDALDriver* poGTiffDrv = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
1055     if( !poGTiffDrv )
1056         return nullptr;
1057     void* pScaledProgress = GDALCreateScaledProgress(
1058             dfCurPixels / dfTotalPixelsToProcess,
1059             1.0,
1060             pfnProgress, pProgressData );
1061 
1062     CPLConfigOptionSetter oSetterInternalMask(
1063         "GDAL_TIFF_INTERNAL_MASK", "YES", false);
1064 
1065     CPLDebug("COG", "Generating final product: start");
1066     auto poRet = poGTiffDrv->CreateCopy(pszFilename, poCurDS, false,
1067                                         aosOptions.List(),
1068                                         GDALScaledProgress, pScaledProgress);
1069 
1070     GDALDestroyScaledProgress(pScaledProgress);
1071 
1072     if( poRet )
1073         poRet->FlushCache();
1074 
1075     CPLDebug("COG", "Generating final product: end");
1076     return poRet;
1077 }
1078 
1079 /************************************************************************/
1080 /*                            COGCreateCopy()                           */
1081 /************************************************************************/
1082 
COGCreateCopy(const char * pszFilename,GDALDataset * poSrcDS,int,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)1083 static GDALDataset* COGCreateCopy( const char * pszFilename,
1084                                    GDALDataset *poSrcDS,
1085                                    int /*bStrict*/, char ** papszOptions,
1086                                    GDALProgressFunc pfnProgress,
1087                                    void * pProgressData )
1088 {
1089     return GDALCOGCreator().Create(pszFilename, poSrcDS,
1090                                    papszOptions, pfnProgress, pProgressData);
1091 }
1092 
1093 /************************************************************************/
1094 /*                          GDALRegister_COG()                          */
1095 /************************************************************************/
1096 
1097 class GDALCOGDriver final: public GDALDriver
1098 {
1099         bool m_bInitialized = false;
1100 
1101         bool bHasLZW = false;
1102         bool bHasDEFLATE = false;
1103         bool bHasLZMA = false;
1104         bool bHasZSTD = false;
1105         bool bHasJPEG = false;
1106         bool bHasWebP = false;
1107         bool bHasLERC = false;
1108         std::string osCompressValues{};
1109 
1110         void InitializeCreationOptionList();
1111 
1112     public:
1113         GDALCOGDriver();
1114 
GetMetadataItem(const char * pszName,const char * pszDomain)1115         const char* GetMetadataItem(const char* pszName, const char* pszDomain) override
1116         {
1117             if( EQUAL(pszName, GDAL_DMD_CREATIONOPTIONLIST) )
1118             {
1119                 InitializeCreationOptionList();
1120             }
1121             return GDALDriver::GetMetadataItem(pszName, pszDomain);
1122         }
1123 
GetMetadata(const char * pszDomain)1124         char** GetMetadata(const char* pszDomain) override
1125         {
1126             InitializeCreationOptionList();
1127             return GDALDriver::GetMetadata(pszDomain);
1128         }
1129 };
1130 
GDALCOGDriver()1131 GDALCOGDriver::GDALCOGDriver()
1132 {
1133     // We could defer this in InitializeCreationOptionList() but with currently
1134     // released libtiff versions where there was a bug (now fixed) in
1135     // TIFFGetConfiguredCODECs(), this wouldn't work properly if the LERC codec
1136     // had been registered in between
1137     osCompressValues = GTiffGetCompressValues(
1138         bHasLZW, bHasDEFLATE, bHasLZMA, bHasZSTD, bHasJPEG, bHasWebP, bHasLERC,
1139         true /* bForCOG */);
1140 }
1141 
InitializeCreationOptionList()1142 void GDALCOGDriver::InitializeCreationOptionList()
1143 {
1144     if( m_bInitialized )
1145         return;
1146     m_bInitialized = true;
1147 
1148     CPLString osOptions;
1149     osOptions = "<CreationOptionList>"
1150                 "   <Option name='COMPRESS' type='string-select'>";
1151     osOptions += osCompressValues;
1152     osOptions += "   </Option>";
1153 
1154     osOptions += "   <Option name='OVERVIEW_COMPRESS' type='string-select'>";
1155     osOptions += osCompressValues;
1156     osOptions += "   </Option>";
1157 
1158     if( bHasLZW || bHasDEFLATE || bHasZSTD )
1159     {
1160         const char* osPredictorOptions =  "     <Value>YES</Value>"
1161                      "     <Value>NO</Value>"
1162                      "     <Value alias='2'>STANDARD</Value>"
1163                      "     <Value alias='3'>FLOATING_POINT</Value>";
1164 
1165         osOptions += "   <Option name='LEVEL' type='int' "
1166             "description='DEFLATE/ZSTD compression level: 1 (fastest)'/>";
1167 
1168         osOptions += "   <Option name='PREDICTOR' type='string-select' default='FALSE'>";
1169         osOptions += osPredictorOptions;
1170         osOptions += "   </Option>"
1171                      "   <Option name='OVERVIEW_PREDICTOR' type='string-select' default='FALSE'>";
1172         osOptions += osPredictorOptions;
1173         osOptions += "   </Option>";
1174     }
1175     if( bHasJPEG || bHasWebP )
1176     {
1177         osOptions += "   <Option name='QUALITY' type='int' "
1178                      "description='JPEG/WEBP quality 1-100' default='75'/>"
1179                      "   <Option name='OVERVIEW_QUALITY' type='int' "
1180                      "description='Overview JPEG/WEBP quality 1-100' default='75'/>";
1181     }
1182     if( bHasLERC )
1183     {
1184         osOptions += ""
1185 "   <Option name='MAX_Z_ERROR' type='float' description='Maximum error for LERC compression' default='0'/>";
1186     }
1187     osOptions +=
1188 "   <Option name='NUM_THREADS' type='string' "
1189         "description='Number of worker threads for compression. "
1190         "Can be set to ALL_CPUS' default='1'/>"
1191 "   <Option name='BLOCKSIZE' type='int' "
1192         "description='Tile size in pixels' min='128' default='512'/>"
1193 "   <Option name='BIGTIFF' type='string-select' description='"
1194         "Force creation of BigTIFF file'>"
1195 "     <Value>YES</Value>"
1196 "     <Value>NO</Value>"
1197 "     <Value>IF_NEEDED</Value>"
1198 "     <Value>IF_SAFER</Value>"
1199 "   </Option>"
1200 "   <Option name='RESAMPLING' type='string' "
1201         "description='Resampling method for overviews or warping'/>"
1202 "   <Option name='OVERVIEW_RESAMPLING' type='string' "
1203         "description='Resampling method for overviews'/>"
1204 "   <Option name='WARP_RESAMPLING' type='string' "
1205         "description='Resampling method for warping'/>"
1206 "   <Option name='OVERVIEWS' type='string-select' description='"
1207         "Behavior regarding overviews'>"
1208 "     <Value>AUTO</Value>"
1209 "     <Value>IGNORE_EXISTING</Value>"
1210 "     <Value>FORCE_USE_EXISTING</Value>"
1211 "     <Value>NONE</Value>"
1212 "   </Option>"
1213 "  <Option name='TILING_SCHEME' type='string' description='"
1214         "Which tiling scheme to use pre-defined value or custom inline/outline "
1215         "JSON definition' default='CUSTOM'>"
1216 "    <Value>CUSTOM</Value>";
1217 
1218     const auto tmsList = gdal::TileMatrixSet::listPredefinedTileMatrixSets();
1219     for( const auto& tmsName: tmsList )
1220     {
1221         const auto poTM = gdal::TileMatrixSet::parse(tmsName.c_str());
1222         if( poTM &&
1223             poTM->haveAllLevelsSameTopLeft() &&
1224             poTM->haveAllLevelsSameTileSize() &&
1225             !poTM->hasVariableMatrixWidth() )
1226         {
1227             osOptions += "    <Value>";
1228             osOptions += tmsName;
1229             osOptions += "</Value>";
1230         }
1231     }
1232 
1233     osOptions +=
1234 "  </Option>"
1235 "  <Option name='ZOOM_LEVEL_STRATEGY' type='string-select' "
1236         "description='Strategy to determine zoom level. "
1237         "Only used for TILING_SCHEME != CUSTOM' default='AUTO'>"
1238 "    <Value>AUTO</Value>"
1239 "    <Value>LOWER</Value>"
1240 "    <Value>UPPER</Value>"
1241 "  </Option>"
1242 "   <Option name='TARGET_SRS' type='string' "
1243         "description='Target SRS as EPSG:XXXX, WKT or PROJ string for reprojection'/>"
1244 "  <Option name='RES' type='float' description='"
1245         "Target resolution for reprojection'/>"
1246 "  <Option name='EXTENT' type='string' description='"
1247         "Target extent as minx,miny,maxx,maxy for reprojection'/>"
1248 "  <Option name='ALIGNED_LEVELS' type='int' description='"
1249         "Number of overview levels for which the tiles from GeoTIFF and the "
1250         "specified tiling scheme match'/>"
1251 "  <Option name='ADD_ALPHA' type='boolean' description='Can be set to NO to "
1252         "disable the addition of an alpha band in case of reprojection' default='YES'/>"
1253 #if LIBGEOTIFF_VERSION >= 1600
1254 "   <Option name='GEOTIFF_VERSION' type='string-select' default='AUTO' description='Which version of GeoTIFF must be used'>"
1255 "       <Value>AUTO</Value>"
1256 "       <Value>1.0</Value>"
1257 "       <Value>1.1</Value>"
1258 "   </Option>"
1259 #endif
1260 "   <Option name='SPARSE_OK' type='boolean' description='Should empty blocks be omitted on disk?' default='FALSE'/>"
1261 "</CreationOptionList>";
1262 
1263     SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST, osOptions.c_str());
1264 }
1265 
GDALRegister_COG()1266 void GDALRegister_COG()
1267 
1268 {
1269     if( GDALGetDriverByName( "COG" ) != nullptr )
1270         return;
1271 
1272     auto poDriver = new GDALCOGDriver();
1273     poDriver->SetDescription( "COG" );
1274     poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
1275     poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Cloud optimized GeoTIFF generator" );
1276     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/cog.html" );
1277 
1278     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1279                                "Byte UInt16 Int16 UInt32 Int32 Float32 "
1280                                "Float64 CInt16 CInt32 CFloat32 CFloat64" );
1281 
1282     poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1283 
1284     poDriver->pfnCreateCopy = COGCreateCopy;
1285 
1286     GetGDALDriverManager()->RegisterDriver( poDriver );
1287 }
1288