Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "mdim mosaic" 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_mdim_mosaic.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "cpl_vsi_virtual.h"
17 : #include "gdal_priv.h"
18 : #include "vrtdataset.h"
19 :
20 : #include <algorithm>
21 : #include <cmath>
22 : #include <optional>
23 :
24 : //! @cond Doxygen_Suppress
25 :
26 : #ifndef _
27 : #define _(x) (x)
28 : #endif
29 :
30 : /************************************************************************/
31 : /* GDALMdimMosaicAlgorithm::GDALMdimMosaicAlgorithm() */
32 : /************************************************************************/
33 :
34 30 : GDALMdimMosaicAlgorithm::GDALMdimMosaicAlgorithm()
35 30 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
36 : {
37 30 : AddProgressArg();
38 30 : AddOutputFormatArg(&m_outputFormat)
39 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES,
40 60 : {GDAL_DCAP_CREATE_MULTIDIMENSIONAL});
41 30 : AddOpenOptionsArg(&m_openOptions);
42 30 : AddInputFormatsArg(&m_inputFormats)
43 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES,
44 60 : {GDAL_ALG_DCAP_RASTER_OR_MULTIDIM_RASTER});
45 30 : AddInputDatasetArg(&m_inputDatasets, GDAL_OF_MULTIDIM_RASTER)
46 30 : .SetDatasetInputFlags(GADV_NAME)
47 30 : .SetDatasetOutputFlags(0)
48 30 : .SetAutoOpenDataset(false)
49 30 : .SetMinCount(1);
50 30 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_MULTIDIM_RASTER);
51 30 : AddCreationOptionsArg(&m_creationOptions);
52 30 : AddOverwriteArg(&m_overwrite);
53 30 : AddArrayNameArg(&m_array, _("Name of array(s) to mosaic."));
54 30 : }
55 :
56 : /************************************************************************/
57 : /* GetDimensionDesc() */
58 : /************************************************************************/
59 :
60 : std::optional<GDALMdimMosaicAlgorithm::DimensionDesc>
61 59 : GDALMdimMosaicAlgorithm::GetDimensionDesc(
62 : const std::string &osDSName,
63 : const std::shared_ptr<GDALDimension> &poDim) const
64 : {
65 118 : auto poVar = poDim->GetIndexingVariable();
66 59 : if (!poVar)
67 : {
68 1 : ReportError(CE_Failure, CPLE_AppDefined,
69 : "Dataset %s: dimension %s lacks an indexing variable",
70 1 : osDSName.c_str(), poDim->GetName().c_str());
71 1 : return std::nullopt;
72 : }
73 58 : if (poVar->GetDimensionCount() != 1)
74 : {
75 0 : ReportError(
76 : CE_Failure, CPLE_AppDefined,
77 : "Dataset %s: indexing variable %s of dimension %s is not 1D",
78 0 : osDSName.c_str(), poVar->GetName().c_str(),
79 0 : poDim->GetName().c_str());
80 0 : return std::nullopt;
81 : }
82 58 : if (poVar->GetDataType().GetClass() != GEDTC_NUMERIC)
83 : {
84 2 : ReportError(CE_Failure, CPLE_AppDefined,
85 : "Dataset %s: indexing variable %s of dimension %s has a "
86 : "non-numeric data type",
87 1 : osDSName.c_str(), poVar->GetName().c_str(),
88 1 : poDim->GetName().c_str());
89 1 : return std::nullopt;
90 : }
91 114 : DimensionDesc desc;
92 57 : desc.osName = poDim->GetName();
93 57 : desc.osType = poDim->GetType();
94 57 : desc.osDirection = poDim->GetDirection();
95 57 : const auto nSize = poVar->GetDimensions()[0]->GetSize();
96 57 : desc.attributes = poVar->GetAttributes();
97 57 : CPLAssert(nSize > 0);
98 57 : desc.nSize = nSize;
99 57 : if (nSize <= 2 || !poVar->IsRegularlySpaced(desc.dfStart, desc.dfIncrement))
100 : {
101 43 : constexpr uint64_t LIMIT = 100 * 1000 * 1000;
102 43 : if (nSize > LIMIT)
103 : {
104 0 : ReportError(CE_Failure, CPLE_AppDefined,
105 : "Dataset %s: indexing variable %s of dimension %s has "
106 : "too large size",
107 0 : osDSName.c_str(), poVar->GetName().c_str(),
108 : desc.osName.c_str());
109 0 : return std::nullopt;
110 : }
111 43 : std::vector<double> adfValues(static_cast<size_t>(nSize));
112 43 : GUInt64 arrayStartIdx[] = {0};
113 43 : size_t anCount[] = {adfValues.size()};
114 43 : GInt64 arrayStep[] = {1};
115 43 : GPtrDiff_t bufferStride[] = {1};
116 86 : if (!poVar->Read(arrayStartIdx, anCount, arrayStep, bufferStride,
117 86 : GDALExtendedDataType::Create(GDT_Float64),
118 43 : adfValues.data()))
119 : {
120 0 : return std::nullopt;
121 : }
122 43 : if (std::isnan(adfValues[0]))
123 : {
124 0 : ReportError(CE_Failure, CPLE_AppDefined,
125 : "Dataset %s: indexing variable %s of dimension %s has "
126 : "NaN values",
127 0 : osDSName.c_str(), poVar->GetName().c_str(),
128 : desc.osName.c_str());
129 0 : return std::nullopt;
130 : }
131 43 : if (nSize > 1)
132 : {
133 13 : const int nSign = adfValues[1] > adfValues[0] ? 1 : -1;
134 13 : if (std::isnan(adfValues[1]))
135 : {
136 0 : ReportError(CE_Failure, CPLE_AppDefined,
137 : "Dataset %s: indexing variable %s of dimension %s "
138 : "has NaN values",
139 0 : osDSName.c_str(), poVar->GetName().c_str(),
140 : desc.osName.c_str());
141 0 : return std::nullopt;
142 : }
143 25 : for (size_t i = 2; i < nSize; ++i)
144 : {
145 12 : if (std::isnan(adfValues[i]))
146 : {
147 0 : ReportError(CE_Failure, CPLE_AppDefined,
148 : "Dataset %s: indexing variable %s of dimension "
149 : "%s has NaN values",
150 0 : osDSName.c_str(), poVar->GetName().c_str(),
151 : desc.osName.c_str());
152 0 : return std::nullopt;
153 : }
154 12 : const int nSign2 = adfValues[i] > adfValues[i - 1] ? 1 : -1;
155 12 : if (nSign != nSign2)
156 : {
157 0 : ReportError(CE_Failure, CPLE_AppDefined,
158 : "Dataset %s: indexing variable %s of dimension "
159 : "%s is not strictly increasing or decreasing",
160 0 : osDSName.c_str(), poVar->GetName().c_str(),
161 : desc.osName.c_str());
162 0 : return std::nullopt;
163 : }
164 : }
165 13 : desc.nProgressionSign = nSign;
166 : }
167 43 : std::sort(adfValues.begin(), adfValues.end());
168 43 : desc.aaValues.push_back(std::move(adfValues));
169 : }
170 57 : return desc;
171 : }
172 :
173 : /************************************************************************/
174 : /* GDALMdimMosaicAlgorithm::BuildArrayParameters() */
175 : /************************************************************************/
176 :
177 28 : bool GDALMdimMosaicAlgorithm::BuildArrayParameters(
178 : const CPLStringList &aosInputDatasetNames,
179 : std::vector<ArrayParameters> &aoArrayParameters)
180 : {
181 59 : for (const char *pszDatasetName : cpl::Iterate(aosInputDatasetNames))
182 : {
183 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
184 : pszDatasetName, GDAL_OF_MULTIDIM_RASTER | GDAL_OF_VERBOSE_ERROR,
185 102 : CPLStringList(m_inputFormats).List(),
186 102 : CPLStringList(m_openOptions).List(), nullptr));
187 51 : if (!poDS)
188 0 : return false;
189 51 : auto poRG = poDS->GetRootGroup();
190 51 : if (!poRG)
191 : {
192 0 : ReportError(CE_Failure, CPLE_AppDefined,
193 : "Cannot get root group for dataset %s", pszDatasetName);
194 0 : return false;
195 : }
196 51 : std::vector<std::shared_ptr<GDALMDArray>> apoArrays;
197 51 : if (!m_array.empty())
198 : {
199 94 : for (const auto &array : m_array)
200 : {
201 50 : auto poArray = poRG->OpenMDArrayFromFullname(array);
202 50 : if (!poArray)
203 : {
204 2 : ReportError(CE_Failure, CPLE_AppDefined,
205 : "Cannot find array %s in dataset %s",
206 : array.c_str(), pszDatasetName);
207 2 : return false;
208 : }
209 48 : apoArrays.push_back(std::move(poArray));
210 : }
211 : }
212 : else
213 : {
214 15 : for (const std::string &arrayName :
215 35 : poRG->GetMDArrayFullNamesRecursive())
216 : {
217 15 : auto poArray = poRG->OpenMDArrayFromFullname(arrayName);
218 15 : if (!poArray)
219 : {
220 0 : ReportError(CE_Failure, CPLE_AppDefined,
221 : "Cannot open array %s of dataset %s",
222 : arrayName.c_str(), pszDatasetName);
223 0 : return false;
224 : }
225 15 : if (poArray->GetDimensionCount() < 2)
226 10 : continue;
227 5 : m_array.push_back(arrayName);
228 5 : apoArrays.push_back(std::move(poArray));
229 : }
230 5 : if (apoArrays.empty())
231 : {
232 1 : ReportError(
233 : CE_Failure, CPLE_AppDefined,
234 : "No array of dimension count >= 2 found in dataset %s",
235 : pszDatasetName);
236 1 : return false;
237 : }
238 : }
239 :
240 48 : if (aoArrayParameters.empty())
241 25 : aoArrayParameters.resize(apoArrays.size());
242 :
243 83 : for (size_t iArray = 0; iArray < apoArrays.size(); ++iArray)
244 : {
245 52 : auto &arrayParameters = aoArrayParameters[iArray];
246 52 : auto &poArray = apoArrays[iArray];
247 52 : if (poArray->GetDimensionCount() == 0)
248 : {
249 1 : ReportError(CE_Failure, CPLE_AppDefined,
250 : "Array %s of dataset %s has no dimensions",
251 1 : poArray->GetName().c_str(), pszDatasetName);
252 17 : return false;
253 : }
254 :
255 51 : std::vector<SourceShortDimDesc> aoSourceShortDimDesc;
256 : const auto AddToSourceShortDimDesc =
257 57 : [&aoSourceShortDimDesc](const DimensionDesc &dimDesc,
258 57 : uint64_t nSize)
259 : {
260 57 : SourceShortDimDesc sourceDesc;
261 57 : sourceDesc.nSize = nSize;
262 57 : sourceDesc.bIsRegularlySpaced = dimDesc.aaValues.empty();
263 57 : if (sourceDesc.bIsRegularlySpaced)
264 14 : sourceDesc.dfStart = dimDesc.dfStart;
265 : else
266 43 : sourceDesc.dfStart = dimDesc.aaValues[0][0];
267 57 : aoSourceShortDimDesc.push_back(std::move(sourceDesc));
268 57 : };
269 :
270 51 : const auto anBlockSize = poArray->GetBlockSize();
271 51 : CPLAssert(anBlockSize.size() == poArray->GetDimensionCount());
272 :
273 51 : if (!arrayParameters.poFirstSourceArray)
274 : {
275 26 : arrayParameters.poFirstSourceArray = poArray;
276 26 : CPLAssert(arrayParameters.mosaicDimensions.empty());
277 26 : size_t iDim = 0;
278 57 : for (auto &poDim : poArray->GetDimensions())
279 : {
280 66 : auto descOpt = GetDimensionDesc(pszDatasetName, poDim);
281 33 : if (!descOpt.has_value())
282 2 : return false;
283 31 : auto &desc = descOpt.value();
284 31 : AddToSourceShortDimDesc(desc, poDim->GetSize());
285 31 : desc.nBlockSize = anBlockSize[iDim];
286 31 : arrayParameters.mosaicDimensions.push_back(std::move(desc));
287 31 : ++iDim;
288 : }
289 : }
290 : else
291 : {
292 50 : if (poArray->GetDimensionCount() !=
293 25 : arrayParameters.mosaicDimensions.size())
294 : {
295 0 : ReportError(CE_Failure, CPLE_AppDefined,
296 : "Array %s of dataset %s does not have the same "
297 : "number of dimensions as in other datasets",
298 0 : poArray->GetName().c_str(), pszDatasetName);
299 14 : return false;
300 : }
301 25 : if (poArray->GetDataType() !=
302 25 : arrayParameters.poFirstSourceArray->GetDataType())
303 : {
304 1 : ReportError(CE_Failure, CPLE_AppDefined,
305 : "Array %s of dataset %s does not have the same "
306 : "data type as in other datasets",
307 1 : poArray->GetName().c_str(), pszDatasetName);
308 1 : return false;
309 : }
310 : const void *poFirstRawNoData =
311 24 : arrayParameters.poFirstSourceArray->GetRawNoDataValue();
312 24 : const void *poRawNoData = poArray->GetRawNoDataValue();
313 26 : if (!((!poFirstRawNoData && !poRawNoData) ||
314 2 : (poFirstRawNoData && poRawNoData &&
315 1 : memcmp(poFirstRawNoData, poRawNoData,
316 1 : poArray->GetDataType().GetSize()) == 0)))
317 : {
318 3 : ReportError(CE_Failure, CPLE_AppDefined,
319 : "Array %s of dataset %s does not have the same "
320 : "nodata value as in other datasets",
321 3 : poArray->GetName().c_str(), pszDatasetName);
322 3 : return false;
323 : }
324 : const std::vector<std::shared_ptr<GDALDimension>> apoDims =
325 21 : poArray->GetDimensions();
326 38 : for (size_t iDim = 0;
327 38 : iDim < arrayParameters.mosaicDimensions.size(); ++iDim)
328 : {
329 : DimensionDesc &desc =
330 27 : arrayParameters.mosaicDimensions[iDim];
331 27 : auto &poDim = apoDims[iDim];
332 27 : if (poDim->GetName() != desc.osName)
333 : {
334 1 : ReportError(
335 : CE_Failure, CPLE_AppDefined,
336 : "Dimension %d of array %s of dataset %s does "
337 : "not have the same name as in other datasets",
338 1 : static_cast<int>(iDim), poArray->GetName().c_str(),
339 : pszDatasetName);
340 10 : return false;
341 : }
342 26 : if (desc.nBlockSize != anBlockSize[iDim])
343 1 : desc.nBlockSize = 0;
344 :
345 : auto descThisDatasetOpt =
346 52 : GetDimensionDesc(pszDatasetName, poDim);
347 26 : if (!descThisDatasetOpt.has_value())
348 0 : return false;
349 26 : auto &descThisDataset = descThisDatasetOpt.value();
350 26 : AddToSourceShortDimDesc(descThisDataset, poDim->GetSize());
351 26 : if (desc.aaValues.empty())
352 : {
353 6 : if (!descThisDataset.aaValues.empty())
354 : {
355 2 : ReportError(
356 : CE_Failure, CPLE_AppDefined,
357 : "Dimension %s of array %s of dataset %s "
358 : "has irregularly-spaced values, contrary "
359 : "to other datasets",
360 1 : poDim->GetName().c_str(),
361 1 : poArray->GetName().c_str(), pszDatasetName);
362 1 : return false;
363 : }
364 5 : if (std::fabs(descThisDataset.dfIncrement -
365 5 : desc.dfIncrement) >
366 5 : 1e-10 * std::fabs(desc.dfIncrement))
367 : {
368 2 : ReportError(
369 : CE_Failure, CPLE_AppDefined,
370 : "Dimension %s of array %s of dataset %s is "
371 : "indexed by a variable with spacing %g, "
372 : "whereas it is %g in other datasets",
373 1 : poDim->GetName().c_str(),
374 1 : poArray->GetName().c_str(), pszDatasetName,
375 : descThisDataset.dfIncrement, desc.dfIncrement);
376 1 : return false;
377 : }
378 4 : const double dfPos =
379 4 : (descThisDataset.dfStart - desc.dfStart) /
380 4 : desc.dfIncrement;
381 4 : if (std::fabs(std::round(dfPos) - dfPos) > 1e-3)
382 : {
383 2 : ReportError(CE_Failure, CPLE_AppDefined,
384 : "Dimension %s of array %s of dataset "
385 : "%s is indexed "
386 : "by a variable whose start value is "
387 : "not aligned "
388 : "with the one of other datasets",
389 1 : poDim->GetName().c_str(),
390 1 : poArray->GetName().c_str(),
391 : pszDatasetName);
392 1 : return false;
393 : }
394 3 : desc.dfStart =
395 3 : std::min(desc.dfStart, descThisDataset.dfStart);
396 : const double dfEnd = std::max(
397 6 : desc.dfStart + static_cast<double>(desc.nSize) *
398 3 : desc.dfIncrement,
399 6 : descThisDataset.dfStart +
400 3 : static_cast<double>(descThisDataset.nSize) *
401 3 : descThisDataset.dfIncrement);
402 3 : const double dfSize =
403 3 : (dfEnd - desc.dfStart) / desc.dfIncrement;
404 3 : constexpr double MAX_INTEGER_REPRESENTABLE =
405 : static_cast<double>(1ULL << 53);
406 3 : if (dfSize > MAX_INTEGER_REPRESENTABLE)
407 : {
408 0 : ReportError(
409 : CE_Failure, CPLE_AppDefined,
410 : "Dimension %s of array %s of dataset %s "
411 : "would be too large if merged.",
412 0 : poDim->GetName().c_str(),
413 0 : poArray->GetName().c_str(), pszDatasetName);
414 0 : return false;
415 : }
416 3 : desc.nSize = static_cast<uint64_t>(dfSize + 0.5);
417 : }
418 : else
419 : {
420 20 : if (descThisDataset.aaValues.empty())
421 : {
422 2 : ReportError(
423 : CE_Failure, CPLE_AppDefined,
424 : "Dimension %s of array %s of dataset %s "
425 : "has regularly spaced labels, contrary to "
426 : "other datasets",
427 1 : poDim->GetName().c_str(),
428 1 : poArray->GetName().c_str(), pszDatasetName);
429 1 : return false;
430 : }
431 19 : if (descThisDataset.nProgressionSign !=
432 19 : desc.nProgressionSign)
433 : {
434 2 : ReportError(
435 : CE_Failure, CPLE_AppDefined,
436 : "Dataset %s: values in indexing variable %s of "
437 : "dimension %s must be either increasing or "
438 : "decreasing in all input datasets.",
439 : pszDatasetName,
440 2 : poDim->GetIndexingVariable()->GetName().c_str(),
441 : desc.osName.c_str());
442 1 : return false;
443 : }
444 18 : CPLAssert(descThisDataset.aaValues.size() == 1);
445 36 : if (descThisDataset.aaValues[0][0] <
446 18 : desc.aaValues[0][0])
447 : {
448 6 : if (descThisDataset.aaValues[0].back() >=
449 3 : desc.aaValues[0][0])
450 : {
451 2 : ReportError(
452 : CE_Failure, CPLE_AppDefined,
453 : "Dataset %s: values in indexing variable "
454 : "%s of "
455 : "dimension %s are not the same as in other "
456 : "datasets",
457 : pszDatasetName,
458 2 : poDim->GetIndexingVariable()
459 1 : ->GetName()
460 : .c_str(),
461 : desc.osName.c_str());
462 1 : return false;
463 : }
464 : desc.aaValues.insert(
465 2 : desc.aaValues.begin(),
466 4 : std::move(descThisDataset.aaValues[0]));
467 : }
468 : else
469 : {
470 20 : for (size_t i = 0; i < desc.aaValues.size(); ++i)
471 : {
472 32 : if (descThisDataset.aaValues[0][0] ==
473 16 : desc.aaValues[i][0])
474 : {
475 5 : if (descThisDataset.aaValues[0] !=
476 5 : desc.aaValues[i])
477 : {
478 2 : ReportError(
479 : CE_Failure, CPLE_AppDefined,
480 : "Dataset %s: values in indexing "
481 : "variable %s of dimension %s are "
482 : "not "
483 : "the same as in other datasets",
484 : pszDatasetName,
485 2 : poDim->GetIndexingVariable()
486 1 : ->GetName()
487 : .c_str(),
488 : desc.osName.c_str());
489 1 : return false;
490 : }
491 : }
492 11 : else if (descThisDataset.aaValues[0][0] >
493 24 : desc.aaValues[i][0] &&
494 11 : (i + 1 == desc.aaValues.size() ||
495 4 : descThisDataset.aaValues[0][0] <
496 2 : desc.aaValues[i + 1][0]))
497 : {
498 10 : if (descThisDataset.aaValues[0][0] <=
499 10 : desc.aaValues[i].back())
500 : {
501 2 : ReportError(
502 : CE_Failure, CPLE_AppDefined,
503 : "Dataset %s: values in indexing "
504 : "variable %s of dimension %s are "
505 : "overlapping with the ones of "
506 : "other "
507 : "datasets",
508 : pszDatasetName,
509 2 : poDim->GetIndexingVariable()
510 1 : ->GetName()
511 : .c_str(),
512 : desc.osName.c_str());
513 1 : return false;
514 : }
515 10 : if (i + 1 < desc.aaValues.size() &&
516 2 : descThisDataset.aaValues[0].back() >=
517 1 : desc.aaValues[i + 1][0])
518 : {
519 2 : ReportError(
520 : CE_Failure, CPLE_AppDefined,
521 : "Dataset %s: values in indexing "
522 : "variable %s of dimension %s are "
523 : "overlapping with the ones of "
524 : "other "
525 : "datasets",
526 : pszDatasetName,
527 2 : poDim->GetIndexingVariable()
528 1 : ->GetName()
529 : .c_str(),
530 : desc.osName.c_str());
531 1 : return false;
532 : }
533 : desc.aaValues.insert(
534 8 : desc.aaValues.begin() + i + 1,
535 16 : std::move(descThisDataset.aaValues[0]));
536 8 : break;
537 : }
538 : }
539 : }
540 : }
541 : }
542 : }
543 :
544 35 : arrayParameters.aaoSourceShortDimDesc.push_back(
545 35 : std::move(aoSourceShortDimDesc));
546 : }
547 : }
548 :
549 8 : return true;
550 : }
551 :
552 : /************************************************************************/
553 : /* GDALMdimMosaicAlgorithm::GetInputDatasetNames() */
554 : /************************************************************************/
555 :
556 28 : bool GDALMdimMosaicAlgorithm::GetInputDatasetNames(
557 : GDALProgressFunc pfnProgress, void *pProgressData,
558 : CPLStringList &aosInputDatasetNames) const
559 : {
560 73 : for (auto &ds : m_inputDatasets)
561 : {
562 45 : if (ds.GetName()[0] == '@')
563 : {
564 : auto f = VSIVirtualHandleUniquePtr(
565 1 : VSIFOpenL(ds.GetName().c_str() + 1, "r"));
566 1 : if (!f)
567 : {
568 0 : ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s",
569 0 : ds.GetName().c_str() + 1);
570 0 : return false;
571 : }
572 3 : while (const char *filename = CPLReadLineL(f.get()))
573 : {
574 2 : aosInputDatasetNames.push_back(filename);
575 2 : }
576 : }
577 44 : else if (ds.GetName().find_first_of("*?[") != std::string::npos)
578 : {
579 5 : CPLStringList aosMatches(VSIGlob(ds.GetName().c_str(), nullptr,
580 10 : pfnProgress, pProgressData));
581 16 : for (const char *pszStr : aosMatches)
582 : {
583 11 : aosInputDatasetNames.push_back(pszStr);
584 : }
585 : }
586 : else
587 : {
588 78 : std::string osDatasetName = ds.GetName();
589 39 : if (!GetReferencePathForRelativePaths().empty())
590 : {
591 0 : osDatasetName = GDALDataset::BuildFilename(
592 : osDatasetName.c_str(),
593 0 : GetReferencePathForRelativePaths().c_str(), true);
594 : }
595 39 : aosInputDatasetNames.push_back(osDatasetName.c_str());
596 : }
597 : }
598 28 : return true;
599 : }
600 :
601 : /************************************************************************/
602 : /* GDALMdimMosaicAlgorithm::RunImpl() */
603 : /************************************************************************/
604 :
605 28 : bool GDALMdimMosaicAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
606 : void *pProgressData)
607 : {
608 28 : CPLAssert(!m_outputDataset.GetDatasetRef());
609 :
610 28 : if (m_outputFormat.empty())
611 : {
612 : const auto aosFormats =
613 : CPLStringList(GDALGetOutputDriversForDatasetName(
614 3 : m_outputDataset.GetName().c_str(), GDAL_OF_MULTIDIM_RASTER,
615 : /* bSingleMatch = */ true,
616 3 : /* bWarn = */ true));
617 3 : if (aosFormats.size() != 1)
618 : {
619 0 : ReportError(CE_Failure, CPLE_AppDefined,
620 : "Cannot guess driver for %s",
621 0 : m_outputDataset.GetName().c_str());
622 0 : return false;
623 : }
624 3 : m_outputFormat = aosFormats[0];
625 : }
626 : auto poOutDrv =
627 28 : GetGDALDriverManager()->GetDriverByName(m_outputFormat.c_str());
628 28 : if (!poOutDrv)
629 : {
630 0 : ReportError(CE_Failure, CPLE_AppDefined, "Driver %s does not exist",
631 : m_outputFormat.c_str());
632 0 : return false;
633 : }
634 :
635 28 : const bool bIsVRT = EQUAL(m_outputFormat.c_str(), "VRT");
636 :
637 56 : CPLStringList aosInputDatasetNames;
638 28 : const double dfIntermediatePercentage = bIsVRT ? 1.0 : 0.1;
639 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)> pScaledData(
640 : GDALCreateScaledProgress(0.0, dfIntermediatePercentage, pfnProgress,
641 : pProgressData),
642 56 : GDALDestroyScaledProgress);
643 28 : if (!GetInputDatasetNames(GDALScaledProgress, pScaledData.get(),
644 : aosInputDatasetNames))
645 0 : return false;
646 :
647 53 : for (std::string &name : m_array)
648 : {
649 25 : if (!name.empty() && name[0] != '/')
650 25 : name = "/" + name;
651 : }
652 :
653 56 : std::vector<ArrayParameters> aoArrayParameters;
654 28 : if (!BuildArrayParameters(aosInputDatasetNames, aoArrayParameters))
655 : {
656 20 : return false;
657 : }
658 :
659 16 : auto poVRTDS = VRTDataset::CreateVRTMultiDimensional("", nullptr, nullptr);
660 8 : CPLAssert(poVRTDS);
661 :
662 16 : auto poDstGroup = poVRTDS->GetRootVRTGroup();
663 8 : CPLAssert(poDstGroup);
664 :
665 : std::map<std::string, std::shared_ptr<GDALDimension>>
666 16 : oMapAlreadyCreatedDims;
667 :
668 : // Iterate over arrays
669 18 : for (auto &arrayParameters : aoArrayParameters)
670 : {
671 10 : auto &poFirstSourceArray = arrayParameters.poFirstSourceArray;
672 10 : CPLAssert(poFirstSourceArray);
673 10 : auto &aaoSourceShortDimDesc = arrayParameters.aaoSourceShortDimDesc;
674 10 : CPLAssert(aaoSourceShortDimDesc.size() ==
675 : static_cast<size_t>(aosInputDatasetNames.size()));
676 :
677 : // Create mosaic array dimensions
678 10 : std::vector<std::shared_ptr<GDALDimension>> apoDstDims;
679 27 : for (auto &desc : arrayParameters.mosaicDimensions)
680 : {
681 17 : auto oIterCreatedDims = oMapAlreadyCreatedDims.find(desc.osName);
682 17 : if (oIterCreatedDims != oMapAlreadyCreatedDims.end())
683 : {
684 4 : apoDstDims.push_back(oIterCreatedDims->second);
685 : }
686 : else
687 : {
688 13 : uint64_t nDimSize64 = desc.nSize;
689 13 : if (!desc.aaValues.empty())
690 : {
691 8 : nDimSize64 = 0;
692 23 : for (const auto &aValues : desc.aaValues)
693 15 : nDimSize64 += aValues.size();
694 : }
695 : auto dstDim = poDstGroup->CreateDimension(
696 13 : desc.osName, desc.osType, desc.osDirection, nDimSize64);
697 13 : if (!dstDim)
698 0 : return false;
699 :
700 : auto var = poDstGroup->CreateVRTMDArray(
701 13 : desc.osName, {dstDim},
702 52 : GDALExtendedDataType::Create(GDT_Float64));
703 13 : if (!var)
704 0 : return false;
705 :
706 31 : for (const auto &attr : desc.attributes)
707 : {
708 36 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
709 : auto dstAttr = var->CreateAttribute(
710 18 : attr->GetName(), attr->GetDimensionsSize(),
711 54 : attr->GetDataType());
712 18 : if (dstAttr)
713 : {
714 18 : auto raw(attr->ReadAsRaw());
715 18 : CPL_IGNORE_RET_VAL(
716 18 : dstAttr->Write(raw.data(), raw.size()));
717 : }
718 : }
719 :
720 13 : if (desc.aaValues.empty())
721 : {
722 : auto poSource =
723 : std::make_unique<VRTMDArraySourceRegularlySpaced>(
724 5 : desc.dfStart, desc.dfIncrement);
725 5 : var->AddSource(std::move(poSource));
726 : }
727 : else
728 : {
729 8 : const size_t nDimSize = static_cast<size_t>(nDimSize64);
730 16 : std::vector<GUInt64> anOffset = {0};
731 16 : std::vector<size_t> anCount = {nDimSize};
732 16 : std::vector<double> adfValues;
733 8 : adfValues.reserve(nDimSize);
734 8 : if (desc.nProgressionSign >= 0)
735 : {
736 23 : for (const auto &aValues : desc.aaValues)
737 15 : adfValues.insert(adfValues.end(), aValues.begin(),
738 30 : aValues.end());
739 : }
740 : else
741 : {
742 0 : for (auto oIter = desc.aaValues.rbegin();
743 0 : oIter != desc.aaValues.rend(); ++oIter)
744 : {
745 0 : adfValues.insert(adfValues.end(), oIter->rbegin(),
746 0 : oIter->rend());
747 : }
748 : }
749 16 : std::vector<GByte> abyValues(nDimSize * sizeof(double));
750 8 : memcpy(abyValues.data(), adfValues.data(),
751 : nDimSize * sizeof(double));
752 : auto poSource =
753 : std::make_unique<VRTMDArraySourceInlinedValues>(
754 0 : var.get(),
755 8 : /* bIsConstantValue = */ false, std::move(anOffset),
756 16 : std::move(anCount), std::move(abyValues));
757 8 : var->AddSource(std::move(poSource));
758 : }
759 13 : dstDim->SetIndexingVariable(std::move(var));
760 13 : oMapAlreadyCreatedDims[dstDim->GetName()] = dstDim;
761 13 : apoDstDims.push_back(std::move(dstDim));
762 : }
763 : }
764 :
765 : // Create mosaic array
766 10 : CPLStringList aosArrayCO;
767 10 : std::string osBlockSize;
768 27 : for (size_t i = 0; i < apoDstDims.size(); ++i)
769 : {
770 17 : if (i > 0)
771 7 : osBlockSize += ',';
772 : osBlockSize +=
773 17 : std::to_string(arrayParameters.mosaicDimensions[i].nBlockSize);
774 : }
775 10 : if (!osBlockSize.empty())
776 10 : aosArrayCO.SetNameValue("BLOCKSIZE", osBlockSize.c_str());
777 :
778 : auto poDstArray = poDstGroup->CreateVRTMDArray(
779 10 : CPLGetFilename(poFirstSourceArray->GetName().c_str()), apoDstDims,
780 30 : poFirstSourceArray->GetDataType(), aosArrayCO);
781 10 : if (!poDstArray)
782 0 : return false;
783 :
784 10 : GUInt64 nCurCost = 0;
785 10 : poDstArray->CopyFromAllExceptValues(poFirstSourceArray.get(), false,
786 : nCurCost, 0, nullptr, nullptr);
787 :
788 : // Add sources to mosaic array
789 30 : for (int iDS = 0; iDS < aosInputDatasetNames.size(); ++iDS)
790 : {
791 20 : const auto &aoSourceShortDimDesc = aaoSourceShortDimDesc[iDS];
792 :
793 20 : const auto dimCount = arrayParameters.mosaicDimensions.size();
794 40 : std::vector<GUInt64> anSrcOffset(dimCount);
795 40 : std::vector<GUInt64> anCount(dimCount);
796 40 : std::vector<GUInt64> anDstOffset;
797 20 : CPLAssert(aoSourceShortDimDesc.size() == dimCount);
798 :
799 53 : for (size_t iDim = 0; iDim < dimCount; ++iDim)
800 : {
801 : const DimensionDesc &desc =
802 33 : arrayParameters.mosaicDimensions[iDim];
803 : const SourceShortDimDesc &sourceDesc =
804 33 : aoSourceShortDimDesc[iDim];
805 33 : if (sourceDesc.bIsRegularlySpaced)
806 : {
807 8 : const double dfPos =
808 8 : (sourceDesc.dfStart - desc.dfStart) / desc.dfIncrement;
809 8 : anDstOffset.push_back(static_cast<uint64_t>(dfPos + 0.5));
810 : }
811 : else
812 : {
813 25 : uint64_t nPos = 0;
814 25 : bool bFound = false;
815 35 : for (size_t i = 0; i < desc.aaValues.size(); ++i)
816 : {
817 35 : if (sourceDesc.dfStart == desc.aaValues[i][0])
818 : {
819 25 : bFound = true;
820 25 : break;
821 : }
822 : else
823 : {
824 10 : nPos += desc.aaValues[i].size();
825 : }
826 : }
827 25 : CPLAssert(bFound);
828 25 : CPL_IGNORE_RET_VAL(bFound);
829 :
830 25 : anDstOffset.push_back(nPos);
831 : }
832 :
833 33 : anCount[iDim] = sourceDesc.nSize;
834 : }
835 :
836 40 : std::vector<GUInt64> anStep(dimCount, 1);
837 40 : std::vector<int> anTransposedAxis;
838 : auto poSource = std::make_unique<VRTMDArraySourceFromArray>(
839 40 : poDstArray.get(), false, false, aosInputDatasetNames[iDS],
840 40 : poFirstSourceArray->GetFullName(), std::string(),
841 20 : std::move(anTransposedAxis),
842 20 : std::string(), // viewExpr
843 20 : std::move(anSrcOffset), std::move(anCount), std::move(anStep),
844 40 : std::move(anDstOffset));
845 20 : poDstArray->AddSource(std::move(poSource));
846 : }
847 : }
848 :
849 8 : pScaledData.reset(GDALCreateScaledProgress(dfIntermediatePercentage, 1.0,
850 : pfnProgress, pProgressData));
851 : auto poOutDS = std::unique_ptr<GDALDataset>(
852 8 : poOutDrv->CreateCopy(m_outputDataset.GetName().c_str(), poVRTDS.get(),
853 8 : false, CPLStringList(m_creationOptions).List(),
854 16 : GDALScaledProgress, pScaledData.get()));
855 :
856 8 : if (poOutDS)
857 8 : m_outputDataset.Set(std::move(poOutDS));
858 :
859 8 : return m_outputDataset.GetDatasetRef() != nullptr;
860 : }
861 :
862 : //! @endcond
|
