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