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 "gdal_priv.h"
30 #include "gtiff.h"
31 #include "gt_overview.h"
32 #include "gdal_utils.h"
33 #include "gdalwarper.h"
34 #include "cogdriver.h"
35 #include "geotiff.h"
36
37 #include "tilematrixset.hpp"
38
39 #include <algorithm>
40 #include <memory>
41 #include <vector>
42
43 extern "C" CPL_DLL void GDALRegister_COG();
44
45 /************************************************************************/
46 /* HasZSTDCompression() */
47 /************************************************************************/
48
HasZSTDCompression()49 static bool HasZSTDCompression()
50 {
51 TIFFCodec *codecs = TIFFGetConfiguredCODECs();
52 bool bHasZSTD = false;
53 for( TIFFCodec *c = codecs; c->name; ++c )
54 {
55 if( c->scheme == COMPRESSION_ZSTD )
56 {
57 bHasZSTD = true;
58 break;
59 }
60 }
61 _TIFFfree( codecs );
62 return bHasZSTD;
63 }
64
65 /************************************************************************/
66 /* GetTmpFilename() */
67 /************************************************************************/
68
GetTmpFilename(const char * pszFilename,const char * pszExt)69 static CPLString GetTmpFilename(const char* pszFilename,
70 const char* pszExt)
71 {
72 CPLString osTmpFilename;
73 osTmpFilename.Printf("%s.%s", pszFilename, pszExt);
74 VSIUnlink(osTmpFilename);
75 return osTmpFilename;
76 }
77
78 /************************************************************************/
79 /* GetResampling() */
80 /************************************************************************/
81
GetResampling(GDALDataset * poSrcDS)82 static const char* GetResampling(GDALDataset* poSrcDS)
83 {
84 return poSrcDS->GetRasterBand(1)->GetColorTable() ? "NEAREST" : "CUBIC";
85 }
86
87 /************************************************************************/
88 /* GetPredictor() */
89 /************************************************************************/
GetPredictor(GDALDataset * poSrcDS,const char * pszPredictor)90 static const char* GetPredictor(GDALDataset* poSrcDS,
91 const char* pszPredictor)
92 {
93 if (pszPredictor == nullptr) return nullptr;
94
95 if( EQUAL(pszPredictor, "YES") || EQUAL(pszPredictor, "ON") || EQUAL(pszPredictor, "TRUE") )
96 {
97 if( GDALDataTypeIsFloating(poSrcDS->GetRasterBand(1)->GetRasterDataType()) )
98 return "3";
99 else
100 return "2";
101 }
102 else if( EQUAL(pszPredictor, "STANDARD") || EQUAL(pszPredictor, "2") )
103 {
104 return "2";
105 }
106 else if( EQUAL(pszPredictor, "FLOATING_POINT") || EQUAL(pszPredictor, "3") )
107 {
108 return "3";
109 }
110 return nullptr;
111 }
112
113 /************************************************************************/
114 /* COGGetWarpingCharacteristics() */
115 /************************************************************************/
116
117 static
COGGetWarpingCharacteristics(GDALDataset * poSrcDS,const char * const * papszOptions,CPLString & osResampling,CPLString & osTargetSRS,int & nXSize,int & nYSize,double & dfMinX,double & dfMinY,double & dfMaxX,double & dfMaxY,std::unique_ptr<gdal::TileMatrixSet> & poTM,int & nZoomLevel,int & nAlignedLevels)118 bool COGGetWarpingCharacteristics(GDALDataset* poSrcDS,
119 const char * const* papszOptions,
120 CPLString& osResampling,
121 CPLString& osTargetSRS,
122 int& nXSize,
123 int& nYSize,
124 double& dfMinX,
125 double& dfMinY,
126 double& dfMaxX,
127 double& dfMaxY,
128 std::unique_ptr<gdal::TileMatrixSet>& poTM,
129 int& nZoomLevel,
130 int& nAlignedLevels)
131 {
132 osTargetSRS = CSLFetchNameValueDef(papszOptions, "TARGET_SRS", "");
133 CPLString osTilingScheme(CSLFetchNameValueDef(papszOptions,
134 "TILING_SCHEME", "CUSTOM"));
135 if( EQUAL(osTargetSRS, "") && EQUAL(osTilingScheme, "CUSTOM") )
136 return false;
137
138 CPLString osExtent(CSLFetchNameValueDef(papszOptions, "EXTENT", ""));
139 CPLString osRes(CSLFetchNameValueDef(papszOptions, "RES", ""));
140 if( !EQUAL(osTilingScheme, "CUSTOM") )
141 {
142 poTM = gdal::TileMatrixSet::parse(osTilingScheme);
143 if( poTM == nullptr )
144 return false;
145 if( !poTM->haveAllLevelsSameTopLeft() )
146 {
147 CPLError(CE_Failure, CPLE_NotSupported,
148 "Unsupported tiling scheme: not all zoom levels have same top left corner");
149 return false;
150 }
151 if( !poTM->haveAllLevelsSameTileSize() )
152 {
153 CPLError(CE_Failure, CPLE_NotSupported,
154 "Unsupported tiling scheme: not all zoom levels have same tile size");
155 return false;
156 }
157 if( poTM->hasVariableMatrixWidth() )
158 {
159 CPLError(CE_Failure, CPLE_NotSupported,
160 "Unsupported tiling scheme: some levels have variable matrix width");
161 return false;
162 }
163 osTargetSRS = poTM->crs();
164
165 // "Normalize" SRS as AUTH:CODE
166 OGRSpatialReference oTargetSRS;
167 oTargetSRS.SetFromUserInput(osTargetSRS);
168 const char* pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
169 const char* pszAuthName = oTargetSRS.GetAuthorityName(nullptr);
170 if( pszAuthName && pszAuthCode )
171 {
172 osTargetSRS = pszAuthName;
173 osTargetSRS += ':';
174 osTargetSRS += pszAuthCode;
175 }
176 }
177
178 CPLStringList aosTO;
179 aosTO.SetNameValue( "DST_SRS", osTargetSRS );
180 void* hTransformArg = nullptr;
181
182 OGRSpatialReference oTargetSRS;
183 oTargetSRS.SetFromUserInput(osTargetSRS);
184 const char* pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
185 const int nEPSGCode = pszAuthCode ? atoi(pszAuthCode) : 0;
186
187 // Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
188 // ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
189 // EPSG:3857.
190 double adfSrcGeoTransform[6];
191 std::unique_ptr<GDALDataset> poTmpDS;
192 if( nEPSGCode == 3857 &&
193 poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None &&
194 adfSrcGeoTransform[2] == 0 &&
195 adfSrcGeoTransform[4] == 0 &&
196 adfSrcGeoTransform[5] < 0 )
197 {
198 const auto poSrcSRS = poSrcDS->GetSpatialRef();
199 if( poSrcSRS && poSrcSRS->IsGeographic() )
200 {
201 double maxLat = adfSrcGeoTransform[3];
202 double minLat = adfSrcGeoTransform[3] +
203 poSrcDS->GetRasterYSize() *
204 adfSrcGeoTransform[5];
205 // Corresponds to the latitude of below MAX_GM
206 constexpr double MAX_LAT = 85.0511287798066;
207 bool bModified = false;
208 if( maxLat > MAX_LAT )
209 {
210 maxLat = MAX_LAT;
211 bModified = true;
212 }
213 if( minLat < -MAX_LAT )
214 {
215 minLat = -MAX_LAT;
216 bModified = true;
217 }
218 if( bModified )
219 {
220 CPLStringList aosOptions;
221 aosOptions.AddString("-of");
222 aosOptions.AddString("VRT");
223 aosOptions.AddString("-projwin");
224 aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0]));
225 aosOptions.AddString(CPLSPrintf("%.18g", maxLat));
226 aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0] + poSrcDS->GetRasterXSize() * adfSrcGeoTransform[1]));
227 aosOptions.AddString(CPLSPrintf("%.18g", minLat));
228 auto psOptions = GDALTranslateOptionsNew(aosOptions.List(), nullptr);
229 poTmpDS.reset(GDALDataset::FromHandle(
230 GDALTranslate("", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
231 GDALTranslateOptionsFree(psOptions);
232 if( poTmpDS )
233 {
234 hTransformArg = GDALCreateGenImgProjTransformer2(
235 GDALDataset::FromHandle(poTmpDS.get()), nullptr, aosTO.List() );
236 if( hTransformArg == nullptr )
237 {
238 return false;
239 }
240 }
241 }
242 }
243 }
244 if( hTransformArg == nullptr )
245 {
246 hTransformArg =
247 GDALCreateGenImgProjTransformer2( poSrcDS, nullptr, aosTO.List() );
248 if( hTransformArg == nullptr )
249 {
250 return false;
251 }
252 }
253
254 GDALTransformerInfo* psInfo = static_cast<GDALTransformerInfo*>(hTransformArg);
255 double adfGeoTransform[6];
256 double adfExtent[4];
257
258 if ( GDALSuggestedWarpOutput2( poSrcDS,
259 psInfo->pfnTransform, hTransformArg,
260 adfGeoTransform,
261 &nXSize, &nYSize,
262 adfExtent, 0 ) != CE_None )
263 {
264 GDALDestroyGenImgProjTransformer( hTransformArg );
265 return false;
266 }
267
268 GDALDestroyGenImgProjTransformer( hTransformArg );
269 hTransformArg = nullptr;
270 poTmpDS.reset();
271
272 dfMinX = adfExtent[0];
273 dfMinY = adfExtent[1];
274 dfMaxX = adfExtent[2];
275 dfMaxY = adfExtent[3];
276 double dfRes = adfGeoTransform[1];
277
278 if( poTM )
279 {
280 const bool bInvertAxis =
281 oTargetSRS.EPSGTreatsAsLatLong() != FALSE ||
282 oTargetSRS.EPSGTreatsAsNorthingEasting() != FALSE;
283
284 const auto& bbox = poTM->bbox();
285 if( bbox.mCrs == poTM->crs() )
286 {
287 if( dfMaxX < (bInvertAxis ? bbox.mLowerCornerY : bbox.mLowerCornerX) ||
288 dfMinX > (bInvertAxis ? bbox.mUpperCornerY : bbox.mUpperCornerX) ||
289 dfMaxY < (bInvertAxis ? bbox.mLowerCornerX : bbox.mLowerCornerY) ||
290 dfMinY > (bInvertAxis ? bbox.mUpperCornerX : bbox.mUpperCornerY) )
291 {
292 CPLError(CE_Failure, CPLE_AppDefined,
293 "Raster extent completely outside of tile matrix set bounding box");
294 return false;
295 }
296 }
297
298 const auto& tmList = poTM->tileMatrixList();
299 const int nBlockSize = atoi(CSLFetchNameValueDef(
300 papszOptions, "BLOCKSIZE", CPLSPrintf("%d", tmList[0].mTileWidth)));
301 const double dfOriX = bInvertAxis ? tmList[0].mTopLeftY : tmList[0].mTopLeftX;
302 const double dfOriY = bInvertAxis ? tmList[0].mTopLeftX : tmList[0].mTopLeftY;
303 double dfComputedRes = adfGeoTransform[1];
304 double dfPrevRes = 0.0;
305 dfRes = 0.0;
306 for( ; nZoomLevel < static_cast<int>(tmList.size()); nZoomLevel++ )
307 {
308 dfRes = tmList[nZoomLevel].mResX * tmList[0].mTileWidth / nBlockSize;
309 if( dfComputedRes > dfRes || fabs( dfComputedRes - dfRes ) / dfRes <= 1e-8 )
310 break;
311 dfPrevRes = dfRes;
312 }
313 if( nZoomLevel == static_cast<int>(tmList.size()) )
314 {
315 CPLError(CE_Failure, CPLE_AppDefined,
316 "Could not find an appropriate zoom level");
317 return false;
318 }
319
320 if( nZoomLevel > 0 && fabs( dfComputedRes - dfRes ) / dfRes > 1e-8 )
321 {
322 const char* pszZoomLevelStrategy = CSLFetchNameValueDef(papszOptions,
323 "ZOOM_LEVEL_STRATEGY",
324 "AUTO");
325 if( EQUAL(pszZoomLevelStrategy, "LOWER") )
326 {
327 nZoomLevel --;
328 }
329 else if( EQUAL(pszZoomLevelStrategy, "UPPER") )
330 {
331 /* do nothing */
332 }
333 else
334 {
335 if( dfPrevRes / dfComputedRes < dfComputedRes / dfRes )
336 nZoomLevel --;
337 }
338 dfRes = tmList[nZoomLevel].mResX * tmList[0].mTileWidth / nBlockSize;
339 }
340
341 CPLDebug("COG", "Using ZOOM_LEVEL %d", nZoomLevel);
342
343 const double dfTileExtent = dfRes * nBlockSize;
344 int nTLTileX = static_cast<int>(std::floor((dfMinX - dfOriX) / dfTileExtent + 1e-10));
345 int nTLTileY = static_cast<int>(std::floor((dfOriY - dfMaxY) / dfTileExtent + 1e-10));
346 int nBRTileX = static_cast<int>(std::ceil((dfMaxX - dfOriX) / dfTileExtent - 1e-10));
347 int nBRTileY = static_cast<int>(std::ceil((dfOriY - dfMinY) / dfTileExtent - 1e-10));
348
349 nAlignedLevels = std::min(std::min(10, atoi(
350 CSLFetchNameValueDef(papszOptions, "ALIGNED_LEVELS", "0"))), nZoomLevel);
351 int nAccDivisor = 1;
352 for( int i = 0; i < nAlignedLevels - 1; i++ )
353 {
354 const int nCurLevel = nZoomLevel - i;
355 const double dfResRatio =
356 tmList[nCurLevel-1].mResX / tmList[nCurLevel].mResX;
357 // Magical number that has a great number of divisors
358 // For example if previous scale denom was 50K and current one
359 // is 20K, then dfResRatio = 2.5 and dfScaledInvResRatio = 24
360 // We must then simplify 60 / 24 as 5 / 2, and make sure to
361 // align tile coordinates on multiple of the 5 numerator
362 constexpr int MAGICAL = 60;
363 const double dfScaledInvResRatio = MAGICAL / dfResRatio;
364 if( dfScaledInvResRatio < 1 || dfScaledInvResRatio > 60 ||
365 std::abs(std::round(dfScaledInvResRatio) - dfScaledInvResRatio) > 1e-10 )
366 {
367 CPLError(CE_Failure, CPLE_AppDefined,
368 "Unsupported ratio of resolution for "
369 "ALIGNED_LEVELS between zoom level %d and %d = %g",
370 nCurLevel-1, nCurLevel, dfResRatio);
371 return false;
372 }
373 const int nScaledInvResRatio = static_cast<int>(
374 std::round(dfScaledInvResRatio));
375 int nNumerator = 0;
376 for( int nDivisor = nScaledInvResRatio; nDivisor >= 2; --nDivisor )
377 {
378 if( (MAGICAL % nDivisor) == 0 && (nScaledInvResRatio % nDivisor) == 0 )
379 {
380 nNumerator = MAGICAL / nDivisor;
381 break;
382 }
383 }
384 if( nNumerator == 0 )
385 {
386 CPLError(CE_Failure, CPLE_AppDefined,
387 "Unsupported ratio of resolution for "
388 "ALIGNED_LEVELS between zoom level %d and %d = %g",
389 nCurLevel-1, nCurLevel, dfResRatio);
390 return false;
391 }
392 nAccDivisor *= nNumerator;
393 }
394 if( nAccDivisor > 1 )
395 {
396 nTLTileX = (nTLTileX / nAccDivisor) * nAccDivisor;
397 nTLTileY = (nTLTileY / nAccDivisor) * nAccDivisor;
398 nBRTileY = ((nBRTileY + nAccDivisor - 1) / nAccDivisor) * nAccDivisor;
399 nBRTileX = ((nBRTileX + nAccDivisor - 1) / nAccDivisor) * nAccDivisor;
400 }
401
402 if( nTLTileX < 0 || nTLTileY < 0 ||
403 nBRTileX > tmList[nZoomLevel].mMatrixWidth ||
404 nBRTileY > tmList[nZoomLevel].mMatrixHeight )
405 {
406 CPLError(CE_Warning, CPLE_AppDefined,
407 "Raster extent partially outside of tile matrix "
408 "bounding box. Clamping it to it");
409 }
410 nTLTileX = std::max(0, nTLTileX);
411 nTLTileY = std::max(0, nTLTileY);
412 nBRTileX = std::min(tmList[nZoomLevel].mMatrixWidth, nBRTileX);
413 nBRTileY = std::min(tmList[nZoomLevel].mMatrixHeight, nBRTileY);
414
415 dfMinX = dfOriX + nTLTileX * dfTileExtent;
416 dfMinY = dfOriY - nBRTileY * dfTileExtent;
417 dfMaxX = dfOriX + nBRTileX * dfTileExtent;
418 dfMaxY = dfOriY - nTLTileY * dfTileExtent;
419 nXSize = static_cast<int>(std::round((dfMaxX - dfMinX) / dfRes));
420 nYSize = static_cast<int>(std::round((dfMaxY - dfMinY) / dfRes));
421 }
422 else if( !osExtent.empty() || !osRes.empty() )
423 {
424 CPLStringList aosTokens(CSLTokenizeString2(osExtent, ",", 0));
425 if( aosTokens.size() != 4 )
426 {
427 CPLError(CE_Failure, CPLE_AppDefined,
428 "Invalid value for EXTENT");
429 return false;
430 }
431 dfMinX = CPLAtof(aosTokens[0]);
432 dfMinY = CPLAtof(aosTokens[1]);
433 dfMaxX = CPLAtof(aosTokens[2]);
434 dfMaxY = CPLAtof(aosTokens[3]);
435 if( !osRes.empty() )
436 dfRes = CPLAtof(osRes);
437 }
438
439 nXSize = static_cast<int>(std::round((dfMaxX - dfMinX) / dfRes));
440 nYSize = static_cast<int>(std::round((dfMaxY - dfMinY) / dfRes));
441
442 osResampling = CSLFetchNameValueDef(papszOptions,
443 "WARP_RESAMPLING",
444 CSLFetchNameValueDef(papszOptions,
445 "RESAMPLING", GetResampling(poSrcDS)));
446
447 return true;
448 }
449
450 // Used by gdalwarp
COGGetWarpingCharacteristics(GDALDataset * poSrcDS,const char * const * papszOptions,CPLString & osResampling,CPLString & osTargetSRS,int & nXSize,int & nYSize,double & dfMinX,double & dfMinY,double & dfMaxX,double & dfMaxY)451 bool COGGetWarpingCharacteristics(GDALDataset* poSrcDS,
452 const char * const* papszOptions,
453 CPLString& osResampling,
454 CPLString& osTargetSRS,
455 int& nXSize,
456 int& nYSize,
457 double& dfMinX,
458 double& dfMinY,
459 double& dfMaxX,
460 double& dfMaxY)
461 {
462 std::unique_ptr<gdal::TileMatrixSet> poTM;
463 int nZoomLevel = 0;
464 int nAlignedLevels = 0;
465 return COGGetWarpingCharacteristics(poSrcDS,
466 papszOptions,
467 osResampling,
468 osTargetSRS,
469 nXSize, nYSize,
470 dfMinX, dfMinY, dfMaxX, dfMaxY,
471 poTM, nZoomLevel, nAlignedLevels);
472 }
473
474 /************************************************************************/
475 /* COGHasWarpingOptions() */
476 /************************************************************************/
477
COGHasWarpingOptions(CSLConstList papszOptions)478 bool COGHasWarpingOptions(CSLConstList papszOptions)
479 {
480 return CSLFetchNameValue(papszOptions, "TARGET_SRS") != nullptr ||
481 !EQUAL(CSLFetchNameValueDef(papszOptions, "TILING_SCHEME", "CUSTOM"),
482 "CUSTOM");
483 }
484
485 /************************************************************************/
486 /* COGRemoveWarpingOptions() */
487 /************************************************************************/
488
COGRemoveWarpingOptions(CPLStringList & aosOptions)489 void COGRemoveWarpingOptions(CPLStringList& aosOptions)
490 {
491 aosOptions.SetNameValue("TARGET_SRS", nullptr);
492 aosOptions.SetNameValue("TILING_SCHEME", nullptr);
493 aosOptions.SetNameValue("EXTENT", nullptr);
494 aosOptions.SetNameValue("RES", nullptr);
495 aosOptions.SetNameValue("ALIGNED_LEVELS", nullptr);
496 aosOptions.SetNameValue("ZOOM_LEVEL_STRATEGY", nullptr);
497 }
498
499 /************************************************************************/
500 /* CreateReprojectedDS() */
501 /************************************************************************/
502
CreateReprojectedDS(const char * pszDstFilename,GDALDataset * poSrcDS,const char * const * papszOptions,const CPLString & osResampling,const CPLString & osTargetSRS,const int nXSize,const int nYSize,const double dfMinX,const double dfMinY,const double dfMaxX,const double dfMaxY,GDALProgressFunc pfnProgress,void * pProgressData,double & dfCurPixels,double & dfTotalPixelsToProcess)503 static std::unique_ptr<GDALDataset> CreateReprojectedDS(
504 const char* pszDstFilename,
505 GDALDataset *poSrcDS,
506 const char * const* papszOptions,
507 const CPLString& osResampling,
508 const CPLString& osTargetSRS,
509 const int nXSize,
510 const int nYSize,
511 const double dfMinX,
512 const double dfMinY,
513 const double dfMaxX,
514 const double dfMaxY,
515 GDALProgressFunc pfnProgress,
516 void * pProgressData,
517 double& dfCurPixels,
518 double& dfTotalPixelsToProcess)
519 {
520 char** papszArg = nullptr;
521 // We could have done a warped VRT, but overview building on it might be
522 // slow, so materialize as GTiff
523 papszArg = CSLAddString(papszArg, "-of");
524 papszArg = CSLAddString(papszArg, "GTiff");
525 papszArg = CSLAddString(papszArg, "-co");
526 papszArg = CSLAddString(papszArg, "TILED=YES");
527 papszArg = CSLAddString(papszArg, "-co");
528 papszArg = CSLAddString(papszArg, "SPARSE_OK=YES");
529 const char* pszBIGTIFF = CSLFetchNameValue(papszOptions, "BIGTIFF");
530 if( pszBIGTIFF )
531 {
532 papszArg = CSLAddString(papszArg, "-co");
533 papszArg = CSLAddString(papszArg, (CPLString("BIGTIFF=") + pszBIGTIFF).c_str());
534 }
535 papszArg = CSLAddString(papszArg, "-co");
536 papszArg = CSLAddString(papszArg,
537 HasZSTDCompression() ? "COMPRESS=ZSTD" : "COMPRESS=LZW");
538 papszArg = CSLAddString(papszArg, "-t_srs");
539 papszArg = CSLAddString(papszArg, osTargetSRS);
540 papszArg = CSLAddString(papszArg, "-te");
541 papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMinX));
542 papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMinY));
543 papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMaxX));
544 papszArg = CSLAddString(papszArg, CPLSPrintf("%.18g,", dfMaxY));
545 papszArg = CSLAddString(papszArg, "-ts");
546 papszArg = CSLAddString(papszArg, CPLSPrintf("%d", nXSize));
547 papszArg = CSLAddString(papszArg, CPLSPrintf("%d", nYSize));
548 int bHasNoData = FALSE;
549 poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoData);
550 if( !bHasNoData && CPLTestBool(CSLFetchNameValueDef(
551 papszOptions, "ADD_ALPHA", "YES")) )
552 {
553 papszArg = CSLAddString(papszArg, "-dstalpha");
554 }
555 papszArg = CSLAddString(papszArg, "-r");
556 papszArg = CSLAddString(papszArg, osResampling);
557 papszArg = CSLAddString(papszArg, "-wo");
558 papszArg = CSLAddString(papszArg, "SAMPLE_GRID=YES");
559 const char* pszNumThreads = CSLFetchNameValue(papszOptions, "NUM_THREADS");
560 if( pszNumThreads )
561 {
562 papszArg = CSLAddString(papszArg, "-wo");
563 papszArg = CSLAddString(papszArg, (CPLString("NUM_THREADS=") + pszNumThreads).c_str());
564 }
565
566 const auto poFirstBand = poSrcDS->GetRasterBand(1);
567 const bool bHasMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
568
569 const int nBands = poSrcDS->GetRasterCount();
570 const char* pszOverviews = CSLFetchNameValueDef(
571 papszOptions, "OVERVIEWS", "AUTO");
572 const bool bRecreateOvr = EQUAL(pszOverviews, "FORCE_USE_EXISTING") ||
573 EQUAL(pszOverviews, "NONE");
574 dfTotalPixelsToProcess =
575 double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) +
576 ((bHasMask && !bRecreateOvr) ? double(nXSize) * nYSize / 3 : 0) +
577 (!bRecreateOvr ? double(nXSize) * nYSize * nBands / 3: 0) +
578 double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) * 4. / 3;
579
580 auto psOptions = GDALWarpAppOptionsNew(papszArg, nullptr);
581 CSLDestroy(papszArg);
582 if( psOptions == nullptr )
583 return nullptr;
584
585 const double dfNextPixels =
586 double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0));
587 void* pScaledProgress = GDALCreateScaledProgress(
588 dfCurPixels / dfTotalPixelsToProcess,
589 dfNextPixels / dfTotalPixelsToProcess,
590 pfnProgress, pProgressData );
591 dfCurPixels = dfNextPixels;
592
593 CPLDebug("COG", "Reprojecting source dataset: start");
594 GDALWarpAppOptionsSetProgress(psOptions, GDALScaledProgress, pScaledProgress );
595 CPLString osTmpFile(GetTmpFilename(pszDstFilename, "warped.tif.tmp"));
596 auto hSrcDS = GDALDataset::ToHandle(poSrcDS);
597 auto hRet = GDALWarp( osTmpFile, nullptr,
598 1, &hSrcDS,
599 psOptions, nullptr);
600 GDALWarpAppOptionsFree(psOptions);
601 CPLDebug("COG", "Reprojecting source dataset: end");
602
603 GDALDestroyScaledProgress(pScaledProgress);
604
605 return std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(hRet));
606 }
607
608 /************************************************************************/
609 /* GDALCOGCreator */
610 /************************************************************************/
611
612 struct GDALCOGCreator final
613 {
614 std::unique_ptr<GDALDataset> m_poReprojectedDS{};
615 std::unique_ptr<GDALDataset> m_poRGBMaskDS{};
616 CPLString m_osTmpOverviewFilename{};
617 CPLString m_osTmpMskOverviewFilename{};
618
619 ~GDALCOGCreator();
620
621 GDALDataset* Create(const char * pszFilename,
622 GDALDataset * const poSrcDS,
623 char ** papszOptions,
624 GDALProgressFunc pfnProgress,
625 void * pProgressData );
626 };
627
628 /************************************************************************/
629 /* GDALCOGCreator::~GDALCOGCreator() */
630 /************************************************************************/
631
~GDALCOGCreator()632 GDALCOGCreator::~GDALCOGCreator()
633 {
634 if( m_poReprojectedDS )
635 {
636 CPLString osProjectedDSName(m_poReprojectedDS->GetDescription());
637 // Destroy m_poRGBMaskDS before m_poReprojectedDS since the former
638 // references the later
639 m_poRGBMaskDS.reset();
640 m_poReprojectedDS.reset();
641 VSIUnlink(osProjectedDSName);
642 }
643 if( !m_osTmpOverviewFilename.empty() )
644 {
645 VSIUnlink(m_osTmpOverviewFilename);
646 }
647 if( !m_osTmpMskOverviewFilename.empty() )
648 {
649 VSIUnlink(m_osTmpMskOverviewFilename);
650 }
651 }
652
653 /************************************************************************/
654 /* GDALCOGCreator::Create() */
655 /************************************************************************/
656
Create(const char * pszFilename,GDALDataset * const poSrcDS,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)657 GDALDataset* GDALCOGCreator::Create(const char * pszFilename,
658 GDALDataset * const poSrcDS,
659 char ** papszOptions,
660 GDALProgressFunc pfnProgress,
661 void * pProgressData )
662 {
663 if( pfnProgress == nullptr )
664 pfnProgress = GDALDummyProgress;
665
666 if( poSrcDS->GetRasterCount() == 0 )
667 {
668 CPLError(CE_Failure, CPLE_NotSupported,
669 "COG driver does not support 0-band source raster");
670 return nullptr;
671 }
672
673 CPLConfigOptionSetter oSetterReportDirtyBlockFlushing(
674 "GDAL_REPORT_DIRTY_BLOCK_FLUSHING", "NO", true);
675
676 double dfCurPixels = 0;
677 double dfTotalPixelsToProcess = 0;
678 GDALDataset* poCurDS = poSrcDS;
679
680 std::unique_ptr<gdal::TileMatrixSet> poTM;
681 int nZoomLevel = 0;
682 int nAlignedLevels = 0;
683 if( COGHasWarpingOptions(papszOptions) )
684 {
685 CPLString osTargetResampling;
686 CPLString osTargetSRS;
687 int nTargetXSize = 0;
688 int nTargetYSize = 0;
689 double dfTargetMinX = 0;
690 double dfTargetMinY = 0;
691 double dfTargetMaxX = 0;
692 double dfTargetMaxY = 0;
693 if( !COGGetWarpingCharacteristics(poCurDS, papszOptions,
694 osTargetResampling,
695 osTargetSRS,
696 nTargetXSize, nTargetYSize,
697 dfTargetMinX, dfTargetMinY,
698 dfTargetMaxX, dfTargetMaxY,
699 poTM, nZoomLevel, nAlignedLevels) )
700 {
701 return nullptr;
702 }
703
704 // Collect information on source dataset to see if it already
705 // matches the warping specifications
706 CPLString osSrcSRS;
707 const auto poSrcSRS = poCurDS->GetSpatialRef();
708 if( poSrcSRS )
709 {
710 const char* pszAuthName = poSrcSRS->GetAuthorityName(nullptr);
711 const char* pszAuthCode = poSrcSRS->GetAuthorityCode(nullptr);
712 if( pszAuthName && pszAuthCode )
713 {
714 osSrcSRS = pszAuthName;
715 osSrcSRS += ':';
716 osSrcSRS += pszAuthCode;
717 }
718 }
719 double dfSrcMinX = 0;
720 double dfSrcMinY = 0;
721 double dfSrcMaxX = 0;
722 double dfSrcMaxY = 0;
723 double adfSrcGT[6];
724 const int nSrcXSize = poCurDS->GetRasterXSize();
725 const int nSrcYSize = poCurDS->GetRasterYSize();
726 if( poCurDS->GetGeoTransform(adfSrcGT) == CE_None )
727 {
728 dfSrcMinX = adfSrcGT[0];
729 dfSrcMaxY = adfSrcGT[3];
730 dfSrcMaxX = adfSrcGT[0] + nSrcXSize * adfSrcGT[1];
731 dfSrcMinY = adfSrcGT[3] + nSrcYSize * adfSrcGT[5];
732 }
733
734 if( nTargetXSize == nSrcXSize &&
735 nTargetYSize == nSrcYSize &&
736 osTargetSRS == osSrcSRS &&
737 fabs(dfSrcMinX - dfTargetMinX) < 1e-10 * fabs(dfSrcMinX) &&
738 fabs(dfSrcMinY - dfTargetMinY) < 1e-10 * fabs(dfSrcMinY) &&
739 fabs(dfSrcMaxX - dfTargetMaxX) < 1e-10 * fabs(dfSrcMaxX) &&
740 fabs(dfSrcMaxY - dfTargetMaxY) < 1e-10 * fabs(dfSrcMaxY) )
741 {
742 CPLDebug("COG", "Skipping reprojection step: "
743 "source dataset matches reprojection specifications");
744 }
745 else
746 {
747 m_poReprojectedDS =
748 CreateReprojectedDS(pszFilename, poCurDS,
749 papszOptions,
750 osTargetResampling,
751 osTargetSRS,
752 nTargetXSize, nTargetYSize,
753 dfTargetMinX, dfTargetMinY,
754 dfTargetMaxX, dfTargetMaxY,
755 pfnProgress, pProgressData,
756 dfCurPixels, dfTotalPixelsToProcess);
757 if( !m_poReprojectedDS )
758 return nullptr;
759 poCurDS = m_poReprojectedDS.get();
760 }
761 }
762
763 CPLString osCompress = CSLFetchNameValueDef(papszOptions, "COMPRESS", "NONE");
764 if( EQUAL(osCompress, "JPEG") &&
765 poCurDS->GetRasterCount() == 4 )
766 {
767 char** papszArg = nullptr;
768 papszArg = CSLAddString(papszArg, "-of");
769 papszArg = CSLAddString(papszArg, "VRT");
770 papszArg = CSLAddString(papszArg, "-b");
771 papszArg = CSLAddString(papszArg, "1");
772 papszArg = CSLAddString(papszArg, "-b");
773 papszArg = CSLAddString(papszArg, "2");
774 papszArg = CSLAddString(papszArg, "-b");
775 papszArg = CSLAddString(papszArg, "3");
776 papszArg = CSLAddString(papszArg, "-mask");
777 papszArg = CSLAddString(papszArg, "4");
778 GDALTranslateOptions* psOptions = GDALTranslateOptionsNew(papszArg, nullptr);
779 CSLDestroy(papszArg);
780 GDALDatasetH hRGBMaskDS = GDALTranslate("",
781 GDALDataset::ToHandle(poCurDS),
782 psOptions,
783 nullptr);
784 GDALTranslateOptionsFree(psOptions);
785 if( !hRGBMaskDS )
786 {
787 return nullptr;
788 }
789 m_poRGBMaskDS.reset( GDALDataset::FromHandle(hRGBMaskDS) );
790 poCurDS = m_poRGBMaskDS.get();
791 }
792
793 const int nBands = poCurDS->GetRasterCount();
794 const int nXSize = poCurDS->GetRasterXSize();
795 const int nYSize = poCurDS->GetRasterYSize();
796
797 CPLString osBlockSize(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", ""));
798 if( osBlockSize.empty() )
799 {
800 if( poTM )
801 {
802 osBlockSize.Printf("%d", poTM->tileMatrixList()[0].mTileWidth);
803 }
804 else
805 {
806 osBlockSize = "512";
807 }
808 }
809
810 const int nOvrThresholdSize = atoi(osBlockSize);
811
812 const auto poFirstBand = poCurDS->GetRasterBand(1);
813 const bool bHasMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
814
815 CPLString osOverviews = CSLFetchNameValueDef(
816 papszOptions, "OVERVIEWS", "AUTO");
817 const bool bRecreateOvr = EQUAL(osOverviews, "FORCE_USE_EXISTING") ||
818 EQUAL(osOverviews, "NONE");
819 const bool bGenerateMskOvr =
820 !bRecreateOvr &&
821 bHasMask &&
822 (nXSize > nOvrThresholdSize || nYSize > nOvrThresholdSize) &&
823 (EQUAL(osOverviews, "IGNORE_EXISTING") ||
824 poFirstBand->GetMaskBand()->GetOverviewCount() == 0);
825 const bool bGenerateOvr =
826 !bRecreateOvr &&
827 (nXSize > nOvrThresholdSize || nYSize > nOvrThresholdSize) &&
828 (EQUAL(osOverviews, "IGNORE_EXISTING") ||
829 poFirstBand->GetOverviewCount() == 0);
830
831 std::vector<std::pair<int, int>> asOverviewDims;
832 int nTmpXSize = nXSize;
833 int nTmpYSize = nYSize;
834 if( poTM )
835 {
836 const auto& tmList = poTM->tileMatrixList();
837 int nCurLevel = nZoomLevel;
838 while( nTmpXSize > nOvrThresholdSize || nTmpYSize > nOvrThresholdSize )
839 {
840 const double dfResRatio = (nCurLevel >= 1) ?
841 tmList[nCurLevel-1].mResX / tmList[nCurLevel].mResX : 2;
842 nTmpXSize = static_cast<int>(nTmpXSize / dfResRatio + 0.5);
843 nTmpYSize = static_cast<int>(nTmpYSize / dfResRatio + 0.5);
844 asOverviewDims.push_back(std::pair<int,int>(nTmpXSize, nTmpYSize));
845 nCurLevel --;
846 }
847 }
848 else
849 {
850 while( nTmpXSize > nOvrThresholdSize || nTmpYSize > nOvrThresholdSize )
851 {
852 nTmpXSize /= 2;
853 nTmpYSize /= 2;
854 asOverviewDims.push_back(std::pair<int,int>(nTmpXSize, nTmpYSize));
855 }
856 }
857
858 if( dfTotalPixelsToProcess == 0.0 )
859 {
860 dfTotalPixelsToProcess =
861 (bGenerateMskOvr ? double(nXSize) * nYSize / 3 : 0) +
862 (bGenerateOvr ? double(nXSize) * nYSize * nBands / 3: 0) +
863 double(nXSize) * nYSize * (nBands + (bHasMask ? 1 : 0)) * 4. / 3;
864 }
865
866 CPLStringList aosOverviewOptions;
867 aosOverviewOptions.SetNameValue("COMPRESS",
868 CPLGetConfigOption("COG_TMP_COMPRESSION", // only for debug purposes
869 HasZSTDCompression() ? "ZSTD" : "LZW"));
870 aosOverviewOptions.SetNameValue("NUM_THREADS",
871 CSLFetchNameValue(papszOptions, "NUM_THREADS"));
872 aosOverviewOptions.SetNameValue("BIGTIFF", "YES");
873 aosOverviewOptions.SetNameValue("SPARSE_OK", "YES");
874
875 if( bGenerateMskOvr )
876 {
877 CPLDebug("COG", "Generating overviews of the mask: start");
878 m_osTmpMskOverviewFilename = GetTmpFilename(pszFilename, "msk.ovr.tmp");
879 GDALRasterBand* poSrcMask = poFirstBand->GetMaskBand();
880 const char* pszResampling = CSLFetchNameValueDef(papszOptions,
881 "OVERVIEW_RESAMPLING",
882 CSLFetchNameValueDef(papszOptions,
883 "RESAMPLING", GetResampling(poSrcDS)));
884
885 double dfNextPixels = dfCurPixels + double(nXSize) * nYSize / 3;
886 void* pScaledProgress = GDALCreateScaledProgress(
887 dfCurPixels / dfTotalPixelsToProcess,
888 dfNextPixels / dfTotalPixelsToProcess,
889 pfnProgress, pProgressData );
890 dfCurPixels = dfNextPixels;
891
892 // Used by GDALRegenerateOverviews() and GDALRegenerateOverviewsMultiBand()
893 CPLConfigOptionSetter oSetterRegeneratedBandIsMask(
894 "GDAL_REGENERATED_BAND_IS_MASK", "YES", true);
895
896 CPLErr eErr = GTIFFBuildOverviewsEx(
897 m_osTmpMskOverviewFilename,
898 1, &poSrcMask,
899 static_cast<int>(asOverviewDims.size()),
900 nullptr, asOverviewDims.data(),
901 pszResampling,
902 aosOverviewOptions.List(),
903 GDALScaledProgress, pScaledProgress );
904 CPLDebug("COG", "Generating overviews of the mask: end");
905
906 GDALDestroyScaledProgress(pScaledProgress);
907 if( eErr != CE_None )
908 {
909 return nullptr;
910 }
911 }
912
913 if( bGenerateOvr )
914 {
915 CPLDebug("COG", "Generating overviews of the imagery: start");
916 m_osTmpOverviewFilename = GetTmpFilename(pszFilename, "ovr.tmp");
917 std::vector<GDALRasterBand*> apoSrcBands;
918 for( int i = 0; i < nBands; i++ )
919 apoSrcBands.push_back( poCurDS->GetRasterBand(i+1) );
920 const char* pszResampling = CSLFetchNameValueDef(papszOptions,
921 "OVERVIEW_RESAMPLING",
922 CSLFetchNameValueDef(papszOptions,
923 "RESAMPLING", GetResampling(poSrcDS)));
924
925 double dfNextPixels = dfCurPixels + double(nXSize) * nYSize * nBands / 3;
926 void* pScaledProgress = GDALCreateScaledProgress(
927 dfCurPixels / dfTotalPixelsToProcess,
928 dfNextPixels / dfTotalPixelsToProcess,
929 pfnProgress, pProgressData );
930 dfCurPixels = dfNextPixels;
931
932 if( nBands > 1 )
933 {
934 aosOverviewOptions.SetNameValue("INTERLEAVE", "PIXEL");
935 }
936 if( !m_osTmpMskOverviewFilename.empty() )
937 {
938 aosOverviewOptions.SetNameValue("MASK_OVERVIEW_DATASET",
939 m_osTmpMskOverviewFilename);
940 }
941 CPLErr eErr = GTIFFBuildOverviewsEx(
942 m_osTmpOverviewFilename,
943 nBands, &apoSrcBands[0],
944 static_cast<int>(asOverviewDims.size()),
945 nullptr, asOverviewDims.data(),
946 pszResampling,
947 aosOverviewOptions.List(),
948 GDALScaledProgress, pScaledProgress );
949 CPLDebug("COG", "Generating overviews of the imagery: end");
950
951 GDALDestroyScaledProgress(pScaledProgress);
952 if( eErr != CE_None )
953 {
954 return nullptr;
955 }
956 }
957
958 CPLStringList aosOptions;
959 aosOptions.SetNameValue("COPY_SRC_OVERVIEWS", "YES");
960 aosOptions.SetNameValue("COMPRESS", osCompress);
961 aosOptions.SetNameValue("TILED", "YES");
962 aosOptions.SetNameValue("BLOCKXSIZE", osBlockSize);
963 aosOptions.SetNameValue("BLOCKYSIZE", osBlockSize);
964 const char* pszPredictor = CSLFetchNameValueDef(papszOptions, "PREDICTOR", "FALSE");
965 const char* pszPredictorValue = GetPredictor(poSrcDS, pszPredictor);
966 if (pszPredictorValue != nullptr)
967 {
968 aosOptions.SetNameValue("PREDICTOR", pszPredictorValue);
969 }
970
971 const char* pszQuality = CSLFetchNameValue(papszOptions, "QUALITY");
972 if( EQUAL(osCompress, "JPEG") )
973 {
974 aosOptions.SetNameValue("JPEG_QUALITY", pszQuality);
975 if( nBands == 3 )
976 aosOptions.SetNameValue("PHOTOMETRIC", "YCBCR");
977 }
978 else if( EQUAL(osCompress, "WEBP") )
979 {
980 if( pszQuality && atoi(pszQuality) == 100 )
981 aosOptions.SetNameValue("WEBP_LOSSLESS", "YES");
982 aosOptions.SetNameValue("WEBP_LEVEL", pszQuality);
983 }
984 else if( EQUAL(osCompress, "DEFLATE") || EQUAL(osCompress, "LERC_DEFLATE") )
985 {
986 aosOptions.SetNameValue("ZLEVEL",
987 CSLFetchNameValue(papszOptions, "LEVEL"));
988 }
989 else if( EQUAL(osCompress, "ZSTD") || EQUAL(osCompress, "LERC_ZSTD") )
990 {
991 aosOptions.SetNameValue("ZSTD_LEVEL",
992 CSLFetchNameValue(papszOptions, "LEVEL"));
993 }
994
995 if( STARTS_WITH_CI(osCompress, "LERC") )
996 {
997 aosOptions.SetNameValue("MAX_Z_ERROR",
998 CSLFetchNameValue(papszOptions, "MAX_Z_ERROR"));
999 }
1000 aosOptions.SetNameValue("BIGTIFF",
1001 CSLFetchNameValue(papszOptions, "BIGTIFF"));
1002 aosOptions.SetNameValue("NUM_THREADS",
1003 CSLFetchNameValue(papszOptions, "NUM_THREADS"));
1004 aosOptions.SetNameValue("GEOTIFF_VERSION",
1005 CSLFetchNameValue(papszOptions, "GEOTIFF_VERSION"));
1006 aosOptions.SetNameValue("SPARSE_OK",
1007 CSLFetchNameValue(papszOptions, "SPARSE_OK"));
1008
1009 if( EQUAL( osOverviews, "NONE") )
1010 {
1011 aosOptions.SetNameValue("@OVERVIEW_DATASET", "");
1012 }
1013 else
1014 {
1015 if( !m_osTmpOverviewFilename.empty() )
1016 {
1017 aosOptions.SetNameValue("@OVERVIEW_DATASET", m_osTmpOverviewFilename);
1018 }
1019 if( !m_osTmpMskOverviewFilename.empty() )
1020 {
1021 aosOptions.SetNameValue("@MASK_OVERVIEW_DATASET", m_osTmpMskOverviewFilename);
1022 }
1023 }
1024
1025 const CPLString osTilingScheme(CSLFetchNameValueDef(papszOptions,
1026 "TILING_SCHEME", "CUSTOM"));
1027 if( osTilingScheme != "CUSTOM" )
1028 {
1029 aosOptions.SetNameValue("@TILING_SCHEME_NAME", osTilingScheme);
1030 aosOptions.SetNameValue("@TILING_SCHEME_ZOOM_LEVEL",
1031 CPLSPrintf("%d", nZoomLevel));
1032 if( nAlignedLevels > 0 )
1033 {
1034 aosOptions.SetNameValue("@TILING_SCHEME_ALIGNED_LEVELS",
1035 CPLSPrintf("%d", nAlignedLevels));
1036 }
1037 }
1038 const char* pszOverviewCompress = CSLFetchNameValue(papszOptions, "OVERVIEW_COMPRESS");
1039
1040 CPLConfigOptionSetter ovrCompressSetter("COMPRESS_OVERVIEW", pszOverviewCompress, true);
1041 CPLConfigOptionSetter ovrQualityJpegSetter("JPEG_QUALITY_OVERVIEW", CSLFetchNameValue(papszOptions, "OVERVIEW_QUALITY"), true);
1042 CPLConfigOptionSetter ovrQualityWebpSetter("WEBP_LEVEL_OVERVIEW", CSLFetchNameValue(papszOptions, "OVERVIEW_QUALITY"), true);
1043
1044 std::unique_ptr<CPLConfigOptionSetter> poPhotometricSetter;
1045 if (pszOverviewCompress != nullptr && nBands == 3 && EQUAL(pszOverviewCompress, "JPEG") )
1046 {
1047 poPhotometricSetter.reset(new CPLConfigOptionSetter("PHOTOMETRIC_OVERVIEW", "YCBCR", true));
1048 }
1049
1050 const char* osOvrPredictor = CSLFetchNameValueDef(papszOptions, "OVERVIEW_PREDICTOR", "FALSE");
1051 const char* pszOvrPredictorValue = GetPredictor(poSrcDS, osOvrPredictor);
1052 CPLConfigOptionSetter ovrPredictorSetter("PREDICTOR_OVERVIEW", pszOvrPredictorValue, true);
1053
1054 GDALDriver* poGTiffDrv = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
1055 if( !poGTiffDrv )
1056 return nullptr;
1057 void* pScaledProgress = GDALCreateScaledProgress(
1058 dfCurPixels / dfTotalPixelsToProcess,
1059 1.0,
1060 pfnProgress, pProgressData );
1061
1062 CPLConfigOptionSetter oSetterInternalMask(
1063 "GDAL_TIFF_INTERNAL_MASK", "YES", false);
1064
1065 CPLDebug("COG", "Generating final product: start");
1066 auto poRet = poGTiffDrv->CreateCopy(pszFilename, poCurDS, false,
1067 aosOptions.List(),
1068 GDALScaledProgress, pScaledProgress);
1069
1070 GDALDestroyScaledProgress(pScaledProgress);
1071
1072 if( poRet )
1073 poRet->FlushCache();
1074
1075 CPLDebug("COG", "Generating final product: end");
1076 return poRet;
1077 }
1078
1079 /************************************************************************/
1080 /* COGCreateCopy() */
1081 /************************************************************************/
1082
COGCreateCopy(const char * pszFilename,GDALDataset * poSrcDS,int,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)1083 static GDALDataset* COGCreateCopy( const char * pszFilename,
1084 GDALDataset *poSrcDS,
1085 int /*bStrict*/, char ** papszOptions,
1086 GDALProgressFunc pfnProgress,
1087 void * pProgressData )
1088 {
1089 return GDALCOGCreator().Create(pszFilename, poSrcDS,
1090 papszOptions, pfnProgress, pProgressData);
1091 }
1092
1093 /************************************************************************/
1094 /* GDALRegister_COG() */
1095 /************************************************************************/
1096
1097 class GDALCOGDriver final: public GDALDriver
1098 {
1099 bool m_bInitialized = false;
1100
1101 bool bHasLZW = false;
1102 bool bHasDEFLATE = false;
1103 bool bHasLZMA = false;
1104 bool bHasZSTD = false;
1105 bool bHasJPEG = false;
1106 bool bHasWebP = false;
1107 bool bHasLERC = false;
1108 std::string osCompressValues{};
1109
1110 void InitializeCreationOptionList();
1111
1112 public:
1113 GDALCOGDriver();
1114
GetMetadataItem(const char * pszName,const char * pszDomain)1115 const char* GetMetadataItem(const char* pszName, const char* pszDomain) override
1116 {
1117 if( EQUAL(pszName, GDAL_DMD_CREATIONOPTIONLIST) )
1118 {
1119 InitializeCreationOptionList();
1120 }
1121 return GDALDriver::GetMetadataItem(pszName, pszDomain);
1122 }
1123
GetMetadata(const char * pszDomain)1124 char** GetMetadata(const char* pszDomain) override
1125 {
1126 InitializeCreationOptionList();
1127 return GDALDriver::GetMetadata(pszDomain);
1128 }
1129 };
1130
GDALCOGDriver()1131 GDALCOGDriver::GDALCOGDriver()
1132 {
1133 // We could defer this in InitializeCreationOptionList() but with currently
1134 // released libtiff versions where there was a bug (now fixed) in
1135 // TIFFGetConfiguredCODECs(), this wouldn't work properly if the LERC codec
1136 // had been registered in between
1137 osCompressValues = GTiffGetCompressValues(
1138 bHasLZW, bHasDEFLATE, bHasLZMA, bHasZSTD, bHasJPEG, bHasWebP, bHasLERC,
1139 true /* bForCOG */);
1140 }
1141
InitializeCreationOptionList()1142 void GDALCOGDriver::InitializeCreationOptionList()
1143 {
1144 if( m_bInitialized )
1145 return;
1146 m_bInitialized = true;
1147
1148 CPLString osOptions;
1149 osOptions = "<CreationOptionList>"
1150 " <Option name='COMPRESS' type='string-select'>";
1151 osOptions += osCompressValues;
1152 osOptions += " </Option>";
1153
1154 osOptions += " <Option name='OVERVIEW_COMPRESS' type='string-select'>";
1155 osOptions += osCompressValues;
1156 osOptions += " </Option>";
1157
1158 if( bHasLZW || bHasDEFLATE || bHasZSTD )
1159 {
1160 const char* osPredictorOptions = " <Value>YES</Value>"
1161 " <Value>NO</Value>"
1162 " <Value alias='2'>STANDARD</Value>"
1163 " <Value alias='3'>FLOATING_POINT</Value>";
1164
1165 osOptions += " <Option name='LEVEL' type='int' "
1166 "description='DEFLATE/ZSTD compression level: 1 (fastest)'/>";
1167
1168 osOptions += " <Option name='PREDICTOR' type='string-select' default='FALSE'>";
1169 osOptions += osPredictorOptions;
1170 osOptions += " </Option>"
1171 " <Option name='OVERVIEW_PREDICTOR' type='string-select' default='FALSE'>";
1172 osOptions += osPredictorOptions;
1173 osOptions += " </Option>";
1174 }
1175 if( bHasJPEG || bHasWebP )
1176 {
1177 osOptions += " <Option name='QUALITY' type='int' "
1178 "description='JPEG/WEBP quality 1-100' default='75'/>"
1179 " <Option name='OVERVIEW_QUALITY' type='int' "
1180 "description='Overview JPEG/WEBP quality 1-100' default='75'/>";
1181 }
1182 if( bHasLERC )
1183 {
1184 osOptions += ""
1185 " <Option name='MAX_Z_ERROR' type='float' description='Maximum error for LERC compression' default='0'/>";
1186 }
1187 osOptions +=
1188 " <Option name='NUM_THREADS' type='string' "
1189 "description='Number of worker threads for compression. "
1190 "Can be set to ALL_CPUS' default='1'/>"
1191 " <Option name='BLOCKSIZE' type='int' "
1192 "description='Tile size in pixels' min='128' default='512'/>"
1193 " <Option name='BIGTIFF' type='string-select' description='"
1194 "Force creation of BigTIFF file'>"
1195 " <Value>YES</Value>"
1196 " <Value>NO</Value>"
1197 " <Value>IF_NEEDED</Value>"
1198 " <Value>IF_SAFER</Value>"
1199 " </Option>"
1200 " <Option name='RESAMPLING' type='string' "
1201 "description='Resampling method for overviews or warping'/>"
1202 " <Option name='OVERVIEW_RESAMPLING' type='string' "
1203 "description='Resampling method for overviews'/>"
1204 " <Option name='WARP_RESAMPLING' type='string' "
1205 "description='Resampling method for warping'/>"
1206 " <Option name='OVERVIEWS' type='string-select' description='"
1207 "Behavior regarding overviews'>"
1208 " <Value>AUTO</Value>"
1209 " <Value>IGNORE_EXISTING</Value>"
1210 " <Value>FORCE_USE_EXISTING</Value>"
1211 " <Value>NONE</Value>"
1212 " </Option>"
1213 " <Option name='TILING_SCHEME' type='string' description='"
1214 "Which tiling scheme to use pre-defined value or custom inline/outline "
1215 "JSON definition' default='CUSTOM'>"
1216 " <Value>CUSTOM</Value>";
1217
1218 const auto tmsList = gdal::TileMatrixSet::listPredefinedTileMatrixSets();
1219 for( const auto& tmsName: tmsList )
1220 {
1221 const auto poTM = gdal::TileMatrixSet::parse(tmsName.c_str());
1222 if( poTM &&
1223 poTM->haveAllLevelsSameTopLeft() &&
1224 poTM->haveAllLevelsSameTileSize() &&
1225 !poTM->hasVariableMatrixWidth() )
1226 {
1227 osOptions += " <Value>";
1228 osOptions += tmsName;
1229 osOptions += "</Value>";
1230 }
1231 }
1232
1233 osOptions +=
1234 " </Option>"
1235 " <Option name='ZOOM_LEVEL_STRATEGY' type='string-select' "
1236 "description='Strategy to determine zoom level. "
1237 "Only used for TILING_SCHEME != CUSTOM' default='AUTO'>"
1238 " <Value>AUTO</Value>"
1239 " <Value>LOWER</Value>"
1240 " <Value>UPPER</Value>"
1241 " </Option>"
1242 " <Option name='TARGET_SRS' type='string' "
1243 "description='Target SRS as EPSG:XXXX, WKT or PROJ string for reprojection'/>"
1244 " <Option name='RES' type='float' description='"
1245 "Target resolution for reprojection'/>"
1246 " <Option name='EXTENT' type='string' description='"
1247 "Target extent as minx,miny,maxx,maxy for reprojection'/>"
1248 " <Option name='ALIGNED_LEVELS' type='int' description='"
1249 "Number of overview levels for which the tiles from GeoTIFF and the "
1250 "specified tiling scheme match'/>"
1251 " <Option name='ADD_ALPHA' type='boolean' description='Can be set to NO to "
1252 "disable the addition of an alpha band in case of reprojection' default='YES'/>"
1253 #if LIBGEOTIFF_VERSION >= 1600
1254 " <Option name='GEOTIFF_VERSION' type='string-select' default='AUTO' description='Which version of GeoTIFF must be used'>"
1255 " <Value>AUTO</Value>"
1256 " <Value>1.0</Value>"
1257 " <Value>1.1</Value>"
1258 " </Option>"
1259 #endif
1260 " <Option name='SPARSE_OK' type='boolean' description='Should empty blocks be omitted on disk?' default='FALSE'/>"
1261 "</CreationOptionList>";
1262
1263 SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST, osOptions.c_str());
1264 }
1265
GDALRegister_COG()1266 void GDALRegister_COG()
1267
1268 {
1269 if( GDALGetDriverByName( "COG" ) != nullptr )
1270 return;
1271
1272 auto poDriver = new GDALCOGDriver();
1273 poDriver->SetDescription( "COG" );
1274 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
1275 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Cloud optimized GeoTIFF generator" );
1276 poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/cog.html" );
1277
1278 poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1279 "Byte UInt16 Int16 UInt32 Int32 Float32 "
1280 "Float64 CInt16 CInt32 CFloat32 CFloat64" );
1281
1282 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1283
1284 poDriver->pfnCreateCopy = COGCreateCopy;
1285
1286 GetGDALDriverManager()->RegisterDriver( poDriver );
1287 }
1288