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