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