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