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