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