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 43989 : const char *GDALCOGDriver::GetMetadataItem(const char *pszName,
1443 : const char *pszDomain)
1444 : {
1445 87978 : std::lock_guard oLock(m_oMutex);
1446 43989 : if (EQUAL(pszName, GDAL_DMD_CREATIONOPTIONLIST))
1447 : {
1448 384 : InitializeCreationOptionList();
1449 : }
1450 87978 : 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 : }
|
