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