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