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