Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster tile" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_tile.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "cpl_mem_cache.h"
17 : #include "cpl_worker_thread_pool.h"
18 : #include "gdal_alg_priv.h"
19 : #include "gdal_priv.h"
20 : #include "gdalwarper.h"
21 : #include "gdal_utils.h"
22 : #include "ogr_spatialref.h"
23 : #include "memdataset.h"
24 : #include "tilematrixset.hpp"
25 :
26 : #include <algorithm>
27 : #include <array>
28 : #include <atomic>
29 : #include <cmath>
30 : #include <mutex>
31 :
32 : //! @cond Doxygen_Suppress
33 :
34 : #ifndef _
35 : #define _(x) (x)
36 : #endif
37 :
38 : /************************************************************************/
39 : /* GDALRasterTileAlgorithm::GDALRasterTileAlgorithm() */
40 : /************************************************************************/
41 :
42 98 : GDALRasterTileAlgorithm::GDALRasterTileAlgorithm()
43 98 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
44 : {
45 98 : AddProgressArg();
46 98 : AddOpenOptionsArg(&m_openOptions);
47 98 : AddInputFormatsArg(&m_inputFormats)
48 196 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER});
49 98 : AddInputDatasetArg(&m_dataset, GDAL_OF_RASTER);
50 98 : AddOutputFormatArg(&m_outputFormat)
51 98 : .SetDefault(m_outputFormat)
52 : .AddMetadataItem(
53 : GAAMDI_REQUIRED_CAPABILITIES,
54 490 : {GDAL_DCAP_RASTER, GDAL_DCAP_CREATECOPY, GDAL_DMD_EXTENSIONS})
55 196 : .AddMetadataItem(GAAMDI_VRT_COMPATIBLE, {"false"});
56 98 : AddCreationOptionsArg(&m_creationOptions);
57 :
58 196 : AddArg("output", 'o', _("Output directory"), &m_outputDirectory)
59 98 : .SetRequired()
60 98 : .SetMinCharCount(1)
61 98 : .SetPositional();
62 :
63 392 : std::vector<std::string> tilingSchemes{"raster"};
64 882 : for (const std::string &scheme :
65 1862 : gdal::TileMatrixSet::listPredefinedTileMatrixSets())
66 : {
67 1764 : auto poTMS = gdal::TileMatrixSet::parse(scheme.c_str());
68 1764 : OGRSpatialReference oSRS_TMS;
69 1764 : if (poTMS && !poTMS->hasVariableMatrixWidth() &&
70 882 : oSRS_TMS.SetFromUserInput(poTMS->crs().c_str()) == OGRERR_NONE)
71 : {
72 882 : std::string identifier = scheme == "GoogleMapsCompatible"
73 : ? "WebMercatorQuad"
74 1764 : : poTMS->identifier();
75 882 : m_mapTileMatrixIdentifierToScheme[identifier] = scheme;
76 882 : tilingSchemes.push_back(std::move(identifier));
77 : }
78 : }
79 196 : AddArg("tiling-scheme", 0, _("Tiling scheme"), &m_tilingScheme)
80 98 : .SetDefault("WebMercatorQuad")
81 98 : .SetChoices(tilingSchemes)
82 : .SetHiddenChoices(
83 : "GoogleMapsCompatible", // equivalent of WebMercatorQuad
84 : "mercator", // gdal2tiles equivalent of WebMercatorQuad
85 : "geodetic" // gdal2tiles (not totally) equivalent of WorldCRS84Quad
86 98 : );
87 :
88 196 : AddArg("min-zoom", 0, _("Minimum zoom level"), &m_minZoomLevel)
89 98 : .SetMinValueIncluded(0);
90 196 : AddArg("max-zoom", 0, _("Maximum zoom level"), &m_maxZoomLevel)
91 98 : .SetMinValueIncluded(0);
92 :
93 196 : AddArg("min-x", 0, _("Minimum tile X coordinate"), &m_minTileX)
94 98 : .SetMinValueIncluded(0);
95 196 : AddArg("max-x", 0, _("Maximum tile X coordinate"), &m_maxTileX)
96 98 : .SetMinValueIncluded(0);
97 196 : AddArg("min-y", 0, _("Minimum tile Y coordinate"), &m_minTileY)
98 98 : .SetMinValueIncluded(0);
99 196 : AddArg("max-y", 0, _("Maximum tile Y coordinate"), &m_maxTileY)
100 98 : .SetMinValueIncluded(0);
101 : AddArg("no-intersection-ok", 0,
102 : _("Whether dataset extent not intersecting tile matrix is only a "
103 : "warning"),
104 98 : &m_noIntersectionIsOK);
105 :
106 : AddArg("resampling", 'r', _("Resampling method for max zoom"),
107 196 : &m_resampling)
108 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
109 : "average", "rms", "mode", "min", "max", "med", "q1", "q3",
110 98 : "sum")
111 98 : .SetDefault("cubic")
112 98 : .SetHiddenChoices("near");
113 : AddArg("overview-resampling", 0, _("Resampling method for overviews"),
114 196 : &m_overviewResampling)
115 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
116 : "average", "rms", "mode", "min", "max", "med", "q1", "q3",
117 98 : "sum")
118 98 : .SetHiddenChoices("near");
119 :
120 : AddArg("convention", 0,
121 : _("Tile numbering convention: xyz (from top) or tms (from bottom)"),
122 196 : &m_convention)
123 98 : .SetDefault(m_convention)
124 98 : .SetChoices("xyz", "tms");
125 196 : AddArg("tile-size", 0, _("Override default tile size"), &m_tileSize)
126 98 : .SetMinValueIncluded(64)
127 98 : .SetMaxValueIncluded(32768);
128 : AddArg("add-alpha", 0, _("Whether to force adding an alpha channel"),
129 196 : &m_addalpha)
130 98 : .SetMutualExclusionGroup("alpha");
131 : AddArg("no-alpha", 0, _("Whether to disable adding an alpha channel"),
132 196 : &m_noalpha)
133 98 : .SetMutualExclusionGroup("alpha");
134 : auto &dstNoDataArg =
135 98 : AddArg("dst-nodata", 0, _("Destination nodata value"), &m_dstNoData);
136 98 : AddArg("skip-blank", 0, _("Do not generate blank tiles"), &m_skipBlank);
137 :
138 : {
139 : auto &arg = AddArg("metadata", 0,
140 196 : _("Add metadata item to output tiles"), &m_metadata)
141 196 : .SetMetaVar("<KEY>=<VALUE>")
142 98 : .SetPackedValuesAllowed(false);
143 4 : arg.AddValidationAction([this, &arg]()
144 102 : { return ParseAndValidateKeyValue(arg); });
145 98 : arg.AddHiddenAlias("mo");
146 : }
147 : AddArg("copy-src-metadata", 0,
148 : _("Whether to copy metadata from source dataset"),
149 98 : &m_copySrcMetadata);
150 :
151 : AddArg("aux-xml", 0, _("Generate .aux.xml sidecar files when needed"),
152 98 : &m_auxXML);
153 98 : AddArg("kml", 0, _("Generate KML files"), &m_kml);
154 98 : AddArg("resume", 0, _("Generate only missing files"), &m_resume);
155 :
156 98 : AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
157 :
158 98 : constexpr const char *ADVANCED_RESAMPLING_CATEGORY = "Advanced Resampling";
159 : auto &excludedValuesArg =
160 : AddArg("excluded-values", 0,
161 : _("Tuples of values (e.g. <R>,<G>,<B> or (<R1>,<G1>,<B1>),"
162 : "(<R2>,<G2>,<B2>)) that must beignored as contributing source "
163 : "pixels during (average) resampling"),
164 196 : &m_excludedValues)
165 98 : .SetCategory(ADVANCED_RESAMPLING_CATEGORY);
166 : auto &excludedValuesPctThresholdArg =
167 : AddArg(
168 : "excluded-values-pct-threshold", 0,
169 : _("Minimum percentage of source pixels that must be set at one of "
170 : "the --excluded-values to cause the excluded value to be used as "
171 : "the target pixel value"),
172 196 : &m_excludedValuesPctThreshold)
173 98 : .SetDefault(m_excludedValuesPctThreshold)
174 98 : .SetMinValueIncluded(0)
175 98 : .SetMaxValueIncluded(100)
176 98 : .SetCategory(ADVANCED_RESAMPLING_CATEGORY);
177 : auto &nodataValuesPctThresholdArg =
178 : AddArg(
179 : "nodata-values-pct-threshold", 0,
180 : _("Minimum percentage of source pixels that must be set at one of "
181 : "nodata (or alpha=0 or any other way to express transparent pixel"
182 : "to cause the target pixel value to be transparent"),
183 196 : &m_nodataValuesPctThreshold)
184 98 : .SetDefault(m_nodataValuesPctThreshold)
185 98 : .SetMinValueIncluded(0)
186 98 : .SetMaxValueIncluded(100)
187 98 : .SetCategory(ADVANCED_RESAMPLING_CATEGORY);
188 :
189 98 : constexpr const char *PUBLICATION_CATEGORY = "Publication";
190 196 : AddArg("webviewer", 0, _("Web viewer to generate"), &m_webviewers)
191 98 : .SetDefault("all")
192 98 : .SetChoices("none", "all", "leaflet", "openlayers", "mapml")
193 98 : .SetCategory(PUBLICATION_CATEGORY);
194 : AddArg("url", 0,
195 : _("URL address where the generated tiles are going to be published"),
196 196 : &m_url)
197 98 : .SetCategory(PUBLICATION_CATEGORY);
198 196 : AddArg("title", 0, _("Title of the map"), &m_title)
199 98 : .SetCategory(PUBLICATION_CATEGORY);
200 196 : AddArg("copyright", 0, _("Copyright for the map"), &m_copyright)
201 98 : .SetCategory(PUBLICATION_CATEGORY);
202 : AddArg("mapml-template", 0,
203 : _("Filename of a template mapml file where variables will be "
204 : "substituted"),
205 196 : &m_mapmlTemplate)
206 98 : .SetMinCharCount(1)
207 98 : .SetCategory(PUBLICATION_CATEGORY);
208 :
209 98 : AddValidationAction(
210 93 : [this, &dstNoDataArg, &excludedValuesArg,
211 564 : &excludedValuesPctThresholdArg, &nodataValuesPctThresholdArg]()
212 : {
213 93 : if (m_minTileX >= 0 && m_maxTileX >= 0 && m_minTileX > m_maxTileX)
214 : {
215 1 : ReportError(CE_Failure, CPLE_IllegalArg,
216 : "'min-x' must be lesser or equal to 'max-x'");
217 1 : return false;
218 : }
219 :
220 92 : if (m_minTileY >= 0 && m_maxTileY >= 0 && m_minTileY > m_maxTileY)
221 : {
222 1 : ReportError(CE_Failure, CPLE_IllegalArg,
223 : "'min-y' must be lesser or equal to 'max-y'");
224 1 : return false;
225 : }
226 :
227 91 : if (m_minZoomLevel >= 0 && m_maxZoomLevel >= 0 &&
228 15 : m_minZoomLevel > m_maxZoomLevel)
229 : {
230 1 : ReportError(CE_Failure, CPLE_IllegalArg,
231 : "'min-zoom' must be lesser or equal to 'max-zoom'");
232 1 : return false;
233 : }
234 :
235 90 : if (m_addalpha && dstNoDataArg.IsExplicitlySet())
236 : {
237 1 : ReportError(
238 : CE_Failure, CPLE_IllegalArg,
239 : "'add-alpha' and 'dst-nodata' are mutually exclusive");
240 1 : return false;
241 : }
242 :
243 261 : for (const auto *arg :
244 : {&excludedValuesArg, &excludedValuesPctThresholdArg,
245 350 : &nodataValuesPctThresholdArg})
246 : {
247 263 : if (arg->IsExplicitlySet() && m_resampling != "average")
248 : {
249 1 : ReportError(CE_Failure, CPLE_AppDefined,
250 : "'%s' can only be specified if 'resampling' is "
251 : "set to 'average'",
252 1 : arg->GetName().c_str());
253 2 : return false;
254 : }
255 263 : if (arg->IsExplicitlySet() && !m_overviewResampling.empty() &&
256 1 : m_overviewResampling != "average")
257 : {
258 1 : ReportError(CE_Failure, CPLE_AppDefined,
259 : "'%s' can only be specified if "
260 : "'overview-resampling' is set to 'average'",
261 1 : arg->GetName().c_str());
262 1 : return false;
263 : }
264 : }
265 :
266 87 : return true;
267 : });
268 98 : }
269 :
270 : /************************************************************************/
271 : /* GetTileIndices() */
272 : /************************************************************************/
273 :
274 177 : static bool GetTileIndices(gdal::TileMatrixSet::TileMatrix &tileMatrix,
275 : bool bInvertAxisTMS, int tileSize,
276 : const double adfExtent[4], int &nMinTileX,
277 : int &nMinTileY, int &nMaxTileX, int &nMaxTileY,
278 : bool noIntersectionIsOK, bool &bIntersects,
279 : bool checkRasterOverflow = true)
280 : {
281 177 : if (tileSize > 0)
282 : {
283 18 : tileMatrix.mResX *=
284 18 : static_cast<double>(tileMatrix.mTileWidth) / tileSize;
285 18 : tileMatrix.mResY *=
286 18 : static_cast<double>(tileMatrix.mTileHeight) / tileSize;
287 18 : tileMatrix.mTileWidth = tileSize;
288 18 : tileMatrix.mTileHeight = tileSize;
289 : }
290 :
291 177 : if (bInvertAxisTMS)
292 2 : std::swap(tileMatrix.mTopLeftX, tileMatrix.mTopLeftY);
293 :
294 177 : const double dfTileWidth = tileMatrix.mResX * tileMatrix.mTileWidth;
295 177 : const double dfTileHeight = tileMatrix.mResY * tileMatrix.mTileHeight;
296 :
297 177 : constexpr double EPSILON = 1e-3;
298 177 : const double dfMinTileX =
299 177 : (adfExtent[0] - tileMatrix.mTopLeftX) / dfTileWidth;
300 177 : nMinTileX = static_cast<int>(
301 354 : std::clamp(std::floor(dfMinTileX + EPSILON), 0.0,
302 177 : static_cast<double>(tileMatrix.mMatrixWidth - 1)));
303 177 : const double dfMinTileY =
304 177 : (tileMatrix.mTopLeftY - adfExtent[3]) / dfTileHeight;
305 177 : nMinTileY = static_cast<int>(
306 354 : std::clamp(std::floor(dfMinTileY + EPSILON), 0.0,
307 177 : static_cast<double>(tileMatrix.mMatrixHeight - 1)));
308 177 : const double dfMaxTileX =
309 177 : (adfExtent[2] - tileMatrix.mTopLeftX) / dfTileWidth;
310 177 : nMaxTileX = static_cast<int>(
311 354 : std::clamp(std::floor(dfMaxTileX + EPSILON), 0.0,
312 177 : static_cast<double>(tileMatrix.mMatrixWidth - 1)));
313 177 : const double dfMaxTileY =
314 177 : (tileMatrix.mTopLeftY - adfExtent[1]) / dfTileHeight;
315 177 : nMaxTileY = static_cast<int>(
316 354 : std::clamp(std::floor(dfMaxTileY + EPSILON), 0.0,
317 177 : static_cast<double>(tileMatrix.mMatrixHeight - 1)));
318 :
319 177 : bIntersects = (dfMinTileX <= tileMatrix.mMatrixWidth && dfMaxTileX >= 0 &&
320 354 : dfMinTileY <= tileMatrix.mMatrixHeight && dfMaxTileY >= 0);
321 177 : if (!bIntersects)
322 : {
323 2 : CPLDebug("gdal_raster_tile",
324 : "dfMinTileX=%g dfMinTileY=%g dfMaxTileX=%g dfMaxTileY=%g",
325 : dfMinTileX, dfMinTileY, dfMaxTileX, dfMaxTileY);
326 2 : CPLError(noIntersectionIsOK ? CE_Warning : CE_Failure, CPLE_AppDefined,
327 : "Extent of source dataset is not compatible with extent of "
328 : "tile matrix %s",
329 : tileMatrix.mId.c_str());
330 2 : return noIntersectionIsOK;
331 : }
332 175 : if (checkRasterOverflow)
333 : {
334 107 : if (nMaxTileX - nMinTileX + 1 > INT_MAX / tileMatrix.mTileWidth ||
335 107 : nMaxTileY - nMinTileY + 1 > INT_MAX / tileMatrix.mTileHeight)
336 : {
337 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too large zoom level");
338 0 : return false;
339 : }
340 : }
341 175 : return true;
342 : }
343 :
344 : /************************************************************************/
345 : /* GetFileY() */
346 : /************************************************************************/
347 :
348 1981 : static int GetFileY(int iY, const gdal::TileMatrixSet::TileMatrix &tileMatrix,
349 : const std::string &convention)
350 : {
351 1981 : return convention == "xyz" ? iY : tileMatrix.mMatrixHeight - 1 - iY;
352 : }
353 :
354 : /************************************************************************/
355 : /* GenerateTile() */
356 : /************************************************************************/
357 :
358 263 : static bool GenerateTile(
359 : GDALDataset *poSrcDS, GDALDriver *poDstDriver, const char *pszExtension,
360 : CSLConstList creationOptions, GDALWarpOperation &oWO,
361 : const OGRSpatialReference &oSRS_TMS, GDALDataType eWorkingDataType,
362 : const gdal::TileMatrixSet::TileMatrix &tileMatrix,
363 : const std::string &outputDirectory, int nBands, const double *pdfDstNoData,
364 : int nZoomLevel, int iX, int iY, const std::string &convention,
365 : int nMinTileX, int nMinTileY, bool bSkipBlank, bool bUserAskedForAlpha,
366 : bool bAuxXML, bool bResume, const std::vector<std::string> &metadata,
367 : const GDALColorTable *poColorTable, std::vector<GByte> &dstBuffer)
368 : {
369 : const std::string osDirZ = CPLFormFilenameSafe(
370 526 : outputDirectory.c_str(), CPLSPrintf("%d", nZoomLevel), nullptr);
371 : const std::string osDirX =
372 526 : CPLFormFilenameSafe(osDirZ.c_str(), CPLSPrintf("%d", iX), nullptr);
373 263 : const int iFileY = GetFileY(iY, tileMatrix, convention);
374 : const std::string osFilename = CPLFormFilenameSafe(
375 526 : osDirX.c_str(), CPLSPrintf("%d", iFileY), pszExtension);
376 :
377 263 : if (bResume)
378 : {
379 : VSIStatBufL sStat;
380 5 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
381 5 : return true;
382 : }
383 :
384 258 : const int nDstXOff = (iX - nMinTileX) * tileMatrix.mTileWidth;
385 258 : const int nDstYOff = (iY - nMinTileY) * tileMatrix.mTileHeight;
386 258 : memset(dstBuffer.data(), 0, dstBuffer.size());
387 516 : const CPLErr eErr = oWO.WarpRegionToBuffer(
388 258 : nDstXOff, nDstYOff, tileMatrix.mTileWidth, tileMatrix.mTileHeight,
389 258 : dstBuffer.data(), eWorkingDataType);
390 258 : if (eErr != CE_None)
391 0 : return false;
392 :
393 : const bool bDstHasAlpha =
394 287 : nBands > poSrcDS->GetRasterCount() ||
395 29 : (nBands == poSrcDS->GetRasterCount() &&
396 28 : poSrcDS->GetRasterBand(nBands)->GetColorInterpretation() ==
397 258 : GCI_AlphaBand);
398 258 : const size_t nBytesPerBand = static_cast<size_t>(tileMatrix.mTileWidth) *
399 258 : tileMatrix.mTileHeight *
400 258 : GDALGetDataTypeSizeBytes(eWorkingDataType);
401 258 : if (bDstHasAlpha && bSkipBlank)
402 : {
403 9 : bool bBlank = true;
404 321952 : for (size_t i = 0; i < nBytesPerBand && bBlank; ++i)
405 : {
406 322022 : bBlank = (dstBuffer[(nBands - 1) * nBytesPerBand + i] == 0);
407 : }
408 0 : if (bBlank)
409 5 : return true;
410 : }
411 174 : if (bDstHasAlpha && !bUserAskedForAlpha)
412 : {
413 238 : bool bAllOpaque = true;
414 11162700 : for (size_t i = 0; i < nBytesPerBand && bAllOpaque; ++i)
415 : {
416 11174300 : bAllOpaque = (dstBuffer[(nBands - 1) * nBytesPerBand + i] == 255);
417 : }
418 0 : if (bAllOpaque)
419 163 : nBands--;
420 : }
421 :
422 0 : VSIMkdir(osDirZ.c_str(), 0755);
423 253 : VSIMkdir(osDirX.c_str(), 0755);
424 :
425 : auto memDS = std::unique_ptr<GDALDataset>(
426 253 : MEMDataset::Create("", tileMatrix.mTileWidth, tileMatrix.mTileHeight, 0,
427 506 : eWorkingDataType, nullptr));
428 1023 : for (int i = 0; i < nBands; ++i)
429 : {
430 770 : char szBuffer[32] = {'\0'};
431 1540 : int nRet = CPLPrintPointer(
432 770 : szBuffer, dstBuffer.data() + i * nBytesPerBand, sizeof(szBuffer));
433 770 : szBuffer[nRet] = 0;
434 :
435 770 : char szOption[64] = {'\0'};
436 770 : snprintf(szOption, sizeof(szOption), "DATAPOINTER=%s", szBuffer);
437 :
438 770 : char *apszOptions[] = {szOption, nullptr};
439 :
440 770 : memDS->AddBand(eWorkingDataType, apszOptions);
441 770 : auto poDstBand = memDS->GetRasterBand(i + 1);
442 770 : if (i + 1 <= poSrcDS->GetRasterCount())
443 695 : poDstBand->SetColorInterpretation(
444 695 : poSrcDS->GetRasterBand(i + 1)->GetColorInterpretation());
445 : else
446 75 : poDstBand->SetColorInterpretation(GCI_AlphaBand);
447 770 : if (pdfDstNoData)
448 9 : poDstBand->SetNoDataValue(*pdfDstNoData);
449 770 : if (i == 0 && poColorTable)
450 1 : poDstBand->SetColorTable(
451 1 : const_cast<GDALColorTable *>(poColorTable));
452 : }
453 506 : const CPLStringList aosMD(metadata);
454 261 : for (const auto [key, value] : cpl::IterateNameValue(aosMD))
455 : {
456 8 : memDS->SetMetadataItem(key, value);
457 : }
458 :
459 : double adfGT[6];
460 252 : adfGT[0] =
461 252 : tileMatrix.mTopLeftX + iX * tileMatrix.mResX * tileMatrix.mTileWidth;
462 252 : adfGT[1] = tileMatrix.mResX;
463 252 : adfGT[2] = 0;
464 252 : adfGT[3] =
465 252 : tileMatrix.mTopLeftY - iY * tileMatrix.mResY * tileMatrix.mTileHeight;
466 252 : adfGT[4] = 0;
467 252 : adfGT[5] = -tileMatrix.mResY;
468 252 : memDS->SetGeoTransform(adfGT);
469 :
470 253 : memDS->SetSpatialRef(&oSRS_TMS);
471 :
472 : CPLConfigOptionSetter oSetter("GDAL_PAM_ENABLED", bAuxXML ? "YES" : "NO",
473 506 : false);
474 :
475 505 : const std::string osTmpFilename = osFilename + ".tmp." + pszExtension;
476 :
477 : std::unique_ptr<GDALDataset> poOutDS(
478 : poDstDriver->CreateCopy(osTmpFilename.c_str(), memDS.get(), false,
479 253 : creationOptions, nullptr, nullptr));
480 253 : bool bRet = poOutDS && poOutDS->Close() == CE_None;
481 253 : poOutDS.reset();
482 253 : if (bRet)
483 : {
484 236 : bRet = VSIRename(osTmpFilename.c_str(), osFilename.c_str()) == 0;
485 236 : if (bAuxXML)
486 : {
487 4 : VSIRename((osTmpFilename + ".aux.xml").c_str(),
488 8 : (osFilename + ".aux.xml").c_str());
489 : }
490 : }
491 : else
492 : {
493 17 : VSIUnlink(osTmpFilename.c_str());
494 : }
495 253 : return bRet;
496 : }
497 :
498 : /************************************************************************/
499 : /* GenerateOverviewTile() */
500 : /************************************************************************/
501 :
502 : static bool
503 71 : GenerateOverviewTile(GDALDataset &oSrcDS, GDALDriver *poDstDriver,
504 : const std::string &outputFormat, const char *pszExtension,
505 : CSLConstList creationOptions,
506 : CSLConstList papszWarpOptions,
507 : const std::string &resampling,
508 : const gdal::TileMatrixSet::TileMatrix &tileMatrix,
509 : const std::string &outputDirectory, int nZoomLevel, int iX,
510 : int iY, const std::string &convention, bool bSkipBlank,
511 : bool bUserAskedForAlpha, bool bAuxXML, bool bResume)
512 : {
513 : const std::string osDirZ = CPLFormFilenameSafe(
514 142 : outputDirectory.c_str(), CPLSPrintf("%d", nZoomLevel), nullptr);
515 : const std::string osDirX =
516 142 : CPLFormFilenameSafe(osDirZ.c_str(), CPLSPrintf("%d", iX), nullptr);
517 :
518 71 : const int iFileY = GetFileY(iY, tileMatrix, convention);
519 : const std::string osFilename = CPLFormFilenameSafe(
520 142 : osDirX.c_str(), CPLSPrintf("%d", iFileY), pszExtension);
521 :
522 71 : if (bResume)
523 : {
524 : VSIStatBufL sStat;
525 2 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
526 1 : return true;
527 : }
528 :
529 70 : VSIMkdir(osDirZ.c_str(), 0755);
530 70 : VSIMkdir(osDirX.c_str(), 0755);
531 :
532 140 : CPLStringList aosOptions;
533 :
534 70 : aosOptions.AddString("-of");
535 70 : aosOptions.AddString(outputFormat.c_str());
536 :
537 94 : for (const char *pszCO : cpl::Iterate(creationOptions))
538 : {
539 24 : aosOptions.AddString("-co");
540 24 : aosOptions.AddString(pszCO);
541 : }
542 : CPLConfigOptionSetter oSetter("GDAL_PAM_ENABLED", bAuxXML ? "YES" : "NO",
543 140 : false);
544 :
545 70 : aosOptions.AddString("-r");
546 70 : aosOptions.AddString(resampling.c_str());
547 :
548 70 : std::unique_ptr<GDALDataset> poOutDS;
549 70 : const double dfMinX =
550 70 : tileMatrix.mTopLeftX + iX * tileMatrix.mResX * tileMatrix.mTileWidth;
551 70 : const double dfMaxY =
552 70 : tileMatrix.mTopLeftY - iY * tileMatrix.mResY * tileMatrix.mTileHeight;
553 70 : const double dfMaxX = dfMinX + tileMatrix.mResX * tileMatrix.mTileWidth;
554 70 : const double dfMinY = dfMaxY - tileMatrix.mResY * tileMatrix.mTileHeight;
555 :
556 : const bool resamplingCompatibleOfTranslate =
557 199 : papszWarpOptions == nullptr &&
558 190 : (resampling == "nearest" || resampling == "average" ||
559 125 : resampling == "bilinear" || resampling == "cubic" ||
560 9 : resampling == "cubicspline" || resampling == "lanczos" ||
561 3 : resampling == "mode");
562 :
563 140 : const std::string osTmpFilename = osFilename + ".tmp." + pszExtension;
564 :
565 70 : if (resamplingCompatibleOfTranslate)
566 : {
567 : double adfUpperGT[6];
568 64 : oSrcDS.GetGeoTransform(adfUpperGT);
569 64 : const double dfMinXUpper = adfUpperGT[0];
570 : const double dfMaxXUpper =
571 64 : dfMinXUpper + adfUpperGT[1] * oSrcDS.GetRasterXSize();
572 64 : const double dfMaxYUpper = adfUpperGT[3];
573 : const double dfMinYUpper =
574 64 : dfMaxYUpper + adfUpperGT[5] * oSrcDS.GetRasterYSize();
575 64 : if (dfMinX >= dfMinXUpper && dfMaxX <= dfMaxXUpper &&
576 48 : dfMinY >= dfMinYUpper && dfMaxY <= dfMaxYUpper)
577 : {
578 : // If the overview tile is fully within the extent of the
579 : // upper zoom level, we can use GDALDataset::RasterIO() directly.
580 :
581 48 : const auto eDT = oSrcDS.GetRasterBand(1)->GetRasterDataType();
582 : const size_t nBytesPerBand =
583 48 : static_cast<size_t>(tileMatrix.mTileWidth) *
584 48 : tileMatrix.mTileHeight * GDALGetDataTypeSizeBytes(eDT);
585 : std::vector<GByte> dstBuffer(nBytesPerBand *
586 48 : oSrcDS.GetRasterCount());
587 :
588 48 : const double dfXOff = (dfMinX - dfMinXUpper) / adfUpperGT[1];
589 48 : const double dfYOff = (dfMaxYUpper - dfMaxY) / -adfUpperGT[5];
590 48 : const double dfXSize = (dfMaxX - dfMinX) / adfUpperGT[1];
591 48 : const double dfYSize = (dfMaxY - dfMinY) / -adfUpperGT[5];
592 : GDALRasterIOExtraArg sExtraArg;
593 48 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
594 48 : CPL_IGNORE_RET_VAL(sExtraArg.eResampleAlg);
595 48 : sExtraArg.eResampleAlg =
596 48 : GDALRasterIOGetResampleAlg(resampling.c_str());
597 48 : sExtraArg.dfXOff = dfXOff;
598 48 : sExtraArg.dfYOff = dfYOff;
599 48 : sExtraArg.dfXSize = dfXSize;
600 48 : sExtraArg.dfYSize = dfYSize;
601 48 : sExtraArg.bFloatingPointWindowValidity =
602 48 : sExtraArg.eResampleAlg != GRIORA_NearestNeighbour;
603 48 : constexpr double EPSILON = 1e-3;
604 48 : if (oSrcDS.RasterIO(GF_Read, static_cast<int>(dfXOff + EPSILON),
605 48 : static_cast<int>(dfYOff + EPSILON),
606 48 : static_cast<int>(dfXSize + 0.5),
607 48 : static_cast<int>(dfYSize + 0.5),
608 48 : dstBuffer.data(), tileMatrix.mTileWidth,
609 48 : tileMatrix.mTileHeight, eDT,
610 : oSrcDS.GetRasterCount(), nullptr, 0, 0, 0,
611 48 : &sExtraArg) == CE_None)
612 : {
613 47 : int nDstBands = oSrcDS.GetRasterCount();
614 : const bool bDstHasAlpha =
615 47 : oSrcDS.GetRasterBand(nDstBands)->GetColorInterpretation() ==
616 47 : GCI_AlphaBand;
617 47 : if (bDstHasAlpha && bSkipBlank)
618 : {
619 2 : bool bBlank = true;
620 65545 : for (size_t i = 0; i < nBytesPerBand && bBlank; ++i)
621 : {
622 65543 : bBlank =
623 65543 : (dstBuffer[(nDstBands - 1) * nBytesPerBand + i] ==
624 : 0);
625 : }
626 2 : if (bBlank)
627 1 : return true;
628 1 : bSkipBlank = false;
629 : }
630 46 : if (bDstHasAlpha && !bUserAskedForAlpha)
631 : {
632 46 : bool bAllOpaque = true;
633 2818100 : for (size_t i = 0; i < nBytesPerBand && bAllOpaque; ++i)
634 : {
635 2818050 : bAllOpaque =
636 2818050 : (dstBuffer[(nDstBands - 1) * nBytesPerBand + i] ==
637 : 255);
638 : }
639 46 : if (bAllOpaque)
640 43 : nDstBands--;
641 : }
642 :
643 92 : auto memDS = std::unique_ptr<GDALDataset>(MEMDataset::Create(
644 46 : "", tileMatrix.mTileWidth, tileMatrix.mTileHeight, 0, eDT,
645 92 : nullptr));
646 187 : for (int i = 0; i < nDstBands; ++i)
647 : {
648 141 : char szBuffer[32] = {'\0'};
649 282 : int nRet = CPLPrintPointer(
650 141 : szBuffer, dstBuffer.data() + i * nBytesPerBand,
651 : sizeof(szBuffer));
652 141 : szBuffer[nRet] = 0;
653 :
654 141 : char szOption[64] = {'\0'};
655 141 : snprintf(szOption, sizeof(szOption), "DATAPOINTER=%s",
656 : szBuffer);
657 :
658 141 : char *apszOptions[] = {szOption, nullptr};
659 :
660 141 : memDS->AddBand(eDT, apszOptions);
661 141 : auto poSrcBand = oSrcDS.GetRasterBand(i + 1);
662 141 : auto poDstBand = memDS->GetRasterBand(i + 1);
663 141 : poDstBand->SetColorInterpretation(
664 141 : poSrcBand->GetColorInterpretation());
665 141 : int bHasNoData = false;
666 : const double dfNoData =
667 141 : poSrcBand->GetNoDataValue(&bHasNoData);
668 141 : if (bHasNoData)
669 0 : poDstBand->SetNoDataValue(dfNoData);
670 141 : if (auto poCT = poSrcBand->GetColorTable())
671 0 : poDstBand->SetColorTable(poCT);
672 : }
673 46 : memDS->SetMetadata(oSrcDS.GetMetadata());
674 : double adfGT[6];
675 46 : adfGT[0] = dfMinX;
676 46 : adfGT[1] = tileMatrix.mResX;
677 46 : adfGT[2] = 0;
678 46 : adfGT[3] = dfMaxY;
679 46 : adfGT[4] = 0;
680 46 : adfGT[5] = -tileMatrix.mResY;
681 46 : memDS->SetGeoTransform(adfGT);
682 :
683 46 : memDS->SetSpatialRef(oSrcDS.GetSpatialRef());
684 :
685 46 : poOutDS.reset(poDstDriver->CreateCopy(
686 : osTmpFilename.c_str(), memDS.get(), false, creationOptions,
687 : nullptr, nullptr));
688 47 : }
689 : }
690 : else
691 : {
692 : // If the overview tile is not fully within the extent of the
693 : // upper zoom level, use GDALTranslate() to use VRT padding
694 :
695 16 : aosOptions.AddString("-q");
696 :
697 16 : aosOptions.AddString("-projwin");
698 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
699 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY));
700 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
701 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinY));
702 :
703 16 : aosOptions.AddString("-outsize");
704 16 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileWidth));
705 16 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileHeight));
706 :
707 : GDALTranslateOptions *psOptions =
708 16 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
709 16 : poOutDS.reset(GDALDataset::FromHandle(GDALTranslate(
710 : osTmpFilename.c_str(), GDALDataset::ToHandle(&oSrcDS),
711 : psOptions, nullptr)));
712 16 : GDALTranslateOptionsFree(psOptions);
713 : }
714 : }
715 : else
716 : {
717 6 : aosOptions.AddString("-te");
718 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
719 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinY));
720 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
721 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY));
722 :
723 6 : aosOptions.AddString("-ts");
724 6 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileWidth));
725 6 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileHeight));
726 :
727 11 : for (int i = 0; papszWarpOptions && papszWarpOptions[i]; ++i)
728 : {
729 5 : aosOptions.AddString("-wo");
730 5 : aosOptions.AddString(papszWarpOptions[i]);
731 : }
732 :
733 : GDALWarpAppOptions *psOptions =
734 6 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
735 6 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(&oSrcDS);
736 6 : poOutDS.reset(GDALDataset::FromHandle(GDALWarp(
737 : osTmpFilename.c_str(), nullptr, 1, &hSrcDS, psOptions, nullptr)));
738 6 : GDALWarpAppOptionsFree(psOptions);
739 : }
740 :
741 69 : bool bRet = poOutDS != nullptr;
742 69 : if (bRet && bSkipBlank)
743 : {
744 10 : auto poLastBand = poOutDS->GetRasterBand(poOutDS->GetRasterCount());
745 10 : if (poLastBand->GetColorInterpretation() == GCI_AlphaBand)
746 : {
747 : std::vector<GByte> buffer(
748 1 : static_cast<size_t>(tileMatrix.mTileWidth) *
749 1 : tileMatrix.mTileHeight *
750 1 : GDALGetDataTypeSizeBytes(poLastBand->GetRasterDataType()));
751 2 : CPL_IGNORE_RET_VAL(poLastBand->RasterIO(
752 1 : GF_Read, 0, 0, tileMatrix.mTileWidth, tileMatrix.mTileHeight,
753 1 : buffer.data(), tileMatrix.mTileWidth, tileMatrix.mTileHeight,
754 : poLastBand->GetRasterDataType(), 0, 0, nullptr));
755 1 : bool bBlank = true;
756 65537 : for (size_t i = 0; i < buffer.size() && bBlank; ++i)
757 : {
758 65536 : bBlank = (buffer[i] == 0);
759 : }
760 1 : if (bBlank)
761 : {
762 1 : poOutDS.reset();
763 1 : VSIUnlink(osTmpFilename.c_str());
764 1 : if (bAuxXML)
765 0 : VSIUnlink((osTmpFilename + ".aux.xml").c_str());
766 1 : return true;
767 : }
768 : }
769 : }
770 68 : bRet = bRet && poOutDS->Close() == CE_None;
771 68 : poOutDS.reset();
772 68 : if (bRet)
773 : {
774 67 : bRet = VSIRename(osTmpFilename.c_str(), osFilename.c_str()) == 0;
775 67 : if (bAuxXML)
776 : {
777 4 : VSIRename((osTmpFilename + ".aux.xml").c_str(),
778 8 : (osFilename + ".aux.xml").c_str());
779 : }
780 : }
781 : else
782 : {
783 1 : VSIUnlink(osTmpFilename.c_str());
784 : }
785 68 : return bRet;
786 : }
787 :
788 : namespace
789 : {
790 :
791 : /************************************************************************/
792 : /* FakeMaxZoomRasterBand */
793 : /************************************************************************/
794 :
795 : class FakeMaxZoomRasterBand : public GDALRasterBand
796 : {
797 : void *m_pDstBuffer = nullptr;
798 : CPL_DISALLOW_COPY_ASSIGN(FakeMaxZoomRasterBand)
799 :
800 : public:
801 601 : FakeMaxZoomRasterBand(int nBandIn, int nWidth, int nHeight,
802 : int nBlockXSizeIn, int nBlockYSizeIn,
803 : GDALDataType eDT, void *pDstBuffer)
804 601 : : m_pDstBuffer(pDstBuffer)
805 : {
806 601 : nBand = nBandIn;
807 601 : nRasterXSize = nWidth;
808 601 : nRasterYSize = nHeight;
809 601 : nBlockXSize = nBlockXSizeIn;
810 601 : nBlockYSize = nBlockYSizeIn;
811 601 : eDataType = eDT;
812 601 : }
813 :
814 0 : CPLErr IReadBlock(int, int, void *) override
815 : {
816 0 : CPLAssert(false);
817 : return CE_Failure;
818 : }
819 :
820 : #ifdef DEBUG
821 0 : CPLErr IWriteBlock(int, int, void *) override
822 : {
823 0 : CPLAssert(false);
824 : return CE_Failure;
825 : }
826 : #endif
827 :
828 490 : CPLErr IRasterIO(GDALRWFlag eRWFlag, [[maybe_unused]] int nXOff,
829 : [[maybe_unused]] int nYOff, [[maybe_unused]] int nXSize,
830 : [[maybe_unused]] int nYSize, void *pData,
831 : [[maybe_unused]] int nBufXSize,
832 : [[maybe_unused]] int nBufYSize, GDALDataType eBufType,
833 : GSpacing nPixelSpace, [[maybe_unused]] GSpacing nLineSpace,
834 : GDALRasterIOExtraArg *) override
835 : {
836 : // For sake of implementation simplicity, check various assumptions of
837 : // how GDALAlphaMask code does I/O
838 490 : CPLAssert((nXOff % nBlockXSize) == 0);
839 490 : CPLAssert((nYOff % nBlockYSize) == 0);
840 490 : CPLAssert(nXSize == nBufXSize);
841 490 : CPLAssert(nXSize == nBlockXSize);
842 490 : CPLAssert(nYSize == nBufYSize);
843 490 : CPLAssert(nYSize == nBlockYSize);
844 490 : CPLAssert(nLineSpace == nBlockXSize * nPixelSpace);
845 490 : CPLAssert(
846 : nBand ==
847 : poDS->GetRasterCount()); // only alpha band is accessed this way
848 490 : if (eRWFlag == GF_Read)
849 : {
850 245 : double dfZero = 0;
851 245 : GDALCopyWords64(&dfZero, GDT_Float64, 0, pData, eBufType,
852 : static_cast<int>(nPixelSpace),
853 245 : static_cast<size_t>(nBlockXSize) * nBlockYSize);
854 : }
855 : else
856 : {
857 245 : GDALCopyWords64(pData, eBufType, static_cast<int>(nPixelSpace),
858 : m_pDstBuffer, eDataType,
859 : GDALGetDataTypeSizeBytes(eDataType),
860 245 : static_cast<size_t>(nBlockXSize) * nBlockYSize);
861 : }
862 490 : return CE_None;
863 : }
864 : };
865 :
866 : /************************************************************************/
867 : /* FakeMaxZoomDataset */
868 : /************************************************************************/
869 :
870 : // This class is used to create a fake output dataset for GDALWarpOperation.
871 : // In particular we need to implement GDALRasterBand::IRasterIO(GF_Write, ...)
872 : // to catch writes (of one single tile) to the alpha band and redirect them
873 : // to the dstBuffer passed to FakeMaxZoomDataset constructor.
874 :
875 : class FakeMaxZoomDataset : public GDALDataset
876 : {
877 : const int m_nBlockXSize;
878 : const int m_nBlockYSize;
879 : const OGRSpatialReference m_oSRS;
880 : double m_adfGT[6];
881 :
882 : public:
883 172 : FakeMaxZoomDataset(int nWidth, int nHeight, int nBandsIn, int nBlockXSize,
884 : int nBlockYSize, GDALDataType eDT, const double adfGT[6],
885 : const OGRSpatialReference &oSRS,
886 : std::vector<GByte> &dstBuffer)
887 172 : : m_nBlockXSize(nBlockXSize), m_nBlockYSize(nBlockYSize), m_oSRS(oSRS)
888 : {
889 172 : eAccess = GA_Update;
890 172 : nRasterXSize = nWidth;
891 172 : nRasterYSize = nHeight;
892 172 : memcpy(m_adfGT, adfGT, sizeof(double) * 6);
893 773 : for (int i = 1; i <= nBandsIn; ++i)
894 : {
895 601 : SetBand(i,
896 : new FakeMaxZoomRasterBand(
897 : i, nWidth, nHeight, nBlockXSize, nBlockYSize, eDT,
898 601 : dstBuffer.data() + static_cast<size_t>(i - 1) *
899 1202 : nBlockXSize * nBlockYSize *
900 601 : GDALGetDataTypeSizeBytes(eDT)));
901 : }
902 172 : }
903 :
904 296 : const OGRSpatialReference *GetSpatialRef() const override
905 : {
906 296 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
907 : }
908 :
909 60 : CPLErr GetGeoTransform(double *padfGT) override
910 : {
911 60 : memcpy(padfGT, m_adfGT, sizeof(double) * 6);
912 60 : return CE_None;
913 : }
914 :
915 : using GDALDataset::Clone;
916 :
917 : std::unique_ptr<FakeMaxZoomDataset>
918 112 : Clone(std::vector<GByte> &dstBuffer) const
919 : {
920 : return std::make_unique<FakeMaxZoomDataset>(
921 112 : nRasterXSize, nRasterYSize, nBands, m_nBlockXSize, m_nBlockYSize,
922 224 : GetRasterBand(1)->GetRasterDataType(), m_adfGT, m_oSRS, dstBuffer);
923 : }
924 : };
925 :
926 : /************************************************************************/
927 : /* MosaicRasterBand */
928 : /************************************************************************/
929 :
930 : class MosaicRasterBand : public GDALRasterBand
931 : {
932 : const int m_tileMinX;
933 : const int m_tileMinY;
934 : const GDALColorInterp m_eColorInterp;
935 : const gdal::TileMatrixSet::TileMatrix m_oTM;
936 : const std::string m_convention;
937 : const std::string m_directory;
938 : const std::string m_extension;
939 : const bool m_hasNoData;
940 : const double m_noData;
941 : std::unique_ptr<GDALColorTable> m_poColorTable{};
942 :
943 : public:
944 158 : MosaicRasterBand(GDALDataset *poDSIn, int nBandIn, int nWidth, int nHeight,
945 : int nBlockXSizeIn, int nBlockYSizeIn, GDALDataType eDT,
946 : GDALColorInterp eColorInterp, int nTileMinX, int nTileMinY,
947 : const gdal::TileMatrixSet::TileMatrix &oTM,
948 : const std::string &convention,
949 : const std::string &directory, const std::string &extension,
950 : const double *pdfDstNoData,
951 : const GDALColorTable *poColorTable)
952 158 : : m_tileMinX(nTileMinX), m_tileMinY(nTileMinY),
953 : m_eColorInterp(eColorInterp), m_oTM(oTM), m_convention(convention),
954 : m_directory(directory), m_extension(extension),
955 158 : m_hasNoData(pdfDstNoData != nullptr),
956 158 : m_noData(pdfDstNoData ? *pdfDstNoData : 0),
957 316 : m_poColorTable(poColorTable ? poColorTable->Clone() : nullptr)
958 : {
959 158 : poDS = poDSIn;
960 158 : nBand = nBandIn;
961 158 : nRasterXSize = nWidth;
962 158 : nRasterYSize = nHeight;
963 158 : nBlockXSize = nBlockXSizeIn;
964 158 : nBlockYSize = nBlockYSizeIn;
965 158 : eDataType = eDT;
966 158 : }
967 :
968 : CPLErr IReadBlock(int nXBlock, int nYBlock, void *pData) override;
969 :
970 2024 : GDALColorTable *GetColorTable() override
971 : {
972 2024 : return m_poColorTable.get();
973 : }
974 :
975 325 : GDALColorInterp GetColorInterpretation() override
976 : {
977 325 : return m_eColorInterp;
978 : }
979 :
980 1635 : double GetNoDataValue(int *pbHasNoData) override
981 : {
982 1635 : if (pbHasNoData)
983 1626 : *pbHasNoData = m_hasNoData;
984 1635 : return m_noData;
985 : }
986 : };
987 :
988 : /************************************************************************/
989 : /* MosaicDataset */
990 : /************************************************************************/
991 :
992 : // This class is to expose the tiles of a given level as a mosaic that
993 : // can be used as a source to generate the immediately below zoom level.
994 :
995 : class MosaicDataset : public GDALDataset
996 : {
997 : friend class MosaicRasterBand;
998 :
999 : const std::string m_directory;
1000 : const std::string m_extension;
1001 : const std::string m_format;
1002 : GDALDataset *const m_poSrcDS;
1003 : const gdal::TileMatrixSet::TileMatrix &m_oTM;
1004 : const OGRSpatialReference m_oSRS;
1005 : const int m_nTileMinX;
1006 : const int m_nTileMinY;
1007 : const int m_nTileMaxX;
1008 : const int m_nTileMaxY;
1009 : const std::string m_convention;
1010 : const GDALDataType m_eDT;
1011 : const double *const m_pdfDstNoData;
1012 : const std::vector<std::string> &m_metadata;
1013 : const GDALColorTable *const m_poCT;
1014 :
1015 : double m_adfGT[6];
1016 : lru11::Cache<std::string, std::shared_ptr<GDALDataset>> m_oCacheTile{};
1017 :
1018 : CPL_DISALLOW_COPY_ASSIGN(MosaicDataset)
1019 :
1020 : public:
1021 51 : MosaicDataset(const std::string &directory, const std::string &extension,
1022 : const std::string &format, GDALDataset *poSrcDS,
1023 : const gdal::TileMatrixSet::TileMatrix &oTM,
1024 : const OGRSpatialReference &oSRS, int nTileMinX, int nTileMinY,
1025 : int nTileMaxX, int nTileMaxY, const std::string &convention,
1026 : int nBandsIn, GDALDataType eDT, const double *pdfDstNoData,
1027 : const std::vector<std::string> &metadata,
1028 : const GDALColorTable *poCT)
1029 51 : : m_directory(directory), m_extension(extension), m_format(format),
1030 : m_poSrcDS(poSrcDS), m_oTM(oTM), m_oSRS(oSRS), m_nTileMinX(nTileMinX),
1031 : m_nTileMinY(nTileMinY), m_nTileMaxX(nTileMaxX),
1032 : m_nTileMaxY(nTileMaxY), m_convention(convention), m_eDT(eDT),
1033 51 : m_pdfDstNoData(pdfDstNoData), m_metadata(metadata), m_poCT(poCT)
1034 : {
1035 51 : nRasterXSize = (nTileMaxX - nTileMinX + 1) * oTM.mTileWidth;
1036 51 : nRasterYSize = (nTileMaxY - nTileMinY + 1) * oTM.mTileHeight;
1037 51 : m_adfGT[0] = oTM.mTopLeftX + nTileMinX * oTM.mResX * oTM.mTileWidth;
1038 51 : m_adfGT[1] = oTM.mResX;
1039 51 : m_adfGT[2] = 0;
1040 51 : m_adfGT[3] = oTM.mTopLeftY - nTileMinY * oTM.mResY * oTM.mTileHeight;
1041 51 : m_adfGT[4] = 0;
1042 51 : m_adfGT[5] = -oTM.mResY;
1043 209 : for (int i = 1; i <= nBandsIn; ++i)
1044 : {
1045 : const GDALColorInterp eColorInterp =
1046 158 : (i <= poSrcDS->GetRasterCount())
1047 158 : ? poSrcDS->GetRasterBand(i)->GetColorInterpretation()
1048 158 : : GCI_AlphaBand;
1049 158 : SetBand(i, new MosaicRasterBand(
1050 158 : this, i, nRasterXSize, nRasterYSize, oTM.mTileWidth,
1051 158 : oTM.mTileHeight, eDT, eColorInterp, nTileMinX,
1052 : nTileMinY, oTM, convention, directory, extension,
1053 158 : pdfDstNoData, poCT));
1054 : }
1055 51 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1056 102 : const CPLStringList aosMD(metadata);
1057 59 : for (const auto [key, value] : cpl::IterateNameValue(aosMD))
1058 : {
1059 8 : SetMetadataItem(key, value);
1060 : }
1061 51 : }
1062 :
1063 92 : const OGRSpatialReference *GetSpatialRef() const override
1064 : {
1065 92 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
1066 : }
1067 :
1068 102 : CPLErr GetGeoTransform(double *padfGT) override
1069 : {
1070 102 : memcpy(padfGT, m_adfGT, sizeof(double) * 6);
1071 102 : return CE_None;
1072 : }
1073 :
1074 : using GDALDataset::Clone;
1075 :
1076 16 : std::unique_ptr<MosaicDataset> Clone() const
1077 : {
1078 : return std::make_unique<MosaicDataset>(
1079 16 : m_directory, m_extension, m_format, m_poSrcDS, m_oTM, m_oSRS,
1080 16 : m_nTileMinX, m_nTileMinY, m_nTileMaxX, m_nTileMaxY, m_convention,
1081 16 : nBands, m_eDT, m_pdfDstNoData, m_metadata, m_poCT);
1082 : }
1083 : };
1084 :
1085 : /************************************************************************/
1086 : /* MosaicRasterBand::IReadBlock() */
1087 : /************************************************************************/
1088 :
1089 1632 : CPLErr MosaicRasterBand::IReadBlock(int nXBlock, int nYBlock, void *pData)
1090 : {
1091 1632 : auto poThisDS = cpl::down_cast<MosaicDataset *>(poDS);
1092 : std::string filename = CPLFormFilenameSafe(
1093 3263 : m_directory.c_str(), CPLSPrintf("%d", m_tileMinX + nXBlock), nullptr);
1094 1631 : const int iFileY = GetFileY(m_tileMinY + nYBlock, m_oTM, m_convention);
1095 3264 : filename = CPLFormFilenameSafe(filename.c_str(), CPLSPrintf("%d", iFileY),
1096 1632 : m_extension.c_str());
1097 :
1098 1632 : std::shared_ptr<GDALDataset> poTileDS;
1099 1632 : if (!poThisDS->m_oCacheTile.tryGet(filename, poTileDS))
1100 : {
1101 422 : const char *const apszAllowedDrivers[] = {poThisDS->m_format.c_str(),
1102 422 : nullptr};
1103 422 : const char *const apszAllowedDriversForCOG[] = {"GTiff", "LIBERTIFF",
1104 : nullptr};
1105 422 : poTileDS.reset(GDALDataset::Open(
1106 : filename.c_str(), GDAL_OF_RASTER | GDAL_OF_INTERNAL,
1107 422 : EQUAL(poThisDS->m_format.c_str(), "COG") ? apszAllowedDriversForCOG
1108 : : apszAllowedDrivers));
1109 423 : if (!poTileDS)
1110 : {
1111 : VSIStatBufL sStat;
1112 6 : if (VSIStatL(filename.c_str(), &sStat) == 0)
1113 : {
1114 1 : CPLError(CE_Failure, CPLE_AppDefined,
1115 : "File %s exists but cannot be opened with %s driver",
1116 : filename.c_str(), poThisDS->m_format.c_str());
1117 1 : return CE_Failure;
1118 : }
1119 : }
1120 422 : poThisDS->m_oCacheTile.insert(filename, poTileDS);
1121 : }
1122 1631 : if (!poTileDS || nBand > poTileDS->GetRasterCount())
1123 : {
1124 767 : memset(pData,
1125 381 : (poTileDS && (nBand == poTileDS->GetRasterCount() + 1)) ? 255
1126 : : 0,
1127 386 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
1128 386 : GDALGetDataTypeSizeBytes(eDataType));
1129 386 : return CE_None;
1130 : }
1131 : else
1132 : {
1133 1244 : return poTileDS->GetRasterBand(nBand)->RasterIO(
1134 : GF_Read, 0, 0, nBlockXSize, nBlockYSize, pData, nBlockXSize,
1135 1245 : nBlockYSize, eDataType, 0, 0, nullptr);
1136 : }
1137 : }
1138 :
1139 : } // namespace
1140 :
1141 : /************************************************************************/
1142 : /* ApplySubstitutions() */
1143 : /************************************************************************/
1144 :
1145 59 : static void ApplySubstitutions(CPLString &s,
1146 : const std::map<std::string, std::string> &substs)
1147 : {
1148 932 : for (const auto &[key, value] : substs)
1149 : {
1150 873 : s.replaceAll("%(" + key + ")s", value);
1151 873 : s.replaceAll("%(" + key + ")d", value);
1152 873 : s.replaceAll("%(" + key + ")f", value);
1153 873 : s.replaceAll("${" + key + "}", value);
1154 : }
1155 59 : }
1156 :
1157 : /************************************************************************/
1158 : /* GenerateLeaflet() */
1159 : /************************************************************************/
1160 :
1161 13 : static void GenerateLeaflet(const std::string &osDirectory,
1162 : const std::string &osTitle, double dfSouthLat,
1163 : double dfWestLon, double dfNorthLat,
1164 : double dfEastLon, int nMinZoom, int nMaxZoom,
1165 : int nTileSize, const std::string &osExtension,
1166 : const std::string &osURL,
1167 : const std::string &osCopyright, bool bXYZ)
1168 : {
1169 13 : if (const char *pszTemplate = CPLFindFile("gdal", "leaflet_template.html"))
1170 : {
1171 26 : const std::string osFilename(pszTemplate);
1172 26 : std::map<std::string, std::string> substs;
1173 :
1174 : // For tests
1175 : const char *pszFmt =
1176 13 : atoi(CPLGetConfigOption("GDAL_RASTER_TILE_HTML_PREC", "17")) == 10
1177 : ? "%.10g"
1178 13 : : "%.17g";
1179 :
1180 26 : substs["double_quote_escaped_title"] =
1181 39 : CPLString(osTitle).replaceAll('"', "\\\"");
1182 13 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1183 13 : substs["xml_escaped_title"] = pszStr;
1184 13 : CPLFree(pszStr);
1185 13 : substs["south"] = CPLSPrintf(pszFmt, dfSouthLat);
1186 13 : substs["west"] = CPLSPrintf(pszFmt, dfWestLon);
1187 13 : substs["north"] = CPLSPrintf(pszFmt, dfNorthLat);
1188 13 : substs["east"] = CPLSPrintf(pszFmt, dfEastLon);
1189 13 : substs["centerlon"] = CPLSPrintf(pszFmt, (dfNorthLat + dfSouthLat) / 2);
1190 13 : substs["centerlat"] = CPLSPrintf(pszFmt, (dfWestLon + dfEastLon) / 2);
1191 13 : substs["minzoom"] = CPLSPrintf("%d", nMinZoom);
1192 13 : substs["maxzoom"] = CPLSPrintf("%d", nMaxZoom);
1193 13 : substs["beginzoom"] = CPLSPrintf("%d", nMaxZoom);
1194 13 : substs["tile_size"] = CPLSPrintf("%d", nTileSize); // not used
1195 13 : substs["tileformat"] = osExtension;
1196 13 : substs["publishurl"] = osURL; // not used
1197 13 : substs["copyright"] = CPLString(osCopyright).replaceAll('"', "\\\"");
1198 13 : substs["tms"] = bXYZ ? "0" : "1";
1199 :
1200 13 : GByte *pabyRet = nullptr;
1201 13 : CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, osFilename.c_str(), &pabyRet,
1202 : nullptr, 10 * 1024 * 1024));
1203 13 : if (pabyRet)
1204 : {
1205 26 : CPLString osHTML(reinterpret_cast<char *>(pabyRet));
1206 13 : CPLFree(pabyRet);
1207 :
1208 13 : ApplySubstitutions(osHTML, substs);
1209 :
1210 13 : VSILFILE *f = VSIFOpenL(CPLFormFilenameSafe(osDirectory.c_str(),
1211 : "leaflet.html", nullptr)
1212 : .c_str(),
1213 : "wb");
1214 13 : if (f)
1215 : {
1216 13 : VSIFWriteL(osHTML.data(), 1, osHTML.size(), f);
1217 13 : VSIFCloseL(f);
1218 : }
1219 : }
1220 : }
1221 13 : }
1222 :
1223 : /************************************************************************/
1224 : /* GenerateMapML() */
1225 : /************************************************************************/
1226 :
1227 : static void
1228 14 : GenerateMapML(const std::string &osDirectory, const std::string &mapmlTemplate,
1229 : const std::string &osTitle, int nMinTileX, int nMinTileY,
1230 : int nMaxTileX, int nMaxTileY, int nMinZoom, int nMaxZoom,
1231 : const std::string &osExtension, const std::string &osURL,
1232 : const std::string &osCopyright, const gdal::TileMatrixSet &tms)
1233 : {
1234 14 : if (const char *pszTemplate =
1235 14 : (mapmlTemplate.empty() ? CPLFindFile("gdal", "template_tiles.mapml")
1236 14 : : mapmlTemplate.c_str()))
1237 : {
1238 28 : const std::string osFilename(pszTemplate);
1239 28 : std::map<std::string, std::string> substs;
1240 :
1241 14 : if (tms.identifier() == "GoogleMapsCompatible")
1242 13 : substs["TILING_SCHEME"] = "OSMTILE";
1243 1 : else if (tms.identifier() == "WorldCRS84Quad")
1244 1 : substs["TILING_SCHEME"] = "WGS84";
1245 : else
1246 0 : substs["TILING_SCHEME"] = tms.identifier();
1247 :
1248 14 : substs["URL"] = osURL.empty() ? "./" : osURL;
1249 14 : substs["MINTILEX"] = CPLSPrintf("%d", nMinTileX);
1250 14 : substs["MINTILEY"] = CPLSPrintf("%d", nMinTileY);
1251 14 : substs["MAXTILEX"] = CPLSPrintf("%d", nMaxTileX);
1252 14 : substs["MAXTILEY"] = CPLSPrintf("%d", nMaxTileY);
1253 14 : substs["CURZOOM"] = CPLSPrintf("%d", nMaxZoom);
1254 14 : substs["MINZOOM"] = CPLSPrintf("%d", nMinZoom);
1255 14 : substs["MAXZOOM"] = CPLSPrintf("%d", nMaxZoom);
1256 14 : substs["TILEEXT"] = osExtension;
1257 14 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1258 14 : substs["TITLE"] = pszStr;
1259 14 : CPLFree(pszStr);
1260 14 : substs["COPYRIGHT"] = osCopyright;
1261 :
1262 14 : GByte *pabyRet = nullptr;
1263 14 : CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, osFilename.c_str(), &pabyRet,
1264 : nullptr, 10 * 1024 * 1024));
1265 14 : if (pabyRet)
1266 : {
1267 28 : CPLString osMAPML(reinterpret_cast<char *>(pabyRet));
1268 14 : CPLFree(pabyRet);
1269 :
1270 14 : ApplySubstitutions(osMAPML, substs);
1271 :
1272 14 : VSILFILE *f = VSIFOpenL(
1273 28 : CPLFormFilenameSafe(osDirectory.c_str(), "mapml.mapml", nullptr)
1274 : .c_str(),
1275 : "wb");
1276 14 : if (f)
1277 : {
1278 14 : VSIFWriteL(osMAPML.data(), 1, osMAPML.size(), f);
1279 14 : VSIFCloseL(f);
1280 : }
1281 : }
1282 : }
1283 14 : }
1284 :
1285 : /************************************************************************/
1286 : /* GenerateOpenLayers() */
1287 : /************************************************************************/
1288 :
1289 21 : static void GenerateOpenLayers(
1290 : const std::string &osDirectory, const std::string &osTitle, double dfMinX,
1291 : double dfMinY, double dfMaxX, double dfMaxY, int nMinZoom, int nMaxZoom,
1292 : int nTileSize, const std::string &osExtension, const std::string &osURL,
1293 : const std::string &osCopyright, const gdal::TileMatrixSet &tms,
1294 : bool bInvertAxisTMS, const OGRSpatialReference &oSRS_TMS, bool bXYZ)
1295 : {
1296 42 : std::map<std::string, std::string> substs;
1297 :
1298 : // For tests
1299 : const char *pszFmt =
1300 21 : atoi(CPLGetConfigOption("GDAL_RASTER_TILE_HTML_PREC", "17")) == 10
1301 : ? "%.10g"
1302 21 : : "%.17g";
1303 :
1304 21 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1305 21 : substs["xml_escaped_title"] = pszStr;
1306 21 : CPLFree(pszStr);
1307 21 : substs["ominx"] = CPLSPrintf(pszFmt, dfMinX);
1308 21 : substs["ominy"] = CPLSPrintf(pszFmt, dfMinY);
1309 21 : substs["omaxx"] = CPLSPrintf(pszFmt, dfMaxX);
1310 21 : substs["omaxy"] = CPLSPrintf(pszFmt, dfMaxY);
1311 21 : substs["center_x"] = CPLSPrintf(pszFmt, (dfMinX + dfMaxX) / 2);
1312 21 : substs["center_y"] = CPLSPrintf(pszFmt, (dfMinY + dfMaxY) / 2);
1313 21 : substs["minzoom"] = CPLSPrintf("%d", nMinZoom);
1314 21 : substs["maxzoom"] = CPLSPrintf("%d", nMaxZoom);
1315 21 : substs["tile_size"] = CPLSPrintf("%d", nTileSize);
1316 21 : substs["tileformat"] = osExtension;
1317 21 : substs["publishurl"] = osURL;
1318 21 : substs["copyright"] = osCopyright;
1319 21 : substs["sign_y"] = bXYZ ? "" : "-";
1320 :
1321 : CPLString s(R"raw(<!DOCTYPE html>
1322 : <html>
1323 : <head>
1324 : <title>%(xml_escaped_title)s</title>
1325 : <meta http-equiv="content-type" content="text/html; charset=utf-8"/>
1326 : <meta http-equiv='imagetoolbar' content='no'/>
1327 : <style type="text/css"> v\:* {behavior:url(#default#VML);}
1328 : html, body { overflow: hidden; padding: 0; height: 100%; width: 100%; font-family: 'Lucida Grande',Geneva,Arial,Verdana,sans-serif; }
1329 : body { margin: 10px; background: #fff; }
1330 : h1 { margin: 0; padding: 6px; border:0; font-size: 20pt; }
1331 : #header { height: 43px; padding: 0; background-color: #eee; border: 1px solid #888; }
1332 : #subheader { height: 12px; text-align: right; font-size: 10px; color: #555;}
1333 : #map { height: 90%; border: 1px solid #888; }
1334 : </style>
1335 : <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@main/dist/en/v7.0.0/legacy/ol.css" type="text/css">
1336 : <script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@main/dist/en/v7.0.0/legacy/ol.js"></script>
1337 : <script src="https://unpkg.com/ol-layerswitcher@4.1.1"></script>
1338 : <link rel="stylesheet" href="https://unpkg.com/ol-layerswitcher@4.1.1/src/ol-layerswitcher.css" />
1339 : </head>
1340 : <body>
1341 : <div id="header"><h1>%(xml_escaped_title)s</h1></div>
1342 : <div id="subheader">Generated by <a href="https://gdal.org/programs/gdal_raster_tile.html">gdal raster tile</a> </div>
1343 : <div id="map" class="map"></div>
1344 : <div id="mouse-position"></div>
1345 : <script type="text/javascript">
1346 : var mousePositionControl = new ol.control.MousePosition({
1347 : className: 'custom-mouse-position',
1348 : target: document.getElementById('mouse-position'),
1349 : undefinedHTML: ' '
1350 : });
1351 : var map = new ol.Map({
1352 : controls: ol.control.defaults.defaults().extend([mousePositionControl]),
1353 42 : target: 'map',)raw");
1354 :
1355 29 : if (tms.identifier() == "GoogleMapsCompatible" ||
1356 8 : tms.identifier() == "WorldCRS84Quad")
1357 : {
1358 15 : s += R"raw(
1359 : layers: [
1360 : new ol.layer.Group({
1361 : title: 'Base maps',
1362 : layers: [
1363 : new ol.layer.Tile({
1364 : title: 'OpenStreetMap',
1365 : type: 'base',
1366 : visible: true,
1367 : source: new ol.source.OSM()
1368 : }),
1369 : ]
1370 : }),)raw";
1371 : }
1372 :
1373 21 : if (tms.identifier() == "GoogleMapsCompatible")
1374 : {
1375 13 : s += R"raw(new ol.layer.Group({
1376 : title: 'Overlay',
1377 : layers: [
1378 : new ol.layer.Tile({
1379 : title: 'Overlay',
1380 : // opacity: 0.7,
1381 : extent: [%(ominx)f, %(ominy)f,%(omaxx)f, %(omaxy)f],
1382 : source: new ol.source.XYZ({
1383 : attributions: '%(copyright)s',
1384 : minZoom: %(minzoom)d,
1385 : maxZoom: %(maxzoom)d,
1386 : url: './{z}/{x}/{%(sign_y)sy}.%(tileformat)s',
1387 : tileSize: [%(tile_size)d, %(tile_size)d]
1388 : })
1389 : }),
1390 : ]
1391 : }),)raw";
1392 : }
1393 8 : else if (tms.identifier() == "WorldCRS84Quad")
1394 : {
1395 2 : const double base_res = 180.0 / nTileSize;
1396 4 : std::string resolutions = "[";
1397 4 : for (int i = 0; i <= nMaxZoom; ++i)
1398 : {
1399 2 : if (i > 0)
1400 0 : resolutions += ",";
1401 2 : resolutions += CPLSPrintf(pszFmt, base_res / (1 << i));
1402 : }
1403 2 : resolutions += "]";
1404 2 : substs["resolutions"] = std::move(resolutions);
1405 :
1406 2 : if (bXYZ)
1407 : {
1408 1 : substs["origin"] = "[-180,90]";
1409 1 : substs["y_formula"] = "tileCoord[2]";
1410 : }
1411 : else
1412 : {
1413 1 : substs["origin"] = "[-180,-90]";
1414 1 : substs["y_formula"] = "- 1 - tileCoord[2]";
1415 : }
1416 :
1417 2 : s += R"raw(
1418 : new ol.layer.Group({
1419 : title: 'Overlay',
1420 : layers: [
1421 : new ol.layer.Tile({
1422 : title: 'Overlay',
1423 : // opacity: 0.7,
1424 : extent: [%(ominx)f, %(ominy)f,%(omaxx)f, %(omaxy)f],
1425 : source: new ol.source.TileImage({
1426 : attributions: '%(copyright)s',
1427 : projection: 'EPSG:4326',
1428 : minZoom: %(minzoom)d,
1429 : maxZoom: %(maxzoom)d,
1430 : tileGrid: new ol.tilegrid.TileGrid({
1431 : extent: [-180,-90,180,90],
1432 : origin: %(origin)s,
1433 : resolutions: %(resolutions)s,
1434 : tileSize: [%(tile_size)d, %(tile_size)d]
1435 : }),
1436 : tileUrlFunction: function(tileCoord) {
1437 : return ('./{z}/{x}/{y}.%(tileformat)s'
1438 : .replace('{z}', String(tileCoord[0]))
1439 : .replace('{x}', String(tileCoord[1]))
1440 : .replace('{y}', String(%(y_formula)s)));
1441 : },
1442 : })
1443 : }),
1444 : ]
1445 : }),)raw";
1446 : }
1447 : else
1448 : {
1449 12 : substs["maxres"] =
1450 12 : CPLSPrintf(pszFmt, tms.tileMatrixList()[nMinZoom].mResX);
1451 12 : std::string resolutions = "[";
1452 16 : for (int i = 0; i <= nMaxZoom; ++i)
1453 : {
1454 10 : if (i > 0)
1455 4 : resolutions += ",";
1456 10 : resolutions += CPLSPrintf(pszFmt, tms.tileMatrixList()[i].mResX);
1457 : }
1458 6 : resolutions += "]";
1459 6 : substs["resolutions"] = std::move(resolutions);
1460 :
1461 12 : std::string matrixsizes = "[";
1462 16 : for (int i = 0; i <= nMaxZoom; ++i)
1463 : {
1464 10 : if (i > 0)
1465 4 : matrixsizes += ",";
1466 : matrixsizes +=
1467 10 : CPLSPrintf("[%d,%d]", tms.tileMatrixList()[i].mMatrixWidth,
1468 20 : tms.tileMatrixList()[i].mMatrixHeight);
1469 : }
1470 6 : matrixsizes += "]";
1471 6 : substs["matrixsizes"] = std::move(matrixsizes);
1472 :
1473 6 : double dfTopLeftX = tms.tileMatrixList()[0].mTopLeftX;
1474 6 : double dfTopLeftY = tms.tileMatrixList()[0].mTopLeftY;
1475 6 : if (bInvertAxisTMS)
1476 0 : std::swap(dfTopLeftX, dfTopLeftY);
1477 :
1478 6 : if (bXYZ)
1479 : {
1480 10 : substs["origin"] =
1481 10 : CPLSPrintf("[%.17g,%.17g]", dfTopLeftX, dfTopLeftY);
1482 5 : substs["y_formula"] = "tileCoord[2]";
1483 : }
1484 : else
1485 : {
1486 2 : substs["origin"] = CPLSPrintf(
1487 : "[%.17g,%.17g]", dfTopLeftX,
1488 1 : dfTopLeftY - tms.tileMatrixList()[0].mResY *
1489 3 : tms.tileMatrixList()[0].mTileHeight);
1490 1 : substs["y_formula"] = "- 1 - tileCoord[2]";
1491 : }
1492 :
1493 12 : substs["tilegrid_extent"] =
1494 : CPLSPrintf("[%.17g,%.17g,%.17g,%.17g]", dfTopLeftX,
1495 6 : dfTopLeftY - tms.tileMatrixList()[0].mMatrixHeight *
1496 6 : tms.tileMatrixList()[0].mResY *
1497 6 : tms.tileMatrixList()[0].mTileHeight,
1498 6 : dfTopLeftX + tms.tileMatrixList()[0].mMatrixWidth *
1499 6 : tms.tileMatrixList()[0].mResX *
1500 6 : tms.tileMatrixList()[0].mTileWidth,
1501 24 : dfTopLeftY);
1502 :
1503 6 : s += R"raw(
1504 : layers: [
1505 : new ol.layer.Group({
1506 : title: 'Overlay',
1507 : layers: [
1508 : new ol.layer.Tile({
1509 : title: 'Overlay',
1510 : // opacity: 0.7,
1511 : extent: [%(ominx)f, %(ominy)f,%(omaxx)f, %(omaxy)f],
1512 : source: new ol.source.TileImage({
1513 : attributions: '%(copyright)s',
1514 : minZoom: %(minzoom)d,
1515 : maxZoom: %(maxzoom)d,
1516 : tileGrid: new ol.tilegrid.TileGrid({
1517 : extent: %(tilegrid_extent)s,
1518 : origin: %(origin)s,
1519 : resolutions: %(resolutions)s,
1520 : sizes: %(matrixsizes)s,
1521 : tileSize: [%(tile_size)d, %(tile_size)d]
1522 : }),
1523 : tileUrlFunction: function(tileCoord) {
1524 : return ('./{z}/{x}/{y}.%(tileformat)s'
1525 : .replace('{z}', String(tileCoord[0]))
1526 : .replace('{x}', String(tileCoord[1]))
1527 : .replace('{y}', String(%(y_formula)s)));
1528 : },
1529 : })
1530 : }),
1531 : ]
1532 : }),)raw";
1533 : }
1534 :
1535 21 : s += R"raw(
1536 : ],
1537 : view: new ol.View({
1538 : center: [%(center_x)f, %(center_y)f],)raw";
1539 :
1540 29 : if (tms.identifier() == "GoogleMapsCompatible" ||
1541 8 : tms.identifier() == "WorldCRS84Quad")
1542 : {
1543 15 : substs["view_zoom"] = substs["minzoom"];
1544 15 : if (tms.identifier() == "WorldCRS84Quad")
1545 : {
1546 2 : substs["view_zoom"] = CPLSPrintf("%d", nMinZoom + 1);
1547 : }
1548 :
1549 15 : s += R"raw(
1550 : zoom: %(view_zoom)d,)raw";
1551 : }
1552 : else
1553 : {
1554 6 : s += R"raw(
1555 : resolution: %(maxres)f,)raw";
1556 : }
1557 :
1558 21 : if (tms.identifier() == "WorldCRS84Quad")
1559 : {
1560 2 : s += R"raw(
1561 : projection: 'EPSG:4326',)raw";
1562 : }
1563 19 : else if (!oSRS_TMS.IsEmpty() && tms.identifier() != "GoogleMapsCompatible")
1564 : {
1565 5 : const char *pszAuthName = oSRS_TMS.GetAuthorityName(nullptr);
1566 5 : const char *pszAuthCode = oSRS_TMS.GetAuthorityCode(nullptr);
1567 5 : if (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
1568 : {
1569 3 : substs["epsg_code"] = pszAuthCode;
1570 3 : if (oSRS_TMS.IsGeographic())
1571 : {
1572 1 : substs["units"] = "deg";
1573 : }
1574 : else
1575 : {
1576 2 : const char *pszUnits = "";
1577 2 : if (oSRS_TMS.GetLinearUnits(&pszUnits) == 1.0)
1578 2 : substs["units"] = "m";
1579 : else
1580 0 : substs["units"] = pszUnits;
1581 : }
1582 3 : s += R"raw(
1583 : projection: new ol.proj.Projection({code: 'EPSG:%(epsg_code)s', units:'%(units)s'}),)raw";
1584 : }
1585 : }
1586 :
1587 21 : s += R"raw(
1588 : })
1589 : });)raw";
1590 :
1591 29 : if (tms.identifier() == "GoogleMapsCompatible" ||
1592 8 : tms.identifier() == "WorldCRS84Quad")
1593 : {
1594 15 : s += R"raw(
1595 : map.addControl(new ol.control.LayerSwitcher());)raw";
1596 : }
1597 :
1598 21 : s += R"raw(
1599 : </script>
1600 : </body>
1601 : </html>)raw";
1602 :
1603 21 : ApplySubstitutions(s, substs);
1604 :
1605 21 : VSILFILE *f = VSIFOpenL(
1606 42 : CPLFormFilenameSafe(osDirectory.c_str(), "openlayers.html", nullptr)
1607 : .c_str(),
1608 : "wb");
1609 21 : if (f)
1610 : {
1611 21 : VSIFWriteL(s.data(), 1, s.size(), f);
1612 21 : VSIFCloseL(f);
1613 : }
1614 21 : }
1615 :
1616 : /************************************************************************/
1617 : /* GetTileBoundingBox() */
1618 : /************************************************************************/
1619 :
1620 6 : static void GetTileBoundingBox(int nTileX, int nTileY, int nTileZ,
1621 : const gdal::TileMatrixSet *poTMS,
1622 : bool bInvertAxisTMS,
1623 : OGRCoordinateTransformation *poCTToWGS84,
1624 : double &dfTLX, double &dfTLY, double &dfTRX,
1625 : double &dfTRY, double &dfLLX, double &dfLLY,
1626 : double &dfLRX, double &dfLRY)
1627 : {
1628 : gdal::TileMatrixSet::TileMatrix tileMatrix =
1629 12 : poTMS->tileMatrixList()[nTileZ];
1630 6 : if (bInvertAxisTMS)
1631 0 : std::swap(tileMatrix.mTopLeftX, tileMatrix.mTopLeftY);
1632 :
1633 6 : dfTLX = tileMatrix.mTopLeftX +
1634 6 : nTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
1635 6 : dfTLY = tileMatrix.mTopLeftY -
1636 6 : nTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
1637 6 : poCTToWGS84->Transform(1, &dfTLX, &dfTLY);
1638 :
1639 6 : dfTRX = tileMatrix.mTopLeftX +
1640 6 : (nTileX + 1) * tileMatrix.mResX * tileMatrix.mTileWidth;
1641 6 : dfTRY = tileMatrix.mTopLeftY -
1642 6 : nTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
1643 6 : poCTToWGS84->Transform(1, &dfTRX, &dfTRY);
1644 :
1645 6 : dfLLX = tileMatrix.mTopLeftX +
1646 6 : nTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
1647 6 : dfLLY = tileMatrix.mTopLeftY -
1648 6 : (nTileY + 1) * tileMatrix.mResY * tileMatrix.mTileHeight;
1649 6 : poCTToWGS84->Transform(1, &dfLLX, &dfLLY);
1650 :
1651 6 : dfLRX = tileMatrix.mTopLeftX +
1652 6 : (nTileX + 1) * tileMatrix.mResX * tileMatrix.mTileWidth;
1653 6 : dfLRY = tileMatrix.mTopLeftY -
1654 6 : (nTileY + 1) * tileMatrix.mResY * tileMatrix.mTileHeight;
1655 6 : poCTToWGS84->Transform(1, &dfLRX, &dfLRY);
1656 6 : }
1657 :
1658 : /************************************************************************/
1659 : /* GenerateKML() */
1660 : /************************************************************************/
1661 :
1662 : namespace
1663 : {
1664 : struct TileCoordinates
1665 : {
1666 : int nTileX = 0;
1667 : int nTileY = 0;
1668 : int nTileZ = 0;
1669 : };
1670 : } // namespace
1671 :
1672 5 : static void GenerateKML(const std::string &osDirectory,
1673 : const std::string &osTitle, int nTileX, int nTileY,
1674 : int nTileZ, int nTileSize,
1675 : const std::string &osExtension,
1676 : const std::string &osURL,
1677 : const gdal::TileMatrixSet *poTMS, bool bInvertAxisTMS,
1678 : const std::string &convention,
1679 : OGRCoordinateTransformation *poCTToWGS84,
1680 : const std::vector<TileCoordinates> &children)
1681 : {
1682 10 : std::map<std::string, std::string> substs;
1683 :
1684 5 : const bool bIsTileKML = nTileX >= 0;
1685 :
1686 : // For tests
1687 : const char *pszFmt =
1688 5 : atoi(CPLGetConfigOption("GDAL_RASTER_TILE_KML_PREC", "14")) == 10
1689 : ? "%.10f"
1690 5 : : "%.14f";
1691 :
1692 5 : substs["tx"] = CPLSPrintf("%d", nTileX);
1693 5 : substs["tz"] = CPLSPrintf("%d", nTileZ);
1694 5 : substs["tileformat"] = osExtension;
1695 5 : substs["minlodpixels"] = CPLSPrintf("%d", nTileSize / 2);
1696 10 : substs["maxlodpixels"] =
1697 10 : children.empty() ? "-1" : CPLSPrintf("%d", nTileSize * 8);
1698 :
1699 5 : double dfTLX = 0;
1700 5 : double dfTLY = 0;
1701 5 : double dfTRX = 0;
1702 5 : double dfTRY = 0;
1703 5 : double dfLLX = 0;
1704 5 : double dfLLY = 0;
1705 5 : double dfLRX = 0;
1706 5 : double dfLRY = 0;
1707 :
1708 5 : int nFileY = -1;
1709 5 : if (!bIsTileKML)
1710 : {
1711 2 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1712 2 : substs["xml_escaped_title"] = pszStr;
1713 2 : CPLFree(pszStr);
1714 : }
1715 : else
1716 : {
1717 3 : nFileY = GetFileY(nTileY, poTMS->tileMatrixList()[nTileZ], convention);
1718 3 : substs["realtiley"] = CPLSPrintf("%d", nFileY);
1719 6 : substs["xml_escaped_title"] =
1720 6 : CPLSPrintf("%d/%d/%d.kml", nTileZ, nTileX, nFileY);
1721 :
1722 3 : GetTileBoundingBox(nTileX, nTileY, nTileZ, poTMS, bInvertAxisTMS,
1723 : poCTToWGS84, dfTLX, dfTLY, dfTRX, dfTRY, dfLLX,
1724 : dfLLY, dfLRX, dfLRY);
1725 : }
1726 :
1727 10 : substs["drawOrder"] = CPLSPrintf("%d", nTileX == 0 ? 2 * nTileZ + 1
1728 4 : : nTileX > 0 ? 2 * nTileZ
1729 14 : : 0);
1730 :
1731 5 : substs["url"] = osURL.empty() && bIsTileKML ? "../../" : "";
1732 :
1733 5 : const bool bIsRectangle =
1734 5 : (dfTLX == dfLLX && dfTRX == dfLRX && dfTLY == dfTRY && dfLLY == dfLRY);
1735 5 : const bool bUseGXNamespace = bIsTileKML && !bIsRectangle;
1736 :
1737 10 : substs["xmlns_gx"] = bUseGXNamespace
1738 : ? " xmlns:gx=\"http://www.google.com/kml/ext/2.2\""
1739 10 : : "";
1740 :
1741 : CPLString s(R"raw(<?xml version="1.0" encoding="utf-8"?>
1742 : <kml xmlns="http://www.opengis.net/kml/2.2"%(xmlns_gx)s>
1743 : <Document>
1744 : <name>%(xml_escaped_title)s</name>
1745 : <description></description>
1746 : <Style>
1747 : <ListStyle id="hideChildren">
1748 : <listItemType>checkHideChildren</listItemType>
1749 : </ListStyle>
1750 : </Style>
1751 10 : )raw");
1752 5 : ApplySubstitutions(s, substs);
1753 :
1754 5 : if (bIsTileKML)
1755 : {
1756 : CPLString s2(R"raw( <Region>
1757 : <LatLonAltBox>
1758 : <north>%(north)f</north>
1759 : <south>%(south)f</south>
1760 : <east>%(east)f</east>
1761 : <west>%(west)f</west>
1762 : </LatLonAltBox>
1763 : <Lod>
1764 : <minLodPixels>%(minlodpixels)d</minLodPixels>
1765 : <maxLodPixels>%(maxlodpixels)d</maxLodPixels>
1766 : </Lod>
1767 : </Region>
1768 : <GroundOverlay>
1769 : <drawOrder>%(drawOrder)d</drawOrder>
1770 : <Icon>
1771 : <href>%(realtiley)d.%(tileformat)s</href>
1772 : </Icon>
1773 : <LatLonBox>
1774 : <north>%(north)f</north>
1775 : <south>%(south)f</south>
1776 : <east>%(east)f</east>
1777 : <west>%(west)f</west>
1778 : </LatLonBox>
1779 6 : )raw");
1780 :
1781 3 : if (!bIsRectangle)
1782 : {
1783 : s2 +=
1784 1 : R"raw( <gx:LatLonQuad><coordinates>%(LLX)f,%(LLY)f %(LRX)f,%(LRY)f %(TRX)f,%(TRY)f %(TLX)f,%(TLY)f</coordinates></gx:LatLonQuad>
1785 : )raw";
1786 : }
1787 :
1788 3 : s2 += R"raw( </GroundOverlay>
1789 : )raw";
1790 3 : substs["north"] = CPLSPrintf(pszFmt, std::max(dfTLY, dfTRY));
1791 3 : substs["south"] = CPLSPrintf(pszFmt, std::min(dfLLY, dfLRY));
1792 3 : substs["east"] = CPLSPrintf(pszFmt, std::max(dfTRX, dfLRX));
1793 3 : substs["west"] = CPLSPrintf(pszFmt, std::min(dfLLX, dfTLX));
1794 :
1795 3 : if (!bIsRectangle)
1796 : {
1797 1 : substs["TLX"] = CPLSPrintf(pszFmt, dfTLX);
1798 1 : substs["TLY"] = CPLSPrintf(pszFmt, dfTLY);
1799 1 : substs["TRX"] = CPLSPrintf(pszFmt, dfTRX);
1800 1 : substs["TRY"] = CPLSPrintf(pszFmt, dfTRY);
1801 1 : substs["LRX"] = CPLSPrintf(pszFmt, dfLRX);
1802 1 : substs["LRY"] = CPLSPrintf(pszFmt, dfLRY);
1803 1 : substs["LLX"] = CPLSPrintf(pszFmt, dfLLX);
1804 1 : substs["LLY"] = CPLSPrintf(pszFmt, dfLLY);
1805 : }
1806 :
1807 3 : ApplySubstitutions(s2, substs);
1808 3 : s += s2;
1809 : }
1810 :
1811 8 : for (const auto &child : children)
1812 : {
1813 3 : substs["tx"] = CPLSPrintf("%d", child.nTileX);
1814 3 : substs["tz"] = CPLSPrintf("%d", child.nTileZ);
1815 6 : substs["realtiley"] = CPLSPrintf(
1816 3 : "%d", GetFileY(child.nTileY, poTMS->tileMatrixList()[child.nTileZ],
1817 6 : convention));
1818 :
1819 3 : GetTileBoundingBox(child.nTileX, child.nTileY, child.nTileZ, poTMS,
1820 : bInvertAxisTMS, poCTToWGS84, dfTLX, dfTLY, dfTRX,
1821 : dfTRY, dfLLX, dfLLY, dfLRX, dfLRY);
1822 :
1823 : CPLString s2(R"raw( <NetworkLink>
1824 : <name>%(tz)d/%(tx)d/%(realtiley)d.%(tileformat)s</name>
1825 : <Region>
1826 : <LatLonAltBox>
1827 : <north>%(north)f</north>
1828 : <south>%(south)f</south>
1829 : <east>%(east)f</east>
1830 : <west>%(west)f</west>
1831 : </LatLonAltBox>
1832 : <Lod>
1833 : <minLodPixels>%(minlodpixels)d</minLodPixels>
1834 : <maxLodPixels>-1</maxLodPixels>
1835 : </Lod>
1836 : </Region>
1837 : <Link>
1838 : <href>%(url)s%(tz)d/%(tx)d/%(realtiley)d.kml</href>
1839 : <viewRefreshMode>onRegion</viewRefreshMode>
1840 : <viewFormat/>
1841 : </Link>
1842 : </NetworkLink>
1843 6 : )raw");
1844 3 : substs["north"] = CPLSPrintf(pszFmt, std::max(dfTLY, dfTRY));
1845 3 : substs["south"] = CPLSPrintf(pszFmt, std::min(dfLLY, dfLRY));
1846 3 : substs["east"] = CPLSPrintf(pszFmt, std::max(dfTRX, dfLRX));
1847 3 : substs["west"] = CPLSPrintf(pszFmt, std::min(dfLLX, dfTLX));
1848 3 : ApplySubstitutions(s2, substs);
1849 3 : s += s2;
1850 : }
1851 :
1852 5 : s += R"raw(</Document>
1853 : </kml>)raw";
1854 :
1855 10 : std::string osFilename(osDirectory);
1856 5 : if (!bIsTileKML)
1857 : {
1858 : osFilename =
1859 2 : CPLFormFilenameSafe(osFilename.c_str(), "doc.kml", nullptr);
1860 : }
1861 : else
1862 : {
1863 6 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
1864 3 : CPLSPrintf("%d", nTileZ), nullptr);
1865 6 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
1866 3 : CPLSPrintf("%d", nTileX), nullptr);
1867 6 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
1868 3 : CPLSPrintf("%d.kml", nFileY), nullptr);
1869 : }
1870 :
1871 5 : VSILFILE *f = VSIFOpenL(osFilename.c_str(), "wb");
1872 5 : if (f)
1873 : {
1874 5 : VSIFWriteL(s.data(), 1, s.size(), f);
1875 5 : VSIFCloseL(f);
1876 : }
1877 5 : }
1878 :
1879 : namespace
1880 : {
1881 :
1882 : /************************************************************************/
1883 : /* ResourceManager */
1884 : /************************************************************************/
1885 :
1886 : // Generic cache managing resources
1887 : template <class Resource> class ResourceManager /* non final */
1888 : {
1889 : public:
1890 94 : virtual ~ResourceManager() = default;
1891 :
1892 270 : std::unique_ptr<Resource> AcquireResources()
1893 : {
1894 542 : std::lock_guard oLock(m_oMutex);
1895 272 : if (!m_oResources.empty())
1896 : {
1897 288 : auto ret = std::move(m_oResources.back());
1898 144 : m_oResources.pop_back();
1899 144 : return ret;
1900 : }
1901 :
1902 128 : return CreateResources();
1903 : }
1904 :
1905 256 : void ReleaseResources(std::unique_ptr<Resource> resources)
1906 : {
1907 512 : std::lock_guard oLock(m_oMutex);
1908 256 : m_oResources.push_back(std::move(resources));
1909 256 : }
1910 :
1911 16 : void SetError()
1912 : {
1913 32 : std::lock_guard oLock(m_oMutex);
1914 16 : if (m_errorMsg.empty())
1915 1 : m_errorMsg = CPLGetLastErrorMsg();
1916 16 : }
1917 :
1918 33 : const std::string &GetErrorMsg() const
1919 : {
1920 33 : std::lock_guard oLock(m_oMutex);
1921 66 : return m_errorMsg;
1922 : }
1923 :
1924 : protected:
1925 : virtual std::unique_ptr<Resource> CreateResources() = 0;
1926 :
1927 : private:
1928 : mutable std::mutex m_oMutex{};
1929 : std::vector<std::unique_ptr<Resource>> m_oResources{};
1930 : std::string m_errorMsg{};
1931 : };
1932 :
1933 : /************************************************************************/
1934 : /* PerThreadMaxZoomResources */
1935 : /************************************************************************/
1936 :
1937 : // Per-thread resources for generation of tiles at full resolution
1938 : struct PerThreadMaxZoomResources
1939 : {
1940 : struct GDALDatasetReleaser
1941 : {
1942 112 : void operator()(GDALDataset *poDS)
1943 : {
1944 112 : if (poDS)
1945 112 : poDS->ReleaseRef();
1946 112 : }
1947 : };
1948 :
1949 : std::unique_ptr<GDALDataset, GDALDatasetReleaser> poSrcDS{};
1950 : std::vector<GByte> dstBuffer{};
1951 : std::unique_ptr<FakeMaxZoomDataset> poFakeMaxZoomDS{};
1952 : std::unique_ptr<void, decltype(&GDALDestroyTransformer)> poTransformer{
1953 : nullptr, GDALDestroyTransformer};
1954 : std::unique_ptr<GDALWarpOperation> poWO{};
1955 : };
1956 :
1957 : /************************************************************************/
1958 : /* PerThreadMaxZoomResourceManager */
1959 : /************************************************************************/
1960 :
1961 : // Manage a cache of PerThreadMaxZoomResources instances
1962 : class PerThreadMaxZoomResourceManager final
1963 : : public ResourceManager<PerThreadMaxZoomResources>
1964 : {
1965 : public:
1966 59 : PerThreadMaxZoomResourceManager(GDALDataset *poSrcDS,
1967 : const GDALWarpOptions *psWO,
1968 : void *pTransformerArg,
1969 : const FakeMaxZoomDataset &oFakeMaxZoomDS,
1970 : size_t nBufferSize)
1971 59 : : m_poSrcDS(poSrcDS), m_psWOSource(psWO),
1972 : m_pTransformerArg(pTransformerArg), m_oFakeMaxZoomDS(oFakeMaxZoomDS),
1973 59 : m_nBufferSize(nBufferSize)
1974 : {
1975 59 : }
1976 :
1977 : protected:
1978 112 : std::unique_ptr<PerThreadMaxZoomResources> CreateResources() override
1979 : {
1980 224 : auto ret = std::make_unique<PerThreadMaxZoomResources>();
1981 :
1982 112 : ret->poSrcDS.reset(GDALGetThreadSafeDataset(m_poSrcDS, GDAL_OF_RASTER));
1983 112 : if (!ret->poSrcDS)
1984 0 : return nullptr;
1985 :
1986 : try
1987 : {
1988 112 : ret->dstBuffer.resize(m_nBufferSize);
1989 : }
1990 0 : catch (const std::exception &)
1991 : {
1992 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
1993 : "Out of memory allocating temporary buffer");
1994 0 : return nullptr;
1995 : }
1996 :
1997 112 : ret->poFakeMaxZoomDS = m_oFakeMaxZoomDS.Clone(ret->dstBuffer);
1998 :
1999 112 : ret->poTransformer.reset(GDALCloneTransformer(m_pTransformerArg));
2000 112 : if (!ret->poTransformer)
2001 0 : return nullptr;
2002 :
2003 : auto psWO =
2004 : std::unique_ptr<GDALWarpOptions, decltype(&GDALDestroyWarpOptions)>(
2005 224 : GDALCloneWarpOptions(m_psWOSource), GDALDestroyWarpOptions);
2006 112 : if (!psWO)
2007 0 : return nullptr;
2008 :
2009 112 : psWO->hSrcDS = GDALDataset::ToHandle(ret->poSrcDS.get());
2010 112 : psWO->hDstDS = GDALDataset::ToHandle(ret->poFakeMaxZoomDS.get());
2011 112 : psWO->pTransformerArg = ret->poTransformer.get();
2012 112 : psWO->pfnTransformer = m_psWOSource->pfnTransformer;
2013 :
2014 112 : ret->poWO = std::make_unique<GDALWarpOperation>();
2015 112 : if (ret->poWO->Initialize(psWO.get()) != CE_None)
2016 0 : return nullptr;
2017 :
2018 112 : return ret;
2019 : }
2020 :
2021 : private:
2022 : GDALDataset *const m_poSrcDS;
2023 : const GDALWarpOptions *const m_psWOSource;
2024 : void *const m_pTransformerArg;
2025 : const FakeMaxZoomDataset &m_oFakeMaxZoomDS;
2026 : const size_t m_nBufferSize;
2027 :
2028 : CPL_DISALLOW_COPY_ASSIGN(PerThreadMaxZoomResourceManager)
2029 : };
2030 :
2031 : /************************************************************************/
2032 : /* PerThreadLowerZoomResources */
2033 : /************************************************************************/
2034 :
2035 : // Per-thread resources for generation of tiles at zoom level < max
2036 : struct PerThreadLowerZoomResources
2037 : {
2038 : std::unique_ptr<GDALDataset> poSrcDS{};
2039 : };
2040 :
2041 : /************************************************************************/
2042 : /* PerThreadLowerZoomResourceManager */
2043 : /************************************************************************/
2044 :
2045 : // Manage a cache of PerThreadLowerZoomResources instances
2046 : class PerThreadLowerZoomResourceManager final
2047 : : public ResourceManager<PerThreadLowerZoomResources>
2048 : {
2049 : public:
2050 35 : explicit PerThreadLowerZoomResourceManager(const MosaicDataset &oSrcDS)
2051 35 : : m_oSrcDS(oSrcDS)
2052 : {
2053 35 : }
2054 :
2055 : protected:
2056 16 : std::unique_ptr<PerThreadLowerZoomResources> CreateResources() override
2057 : {
2058 16 : auto ret = std::make_unique<PerThreadLowerZoomResources>();
2059 16 : ret->poSrcDS = m_oSrcDS.Clone();
2060 16 : return ret;
2061 : }
2062 :
2063 : private:
2064 : const MosaicDataset &m_oSrcDS;
2065 : };
2066 :
2067 : } // namespace
2068 :
2069 : /************************************************************************/
2070 : /* GDALRasterTileAlgorithm::RunImpl() */
2071 : /************************************************************************/
2072 :
2073 87 : bool GDALRasterTileAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
2074 : void *pProgressData)
2075 : {
2076 87 : auto poSrcDS = m_dataset.GetDatasetRef();
2077 87 : CPLAssert(poSrcDS);
2078 87 : const int nSrcWidth = poSrcDS->GetRasterXSize();
2079 87 : const int nSrcHeight = poSrcDS->GetRasterYSize();
2080 87 : if (poSrcDS->GetRasterCount() == 0 || nSrcWidth == 0 || nSrcHeight == 0)
2081 : {
2082 1 : ReportError(CE_Failure, CPLE_AppDefined, "Invalid source dataset");
2083 1 : return false;
2084 : }
2085 :
2086 86 : if (m_resampling == "near")
2087 1 : m_resampling = "nearest";
2088 86 : if (m_overviewResampling == "near")
2089 1 : m_overviewResampling = "nearest";
2090 85 : else if (m_overviewResampling.empty())
2091 80 : m_overviewResampling = m_resampling;
2092 :
2093 172 : CPLStringList aosWarpOptions;
2094 86 : if (!m_excludedValues.empty() || m_nodataValuesPctThreshold < 100)
2095 : {
2096 : aosWarpOptions.SetNameValue(
2097 : "NODATA_VALUES_PCT_THRESHOLD",
2098 3 : CPLSPrintf("%g", m_nodataValuesPctThreshold));
2099 3 : if (!m_excludedValues.empty())
2100 : {
2101 : aosWarpOptions.SetNameValue("EXCLUDED_VALUES",
2102 1 : m_excludedValues.c_str());
2103 : aosWarpOptions.SetNameValue(
2104 : "EXCLUDED_VALUES_PCT_THRESHOLD",
2105 1 : CPLSPrintf("%g", m_excludedValuesPctThreshold));
2106 : }
2107 : }
2108 :
2109 86 : if (poSrcDS->GetRasterBand(1)->GetColorInterpretation() ==
2110 90 : GCI_PaletteIndex &&
2111 4 : ((m_resampling != "nearest" && m_resampling != "mode") ||
2112 1 : (m_overviewResampling != "nearest" && m_overviewResampling != "mode")))
2113 : {
2114 1 : ReportError(CE_Failure, CPLE_NotSupported,
2115 : "Datasets with color table not supported with non-nearest "
2116 : "or non-mode resampling. Run 'gdal raster "
2117 : "color-map' before or set the 'resampling' argument to "
2118 : "'nearest' or 'mode'.");
2119 1 : return false;
2120 : }
2121 :
2122 85 : const auto eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
2123 : auto poDstDriver =
2124 85 : GetGDALDriverManager()->GetDriverByName(m_outputFormat.c_str());
2125 85 : if (!poDstDriver)
2126 : {
2127 1 : ReportError(CE_Failure, CPLE_AppDefined,
2128 : "Invalid value for argument 'output-format'. Driver '%s' "
2129 : "does not exist",
2130 : m_outputFormat.c_str());
2131 1 : return false;
2132 : }
2133 :
2134 84 : if (m_outputFormat == "PNG")
2135 : {
2136 60 : if (poSrcDS->GetRasterCount() > 4)
2137 : {
2138 1 : ReportError(CE_Failure, CPLE_NotSupported,
2139 : "Only up to 4 bands supported for PNG.");
2140 1 : return false;
2141 : }
2142 59 : if (eSrcDT != GDT_Byte && eSrcDT != GDT_UInt16)
2143 : {
2144 1 : ReportError(CE_Failure, CPLE_NotSupported,
2145 : "Only Byte and UInt16 data types supported for PNG.");
2146 1 : return false;
2147 : }
2148 : }
2149 24 : else if (m_outputFormat == "JPEG")
2150 : {
2151 5 : if (poSrcDS->GetRasterCount() > 4)
2152 : {
2153 1 : ReportError(
2154 : CE_Failure, CPLE_NotSupported,
2155 : "Only up to 4 bands supported for JPEG (with alpha ignored).");
2156 1 : return false;
2157 : }
2158 8 : const bool bUInt16Supported = strstr(
2159 4 : poDstDriver->GetMetadataItem(GDAL_DMD_CREATIONDATATYPES), "UInt16");
2160 4 : if (eSrcDT != GDT_Byte && !(eSrcDT == GDT_UInt16 && bUInt16Supported))
2161 : {
2162 1 : ReportError(
2163 : CE_Failure, CPLE_NotSupported,
2164 : bUInt16Supported
2165 : ? "Only Byte and UInt16 data types supported for JPEG."
2166 : : "Only Byte data type supported for JPEG.");
2167 1 : return false;
2168 : }
2169 3 : if (eSrcDT == GDT_UInt16)
2170 : {
2171 3 : if (const char *pszNBITS =
2172 6 : poSrcDS->GetRasterBand(1)->GetMetadataItem(
2173 3 : "NBITS", "IMAGE_STRUCTURE"))
2174 : {
2175 1 : if (atoi(pszNBITS) > 12)
2176 : {
2177 1 : ReportError(CE_Failure, CPLE_NotSupported,
2178 : "JPEG output only supported up to 12 bits");
2179 1 : return false;
2180 : }
2181 : }
2182 : else
2183 : {
2184 2 : double adfMinMax[2] = {0, 0};
2185 2 : poSrcDS->GetRasterBand(1)->ComputeRasterMinMax(
2186 2 : /* bApproxOK = */ true, adfMinMax);
2187 2 : if (adfMinMax[1] >= (1 << 12))
2188 : {
2189 1 : ReportError(CE_Failure, CPLE_NotSupported,
2190 : "JPEG output only supported up to 12 bits");
2191 1 : return false;
2192 : }
2193 : }
2194 : }
2195 : }
2196 19 : else if (m_outputFormat == "WEBP")
2197 : {
2198 2 : if (poSrcDS->GetRasterCount() != 3 && poSrcDS->GetRasterCount() != 4)
2199 : {
2200 1 : ReportError(CE_Failure, CPLE_NotSupported,
2201 : "Only 3 or 4 bands supported for WEBP.");
2202 1 : return false;
2203 : }
2204 1 : if (eSrcDT != GDT_Byte)
2205 : {
2206 1 : ReportError(CE_Failure, CPLE_NotSupported,
2207 : "Only Byte data type supported for WEBP.");
2208 1 : return false;
2209 : }
2210 : }
2211 :
2212 : const char *pszExtensions =
2213 76 : poDstDriver->GetMetadataItem(GDAL_DMD_EXTENSIONS);
2214 76 : CPLAssert(pszExtensions && pszExtensions[0] != 0);
2215 : const CPLStringList aosExtensions(
2216 152 : CSLTokenizeString2(pszExtensions, " ", 0));
2217 76 : const char *pszExtension = aosExtensions[0];
2218 : std::array<double, 6> adfSrcGeoTransform;
2219 : const bool bHasSrcGT =
2220 76 : poSrcDS->GetGeoTransform(adfSrcGeoTransform.data()) == CE_None;
2221 74 : const bool bHasNorthUpSrcGT = bHasSrcGT && adfSrcGeoTransform[2] == 0 &&
2222 224 : adfSrcGeoTransform[4] == 0 &&
2223 74 : adfSrcGeoTransform[5] < 0;
2224 152 : OGRSpatialReference oSRS_TMS;
2225 :
2226 76 : if (m_tilingScheme == "raster")
2227 : {
2228 4 : if (const auto poSRS = poSrcDS->GetSpatialRef())
2229 3 : oSRS_TMS = *poSRS;
2230 : }
2231 : else
2232 : {
2233 1 : if (!bHasSrcGT && poSrcDS->GetGCPCount() == 0 &&
2234 74 : poSrcDS->GetMetadata("GEOLOCATION") == nullptr &&
2235 1 : poSrcDS->GetMetadata("RPC") == nullptr)
2236 : {
2237 1 : ReportError(CE_Failure, CPLE_NotSupported,
2238 : "Ungeoreferenced datasets are not supported, unless "
2239 : "'tiling-scheme' is set to 'raster'");
2240 1 : return false;
2241 : }
2242 :
2243 71 : if (poSrcDS->GetMetadata("GEOLOCATION") == nullptr &&
2244 71 : poSrcDS->GetMetadata("RPC") == nullptr &&
2245 143 : poSrcDS->GetSpatialRef() == nullptr &&
2246 1 : poSrcDS->GetGCPSpatialRef() == nullptr)
2247 : {
2248 1 : ReportError(CE_Failure, CPLE_NotSupported,
2249 : "Ungeoreferenced datasets are not supported, unless "
2250 : "'tiling-scheme' is set to 'raster'");
2251 1 : return false;
2252 : }
2253 : }
2254 :
2255 74 : if (m_copySrcMetadata)
2256 : {
2257 8 : CPLStringList aosMD(CSLDuplicate(poSrcDS->GetMetadata()));
2258 4 : const CPLStringList aosNewMD(m_metadata);
2259 8 : for (const auto [key, value] : cpl::IterateNameValue(aosNewMD))
2260 : {
2261 4 : aosMD.SetNameValue(key, value);
2262 : }
2263 4 : m_metadata = aosMD;
2264 : }
2265 :
2266 74 : std::array<double, 6> adfSrcGeoTransformModif{0, 1, 0, 0, 0, -1};
2267 :
2268 74 : if (m_tilingScheme == "mercator")
2269 1 : m_tilingScheme = "WebMercatorQuad";
2270 73 : else if (m_tilingScheme == "geodetic")
2271 1 : m_tilingScheme = "WorldCRS84Quad";
2272 72 : else if (m_tilingScheme == "raster")
2273 : {
2274 4 : if (m_tileSize == 0)
2275 4 : m_tileSize = 256;
2276 4 : if (m_maxZoomLevel < 0)
2277 : {
2278 3 : m_maxZoomLevel = static_cast<int>(std::ceil(std::log2(
2279 6 : std::max(1, std::max(nSrcWidth, nSrcHeight) / m_tileSize))));
2280 : }
2281 4 : if (bHasNorthUpSrcGT)
2282 : {
2283 3 : adfSrcGeoTransformModif = adfSrcGeoTransform;
2284 : }
2285 : }
2286 :
2287 : auto poTMS =
2288 74 : m_tilingScheme == "raster"
2289 : ? gdal::TileMatrixSet::createRaster(
2290 4 : nSrcWidth, nSrcHeight, m_tileSize, 1 + m_maxZoomLevel,
2291 8 : adfSrcGeoTransformModif[0], adfSrcGeoTransformModif[3],
2292 8 : adfSrcGeoTransformModif[1], -adfSrcGeoTransformModif[5],
2293 78 : oSRS_TMS.IsEmpty() ? std::string() : oSRS_TMS.exportToWkt())
2294 : : gdal::TileMatrixSet::parse(
2295 152 : m_mapTileMatrixIdentifierToScheme[m_tilingScheme].c_str());
2296 : // Enforced by SetChoices() on the m_tilingScheme argument
2297 74 : CPLAssert(poTMS && !poTMS->hasVariableMatrixWidth());
2298 :
2299 148 : CPLStringList aosTO;
2300 74 : if (m_tilingScheme == "raster")
2301 : {
2302 4 : aosTO.SetNameValue("SRC_METHOD", "GEOTRANSFORM");
2303 : }
2304 : else
2305 : {
2306 70 : CPL_IGNORE_RET_VAL(oSRS_TMS.SetFromUserInput(poTMS->crs().c_str()));
2307 70 : aosTO.SetNameValue("DST_SRS", oSRS_TMS.exportToWkt().c_str());
2308 : }
2309 :
2310 74 : const char *pszAuthName = oSRS_TMS.GetAuthorityName(nullptr);
2311 74 : const char *pszAuthCode = oSRS_TMS.GetAuthorityCode(nullptr);
2312 74 : const int nEPSGCode =
2313 73 : (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
2314 147 : ? atoi(pszAuthCode)
2315 : : 0;
2316 :
2317 : const bool bInvertAxisTMS =
2318 144 : m_tilingScheme != "raster" &&
2319 70 : (oSRS_TMS.EPSGTreatsAsLatLong() != FALSE ||
2320 70 : oSRS_TMS.EPSGTreatsAsNorthingEasting() != FALSE);
2321 :
2322 74 : oSRS_TMS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2323 :
2324 : std::unique_ptr<void, decltype(&GDALDestroyTransformer)> hTransformArg(
2325 148 : nullptr, GDALDestroyTransformer);
2326 :
2327 : // Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
2328 : // ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
2329 : // EPSG:3857.
2330 74 : std::unique_ptr<GDALDataset> poTmpDS;
2331 74 : bool bEPSG3857Adjust = false;
2332 74 : if (nEPSGCode == 3857 && bHasNorthUpSrcGT)
2333 : {
2334 48 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
2335 48 : if (poSrcSRS && poSrcSRS->IsGeographic())
2336 : {
2337 17 : double maxLat = adfSrcGeoTransform[3];
2338 : double minLat =
2339 17 : adfSrcGeoTransform[3] + nSrcHeight * adfSrcGeoTransform[5];
2340 : // Corresponds to the latitude of below MAX_GM
2341 17 : constexpr double MAX_LAT = 85.0511287798066;
2342 17 : bool bModified = false;
2343 17 : if (maxLat > MAX_LAT)
2344 : {
2345 16 : maxLat = MAX_LAT;
2346 16 : bModified = true;
2347 : }
2348 17 : if (minLat < -MAX_LAT)
2349 : {
2350 16 : minLat = -MAX_LAT;
2351 16 : bModified = true;
2352 : }
2353 17 : if (bModified)
2354 : {
2355 32 : CPLStringList aosOptions;
2356 16 : aosOptions.AddString("-of");
2357 16 : aosOptions.AddString("VRT");
2358 16 : aosOptions.AddString("-projwin");
2359 : aosOptions.AddString(
2360 16 : CPLSPrintf("%.17g", adfSrcGeoTransform[0]));
2361 16 : aosOptions.AddString(CPLSPrintf("%.17g", maxLat));
2362 : aosOptions.AddString(
2363 16 : CPLSPrintf("%.17g", adfSrcGeoTransform[0] +
2364 16 : nSrcWidth * adfSrcGeoTransform[1]));
2365 16 : aosOptions.AddString(CPLSPrintf("%.17g", minLat));
2366 : auto psOptions =
2367 16 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
2368 16 : poTmpDS.reset(GDALDataset::FromHandle(GDALTranslate(
2369 : "", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
2370 16 : GDALTranslateOptionsFree(psOptions);
2371 16 : if (poTmpDS)
2372 : {
2373 16 : bEPSG3857Adjust = true;
2374 16 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
2375 16 : GDALDataset::FromHandle(poTmpDS.get()), nullptr,
2376 16 : aosTO.List()));
2377 : }
2378 : }
2379 : }
2380 : }
2381 :
2382 : std::array<double, 6> adfDstGeoTransform;
2383 : double adfExtent[4];
2384 : int nXSize, nYSize;
2385 :
2386 : bool bSuggestOK;
2387 74 : if (m_tilingScheme == "raster")
2388 : {
2389 4 : bSuggestOK = true;
2390 4 : nXSize = nSrcWidth;
2391 4 : nYSize = nSrcHeight;
2392 4 : adfDstGeoTransform = adfSrcGeoTransformModif;
2393 4 : adfExtent[0] = adfDstGeoTransform[0];
2394 4 : adfExtent[1] =
2395 4 : adfDstGeoTransform[3] + nSrcHeight * adfDstGeoTransform[5];
2396 4 : adfExtent[2] =
2397 4 : adfDstGeoTransform[0] + nSrcWidth * adfDstGeoTransform[1];
2398 4 : adfExtent[3] = adfDstGeoTransform[3];
2399 : }
2400 : else
2401 : {
2402 70 : if (!hTransformArg)
2403 : {
2404 54 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
2405 54 : poSrcDS, nullptr, aosTO.List()));
2406 : }
2407 70 : if (!hTransformArg)
2408 : {
2409 1 : return false;
2410 : }
2411 69 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
2412 69 : bSuggestOK =
2413 138 : (GDALSuggestedWarpOutput2(
2414 : poSrcDS,
2415 69 : static_cast<GDALTransformerInfo *>(hTransformArg.get())
2416 : ->pfnTransform,
2417 : hTransformArg.get(), adfDstGeoTransform.data(), &nXSize,
2418 : &nYSize, adfExtent, 0) == CE_None);
2419 : }
2420 73 : if (!bSuggestOK)
2421 : {
2422 1 : ReportError(CE_Failure, CPLE_AppDefined,
2423 : "Cannot determine extent of raster in target CRS");
2424 1 : return false;
2425 : }
2426 :
2427 72 : poTmpDS.reset();
2428 :
2429 72 : if (bEPSG3857Adjust)
2430 : {
2431 16 : constexpr double SPHERICAL_RADIUS = 6378137.0;
2432 16 : constexpr double MAX_GM =
2433 : SPHERICAL_RADIUS * M_PI; // 20037508.342789244
2434 16 : double maxNorthing = adfDstGeoTransform[3];
2435 : double minNorthing =
2436 16 : adfDstGeoTransform[3] + adfDstGeoTransform[5] * nYSize;
2437 16 : bool bChanged = false;
2438 16 : if (maxNorthing > MAX_GM)
2439 : {
2440 13 : bChanged = true;
2441 13 : maxNorthing = MAX_GM;
2442 : }
2443 16 : if (minNorthing < -MAX_GM)
2444 : {
2445 13 : bChanged = true;
2446 13 : minNorthing = -MAX_GM;
2447 : }
2448 16 : if (bChanged)
2449 : {
2450 13 : adfDstGeoTransform[3] = maxNorthing;
2451 13 : nYSize = int(
2452 13 : (maxNorthing - minNorthing) / (-adfDstGeoTransform[5]) + 0.5);
2453 13 : adfExtent[1] = maxNorthing + nYSize * adfDstGeoTransform[5];
2454 13 : adfExtent[3] = maxNorthing;
2455 : }
2456 : }
2457 :
2458 72 : const auto &tileMatrixList = poTMS->tileMatrixList();
2459 72 : if (m_maxZoomLevel >= 0)
2460 : {
2461 23 : if (m_maxZoomLevel >= static_cast<int>(tileMatrixList.size()))
2462 : {
2463 1 : ReportError(CE_Failure, CPLE_AppDefined,
2464 : "max-zoom = %d is invalid. It must be in [0,%d] range",
2465 : m_maxZoomLevel,
2466 1 : static_cast<int>(tileMatrixList.size()) - 1);
2467 1 : return false;
2468 : }
2469 : }
2470 : else
2471 : {
2472 49 : const double dfComputedRes = adfDstGeoTransform[1];
2473 49 : double dfPrevRes = 0.0;
2474 49 : double dfRes = 0.0;
2475 49 : constexpr double EPSILON = 1e-8;
2476 :
2477 49 : if (m_minZoomLevel >= 0)
2478 17 : m_maxZoomLevel = m_minZoomLevel;
2479 : else
2480 32 : m_maxZoomLevel = 0;
2481 :
2482 164 : for (; m_maxZoomLevel < static_cast<int>(tileMatrixList.size());
2483 115 : m_maxZoomLevel++)
2484 : {
2485 163 : dfRes = tileMatrixList[m_maxZoomLevel].mResX;
2486 163 : if (dfComputedRes > dfRes ||
2487 115 : fabs(dfComputedRes - dfRes) / dfRes <= EPSILON)
2488 : break;
2489 115 : dfPrevRes = dfRes;
2490 : }
2491 49 : if (m_maxZoomLevel >= static_cast<int>(tileMatrixList.size()))
2492 : {
2493 1 : ReportError(CE_Failure, CPLE_AppDefined,
2494 : "Could not find an appropriate zoom level. Perhaps "
2495 : "min-zoom is too large?");
2496 1 : return false;
2497 : }
2498 :
2499 48 : if (m_maxZoomLevel > 0 && fabs(dfComputedRes - dfRes) / dfRes > EPSILON)
2500 : {
2501 : // Round to closest resolution
2502 43 : if (dfPrevRes / dfComputedRes < dfComputedRes / dfRes)
2503 23 : m_maxZoomLevel--;
2504 : }
2505 : }
2506 70 : if (m_minZoomLevel < 0)
2507 40 : m_minZoomLevel = m_maxZoomLevel;
2508 :
2509 140 : auto tileMatrix = tileMatrixList[m_maxZoomLevel];
2510 70 : int nMinTileX = 0;
2511 70 : int nMinTileY = 0;
2512 70 : int nMaxTileX = 0;
2513 70 : int nMaxTileY = 0;
2514 70 : bool bIntersects = false;
2515 70 : if (!GetTileIndices(tileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
2516 : nMinTileX, nMinTileY, nMaxTileX, nMaxTileY,
2517 70 : m_noIntersectionIsOK, bIntersects,
2518 : /* checkRasterOverflow = */ false))
2519 : {
2520 1 : return false;
2521 : }
2522 69 : if (!bIntersects)
2523 1 : return true;
2524 :
2525 : // Potentially restrict tiling to user specified coordinates
2526 68 : if (m_minTileX >= tileMatrix.mMatrixWidth)
2527 : {
2528 1 : ReportError(CE_Failure, CPLE_IllegalArg,
2529 : "'min-x' value must be in [0,%d] range",
2530 1 : tileMatrix.mMatrixWidth - 1);
2531 1 : return false;
2532 : }
2533 67 : if (m_maxTileX >= tileMatrix.mMatrixWidth)
2534 : {
2535 1 : ReportError(CE_Failure, CPLE_IllegalArg,
2536 : "'max-x' value must be in [0,%d] range",
2537 1 : tileMatrix.mMatrixWidth - 1);
2538 1 : return false;
2539 : }
2540 66 : if (m_minTileY >= tileMatrix.mMatrixHeight)
2541 : {
2542 1 : ReportError(CE_Failure, CPLE_IllegalArg,
2543 : "'min-y' value must be in [0,%d] range",
2544 1 : tileMatrix.mMatrixHeight - 1);
2545 1 : return false;
2546 : }
2547 65 : if (m_maxTileY >= tileMatrix.mMatrixHeight)
2548 : {
2549 1 : ReportError(CE_Failure, CPLE_IllegalArg,
2550 : "'max-y' value must be in [0,%d] range",
2551 1 : tileMatrix.mMatrixHeight - 1);
2552 1 : return false;
2553 : }
2554 :
2555 64 : if ((m_minTileX >= 0 && m_minTileX > nMaxTileX) ||
2556 62 : (m_minTileY >= 0 && m_minTileY > nMaxTileY) ||
2557 62 : (m_maxTileX >= 0 && m_maxTileX < nMinTileX) ||
2558 62 : (m_maxTileY >= 0 && m_maxTileY < nMinTileY))
2559 : {
2560 2 : ReportError(
2561 2 : m_noIntersectionIsOK ? CE_Warning : CE_Failure, CPLE_AppDefined,
2562 : "Dataset extent not intersecting specified min/max X/Y tile "
2563 : "coordinates");
2564 2 : return m_noIntersectionIsOK;
2565 : }
2566 62 : if (m_minTileX >= 0 && m_minTileX > nMinTileX)
2567 : {
2568 2 : nMinTileX = m_minTileX;
2569 2 : adfExtent[0] = tileMatrix.mTopLeftX +
2570 2 : nMinTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
2571 : }
2572 62 : if (m_minTileY >= 0 && m_minTileY > nMinTileY)
2573 : {
2574 1 : nMinTileY = m_minTileY;
2575 1 : adfExtent[3] = tileMatrix.mTopLeftY -
2576 1 : nMinTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
2577 : }
2578 62 : if (m_maxTileX >= 0 && m_maxTileX < nMaxTileX)
2579 : {
2580 2 : nMaxTileX = m_maxTileX;
2581 2 : adfExtent[2] = tileMatrix.mTopLeftX + (nMaxTileX + 1) *
2582 2 : tileMatrix.mResX *
2583 2 : tileMatrix.mTileWidth;
2584 : }
2585 62 : if (m_maxTileY >= 0 && m_maxTileY < nMaxTileY)
2586 : {
2587 2 : nMaxTileY = m_maxTileY;
2588 2 : adfExtent[1] = tileMatrix.mTopLeftY - (nMaxTileY + 1) *
2589 2 : tileMatrix.mResY *
2590 2 : tileMatrix.mTileHeight;
2591 : }
2592 :
2593 62 : if (nMaxTileX - nMinTileX + 1 > INT_MAX / tileMatrix.mTileWidth ||
2594 61 : nMaxTileY - nMinTileY + 1 > INT_MAX / tileMatrix.mTileHeight)
2595 : {
2596 1 : ReportError(CE_Failure, CPLE_AppDefined, "Too large zoom level");
2597 1 : return false;
2598 : }
2599 :
2600 122 : adfDstGeoTransform[0] = tileMatrix.mTopLeftX + nMinTileX *
2601 61 : tileMatrix.mResX *
2602 61 : tileMatrix.mTileWidth;
2603 61 : adfDstGeoTransform[1] = tileMatrix.mResX;
2604 61 : adfDstGeoTransform[2] = 0;
2605 122 : adfDstGeoTransform[3] = tileMatrix.mTopLeftY - nMinTileY *
2606 61 : tileMatrix.mResY *
2607 61 : tileMatrix.mTileHeight;
2608 61 : adfDstGeoTransform[4] = 0;
2609 61 : adfDstGeoTransform[5] = -tileMatrix.mResY;
2610 :
2611 : /* -------------------------------------------------------------------- */
2612 : /* Setup warp options. */
2613 : /* -------------------------------------------------------------------- */
2614 : std::unique_ptr<GDALWarpOptions, decltype(&GDALDestroyWarpOptions)> psWO(
2615 122 : GDALCreateWarpOptions(), GDALDestroyWarpOptions);
2616 :
2617 61 : psWO->papszWarpOptions = CSLSetNameValue(nullptr, "OPTIMIZE_SIZE", "YES");
2618 122 : psWO->papszWarpOptions =
2619 61 : CSLSetNameValue(psWO->papszWarpOptions, "SAMPLE_GRID", "YES");
2620 122 : psWO->papszWarpOptions =
2621 61 : CSLMerge(psWO->papszWarpOptions, aosWarpOptions.List());
2622 :
2623 61 : int bHasSrcNoData = false;
2624 : const double dfSrcNoDataValue =
2625 61 : poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasSrcNoData);
2626 :
2627 : const bool bLastSrcBandIsAlpha =
2628 93 : (poSrcDS->GetRasterCount() > 1 &&
2629 32 : poSrcDS->GetRasterBand(poSrcDS->GetRasterCount())
2630 32 : ->GetColorInterpretation() == GCI_AlphaBand);
2631 :
2632 61 : const bool bOutputSupportsAlpha = !EQUAL(m_outputFormat.c_str(), "JPEG");
2633 61 : const bool bOutputSupportsNoData = EQUAL(m_outputFormat.c_str(), "GTiff");
2634 61 : const bool bDstNoDataSpecified = GetArg("dst-nodata")->IsExplicitlySet();
2635 : const GDALColorTable *poColorTable =
2636 61 : poSrcDS->GetRasterBand(1)->GetColorTable();
2637 :
2638 61 : const bool bUserAskedForAlpha = m_addalpha;
2639 61 : if (!m_noalpha && !m_addalpha)
2640 : {
2641 71 : m_addalpha = !(bHasSrcNoData && bOutputSupportsNoData) &&
2642 71 : !bDstNoDataSpecified && poColorTable == nullptr;
2643 : }
2644 61 : m_addalpha &= bOutputSupportsAlpha;
2645 :
2646 61 : psWO->nBandCount = poSrcDS->GetRasterCount();
2647 61 : if (bLastSrcBandIsAlpha)
2648 : {
2649 5 : --psWO->nBandCount;
2650 5 : psWO->nSrcAlphaBand = poSrcDS->GetRasterCount();
2651 : }
2652 :
2653 61 : if (bHasSrcNoData)
2654 : {
2655 26 : psWO->padfSrcNoDataReal =
2656 13 : static_cast<double *>(CPLCalloc(psWO->nBandCount, sizeof(double)));
2657 32 : for (int i = 0; i < psWO->nBandCount; ++i)
2658 : {
2659 19 : psWO->padfSrcNoDataReal[i] = dfSrcNoDataValue;
2660 : }
2661 : }
2662 :
2663 61 : if ((bHasSrcNoData && !m_addalpha && bOutputSupportsNoData) ||
2664 : bDstNoDataSpecified)
2665 : {
2666 18 : psWO->padfDstNoDataReal =
2667 9 : static_cast<double *>(CPLCalloc(psWO->nBandCount, sizeof(double)));
2668 18 : for (int i = 0; i < psWO->nBandCount; ++i)
2669 : {
2670 9 : psWO->padfDstNoDataReal[i] =
2671 9 : bDstNoDataSpecified ? m_dstNoData : dfSrcNoDataValue;
2672 : }
2673 : }
2674 :
2675 61 : psWO->eWorkingDataType = eSrcDT;
2676 :
2677 61 : GDALGetWarpResampleAlg(m_resampling.c_str(), psWO->eResampleAlg);
2678 :
2679 : /* -------------------------------------------------------------------- */
2680 : /* Setup band mapping. */
2681 : /* -------------------------------------------------------------------- */
2682 :
2683 122 : psWO->panSrcBands =
2684 61 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
2685 122 : psWO->panDstBands =
2686 61 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
2687 :
2688 65719 : for (int i = 0; i < psWO->nBandCount; i++)
2689 : {
2690 65658 : psWO->panSrcBands[i] = i + 1;
2691 65658 : psWO->panDstBands[i] = i + 1;
2692 : }
2693 :
2694 61 : if (m_addalpha)
2695 48 : psWO->nDstAlphaBand = psWO->nBandCount + 1;
2696 :
2697 : const int nDstBands =
2698 61 : psWO->nDstAlphaBand ? psWO->nDstAlphaBand : psWO->nBandCount;
2699 :
2700 122 : std::vector<GByte> dstBuffer;
2701 : const uint64_t dstBufferSize =
2702 61 : static_cast<uint64_t>(tileMatrix.mTileWidth) * tileMatrix.mTileHeight *
2703 61 : nDstBands * GDALGetDataTypeSizeBytes(psWO->eWorkingDataType);
2704 : const uint64_t nUsableRAM =
2705 61 : std::min<uint64_t>(INT_MAX, CPLGetUsablePhysicalRAM() / 4);
2706 61 : if (dstBufferSize <=
2707 61 : (nUsableRAM ? nUsableRAM : static_cast<uint64_t>(INT_MAX)))
2708 : {
2709 : try
2710 : {
2711 60 : dstBuffer.resize(static_cast<size_t>(dstBufferSize));
2712 : }
2713 0 : catch (const std::exception &)
2714 : {
2715 : }
2716 : }
2717 61 : if (dstBuffer.size() < dstBufferSize)
2718 : {
2719 1 : ReportError(CE_Failure, CPLE_AppDefined,
2720 : "Tile size and/or number of bands too large compared to "
2721 : "available RAM");
2722 1 : return false;
2723 : }
2724 :
2725 : FakeMaxZoomDataset oFakeMaxZoomDS(
2726 60 : (nMaxTileX - nMinTileX + 1) * tileMatrix.mTileWidth,
2727 60 : (nMaxTileY - nMinTileY + 1) * tileMatrix.mTileHeight, nDstBands,
2728 60 : tileMatrix.mTileWidth, tileMatrix.mTileHeight, psWO->eWorkingDataType,
2729 180 : adfDstGeoTransform.data(), oSRS_TMS, dstBuffer);
2730 60 : CPL_IGNORE_RET_VAL(oFakeMaxZoomDS.GetSpatialRef());
2731 :
2732 60 : psWO->hSrcDS = GDALDataset::ToHandle(poSrcDS);
2733 60 : psWO->hDstDS = GDALDataset::ToHandle(&oFakeMaxZoomDS);
2734 :
2735 60 : std::unique_ptr<GDALDataset> tmpSrcDS;
2736 60 : if (m_tilingScheme == "raster" && !bHasNorthUpSrcGT)
2737 : {
2738 1 : CPLStringList aosOptions;
2739 1 : aosOptions.AddString("-of");
2740 1 : aosOptions.AddString("VRT");
2741 1 : aosOptions.AddString("-a_ullr");
2742 1 : aosOptions.AddString(CPLSPrintf("%.17g", adfSrcGeoTransformModif[0]));
2743 1 : aosOptions.AddString(CPLSPrintf("%.17g", adfSrcGeoTransformModif[3]));
2744 : aosOptions.AddString(
2745 1 : CPLSPrintf("%.17g", adfSrcGeoTransformModif[0] +
2746 1 : nSrcWidth * adfSrcGeoTransformModif[1]));
2747 : aosOptions.AddString(
2748 1 : CPLSPrintf("%.17g", adfSrcGeoTransformModif[3] +
2749 1 : nSrcHeight * adfSrcGeoTransformModif[5]));
2750 1 : if (oSRS_TMS.IsEmpty())
2751 : {
2752 1 : aosOptions.AddString("-a_srs");
2753 1 : aosOptions.AddString("none");
2754 : }
2755 :
2756 : GDALTranslateOptions *psOptions =
2757 1 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
2758 :
2759 1 : tmpSrcDS.reset(GDALDataset::FromHandle(GDALTranslate(
2760 : "", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
2761 1 : GDALTranslateOptionsFree(psOptions);
2762 1 : if (!tmpSrcDS)
2763 0 : return false;
2764 : }
2765 61 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
2766 61 : tmpSrcDS ? tmpSrcDS.get() : poSrcDS, &oFakeMaxZoomDS, aosTO.List()));
2767 60 : CPLAssert(hTransformArg);
2768 :
2769 : /* -------------------------------------------------------------------- */
2770 : /* Warp the transformer with a linear approximator */
2771 : /* -------------------------------------------------------------------- */
2772 60 : hTransformArg.reset(GDALCreateApproxTransformer(
2773 : GDALGenImgProjTransform, hTransformArg.release(), 0.125));
2774 60 : GDALApproxTransformerOwnsSubtransformer(hTransformArg.get(), TRUE);
2775 :
2776 60 : psWO->pfnTransformer = GDALApproxTransform;
2777 60 : psWO->pTransformerArg = hTransformArg.get();
2778 :
2779 : /* -------------------------------------------------------------------- */
2780 : /* Determine total number of tiles */
2781 : /* -------------------------------------------------------------------- */
2782 60 : uint64_t nTotalTiles = static_cast<uint64_t>(nMaxTileY - nMinTileY + 1) *
2783 60 : (nMaxTileX - nMinTileX + 1);
2784 60 : const uint64_t nBaseTiles = nTotalTiles;
2785 60 : std::atomic<uint64_t> nCurTile = 0;
2786 60 : bool bRet = true;
2787 :
2788 95 : for (int iZ = m_maxZoomLevel - 1;
2789 95 : bRet && bIntersects && iZ >= m_minZoomLevel; --iZ)
2790 : {
2791 70 : auto ovrTileMatrix = tileMatrixList[iZ];
2792 35 : int nOvrMinTileX = 0;
2793 35 : int nOvrMinTileY = 0;
2794 35 : int nOvrMaxTileX = 0;
2795 35 : int nOvrMaxTileY = 0;
2796 : bRet =
2797 70 : GetTileIndices(ovrTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
2798 : nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX,
2799 35 : nOvrMaxTileY, m_noIntersectionIsOK, bIntersects);
2800 35 : if (bIntersects)
2801 : {
2802 35 : nTotalTiles +=
2803 35 : static_cast<uint64_t>(nOvrMaxTileY - nOvrMinTileY + 1) *
2804 35 : (nOvrMaxTileX - nOvrMinTileX + 1);
2805 : }
2806 : }
2807 :
2808 : /* -------------------------------------------------------------------- */
2809 : /* Generate tiles at max zoom level */
2810 : /* -------------------------------------------------------------------- */
2811 120 : GDALWarpOperation oWO;
2812 :
2813 60 : bRet = oWO.Initialize(psWO.get()) == CE_None && bRet;
2814 :
2815 : const auto GetUpdatedCreationOptions =
2816 260 : [this](const gdal::TileMatrixSet::TileMatrix &oTM)
2817 : {
2818 94 : CPLStringList aosCreationOptions(m_creationOptions);
2819 94 : if (m_outputFormat == "GTiff")
2820 : {
2821 48 : if (aosCreationOptions.FetchNameValue("TILED") == nullptr &&
2822 24 : aosCreationOptions.FetchNameValue("BLOCKYSIZE") == nullptr)
2823 : {
2824 24 : if (oTM.mTileWidth <= 512 && oTM.mTileHeight <= 512)
2825 : {
2826 22 : aosCreationOptions.SetNameValue(
2827 22 : "BLOCKYSIZE", CPLSPrintf("%d", oTM.mTileHeight));
2828 : }
2829 : else
2830 : {
2831 2 : aosCreationOptions.SetNameValue("TILED", "YES");
2832 : }
2833 : }
2834 24 : if (aosCreationOptions.FetchNameValue("COMPRESS") == nullptr)
2835 24 : aosCreationOptions.SetNameValue("COMPRESS", "LZW");
2836 : }
2837 70 : else if (m_outputFormat == "COG")
2838 : {
2839 2 : if (aosCreationOptions.FetchNameValue("OVERVIEW_RESAMPLING") ==
2840 : nullptr)
2841 : {
2842 : aosCreationOptions.SetNameValue("OVERVIEW_RESAMPLING",
2843 2 : m_overviewResampling.c_str());
2844 : }
2845 2 : if (aosCreationOptions.FetchNameValue("BLOCKSIZE") == nullptr &&
2846 2 : oTM.mTileWidth <= 512 && oTM.mTileWidth == oTM.mTileHeight)
2847 : {
2848 : aosCreationOptions.SetNameValue(
2849 2 : "BLOCKSIZE", CPLSPrintf("%d", oTM.mTileWidth));
2850 : }
2851 : }
2852 94 : return aosCreationOptions;
2853 60 : };
2854 :
2855 60 : VSIMkdir(m_outputDirectory.c_str(), 0755);
2856 : VSIStatBufL sStat;
2857 119 : if (VSIStatL(m_outputDirectory.c_str(), &sStat) != 0 ||
2858 59 : !VSI_ISDIR(sStat.st_mode))
2859 : {
2860 1 : ReportError(CE_Failure, CPLE_FileIO,
2861 : "Cannot create output directory %s",
2862 : m_outputDirectory.c_str());
2863 1 : return false;
2864 : }
2865 :
2866 118 : OGRSpatialReference oWGS84;
2867 59 : oWGS84.importFromEPSG(4326);
2868 59 : oWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2869 :
2870 59 : std::unique_ptr<OGRCoordinateTransformation> poCTToWGS84;
2871 59 : if (!oSRS_TMS.IsEmpty())
2872 : {
2873 116 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
2874 58 : poCTToWGS84.reset(
2875 : OGRCreateCoordinateTransformation(&oSRS_TMS, &oWGS84));
2876 : }
2877 :
2878 61 : const bool kmlCompatible = m_kml &&
2879 17 : [this, &poTMS, &poCTToWGS84, bInvertAxisTMS]()
2880 : {
2881 2 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
2882 2 : double dfX = poTMS->tileMatrixList()[0].mTopLeftX;
2883 2 : double dfY = poTMS->tileMatrixList()[0].mTopLeftY;
2884 2 : if (bInvertAxisTMS)
2885 0 : std::swap(dfX, dfY);
2886 3 : return (m_minZoomLevel == m_maxZoomLevel ||
2887 1 : (poTMS->haveAllLevelsSameTopLeft() &&
2888 1 : poTMS->haveAllLevelsSameTileSize() &&
2889 3 : poTMS->hasOnlyPowerOfTwoVaryingScales())) &&
2890 7 : poCTToWGS84 && poCTToWGS84->Transform(1, &dfX, &dfY);
2891 2 : }();
2892 : const int kmlTileSize =
2893 59 : m_tileSize > 0 ? m_tileSize : poTMS->tileMatrixList()[0].mTileWidth;
2894 59 : if (m_kml && !kmlCompatible)
2895 : {
2896 0 : ReportError(CE_Failure, CPLE_NotSupported,
2897 : "Tiling scheme not compatible with KML output");
2898 0 : return false;
2899 : }
2900 :
2901 59 : if (m_title.empty())
2902 57 : m_title = CPLGetFilename(m_dataset.GetName().c_str());
2903 :
2904 59 : if (!m_url.empty())
2905 : {
2906 2 : if (m_url.back() != '/')
2907 2 : m_url += '/';
2908 4 : std::string out_path = m_outputDirectory;
2909 2 : if (m_outputDirectory.back() == '/')
2910 0 : out_path.pop_back();
2911 2 : m_url += CPLGetFilename(out_path.c_str());
2912 : }
2913 :
2914 59 : CPLWorkerThreadPool oThreadPool;
2915 :
2916 : {
2917 : PerThreadMaxZoomResourceManager oResourceManager(
2918 59 : poSrcDS, psWO.get(), hTransformArg.get(), oFakeMaxZoomDS,
2919 177 : dstBuffer.size());
2920 :
2921 : const CPLStringList aosCreationOptions(
2922 118 : GetUpdatedCreationOptions(tileMatrix));
2923 :
2924 59 : CPLDebug("gdal_raster_tile",
2925 : "Generating tiles z=%d, y=%d...%d, x=%d...%d", m_maxZoomLevel,
2926 : nMinTileY, nMaxTileY, nMinTileX, nMaxTileX);
2927 :
2928 59 : if (static_cast<uint64_t>(m_numThreads) > nBaseTiles)
2929 37 : m_numThreads = static_cast<int>(nBaseTiles);
2930 :
2931 59 : if (bRet && m_numThreads > 1)
2932 : {
2933 28 : CPLDebug("gdal_raster_tile", "Using %d threads", m_numThreads);
2934 28 : bRet = oThreadPool.Setup(m_numThreads, nullptr, nullptr);
2935 : }
2936 :
2937 59 : std::atomic<bool> bFailure = false;
2938 59 : std::atomic<int> nQueuedJobs = 0;
2939 :
2940 156 : for (int iY = nMinTileY; bRet && iY <= nMaxTileY; ++iY)
2941 : {
2942 360 : for (int iX = nMinTileX; bRet && iX <= nMaxTileX; ++iX)
2943 : {
2944 263 : if (m_numThreads > 1)
2945 : {
2946 232 : auto job = [this, &oResourceManager, &bFailure, &nCurTile,
2947 : &nQueuedJobs, poDstDriver, pszExtension,
2948 : &aosCreationOptions, &psWO, &tileMatrix,
2949 : nDstBands, iX, iY, nMinTileX, nMinTileY,
2950 1870 : poColorTable, bUserAskedForAlpha]()
2951 : {
2952 464 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
2953 :
2954 231 : --nQueuedJobs;
2955 463 : auto resources = oResourceManager.AcquireResources();
2956 464 : if (resources &&
2957 696 : GenerateTile(
2958 232 : resources->poSrcDS.get(), poDstDriver,
2959 : pszExtension, aosCreationOptions.List(),
2960 232 : *(resources->poWO.get()),
2961 232 : *(resources->poFakeMaxZoomDS->GetSpatialRef()),
2962 232 : psWO->eWorkingDataType, tileMatrix,
2963 232 : m_outputDirectory, nDstBands,
2964 232 : psWO->padfDstNoDataReal
2965 0 : ? &(psWO->padfDstNoDataReal[0])
2966 : : nullptr,
2967 232 : m_maxZoomLevel, iX, iY, m_convention, nMinTileX,
2968 232 : nMinTileY, m_skipBlank, bUserAskedForAlpha,
2969 232 : m_auxXML, m_resume, m_metadata, poColorTable,
2970 464 : resources->dstBuffer))
2971 : {
2972 216 : oResourceManager.ReleaseResources(
2973 216 : std::move(resources));
2974 : }
2975 : else
2976 : {
2977 16 : oResourceManager.SetError();
2978 16 : bFailure = true;
2979 : }
2980 232 : ++nCurTile;
2981 464 : };
2982 :
2983 : // Avoid queueing too many jobs at once
2984 253 : while (bRet && nQueuedJobs / 10 > m_numThreads)
2985 : {
2986 21 : oThreadPool.WaitEvent();
2987 :
2988 21 : bRet &=
2989 32 : !bFailure &&
2990 11 : (!pfnProgress ||
2991 22 : pfnProgress(static_cast<double>(nCurTile) /
2992 11 : static_cast<double>(nTotalTiles),
2993 21 : "", pProgressData));
2994 : }
2995 :
2996 232 : ++nQueuedJobs;
2997 232 : oThreadPool.SubmitJob(std::move(job));
2998 : }
2999 : else
3000 : {
3001 62 : bRet = GenerateTile(
3002 : poSrcDS, poDstDriver, pszExtension,
3003 : aosCreationOptions.List(), oWO, oSRS_TMS,
3004 31 : psWO->eWorkingDataType, tileMatrix, m_outputDirectory,
3005 : nDstBands,
3006 31 : psWO->padfDstNoDataReal ? &(psWO->padfDstNoDataReal[0])
3007 : : nullptr,
3008 31 : m_maxZoomLevel, iX, iY, m_convention, nMinTileX,
3009 31 : nMinTileY, m_skipBlank, bUserAskedForAlpha, m_auxXML,
3010 31 : m_resume, m_metadata, poColorTable, dstBuffer);
3011 :
3012 31 : ++nCurTile;
3013 33 : bRet &= (!pfnProgress ||
3014 4 : pfnProgress(static_cast<double>(nCurTile) /
3015 2 : static_cast<double>(nTotalTiles),
3016 31 : "", pProgressData));
3017 : }
3018 : }
3019 : }
3020 :
3021 59 : if (m_numThreads > 1)
3022 : {
3023 : // Wait for completion of all jobs
3024 154 : while (bRet && nQueuedJobs > 0)
3025 : {
3026 126 : oThreadPool.WaitEvent();
3027 175 : bRet &= !bFailure &&
3028 49 : (!pfnProgress ||
3029 98 : pfnProgress(static_cast<double>(nCurTile) /
3030 49 : static_cast<double>(nTotalTiles),
3031 126 : "", pProgressData));
3032 : }
3033 28 : oThreadPool.WaitCompletion();
3034 28 : bRet &=
3035 29 : !bFailure && (!pfnProgress ||
3036 2 : pfnProgress(static_cast<double>(nCurTile) /
3037 1 : static_cast<double>(nTotalTiles),
3038 28 : "", pProgressData));
3039 :
3040 28 : if (!oResourceManager.GetErrorMsg().empty())
3041 : {
3042 : // Re-emit error message from worker thread to main thread
3043 1 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
3044 1 : oResourceManager.GetErrorMsg().c_str());
3045 : }
3046 : }
3047 :
3048 59 : if (m_kml && bRet)
3049 : {
3050 4 : for (int iY = nMinTileY; iY <= nMaxTileY; ++iY)
3051 : {
3052 4 : for (int iX = nMinTileX; iX <= nMaxTileX; ++iX)
3053 : {
3054 : const int nFileY =
3055 2 : GetFileY(iY, poTMS->tileMatrixList()[m_maxZoomLevel],
3056 2 : m_convention);
3057 : std::string osFilename = CPLFormFilenameSafe(
3058 : m_outputDirectory.c_str(),
3059 4 : CPLSPrintf("%d", m_maxZoomLevel), nullptr);
3060 4 : osFilename = CPLFormFilenameSafe(
3061 2 : osFilename.c_str(), CPLSPrintf("%d", iX), nullptr);
3062 4 : osFilename = CPLFormFilenameSafe(
3063 : osFilename.c_str(),
3064 2 : CPLSPrintf("%d.%s", nFileY, pszExtension), nullptr);
3065 2 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
3066 : {
3067 4 : GenerateKML(m_outputDirectory, m_title, iX, iY,
3068 : m_maxZoomLevel, kmlTileSize, pszExtension,
3069 2 : m_url, poTMS.get(), bInvertAxisTMS,
3070 2 : m_convention, poCTToWGS84.get(), {});
3071 : }
3072 : }
3073 : }
3074 : }
3075 : }
3076 :
3077 : /* -------------------------------------------------------------------- */
3078 : /* Generate tiles at lower zoom levels */
3079 : /* -------------------------------------------------------------------- */
3080 94 : for (int iZ = m_maxZoomLevel - 1; bRet && iZ >= m_minZoomLevel; --iZ)
3081 : {
3082 70 : auto srcTileMatrix = tileMatrixList[iZ + 1];
3083 35 : int nSrcMinTileX = 0;
3084 35 : int nSrcMinTileY = 0;
3085 35 : int nSrcMaxTileX = 0;
3086 35 : int nSrcMaxTileY = 0;
3087 :
3088 35 : CPL_IGNORE_RET_VAL(
3089 35 : GetTileIndices(srcTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
3090 : nSrcMinTileX, nSrcMinTileY, nSrcMaxTileX,
3091 35 : nSrcMaxTileY, m_noIntersectionIsOK, bIntersects));
3092 :
3093 : MosaicDataset oSrcDS(
3094 70 : CPLFormFilenameSafe(m_outputDirectory.c_str(),
3095 : CPLSPrintf("%d", iZ + 1), nullptr),
3096 35 : pszExtension, m_outputFormat, poSrcDS, srcTileMatrix, oSRS_TMS,
3097 : nSrcMinTileX, nSrcMinTileY, nSrcMaxTileX, nSrcMaxTileY,
3098 35 : m_convention, nDstBands, psWO->eWorkingDataType,
3099 9 : psWO->padfDstNoDataReal ? &(psWO->padfDstNoDataReal[0]) : nullptr,
3100 184 : m_metadata, poColorTable);
3101 :
3102 70 : auto ovrTileMatrix = tileMatrixList[iZ];
3103 35 : int nOvrMinTileX = 0;
3104 35 : int nOvrMinTileY = 0;
3105 35 : int nOvrMaxTileX = 0;
3106 35 : int nOvrMaxTileY = 0;
3107 35 : CPL_IGNORE_RET_VAL(
3108 35 : GetTileIndices(ovrTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
3109 : nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX,
3110 35 : nOvrMaxTileY, m_noIntersectionIsOK, bIntersects));
3111 35 : bRet = bIntersects;
3112 :
3113 35 : if (bRet)
3114 : {
3115 35 : CPLDebug("gdal_raster_tile",
3116 : "Generating overview tiles z=%d, y=%d...%d, x=%d...%d", iZ,
3117 : nOvrMinTileY, nOvrMaxTileY, nOvrMinTileX, nOvrMaxTileX);
3118 : }
3119 :
3120 : const CPLStringList aosCreationOptions(
3121 70 : GetUpdatedCreationOptions(ovrTileMatrix));
3122 :
3123 70 : PerThreadLowerZoomResourceManager oResourceManager(oSrcDS);
3124 35 : std::atomic<bool> bFailure = false;
3125 35 : std::atomic<int> nQueuedJobs = 0;
3126 :
3127 35 : const bool bUseThreads =
3128 52 : m_numThreads > 1 &&
3129 17 : (nOvrMaxTileY > nOvrMinTileY || nOvrMaxTileX > nOvrMinTileX);
3130 :
3131 78 : for (int iY = nOvrMinTileY; bRet && iY <= nOvrMaxTileY; ++iY)
3132 : {
3133 114 : for (int iX = nOvrMinTileX; bRet && iX <= nOvrMaxTileX; ++iX)
3134 : {
3135 71 : if (bUseThreads)
3136 : {
3137 40 : auto job = [this, &oResourceManager, poDstDriver, &bFailure,
3138 : &nCurTile, &nQueuedJobs, pszExtension,
3139 : &aosCreationOptions, &aosWarpOptions,
3140 : &ovrTileMatrix, iZ, iX, iY,
3141 280 : bUserAskedForAlpha]()
3142 : {
3143 80 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3144 :
3145 40 : --nQueuedJobs;
3146 80 : auto resources = oResourceManager.AcquireResources();
3147 80 : if (resources &&
3148 80 : GenerateOverviewTile(
3149 40 : *(resources->poSrcDS.get()), poDstDriver,
3150 40 : m_outputFormat, pszExtension,
3151 : aosCreationOptions.List(),
3152 40 : aosWarpOptions.List(), m_overviewResampling,
3153 40 : ovrTileMatrix, m_outputDirectory, iZ, iX, iY,
3154 40 : m_convention, m_skipBlank, bUserAskedForAlpha,
3155 80 : m_auxXML, m_resume))
3156 : {
3157 40 : oResourceManager.ReleaseResources(
3158 40 : std::move(resources));
3159 : }
3160 : else
3161 : {
3162 0 : oResourceManager.SetError();
3163 0 : bFailure = true;
3164 : }
3165 40 : ++nCurTile;
3166 80 : };
3167 :
3168 : // Avoid queueing too many jobs at once
3169 40 : while (bRet && nQueuedJobs / 10 > m_numThreads)
3170 : {
3171 0 : oThreadPool.WaitEvent();
3172 :
3173 0 : bRet &=
3174 0 : !bFailure &&
3175 0 : (!pfnProgress ||
3176 0 : pfnProgress(static_cast<double>(nCurTile) /
3177 0 : static_cast<double>(nTotalTiles),
3178 0 : "", pProgressData));
3179 : }
3180 :
3181 40 : ++nQueuedJobs;
3182 40 : oThreadPool.SubmitJob(std::move(job));
3183 : }
3184 : else
3185 : {
3186 31 : bRet = GenerateOverviewTile(
3187 31 : oSrcDS, poDstDriver, m_outputFormat, pszExtension,
3188 31 : aosCreationOptions.List(), aosWarpOptions.List(),
3189 31 : m_overviewResampling, ovrTileMatrix, m_outputDirectory,
3190 31 : iZ, iX, iY, m_convention, m_skipBlank,
3191 31 : bUserAskedForAlpha, m_auxXML, m_resume);
3192 :
3193 31 : ++nCurTile;
3194 32 : bRet &= (!pfnProgress ||
3195 2 : pfnProgress(static_cast<double>(nCurTile) /
3196 1 : static_cast<double>(nTotalTiles),
3197 31 : "", pProgressData));
3198 : }
3199 : }
3200 : }
3201 :
3202 35 : if (bUseThreads)
3203 : {
3204 : // Wait for completion of all jobs
3205 30 : while (bRet && nQueuedJobs > 0)
3206 : {
3207 26 : oThreadPool.WaitEvent();
3208 39 : bRet &= !bFailure &&
3209 13 : (!pfnProgress ||
3210 26 : pfnProgress(static_cast<double>(nCurTile) /
3211 13 : static_cast<double>(nTotalTiles),
3212 26 : "", pProgressData));
3213 : }
3214 4 : oThreadPool.WaitCompletion();
3215 4 : bRet &=
3216 6 : !bFailure && (!pfnProgress ||
3217 4 : pfnProgress(static_cast<double>(nCurTile) /
3218 2 : static_cast<double>(nTotalTiles),
3219 4 : "", pProgressData));
3220 :
3221 4 : if (!oResourceManager.GetErrorMsg().empty())
3222 : {
3223 : // Re-emit error message from worker thread to main thread
3224 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
3225 0 : oResourceManager.GetErrorMsg().c_str());
3226 : }
3227 : }
3228 :
3229 35 : if (m_kml && bRet)
3230 : {
3231 2 : for (int iY = nOvrMinTileY; bRet && iY <= nOvrMaxTileY; ++iY)
3232 : {
3233 2 : for (int iX = nOvrMinTileX; bRet && iX <= nOvrMaxTileX; ++iX)
3234 : {
3235 : int nFileY =
3236 1 : GetFileY(iY, poTMS->tileMatrixList()[iZ], m_convention);
3237 : std::string osFilename =
3238 : CPLFormFilenameSafe(m_outputDirectory.c_str(),
3239 2 : CPLSPrintf("%d", iZ), nullptr);
3240 2 : osFilename = CPLFormFilenameSafe(
3241 1 : osFilename.c_str(), CPLSPrintf("%d", iX), nullptr);
3242 2 : osFilename = CPLFormFilenameSafe(
3243 : osFilename.c_str(),
3244 1 : CPLSPrintf("%d.%s", nFileY, pszExtension), nullptr);
3245 1 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
3246 : {
3247 1 : std::vector<TileCoordinates> children;
3248 :
3249 3 : for (int iChildY = 0; iChildY <= 1; ++iChildY)
3250 : {
3251 6 : for (int iChildX = 0; iChildX <= 1; ++iChildX)
3252 : {
3253 : nFileY =
3254 4 : GetFileY(iY * 2 + iChildY,
3255 4 : poTMS->tileMatrixList()[iZ + 1],
3256 4 : m_convention);
3257 8 : osFilename = CPLFormFilenameSafe(
3258 : m_outputDirectory.c_str(),
3259 4 : CPLSPrintf("%d", iZ + 1), nullptr);
3260 8 : osFilename = CPLFormFilenameSafe(
3261 : osFilename.c_str(),
3262 4 : CPLSPrintf("%d", iX * 2 + iChildX),
3263 4 : nullptr);
3264 8 : osFilename = CPLFormFilenameSafe(
3265 : osFilename.c_str(),
3266 : CPLSPrintf("%d.%s", nFileY, pszExtension),
3267 4 : nullptr);
3268 4 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
3269 : {
3270 1 : TileCoordinates tc;
3271 1 : tc.nTileX = iX * 2 + iChildX;
3272 1 : tc.nTileY = iY * 2 + iChildY;
3273 1 : tc.nTileZ = iZ + 1;
3274 1 : children.push_back(std::move(tc));
3275 : }
3276 : }
3277 : }
3278 :
3279 2 : GenerateKML(m_outputDirectory, m_title, iX, iY, iZ,
3280 1 : kmlTileSize, pszExtension, m_url,
3281 1 : poTMS.get(), bInvertAxisTMS, m_convention,
3282 : poCTToWGS84.get(), children);
3283 : }
3284 : }
3285 : }
3286 : }
3287 : }
3288 :
3289 450 : const auto IsWebViewerEnabled = [this](const char *name)
3290 : {
3291 150 : return std::find_if(m_webviewers.begin(), m_webviewers.end(),
3292 253 : [name](const std::string &s) {
3293 152 : return s == "all" || s == name;
3294 150 : }) != m_webviewers.end();
3295 59 : };
3296 :
3297 97 : if (bRet && poTMS->identifier() == "GoogleMapsCompatible" &&
3298 38 : IsWebViewerEnabled("leaflet"))
3299 : {
3300 13 : double dfSouthLat = -90;
3301 13 : double dfWestLon = -180;
3302 13 : double dfNorthLat = 90;
3303 13 : double dfEastLon = 180;
3304 :
3305 13 : if (poCTToWGS84)
3306 : {
3307 13 : poCTToWGS84->TransformBounds(
3308 : adfExtent[0], adfExtent[1], adfExtent[2], adfExtent[3],
3309 13 : &dfWestLon, &dfSouthLat, &dfEastLon, &dfNorthLat, 21);
3310 : }
3311 :
3312 13 : GenerateLeaflet(m_outputDirectory, m_title, dfSouthLat, dfWestLon,
3313 : dfNorthLat, dfEastLon, m_minZoomLevel, m_maxZoomLevel,
3314 13 : tileMatrix.mTileWidth, pszExtension, m_url, m_copyright,
3315 13 : m_convention == "xyz");
3316 : }
3317 :
3318 59 : if (bRet && IsWebViewerEnabled("openlayers"))
3319 : {
3320 42 : GenerateOpenLayers(
3321 21 : m_outputDirectory, m_title, adfExtent[0], adfExtent[1],
3322 : adfExtent[2], adfExtent[3], m_minZoomLevel, m_maxZoomLevel,
3323 21 : tileMatrix.mTileWidth, pszExtension, m_url, m_copyright,
3324 21 : *(poTMS.get()), bInvertAxisTMS, oSRS_TMS, m_convention == "xyz");
3325 : }
3326 :
3327 75 : if (bRet && IsWebViewerEnabled("mapml") &&
3328 134 : poTMS->identifier() != "raster" && m_convention == "xyz")
3329 : {
3330 14 : GenerateMapML(m_outputDirectory, m_mapmlTemplate, m_title, nMinTileX,
3331 : nMinTileY, nMaxTileX, nMaxTileY, m_minZoomLevel,
3332 14 : m_maxZoomLevel, pszExtension, m_url, m_copyright,
3333 14 : *(poTMS.get()));
3334 : }
3335 :
3336 59 : if (bRet && m_kml)
3337 : {
3338 4 : std::vector<TileCoordinates> children;
3339 :
3340 2 : auto ovrTileMatrix = tileMatrixList[m_minZoomLevel];
3341 2 : int nOvrMinTileX = 0;
3342 2 : int nOvrMinTileY = 0;
3343 2 : int nOvrMaxTileX = 0;
3344 2 : int nOvrMaxTileY = 0;
3345 2 : CPL_IGNORE_RET_VAL(
3346 2 : GetTileIndices(ovrTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
3347 : nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX,
3348 2 : nOvrMaxTileY, m_noIntersectionIsOK, bIntersects));
3349 :
3350 4 : for (int iY = nOvrMinTileY; bRet && iY <= nOvrMaxTileY; ++iY)
3351 : {
3352 4 : for (int iX = nOvrMinTileX; bRet && iX <= nOvrMaxTileX; ++iX)
3353 : {
3354 2 : int nFileY = GetFileY(
3355 2 : iY, poTMS->tileMatrixList()[m_minZoomLevel], m_convention);
3356 : std::string osFilename = CPLFormFilenameSafe(
3357 : m_outputDirectory.c_str(), CPLSPrintf("%d", m_minZoomLevel),
3358 4 : nullptr);
3359 4 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
3360 2 : CPLSPrintf("%d", iX), nullptr);
3361 4 : osFilename = CPLFormFilenameSafe(
3362 : osFilename.c_str(),
3363 2 : CPLSPrintf("%d.%s", nFileY, pszExtension), nullptr);
3364 2 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
3365 : {
3366 2 : TileCoordinates tc;
3367 2 : tc.nTileX = iX;
3368 2 : tc.nTileY = iY;
3369 2 : tc.nTileZ = m_minZoomLevel;
3370 2 : children.push_back(std::move(tc));
3371 : }
3372 : }
3373 : }
3374 4 : GenerateKML(m_outputDirectory, m_title, -1, -1, -1, kmlTileSize,
3375 2 : pszExtension, m_url, poTMS.get(), bInvertAxisTMS,
3376 2 : m_convention, poCTToWGS84.get(), children);
3377 : }
3378 :
3379 59 : return bRet;
3380 : }
3381 :
3382 : //! @endcond
|
