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