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_spawn.h"
18 : #include "cpl_worker_thread_pool.h"
19 : #include "gdal_alg_priv.h"
20 : #include "gdal_priv.h"
21 : #include "gdalgetgdalpath.h"
22 : #include "gdalwarper.h"
23 : #include "gdal_utils.h"
24 : #include "ogr_spatialref.h"
25 : #include "memdataset.h"
26 : #include "tilematrixset.hpp"
27 :
28 : #include <algorithm>
29 : #include <array>
30 : #include <atomic>
31 : #include <cinttypes>
32 : #include <cmath>
33 : #include <mutex>
34 :
35 : //! @cond Doxygen_Suppress
36 :
37 : #ifndef _
38 : #define _(x) (x)
39 : #endif
40 :
41 : // Unlikely substring to appear in stdout. We do that in case some GDAL
42 : // driver would output on stdout.
43 : constexpr const char PROGRESS_MARKER[] = {'!', '.', 'x'};
44 :
45 : /************************************************************************/
46 : /* GetThresholdMinTilesPerJob() */
47 : /************************************************************************/
48 :
49 6 : static int GetThresholdMinThreadsForSpawn()
50 : {
51 : // Minimum number of threads for automatic switch to spawning
52 6 : constexpr int THRESHOLD_MIN_THREADS_FOR_SPAWN = 8;
53 :
54 : // Config option for test only
55 6 : return std::max(1, atoi(CPLGetConfigOption(
56 : "GDAL_THRESHOLD_MIN_THREADS_FOR_SPAWN",
57 6 : CPLSPrintf("%d", THRESHOLD_MIN_THREADS_FOR_SPAWN))));
58 : }
59 :
60 : /************************************************************************/
61 : /* GetThresholdMinTilesPerJob() */
62 : /************************************************************************/
63 :
64 129 : static int GetThresholdMinTilesPerJob()
65 : {
66 : // Minimum number of tiles per job to decide for automatic switch to spawning
67 129 : constexpr int THRESHOLD_TILES_PER_JOB = 100;
68 :
69 : // Config option for test only
70 : return std::max(
71 129 : 1, atoi(CPLGetConfigOption("GDAL_THRESHOLD_MIN_TILES_PER_JOB",
72 129 : CPLSPrintf("%d", THRESHOLD_TILES_PER_JOB))));
73 : }
74 :
75 : /************************************************************************/
76 : /* GDALRasterTileAlgorithm::GDALRasterTileAlgorithm() */
77 : /************************************************************************/
78 :
79 118 : GDALRasterTileAlgorithm::GDALRasterTileAlgorithm()
80 118 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
81 : {
82 118 : AddProgressArg();
83 : AddArg("progress-forked", 0, _("Report progress as a forked child"),
84 236 : &m_progressForked)
85 118 : .SetHidden(); // Used in spawn mode
86 236 : AddArg("config-options-in-stdin", 0, _(""), &m_dummy)
87 118 : .SetHidden(); // Used in spawn mode
88 : AddArg("ovr-zoom-level", 0, _("Overview zoom level to compute"),
89 236 : &m_ovrZoomLevel)
90 118 : .SetMinValueIncluded(0)
91 118 : .SetHidden(); // Used in spawn mode
92 236 : AddArg("ovr-min-x", 0, _("Minimum tile X coordinate"), &m_minOvrTileX)
93 118 : .SetMinValueIncluded(0)
94 118 : .SetHidden(); // Used in spawn mode
95 236 : AddArg("ovr-max-x", 0, _("Maximum tile X coordinate"), &m_maxOvrTileX)
96 118 : .SetMinValueIncluded(0)
97 118 : .SetHidden(); // Used in spawn mode
98 236 : AddArg("ovr-min-y", 0, _("Minimum tile Y coordinate"), &m_minOvrTileY)
99 118 : .SetMinValueIncluded(0)
100 118 : .SetHidden(); // Used in spawn mode
101 236 : AddArg("ovr-max-y", 0, _("Maximum tile Y coordinate"), &m_maxOvrTileY)
102 118 : .SetMinValueIncluded(0)
103 118 : .SetHidden(); // Used in spawn mode
104 :
105 118 : AddOpenOptionsArg(&m_openOptions);
106 118 : AddInputFormatsArg(&m_inputFormats)
107 236 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER});
108 118 : AddInputDatasetArg(&m_dataset, GDAL_OF_RASTER);
109 118 : AddOutputFormatArg(&m_outputFormat)
110 118 : .SetDefault(m_outputFormat)
111 : .AddMetadataItem(
112 : GAAMDI_REQUIRED_CAPABILITIES,
113 590 : {GDAL_DCAP_RASTER, GDAL_DCAP_CREATECOPY, GDAL_DMD_EXTENSIONS})
114 236 : .AddMetadataItem(GAAMDI_VRT_COMPATIBLE, {"false"});
115 118 : AddCreationOptionsArg(&m_creationOptions);
116 :
117 236 : AddArg("output", 'o', _("Output directory"), &m_outputDirectory)
118 118 : .SetRequired()
119 118 : .SetMinCharCount(1)
120 118 : .SetPositional();
121 :
122 472 : std::vector<std::string> tilingSchemes{"raster"};
123 1062 : for (const std::string &scheme :
124 2242 : gdal::TileMatrixSet::listPredefinedTileMatrixSets())
125 : {
126 2124 : auto poTMS = gdal::TileMatrixSet::parse(scheme.c_str());
127 2124 : OGRSpatialReference oSRS_TMS;
128 2124 : if (poTMS && !poTMS->hasVariableMatrixWidth() &&
129 1062 : oSRS_TMS.SetFromUserInput(poTMS->crs().c_str()) == OGRERR_NONE)
130 : {
131 1062 : std::string identifier = scheme == "GoogleMapsCompatible"
132 : ? "WebMercatorQuad"
133 2124 : : poTMS->identifier();
134 1062 : m_mapTileMatrixIdentifierToScheme[identifier] = scheme;
135 1062 : tilingSchemes.push_back(std::move(identifier));
136 : }
137 : }
138 236 : AddArg("tiling-scheme", 0, _("Tiling scheme"), &m_tilingScheme)
139 118 : .SetDefault("WebMercatorQuad")
140 118 : .SetChoices(tilingSchemes)
141 : .SetHiddenChoices(
142 : "GoogleMapsCompatible", // equivalent of WebMercatorQuad
143 : "mercator", // gdal2tiles equivalent of WebMercatorQuad
144 : "geodetic" // gdal2tiles (not totally) equivalent of WorldCRS84Quad
145 118 : );
146 :
147 236 : AddArg("min-zoom", 0, _("Minimum zoom level"), &m_minZoomLevel)
148 118 : .SetMinValueIncluded(0);
149 236 : AddArg("max-zoom", 0, _("Maximum zoom level"), &m_maxZoomLevel)
150 118 : .SetMinValueIncluded(0);
151 :
152 236 : AddArg("min-x", 0, _("Minimum tile X coordinate"), &m_minTileX)
153 118 : .SetMinValueIncluded(0);
154 236 : AddArg("max-x", 0, _("Maximum tile X coordinate"), &m_maxTileX)
155 118 : .SetMinValueIncluded(0);
156 236 : AddArg("min-y", 0, _("Minimum tile Y coordinate"), &m_minTileY)
157 118 : .SetMinValueIncluded(0);
158 236 : AddArg("max-y", 0, _("Maximum tile Y coordinate"), &m_maxTileY)
159 118 : .SetMinValueIncluded(0);
160 : AddArg("no-intersection-ok", 0,
161 : _("Whether dataset extent not intersecting tile matrix is only a "
162 : "warning"),
163 118 : &m_noIntersectionIsOK);
164 :
165 : AddArg("resampling", 'r', _("Resampling method for max zoom"),
166 236 : &m_resampling)
167 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
168 : "average", "rms", "mode", "min", "max", "med", "q1", "q3",
169 118 : "sum")
170 118 : .SetDefault("cubic")
171 118 : .SetHiddenChoices("near");
172 : AddArg("overview-resampling", 0, _("Resampling method for overviews"),
173 236 : &m_overviewResampling)
174 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
175 : "average", "rms", "mode", "min", "max", "med", "q1", "q3",
176 118 : "sum")
177 118 : .SetHiddenChoices("near");
178 :
179 : AddArg("convention", 0,
180 : _("Tile numbering convention: xyz (from top) or tms (from bottom)"),
181 236 : &m_convention)
182 118 : .SetDefault(m_convention)
183 118 : .SetChoices("xyz", "tms");
184 236 : AddArg("tile-size", 0, _("Override default tile size"), &m_tileSize)
185 118 : .SetMinValueIncluded(64)
186 118 : .SetMaxValueIncluded(32768);
187 : AddArg("add-alpha", 0, _("Whether to force adding an alpha channel"),
188 236 : &m_addalpha)
189 118 : .SetMutualExclusionGroup("alpha");
190 : AddArg("no-alpha", 0, _("Whether to disable adding an alpha channel"),
191 236 : &m_noalpha)
192 118 : .SetMutualExclusionGroup("alpha");
193 : auto &dstNoDataArg =
194 118 : AddArg("dst-nodata", 0, _("Destination nodata value"), &m_dstNoData);
195 118 : AddArg("skip-blank", 0, _("Do not generate blank tiles"), &m_skipBlank);
196 :
197 : {
198 : auto &arg = AddArg("metadata", 0,
199 236 : _("Add metadata item to output tiles"), &m_metadata)
200 236 : .SetMetaVar("<KEY>=<VALUE>")
201 118 : .SetPackedValuesAllowed(false);
202 17 : arg.AddValidationAction([this, &arg]()
203 135 : { return ParseAndValidateKeyValue(arg); });
204 118 : arg.AddHiddenAlias("mo");
205 : }
206 : AddArg("copy-src-metadata", 0,
207 : _("Whether to copy metadata from source dataset"),
208 118 : &m_copySrcMetadata);
209 :
210 : AddArg("aux-xml", 0, _("Generate .aux.xml sidecar files when needed"),
211 118 : &m_auxXML);
212 118 : AddArg("kml", 0, _("Generate KML files"), &m_kml);
213 118 : AddArg("resume", 0, _("Generate only missing files"), &m_resume);
214 :
215 118 : AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
216 : AddArg("parallel-method", 0, _("Parallelization method (thread / spawn)"),
217 236 : &m_parallelMethod)
218 118 : .SetChoices("thread", "spawn");
219 :
220 118 : constexpr const char *ADVANCED_RESAMPLING_CATEGORY = "Advanced Resampling";
221 : auto &excludedValuesArg =
222 : AddArg("excluded-values", 0,
223 : _("Tuples of values (e.g. <R>,<G>,<B> or (<R1>,<G1>,<B1>),"
224 : "(<R2>,<G2>,<B2>)) that must beignored as contributing source "
225 : "pixels during (average) resampling"),
226 236 : &m_excludedValues)
227 118 : .SetCategory(ADVANCED_RESAMPLING_CATEGORY);
228 : auto &excludedValuesPctThresholdArg =
229 : AddArg(
230 : "excluded-values-pct-threshold", 0,
231 : _("Minimum percentage of source pixels that must be set at one of "
232 : "the --excluded-values to cause the excluded value to be used as "
233 : "the target pixel value"),
234 236 : &m_excludedValuesPctThreshold)
235 118 : .SetDefault(m_excludedValuesPctThreshold)
236 118 : .SetMinValueIncluded(0)
237 118 : .SetMaxValueIncluded(100)
238 118 : .SetCategory(ADVANCED_RESAMPLING_CATEGORY);
239 : auto &nodataValuesPctThresholdArg =
240 : AddArg(
241 : "nodata-values-pct-threshold", 0,
242 : _("Minimum percentage of source pixels that must be set at one of "
243 : "nodata (or alpha=0 or any other way to express transparent pixel"
244 : "to cause the target pixel value to be transparent"),
245 236 : &m_nodataValuesPctThreshold)
246 118 : .SetDefault(m_nodataValuesPctThreshold)
247 118 : .SetMinValueIncluded(0)
248 118 : .SetMaxValueIncluded(100)
249 118 : .SetCategory(ADVANCED_RESAMPLING_CATEGORY);
250 :
251 118 : constexpr const char *PUBLICATION_CATEGORY = "Publication";
252 236 : AddArg("webviewer", 0, _("Web viewer to generate"), &m_webviewers)
253 118 : .SetDefault("all")
254 118 : .SetChoices("none", "all", "leaflet", "openlayers", "mapml")
255 118 : .SetCategory(PUBLICATION_CATEGORY);
256 : AddArg("url", 0,
257 : _("URL address where the generated tiles are going to be published"),
258 236 : &m_url)
259 118 : .SetCategory(PUBLICATION_CATEGORY);
260 236 : AddArg("title", 0, _("Title of the map"), &m_title)
261 118 : .SetCategory(PUBLICATION_CATEGORY);
262 236 : AddArg("copyright", 0, _("Copyright for the map"), &m_copyright)
263 118 : .SetCategory(PUBLICATION_CATEGORY);
264 : AddArg("mapml-template", 0,
265 : _("Filename of a template mapml file where variables will be "
266 : "substituted"),
267 236 : &m_mapmlTemplate)
268 118 : .SetMinCharCount(1)
269 118 : .SetCategory(PUBLICATION_CATEGORY);
270 :
271 118 : AddValidationAction(
272 127 : [this, &dstNoDataArg, &excludedValuesArg,
273 818 : &excludedValuesPctThresholdArg, &nodataValuesPctThresholdArg]()
274 : {
275 127 : if (m_minTileX >= 0 && m_maxTileX >= 0 && m_minTileX > m_maxTileX)
276 : {
277 1 : ReportError(CE_Failure, CPLE_IllegalArg,
278 : "'min-x' must be lesser or equal to 'max-x'");
279 1 : return false;
280 : }
281 :
282 126 : if (m_minTileY >= 0 && m_maxTileY >= 0 && m_minTileY > m_maxTileY)
283 : {
284 1 : ReportError(CE_Failure, CPLE_IllegalArg,
285 : "'min-y' must be lesser or equal to 'max-y'");
286 1 : return false;
287 : }
288 :
289 125 : if (m_minZoomLevel >= 0 && m_maxZoomLevel >= 0 &&
290 33 : m_minZoomLevel > m_maxZoomLevel)
291 : {
292 1 : ReportError(CE_Failure, CPLE_IllegalArg,
293 : "'min-zoom' must be lesser or equal to 'max-zoom'");
294 1 : return false;
295 : }
296 :
297 124 : if (m_addalpha && dstNoDataArg.IsExplicitlySet())
298 : {
299 1 : ReportError(
300 : CE_Failure, CPLE_IllegalArg,
301 : "'add-alpha' and 'dst-nodata' are mutually exclusive");
302 1 : return false;
303 : }
304 :
305 363 : for (const auto *arg :
306 : {&excludedValuesArg, &excludedValuesPctThresholdArg,
307 486 : &nodataValuesPctThresholdArg})
308 : {
309 365 : if (arg->IsExplicitlySet() && m_resampling != "average")
310 : {
311 1 : ReportError(CE_Failure, CPLE_AppDefined,
312 : "'%s' can only be specified if 'resampling' is "
313 : "set to 'average'",
314 1 : arg->GetName().c_str());
315 2 : return false;
316 : }
317 365 : if (arg->IsExplicitlySet() && !m_overviewResampling.empty() &&
318 1 : m_overviewResampling != "average")
319 : {
320 1 : ReportError(CE_Failure, CPLE_AppDefined,
321 : "'%s' can only be specified if "
322 : "'overview-resampling' is set to 'average'",
323 1 : arg->GetName().c_str());
324 1 : return false;
325 : }
326 : }
327 :
328 121 : return true;
329 : });
330 118 : }
331 :
332 : /************************************************************************/
333 : /* GetTileIndices() */
334 : /************************************************************************/
335 :
336 250 : static bool GetTileIndices(gdal::TileMatrixSet::TileMatrix &tileMatrix,
337 : bool bInvertAxisTMS, int tileSize,
338 : const double adfExtent[4], int &nMinTileX,
339 : int &nMinTileY, int &nMaxTileX, int &nMaxTileY,
340 : bool noIntersectionIsOK, bool &bIntersects,
341 : bool checkRasterOverflow = true)
342 : {
343 250 : if (tileSize > 0)
344 : {
345 18 : tileMatrix.mResX *=
346 18 : static_cast<double>(tileMatrix.mTileWidth) / tileSize;
347 18 : tileMatrix.mResY *=
348 18 : static_cast<double>(tileMatrix.mTileHeight) / tileSize;
349 18 : tileMatrix.mTileWidth = tileSize;
350 18 : tileMatrix.mTileHeight = tileSize;
351 : }
352 :
353 250 : if (bInvertAxisTMS)
354 2 : std::swap(tileMatrix.mTopLeftX, tileMatrix.mTopLeftY);
355 :
356 250 : const double dfTileWidth = tileMatrix.mResX * tileMatrix.mTileWidth;
357 250 : const double dfTileHeight = tileMatrix.mResY * tileMatrix.mTileHeight;
358 :
359 250 : constexpr double EPSILON = 1e-3;
360 250 : const double dfMinTileX =
361 250 : (adfExtent[0] - tileMatrix.mTopLeftX) / dfTileWidth;
362 250 : nMinTileX = static_cast<int>(
363 500 : std::clamp(std::floor(dfMinTileX + EPSILON), 0.0,
364 250 : static_cast<double>(tileMatrix.mMatrixWidth - 1)));
365 250 : const double dfMinTileY =
366 250 : (tileMatrix.mTopLeftY - adfExtent[3]) / dfTileHeight;
367 250 : nMinTileY = static_cast<int>(
368 500 : std::clamp(std::floor(dfMinTileY + EPSILON), 0.0,
369 250 : static_cast<double>(tileMatrix.mMatrixHeight - 1)));
370 250 : const double dfMaxTileX =
371 250 : (adfExtent[2] - tileMatrix.mTopLeftX) / dfTileWidth;
372 250 : nMaxTileX = static_cast<int>(
373 500 : std::clamp(std::floor(dfMaxTileX + EPSILON), 0.0,
374 250 : static_cast<double>(tileMatrix.mMatrixWidth - 1)));
375 250 : const double dfMaxTileY =
376 250 : (tileMatrix.mTopLeftY - adfExtent[1]) / dfTileHeight;
377 250 : nMaxTileY = static_cast<int>(
378 500 : std::clamp(std::floor(dfMaxTileY + EPSILON), 0.0,
379 250 : static_cast<double>(tileMatrix.mMatrixHeight - 1)));
380 :
381 250 : bIntersects = (dfMinTileX <= tileMatrix.mMatrixWidth && dfMaxTileX >= 0 &&
382 500 : dfMinTileY <= tileMatrix.mMatrixHeight && dfMaxTileY >= 0);
383 250 : if (!bIntersects)
384 : {
385 2 : CPLDebug("gdal_raster_tile",
386 : "dfMinTileX=%g dfMinTileY=%g dfMaxTileX=%g dfMaxTileY=%g",
387 : dfMinTileX, dfMinTileY, dfMaxTileX, dfMaxTileY);
388 2 : CPLError(noIntersectionIsOK ? CE_Warning : CE_Failure, CPLE_AppDefined,
389 : "Extent of source dataset is not compatible with extent of "
390 : "tile matrix %s",
391 : tileMatrix.mId.c_str());
392 2 : return noIntersectionIsOK;
393 : }
394 248 : if (checkRasterOverflow)
395 : {
396 163 : if (nMaxTileX - nMinTileX + 1 > INT_MAX / tileMatrix.mTileWidth ||
397 163 : nMaxTileY - nMinTileY + 1 > INT_MAX / tileMatrix.mTileHeight)
398 : {
399 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too large zoom level");
400 0 : return false;
401 : }
402 : }
403 248 : return true;
404 : }
405 :
406 : /************************************************************************/
407 : /* GetFileY() */
408 : /************************************************************************/
409 :
410 2818 : static int GetFileY(int iY, const gdal::TileMatrixSet::TileMatrix &tileMatrix,
411 : const std::string &convention)
412 : {
413 2818 : return convention == "xyz" ? iY : tileMatrix.mMatrixHeight - 1 - iY;
414 : }
415 :
416 : /************************************************************************/
417 : /* GenerateTile() */
418 : /************************************************************************/
419 :
420 378 : static bool GenerateTile(
421 : GDALDataset *poSrcDS, GDALDriver *m_poDstDriver, const char *pszExtension,
422 : CSLConstList creationOptions, GDALWarpOperation &oWO,
423 : const OGRSpatialReference &oSRS_TMS, GDALDataType eWorkingDataType,
424 : const gdal::TileMatrixSet::TileMatrix &tileMatrix,
425 : const std::string &outputDirectory, int nBands, const double *pdfDstNoData,
426 : int nZoomLevel, int iX, int iY, const std::string &convention,
427 : int nMinTileX, int nMinTileY, bool bSkipBlank, bool bUserAskedForAlpha,
428 : bool bAuxXML, bool bResume, const std::vector<std::string> &metadata,
429 : const GDALColorTable *poColorTable, std::vector<GByte> &dstBuffer)
430 : {
431 : const std::string osDirZ = CPLFormFilenameSafe(
432 756 : outputDirectory.c_str(), CPLSPrintf("%d", nZoomLevel), nullptr);
433 : const std::string osDirX =
434 756 : CPLFormFilenameSafe(osDirZ.c_str(), CPLSPrintf("%d", iX), nullptr);
435 378 : const int iFileY = GetFileY(iY, tileMatrix, convention);
436 : const std::string osFilename = CPLFormFilenameSafe(
437 756 : osDirX.c_str(), CPLSPrintf("%d", iFileY), pszExtension);
438 :
439 378 : if (bResume)
440 : {
441 : VSIStatBufL sStat;
442 5 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
443 5 : return true;
444 : }
445 :
446 373 : const int nDstXOff = (iX - nMinTileX) * tileMatrix.mTileWidth;
447 373 : const int nDstYOff = (iY - nMinTileY) * tileMatrix.mTileHeight;
448 373 : memset(dstBuffer.data(), 0, dstBuffer.size());
449 746 : const CPLErr eErr = oWO.WarpRegionToBuffer(
450 373 : nDstXOff, nDstYOff, tileMatrix.mTileWidth, tileMatrix.mTileHeight,
451 373 : dstBuffer.data(), eWorkingDataType);
452 373 : if (eErr != CE_None)
453 2 : return false;
454 :
455 : const bool bDstHasAlpha =
456 400 : nBands > poSrcDS->GetRasterCount() ||
457 29 : (nBands == poSrcDS->GetRasterCount() &&
458 28 : poSrcDS->GetRasterBand(nBands)->GetColorInterpretation() ==
459 371 : GCI_AlphaBand);
460 371 : const size_t nBytesPerBand = static_cast<size_t>(tileMatrix.mTileWidth) *
461 371 : tileMatrix.mTileHeight *
462 371 : GDALGetDataTypeSizeBytes(eWorkingDataType);
463 371 : if (bDstHasAlpha && bSkipBlank)
464 : {
465 9 : bool bBlank = true;
466 327706 : for (size_t i = 0; i < nBytesPerBand && bBlank; ++i)
467 : {
468 327697 : bBlank = (dstBuffer[(nBands - 1) * nBytesPerBand + i] == 0);
469 : }
470 9 : if (bBlank)
471 5 : return true;
472 : }
473 366 : if (bDstHasAlpha && !bUserAskedForAlpha)
474 : {
475 351 : bool bAllOpaque = true;
476 18527600 : for (size_t i = 0; i < nBytesPerBand && bAllOpaque; ++i)
477 : {
478 18636200 : bAllOpaque = (dstBuffer[(nBands - 1) * nBytesPerBand + i] == 255);
479 : }
480 0 : if (bAllOpaque)
481 276 : nBands--;
482 : }
483 :
484 0 : VSIMkdir(osDirZ.c_str(), 0755);
485 366 : VSIMkdir(osDirX.c_str(), 0755);
486 :
487 : auto memDS = std::unique_ptr<GDALDataset>(
488 366 : MEMDataset::Create("", tileMatrix.mTileWidth, tileMatrix.mTileHeight, 0,
489 732 : eWorkingDataType, nullptr));
490 1475 : for (int i = 0; i < nBands; ++i)
491 : {
492 1109 : char szBuffer[32] = {'\0'};
493 2218 : int nRet = CPLPrintPointer(
494 1109 : szBuffer, dstBuffer.data() + i * nBytesPerBand, sizeof(szBuffer));
495 1109 : szBuffer[nRet] = 0;
496 :
497 1109 : char szOption[64] = {'\0'};
498 1109 : snprintf(szOption, sizeof(szOption), "DATAPOINTER=%s", szBuffer);
499 :
500 1109 : char *apszOptions[] = {szOption, nullptr};
501 :
502 1109 : memDS->AddBand(eWorkingDataType, apszOptions);
503 1109 : auto poDstBand = memDS->GetRasterBand(i + 1);
504 1109 : if (i + 1 <= poSrcDS->GetRasterCount())
505 1034 : poDstBand->SetColorInterpretation(
506 1034 : poSrcDS->GetRasterBand(i + 1)->GetColorInterpretation());
507 : else
508 75 : poDstBand->SetColorInterpretation(GCI_AlphaBand);
509 1109 : if (pdfDstNoData)
510 9 : poDstBand->SetNoDataValue(*pdfDstNoData);
511 1109 : if (i == 0 && poColorTable)
512 1 : poDstBand->SetColorTable(
513 1 : const_cast<GDALColorTable *>(poColorTable));
514 : }
515 732 : const CPLStringList aosMD(metadata);
516 438 : for (const auto [key, value] : cpl::IterateNameValue(aosMD))
517 : {
518 72 : memDS->SetMetadataItem(key, value);
519 : }
520 :
521 366 : GDALGeoTransform gt;
522 732 : gt[0] =
523 366 : tileMatrix.mTopLeftX + iX * tileMatrix.mResX * tileMatrix.mTileWidth;
524 366 : gt[1] = tileMatrix.mResX;
525 366 : gt[2] = 0;
526 732 : gt[3] =
527 366 : tileMatrix.mTopLeftY - iY * tileMatrix.mResY * tileMatrix.mTileHeight;
528 366 : gt[4] = 0;
529 366 : gt[5] = -tileMatrix.mResY;
530 366 : memDS->SetGeoTransform(gt);
531 :
532 366 : memDS->SetSpatialRef(&oSRS_TMS);
533 :
534 : CPLConfigOptionSetter oSetter("GDAL_PAM_ENABLED", bAuxXML ? "YES" : "NO",
535 732 : false);
536 : CPLConfigOptionSetter oSetter2("GDAL_DISABLE_READDIR_ON_OPEN", "YES",
537 732 : false);
538 :
539 366 : std::unique_ptr<CPLConfigOptionSetter> poSetter;
540 : // No need to reopen the dataset at end of CreateCopy() (for PNG
541 : // and JPEG) if we don't need to generate .aux.xml
542 366 : if (!bAuxXML)
543 362 : poSetter = std::make_unique<CPLConfigOptionSetter>(
544 724 : "GDAL_OPEN_AFTER_COPY", "NO", false);
545 366 : CPL_IGNORE_RET_VAL(poSetter);
546 :
547 732 : const std::string osTmpFilename = osFilename + ".tmp." + pszExtension;
548 :
549 : std::unique_ptr<GDALDataset> poOutDS(
550 : m_poDstDriver->CreateCopy(osTmpFilename.c_str(), memDS.get(), false,
551 366 : creationOptions, nullptr, nullptr));
552 366 : bool bRet = poOutDS && poOutDS->Close() == CE_None;
553 366 : poOutDS.reset();
554 366 : if (bRet)
555 : {
556 364 : bRet = VSIRename(osTmpFilename.c_str(), osFilename.c_str()) == 0;
557 364 : if (bAuxXML)
558 : {
559 4 : VSIRename((osTmpFilename + ".aux.xml").c_str(),
560 8 : (osFilename + ".aux.xml").c_str());
561 : }
562 : }
563 : else
564 : {
565 2 : VSIUnlink(osTmpFilename.c_str());
566 : }
567 366 : return bRet;
568 : }
569 :
570 : /************************************************************************/
571 : /* GenerateOverviewTile() */
572 : /************************************************************************/
573 :
574 : static bool
575 113 : GenerateOverviewTile(GDALDataset &oSrcDS, GDALDriver *m_poDstDriver,
576 : const std::string &outputFormat, const char *pszExtension,
577 : CSLConstList creationOptions,
578 : CSLConstList papszWarpOptions,
579 : const std::string &resampling,
580 : const gdal::TileMatrixSet::TileMatrix &tileMatrix,
581 : const std::string &outputDirectory, int nZoomLevel, int iX,
582 : int iY, const std::string &convention, bool bSkipBlank,
583 : bool bUserAskedForAlpha, bool bAuxXML, bool bResume)
584 : {
585 : const std::string osDirZ = CPLFormFilenameSafe(
586 226 : outputDirectory.c_str(), CPLSPrintf("%d", nZoomLevel), nullptr);
587 : const std::string osDirX =
588 226 : CPLFormFilenameSafe(osDirZ.c_str(), CPLSPrintf("%d", iX), nullptr);
589 :
590 113 : const int iFileY = GetFileY(iY, tileMatrix, convention);
591 : const std::string osFilename = CPLFormFilenameSafe(
592 226 : osDirX.c_str(), CPLSPrintf("%d", iFileY), pszExtension);
593 :
594 113 : if (bResume)
595 : {
596 : VSIStatBufL sStat;
597 2 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
598 1 : return true;
599 : }
600 :
601 112 : VSIMkdir(osDirZ.c_str(), 0755);
602 112 : VSIMkdir(osDirX.c_str(), 0755);
603 :
604 224 : CPLStringList aosOptions;
605 :
606 112 : aosOptions.AddString("-of");
607 112 : aosOptions.AddString(outputFormat.c_str());
608 :
609 136 : for (const char *pszCO : cpl::Iterate(creationOptions))
610 : {
611 24 : aosOptions.AddString("-co");
612 24 : aosOptions.AddString(pszCO);
613 : }
614 : CPLConfigOptionSetter oSetter("GDAL_PAM_ENABLED", bAuxXML ? "YES" : "NO",
615 224 : false);
616 : CPLConfigOptionSetter oSetter2("GDAL_DISABLE_READDIR_ON_OPEN", "YES",
617 224 : false);
618 :
619 112 : aosOptions.AddString("-r");
620 112 : aosOptions.AddString(resampling.c_str());
621 :
622 112 : std::unique_ptr<GDALDataset> poOutDS;
623 112 : const double dfMinX =
624 112 : tileMatrix.mTopLeftX + iX * tileMatrix.mResX * tileMatrix.mTileWidth;
625 112 : const double dfMaxY =
626 112 : tileMatrix.mTopLeftY - iY * tileMatrix.mResY * tileMatrix.mTileHeight;
627 112 : const double dfMaxX = dfMinX + tileMatrix.mResX * tileMatrix.mTileWidth;
628 112 : const double dfMinY = dfMaxY - tileMatrix.mResY * tileMatrix.mTileHeight;
629 :
630 : const bool resamplingCompatibleOfTranslate =
631 325 : papszWarpOptions == nullptr &&
632 316 : (resampling == "nearest" || resampling == "average" ||
633 209 : resampling == "bilinear" || resampling == "cubic" ||
634 9 : resampling == "cubicspline" || resampling == "lanczos" ||
635 3 : resampling == "mode");
636 :
637 224 : const std::string osTmpFilename = osFilename + ".tmp." + pszExtension;
638 :
639 112 : if (resamplingCompatibleOfTranslate)
640 : {
641 106 : GDALGeoTransform upperGT;
642 106 : oSrcDS.GetGeoTransform(upperGT);
643 106 : const double dfMinXUpper = upperGT[0];
644 : const double dfMaxXUpper =
645 106 : dfMinXUpper + upperGT[1] * oSrcDS.GetRasterXSize();
646 106 : const double dfMaxYUpper = upperGT[3];
647 : const double dfMinYUpper =
648 106 : dfMaxYUpper + upperGT[5] * oSrcDS.GetRasterYSize();
649 106 : if (dfMinX >= dfMinXUpper && dfMaxX <= dfMaxXUpper &&
650 90 : dfMinY >= dfMinYUpper && dfMaxY <= dfMaxYUpper)
651 : {
652 : // If the overview tile is fully within the extent of the
653 : // upper zoom level, we can use GDALDataset::RasterIO() directly.
654 :
655 90 : const auto eDT = oSrcDS.GetRasterBand(1)->GetRasterDataType();
656 : const size_t nBytesPerBand =
657 90 : static_cast<size_t>(tileMatrix.mTileWidth) *
658 90 : tileMatrix.mTileHeight * GDALGetDataTypeSizeBytes(eDT);
659 : std::vector<GByte> dstBuffer(nBytesPerBand *
660 90 : oSrcDS.GetRasterCount());
661 :
662 90 : const double dfXOff = (dfMinX - dfMinXUpper) / upperGT[1];
663 90 : const double dfYOff = (dfMaxYUpper - dfMaxY) / -upperGT[5];
664 90 : const double dfXSize = (dfMaxX - dfMinX) / upperGT[1];
665 90 : const double dfYSize = (dfMaxY - dfMinY) / -upperGT[5];
666 : GDALRasterIOExtraArg sExtraArg;
667 90 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
668 90 : CPL_IGNORE_RET_VAL(sExtraArg.eResampleAlg);
669 90 : sExtraArg.eResampleAlg =
670 90 : GDALRasterIOGetResampleAlg(resampling.c_str());
671 90 : sExtraArg.dfXOff = dfXOff;
672 90 : sExtraArg.dfYOff = dfYOff;
673 90 : sExtraArg.dfXSize = dfXSize;
674 90 : sExtraArg.dfYSize = dfYSize;
675 90 : sExtraArg.bFloatingPointWindowValidity =
676 90 : sExtraArg.eResampleAlg != GRIORA_NearestNeighbour;
677 90 : constexpr double EPSILON = 1e-3;
678 90 : if (oSrcDS.RasterIO(GF_Read, static_cast<int>(dfXOff + EPSILON),
679 90 : static_cast<int>(dfYOff + EPSILON),
680 90 : static_cast<int>(dfXSize + 0.5),
681 90 : static_cast<int>(dfYSize + 0.5),
682 90 : dstBuffer.data(), tileMatrix.mTileWidth,
683 90 : tileMatrix.mTileHeight, eDT,
684 : oSrcDS.GetRasterCount(), nullptr, 0, 0, 0,
685 90 : &sExtraArg) == CE_None)
686 : {
687 89 : int nDstBands = oSrcDS.GetRasterCount();
688 : const bool bDstHasAlpha =
689 89 : oSrcDS.GetRasterBand(nDstBands)->GetColorInterpretation() ==
690 89 : GCI_AlphaBand;
691 89 : if (bDstHasAlpha && bSkipBlank)
692 : {
693 2 : bool bBlank = true;
694 65545 : for (size_t i = 0; i < nBytesPerBand && bBlank; ++i)
695 : {
696 65543 : bBlank =
697 65543 : (dstBuffer[(nDstBands - 1) * nBytesPerBand + i] ==
698 : 0);
699 : }
700 2 : if (bBlank)
701 1 : return true;
702 1 : bSkipBlank = false;
703 : }
704 88 : if (bDstHasAlpha && !bUserAskedForAlpha)
705 : {
706 88 : bool bAllOpaque = true;
707 5570650 : for (size_t i = 0; i < nBytesPerBand && bAllOpaque; ++i)
708 : {
709 5570560 : bAllOpaque =
710 5570560 : (dstBuffer[(nDstBands - 1) * nBytesPerBand + i] ==
711 : 255);
712 : }
713 88 : if (bAllOpaque)
714 85 : nDstBands--;
715 : }
716 :
717 176 : auto memDS = std::unique_ptr<GDALDataset>(MEMDataset::Create(
718 88 : "", tileMatrix.mTileWidth, tileMatrix.mTileHeight, 0, eDT,
719 176 : nullptr));
720 355 : for (int i = 0; i < nDstBands; ++i)
721 : {
722 267 : char szBuffer[32] = {'\0'};
723 534 : int nRet = CPLPrintPointer(
724 267 : szBuffer, dstBuffer.data() + i * nBytesPerBand,
725 : sizeof(szBuffer));
726 267 : szBuffer[nRet] = 0;
727 :
728 267 : char szOption[64] = {'\0'};
729 267 : snprintf(szOption, sizeof(szOption), "DATAPOINTER=%s",
730 : szBuffer);
731 :
732 267 : char *apszOptions[] = {szOption, nullptr};
733 :
734 267 : memDS->AddBand(eDT, apszOptions);
735 267 : auto poSrcBand = oSrcDS.GetRasterBand(i + 1);
736 267 : auto poDstBand = memDS->GetRasterBand(i + 1);
737 267 : poDstBand->SetColorInterpretation(
738 267 : poSrcBand->GetColorInterpretation());
739 267 : int bHasNoData = false;
740 : const double dfNoData =
741 267 : poSrcBand->GetNoDataValue(&bHasNoData);
742 267 : if (bHasNoData)
743 0 : poDstBand->SetNoDataValue(dfNoData);
744 267 : if (auto poCT = poSrcBand->GetColorTable())
745 0 : poDstBand->SetColorTable(poCT);
746 : }
747 88 : memDS->SetMetadata(oSrcDS.GetMetadata());
748 176 : memDS->SetGeoTransform(GDALGeoTransform(
749 88 : dfMinX, tileMatrix.mResX, 0, dfMaxY, 0, -tileMatrix.mResY));
750 :
751 88 : memDS->SetSpatialRef(oSrcDS.GetSpatialRef());
752 :
753 88 : std::unique_ptr<CPLConfigOptionSetter> poSetter;
754 : // No need to reopen the dataset at end of CreateCopy() (for PNG
755 : // and JPEG) if we don't need to generate .aux.xml
756 88 : if (!bAuxXML)
757 88 : poSetter = std::make_unique<CPLConfigOptionSetter>(
758 176 : "GDAL_OPEN_AFTER_COPY", "NO", false);
759 88 : CPL_IGNORE_RET_VAL(poSetter);
760 :
761 88 : poOutDS.reset(m_poDstDriver->CreateCopy(
762 : osTmpFilename.c_str(), memDS.get(), false, creationOptions,
763 : nullptr, nullptr));
764 89 : }
765 : }
766 : else
767 : {
768 : // If the overview tile is not fully within the extent of the
769 : // upper zoom level, use GDALTranslate() to use VRT padding
770 :
771 16 : aosOptions.AddString("-q");
772 :
773 16 : aosOptions.AddString("-projwin");
774 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
775 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY));
776 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
777 16 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinY));
778 :
779 16 : aosOptions.AddString("-outsize");
780 16 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileWidth));
781 16 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileHeight));
782 :
783 : GDALTranslateOptions *psOptions =
784 16 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
785 16 : poOutDS.reset(GDALDataset::FromHandle(GDALTranslate(
786 : osTmpFilename.c_str(), GDALDataset::ToHandle(&oSrcDS),
787 : psOptions, nullptr)));
788 16 : GDALTranslateOptionsFree(psOptions);
789 : }
790 : }
791 : else
792 : {
793 6 : aosOptions.AddString("-te");
794 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
795 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinY));
796 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
797 6 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY));
798 :
799 6 : aosOptions.AddString("-ts");
800 6 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileWidth));
801 6 : aosOptions.AddString(CPLSPrintf("%d", tileMatrix.mTileHeight));
802 :
803 11 : for (int i = 0; papszWarpOptions && papszWarpOptions[i]; ++i)
804 : {
805 5 : aosOptions.AddString("-wo");
806 5 : aosOptions.AddString(papszWarpOptions[i]);
807 : }
808 :
809 : GDALWarpAppOptions *psOptions =
810 6 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
811 6 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(&oSrcDS);
812 6 : poOutDS.reset(GDALDataset::FromHandle(GDALWarp(
813 : osTmpFilename.c_str(), nullptr, 1, &hSrcDS, psOptions, nullptr)));
814 6 : GDALWarpAppOptionsFree(psOptions);
815 : }
816 :
817 111 : bool bRet = poOutDS != nullptr;
818 111 : if (bRet && bSkipBlank)
819 : {
820 10 : auto poLastBand = poOutDS->GetRasterBand(poOutDS->GetRasterCount());
821 10 : if (poLastBand->GetColorInterpretation() == GCI_AlphaBand)
822 : {
823 : std::vector<GByte> buffer(
824 1 : static_cast<size_t>(tileMatrix.mTileWidth) *
825 1 : tileMatrix.mTileHeight *
826 1 : GDALGetDataTypeSizeBytes(poLastBand->GetRasterDataType()));
827 2 : CPL_IGNORE_RET_VAL(poLastBand->RasterIO(
828 1 : GF_Read, 0, 0, tileMatrix.mTileWidth, tileMatrix.mTileHeight,
829 1 : buffer.data(), tileMatrix.mTileWidth, tileMatrix.mTileHeight,
830 : poLastBand->GetRasterDataType(), 0, 0, nullptr));
831 1 : bool bBlank = true;
832 65537 : for (size_t i = 0; i < buffer.size() && bBlank; ++i)
833 : {
834 65536 : bBlank = (buffer[i] == 0);
835 : }
836 1 : if (bBlank)
837 : {
838 1 : poOutDS.reset();
839 1 : VSIUnlink(osTmpFilename.c_str());
840 1 : if (bAuxXML)
841 0 : VSIUnlink((osTmpFilename + ".aux.xml").c_str());
842 1 : return true;
843 : }
844 : }
845 : }
846 110 : bRet = bRet && poOutDS->Close() == CE_None;
847 110 : poOutDS.reset();
848 110 : if (bRet)
849 : {
850 109 : bRet = VSIRename(osTmpFilename.c_str(), osFilename.c_str()) == 0;
851 109 : if (bAuxXML)
852 : {
853 4 : VSIRename((osTmpFilename + ".aux.xml").c_str(),
854 8 : (osFilename + ".aux.xml").c_str());
855 : }
856 : }
857 : else
858 : {
859 1 : VSIUnlink(osTmpFilename.c_str());
860 : }
861 110 : return bRet;
862 : }
863 :
864 : namespace
865 : {
866 :
867 : /************************************************************************/
868 : /* FakeMaxZoomRasterBand */
869 : /************************************************************************/
870 :
871 : class FakeMaxZoomRasterBand : public GDALRasterBand
872 : {
873 : void *m_pDstBuffer = nullptr;
874 : CPL_DISALLOW_COPY_ASSIGN(FakeMaxZoomRasterBand)
875 :
876 : public:
877 269 : FakeMaxZoomRasterBand(int nBandIn, int nWidth, int nHeight,
878 : int nBlockXSizeIn, int nBlockYSizeIn,
879 : GDALDataType eDT, void *pDstBuffer)
880 269 : : m_pDstBuffer(pDstBuffer)
881 : {
882 269 : nBand = nBandIn;
883 269 : nRasterXSize = nWidth;
884 269 : nRasterYSize = nHeight;
885 269 : nBlockXSize = nBlockXSizeIn;
886 269 : nBlockYSize = nBlockYSizeIn;
887 269 : eDataType = eDT;
888 269 : }
889 :
890 0 : CPLErr IReadBlock(int, int, void *) override
891 : {
892 0 : CPLAssert(false);
893 : return CE_Failure;
894 : }
895 :
896 : #ifdef DEBUG
897 0 : CPLErr IWriteBlock(int, int, void *) override
898 : {
899 0 : CPLAssert(false);
900 : return CE_Failure;
901 : }
902 : #endif
903 :
904 716 : CPLErr IRasterIO(GDALRWFlag eRWFlag, [[maybe_unused]] int nXOff,
905 : [[maybe_unused]] int nYOff, [[maybe_unused]] int nXSize,
906 : [[maybe_unused]] int nYSize, void *pData,
907 : [[maybe_unused]] int nBufXSize,
908 : [[maybe_unused]] int nBufYSize, GDALDataType eBufType,
909 : GSpacing nPixelSpace, [[maybe_unused]] GSpacing nLineSpace,
910 : GDALRasterIOExtraArg *) override
911 : {
912 : // For sake of implementation simplicity, check various assumptions of
913 : // how GDALAlphaMask code does I/O
914 716 : CPLAssert((nXOff % nBlockXSize) == 0);
915 716 : CPLAssert((nYOff % nBlockYSize) == 0);
916 716 : CPLAssert(nXSize == nBufXSize);
917 716 : CPLAssert(nXSize == nBlockXSize);
918 716 : CPLAssert(nYSize == nBufYSize);
919 716 : CPLAssert(nYSize == nBlockYSize);
920 716 : CPLAssert(nLineSpace == nBlockXSize * nPixelSpace);
921 716 : CPLAssert(
922 : nBand ==
923 : poDS->GetRasterCount()); // only alpha band is accessed this way
924 716 : if (eRWFlag == GF_Read)
925 : {
926 358 : double dfZero = 0;
927 358 : GDALCopyWords64(&dfZero, GDT_Float64, 0, pData, eBufType,
928 : static_cast<int>(nPixelSpace),
929 358 : static_cast<size_t>(nBlockXSize) * nBlockYSize);
930 : }
931 : else
932 : {
933 358 : GDALCopyWords64(pData, eBufType, static_cast<int>(nPixelSpace),
934 : m_pDstBuffer, eDataType,
935 : GDALGetDataTypeSizeBytes(eDataType),
936 358 : static_cast<size_t>(nBlockXSize) * nBlockYSize);
937 : }
938 716 : return CE_None;
939 : }
940 : };
941 :
942 : /************************************************************************/
943 : /* FakeMaxZoomDataset */
944 : /************************************************************************/
945 :
946 : // This class is used to create a fake output dataset for GDALWarpOperation.
947 : // In particular we need to implement GDALRasterBand::IRasterIO(GF_Write, ...)
948 : // to catch writes (of one single tile) to the alpha band and redirect them
949 : // to the dstBuffer passed to FakeMaxZoomDataset constructor.
950 :
951 : class FakeMaxZoomDataset : public GDALDataset
952 : {
953 : const int m_nBlockXSize;
954 : const int m_nBlockYSize;
955 : const OGRSpatialReference m_oSRS;
956 : const GDALGeoTransform m_gt{};
957 :
958 : public:
959 85 : FakeMaxZoomDataset(int nWidth, int nHeight, int nBandsIn, int nBlockXSize,
960 : int nBlockYSize, GDALDataType eDT,
961 : const GDALGeoTransform >,
962 : const OGRSpatialReference &oSRS,
963 : std::vector<GByte> &dstBuffer)
964 85 : : m_nBlockXSize(nBlockXSize), m_nBlockYSize(nBlockYSize), m_oSRS(oSRS),
965 85 : m_gt(gt)
966 : {
967 85 : eAccess = GA_Update;
968 85 : nRasterXSize = nWidth;
969 85 : nRasterYSize = nHeight;
970 354 : for (int i = 1; i <= nBandsIn; ++i)
971 : {
972 269 : SetBand(i,
973 : new FakeMaxZoomRasterBand(
974 : i, nWidth, nHeight, nBlockXSize, nBlockYSize, eDT,
975 269 : dstBuffer.data() + static_cast<size_t>(i - 1) *
976 538 : nBlockXSize * nBlockYSize *
977 269 : GDALGetDataTypeSizeBytes(eDT)));
978 : }
979 85 : }
980 :
981 209 : const OGRSpatialReference *GetSpatialRef() const override
982 : {
983 209 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
984 : }
985 :
986 77 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
987 : {
988 77 : gt = m_gt;
989 77 : return CE_None;
990 : }
991 :
992 : using GDALDataset::Clone;
993 :
994 : std::unique_ptr<FakeMaxZoomDataset>
995 8 : Clone(std::vector<GByte> &dstBuffer) const
996 : {
997 : return std::make_unique<FakeMaxZoomDataset>(
998 8 : nRasterXSize, nRasterYSize, nBands, m_nBlockXSize, m_nBlockYSize,
999 16 : GetRasterBand(1)->GetRasterDataType(), m_gt, m_oSRS, dstBuffer);
1000 : }
1001 : };
1002 :
1003 : /************************************************************************/
1004 : /* MosaicRasterBand */
1005 : /************************************************************************/
1006 :
1007 : class MosaicRasterBand : public GDALRasterBand
1008 : {
1009 : const int m_tileMinX;
1010 : const int m_tileMinY;
1011 : const GDALColorInterp m_eColorInterp;
1012 : const gdal::TileMatrixSet::TileMatrix m_oTM;
1013 : const std::string m_convention;
1014 : const std::string m_directory;
1015 : const std::string m_extension;
1016 : const bool m_hasNoData;
1017 : const double m_noData;
1018 : std::unique_ptr<GDALColorTable> m_poColorTable{};
1019 :
1020 : public:
1021 206 : MosaicRasterBand(GDALDataset *poDSIn, int nBandIn, int nWidth, int nHeight,
1022 : int nBlockXSizeIn, int nBlockYSizeIn, GDALDataType eDT,
1023 : GDALColorInterp eColorInterp, int nTileMinX, int nTileMinY,
1024 : const gdal::TileMatrixSet::TileMatrix &oTM,
1025 : const std::string &convention,
1026 : const std::string &directory, const std::string &extension,
1027 : const double *pdfDstNoData,
1028 : const GDALColorTable *poColorTable)
1029 206 : : m_tileMinX(nTileMinX), m_tileMinY(nTileMinY),
1030 : m_eColorInterp(eColorInterp), m_oTM(oTM), m_convention(convention),
1031 : m_directory(directory), m_extension(extension),
1032 206 : m_hasNoData(pdfDstNoData != nullptr),
1033 206 : m_noData(pdfDstNoData ? *pdfDstNoData : 0),
1034 412 : m_poColorTable(poColorTable ? poColorTable->Clone() : nullptr)
1035 : {
1036 206 : poDS = poDSIn;
1037 206 : nBand = nBandIn;
1038 206 : nRasterXSize = nWidth;
1039 206 : nRasterYSize = nHeight;
1040 206 : nBlockXSize = nBlockXSizeIn;
1041 206 : nBlockYSize = nBlockYSizeIn;
1042 206 : eDataType = eDT;
1043 206 : }
1044 :
1045 : CPLErr IReadBlock(int nXBlock, int nYBlock, void *pData) override;
1046 :
1047 2402 : GDALColorTable *GetColorTable() override
1048 : {
1049 2402 : return m_poColorTable.get();
1050 : }
1051 :
1052 529 : GDALColorInterp GetColorInterpretation() override
1053 : {
1054 529 : return m_eColorInterp;
1055 : }
1056 :
1057 1935 : double GetNoDataValue(int *pbHasNoData) override
1058 : {
1059 1935 : if (pbHasNoData)
1060 1926 : *pbHasNoData = m_hasNoData;
1061 1935 : return m_noData;
1062 : }
1063 : };
1064 :
1065 : /************************************************************************/
1066 : /* MosaicDataset */
1067 : /************************************************************************/
1068 :
1069 : // This class is to expose the tiles of a given level as a mosaic that
1070 : // can be used as a source to generate the immediately below zoom level.
1071 :
1072 : class MosaicDataset : public GDALDataset
1073 : {
1074 : friend class MosaicRasterBand;
1075 :
1076 : const std::string m_directory;
1077 : const std::string m_extension;
1078 : const std::string m_format;
1079 : const std::vector<GDALColorInterp> m_aeColorInterp;
1080 : const gdal::TileMatrixSet::TileMatrix &m_oTM;
1081 : const OGRSpatialReference m_oSRS;
1082 : const int m_nTileMinX;
1083 : const int m_nTileMinY;
1084 : const int m_nTileMaxX;
1085 : const int m_nTileMaxY;
1086 : const std::string m_convention;
1087 : const GDALDataType m_eDT;
1088 : const double *const m_pdfDstNoData;
1089 : const std::vector<std::string> &m_metadata;
1090 : const GDALColorTable *const m_poCT;
1091 :
1092 : GDALGeoTransform m_gt{};
1093 : const int m_nMaxCacheTileSize;
1094 : lru11::Cache<std::string, std::shared_ptr<GDALDataset>> m_oCacheTile;
1095 :
1096 : CPL_DISALLOW_COPY_ASSIGN(MosaicDataset)
1097 :
1098 : public:
1099 63 : MosaicDataset(const std::string &directory, const std::string &extension,
1100 : const std::string &format,
1101 : const std::vector<GDALColorInterp> &aeColorInterp,
1102 : const gdal::TileMatrixSet::TileMatrix &oTM,
1103 : const OGRSpatialReference &oSRS, int nTileMinX, int nTileMinY,
1104 : int nTileMaxX, int nTileMaxY, const std::string &convention,
1105 : int nBandsIn, GDALDataType eDT, const double *pdfDstNoData,
1106 : const std::vector<std::string> &metadata,
1107 : const GDALColorTable *poCT, int maxCacheTileSize)
1108 63 : : m_directory(directory), m_extension(extension), m_format(format),
1109 : m_aeColorInterp(aeColorInterp), m_oTM(oTM), m_oSRS(oSRS),
1110 : m_nTileMinX(nTileMinX), m_nTileMinY(nTileMinY),
1111 : m_nTileMaxX(nTileMaxX), m_nTileMaxY(nTileMaxY),
1112 : m_convention(convention), m_eDT(eDT), m_pdfDstNoData(pdfDstNoData),
1113 : m_metadata(metadata), m_poCT(poCT),
1114 : m_nMaxCacheTileSize(maxCacheTileSize),
1115 63 : m_oCacheTile(/* max_size = */ maxCacheTileSize, /* elasticity = */ 0)
1116 : {
1117 63 : nRasterXSize = (nTileMaxX - nTileMinX + 1) * oTM.mTileWidth;
1118 63 : nRasterYSize = (nTileMaxY - nTileMinY + 1) * oTM.mTileHeight;
1119 63 : m_gt[0] = oTM.mTopLeftX + nTileMinX * oTM.mResX * oTM.mTileWidth;
1120 63 : m_gt[1] = oTM.mResX;
1121 63 : m_gt[2] = 0;
1122 63 : m_gt[3] = oTM.mTopLeftY - nTileMinY * oTM.mResY * oTM.mTileHeight;
1123 63 : m_gt[4] = 0;
1124 63 : m_gt[5] = -oTM.mResY;
1125 269 : for (int i = 1; i <= nBandsIn; ++i)
1126 : {
1127 : const GDALColorInterp eColorInterp =
1128 206 : (i <= static_cast<int>(m_aeColorInterp.size()))
1129 206 : ? m_aeColorInterp[i - 1]
1130 206 : : GCI_AlphaBand;
1131 206 : SetBand(i, new MosaicRasterBand(
1132 206 : this, i, nRasterXSize, nRasterYSize, oTM.mTileWidth,
1133 206 : oTM.mTileHeight, eDT, eColorInterp, nTileMinX,
1134 : nTileMinY, oTM, convention, directory, extension,
1135 206 : pdfDstNoData, poCT));
1136 : }
1137 63 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1138 126 : const CPLStringList aosMD(metadata);
1139 80 : for (const auto [key, value] : cpl::IterateNameValue(aosMD))
1140 : {
1141 17 : SetMetadataItem(key, value);
1142 : }
1143 63 : }
1144 :
1145 134 : const OGRSpatialReference *GetSpatialRef() const override
1146 : {
1147 134 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
1148 : }
1149 :
1150 144 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
1151 : {
1152 144 : gt = m_gt;
1153 144 : return CE_None;
1154 : }
1155 :
1156 : using GDALDataset::Clone;
1157 :
1158 16 : std::unique_ptr<MosaicDataset> Clone() const
1159 : {
1160 : return std::make_unique<MosaicDataset>(
1161 16 : m_directory, m_extension, m_format, m_aeColorInterp, m_oTM, m_oSRS,
1162 16 : m_nTileMinX, m_nTileMinY, m_nTileMaxX, m_nTileMaxY, m_convention,
1163 16 : nBands, m_eDT, m_pdfDstNoData, m_metadata, m_poCT,
1164 16 : m_nMaxCacheTileSize);
1165 : }
1166 : };
1167 :
1168 : /************************************************************************/
1169 : /* MosaicRasterBand::IReadBlock() */
1170 : /************************************************************************/
1171 :
1172 2312 : CPLErr MosaicRasterBand::IReadBlock(int nXBlock, int nYBlock, void *pData)
1173 : {
1174 2312 : auto poThisDS = cpl::down_cast<MosaicDataset *>(poDS);
1175 : std::string filename = CPLFormFilenameSafe(
1176 4624 : m_directory.c_str(), CPLSPrintf("%d", m_tileMinX + nXBlock), nullptr);
1177 2312 : const int iFileY = GetFileY(m_tileMinY + nYBlock, m_oTM, m_convention);
1178 4623 : filename = CPLFormFilenameSafe(filename.c_str(), CPLSPrintf("%d", iFileY),
1179 2312 : m_extension.c_str());
1180 :
1181 2312 : std::shared_ptr<GDALDataset> poTileDS;
1182 2312 : if (!poThisDS->m_oCacheTile.tryGet(filename, poTileDS))
1183 : {
1184 665 : const char *const apszAllowedDrivers[] = {poThisDS->m_format.c_str(),
1185 665 : nullptr};
1186 665 : const char *const apszAllowedDriversForCOG[] = {"GTiff", "LIBERTIFF",
1187 : nullptr};
1188 : // CPLDebugOnly("gdal_raster_tile", "Opening %s", filename.c_str());
1189 665 : poTileDS.reset(GDALDataset::Open(
1190 : filename.c_str(), GDAL_OF_RASTER | GDAL_OF_INTERNAL,
1191 665 : EQUAL(poThisDS->m_format.c_str(), "COG") ? apszAllowedDriversForCOG
1192 : : apszAllowedDrivers));
1193 664 : if (!poTileDS)
1194 : {
1195 : VSIStatBufL sStat;
1196 6 : if (VSIStatL(filename.c_str(), &sStat) == 0)
1197 : {
1198 1 : CPLError(CE_Failure, CPLE_AppDefined,
1199 : "File %s exists but cannot be opened with %s driver",
1200 : filename.c_str(), poThisDS->m_format.c_str());
1201 1 : return CE_Failure;
1202 : }
1203 : }
1204 663 : poThisDS->m_oCacheTile.insert(filename, poTileDS);
1205 : }
1206 2311 : if (!poTileDS || nBand > poTileDS->GetRasterCount())
1207 : {
1208 1107 : memset(pData,
1209 551 : (poTileDS && (nBand == poTileDS->GetRasterCount() + 1)) ? 255
1210 : : 0,
1211 556 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
1212 556 : GDALGetDataTypeSizeBytes(eDataType));
1213 556 : return CE_None;
1214 : }
1215 : else
1216 : {
1217 1755 : return poTileDS->GetRasterBand(nBand)->RasterIO(
1218 : GF_Read, 0, 0, nBlockXSize, nBlockYSize, pData, nBlockXSize,
1219 1755 : nBlockYSize, eDataType, 0, 0, nullptr);
1220 : }
1221 : }
1222 :
1223 : } // namespace
1224 :
1225 : /************************************************************************/
1226 : /* ApplySubstitutions() */
1227 : /************************************************************************/
1228 :
1229 65 : static void ApplySubstitutions(CPLString &s,
1230 : const std::map<std::string, std::string> &substs)
1231 : {
1232 1024 : for (const auto &[key, value] : substs)
1233 : {
1234 959 : s.replaceAll("%(" + key + ")s", value);
1235 959 : s.replaceAll("%(" + key + ")d", value);
1236 959 : s.replaceAll("%(" + key + ")f", value);
1237 959 : s.replaceAll("${" + key + "}", value);
1238 : }
1239 65 : }
1240 :
1241 : /************************************************************************/
1242 : /* GenerateLeaflet() */
1243 : /************************************************************************/
1244 :
1245 15 : static void GenerateLeaflet(const std::string &osDirectory,
1246 : const std::string &osTitle, double dfSouthLat,
1247 : double dfWestLon, double dfNorthLat,
1248 : double dfEastLon, int nMinZoom, int nMaxZoom,
1249 : int nTileSize, const std::string &osExtension,
1250 : const std::string &osURL,
1251 : const std::string &osCopyright, bool bXYZ)
1252 : {
1253 15 : if (const char *pszTemplate = CPLFindFile("gdal", "leaflet_template.html"))
1254 : {
1255 30 : const std::string osFilename(pszTemplate);
1256 30 : std::map<std::string, std::string> substs;
1257 :
1258 : // For tests
1259 : const char *pszFmt =
1260 15 : atoi(CPLGetConfigOption("GDAL_RASTER_TILE_HTML_PREC", "17")) == 10
1261 : ? "%.10g"
1262 15 : : "%.17g";
1263 :
1264 30 : substs["double_quote_escaped_title"] =
1265 45 : CPLString(osTitle).replaceAll('"', "\\\"");
1266 15 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1267 15 : substs["xml_escaped_title"] = pszStr;
1268 15 : CPLFree(pszStr);
1269 15 : substs["south"] = CPLSPrintf(pszFmt, dfSouthLat);
1270 15 : substs["west"] = CPLSPrintf(pszFmt, dfWestLon);
1271 15 : substs["north"] = CPLSPrintf(pszFmt, dfNorthLat);
1272 15 : substs["east"] = CPLSPrintf(pszFmt, dfEastLon);
1273 15 : substs["centerlon"] = CPLSPrintf(pszFmt, (dfNorthLat + dfSouthLat) / 2);
1274 15 : substs["centerlat"] = CPLSPrintf(pszFmt, (dfWestLon + dfEastLon) / 2);
1275 15 : substs["minzoom"] = CPLSPrintf("%d", nMinZoom);
1276 15 : substs["maxzoom"] = CPLSPrintf("%d", nMaxZoom);
1277 15 : substs["beginzoom"] = CPLSPrintf("%d", nMaxZoom);
1278 15 : substs["tile_size"] = CPLSPrintf("%d", nTileSize); // not used
1279 15 : substs["tileformat"] = osExtension;
1280 15 : substs["publishurl"] = osURL; // not used
1281 15 : substs["copyright"] = CPLString(osCopyright).replaceAll('"', "\\\"");
1282 15 : substs["tms"] = bXYZ ? "0" : "1";
1283 :
1284 15 : GByte *pabyRet = nullptr;
1285 15 : CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, osFilename.c_str(), &pabyRet,
1286 : nullptr, 10 * 1024 * 1024));
1287 15 : if (pabyRet)
1288 : {
1289 30 : CPLString osHTML(reinterpret_cast<char *>(pabyRet));
1290 15 : CPLFree(pabyRet);
1291 :
1292 15 : ApplySubstitutions(osHTML, substs);
1293 :
1294 15 : VSILFILE *f = VSIFOpenL(CPLFormFilenameSafe(osDirectory.c_str(),
1295 : "leaflet.html", nullptr)
1296 : .c_str(),
1297 : "wb");
1298 15 : if (f)
1299 : {
1300 15 : VSIFWriteL(osHTML.data(), 1, osHTML.size(), f);
1301 15 : VSIFCloseL(f);
1302 : }
1303 : }
1304 : }
1305 15 : }
1306 :
1307 : /************************************************************************/
1308 : /* GenerateMapML() */
1309 : /************************************************************************/
1310 :
1311 : static void
1312 16 : GenerateMapML(const std::string &osDirectory, const std::string &mapmlTemplate,
1313 : const std::string &osTitle, int nMinTileX, int nMinTileY,
1314 : int nMaxTileX, int nMaxTileY, int nMinZoom, int nMaxZoom,
1315 : const std::string &osExtension, const std::string &osURL,
1316 : const std::string &osCopyright, const gdal::TileMatrixSet &tms)
1317 : {
1318 16 : if (const char *pszTemplate =
1319 16 : (mapmlTemplate.empty() ? CPLFindFile("gdal", "template_tiles.mapml")
1320 16 : : mapmlTemplate.c_str()))
1321 : {
1322 32 : const std::string osFilename(pszTemplate);
1323 32 : std::map<std::string, std::string> substs;
1324 :
1325 16 : if (tms.identifier() == "GoogleMapsCompatible")
1326 15 : substs["TILING_SCHEME"] = "OSMTILE";
1327 1 : else if (tms.identifier() == "WorldCRS84Quad")
1328 1 : substs["TILING_SCHEME"] = "WGS84";
1329 : else
1330 0 : substs["TILING_SCHEME"] = tms.identifier();
1331 :
1332 16 : substs["URL"] = osURL.empty() ? "./" : osURL;
1333 16 : substs["MINTILEX"] = CPLSPrintf("%d", nMinTileX);
1334 16 : substs["MINTILEY"] = CPLSPrintf("%d", nMinTileY);
1335 16 : substs["MAXTILEX"] = CPLSPrintf("%d", nMaxTileX);
1336 16 : substs["MAXTILEY"] = CPLSPrintf("%d", nMaxTileY);
1337 16 : substs["CURZOOM"] = CPLSPrintf("%d", nMaxZoom);
1338 16 : substs["MINZOOM"] = CPLSPrintf("%d", nMinZoom);
1339 16 : substs["MAXZOOM"] = CPLSPrintf("%d", nMaxZoom);
1340 16 : substs["TILEEXT"] = osExtension;
1341 16 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1342 16 : substs["TITLE"] = pszStr;
1343 16 : CPLFree(pszStr);
1344 16 : substs["COPYRIGHT"] = osCopyright;
1345 :
1346 16 : GByte *pabyRet = nullptr;
1347 16 : CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, osFilename.c_str(), &pabyRet,
1348 : nullptr, 10 * 1024 * 1024));
1349 16 : if (pabyRet)
1350 : {
1351 32 : CPLString osMAPML(reinterpret_cast<char *>(pabyRet));
1352 16 : CPLFree(pabyRet);
1353 :
1354 16 : ApplySubstitutions(osMAPML, substs);
1355 :
1356 16 : VSILFILE *f = VSIFOpenL(
1357 32 : CPLFormFilenameSafe(osDirectory.c_str(), "mapml.mapml", nullptr)
1358 : .c_str(),
1359 : "wb");
1360 16 : if (f)
1361 : {
1362 16 : VSIFWriteL(osMAPML.data(), 1, osMAPML.size(), f);
1363 16 : VSIFCloseL(f);
1364 : }
1365 : }
1366 : }
1367 16 : }
1368 :
1369 : /************************************************************************/
1370 : /* GenerateOpenLayers() */
1371 : /************************************************************************/
1372 :
1373 23 : static void GenerateOpenLayers(
1374 : const std::string &osDirectory, const std::string &osTitle, double dfMinX,
1375 : double dfMinY, double dfMaxX, double dfMaxY, int nMinZoom, int nMaxZoom,
1376 : int nTileSize, const std::string &osExtension, const std::string &osURL,
1377 : const std::string &osCopyright, const gdal::TileMatrixSet &tms,
1378 : bool bInvertAxisTMS, const OGRSpatialReference &oSRS_TMS, bool bXYZ)
1379 : {
1380 46 : std::map<std::string, std::string> substs;
1381 :
1382 : // For tests
1383 : const char *pszFmt =
1384 23 : atoi(CPLGetConfigOption("GDAL_RASTER_TILE_HTML_PREC", "17")) == 10
1385 : ? "%.10g"
1386 23 : : "%.17g";
1387 :
1388 23 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1389 23 : substs["xml_escaped_title"] = pszStr;
1390 23 : CPLFree(pszStr);
1391 23 : substs["ominx"] = CPLSPrintf(pszFmt, dfMinX);
1392 23 : substs["ominy"] = CPLSPrintf(pszFmt, dfMinY);
1393 23 : substs["omaxx"] = CPLSPrintf(pszFmt, dfMaxX);
1394 23 : substs["omaxy"] = CPLSPrintf(pszFmt, dfMaxY);
1395 23 : substs["center_x"] = CPLSPrintf(pszFmt, (dfMinX + dfMaxX) / 2);
1396 23 : substs["center_y"] = CPLSPrintf(pszFmt, (dfMinY + dfMaxY) / 2);
1397 23 : substs["minzoom"] = CPLSPrintf("%d", nMinZoom);
1398 23 : substs["maxzoom"] = CPLSPrintf("%d", nMaxZoom);
1399 23 : substs["tile_size"] = CPLSPrintf("%d", nTileSize);
1400 23 : substs["tileformat"] = osExtension;
1401 23 : substs["publishurl"] = osURL;
1402 23 : substs["copyright"] = osCopyright;
1403 23 : substs["sign_y"] = bXYZ ? "" : "-";
1404 :
1405 : CPLString s(R"raw(<!DOCTYPE html>
1406 : <html>
1407 : <head>
1408 : <title>%(xml_escaped_title)s</title>
1409 : <meta http-equiv="content-type" content="text/html; charset=utf-8"/>
1410 : <meta http-equiv='imagetoolbar' content='no'/>
1411 : <style type="text/css"> v\:* {behavior:url(#default#VML);}
1412 : html, body { overflow: hidden; padding: 0; height: 100%; width: 100%; font-family: 'Lucida Grande',Geneva,Arial,Verdana,sans-serif; }
1413 : body { margin: 10px; background: #fff; }
1414 : h1 { margin: 0; padding: 6px; border:0; font-size: 20pt; }
1415 : #header { height: 43px; padding: 0; background-color: #eee; border: 1px solid #888; }
1416 : #subheader { height: 12px; text-align: right; font-size: 10px; color: #555;}
1417 : #map { height: 90%; border: 1px solid #888; }
1418 : </style>
1419 : <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">
1420 : <script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@main/dist/en/v7.0.0/legacy/ol.js"></script>
1421 : <script src="https://unpkg.com/ol-layerswitcher@4.1.1"></script>
1422 : <link rel="stylesheet" href="https://unpkg.com/ol-layerswitcher@4.1.1/src/ol-layerswitcher.css" />
1423 : </head>
1424 : <body>
1425 : <div id="header"><h1>%(xml_escaped_title)s</h1></div>
1426 : <div id="subheader">Generated by <a href="https://gdal.org/programs/gdal_raster_tile.html">gdal raster tile</a> </div>
1427 : <div id="map" class="map"></div>
1428 : <div id="mouse-position"></div>
1429 : <script type="text/javascript">
1430 : var mousePositionControl = new ol.control.MousePosition({
1431 : className: 'custom-mouse-position',
1432 : target: document.getElementById('mouse-position'),
1433 : undefinedHTML: ' '
1434 : });
1435 : var map = new ol.Map({
1436 : controls: ol.control.defaults.defaults().extend([mousePositionControl]),
1437 46 : target: 'map',)raw");
1438 :
1439 31 : if (tms.identifier() == "GoogleMapsCompatible" ||
1440 8 : tms.identifier() == "WorldCRS84Quad")
1441 : {
1442 17 : s += R"raw(
1443 : layers: [
1444 : new ol.layer.Group({
1445 : title: 'Base maps',
1446 : layers: [
1447 : new ol.layer.Tile({
1448 : title: 'OpenStreetMap',
1449 : type: 'base',
1450 : visible: true,
1451 : source: new ol.source.OSM()
1452 : }),
1453 : ]
1454 : }),)raw";
1455 : }
1456 :
1457 23 : if (tms.identifier() == "GoogleMapsCompatible")
1458 : {
1459 15 : s += R"raw(new ol.layer.Group({
1460 : title: 'Overlay',
1461 : layers: [
1462 : new ol.layer.Tile({
1463 : title: 'Overlay',
1464 : // opacity: 0.7,
1465 : extent: [%(ominx)f, %(ominy)f,%(omaxx)f, %(omaxy)f],
1466 : source: new ol.source.XYZ({
1467 : attributions: '%(copyright)s',
1468 : minZoom: %(minzoom)d,
1469 : maxZoom: %(maxzoom)d,
1470 : url: './{z}/{x}/{%(sign_y)sy}.%(tileformat)s',
1471 : tileSize: [%(tile_size)d, %(tile_size)d]
1472 : })
1473 : }),
1474 : ]
1475 : }),)raw";
1476 : }
1477 8 : else if (tms.identifier() == "WorldCRS84Quad")
1478 : {
1479 2 : const double base_res = 180.0 / nTileSize;
1480 4 : std::string resolutions = "[";
1481 4 : for (int i = 0; i <= nMaxZoom; ++i)
1482 : {
1483 2 : if (i > 0)
1484 0 : resolutions += ",";
1485 2 : resolutions += CPLSPrintf(pszFmt, base_res / (1 << i));
1486 : }
1487 2 : resolutions += "]";
1488 2 : substs["resolutions"] = std::move(resolutions);
1489 :
1490 2 : if (bXYZ)
1491 : {
1492 1 : substs["origin"] = "[-180,90]";
1493 1 : substs["y_formula"] = "tileCoord[2]";
1494 : }
1495 : else
1496 : {
1497 1 : substs["origin"] = "[-180,-90]";
1498 1 : substs["y_formula"] = "- 1 - tileCoord[2]";
1499 : }
1500 :
1501 2 : s += R"raw(
1502 : new ol.layer.Group({
1503 : title: 'Overlay',
1504 : layers: [
1505 : new ol.layer.Tile({
1506 : title: 'Overlay',
1507 : // opacity: 0.7,
1508 : extent: [%(ominx)f, %(ominy)f,%(omaxx)f, %(omaxy)f],
1509 : source: new ol.source.TileImage({
1510 : attributions: '%(copyright)s',
1511 : projection: 'EPSG:4326',
1512 : minZoom: %(minzoom)d,
1513 : maxZoom: %(maxzoom)d,
1514 : tileGrid: new ol.tilegrid.TileGrid({
1515 : extent: [-180,-90,180,90],
1516 : origin: %(origin)s,
1517 : resolutions: %(resolutions)s,
1518 : tileSize: [%(tile_size)d, %(tile_size)d]
1519 : }),
1520 : tileUrlFunction: function(tileCoord) {
1521 : return ('./{z}/{x}/{y}.%(tileformat)s'
1522 : .replace('{z}', String(tileCoord[0]))
1523 : .replace('{x}', String(tileCoord[1]))
1524 : .replace('{y}', String(%(y_formula)s)));
1525 : },
1526 : })
1527 : }),
1528 : ]
1529 : }),)raw";
1530 : }
1531 : else
1532 : {
1533 12 : substs["maxres"] =
1534 12 : CPLSPrintf(pszFmt, tms.tileMatrixList()[nMinZoom].mResX);
1535 12 : std::string resolutions = "[";
1536 16 : for (int i = 0; i <= nMaxZoom; ++i)
1537 : {
1538 10 : if (i > 0)
1539 4 : resolutions += ",";
1540 10 : resolutions += CPLSPrintf(pszFmt, tms.tileMatrixList()[i].mResX);
1541 : }
1542 6 : resolutions += "]";
1543 6 : substs["resolutions"] = std::move(resolutions);
1544 :
1545 12 : std::string matrixsizes = "[";
1546 16 : for (int i = 0; i <= nMaxZoom; ++i)
1547 : {
1548 10 : if (i > 0)
1549 4 : matrixsizes += ",";
1550 : matrixsizes +=
1551 10 : CPLSPrintf("[%d,%d]", tms.tileMatrixList()[i].mMatrixWidth,
1552 20 : tms.tileMatrixList()[i].mMatrixHeight);
1553 : }
1554 6 : matrixsizes += "]";
1555 6 : substs["matrixsizes"] = std::move(matrixsizes);
1556 :
1557 6 : double dfTopLeftX = tms.tileMatrixList()[0].mTopLeftX;
1558 6 : double dfTopLeftY = tms.tileMatrixList()[0].mTopLeftY;
1559 6 : if (bInvertAxisTMS)
1560 0 : std::swap(dfTopLeftX, dfTopLeftY);
1561 :
1562 6 : if (bXYZ)
1563 : {
1564 10 : substs["origin"] =
1565 10 : CPLSPrintf("[%.17g,%.17g]", dfTopLeftX, dfTopLeftY);
1566 5 : substs["y_formula"] = "tileCoord[2]";
1567 : }
1568 : else
1569 : {
1570 2 : substs["origin"] = CPLSPrintf(
1571 : "[%.17g,%.17g]", dfTopLeftX,
1572 1 : dfTopLeftY - tms.tileMatrixList()[0].mResY *
1573 3 : tms.tileMatrixList()[0].mTileHeight);
1574 1 : substs["y_formula"] = "- 1 - tileCoord[2]";
1575 : }
1576 :
1577 12 : substs["tilegrid_extent"] =
1578 : CPLSPrintf("[%.17g,%.17g,%.17g,%.17g]", dfTopLeftX,
1579 6 : dfTopLeftY - tms.tileMatrixList()[0].mMatrixHeight *
1580 6 : tms.tileMatrixList()[0].mResY *
1581 6 : tms.tileMatrixList()[0].mTileHeight,
1582 6 : dfTopLeftX + tms.tileMatrixList()[0].mMatrixWidth *
1583 6 : tms.tileMatrixList()[0].mResX *
1584 6 : tms.tileMatrixList()[0].mTileWidth,
1585 24 : dfTopLeftY);
1586 :
1587 6 : s += R"raw(
1588 : layers: [
1589 : new ol.layer.Group({
1590 : title: 'Overlay',
1591 : layers: [
1592 : new ol.layer.Tile({
1593 : title: 'Overlay',
1594 : // opacity: 0.7,
1595 : extent: [%(ominx)f, %(ominy)f,%(omaxx)f, %(omaxy)f],
1596 : source: new ol.source.TileImage({
1597 : attributions: '%(copyright)s',
1598 : minZoom: %(minzoom)d,
1599 : maxZoom: %(maxzoom)d,
1600 : tileGrid: new ol.tilegrid.TileGrid({
1601 : extent: %(tilegrid_extent)s,
1602 : origin: %(origin)s,
1603 : resolutions: %(resolutions)s,
1604 : sizes: %(matrixsizes)s,
1605 : tileSize: [%(tile_size)d, %(tile_size)d]
1606 : }),
1607 : tileUrlFunction: function(tileCoord) {
1608 : return ('./{z}/{x}/{y}.%(tileformat)s'
1609 : .replace('{z}', String(tileCoord[0]))
1610 : .replace('{x}', String(tileCoord[1]))
1611 : .replace('{y}', String(%(y_formula)s)));
1612 : },
1613 : })
1614 : }),
1615 : ]
1616 : }),)raw";
1617 : }
1618 :
1619 23 : s += R"raw(
1620 : ],
1621 : view: new ol.View({
1622 : center: [%(center_x)f, %(center_y)f],)raw";
1623 :
1624 31 : if (tms.identifier() == "GoogleMapsCompatible" ||
1625 8 : tms.identifier() == "WorldCRS84Quad")
1626 : {
1627 17 : substs["view_zoom"] = substs["minzoom"];
1628 17 : if (tms.identifier() == "WorldCRS84Quad")
1629 : {
1630 2 : substs["view_zoom"] = CPLSPrintf("%d", nMinZoom + 1);
1631 : }
1632 :
1633 17 : s += R"raw(
1634 : zoom: %(view_zoom)d,)raw";
1635 : }
1636 : else
1637 : {
1638 6 : s += R"raw(
1639 : resolution: %(maxres)f,)raw";
1640 : }
1641 :
1642 23 : if (tms.identifier() == "WorldCRS84Quad")
1643 : {
1644 2 : s += R"raw(
1645 : projection: 'EPSG:4326',)raw";
1646 : }
1647 21 : else if (!oSRS_TMS.IsEmpty() && tms.identifier() != "GoogleMapsCompatible")
1648 : {
1649 5 : const char *pszAuthName = oSRS_TMS.GetAuthorityName(nullptr);
1650 5 : const char *pszAuthCode = oSRS_TMS.GetAuthorityCode(nullptr);
1651 5 : if (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
1652 : {
1653 3 : substs["epsg_code"] = pszAuthCode;
1654 3 : if (oSRS_TMS.IsGeographic())
1655 : {
1656 1 : substs["units"] = "deg";
1657 : }
1658 : else
1659 : {
1660 2 : const char *pszUnits = "";
1661 2 : if (oSRS_TMS.GetLinearUnits(&pszUnits) == 1.0)
1662 2 : substs["units"] = "m";
1663 : else
1664 0 : substs["units"] = pszUnits;
1665 : }
1666 3 : s += R"raw(
1667 : projection: new ol.proj.Projection({code: 'EPSG:%(epsg_code)s', units:'%(units)s'}),)raw";
1668 : }
1669 : }
1670 :
1671 23 : s += R"raw(
1672 : })
1673 : });)raw";
1674 :
1675 31 : if (tms.identifier() == "GoogleMapsCompatible" ||
1676 8 : tms.identifier() == "WorldCRS84Quad")
1677 : {
1678 17 : s += R"raw(
1679 : map.addControl(new ol.control.LayerSwitcher());)raw";
1680 : }
1681 :
1682 23 : s += R"raw(
1683 : </script>
1684 : </body>
1685 : </html>)raw";
1686 :
1687 23 : ApplySubstitutions(s, substs);
1688 :
1689 23 : VSILFILE *f = VSIFOpenL(
1690 46 : CPLFormFilenameSafe(osDirectory.c_str(), "openlayers.html", nullptr)
1691 : .c_str(),
1692 : "wb");
1693 23 : if (f)
1694 : {
1695 23 : VSIFWriteL(s.data(), 1, s.size(), f);
1696 23 : VSIFCloseL(f);
1697 : }
1698 23 : }
1699 :
1700 : /************************************************************************/
1701 : /* GetTileBoundingBox() */
1702 : /************************************************************************/
1703 :
1704 6 : static void GetTileBoundingBox(int nTileX, int nTileY, int nTileZ,
1705 : const gdal::TileMatrixSet *poTMS,
1706 : bool bInvertAxisTMS,
1707 : OGRCoordinateTransformation *poCTToWGS84,
1708 : double &dfTLX, double &dfTLY, double &dfTRX,
1709 : double &dfTRY, double &dfLLX, double &dfLLY,
1710 : double &dfLRX, double &dfLRY)
1711 : {
1712 : gdal::TileMatrixSet::TileMatrix tileMatrix =
1713 12 : poTMS->tileMatrixList()[nTileZ];
1714 6 : if (bInvertAxisTMS)
1715 0 : std::swap(tileMatrix.mTopLeftX, tileMatrix.mTopLeftY);
1716 :
1717 6 : dfTLX = tileMatrix.mTopLeftX +
1718 6 : nTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
1719 6 : dfTLY = tileMatrix.mTopLeftY -
1720 6 : nTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
1721 6 : poCTToWGS84->Transform(1, &dfTLX, &dfTLY);
1722 :
1723 6 : dfTRX = tileMatrix.mTopLeftX +
1724 6 : (nTileX + 1) * tileMatrix.mResX * tileMatrix.mTileWidth;
1725 6 : dfTRY = tileMatrix.mTopLeftY -
1726 6 : nTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
1727 6 : poCTToWGS84->Transform(1, &dfTRX, &dfTRY);
1728 :
1729 6 : dfLLX = tileMatrix.mTopLeftX +
1730 6 : nTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
1731 6 : dfLLY = tileMatrix.mTopLeftY -
1732 6 : (nTileY + 1) * tileMatrix.mResY * tileMatrix.mTileHeight;
1733 6 : poCTToWGS84->Transform(1, &dfLLX, &dfLLY);
1734 :
1735 6 : dfLRX = tileMatrix.mTopLeftX +
1736 6 : (nTileX + 1) * tileMatrix.mResX * tileMatrix.mTileWidth;
1737 6 : dfLRY = tileMatrix.mTopLeftY -
1738 6 : (nTileY + 1) * tileMatrix.mResY * tileMatrix.mTileHeight;
1739 6 : poCTToWGS84->Transform(1, &dfLRX, &dfLRY);
1740 6 : }
1741 :
1742 : /************************************************************************/
1743 : /* GenerateKML() */
1744 : /************************************************************************/
1745 :
1746 : namespace
1747 : {
1748 : struct TileCoordinates
1749 : {
1750 : int nTileX = 0;
1751 : int nTileY = 0;
1752 : int nTileZ = 0;
1753 : };
1754 : } // namespace
1755 :
1756 5 : static void GenerateKML(const std::string &osDirectory,
1757 : const std::string &osTitle, int nTileX, int nTileY,
1758 : int nTileZ, int nTileSize,
1759 : const std::string &osExtension,
1760 : const std::string &osURL,
1761 : const gdal::TileMatrixSet *poTMS, bool bInvertAxisTMS,
1762 : const std::string &convention,
1763 : OGRCoordinateTransformation *poCTToWGS84,
1764 : const std::vector<TileCoordinates> &children)
1765 : {
1766 10 : std::map<std::string, std::string> substs;
1767 :
1768 5 : const bool bIsTileKML = nTileX >= 0;
1769 :
1770 : // For tests
1771 : const char *pszFmt =
1772 5 : atoi(CPLGetConfigOption("GDAL_RASTER_TILE_KML_PREC", "14")) == 10
1773 : ? "%.10f"
1774 5 : : "%.14f";
1775 :
1776 5 : substs["tx"] = CPLSPrintf("%d", nTileX);
1777 5 : substs["tz"] = CPLSPrintf("%d", nTileZ);
1778 5 : substs["tileformat"] = osExtension;
1779 5 : substs["minlodpixels"] = CPLSPrintf("%d", nTileSize / 2);
1780 10 : substs["maxlodpixels"] =
1781 10 : children.empty() ? "-1" : CPLSPrintf("%d", nTileSize * 8);
1782 :
1783 5 : double dfTLX = 0;
1784 5 : double dfTLY = 0;
1785 5 : double dfTRX = 0;
1786 5 : double dfTRY = 0;
1787 5 : double dfLLX = 0;
1788 5 : double dfLLY = 0;
1789 5 : double dfLRX = 0;
1790 5 : double dfLRY = 0;
1791 :
1792 5 : int nFileY = -1;
1793 5 : if (!bIsTileKML)
1794 : {
1795 2 : char *pszStr = CPLEscapeString(osTitle.c_str(), -1, CPLES_XML);
1796 2 : substs["xml_escaped_title"] = pszStr;
1797 2 : CPLFree(pszStr);
1798 : }
1799 : else
1800 : {
1801 3 : nFileY = GetFileY(nTileY, poTMS->tileMatrixList()[nTileZ], convention);
1802 3 : substs["realtiley"] = CPLSPrintf("%d", nFileY);
1803 6 : substs["xml_escaped_title"] =
1804 6 : CPLSPrintf("%d/%d/%d.kml", nTileZ, nTileX, nFileY);
1805 :
1806 3 : GetTileBoundingBox(nTileX, nTileY, nTileZ, poTMS, bInvertAxisTMS,
1807 : poCTToWGS84, dfTLX, dfTLY, dfTRX, dfTRY, dfLLX,
1808 : dfLLY, dfLRX, dfLRY);
1809 : }
1810 :
1811 10 : substs["drawOrder"] = CPLSPrintf("%d", nTileX == 0 ? 2 * nTileZ + 1
1812 4 : : nTileX > 0 ? 2 * nTileZ
1813 14 : : 0);
1814 :
1815 5 : substs["url"] = osURL.empty() && bIsTileKML ? "../../" : "";
1816 :
1817 5 : const bool bIsRectangle =
1818 5 : (dfTLX == dfLLX && dfTRX == dfLRX && dfTLY == dfTRY && dfLLY == dfLRY);
1819 5 : const bool bUseGXNamespace = bIsTileKML && !bIsRectangle;
1820 :
1821 10 : substs["xmlns_gx"] = bUseGXNamespace
1822 : ? " xmlns:gx=\"http://www.google.com/kml/ext/2.2\""
1823 10 : : "";
1824 :
1825 : CPLString s(R"raw(<?xml version="1.0" encoding="utf-8"?>
1826 : <kml xmlns="http://www.opengis.net/kml/2.2"%(xmlns_gx)s>
1827 : <Document>
1828 : <name>%(xml_escaped_title)s</name>
1829 : <description></description>
1830 : <Style>
1831 : <ListStyle id="hideChildren">
1832 : <listItemType>checkHideChildren</listItemType>
1833 : </ListStyle>
1834 : </Style>
1835 10 : )raw");
1836 5 : ApplySubstitutions(s, substs);
1837 :
1838 5 : if (bIsTileKML)
1839 : {
1840 : CPLString s2(R"raw( <Region>
1841 : <LatLonAltBox>
1842 : <north>%(north)f</north>
1843 : <south>%(south)f</south>
1844 : <east>%(east)f</east>
1845 : <west>%(west)f</west>
1846 : </LatLonAltBox>
1847 : <Lod>
1848 : <minLodPixels>%(minlodpixels)d</minLodPixels>
1849 : <maxLodPixels>%(maxlodpixels)d</maxLodPixels>
1850 : </Lod>
1851 : </Region>
1852 : <GroundOverlay>
1853 : <drawOrder>%(drawOrder)d</drawOrder>
1854 : <Icon>
1855 : <href>%(realtiley)d.%(tileformat)s</href>
1856 : </Icon>
1857 : <LatLonBox>
1858 : <north>%(north)f</north>
1859 : <south>%(south)f</south>
1860 : <east>%(east)f</east>
1861 : <west>%(west)f</west>
1862 : </LatLonBox>
1863 6 : )raw");
1864 :
1865 3 : if (!bIsRectangle)
1866 : {
1867 : s2 +=
1868 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>
1869 : )raw";
1870 : }
1871 :
1872 3 : s2 += R"raw( </GroundOverlay>
1873 : )raw";
1874 3 : substs["north"] = CPLSPrintf(pszFmt, std::max(dfTLY, dfTRY));
1875 3 : substs["south"] = CPLSPrintf(pszFmt, std::min(dfLLY, dfLRY));
1876 3 : substs["east"] = CPLSPrintf(pszFmt, std::max(dfTRX, dfLRX));
1877 3 : substs["west"] = CPLSPrintf(pszFmt, std::min(dfLLX, dfTLX));
1878 :
1879 3 : if (!bIsRectangle)
1880 : {
1881 1 : substs["TLX"] = CPLSPrintf(pszFmt, dfTLX);
1882 1 : substs["TLY"] = CPLSPrintf(pszFmt, dfTLY);
1883 1 : substs["TRX"] = CPLSPrintf(pszFmt, dfTRX);
1884 1 : substs["TRY"] = CPLSPrintf(pszFmt, dfTRY);
1885 1 : substs["LRX"] = CPLSPrintf(pszFmt, dfLRX);
1886 1 : substs["LRY"] = CPLSPrintf(pszFmt, dfLRY);
1887 1 : substs["LLX"] = CPLSPrintf(pszFmt, dfLLX);
1888 1 : substs["LLY"] = CPLSPrintf(pszFmt, dfLLY);
1889 : }
1890 :
1891 3 : ApplySubstitutions(s2, substs);
1892 3 : s += s2;
1893 : }
1894 :
1895 8 : for (const auto &child : children)
1896 : {
1897 3 : substs["tx"] = CPLSPrintf("%d", child.nTileX);
1898 3 : substs["tz"] = CPLSPrintf("%d", child.nTileZ);
1899 6 : substs["realtiley"] = CPLSPrintf(
1900 3 : "%d", GetFileY(child.nTileY, poTMS->tileMatrixList()[child.nTileZ],
1901 6 : convention));
1902 :
1903 3 : GetTileBoundingBox(child.nTileX, child.nTileY, child.nTileZ, poTMS,
1904 : bInvertAxisTMS, poCTToWGS84, dfTLX, dfTLY, dfTRX,
1905 : dfTRY, dfLLX, dfLLY, dfLRX, dfLRY);
1906 :
1907 : CPLString s2(R"raw( <NetworkLink>
1908 : <name>%(tz)d/%(tx)d/%(realtiley)d.%(tileformat)s</name>
1909 : <Region>
1910 : <LatLonAltBox>
1911 : <north>%(north)f</north>
1912 : <south>%(south)f</south>
1913 : <east>%(east)f</east>
1914 : <west>%(west)f</west>
1915 : </LatLonAltBox>
1916 : <Lod>
1917 : <minLodPixels>%(minlodpixels)d</minLodPixels>
1918 : <maxLodPixels>-1</maxLodPixels>
1919 : </Lod>
1920 : </Region>
1921 : <Link>
1922 : <href>%(url)s%(tz)d/%(tx)d/%(realtiley)d.kml</href>
1923 : <viewRefreshMode>onRegion</viewRefreshMode>
1924 : <viewFormat/>
1925 : </Link>
1926 : </NetworkLink>
1927 6 : )raw");
1928 3 : substs["north"] = CPLSPrintf(pszFmt, std::max(dfTLY, dfTRY));
1929 3 : substs["south"] = CPLSPrintf(pszFmt, std::min(dfLLY, dfLRY));
1930 3 : substs["east"] = CPLSPrintf(pszFmt, std::max(dfTRX, dfLRX));
1931 3 : substs["west"] = CPLSPrintf(pszFmt, std::min(dfLLX, dfTLX));
1932 3 : ApplySubstitutions(s2, substs);
1933 3 : s += s2;
1934 : }
1935 :
1936 5 : s += R"raw(</Document>
1937 : </kml>)raw";
1938 :
1939 10 : std::string osFilename(osDirectory);
1940 5 : if (!bIsTileKML)
1941 : {
1942 : osFilename =
1943 2 : CPLFormFilenameSafe(osFilename.c_str(), "doc.kml", nullptr);
1944 : }
1945 : else
1946 : {
1947 6 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
1948 3 : CPLSPrintf("%d", nTileZ), nullptr);
1949 6 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
1950 3 : CPLSPrintf("%d", nTileX), nullptr);
1951 6 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
1952 3 : CPLSPrintf("%d.kml", nFileY), nullptr);
1953 : }
1954 :
1955 5 : VSILFILE *f = VSIFOpenL(osFilename.c_str(), "wb");
1956 5 : if (f)
1957 : {
1958 5 : VSIFWriteL(s.data(), 1, s.size(), f);
1959 5 : VSIFCloseL(f);
1960 : }
1961 5 : }
1962 :
1963 : namespace
1964 : {
1965 :
1966 : /************************************************************************/
1967 : /* ResourceManager */
1968 : /************************************************************************/
1969 :
1970 : // Generic cache managing resources
1971 : template <class Resource> class ResourceManager /* non final */
1972 : {
1973 : public:
1974 113 : virtual ~ResourceManager() = default;
1975 :
1976 24 : std::unique_ptr<Resource> AcquireResources()
1977 : {
1978 48 : std::lock_guard oLock(m_oMutex);
1979 24 : if (!m_oResources.empty())
1980 : {
1981 0 : auto ret = std::move(m_oResources.back());
1982 0 : m_oResources.pop_back();
1983 0 : return ret;
1984 : }
1985 :
1986 24 : return CreateResources();
1987 : }
1988 :
1989 24 : void ReleaseResources(std::unique_ptr<Resource> resources)
1990 : {
1991 48 : std::lock_guard oLock(m_oMutex);
1992 24 : m_oResources.push_back(std::move(resources));
1993 24 : }
1994 :
1995 0 : void SetError()
1996 : {
1997 0 : std::lock_guard oLock(m_oMutex);
1998 0 : if (m_errorMsg.empty())
1999 0 : m_errorMsg = CPLGetLastErrorMsg();
2000 0 : }
2001 :
2002 6 : const std::string &GetErrorMsg() const
2003 : {
2004 6 : std::lock_guard oLock(m_oMutex);
2005 12 : return m_errorMsg;
2006 : }
2007 :
2008 : protected:
2009 : virtual std::unique_ptr<Resource> CreateResources() = 0;
2010 :
2011 : private:
2012 : mutable std::mutex m_oMutex{};
2013 : std::vector<std::unique_ptr<Resource>> m_oResources{};
2014 : std::string m_errorMsg{};
2015 : };
2016 :
2017 : /************************************************************************/
2018 : /* PerThreadMaxZoomResources */
2019 : /************************************************************************/
2020 :
2021 : // Per-thread resources for generation of tiles at full resolution
2022 : struct PerThreadMaxZoomResources
2023 : {
2024 : struct GDALDatasetReleaser
2025 : {
2026 8 : void operator()(GDALDataset *poDS)
2027 : {
2028 8 : if (poDS)
2029 8 : poDS->ReleaseRef();
2030 8 : }
2031 : };
2032 :
2033 : std::unique_ptr<GDALDataset, GDALDatasetReleaser> poSrcDS{};
2034 : std::vector<GByte> dstBuffer{};
2035 : std::unique_ptr<FakeMaxZoomDataset> poFakeMaxZoomDS{};
2036 : std::unique_ptr<void, decltype(&GDALDestroyTransformer)> poTransformer{
2037 : nullptr, GDALDestroyTransformer};
2038 : std::unique_ptr<GDALWarpOperation> poWO{};
2039 : };
2040 :
2041 : /************************************************************************/
2042 : /* PerThreadMaxZoomResourceManager */
2043 : /************************************************************************/
2044 :
2045 : // Manage a cache of PerThreadMaxZoomResources instances
2046 : class PerThreadMaxZoomResourceManager final
2047 : : public ResourceManager<PerThreadMaxZoomResources>
2048 : {
2049 : public:
2050 66 : PerThreadMaxZoomResourceManager(GDALDataset *poSrcDS,
2051 : const GDALWarpOptions *psWO,
2052 : void *pTransformerArg,
2053 : const FakeMaxZoomDataset &oFakeMaxZoomDS,
2054 : size_t nBufferSize)
2055 66 : : m_poSrcDS(poSrcDS), m_psWOSource(psWO),
2056 : m_pTransformerArg(pTransformerArg), m_oFakeMaxZoomDS(oFakeMaxZoomDS),
2057 66 : m_nBufferSize(nBufferSize)
2058 : {
2059 66 : }
2060 :
2061 : protected:
2062 8 : std::unique_ptr<PerThreadMaxZoomResources> CreateResources() override
2063 : {
2064 16 : auto ret = std::make_unique<PerThreadMaxZoomResources>();
2065 :
2066 8 : ret->poSrcDS.reset(GDALGetThreadSafeDataset(m_poSrcDS, GDAL_OF_RASTER));
2067 8 : if (!ret->poSrcDS)
2068 0 : return nullptr;
2069 :
2070 : try
2071 : {
2072 8 : ret->dstBuffer.resize(m_nBufferSize);
2073 : }
2074 0 : catch (const std::exception &)
2075 : {
2076 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
2077 : "Out of memory allocating temporary buffer");
2078 0 : return nullptr;
2079 : }
2080 :
2081 8 : ret->poFakeMaxZoomDS = m_oFakeMaxZoomDS.Clone(ret->dstBuffer);
2082 :
2083 8 : ret->poTransformer.reset(GDALCloneTransformer(m_pTransformerArg));
2084 8 : if (!ret->poTransformer)
2085 0 : return nullptr;
2086 :
2087 : auto psWO =
2088 : std::unique_ptr<GDALWarpOptions, decltype(&GDALDestroyWarpOptions)>(
2089 16 : GDALCloneWarpOptions(m_psWOSource), GDALDestroyWarpOptions);
2090 8 : if (!psWO)
2091 0 : return nullptr;
2092 :
2093 8 : psWO->hSrcDS = GDALDataset::ToHandle(ret->poSrcDS.get());
2094 8 : psWO->hDstDS = GDALDataset::ToHandle(ret->poFakeMaxZoomDS.get());
2095 8 : psWO->pTransformerArg = ret->poTransformer.get();
2096 8 : psWO->pfnTransformer = m_psWOSource->pfnTransformer;
2097 :
2098 8 : ret->poWO = std::make_unique<GDALWarpOperation>();
2099 8 : if (ret->poWO->Initialize(psWO.get()) != CE_None)
2100 0 : return nullptr;
2101 :
2102 8 : return ret;
2103 : }
2104 :
2105 : private:
2106 : GDALDataset *const m_poSrcDS;
2107 : const GDALWarpOptions *const m_psWOSource;
2108 : void *const m_pTransformerArg;
2109 : const FakeMaxZoomDataset &m_oFakeMaxZoomDS;
2110 : const size_t m_nBufferSize;
2111 :
2112 : CPL_DISALLOW_COPY_ASSIGN(PerThreadMaxZoomResourceManager)
2113 : };
2114 :
2115 : /************************************************************************/
2116 : /* PerThreadLowerZoomResources */
2117 : /************************************************************************/
2118 :
2119 : // Per-thread resources for generation of tiles at zoom level < max
2120 : struct PerThreadLowerZoomResources
2121 : {
2122 : std::unique_ptr<GDALDataset> poSrcDS{};
2123 : };
2124 :
2125 : /************************************************************************/
2126 : /* PerThreadLowerZoomResourceManager */
2127 : /************************************************************************/
2128 :
2129 : // Manage a cache of PerThreadLowerZoomResources instances
2130 : class PerThreadLowerZoomResourceManager final
2131 : : public ResourceManager<PerThreadLowerZoomResources>
2132 : {
2133 : public:
2134 47 : explicit PerThreadLowerZoomResourceManager(const MosaicDataset &oSrcDS)
2135 47 : : m_oSrcDS(oSrcDS)
2136 : {
2137 47 : }
2138 :
2139 : protected:
2140 16 : std::unique_ptr<PerThreadLowerZoomResources> CreateResources() override
2141 : {
2142 16 : auto ret = std::make_unique<PerThreadLowerZoomResources>();
2143 16 : ret->poSrcDS = m_oSrcDS.Clone();
2144 16 : return ret;
2145 : }
2146 :
2147 : private:
2148 : const MosaicDataset &m_oSrcDS;
2149 : };
2150 :
2151 : } // namespace
2152 :
2153 : /************************************************************************/
2154 : /* GDALRasterTileAlgorithm::ValidateOutputFormat() */
2155 : /************************************************************************/
2156 :
2157 101 : bool GDALRasterTileAlgorithm::ValidateOutputFormat(GDALDataType eSrcDT) const
2158 : {
2159 101 : if (m_outputFormat == "PNG")
2160 : {
2161 77 : if (m_poSrcDS->GetRasterCount() > 4)
2162 : {
2163 1 : ReportError(CE_Failure, CPLE_NotSupported,
2164 : "Only up to 4 bands supported for PNG.");
2165 1 : return false;
2166 : }
2167 76 : if (eSrcDT != GDT_Byte && eSrcDT != GDT_UInt16)
2168 : {
2169 1 : ReportError(CE_Failure, CPLE_NotSupported,
2170 : "Only Byte and UInt16 data types supported for PNG.");
2171 1 : return false;
2172 : }
2173 : }
2174 24 : else if (m_outputFormat == "JPEG")
2175 : {
2176 5 : if (m_poSrcDS->GetRasterCount() > 4)
2177 : {
2178 1 : ReportError(
2179 : CE_Failure, CPLE_NotSupported,
2180 : "Only up to 4 bands supported for JPEG (with alpha ignored).");
2181 1 : return false;
2182 : }
2183 : const bool bUInt16Supported =
2184 4 : strstr(m_poDstDriver->GetMetadataItem(GDAL_DMD_CREATIONDATATYPES),
2185 4 : "UInt16");
2186 4 : if (eSrcDT != GDT_Byte && !(eSrcDT == GDT_UInt16 && bUInt16Supported))
2187 : {
2188 1 : ReportError(
2189 : CE_Failure, CPLE_NotSupported,
2190 : bUInt16Supported
2191 : ? "Only Byte and UInt16 data types supported for JPEG."
2192 : : "Only Byte data type supported for JPEG.");
2193 1 : return false;
2194 : }
2195 3 : if (eSrcDT == GDT_UInt16)
2196 : {
2197 3 : if (const char *pszNBITS =
2198 6 : m_poSrcDS->GetRasterBand(1)->GetMetadataItem(
2199 3 : "NBITS", "IMAGE_STRUCTURE"))
2200 : {
2201 1 : if (atoi(pszNBITS) > 12)
2202 : {
2203 1 : ReportError(CE_Failure, CPLE_NotSupported,
2204 : "JPEG output only supported up to 12 bits");
2205 1 : return false;
2206 : }
2207 : }
2208 : else
2209 : {
2210 2 : double adfMinMax[2] = {0, 0};
2211 2 : m_poSrcDS->GetRasterBand(1)->ComputeRasterMinMax(
2212 2 : /* bApproxOK = */ true, adfMinMax);
2213 2 : if (adfMinMax[1] >= (1 << 12))
2214 : {
2215 1 : ReportError(CE_Failure, CPLE_NotSupported,
2216 : "JPEG output only supported up to 12 bits");
2217 1 : return false;
2218 : }
2219 : }
2220 : }
2221 : }
2222 19 : else if (m_outputFormat == "WEBP")
2223 : {
2224 3 : if (m_poSrcDS->GetRasterCount() != 3 &&
2225 1 : m_poSrcDS->GetRasterCount() != 4)
2226 : {
2227 1 : ReportError(CE_Failure, CPLE_NotSupported,
2228 : "Only 3 or 4 bands supported for WEBP.");
2229 1 : return false;
2230 : }
2231 1 : if (eSrcDT != GDT_Byte)
2232 : {
2233 1 : ReportError(CE_Failure, CPLE_NotSupported,
2234 : "Only Byte data type supported for WEBP.");
2235 1 : return false;
2236 : }
2237 : }
2238 93 : return true;
2239 : }
2240 :
2241 : /************************************************************************/
2242 : /* GDALRasterTileAlgorithm::ComputeJobChunkSize() */
2243 : /************************************************************************/
2244 :
2245 : // Given a number of tiles in the Y dimension being nTilesPerCol and
2246 : // in the X dimension being nTilesPerRow, compute the (upper bound of)
2247 : // number of jobs needed to be nYOuterIterations x nXOuterIterations,
2248 : // with each job processing in average dfTilesYPerJob x dfTilesXPerJob
2249 : // tiles.
2250 : /* static */
2251 10 : void GDALRasterTileAlgorithm::ComputeJobChunkSize(
2252 : int nMaxJobCount, int nTilesPerCol, int nTilesPerRow,
2253 : double &dfTilesYPerJob, int &nYOuterIterations, double &dfTilesXPerJob,
2254 : int &nXOuterIterations)
2255 : {
2256 10 : CPLAssert(nMaxJobCount >= 1);
2257 10 : dfTilesYPerJob = static_cast<double>(nTilesPerCol) / nMaxJobCount;
2258 10 : nYOuterIterations = dfTilesYPerJob >= 1 ? nMaxJobCount : 1;
2259 :
2260 20 : dfTilesXPerJob = dfTilesYPerJob >= 1
2261 10 : ? nTilesPerRow
2262 3 : : static_cast<double>(nTilesPerRow) / nMaxJobCount;
2263 10 : nXOuterIterations = dfTilesYPerJob >= 1 ? 1 : nMaxJobCount;
2264 :
2265 10 : if (dfTilesYPerJob < 1 && dfTilesXPerJob < 1 &&
2266 3 : nTilesPerCol <= nMaxJobCount / nTilesPerRow)
2267 : {
2268 3 : dfTilesYPerJob = 1;
2269 3 : dfTilesXPerJob = 1;
2270 3 : nYOuterIterations = nTilesPerCol;
2271 3 : nXOuterIterations = nTilesPerRow;
2272 : }
2273 10 : }
2274 :
2275 : /************************************************************************/
2276 : /* GDALRasterTileAlgorithm::AddArgToArgv() */
2277 : /************************************************************************/
2278 :
2279 48 : bool GDALRasterTileAlgorithm::AddArgToArgv(const GDALAlgorithmArg *arg,
2280 : CPLStringList &aosArgv) const
2281 : {
2282 48 : aosArgv.push_back(CPLSPrintf("--%s", arg->GetName().c_str()));
2283 48 : if (arg->GetType() == GAAT_STRING)
2284 : {
2285 14 : aosArgv.push_back(arg->Get<std::string>().c_str());
2286 : }
2287 34 : else if (arg->GetType() == GAAT_STRING_LIST)
2288 : {
2289 12 : bool bFirst = true;
2290 24 : for (const std::string &s : arg->Get<std::vector<std::string>>())
2291 : {
2292 12 : if (!bFirst)
2293 : {
2294 0 : aosArgv.push_back(CPLSPrintf("--%s", arg->GetName().c_str()));
2295 : }
2296 12 : bFirst = false;
2297 12 : aosArgv.push_back(s.c_str());
2298 : }
2299 : }
2300 22 : else if (arg->GetType() == GAAT_REAL)
2301 : {
2302 0 : aosArgv.push_back(CPLSPrintf("%.17g", arg->Get<double>()));
2303 : }
2304 22 : else if (arg->GetType() == GAAT_INTEGER)
2305 : {
2306 22 : aosArgv.push_back(CPLSPrintf("%d", arg->Get<int>()));
2307 : }
2308 0 : else if (arg->GetType() != GAAT_BOOLEAN)
2309 : {
2310 0 : ReportError(CE_Failure, CPLE_AppDefined,
2311 : "Bug: argument of type %d not handled "
2312 : "by gdal raster tile!",
2313 0 : static_cast<int>(arg->GetType()));
2314 0 : return false;
2315 : }
2316 48 : return true;
2317 : }
2318 :
2319 : /************************************************************************/
2320 : /* GDALRasterTileAlgorithm::IsCompatibleOfSpawn() */
2321 : /************************************************************************/
2322 :
2323 8 : bool GDALRasterTileAlgorithm::IsCompatibleOfSpawn(const char *&pszErrorMsg)
2324 : {
2325 8 : pszErrorMsg = "";
2326 8 : if (!m_bIsNamedNonMemSrcDS)
2327 : {
2328 1 : pszErrorMsg = "Unnamed or memory dataset sources are not supported "
2329 : "with spawn parallelization method";
2330 1 : return false;
2331 : }
2332 7 : if (cpl::starts_with(m_outputDirectory, "/vsimem/"))
2333 : {
2334 4 : pszErrorMsg = "/vsimem/ output directory not supported with spawn "
2335 : "parallelization method";
2336 4 : return false;
2337 : }
2338 :
2339 3 : if (m_osGDALPath.empty())
2340 3 : m_osGDALPath = GDALGetGDALPath();
2341 3 : return !(m_osGDALPath.empty());
2342 : }
2343 :
2344 : /************************************************************************/
2345 : /* GetProgressForChildProcesses() */
2346 : /************************************************************************/
2347 :
2348 4 : static void GetProgressForChildProcesses(
2349 : bool &bRet, std::vector<CPLSpawnedProcess *> &ahSpawnedProcesses,
2350 : std::vector<uint64_t> &anRemainingTilesForProcess, uint64_t &nCurTile,
2351 : uint64_t nTotalTiles, GDALProgressFunc pfnProgress, void *pProgressData)
2352 : {
2353 8 : std::vector<unsigned int> anProgressState(ahSpawnedProcesses.size(), 0);
2354 :
2355 91 : while (bRet)
2356 : {
2357 91 : size_t iProcess = 0;
2358 91 : size_t nFinished = 0;
2359 447 : for (CPLSpawnedProcess *hSpawnedProcess : ahSpawnedProcesses)
2360 : {
2361 356 : char ch = 0;
2362 700 : if (anRemainingTilesForProcess[iProcess] == 0 ||
2363 344 : !CPLPipeRead(CPLSpawnAsyncGetInputFileHandle(hSpawnedProcess),
2364 : &ch, 1))
2365 : {
2366 14 : ++nFinished;
2367 : }
2368 342 : else if (ch == PROGRESS_MARKER[anProgressState[iProcess]])
2369 : {
2370 258 : ++anProgressState[iProcess];
2371 258 : if (anProgressState[iProcess] == sizeof(PROGRESS_MARKER))
2372 : {
2373 86 : anProgressState[iProcess] = 0;
2374 86 : --anRemainingTilesForProcess[iProcess];
2375 86 : ++nCurTile;
2376 170 : bRet &= (!pfnProgress ||
2377 84 : pfnProgress(static_cast<double>(nCurTile) /
2378 84 : static_cast<double>(nTotalTiles),
2379 86 : "", pProgressData));
2380 : }
2381 : }
2382 : else
2383 : {
2384 84 : CPLErrorOnce(
2385 : CE_Warning, CPLE_AppDefined,
2386 : "Spurious character detected on stdout of child process");
2387 84 : anProgressState[iProcess] = 0;
2388 84 : if (ch == PROGRESS_MARKER[anProgressState[iProcess]])
2389 : {
2390 84 : ++anProgressState[iProcess];
2391 : }
2392 : }
2393 356 : ++iProcess;
2394 : }
2395 91 : if (!bRet || nFinished == ahSpawnedProcesses.size())
2396 4 : break;
2397 : }
2398 4 : }
2399 :
2400 : /************************************************************************/
2401 : /* WaitForSpawnedProcesses() */
2402 : /************************************************************************/
2403 :
2404 4 : void GDALRasterTileAlgorithm::WaitForSpawnedProcesses(
2405 : bool &bRet, const std::vector<std::string> &asCommandLines,
2406 : std::vector<CPLSpawnedProcess *> &ahSpawnedProcesses) const
2407 : {
2408 4 : size_t iProcess = 0;
2409 18 : for (CPLSpawnedProcess *hSpawnedProcess : ahSpawnedProcesses)
2410 : {
2411 14 : CPLSpawnAsyncCloseInputFileHandle(hSpawnedProcess);
2412 :
2413 14 : char ch = 0;
2414 14 : std::string errorMsg;
2415 483 : while (CPLPipeRead(CPLSpawnAsyncGetErrorFileHandle(hSpawnedProcess),
2416 483 : &ch, 1))
2417 : {
2418 469 : if (ch == '\n')
2419 : {
2420 6 : if (!errorMsg.empty())
2421 : {
2422 6 : if (cpl::starts_with(errorMsg, "ERROR "))
2423 : {
2424 6 : const auto nPos = errorMsg.find(": ");
2425 6 : if (nPos != std::string::npos)
2426 6 : errorMsg = errorMsg.substr(nPos + 1);
2427 6 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
2428 : errorMsg.c_str());
2429 : }
2430 : else
2431 : {
2432 0 : std::string osComp = "GDAL";
2433 0 : const auto nPos = errorMsg.find(": ");
2434 0 : if (nPos != std::string::npos)
2435 : {
2436 0 : osComp = errorMsg.substr(0, nPos);
2437 0 : errorMsg = errorMsg.substr(nPos + 1);
2438 : }
2439 0 : CPLDebug(osComp.c_str(), "%s", errorMsg.c_str());
2440 : }
2441 6 : errorMsg.clear();
2442 : }
2443 : }
2444 : else
2445 : {
2446 463 : errorMsg += ch;
2447 : }
2448 : }
2449 14 : CPLSpawnAsyncCloseErrorFileHandle(hSpawnedProcess);
2450 :
2451 14 : if (CPLSpawnAsyncFinish(hSpawnedProcess, /* bWait = */ true,
2452 14 : /* bKill = */ false) != 0)
2453 : {
2454 2 : bRet = false;
2455 2 : ReportError(CE_Failure, CPLE_AppDefined,
2456 : "Child process '%s' failed",
2457 2 : asCommandLines[iProcess].c_str());
2458 : }
2459 14 : ++iProcess;
2460 : }
2461 4 : }
2462 :
2463 : /************************************************************************/
2464 : /* GDALRasterTileAlgorithm::GetMaxChildCount() */
2465 : /**********************************f**************************************/
2466 :
2467 4 : int GDALRasterTileAlgorithm::GetMaxChildCount(int nMaxJobCount) const
2468 : {
2469 : #ifndef _WIN32
2470 : // Limit the number of jobs compared to how many file descriptors we have
2471 : // left
2472 : const int remainingFileDescriptorCount =
2473 4 : CPLGetRemainingFileDescriptorCount();
2474 4 : constexpr int SOME_MARGIN = 3;
2475 4 : constexpr int FD_PER_CHILD = 3; /* stdin, stdout and stderr */
2476 4 : if (FD_PER_CHILD * nMaxJobCount + SOME_MARGIN >
2477 : remainingFileDescriptorCount)
2478 : {
2479 0 : nMaxJobCount = std::max(
2480 0 : 1, (remainingFileDescriptorCount - SOME_MARGIN) / FD_PER_CHILD);
2481 0 : ReportError(
2482 : CE_Warning, CPLE_AppDefined,
2483 : "Limiting the number of child workers to %d (instead of %d), "
2484 : "because there are not enough file descriptors left (%d)",
2485 0 : nMaxJobCount, m_numThreads, remainingFileDescriptorCount);
2486 : }
2487 : #endif
2488 4 : return nMaxJobCount;
2489 : }
2490 :
2491 : /************************************************************************/
2492 : /* SendConfigOptions() */
2493 : /************************************************************************/
2494 :
2495 14 : static void SendConfigOptions(CPLSpawnedProcess *hSpawnedProcess, bool &bRet)
2496 : {
2497 : // Send most config options through pipe, to avoid leaking
2498 : // secrets when listing processes
2499 14 : auto handle = CPLSpawnAsyncGetOutputFileHandle(hSpawnedProcess);
2500 42 : for (auto pfnFunc : {&CPLGetConfigOptions, &CPLGetThreadLocalConfigOptions})
2501 : {
2502 56 : CPLStringList aosConfigOptions((*pfnFunc)());
2503 78 : for (const char *pszNameValue : aosConfigOptions)
2504 : {
2505 50 : if (!STARTS_WITH(pszNameValue, "GDAL_CACHEMAX") &&
2506 50 : !STARTS_WITH(pszNameValue, "GDAL_NUM_THREADS"))
2507 : {
2508 50 : constexpr const char *CONFIG_MARKER = "--config\n";
2509 50 : bRet &= CPL_TO_BOOL(
2510 : CPLPipeWrite(handle, CONFIG_MARKER,
2511 50 : static_cast<int>(strlen(CONFIG_MARKER))));
2512 50 : char *pszEscaped = CPLEscapeString(pszNameValue, -1, CPLES_URL);
2513 50 : bRet &= CPL_TO_BOOL(CPLPipeWrite(
2514 50 : handle, pszEscaped, static_cast<int>(strlen(pszEscaped))));
2515 50 : CPLFree(pszEscaped);
2516 50 : bRet &= CPL_TO_BOOL(CPLPipeWrite(handle, "\n", 1));
2517 : }
2518 : }
2519 : }
2520 14 : constexpr const char *END_MARKER = "END\n";
2521 14 : bRet &= CPL_TO_BOOL(
2522 14 : CPLPipeWrite(handle, END_MARKER, static_cast<int>(strlen(END_MARKER))));
2523 14 : CPLSpawnAsyncCloseOutputFileHandle(hSpawnedProcess);
2524 14 : }
2525 :
2526 : /************************************************************************/
2527 : /* GDALRasterTileAlgorithm::GenerateBaseTilesSpawnMethod() */
2528 : /************************************************************************/
2529 :
2530 2 : bool GDALRasterTileAlgorithm::GenerateBaseTilesSpawnMethod(
2531 : int nBaseTilesPerCol, int nBaseTilesPerRow, int nMinTileX, int nMinTileY,
2532 : int nMaxTileX, int nMaxTileY, uint64_t nTotalTiles, uint64_t nBaseTiles,
2533 : GDALProgressFunc pfnProgress, void *pProgressData)
2534 : {
2535 2 : CPLAssert(!m_osGDALPath.empty());
2536 :
2537 2 : const int nMaxJobCount = GetMaxChildCount(std::max(
2538 0 : 1, static_cast<int>(std::min<uint64_t>(
2539 2 : m_numThreads, nBaseTiles / GetThresholdMinTilesPerJob()))));
2540 :
2541 : double dfTilesYPerJob;
2542 : int nYOuterIterations;
2543 : double dfTilesXPerJob;
2544 : int nXOuterIterations;
2545 2 : ComputeJobChunkSize(nMaxJobCount, nBaseTilesPerCol, nBaseTilesPerRow,
2546 : dfTilesYPerJob, nYOuterIterations, dfTilesXPerJob,
2547 : nXOuterIterations);
2548 :
2549 2 : CPLDebugOnly("gdal_raster_tile",
2550 : "nYOuterIterations=%d, dfTilesYPerJob=%g, "
2551 : "nXOuterIterations=%d, dfTilesXPerJob=%g",
2552 : nYOuterIterations, dfTilesYPerJob, nXOuterIterations,
2553 : dfTilesXPerJob);
2554 :
2555 4 : std::vector<std::string> asCommandLines;
2556 4 : std::vector<CPLSpawnedProcess *> ahSpawnedProcesses;
2557 4 : std::vector<uint64_t> anRemainingTilesForProcess;
2558 :
2559 2 : const uint64_t nCacheMaxPerProcess = GDALGetCacheMax64() / nMaxJobCount;
2560 :
2561 2 : int nLastYEndIncluded = nMinTileY - 1;
2562 :
2563 2 : bool bRet = true;
2564 8 : for (int iYOuterIter = 0; bRet && iYOuterIter < nYOuterIterations &&
2565 6 : nLastYEndIncluded < nMaxTileY;
2566 : ++iYOuterIter)
2567 : {
2568 6 : const int iYStart = nLastYEndIncluded + 1;
2569 : const int iYEndIncluded =
2570 6 : iYOuterIter + 1 == nYOuterIterations
2571 10 : ? nMaxTileY
2572 : : std::max(
2573 : iYStart,
2574 10 : static_cast<int>(std::floor(
2575 4 : nMinTileY + (iYOuterIter + 1) * dfTilesYPerJob - 1)));
2576 :
2577 6 : nLastYEndIncluded = iYEndIncluded;
2578 :
2579 6 : int nLastXEndIncluded = nMinTileX - 1;
2580 12 : for (int iXOuterIter = 0; bRet && iXOuterIter < nXOuterIterations &&
2581 6 : nLastXEndIncluded < nMaxTileX;
2582 : ++iXOuterIter)
2583 : {
2584 6 : const int iXStart = nLastXEndIncluded + 1;
2585 : const int iXEndIncluded =
2586 6 : iXOuterIter + 1 == nXOuterIterations
2587 6 : ? nMaxTileX
2588 : : std::max(iXStart,
2589 6 : static_cast<int>(std::floor(
2590 0 : nMinTileX +
2591 0 : (iXOuterIter + 1) * dfTilesXPerJob - 1)));
2592 :
2593 6 : nLastXEndIncluded = iXEndIncluded;
2594 :
2595 6 : anRemainingTilesForProcess.push_back(
2596 0 : static_cast<uint64_t>(iYEndIncluded - iYStart + 1) *
2597 6 : (iXEndIncluded - iXStart + 1));
2598 :
2599 6 : CPLStringList aosArgv;
2600 6 : aosArgv.push_back(m_osGDALPath.c_str());
2601 6 : aosArgv.push_back("raster");
2602 6 : aosArgv.push_back("tile");
2603 6 : aosArgv.push_back("--config-options-in-stdin");
2604 6 : aosArgv.push_back("--config");
2605 6 : aosArgv.push_back("GDAL_NUM_THREADS=1");
2606 6 : aosArgv.push_back("--config");
2607 6 : aosArgv.push_back(
2608 : CPLSPrintf("GDAL_CACHEMAX=%" PRIu64, nCacheMaxPerProcess));
2609 6 : aosArgv.push_back("--num-threads");
2610 6 : aosArgv.push_back("1");
2611 6 : aosArgv.push_back("--min-x");
2612 6 : aosArgv.push_back(CPLSPrintf("%d", iXStart));
2613 6 : aosArgv.push_back("--max-x");
2614 6 : aosArgv.push_back(CPLSPrintf("%d", iXEndIncluded));
2615 6 : aosArgv.push_back("--min-y");
2616 6 : aosArgv.push_back(CPLSPrintf("%d", iYStart));
2617 6 : aosArgv.push_back("--max-y");
2618 6 : aosArgv.push_back(CPLSPrintf("%d", iYEndIncluded));
2619 6 : aosArgv.push_back("--webviewer");
2620 6 : aosArgv.push_back("none");
2621 6 : aosArgv.push_back("--progress-forked");
2622 6 : aosArgv.push_back("--input");
2623 6 : aosArgv.push_back(m_poSrcDS->GetDescription());
2624 300 : for (const auto &arg : GetArgs())
2625 : {
2626 362 : if (arg->IsExplicitlySet() && arg->GetName() != "min-x" &&
2627 102 : arg->GetName() != "min-y" && arg->GetName() != "max-x" &&
2628 98 : arg->GetName() != "max-y" && arg->GetName() != "min-zoom" &&
2629 60 : arg->GetName() != "progress" &&
2630 60 : arg->GetName() != "progress-forked" &&
2631 54 : arg->GetName() != "input" &&
2632 46 : arg->GetName() != "num-threads" &&
2633 350 : arg->GetName() != "webviewer" &&
2634 22 : arg->GetName() != "parallel-method")
2635 : {
2636 16 : if (!AddArgToArgv(arg.get(), aosArgv))
2637 0 : return false;
2638 : }
2639 : }
2640 :
2641 6 : std::string cmdLine;
2642 176 : for (const char *arg : aosArgv)
2643 : {
2644 170 : if (!cmdLine.empty())
2645 164 : cmdLine += ' ';
2646 170 : cmdLine += arg;
2647 : }
2648 6 : CPLDebugOnly("gdal_raster_tile", "Spawning %s", cmdLine.c_str());
2649 6 : asCommandLines.push_back(std::move(cmdLine));
2650 :
2651 : CPLSpawnedProcess *hSpawnedProcess =
2652 6 : CPLSpawnAsync(nullptr, aosArgv.List(),
2653 : /* bCreateInputPipe = */ true,
2654 : /* bCreateOutputPipe = */ true,
2655 6 : /* bCreateErrorPipe = */ true, nullptr);
2656 6 : if (!hSpawnedProcess)
2657 : {
2658 0 : ReportError(CE_Failure, CPLE_AppDefined,
2659 : "Spawning child gdal process '%s' failed",
2660 0 : asCommandLines.back().c_str());
2661 0 : bRet = false;
2662 0 : break;
2663 : }
2664 :
2665 6 : CPLDebugOnly("gdal_raster_tile",
2666 : "Job for y in [%d,%d] and x in [%d,%d], "
2667 : "run by process %" PRIu64,
2668 : iYStart, iYEndIncluded, iXStart, iXEndIncluded,
2669 : static_cast<uint64_t>(
2670 : CPLSpawnAsyncGetChildProcessId(hSpawnedProcess)));
2671 :
2672 6 : ahSpawnedProcesses.push_back(hSpawnedProcess);
2673 :
2674 6 : SendConfigOptions(hSpawnedProcess, bRet);
2675 6 : if (!bRet)
2676 : {
2677 0 : ReportError(CE_Failure, CPLE_AppDefined,
2678 : "Could not transmit config options to child gdal "
2679 : "process '%s'",
2680 0 : asCommandLines.back().c_str());
2681 0 : break;
2682 : }
2683 : }
2684 : }
2685 :
2686 2 : uint64_t nCurTile = 0;
2687 2 : GetProgressForChildProcesses(bRet, ahSpawnedProcesses,
2688 : anRemainingTilesForProcess, nCurTile,
2689 : nTotalTiles, pfnProgress, pProgressData);
2690 :
2691 2 : WaitForSpawnedProcesses(bRet, asCommandLines, ahSpawnedProcesses);
2692 :
2693 2 : if (bRet && nCurTile != nBaseTiles)
2694 : {
2695 0 : bRet = false;
2696 0 : ReportError(CE_Failure, CPLE_AppDefined,
2697 : "Not all tiles at max zoom level have been "
2698 : "generated. Got %" PRIu64 ", expected %" PRIu64,
2699 : nCurTile, nBaseTiles);
2700 : }
2701 :
2702 2 : return bRet;
2703 : }
2704 :
2705 : /************************************************************************/
2706 : /* GDALRasterTileAlgorithm::GenerateOverviewTilesSpawnMethod() */
2707 : /************************************************************************/
2708 :
2709 2 : bool GDALRasterTileAlgorithm::GenerateOverviewTilesSpawnMethod(
2710 : int iZ, int nOvrMinTileX, int nOvrMinTileY, int nOvrMaxTileX,
2711 : int nOvrMaxTileY, std::atomic<uint64_t> &nCurTile, uint64_t nTotalTiles,
2712 : GDALProgressFunc pfnProgress, void *pProgressData)
2713 : {
2714 2 : CPLAssert(!m_osGDALPath.empty());
2715 :
2716 2 : const int nOvrTilesPerCol = nOvrMaxTileY - nOvrMinTileY + 1;
2717 2 : const int nOvrTilesPerRow = nOvrMaxTileX - nOvrMinTileX + 1;
2718 2 : const uint64_t nExpectedOvrTileCount =
2719 2 : static_cast<uint64_t>(nOvrTilesPerCol) * nOvrTilesPerRow;
2720 :
2721 2 : const int nMaxJobCount = GetMaxChildCount(
2722 0 : std::max(1, static_cast<int>(std::min<uint64_t>(
2723 0 : m_numThreads, nExpectedOvrTileCount /
2724 2 : GetThresholdMinTilesPerJob()))));
2725 :
2726 : double dfTilesYPerJob;
2727 : int nYOuterIterations;
2728 : double dfTilesXPerJob;
2729 : int nXOuterIterations;
2730 2 : ComputeJobChunkSize(nMaxJobCount, nOvrTilesPerCol, nOvrTilesPerRow,
2731 : dfTilesYPerJob, nYOuterIterations, dfTilesXPerJob,
2732 : nXOuterIterations);
2733 :
2734 2 : CPLDebugOnly("gdal_raster_tile",
2735 : "z=%d, nYOuterIterations=%d, dfTilesYPerJob=%g, "
2736 : "nXOuterIterations=%d, dfTilesXPerJob=%g",
2737 : iZ, nYOuterIterations, dfTilesYPerJob, nXOuterIterations,
2738 : dfTilesXPerJob);
2739 :
2740 4 : std::vector<std::string> asCommandLines;
2741 4 : std::vector<CPLSpawnedProcess *> ahSpawnedProcesses;
2742 4 : std::vector<uint64_t> anRemainingTilesForProcess;
2743 :
2744 2 : const uint64_t nCacheMaxPerProcess = GDALGetCacheMax64() / nMaxJobCount;
2745 :
2746 2 : int nLastYEndIncluded = nOvrMinTileY - 1;
2747 2 : bool bRet = true;
2748 8 : for (int iYOuterIter = 0; bRet && iYOuterIter < nYOuterIterations &&
2749 6 : nLastYEndIncluded < nOvrMaxTileY;
2750 : ++iYOuterIter)
2751 : {
2752 6 : const int iYStart = nLastYEndIncluded + 1;
2753 : const int iYEndIncluded =
2754 6 : iYOuterIter + 1 == nYOuterIterations
2755 10 : ? nOvrMaxTileY
2756 : : std::max(iYStart,
2757 10 : static_cast<int>(std::floor(
2758 4 : nOvrMinTileY +
2759 4 : (iYOuterIter + 1) * dfTilesYPerJob - 1)));
2760 :
2761 6 : nLastYEndIncluded = iYEndIncluded;
2762 :
2763 6 : int nLastXEndIncluded = nOvrMinTileX - 1;
2764 14 : for (int iXOuterIter = 0; bRet && iXOuterIter < nXOuterIterations &&
2765 8 : nLastXEndIncluded < nOvrMaxTileX;
2766 : ++iXOuterIter)
2767 : {
2768 8 : const int iXStart = nLastXEndIncluded + 1;
2769 : const int iXEndIncluded =
2770 8 : iXOuterIter + 1 == nXOuterIterations
2771 10 : ? nOvrMaxTileX
2772 : : std::max(iXStart,
2773 10 : static_cast<int>(std::floor(
2774 2 : nOvrMinTileX +
2775 2 : (iXOuterIter + 1) * dfTilesXPerJob - 1)));
2776 :
2777 8 : nLastXEndIncluded = iXEndIncluded;
2778 :
2779 8 : anRemainingTilesForProcess.push_back(
2780 0 : static_cast<uint64_t>(iYEndIncluded - iYStart + 1) *
2781 8 : (iXEndIncluded - iXStart + 1));
2782 :
2783 8 : CPLStringList aosArgv;
2784 8 : aosArgv.push_back(m_osGDALPath.c_str());
2785 8 : aosArgv.push_back("raster");
2786 8 : aosArgv.push_back("tile");
2787 8 : aosArgv.push_back("--config-options-in-stdin");
2788 8 : aosArgv.push_back("--config");
2789 8 : aosArgv.push_back("GDAL_NUM_THREADS=1");
2790 8 : aosArgv.push_back("--config");
2791 8 : aosArgv.push_back(
2792 : CPLSPrintf("GDAL_CACHEMAX=%" PRIu64, nCacheMaxPerProcess));
2793 8 : aosArgv.push_back("--num-threads");
2794 8 : aosArgv.push_back("1");
2795 8 : aosArgv.push_back("--ovr-zoom-level");
2796 8 : aosArgv.push_back(CPLSPrintf("%d", iZ));
2797 8 : aosArgv.push_back("--ovr-min-x");
2798 8 : aosArgv.push_back(CPLSPrintf("%d", iXStart));
2799 8 : aosArgv.push_back("--ovr-max-x");
2800 8 : aosArgv.push_back(CPLSPrintf("%d", iXEndIncluded));
2801 8 : aosArgv.push_back("--ovr-min-y");
2802 8 : aosArgv.push_back(CPLSPrintf("%d", iYStart));
2803 8 : aosArgv.push_back("--ovr-max-y");
2804 8 : aosArgv.push_back(CPLSPrintf("%d", iYEndIncluded));
2805 8 : aosArgv.push_back("--webviewer");
2806 8 : aosArgv.push_back("none");
2807 8 : aosArgv.push_back("--progress-forked");
2808 8 : aosArgv.push_back("--input");
2809 8 : aosArgv.push_back(m_dataset.GetName().c_str());
2810 400 : for (const auto &arg : GetArgs())
2811 : {
2812 488 : if (arg->IsExplicitlySet() && arg->GetName() != "progress" &&
2813 96 : arg->GetName() != "progress-forked" &&
2814 88 : arg->GetName() != "input" &&
2815 80 : arg->GetName() != "num-threads" &&
2816 480 : arg->GetName() != "webviewer" &&
2817 40 : arg->GetName() != "parallel-method")
2818 : {
2819 32 : if (!AddArgToArgv(arg.get(), aosArgv))
2820 0 : return false;
2821 : }
2822 : }
2823 :
2824 8 : std::string cmdLine;
2825 272 : for (const char *arg : aosArgv)
2826 : {
2827 264 : if (!cmdLine.empty())
2828 256 : cmdLine += ' ';
2829 264 : cmdLine += arg;
2830 : }
2831 8 : CPLDebugOnly("gdal_raster_tile", "Spawning %s", cmdLine.c_str());
2832 8 : asCommandLines.push_back(std::move(cmdLine));
2833 :
2834 : CPLSpawnedProcess *hSpawnedProcess =
2835 8 : CPLSpawnAsync(nullptr, aosArgv.List(),
2836 : /* bCreateInputPipe = */ true,
2837 : /* bCreateOutputPipe = */ true,
2838 8 : /* bCreateErrorPipe = */ true, nullptr);
2839 8 : if (!hSpawnedProcess)
2840 : {
2841 0 : ReportError(CE_Failure, CPLE_AppDefined,
2842 : "Spawning child gdal process '%s' failed",
2843 0 : asCommandLines.back().c_str());
2844 0 : bRet = false;
2845 0 : break;
2846 : }
2847 :
2848 8 : CPLDebugOnly("gdal_raster_tile",
2849 : "Job for z = %d, y in [%d,%d] and x in [%d,%d], "
2850 : "run by process %" PRIu64,
2851 : iZ, iYStart, iYEndIncluded, iXStart, iXEndIncluded,
2852 : static_cast<uint64_t>(
2853 : CPLSpawnAsyncGetChildProcessId(hSpawnedProcess)));
2854 :
2855 8 : ahSpawnedProcesses.push_back(hSpawnedProcess);
2856 :
2857 8 : SendConfigOptions(hSpawnedProcess, bRet);
2858 8 : if (!bRet)
2859 : {
2860 0 : ReportError(CE_Failure, CPLE_AppDefined,
2861 : "Could not transmit config options to child gdal "
2862 : "process '%s'",
2863 0 : asCommandLines.back().c_str());
2864 0 : break;
2865 : }
2866 : }
2867 : }
2868 :
2869 2 : uint64_t nCurTileLocal = nCurTile;
2870 2 : GetProgressForChildProcesses(bRet, ahSpawnedProcesses,
2871 : anRemainingTilesForProcess, nCurTileLocal,
2872 : nTotalTiles, pfnProgress, pProgressData);
2873 :
2874 2 : WaitForSpawnedProcesses(bRet, asCommandLines, ahSpawnedProcesses);
2875 :
2876 2 : if (bRet && nCurTileLocal - nCurTile != nExpectedOvrTileCount)
2877 : {
2878 0 : bRet = false;
2879 0 : ReportError(CE_Failure, CPLE_AppDefined,
2880 : "Not all tiles at zoom level %d have been "
2881 : "generated. Got %" PRIu64 ", expected %" PRIu64,
2882 0 : iZ, nCurTileLocal - nCurTile, nExpectedOvrTileCount);
2883 : }
2884 :
2885 2 : nCurTile = nCurTileLocal;
2886 :
2887 2 : return bRet;
2888 : }
2889 :
2890 : /************************************************************************/
2891 : /* GDALRasterTileAlgorithm::RunImpl() */
2892 : /************************************************************************/
2893 :
2894 107 : bool GDALRasterTileAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
2895 : void *pProgressData)
2896 : {
2897 107 : m_poSrcDS = m_dataset.GetDatasetRef();
2898 107 : CPLAssert(m_poSrcDS);
2899 107 : const int nSrcWidth = m_poSrcDS->GetRasterXSize();
2900 107 : const int nSrcHeight = m_poSrcDS->GetRasterYSize();
2901 107 : if (m_poSrcDS->GetRasterCount() == 0 || nSrcWidth == 0 || nSrcHeight == 0)
2902 : {
2903 1 : ReportError(CE_Failure, CPLE_AppDefined, "Invalid source dataset");
2904 1 : return false;
2905 : }
2906 :
2907 106 : auto poSrcDriver = m_poSrcDS->GetDriver();
2908 106 : m_bIsNamedNonMemSrcDS =
2909 157 : (m_poSrcDS->GetDescription()[0] != 0 && poSrcDriver &&
2910 51 : !EQUAL(poSrcDriver->GetDescription(), "MEM"));
2911 :
2912 106 : if (m_parallelMethod == "spawn")
2913 : {
2914 5 : const char *pszErrorMsg = "";
2915 5 : if (!IsCompatibleOfSpawn(pszErrorMsg))
2916 : {
2917 3 : if (pszErrorMsg[0])
2918 2 : ReportError(CE_Failure, CPLE_AppDefined, "%s", pszErrorMsg);
2919 3 : return false;
2920 : }
2921 : }
2922 :
2923 103 : if (m_resampling == "near")
2924 1 : m_resampling = "nearest";
2925 103 : if (m_overviewResampling == "near")
2926 1 : m_overviewResampling = "nearest";
2927 102 : else if (m_overviewResampling.empty())
2928 97 : m_overviewResampling = m_resampling;
2929 :
2930 206 : CPLStringList aosWarpOptions;
2931 103 : if (!m_excludedValues.empty() || m_nodataValuesPctThreshold < 100)
2932 : {
2933 : aosWarpOptions.SetNameValue(
2934 : "NODATA_VALUES_PCT_THRESHOLD",
2935 3 : CPLSPrintf("%g", m_nodataValuesPctThreshold));
2936 3 : if (!m_excludedValues.empty())
2937 : {
2938 : aosWarpOptions.SetNameValue("EXCLUDED_VALUES",
2939 1 : m_excludedValues.c_str());
2940 : aosWarpOptions.SetNameValue(
2941 : "EXCLUDED_VALUES_PCT_THRESHOLD",
2942 1 : CPLSPrintf("%g", m_excludedValuesPctThreshold));
2943 : }
2944 : }
2945 :
2946 103 : if (m_poSrcDS->GetRasterBand(1)->GetColorInterpretation() ==
2947 107 : GCI_PaletteIndex &&
2948 4 : ((m_resampling != "nearest" && m_resampling != "mode") ||
2949 1 : (m_overviewResampling != "nearest" && m_overviewResampling != "mode")))
2950 : {
2951 1 : ReportError(CE_Failure, CPLE_NotSupported,
2952 : "Datasets with color table not supported with non-nearest "
2953 : "or non-mode resampling. Run 'gdal raster "
2954 : "color-map' before or set the 'resampling' argument to "
2955 : "'nearest' or 'mode'.");
2956 1 : return false;
2957 : }
2958 :
2959 102 : const auto eSrcDT = m_poSrcDS->GetRasterBand(1)->GetRasterDataType();
2960 102 : m_poDstDriver =
2961 102 : GetGDALDriverManager()->GetDriverByName(m_outputFormat.c_str());
2962 102 : if (!m_poDstDriver)
2963 : {
2964 1 : ReportError(CE_Failure, CPLE_AppDefined,
2965 : "Invalid value for argument 'output-format'. Driver '%s' "
2966 : "does not exist",
2967 : m_outputFormat.c_str());
2968 1 : return false;
2969 : }
2970 :
2971 101 : if (!ValidateOutputFormat(eSrcDT))
2972 8 : return false;
2973 :
2974 : const char *pszExtensions =
2975 93 : m_poDstDriver->GetMetadataItem(GDAL_DMD_EXTENSIONS);
2976 93 : CPLAssert(pszExtensions && pszExtensions[0] != 0);
2977 : const CPLStringList aosExtensions(
2978 186 : CSLTokenizeString2(pszExtensions, " ", 0));
2979 93 : const char *pszExtension = aosExtensions[0];
2980 93 : GDALGeoTransform srcGT;
2981 93 : const bool bHasSrcGT = m_poSrcDS->GetGeoTransform(srcGT) == CE_None;
2982 : const bool bHasNorthUpSrcGT =
2983 93 : bHasSrcGT && srcGT[2] == 0 && srcGT[4] == 0 && srcGT[5] < 0;
2984 186 : OGRSpatialReference oSRS_TMS;
2985 :
2986 93 : if (m_tilingScheme == "raster")
2987 : {
2988 4 : if (const auto poSRS = m_poSrcDS->GetSpatialRef())
2989 3 : oSRS_TMS = *poSRS;
2990 : }
2991 : else
2992 : {
2993 1 : if (!bHasSrcGT && m_poSrcDS->GetGCPCount() == 0 &&
2994 91 : m_poSrcDS->GetMetadata("GEOLOCATION") == nullptr &&
2995 1 : m_poSrcDS->GetMetadata("RPC") == nullptr)
2996 : {
2997 1 : ReportError(CE_Failure, CPLE_NotSupported,
2998 : "Ungeoreferenced datasets are not supported, unless "
2999 : "'tiling-scheme' is set to 'raster'");
3000 1 : return false;
3001 : }
3002 :
3003 88 : if (m_poSrcDS->GetMetadata("GEOLOCATION") == nullptr &&
3004 88 : m_poSrcDS->GetMetadata("RPC") == nullptr &&
3005 177 : m_poSrcDS->GetSpatialRef() == nullptr &&
3006 1 : m_poSrcDS->GetGCPSpatialRef() == nullptr)
3007 : {
3008 1 : ReportError(CE_Failure, CPLE_NotSupported,
3009 : "Ungeoreferenced datasets are not supported, unless "
3010 : "'tiling-scheme' is set to 'raster'");
3011 1 : return false;
3012 : }
3013 : }
3014 :
3015 91 : if (m_copySrcMetadata)
3016 : {
3017 8 : CPLStringList aosMD(CSLDuplicate(m_poSrcDS->GetMetadata()));
3018 4 : const CPLStringList aosNewMD(m_metadata);
3019 8 : for (const auto [key, value] : cpl::IterateNameValue(aosNewMD))
3020 : {
3021 4 : aosMD.SetNameValue(key, value);
3022 : }
3023 4 : m_metadata = aosMD;
3024 : }
3025 :
3026 91 : GDALGeoTransform srcGTModif{0, 1, 0, 0, 0, -1};
3027 :
3028 91 : if (m_tilingScheme == "mercator")
3029 1 : m_tilingScheme = "WebMercatorQuad";
3030 90 : else if (m_tilingScheme == "geodetic")
3031 1 : m_tilingScheme = "WorldCRS84Quad";
3032 89 : else if (m_tilingScheme == "raster")
3033 : {
3034 4 : if (m_tileSize == 0)
3035 4 : m_tileSize = 256;
3036 4 : if (m_maxZoomLevel < 0)
3037 : {
3038 3 : m_maxZoomLevel = static_cast<int>(std::ceil(std::log2(
3039 6 : std::max(1, std::max(nSrcWidth, nSrcHeight) / m_tileSize))));
3040 : }
3041 4 : if (bHasNorthUpSrcGT)
3042 : {
3043 3 : srcGTModif = srcGT;
3044 : }
3045 : }
3046 :
3047 : auto poTMS =
3048 91 : m_tilingScheme == "raster"
3049 : ? gdal::TileMatrixSet::createRaster(
3050 4 : nSrcWidth, nSrcHeight, m_tileSize, 1 + m_maxZoomLevel,
3051 16 : srcGTModif[0], srcGTModif[3], srcGTModif[1], -srcGTModif[5],
3052 95 : oSRS_TMS.IsEmpty() ? std::string() : oSRS_TMS.exportToWkt())
3053 : : gdal::TileMatrixSet::parse(
3054 186 : m_mapTileMatrixIdentifierToScheme[m_tilingScheme].c_str());
3055 : // Enforced by SetChoices() on the m_tilingScheme argument
3056 91 : CPLAssert(poTMS && !poTMS->hasVariableMatrixWidth());
3057 :
3058 182 : CPLStringList aosTO;
3059 91 : if (m_tilingScheme == "raster")
3060 : {
3061 4 : aosTO.SetNameValue("SRC_METHOD", "GEOTRANSFORM");
3062 : }
3063 : else
3064 : {
3065 87 : CPL_IGNORE_RET_VAL(oSRS_TMS.SetFromUserInput(poTMS->crs().c_str()));
3066 87 : aosTO.SetNameValue("DST_SRS", oSRS_TMS.exportToWkt().c_str());
3067 : }
3068 :
3069 91 : const char *pszAuthName = oSRS_TMS.GetAuthorityName(nullptr);
3070 91 : const char *pszAuthCode = oSRS_TMS.GetAuthorityCode(nullptr);
3071 91 : const int nEPSGCode =
3072 90 : (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
3073 181 : ? atoi(pszAuthCode)
3074 : : 0;
3075 :
3076 : const bool bInvertAxisTMS =
3077 178 : m_tilingScheme != "raster" &&
3078 87 : (oSRS_TMS.EPSGTreatsAsLatLong() != FALSE ||
3079 87 : oSRS_TMS.EPSGTreatsAsNorthingEasting() != FALSE);
3080 :
3081 91 : oSRS_TMS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3082 :
3083 : std::unique_ptr<void, decltype(&GDALDestroyTransformer)> hTransformArg(
3084 182 : nullptr, GDALDestroyTransformer);
3085 :
3086 : // Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
3087 : // ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
3088 : // EPSG:3857.
3089 91 : std::unique_ptr<GDALDataset> poTmpDS;
3090 91 : bool bEPSG3857Adjust = false;
3091 91 : if (nEPSGCode == 3857 && bHasNorthUpSrcGT)
3092 : {
3093 65 : const auto poSrcSRS = m_poSrcDS->GetSpatialRef();
3094 65 : if (poSrcSRS && poSrcSRS->IsGeographic())
3095 : {
3096 34 : double maxLat = srcGT[3];
3097 34 : double minLat = srcGT[3] + nSrcHeight * srcGT[5];
3098 : // Corresponds to the latitude of below MAX_GM
3099 34 : constexpr double MAX_LAT = 85.0511287798066;
3100 34 : bool bModified = false;
3101 34 : if (maxLat > MAX_LAT)
3102 : {
3103 33 : maxLat = MAX_LAT;
3104 33 : bModified = true;
3105 : }
3106 34 : if (minLat < -MAX_LAT)
3107 : {
3108 33 : minLat = -MAX_LAT;
3109 33 : bModified = true;
3110 : }
3111 34 : if (bModified)
3112 : {
3113 66 : CPLStringList aosOptions;
3114 33 : aosOptions.AddString("-of");
3115 33 : aosOptions.AddString("VRT");
3116 33 : aosOptions.AddString("-projwin");
3117 33 : aosOptions.AddString(CPLSPrintf("%.17g", srcGT[0]));
3118 33 : aosOptions.AddString(CPLSPrintf("%.17g", maxLat));
3119 : aosOptions.AddString(
3120 33 : CPLSPrintf("%.17g", srcGT[0] + nSrcWidth * srcGT[1]));
3121 33 : aosOptions.AddString(CPLSPrintf("%.17g", minLat));
3122 : auto psOptions =
3123 33 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
3124 33 : poTmpDS.reset(GDALDataset::FromHandle(GDALTranslate(
3125 : "", GDALDataset::ToHandle(m_poSrcDS), psOptions, nullptr)));
3126 33 : GDALTranslateOptionsFree(psOptions);
3127 33 : if (poTmpDS)
3128 : {
3129 33 : bEPSG3857Adjust = true;
3130 33 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
3131 33 : GDALDataset::FromHandle(poTmpDS.get()), nullptr,
3132 33 : aosTO.List()));
3133 : }
3134 : }
3135 : }
3136 : }
3137 :
3138 91 : GDALGeoTransform dstGT;
3139 : double adfExtent[4];
3140 : int nXSize, nYSize;
3141 :
3142 : bool bSuggestOK;
3143 91 : if (m_tilingScheme == "raster")
3144 : {
3145 4 : bSuggestOK = true;
3146 4 : nXSize = nSrcWidth;
3147 4 : nYSize = nSrcHeight;
3148 4 : dstGT = srcGTModif;
3149 4 : adfExtent[0] = dstGT[0];
3150 4 : adfExtent[1] = dstGT[3] + nSrcHeight * dstGT[5];
3151 4 : adfExtent[2] = dstGT[0] + nSrcWidth * dstGT[1];
3152 4 : adfExtent[3] = dstGT[3];
3153 : }
3154 : else
3155 : {
3156 87 : if (!hTransformArg)
3157 : {
3158 54 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
3159 54 : m_poSrcDS, nullptr, aosTO.List()));
3160 : }
3161 87 : if (!hTransformArg)
3162 : {
3163 1 : return false;
3164 : }
3165 86 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3166 86 : bSuggestOK =
3167 172 : (GDALSuggestedWarpOutput2(
3168 86 : m_poSrcDS,
3169 86 : static_cast<GDALTransformerInfo *>(hTransformArg.get())
3170 : ->pfnTransform,
3171 : hTransformArg.get(), dstGT.data(), &nXSize, &nYSize, adfExtent,
3172 : 0) == CE_None);
3173 : }
3174 90 : if (!bSuggestOK)
3175 : {
3176 1 : ReportError(CE_Failure, CPLE_AppDefined,
3177 : "Cannot determine extent of raster in target CRS");
3178 1 : return false;
3179 : }
3180 :
3181 89 : poTmpDS.reset();
3182 :
3183 89 : if (bEPSG3857Adjust)
3184 : {
3185 33 : constexpr double SPHERICAL_RADIUS = 6378137.0;
3186 33 : constexpr double MAX_GM =
3187 : SPHERICAL_RADIUS * M_PI; // 20037508.342789244
3188 33 : double maxNorthing = dstGT[3];
3189 33 : double minNorthing = dstGT[3] + dstGT[5] * nYSize;
3190 33 : bool bChanged = false;
3191 33 : if (maxNorthing > MAX_GM)
3192 : {
3193 30 : bChanged = true;
3194 30 : maxNorthing = MAX_GM;
3195 : }
3196 33 : if (minNorthing < -MAX_GM)
3197 : {
3198 30 : bChanged = true;
3199 30 : minNorthing = -MAX_GM;
3200 : }
3201 33 : if (bChanged)
3202 : {
3203 30 : dstGT[3] = maxNorthing;
3204 30 : nYSize = int((maxNorthing - minNorthing) / (-dstGT[5]) + 0.5);
3205 30 : adfExtent[1] = maxNorthing + nYSize * dstGT[5];
3206 30 : adfExtent[3] = maxNorthing;
3207 : }
3208 : }
3209 :
3210 89 : const auto &tileMatrixList = poTMS->tileMatrixList();
3211 89 : if (m_maxZoomLevel >= 0)
3212 : {
3213 40 : if (m_maxZoomLevel >= static_cast<int>(tileMatrixList.size()))
3214 : {
3215 1 : ReportError(CE_Failure, CPLE_AppDefined,
3216 : "max-zoom = %d is invalid. It must be in [0,%d] range",
3217 : m_maxZoomLevel,
3218 1 : static_cast<int>(tileMatrixList.size()) - 1);
3219 1 : return false;
3220 : }
3221 : }
3222 : else
3223 : {
3224 49 : const double dfComputedRes = dstGT[1];
3225 49 : double dfPrevRes = 0.0;
3226 49 : double dfRes = 0.0;
3227 49 : constexpr double EPSILON = 1e-8;
3228 :
3229 49 : if (m_minZoomLevel >= 0)
3230 17 : m_maxZoomLevel = m_minZoomLevel;
3231 : else
3232 32 : m_maxZoomLevel = 0;
3233 :
3234 164 : for (; m_maxZoomLevel < static_cast<int>(tileMatrixList.size());
3235 115 : m_maxZoomLevel++)
3236 : {
3237 163 : dfRes = tileMatrixList[m_maxZoomLevel].mResX;
3238 163 : if (dfComputedRes > dfRes ||
3239 115 : fabs(dfComputedRes - dfRes) / dfRes <= EPSILON)
3240 : break;
3241 115 : dfPrevRes = dfRes;
3242 : }
3243 49 : if (m_maxZoomLevel >= static_cast<int>(tileMatrixList.size()))
3244 : {
3245 1 : ReportError(CE_Failure, CPLE_AppDefined,
3246 : "Could not find an appropriate zoom level. Perhaps "
3247 : "min-zoom is too large?");
3248 1 : return false;
3249 : }
3250 :
3251 48 : if (m_maxZoomLevel > 0 && fabs(dfComputedRes - dfRes) / dfRes > EPSILON)
3252 : {
3253 : // Round to closest resolution
3254 43 : if (dfPrevRes / dfComputedRes < dfComputedRes / dfRes)
3255 23 : m_maxZoomLevel--;
3256 : }
3257 : }
3258 87 : if (m_minZoomLevel < 0)
3259 47 : m_minZoomLevel = m_maxZoomLevel;
3260 :
3261 174 : auto tileMatrix = tileMatrixList[m_maxZoomLevel];
3262 87 : int nMinTileX = 0;
3263 87 : int nMinTileY = 0;
3264 87 : int nMaxTileX = 0;
3265 87 : int nMaxTileY = 0;
3266 87 : bool bIntersects = false;
3267 87 : if (!GetTileIndices(tileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
3268 : nMinTileX, nMinTileY, nMaxTileX, nMaxTileY,
3269 87 : m_noIntersectionIsOK, bIntersects,
3270 : /* checkRasterOverflow = */ false))
3271 : {
3272 1 : return false;
3273 : }
3274 86 : if (!bIntersects)
3275 1 : return true;
3276 :
3277 : // Potentially restrict tiling to user specified coordinates
3278 85 : if (m_minTileX >= tileMatrix.mMatrixWidth)
3279 : {
3280 1 : ReportError(CE_Failure, CPLE_IllegalArg,
3281 : "'min-x' value must be in [0,%d] range",
3282 1 : tileMatrix.mMatrixWidth - 1);
3283 1 : return false;
3284 : }
3285 84 : if (m_maxTileX >= tileMatrix.mMatrixWidth)
3286 : {
3287 1 : ReportError(CE_Failure, CPLE_IllegalArg,
3288 : "'max-x' value must be in [0,%d] range",
3289 1 : tileMatrix.mMatrixWidth - 1);
3290 1 : return false;
3291 : }
3292 83 : if (m_minTileY >= tileMatrix.mMatrixHeight)
3293 : {
3294 1 : ReportError(CE_Failure, CPLE_IllegalArg,
3295 : "'min-y' value must be in [0,%d] range",
3296 1 : tileMatrix.mMatrixHeight - 1);
3297 1 : return false;
3298 : }
3299 82 : if (m_maxTileY >= tileMatrix.mMatrixHeight)
3300 : {
3301 1 : ReportError(CE_Failure, CPLE_IllegalArg,
3302 : "'max-y' value must be in [0,%d] range",
3303 1 : tileMatrix.mMatrixHeight - 1);
3304 1 : return false;
3305 : }
3306 :
3307 81 : if ((m_minTileX >= 0 && m_minTileX > nMaxTileX) ||
3308 79 : (m_minTileY >= 0 && m_minTileY > nMaxTileY) ||
3309 79 : (m_maxTileX >= 0 && m_maxTileX < nMinTileX) ||
3310 79 : (m_maxTileY >= 0 && m_maxTileY < nMinTileY))
3311 : {
3312 2 : ReportError(
3313 2 : m_noIntersectionIsOK ? CE_Warning : CE_Failure, CPLE_AppDefined,
3314 : "Dataset extent not intersecting specified min/max X/Y tile "
3315 : "coordinates");
3316 2 : return m_noIntersectionIsOK;
3317 : }
3318 79 : if (m_minTileX >= 0 && m_minTileX > nMinTileX)
3319 : {
3320 2 : nMinTileX = m_minTileX;
3321 2 : adfExtent[0] = tileMatrix.mTopLeftX +
3322 2 : nMinTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
3323 : }
3324 79 : if (m_minTileY >= 0 && m_minTileY > nMinTileY)
3325 : {
3326 5 : nMinTileY = m_minTileY;
3327 5 : adfExtent[3] = tileMatrix.mTopLeftY -
3328 5 : nMinTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
3329 : }
3330 79 : if (m_maxTileX >= 0 && m_maxTileX < nMaxTileX)
3331 : {
3332 2 : nMaxTileX = m_maxTileX;
3333 2 : adfExtent[2] = tileMatrix.mTopLeftX + (nMaxTileX + 1) *
3334 2 : tileMatrix.mResX *
3335 2 : tileMatrix.mTileWidth;
3336 : }
3337 79 : if (m_maxTileY >= 0 && m_maxTileY < nMaxTileY)
3338 : {
3339 6 : nMaxTileY = m_maxTileY;
3340 6 : adfExtent[1] = tileMatrix.mTopLeftY - (nMaxTileY + 1) *
3341 6 : tileMatrix.mResY *
3342 6 : tileMatrix.mTileHeight;
3343 : }
3344 :
3345 79 : if (nMaxTileX - nMinTileX + 1 > INT_MAX / tileMatrix.mTileWidth ||
3346 78 : nMaxTileY - nMinTileY + 1 > INT_MAX / tileMatrix.mTileHeight)
3347 : {
3348 1 : ReportError(CE_Failure, CPLE_AppDefined, "Too large zoom level");
3349 1 : return false;
3350 : }
3351 :
3352 156 : dstGT[0] = tileMatrix.mTopLeftX +
3353 78 : nMinTileX * tileMatrix.mResX * tileMatrix.mTileWidth;
3354 78 : dstGT[1] = tileMatrix.mResX;
3355 78 : dstGT[2] = 0;
3356 156 : dstGT[3] = tileMatrix.mTopLeftY -
3357 78 : nMinTileY * tileMatrix.mResY * tileMatrix.mTileHeight;
3358 78 : dstGT[4] = 0;
3359 78 : dstGT[5] = -tileMatrix.mResY;
3360 :
3361 : /* -------------------------------------------------------------------- */
3362 : /* Setup warp options. */
3363 : /* -------------------------------------------------------------------- */
3364 : std::unique_ptr<GDALWarpOptions, decltype(&GDALDestroyWarpOptions)> psWO(
3365 156 : GDALCreateWarpOptions(), GDALDestroyWarpOptions);
3366 :
3367 78 : psWO->papszWarpOptions = CSLSetNameValue(nullptr, "OPTIMIZE_SIZE", "YES");
3368 156 : psWO->papszWarpOptions =
3369 78 : CSLSetNameValue(psWO->papszWarpOptions, "SAMPLE_GRID", "YES");
3370 156 : psWO->papszWarpOptions =
3371 78 : CSLMerge(psWO->papszWarpOptions, aosWarpOptions.List());
3372 :
3373 78 : int bHasSrcNoData = false;
3374 : const double dfSrcNoDataValue =
3375 78 : m_poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasSrcNoData);
3376 :
3377 : const bool bLastSrcBandIsAlpha =
3378 127 : (m_poSrcDS->GetRasterCount() > 1 &&
3379 49 : m_poSrcDS->GetRasterBand(m_poSrcDS->GetRasterCount())
3380 49 : ->GetColorInterpretation() == GCI_AlphaBand);
3381 :
3382 78 : const bool bOutputSupportsAlpha = !EQUAL(m_outputFormat.c_str(), "JPEG");
3383 78 : const bool bOutputSupportsNoData = EQUAL(m_outputFormat.c_str(), "GTiff");
3384 78 : const bool bDstNoDataSpecified = GetArg("dst-nodata")->IsExplicitlySet();
3385 : auto poColorTable = std::unique_ptr<GDALColorTable>(
3386 78 : [this]()
3387 : {
3388 78 : auto poCT = m_poSrcDS->GetRasterBand(1)->GetColorTable();
3389 78 : return poCT ? poCT->Clone() : nullptr;
3390 156 : }());
3391 :
3392 78 : const bool bUserAskedForAlpha = m_addalpha;
3393 78 : if (!m_noalpha && !m_addalpha)
3394 : {
3395 88 : m_addalpha = !(bHasSrcNoData && bOutputSupportsNoData) &&
3396 88 : !bDstNoDataSpecified && poColorTable == nullptr;
3397 : }
3398 78 : m_addalpha &= bOutputSupportsAlpha;
3399 :
3400 78 : psWO->nBandCount = m_poSrcDS->GetRasterCount();
3401 78 : if (bLastSrcBandIsAlpha)
3402 : {
3403 5 : --psWO->nBandCount;
3404 5 : psWO->nSrcAlphaBand = m_poSrcDS->GetRasterCount();
3405 : }
3406 :
3407 78 : if (bHasSrcNoData)
3408 : {
3409 26 : psWO->padfSrcNoDataReal =
3410 13 : static_cast<double *>(CPLCalloc(psWO->nBandCount, sizeof(double)));
3411 32 : for (int i = 0; i < psWO->nBandCount; ++i)
3412 : {
3413 19 : psWO->padfSrcNoDataReal[i] = dfSrcNoDataValue;
3414 : }
3415 : }
3416 :
3417 78 : if ((bHasSrcNoData && !m_addalpha && bOutputSupportsNoData) ||
3418 : bDstNoDataSpecified)
3419 : {
3420 18 : psWO->padfDstNoDataReal =
3421 9 : static_cast<double *>(CPLCalloc(psWO->nBandCount, sizeof(double)));
3422 18 : for (int i = 0; i < psWO->nBandCount; ++i)
3423 : {
3424 9 : psWO->padfDstNoDataReal[i] =
3425 9 : bDstNoDataSpecified ? m_dstNoData : dfSrcNoDataValue;
3426 : }
3427 : }
3428 :
3429 78 : psWO->eWorkingDataType = eSrcDT;
3430 :
3431 78 : GDALGetWarpResampleAlg(m_resampling.c_str(), psWO->eResampleAlg);
3432 :
3433 : /* -------------------------------------------------------------------- */
3434 : /* Setup band mapping. */
3435 : /* -------------------------------------------------------------------- */
3436 :
3437 156 : psWO->panSrcBands =
3438 78 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
3439 156 : psWO->panDstBands =
3440 78 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
3441 :
3442 65787 : for (int i = 0; i < psWO->nBandCount; i++)
3443 : {
3444 65709 : psWO->panSrcBands[i] = i + 1;
3445 65709 : psWO->panDstBands[i] = i + 1;
3446 : }
3447 :
3448 78 : if (m_addalpha)
3449 65 : psWO->nDstAlphaBand = psWO->nBandCount + 1;
3450 :
3451 : const int nDstBands =
3452 78 : psWO->nDstAlphaBand ? psWO->nDstAlphaBand : psWO->nBandCount;
3453 :
3454 156 : std::vector<GByte> dstBuffer;
3455 : const uint64_t dstBufferSize =
3456 78 : static_cast<uint64_t>(tileMatrix.mTileWidth) * tileMatrix.mTileHeight *
3457 78 : nDstBands * GDALGetDataTypeSizeBytes(psWO->eWorkingDataType);
3458 : const uint64_t nUsableRAM =
3459 78 : std::min<uint64_t>(INT_MAX, CPLGetUsablePhysicalRAM() / 4);
3460 78 : if (dstBufferSize <=
3461 78 : (nUsableRAM ? nUsableRAM : static_cast<uint64_t>(INT_MAX)))
3462 : {
3463 : try
3464 : {
3465 77 : dstBuffer.resize(static_cast<size_t>(dstBufferSize));
3466 : }
3467 0 : catch (const std::exception &)
3468 : {
3469 : }
3470 : }
3471 78 : if (dstBuffer.size() < dstBufferSize)
3472 : {
3473 1 : ReportError(CE_Failure, CPLE_AppDefined,
3474 : "Tile size and/or number of bands too large compared to "
3475 : "available RAM");
3476 1 : return false;
3477 : }
3478 :
3479 : FakeMaxZoomDataset oFakeMaxZoomDS(
3480 77 : (nMaxTileX - nMinTileX + 1) * tileMatrix.mTileWidth,
3481 77 : (nMaxTileY - nMinTileY + 1) * tileMatrix.mTileHeight, nDstBands,
3482 77 : tileMatrix.mTileWidth, tileMatrix.mTileHeight, psWO->eWorkingDataType,
3483 154 : dstGT, oSRS_TMS, dstBuffer);
3484 77 : CPL_IGNORE_RET_VAL(oFakeMaxZoomDS.GetSpatialRef());
3485 :
3486 77 : psWO->hSrcDS = GDALDataset::ToHandle(m_poSrcDS);
3487 77 : psWO->hDstDS = GDALDataset::ToHandle(&oFakeMaxZoomDS);
3488 :
3489 77 : std::unique_ptr<GDALDataset> tmpSrcDS;
3490 77 : if (m_tilingScheme == "raster" && !bHasNorthUpSrcGT)
3491 : {
3492 1 : CPLStringList aosOptions;
3493 1 : aosOptions.AddString("-of");
3494 1 : aosOptions.AddString("VRT");
3495 1 : aosOptions.AddString("-a_ullr");
3496 1 : aosOptions.AddString(CPLSPrintf("%.17g", srcGTModif[0]));
3497 1 : aosOptions.AddString(CPLSPrintf("%.17g", srcGTModif[3]));
3498 : aosOptions.AddString(
3499 1 : CPLSPrintf("%.17g", srcGTModif[0] + nSrcWidth * srcGTModif[1]));
3500 : aosOptions.AddString(
3501 1 : CPLSPrintf("%.17g", srcGTModif[3] + nSrcHeight * srcGTModif[5]));
3502 1 : if (oSRS_TMS.IsEmpty())
3503 : {
3504 1 : aosOptions.AddString("-a_srs");
3505 1 : aosOptions.AddString("none");
3506 : }
3507 :
3508 : GDALTranslateOptions *psOptions =
3509 1 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
3510 :
3511 1 : tmpSrcDS.reset(GDALDataset::FromHandle(GDALTranslate(
3512 : "", GDALDataset::ToHandle(m_poSrcDS), psOptions, nullptr)));
3513 1 : GDALTranslateOptionsFree(psOptions);
3514 1 : if (!tmpSrcDS)
3515 0 : return false;
3516 : }
3517 78 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
3518 78 : tmpSrcDS ? tmpSrcDS.get() : m_poSrcDS, &oFakeMaxZoomDS, aosTO.List()));
3519 77 : CPLAssert(hTransformArg);
3520 :
3521 : /* -------------------------------------------------------------------- */
3522 : /* Warp the transformer with a linear approximator */
3523 : /* -------------------------------------------------------------------- */
3524 77 : hTransformArg.reset(GDALCreateApproxTransformer(
3525 : GDALGenImgProjTransform, hTransformArg.release(), 0.125));
3526 77 : GDALApproxTransformerOwnsSubtransformer(hTransformArg.get(), TRUE);
3527 :
3528 77 : psWO->pfnTransformer = GDALApproxTransform;
3529 77 : psWO->pTransformerArg = hTransformArg.get();
3530 :
3531 : /* -------------------------------------------------------------------- */
3532 : /* Determine total number of tiles */
3533 : /* -------------------------------------------------------------------- */
3534 77 : const int nBaseTilesPerRow = nMaxTileX - nMinTileX + 1;
3535 77 : const int nBaseTilesPerCol = nMaxTileY - nMinTileY + 1;
3536 77 : const uint64_t nBaseTiles =
3537 77 : static_cast<uint64_t>(nBaseTilesPerCol) * nBaseTilesPerRow;
3538 77 : uint64_t nTotalTiles = nBaseTiles;
3539 77 : std::atomic<uint64_t> nCurTile = 0;
3540 77 : bool bRet = true;
3541 :
3542 142 : for (int iZ = m_maxZoomLevel - 1;
3543 142 : bRet && bIntersects && iZ >= m_minZoomLevel; --iZ)
3544 : {
3545 130 : auto ovrTileMatrix = tileMatrixList[iZ];
3546 65 : int nOvrMinTileX = 0;
3547 65 : int nOvrMinTileY = 0;
3548 65 : int nOvrMaxTileX = 0;
3549 65 : int nOvrMaxTileY = 0;
3550 65 : bRet =
3551 130 : GetTileIndices(ovrTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
3552 : nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX,
3553 65 : nOvrMaxTileY, m_noIntersectionIsOK, bIntersects);
3554 65 : if (bIntersects)
3555 : {
3556 65 : nTotalTiles +=
3557 65 : static_cast<uint64_t>(nOvrMaxTileY - nOvrMinTileY + 1) *
3558 65 : (nOvrMaxTileX - nOvrMinTileX + 1);
3559 : }
3560 : }
3561 :
3562 : /* -------------------------------------------------------------------- */
3563 : /* Generate tiles at max zoom level */
3564 : /* -------------------------------------------------------------------- */
3565 154 : GDALWarpOperation oWO;
3566 :
3567 77 : bRet = oWO.Initialize(psWO.get()) == CE_None && bRet;
3568 :
3569 : const auto GetUpdatedCreationOptions =
3570 317 : [this](const gdal::TileMatrixSet::TileMatrix &oTM)
3571 : {
3572 113 : CPLStringList aosCreationOptions(m_creationOptions);
3573 113 : if (m_outputFormat == "GTiff")
3574 : {
3575 48 : if (aosCreationOptions.FetchNameValue("TILED") == nullptr &&
3576 24 : aosCreationOptions.FetchNameValue("BLOCKYSIZE") == nullptr)
3577 : {
3578 24 : if (oTM.mTileWidth <= 512 && oTM.mTileHeight <= 512)
3579 : {
3580 22 : aosCreationOptions.SetNameValue(
3581 22 : "BLOCKYSIZE", CPLSPrintf("%d", oTM.mTileHeight));
3582 : }
3583 : else
3584 : {
3585 2 : aosCreationOptions.SetNameValue("TILED", "YES");
3586 : }
3587 : }
3588 24 : if (aosCreationOptions.FetchNameValue("COMPRESS") == nullptr)
3589 24 : aosCreationOptions.SetNameValue("COMPRESS", "LZW");
3590 : }
3591 89 : else if (m_outputFormat == "COG")
3592 : {
3593 2 : if (aosCreationOptions.FetchNameValue("OVERVIEW_RESAMPLING") ==
3594 : nullptr)
3595 : {
3596 : aosCreationOptions.SetNameValue("OVERVIEW_RESAMPLING",
3597 2 : m_overviewResampling.c_str());
3598 : }
3599 2 : if (aosCreationOptions.FetchNameValue("BLOCKSIZE") == nullptr &&
3600 2 : oTM.mTileWidth <= 512 && oTM.mTileWidth == oTM.mTileHeight)
3601 : {
3602 : aosCreationOptions.SetNameValue(
3603 2 : "BLOCKSIZE", CPLSPrintf("%d", oTM.mTileWidth));
3604 : }
3605 : }
3606 113 : return aosCreationOptions;
3607 77 : };
3608 :
3609 77 : VSIMkdir(m_outputDirectory.c_str(), 0755);
3610 : VSIStatBufL sStat;
3611 153 : if (VSIStatL(m_outputDirectory.c_str(), &sStat) != 0 ||
3612 76 : !VSI_ISDIR(sStat.st_mode))
3613 : {
3614 1 : ReportError(CE_Failure, CPLE_FileIO,
3615 : "Cannot create output directory %s",
3616 : m_outputDirectory.c_str());
3617 1 : return false;
3618 : }
3619 :
3620 152 : OGRSpatialReference oWGS84;
3621 76 : oWGS84.importFromEPSG(4326);
3622 76 : oWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3623 :
3624 76 : std::unique_ptr<OGRCoordinateTransformation> poCTToWGS84;
3625 76 : if (!oSRS_TMS.IsEmpty())
3626 : {
3627 150 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3628 75 : poCTToWGS84.reset(
3629 : OGRCreateCoordinateTransformation(&oSRS_TMS, &oWGS84));
3630 : }
3631 :
3632 78 : const bool kmlCompatible = m_kml &&
3633 17 : [this, &poTMS, &poCTToWGS84, bInvertAxisTMS]()
3634 : {
3635 2 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3636 2 : double dfX = poTMS->tileMatrixList()[0].mTopLeftX;
3637 2 : double dfY = poTMS->tileMatrixList()[0].mTopLeftY;
3638 2 : if (bInvertAxisTMS)
3639 0 : std::swap(dfX, dfY);
3640 3 : return (m_minZoomLevel == m_maxZoomLevel ||
3641 1 : (poTMS->haveAllLevelsSameTopLeft() &&
3642 1 : poTMS->haveAllLevelsSameTileSize() &&
3643 3 : poTMS->hasOnlyPowerOfTwoVaryingScales())) &&
3644 7 : poCTToWGS84 && poCTToWGS84->Transform(1, &dfX, &dfY);
3645 2 : }();
3646 : const int kmlTileSize =
3647 76 : m_tileSize > 0 ? m_tileSize : poTMS->tileMatrixList()[0].mTileWidth;
3648 76 : if (m_kml && !kmlCompatible)
3649 : {
3650 0 : ReportError(CE_Failure, CPLE_NotSupported,
3651 : "Tiling scheme not compatible with KML output");
3652 0 : return false;
3653 : }
3654 :
3655 76 : if (m_title.empty())
3656 74 : m_title = CPLGetFilename(m_dataset.GetName().c_str());
3657 :
3658 76 : if (!m_url.empty())
3659 : {
3660 2 : if (m_url.back() != '/')
3661 2 : m_url += '/';
3662 4 : std::string out_path = m_outputDirectory;
3663 2 : if (m_outputDirectory.back() == '/')
3664 0 : out_path.pop_back();
3665 2 : m_url += CPLGetFilename(out_path.c_str());
3666 : }
3667 :
3668 152 : CPLWorkerThreadPool oThreadPool;
3669 :
3670 76 : bool bThreadPoolInitialized = false;
3671 : const auto InitThreadPool =
3672 380 : [this, &oThreadPool, &bRet, &bThreadPoolInitialized]()
3673 : {
3674 113 : if (!bThreadPoolInitialized)
3675 : {
3676 75 : bThreadPoolInitialized = true;
3677 :
3678 75 : if (bRet && m_numThreads > 1)
3679 : {
3680 2 : CPLDebug("gdal_raster_tile", "Using %d threads", m_numThreads);
3681 2 : bRet = oThreadPool.Setup(m_numThreads, nullptr, nullptr);
3682 : }
3683 : }
3684 :
3685 113 : return bRet;
3686 76 : };
3687 :
3688 : // Just for unit test purposes
3689 76 : const bool bEmitSpuriousCharsOnStdout = CPLTestBool(
3690 : CPLGetConfigOption("GDAL_RASTER_TILE_EMIT_SPURIOUS_CHARS", "NO"));
3691 :
3692 6 : const auto IsCompatibleOfSpawnSilent = [this]()
3693 : {
3694 6 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3695 3 : const char *pszErrorMsg = "";
3696 6 : return IsCompatibleOfSpawn(pszErrorMsg);
3697 76 : };
3698 :
3699 76 : m_numThreads = std::max(
3700 152 : 1, static_cast<int>(std::min<uint64_t>(
3701 76 : m_numThreads, nBaseTiles / GetThresholdMinTilesPerJob())));
3702 :
3703 76 : if (m_ovrZoomLevel >= 0)
3704 : {
3705 : // do not generate base tiles if called as a child process with
3706 : // --ovr-zoom-level
3707 : }
3708 76 : else if (m_numThreads > 1 && nBaseTiles > 1 &&
3709 4 : ((m_parallelMethod.empty() &&
3710 2 : m_numThreads >= GetThresholdMinThreadsForSpawn() &&
3711 5 : IsCompatibleOfSpawnSilent()) ||
3712 4 : m_parallelMethod == "spawn"))
3713 : {
3714 2 : if (!GenerateBaseTilesSpawnMethod(nBaseTilesPerCol, nBaseTilesPerRow,
3715 : nMinTileX, nMinTileY, nMaxTileX,
3716 : nMaxTileY, nTotalTiles, nBaseTiles,
3717 : pfnProgress, pProgressData))
3718 : {
3719 1 : return false;
3720 : }
3721 1 : nCurTile = nBaseTiles;
3722 : }
3723 : else
3724 : {
3725 : // Branch for multi-threaded or single-threaded max zoom level tile
3726 : // generation
3727 :
3728 : PerThreadMaxZoomResourceManager oResourceManager(
3729 66 : m_poSrcDS, psWO.get(), hTransformArg.get(), oFakeMaxZoomDS,
3730 198 : dstBuffer.size());
3731 :
3732 : const CPLStringList aosCreationOptions(
3733 132 : GetUpdatedCreationOptions(tileMatrix));
3734 :
3735 66 : CPLDebug("gdal_raster_tile",
3736 : "Generating tiles z=%d, y=%d...%d, x=%d...%d", m_maxZoomLevel,
3737 : nMinTileY, nMaxTileY, nMinTileX, nMaxTileX);
3738 :
3739 66 : bRet &= InitThreadPool();
3740 :
3741 66 : if (bRet && m_numThreads > 1)
3742 : {
3743 2 : std::atomic<bool> bFailure = false;
3744 2 : std::atomic<int> nQueuedJobs = 0;
3745 :
3746 : double dfTilesYPerJob;
3747 : int nYOuterIterations;
3748 : double dfTilesXPerJob;
3749 : int nXOuterIterations;
3750 2 : ComputeJobChunkSize(m_numThreads, nBaseTilesPerCol,
3751 : nBaseTilesPerRow, dfTilesYPerJob,
3752 : nYOuterIterations, dfTilesXPerJob,
3753 : nXOuterIterations);
3754 :
3755 2 : CPLDebugOnly("gdal_raster_tile",
3756 : "nYOuterIterations=%d, dfTilesYPerJob=%g, "
3757 : "nXOuterIterations=%d, dfTilesXPerJob=%g",
3758 : nYOuterIterations, dfTilesYPerJob, nXOuterIterations,
3759 : dfTilesXPerJob);
3760 :
3761 2 : int nLastYEndIncluded = nMinTileY - 1;
3762 10 : for (int iYOuterIter = 0; bRet && iYOuterIter < nYOuterIterations &&
3763 8 : nLastYEndIncluded < nMaxTileY;
3764 : ++iYOuterIter)
3765 : {
3766 8 : const int iYStart = nLastYEndIncluded + 1;
3767 : const int iYEndIncluded =
3768 8 : iYOuterIter + 1 == nYOuterIterations
3769 14 : ? nMaxTileY
3770 : : std::max(
3771 : iYStart,
3772 14 : static_cast<int>(std::floor(
3773 6 : nMinTileY +
3774 6 : (iYOuterIter + 1) * dfTilesYPerJob - 1)));
3775 :
3776 8 : nLastYEndIncluded = iYEndIncluded;
3777 :
3778 8 : int nLastXEndIncluded = nMinTileX - 1;
3779 8 : for (int iXOuterIter = 0;
3780 16 : bRet && iXOuterIter < nXOuterIterations &&
3781 8 : nLastXEndIncluded < nMaxTileX;
3782 : ++iXOuterIter)
3783 : {
3784 8 : const int iXStart = nLastXEndIncluded + 1;
3785 : const int iXEndIncluded =
3786 8 : iXOuterIter + 1 == nXOuterIterations
3787 8 : ? nMaxTileX
3788 : : std::max(
3789 : iXStart,
3790 8 : static_cast<int>(std::floor(
3791 0 : nMinTileX +
3792 0 : (iXOuterIter + 1) * dfTilesXPerJob - 1)));
3793 :
3794 8 : nLastXEndIncluded = iXEndIncluded;
3795 :
3796 8 : CPLDebugOnly("gdal_raster_tile",
3797 : "Job for y in [%d,%d] and x in [%d,%d]",
3798 : iYStart, iYEndIncluded, iXStart,
3799 : iXEndIncluded);
3800 :
3801 8 : auto job = [this, &oThreadPool, &oResourceManager,
3802 : &bFailure, &nCurTile, &nQueuedJobs,
3803 : pszExtension, &aosCreationOptions, &psWO,
3804 : &tileMatrix, nDstBands, iXStart, iXEndIncluded,
3805 : iYStart, iYEndIncluded, nMinTileX, nMinTileY,
3806 1088 : &poColorTable, bUserAskedForAlpha]()
3807 : {
3808 8 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3809 :
3810 8 : auto resources = oResourceManager.AcquireResources();
3811 8 : if (resources)
3812 : {
3813 24 : for (int iY = iYStart; iY <= iYEndIncluded; ++iY)
3814 : {
3815 144 : for (int iX = iXStart; iX <= iXEndIncluded;
3816 : ++iX)
3817 : {
3818 384 : if (!GenerateTile(
3819 128 : resources->poSrcDS.get(),
3820 : m_poDstDriver, pszExtension,
3821 : aosCreationOptions.List(),
3822 128 : *(resources->poWO.get()),
3823 128 : *(resources->poFakeMaxZoomDS
3824 128 : ->GetSpatialRef()),
3825 128 : psWO->eWorkingDataType, tileMatrix,
3826 128 : m_outputDirectory, nDstBands,
3827 128 : psWO->padfDstNoDataReal
3828 0 : ? &(psWO->padfDstNoDataReal[0])
3829 : : nullptr,
3830 : m_maxZoomLevel, iX, iY,
3831 128 : m_convention, nMinTileX, nMinTileY,
3832 128 : m_skipBlank, bUserAskedForAlpha,
3833 128 : m_auxXML, m_resume, m_metadata,
3834 128 : poColorTable.get(),
3835 128 : resources->dstBuffer))
3836 : {
3837 0 : oResourceManager.SetError();
3838 0 : bFailure = true;
3839 0 : --nQueuedJobs;
3840 0 : return;
3841 : }
3842 128 : ++nCurTile;
3843 128 : oThreadPool.WakeUpWaitEvent();
3844 : }
3845 : }
3846 8 : oResourceManager.ReleaseResources(
3847 8 : std::move(resources));
3848 : }
3849 : else
3850 : {
3851 0 : oResourceManager.SetError();
3852 0 : bFailure = true;
3853 : }
3854 :
3855 8 : --nQueuedJobs;
3856 8 : };
3857 :
3858 8 : ++nQueuedJobs;
3859 8 : oThreadPool.SubmitJob(std::move(job));
3860 : }
3861 : }
3862 :
3863 : // Wait for completion of all jobs
3864 134 : while (bRet && nQueuedJobs > 0)
3865 : {
3866 132 : oThreadPool.WaitEvent();
3867 197 : bRet &= !bFailure &&
3868 65 : (!pfnProgress ||
3869 130 : pfnProgress(static_cast<double>(nCurTile) /
3870 65 : static_cast<double>(nTotalTiles),
3871 132 : "", pProgressData));
3872 : }
3873 2 : oThreadPool.WaitCompletion();
3874 2 : bRet &=
3875 3 : !bFailure && (!pfnProgress ||
3876 2 : pfnProgress(static_cast<double>(nCurTile) /
3877 1 : static_cast<double>(nTotalTiles),
3878 2 : "", pProgressData));
3879 :
3880 2 : if (!oResourceManager.GetErrorMsg().empty())
3881 : {
3882 : // Re-emit error message from worker thread to main thread
3883 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
3884 0 : oResourceManager.GetErrorMsg().c_str());
3885 2 : }
3886 : }
3887 : else
3888 : {
3889 : // Branch for single-thread max zoom level tile generation
3890 160 : for (int iY = nMinTileY; bRet && iY <= nMaxTileY; ++iY)
3891 : {
3892 346 : for (int iX = nMinTileX; bRet && iX <= nMaxTileX; ++iX)
3893 : {
3894 500 : bRet = GenerateTile(
3895 : m_poSrcDS, m_poDstDriver, pszExtension,
3896 : aosCreationOptions.List(), oWO, oSRS_TMS,
3897 250 : psWO->eWorkingDataType, tileMatrix, m_outputDirectory,
3898 : nDstBands,
3899 250 : psWO->padfDstNoDataReal ? &(psWO->padfDstNoDataReal[0])
3900 : : nullptr,
3901 250 : m_maxZoomLevel, iX, iY, m_convention, nMinTileX,
3902 250 : nMinTileY, m_skipBlank, bUserAskedForAlpha, m_auxXML,
3903 250 : m_resume, m_metadata, poColorTable.get(), dstBuffer);
3904 :
3905 250 : if (m_progressForked)
3906 : {
3907 66 : if (bEmitSpuriousCharsOnStdout)
3908 64 : fwrite(&PROGRESS_MARKER[0], 1, 1, stdout);
3909 66 : fwrite(PROGRESS_MARKER, sizeof(PROGRESS_MARKER), 1,
3910 : stdout);
3911 66 : fflush(stdout);
3912 : }
3913 : else
3914 : {
3915 184 : ++nCurTile;
3916 184 : bRet &=
3917 186 : (!pfnProgress ||
3918 4 : pfnProgress(static_cast<double>(nCurTile) /
3919 2 : static_cast<double>(nTotalTiles),
3920 184 : "", pProgressData));
3921 : }
3922 : }
3923 : }
3924 : }
3925 :
3926 66 : if (m_kml && bRet)
3927 : {
3928 4 : for (int iY = nMinTileY; iY <= nMaxTileY; ++iY)
3929 : {
3930 4 : for (int iX = nMinTileX; iX <= nMaxTileX; ++iX)
3931 : {
3932 : const int nFileY =
3933 2 : GetFileY(iY, poTMS->tileMatrixList()[m_maxZoomLevel],
3934 2 : m_convention);
3935 : std::string osFilename = CPLFormFilenameSafe(
3936 : m_outputDirectory.c_str(),
3937 4 : CPLSPrintf("%d", m_maxZoomLevel), nullptr);
3938 4 : osFilename = CPLFormFilenameSafe(
3939 2 : osFilename.c_str(), CPLSPrintf("%d", iX), nullptr);
3940 4 : osFilename = CPLFormFilenameSafe(
3941 : osFilename.c_str(),
3942 2 : CPLSPrintf("%d.%s", nFileY, pszExtension), nullptr);
3943 2 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
3944 : {
3945 4 : GenerateKML(m_outputDirectory, m_title, iX, iY,
3946 : m_maxZoomLevel, kmlTileSize, pszExtension,
3947 2 : m_url, poTMS.get(), bInvertAxisTMS,
3948 2 : m_convention, poCTToWGS84.get(), {});
3949 : }
3950 : }
3951 : }
3952 : }
3953 : }
3954 :
3955 : // Close source dataset if we have opened it (in GDALAlgorithm core code),
3956 : // to free file descriptors, particularly if it is a VRT file.
3957 75 : std::vector<GDALColorInterp> aeColorInterp;
3958 247 : for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
3959 172 : aeColorInterp.push_back(
3960 172 : m_poSrcDS->GetRasterBand(i)->GetColorInterpretation());
3961 75 : if (m_dataset.HasDatasetBeenOpenedByAlgorithm())
3962 : {
3963 37 : m_dataset.Close();
3964 37 : m_poSrcDS = nullptr;
3965 : }
3966 :
3967 : /* -------------------------------------------------------------------- */
3968 : /* Generate tiles at lower zoom levels */
3969 : /* -------------------------------------------------------------------- */
3970 75 : const int iZStart =
3971 75 : m_ovrZoomLevel >= 0 ? m_ovrZoomLevel : m_maxZoomLevel - 1;
3972 75 : const int iZEnd = m_ovrZoomLevel >= 0 ? m_ovrZoomLevel : m_minZoomLevel;
3973 124 : for (int iZ = iZStart; bRet && iZ >= iZEnd; --iZ)
3974 : {
3975 49 : int nOvrMinTileX = 0;
3976 49 : int nOvrMinTileY = 0;
3977 49 : int nOvrMaxTileX = 0;
3978 49 : int nOvrMaxTileY = 0;
3979 :
3980 98 : auto ovrTileMatrix = tileMatrixList[iZ];
3981 49 : CPL_IGNORE_RET_VAL(
3982 49 : GetTileIndices(ovrTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
3983 : nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX,
3984 49 : nOvrMaxTileY, m_noIntersectionIsOK, bIntersects));
3985 :
3986 49 : bRet = bIntersects;
3987 :
3988 49 : if (m_minOvrTileX >= 0)
3989 : {
3990 8 : bRet = true;
3991 8 : nOvrMinTileX = m_minOvrTileX;
3992 8 : nOvrMinTileY = m_minOvrTileY;
3993 8 : nOvrMaxTileX = m_maxOvrTileX;
3994 8 : nOvrMaxTileY = m_maxOvrTileY;
3995 : }
3996 :
3997 49 : if (bRet)
3998 : {
3999 49 : CPLDebug("gdal_raster_tile",
4000 : "Generating overview tiles z=%d, y=%d...%d, x=%d...%d", iZ,
4001 : nOvrMinTileY, nOvrMaxTileY, nOvrMinTileX, nOvrMaxTileX);
4002 : }
4003 :
4004 49 : const int nOvrTilesPerCol = nOvrMaxTileY - nOvrMinTileY + 1;
4005 49 : const int nOvrTilesPerRow = nOvrMaxTileX - nOvrMinTileX + 1;
4006 49 : const uint64_t nOvrTileCount =
4007 49 : static_cast<uint64_t>(nOvrTilesPerCol) * nOvrTilesPerRow;
4008 :
4009 49 : m_numThreads = std::max(
4010 98 : 1,
4011 98 : static_cast<int>(std::min<uint64_t>(
4012 49 : m_numThreads, nOvrTileCount / GetThresholdMinTilesPerJob())));
4013 :
4014 61 : if (m_numThreads > 1 && nOvrTileCount > 1 &&
4015 6 : ((m_parallelMethod.empty() &&
4016 4 : m_numThreads >= GetThresholdMinThreadsForSpawn() &&
4017 8 : IsCompatibleOfSpawnSilent()) ||
4018 6 : m_parallelMethod == "spawn"))
4019 : {
4020 2 : bRet &= GenerateOverviewTilesSpawnMethod(
4021 : iZ, nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX, nOvrMaxTileY,
4022 2 : nCurTile, nTotalTiles, pfnProgress, pProgressData);
4023 : }
4024 : else
4025 : {
4026 47 : bRet &= InitThreadPool();
4027 :
4028 94 : auto srcTileMatrix = tileMatrixList[iZ + 1];
4029 47 : int nSrcMinTileX = 0;
4030 47 : int nSrcMinTileY = 0;
4031 47 : int nSrcMaxTileX = 0;
4032 47 : int nSrcMaxTileY = 0;
4033 :
4034 47 : CPL_IGNORE_RET_VAL(GetTileIndices(
4035 : srcTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
4036 : nSrcMinTileX, nSrcMinTileY, nSrcMaxTileX, nSrcMaxTileY,
4037 47 : m_noIntersectionIsOK, bIntersects));
4038 :
4039 47 : constexpr double EPSILON = 1e-3;
4040 47 : int maxCacheTileSizePerThread = static_cast<int>(
4041 47 : (1 + std::ceil(
4042 47 : (ovrTileMatrix.mResY * ovrTileMatrix.mTileHeight) /
4043 47 : (srcTileMatrix.mResY * srcTileMatrix.mTileHeight) -
4044 47 : EPSILON)) *
4045 47 : (1 + std::ceil(
4046 47 : (ovrTileMatrix.mResX * ovrTileMatrix.mTileWidth) /
4047 47 : (srcTileMatrix.mResX * srcTileMatrix.mTileWidth) -
4048 : EPSILON)));
4049 :
4050 47 : CPLDebugOnly("gdal_raster_tile",
4051 : "Ideal maxCacheTileSizePerThread = %d",
4052 : maxCacheTileSizePerThread);
4053 :
4054 : #ifndef _WIN32
4055 : const int remainingFileDescriptorCount =
4056 47 : CPLGetRemainingFileDescriptorCount();
4057 47 : CPLDebugOnly("gdal_raster_tile",
4058 : "remainingFileDescriptorCount = %d",
4059 : remainingFileDescriptorCount);
4060 47 : if (remainingFileDescriptorCount >= 0 &&
4061 : remainingFileDescriptorCount <
4062 47 : (1 + maxCacheTileSizePerThread) * m_numThreads)
4063 : {
4064 : const int newNumThreads =
4065 0 : std::max(1, remainingFileDescriptorCount /
4066 0 : (1 + maxCacheTileSizePerThread));
4067 0 : if (newNumThreads < m_numThreads)
4068 : {
4069 0 : CPLError(CE_Warning, CPLE_AppDefined,
4070 : "Not enough file descriptors available given the "
4071 : "number of "
4072 : "threads. Reducing the number of threads %d to %d",
4073 : m_numThreads, newNumThreads);
4074 0 : m_numThreads = newNumThreads;
4075 : }
4076 : }
4077 : #endif
4078 :
4079 : MosaicDataset oSrcDS(
4080 94 : CPLFormFilenameSafe(m_outputDirectory.c_str(),
4081 : CPLSPrintf("%d", iZ + 1), nullptr),
4082 47 : pszExtension, m_outputFormat, aeColorInterp, srcTileMatrix,
4083 : oSRS_TMS, nSrcMinTileX, nSrcMinTileY, nSrcMaxTileX,
4084 47 : nSrcMaxTileY, m_convention, nDstBands, psWO->eWorkingDataType,
4085 9 : psWO->padfDstNoDataReal ? &(psWO->padfDstNoDataReal[0])
4086 : : nullptr,
4087 244 : m_metadata, poColorTable.get(), maxCacheTileSizePerThread);
4088 :
4089 : const CPLStringList aosCreationOptions(
4090 94 : GetUpdatedCreationOptions(ovrTileMatrix));
4091 :
4092 94 : PerThreadLowerZoomResourceManager oResourceManager(oSrcDS);
4093 47 : std::atomic<bool> bFailure = false;
4094 47 : std::atomic<int> nQueuedJobs = 0;
4095 :
4096 47 : const bool bUseThreads = m_numThreads > 1 && nOvrTileCount > 1;
4097 :
4098 47 : if (bUseThreads)
4099 : {
4100 : double dfTilesYPerJob;
4101 : int nYOuterIterations;
4102 : double dfTilesXPerJob;
4103 : int nXOuterIterations;
4104 4 : ComputeJobChunkSize(m_numThreads, nOvrTilesPerCol,
4105 : nOvrTilesPerRow, dfTilesYPerJob,
4106 : nYOuterIterations, dfTilesXPerJob,
4107 : nXOuterIterations);
4108 :
4109 4 : CPLDebugOnly("gdal_raster_tile",
4110 : "z=%d, nYOuterIterations=%d, dfTilesYPerJob=%g, "
4111 : "nXOuterIterations=%d, dfTilesXPerJob=%g",
4112 : iZ, nYOuterIterations, dfTilesYPerJob,
4113 : nXOuterIterations, dfTilesXPerJob);
4114 :
4115 4 : int nLastYEndIncluded = nOvrMinTileY - 1;
4116 16 : for (int iYOuterIter = 0;
4117 16 : bRet && iYOuterIter < nYOuterIterations &&
4118 12 : nLastYEndIncluded < nOvrMaxTileY;
4119 : ++iYOuterIter)
4120 : {
4121 12 : const int iYStart = nLastYEndIncluded + 1;
4122 : const int iYEndIncluded =
4123 12 : iYOuterIter + 1 == nYOuterIterations
4124 20 : ? nOvrMaxTileY
4125 : : std::max(
4126 : iYStart,
4127 20 : static_cast<int>(std::floor(
4128 8 : nOvrMinTileY +
4129 8 : (iYOuterIter + 1) * dfTilesYPerJob - 1)));
4130 :
4131 12 : nLastYEndIncluded = iYEndIncluded;
4132 :
4133 12 : int nLastXEndIncluded = nOvrMinTileX - 1;
4134 12 : for (int iXOuterIter = 0;
4135 28 : bRet && iXOuterIter < nXOuterIterations &&
4136 16 : nLastXEndIncluded < nOvrMaxTileX;
4137 : ++iXOuterIter)
4138 : {
4139 16 : const int iXStart = nLastXEndIncluded + 1;
4140 : const int iXEndIncluded =
4141 16 : iXOuterIter + 1 == nXOuterIterations
4142 20 : ? nOvrMaxTileX
4143 20 : : std::max(iXStart, static_cast<int>(std::floor(
4144 4 : nOvrMinTileX +
4145 4 : (iXOuterIter + 1) *
4146 : dfTilesXPerJob -
4147 4 : 1)));
4148 :
4149 16 : nLastXEndIncluded = iXEndIncluded;
4150 :
4151 16 : CPLDebugOnly(
4152 : "gdal_raster_tile",
4153 : "Job for z=%d, y in [%d,%d] and x in [%d,%d]", iZ,
4154 : iYStart, iYEndIncluded, iXStart, iXEndIncluded);
4155 16 : auto job = [this, &oThreadPool, &oResourceManager,
4156 : &bFailure, &nCurTile, &nQueuedJobs,
4157 : pszExtension, &aosCreationOptions,
4158 : &aosWarpOptions, &ovrTileMatrix, iZ,
4159 : iXStart, iXEndIncluded, iYStart,
4160 336 : iYEndIncluded, bUserAskedForAlpha]()
4161 : {
4162 : CPLErrorStateBackuper oBackuper(
4163 16 : CPLQuietErrorHandler);
4164 :
4165 : auto resources =
4166 16 : oResourceManager.AcquireResources();
4167 16 : if (resources)
4168 : {
4169 32 : for (int iY = iYStart; iY <= iYEndIncluded;
4170 : ++iY)
4171 : {
4172 56 : for (int iX = iXStart; iX <= iXEndIncluded;
4173 : ++iX)
4174 : {
4175 80 : if (!GenerateOverviewTile(
4176 40 : *(resources->poSrcDS.get()),
4177 40 : m_poDstDriver, m_outputFormat,
4178 : pszExtension,
4179 : aosCreationOptions.List(),
4180 40 : aosWarpOptions.List(),
4181 40 : m_overviewResampling,
4182 : ovrTileMatrix,
4183 40 : m_outputDirectory, iZ, iX, iY,
4184 40 : m_convention, m_skipBlank,
4185 40 : bUserAskedForAlpha, m_auxXML,
4186 40 : m_resume))
4187 : {
4188 0 : oResourceManager.SetError();
4189 0 : bFailure = true;
4190 0 : --nQueuedJobs;
4191 0 : return;
4192 : }
4193 :
4194 40 : ++nCurTile;
4195 40 : oThreadPool.WakeUpWaitEvent();
4196 : }
4197 : }
4198 16 : oResourceManager.ReleaseResources(
4199 16 : std::move(resources));
4200 : }
4201 : else
4202 : {
4203 0 : oResourceManager.SetError();
4204 0 : bFailure = true;
4205 : }
4206 16 : --nQueuedJobs;
4207 16 : };
4208 :
4209 16 : ++nQueuedJobs;
4210 16 : oThreadPool.SubmitJob(std::move(job));
4211 : }
4212 : }
4213 :
4214 : // Wait for completion of all jobs
4215 51 : while (bRet && nQueuedJobs > 0)
4216 : {
4217 47 : oThreadPool.WaitEvent();
4218 69 : bRet &= !bFailure &&
4219 22 : (!pfnProgress ||
4220 44 : pfnProgress(static_cast<double>(nCurTile) /
4221 22 : static_cast<double>(nTotalTiles),
4222 47 : "", pProgressData));
4223 : }
4224 4 : oThreadPool.WaitCompletion();
4225 6 : bRet &= !bFailure &&
4226 2 : (!pfnProgress ||
4227 4 : pfnProgress(static_cast<double>(nCurTile) /
4228 2 : static_cast<double>(nTotalTiles),
4229 4 : "", pProgressData));
4230 :
4231 4 : if (!oResourceManager.GetErrorMsg().empty())
4232 : {
4233 : // Re-emit error message from worker thread to main thread
4234 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
4235 0 : oResourceManager.GetErrorMsg().c_str());
4236 : }
4237 : }
4238 : else
4239 : {
4240 : // Branch for single-thread overview generation
4241 :
4242 90 : for (int iY = nOvrMinTileY; bRet && iY <= nOvrMaxTileY; ++iY)
4243 : {
4244 120 : for (int iX = nOvrMinTileX; bRet && iX <= nOvrMaxTileX;
4245 : ++iX)
4246 : {
4247 73 : bRet = GenerateOverviewTile(
4248 73 : oSrcDS, m_poDstDriver, m_outputFormat, pszExtension,
4249 73 : aosCreationOptions.List(), aosWarpOptions.List(),
4250 73 : m_overviewResampling, ovrTileMatrix,
4251 73 : m_outputDirectory, iZ, iX, iY, m_convention,
4252 73 : m_skipBlank, bUserAskedForAlpha, m_auxXML,
4253 73 : m_resume);
4254 :
4255 73 : if (m_progressForked)
4256 : {
4257 20 : if (bEmitSpuriousCharsOnStdout)
4258 20 : fwrite(&PROGRESS_MARKER[0], 1, 1, stdout);
4259 20 : fwrite(PROGRESS_MARKER, sizeof(PROGRESS_MARKER), 1,
4260 : stdout);
4261 20 : fflush(stdout);
4262 : }
4263 : else
4264 : {
4265 53 : ++nCurTile;
4266 55 : bRet &= (!pfnProgress ||
4267 2 : pfnProgress(
4268 2 : static_cast<double>(nCurTile) /
4269 2 : static_cast<double>(nTotalTiles),
4270 53 : "", pProgressData));
4271 : }
4272 : }
4273 : }
4274 : }
4275 : }
4276 :
4277 49 : if (m_kml && bRet)
4278 : {
4279 2 : for (int iY = nOvrMinTileY; bRet && iY <= nOvrMaxTileY; ++iY)
4280 : {
4281 2 : for (int iX = nOvrMinTileX; bRet && iX <= nOvrMaxTileX; ++iX)
4282 : {
4283 : int nFileY =
4284 1 : GetFileY(iY, poTMS->tileMatrixList()[iZ], m_convention);
4285 : std::string osFilename =
4286 : CPLFormFilenameSafe(m_outputDirectory.c_str(),
4287 2 : CPLSPrintf("%d", iZ), nullptr);
4288 2 : osFilename = CPLFormFilenameSafe(
4289 1 : osFilename.c_str(), CPLSPrintf("%d", iX), nullptr);
4290 2 : osFilename = CPLFormFilenameSafe(
4291 : osFilename.c_str(),
4292 1 : CPLSPrintf("%d.%s", nFileY, pszExtension), nullptr);
4293 1 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
4294 : {
4295 1 : std::vector<TileCoordinates> children;
4296 :
4297 3 : for (int iChildY = 0; iChildY <= 1; ++iChildY)
4298 : {
4299 6 : for (int iChildX = 0; iChildX <= 1; ++iChildX)
4300 : {
4301 : nFileY =
4302 4 : GetFileY(iY * 2 + iChildY,
4303 4 : poTMS->tileMatrixList()[iZ + 1],
4304 4 : m_convention);
4305 8 : osFilename = CPLFormFilenameSafe(
4306 : m_outputDirectory.c_str(),
4307 4 : CPLSPrintf("%d", iZ + 1), nullptr);
4308 8 : osFilename = CPLFormFilenameSafe(
4309 : osFilename.c_str(),
4310 4 : CPLSPrintf("%d", iX * 2 + iChildX),
4311 4 : nullptr);
4312 8 : osFilename = CPLFormFilenameSafe(
4313 : osFilename.c_str(),
4314 : CPLSPrintf("%d.%s", nFileY, pszExtension),
4315 4 : nullptr);
4316 4 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
4317 : {
4318 1 : TileCoordinates tc;
4319 1 : tc.nTileX = iX * 2 + iChildX;
4320 1 : tc.nTileY = iY * 2 + iChildY;
4321 1 : tc.nTileZ = iZ + 1;
4322 1 : children.push_back(std::move(tc));
4323 : }
4324 : }
4325 : }
4326 :
4327 2 : GenerateKML(m_outputDirectory, m_title, iX, iY, iZ,
4328 1 : kmlTileSize, pszExtension, m_url,
4329 1 : poTMS.get(), bInvertAxisTMS, m_convention,
4330 : poCTToWGS84.get(), children);
4331 : }
4332 : }
4333 : }
4334 : }
4335 : }
4336 :
4337 504 : const auto IsWebViewerEnabled = [this](const char *name)
4338 : {
4339 168 : return std::find_if(m_webviewers.begin(), m_webviewers.end(),
4340 283 : [name](const std::string &s) {
4341 170 : return s == "all" || s == name;
4342 168 : }) != m_webviewers.end();
4343 75 : };
4344 :
4345 129 : if (m_ovrZoomLevel < 0 && bRet &&
4346 204 : poTMS->identifier() == "GoogleMapsCompatible" &&
4347 44 : IsWebViewerEnabled("leaflet"))
4348 : {
4349 15 : double dfSouthLat = -90;
4350 15 : double dfWestLon = -180;
4351 15 : double dfNorthLat = 90;
4352 15 : double dfEastLon = 180;
4353 :
4354 15 : if (poCTToWGS84)
4355 : {
4356 15 : poCTToWGS84->TransformBounds(
4357 : adfExtent[0], adfExtent[1], adfExtent[2], adfExtent[3],
4358 15 : &dfWestLon, &dfSouthLat, &dfEastLon, &dfNorthLat, 21);
4359 : }
4360 :
4361 15 : GenerateLeaflet(m_outputDirectory, m_title, dfSouthLat, dfWestLon,
4362 : dfNorthLat, dfEastLon, m_minZoomLevel, m_maxZoomLevel,
4363 15 : tileMatrix.mTileWidth, pszExtension, m_url, m_copyright,
4364 15 : m_convention == "xyz");
4365 : }
4366 :
4367 75 : if (m_ovrZoomLevel < 0 && bRet && IsWebViewerEnabled("openlayers"))
4368 : {
4369 46 : GenerateOpenLayers(
4370 23 : m_outputDirectory, m_title, adfExtent[0], adfExtent[1],
4371 : adfExtent[2], adfExtent[3], m_minZoomLevel, m_maxZoomLevel,
4372 23 : tileMatrix.mTileWidth, pszExtension, m_url, m_copyright,
4373 23 : *(poTMS.get()), bInvertAxisTMS, oSRS_TMS, m_convention == "xyz");
4374 : }
4375 :
4376 88 : if (m_ovrZoomLevel < 0 && bRet && IsWebViewerEnabled("mapml") &&
4377 163 : poTMS->identifier() != "raster" && m_convention == "xyz")
4378 : {
4379 16 : GenerateMapML(m_outputDirectory, m_mapmlTemplate, m_title, nMinTileX,
4380 : nMinTileY, nMaxTileX, nMaxTileY, m_minZoomLevel,
4381 16 : m_maxZoomLevel, pszExtension, m_url, m_copyright,
4382 16 : *(poTMS.get()));
4383 : }
4384 :
4385 75 : if (m_ovrZoomLevel < 0 && bRet && m_kml)
4386 : {
4387 4 : std::vector<TileCoordinates> children;
4388 :
4389 2 : auto ovrTileMatrix = tileMatrixList[m_minZoomLevel];
4390 2 : int nOvrMinTileX = 0;
4391 2 : int nOvrMinTileY = 0;
4392 2 : int nOvrMaxTileX = 0;
4393 2 : int nOvrMaxTileY = 0;
4394 2 : CPL_IGNORE_RET_VAL(
4395 2 : GetTileIndices(ovrTileMatrix, bInvertAxisTMS, m_tileSize, adfExtent,
4396 : nOvrMinTileX, nOvrMinTileY, nOvrMaxTileX,
4397 2 : nOvrMaxTileY, m_noIntersectionIsOK, bIntersects));
4398 :
4399 4 : for (int iY = nOvrMinTileY; bRet && iY <= nOvrMaxTileY; ++iY)
4400 : {
4401 4 : for (int iX = nOvrMinTileX; bRet && iX <= nOvrMaxTileX; ++iX)
4402 : {
4403 2 : int nFileY = GetFileY(
4404 2 : iY, poTMS->tileMatrixList()[m_minZoomLevel], m_convention);
4405 : std::string osFilename = CPLFormFilenameSafe(
4406 : m_outputDirectory.c_str(), CPLSPrintf("%d", m_minZoomLevel),
4407 4 : nullptr);
4408 4 : osFilename = CPLFormFilenameSafe(osFilename.c_str(),
4409 2 : CPLSPrintf("%d", iX), nullptr);
4410 4 : osFilename = CPLFormFilenameSafe(
4411 : osFilename.c_str(),
4412 2 : CPLSPrintf("%d.%s", nFileY, pszExtension), nullptr);
4413 2 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
4414 : {
4415 2 : TileCoordinates tc;
4416 2 : tc.nTileX = iX;
4417 2 : tc.nTileY = iY;
4418 2 : tc.nTileZ = m_minZoomLevel;
4419 2 : children.push_back(std::move(tc));
4420 : }
4421 : }
4422 : }
4423 4 : GenerateKML(m_outputDirectory, m_title, -1, -1, -1, kmlTileSize,
4424 2 : pszExtension, m_url, poTMS.get(), bInvertAxisTMS,
4425 2 : m_convention, poCTToWGS84.get(), children);
4426 : }
4427 :
4428 75 : if (!bRet && CPLGetLastErrorType() == CE_None)
4429 : {
4430 : // If that happens, this is a programming error
4431 0 : ReportError(CE_Failure, CPLE_AppDefined,
4432 : "Bug: process failed without returning an error message");
4433 : }
4434 :
4435 75 : return bRet;
4436 : }
4437 :
4438 : //! @endcond
|
