Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: gdalmultidim_bridge_classic.cpp
4 : * Project: GDAL Core
5 : * Purpose: Bridge from GDALMDArray to GDALRasterBand/GDALDataset
6 : * Author: Even Rouault <even.rouault at spatialys.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2019, Even Rouault <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdal_multidim.h"
15 : #include "gdal_pam.h"
16 : #include "gdal_pam_multidim.h"
17 :
18 : #include <algorithm>
19 :
20 : //! @cond Doxygen_Suppress
21 :
22 : /************************************************************************/
23 : /* GDALDatasetFromArray() */
24 : /************************************************************************/
25 :
26 : class GDALDatasetFromArray;
27 :
28 : namespace
29 : {
30 : struct MetadataItem
31 : {
32 : std::shared_ptr<GDALAbstractMDArray> poArray{};
33 : std::string osName{};
34 : std::string osDefinition{};
35 : bool bDefinitionUsesPctForG = false;
36 : };
37 :
38 : struct BandImageryMetadata
39 : {
40 : std::shared_ptr<GDALAbstractMDArray> poCentralWavelengthArray{};
41 : double dfCentralWavelengthToMicrometer = 1.0;
42 : std::shared_ptr<GDALAbstractMDArray> poFWHMArray{};
43 : double dfFWHMToMicrometer = 1.0;
44 : };
45 :
46 : } // namespace
47 :
48 : class GDALRasterBandFromArray final : public GDALPamRasterBand
49 : {
50 : std::vector<GUInt64> m_anOffset{};
51 : std::vector<size_t> m_anCount{};
52 : std::vector<GPtrDiff_t> m_anStride{};
53 :
54 : protected:
55 : CPLErr IReadBlock(int, int, void *) override;
56 : CPLErr IWriteBlock(int, int, void *) override;
57 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
58 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
59 : GDALDataType eBufType, GSpacing nPixelSpaceBuf,
60 : GSpacing nLineSpaceBuf,
61 : GDALRasterIOExtraArg *psExtraArg) override;
62 :
63 : public:
64 : explicit GDALRasterBandFromArray(
65 : GDALDatasetFromArray *poDSIn,
66 : const std::vector<GUInt64> &anOtherDimCoord,
67 : const std::vector<std::vector<MetadataItem>>
68 : &aoBandParameterMetadataItems,
69 : const std::vector<BandImageryMetadata> &aoBandImageryMetadata,
70 : double dfDelay, time_t nStartTime, bool &bHasWarned);
71 :
72 : double GetNoDataValue(int *pbHasNoData) override;
73 : int64_t GetNoDataValueAsInt64(int *pbHasNoData) override;
74 : uint64_t GetNoDataValueAsUInt64(int *pbHasNoData) override;
75 : double GetOffset(int *pbHasOffset) override;
76 : double GetScale(int *pbHasScale) override;
77 : const char *GetUnitType() override;
78 : GDALColorInterp GetColorInterpretation() override;
79 : int GetOverviewCount() override;
80 : GDALRasterBand *GetOverview(int idx) override;
81 : CPLErr AdviseRead(int nXOff, int nYOff, int nXSize, int nYSize,
82 : int nBufXSize, int nBufYSize, GDALDataType eBufType,
83 : CSLConstList papszOptions) override;
84 : };
85 :
86 : class GDALDatasetFromArray final : public GDALPamDataset
87 : {
88 : friend class GDALRasterBandFromArray;
89 :
90 : std::shared_ptr<GDALMDArray> m_poArray;
91 : const size_t m_iXDim;
92 : const size_t m_iYDim;
93 : const CPLStringList m_aosOptions;
94 : GDALGeoTransform m_gt{};
95 : bool m_bHasGT = false;
96 : mutable std::shared_ptr<OGRSpatialReference> m_poSRS{};
97 : GDALMultiDomainMetadata m_oMDD{};
98 : std::string m_osOvrFilename{};
99 : bool m_bOverviewsDiscovered = false;
100 : bool m_bPixelInterleaved = false;
101 : std::vector<std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>>
102 : m_apoOverviews{};
103 : std::vector<GUInt64> m_anOffset{};
104 : std::vector<size_t> m_anCount{};
105 : std::vector<GPtrDiff_t> m_anStride{};
106 :
107 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
108 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
109 : GDALDataType eBufType, int nBandCount,
110 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
111 : GSpacing nLineSpace, GSpacing nBandSpace,
112 : GDALRasterIOExtraArg *psExtraArg) override;
113 :
114 : public:
115 476 : GDALDatasetFromArray(const std::shared_ptr<GDALMDArray> &array,
116 : size_t iXDim, size_t iYDim,
117 : const CPLStringList &aosOptions)
118 476 : : m_poArray(array), m_iXDim(iXDim), m_iYDim(iYDim),
119 476 : m_aosOptions(aosOptions)
120 : {
121 476 : const auto nDimCount = m_poArray->GetDimensionCount();
122 476 : m_anOffset.resize(nDimCount);
123 476 : m_anCount.resize(nDimCount, 1);
124 476 : m_anStride.resize(nDimCount);
125 476 : }
126 :
127 : static std::unique_ptr<GDALDatasetFromArray>
128 : Create(const std::shared_ptr<GDALMDArray> &array, size_t iXDim,
129 : size_t iYDim, const std::shared_ptr<GDALGroup> &poRootGroup,
130 : CSLConstList papszOptions);
131 :
132 : ~GDALDatasetFromArray() override;
133 :
134 673 : CPLErr Close(GDALProgressFunc = nullptr, void * = nullptr) override
135 : {
136 673 : CPLErr eErr = CE_None;
137 673 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
138 : {
139 673 : if (GDALDatasetFromArray::FlushCache(/*bAtClosing=*/true) !=
140 : CE_None)
141 0 : eErr = CE_Failure;
142 673 : m_poArray.reset();
143 : }
144 673 : return eErr;
145 : }
146 :
147 83 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
148 : {
149 83 : gt = m_gt;
150 83 : return m_bHasGT ? CE_None : CE_Failure;
151 : }
152 :
153 93 : const OGRSpatialReference *GetSpatialRef() const override
154 : {
155 93 : if (m_poArray->GetDimensionCount() < 2)
156 3 : return nullptr;
157 90 : m_poSRS = m_poArray->GetSpatialRef();
158 90 : if (m_poSRS)
159 : {
160 46 : m_poSRS.reset(m_poSRS->Clone());
161 92 : auto axisMapping = m_poSRS->GetDataAxisToSRSAxisMapping();
162 138 : for (auto &m : axisMapping)
163 : {
164 92 : if (m == static_cast<int>(m_iXDim) + 1)
165 46 : m = 1;
166 46 : else if (m == static_cast<int>(m_iYDim) + 1)
167 46 : m = 2;
168 : }
169 46 : m_poSRS->SetDataAxisToSRSAxisMapping(axisMapping);
170 : }
171 90 : return m_poSRS.get();
172 : }
173 :
174 9 : CPLErr SetMetadata(CSLConstList papszMetadata,
175 : const char *pszDomain) override
176 : {
177 9 : return m_oMDD.SetMetadata(papszMetadata, pszDomain);
178 : }
179 :
180 184 : CSLConstList GetMetadata(const char *pszDomain) override
181 : {
182 184 : return m_oMDD.GetMetadata(pszDomain);
183 : }
184 :
185 312 : const char *GetMetadataItem(const char *pszName,
186 : const char *pszDomain) override
187 : {
188 312 : const std::string &osFilename = m_poArray->GetFilename();
189 553 : if (!osFilename.empty() && pszName && EQUAL(pszName, "OVERVIEW_FILE") &&
190 553 : pszDomain && EQUAL(pszDomain, "OVERVIEWS"))
191 : {
192 57 : if (m_osOvrFilename.empty())
193 : {
194 : // Legacy strategy (pre GDAL 3.13)
195 96 : std::string osOvrFilename = osFilename;
196 48 : osOvrFilename += '.';
197 556 : for (char ch : m_poArray->GetName())
198 : {
199 508 : if ((ch >= '0' && ch <= '9') || (ch >= 'A' && ch <= 'Z') ||
200 449 : (ch >= 'a' && ch <= 'z') || ch == '_')
201 : {
202 425 : osOvrFilename += ch;
203 : }
204 : else
205 : {
206 83 : osOvrFilename += '_';
207 : }
208 : }
209 48 : osOvrFilename += ".ovr";
210 : VSIStatBufL sStatBuf;
211 48 : if (VSIStatL(osOvrFilename.c_str(), &sStatBuf) == 0)
212 : {
213 0 : m_osOvrFilename = std::move(osOvrFilename);
214 : }
215 : else
216 : {
217 48 : auto poPAM = GDALPamMultiDim::GetPAM(m_poArray);
218 48 : if (!poPAM)
219 6 : poPAM = std::make_shared<GDALPamMultiDim>(osFilename);
220 192 : m_osOvrFilename = poPAM->GetOverviewFilename(
221 144 : m_poArray->GetFullName(), m_poArray->GetContext());
222 : }
223 :
224 48 : if (!m_osOvrFilename.empty())
225 2 : oOvManager.Initialize(this, ":::VIRTUAL:::");
226 : }
227 57 : return !m_osOvrFilename.empty() ? m_osOvrFilename.c_str() : nullptr;
228 : }
229 255 : return m_oMDD.GetMetadataItem(pszName, pszDomain);
230 : }
231 :
232 77 : bool HasRegisteredPAMOverviewFile()
233 : {
234 77 : return GetMetadataItem("OVERVIEW_FILE", "OVERVIEWS") != nullptr;
235 : }
236 :
237 3 : CPLErr IBuildOverviews(const char *pszResampling, int nOverviews,
238 : const int *panOverviewList, int nListBands,
239 : const int *panBandList, GDALProgressFunc pfnProgress,
240 : void *pProgressData,
241 : CSLConstList papszOptions) override
242 : {
243 : // Try the multidimensional array path. Use quiet handler to
244 : // suppress the "not supported" error from the base class stub.
245 3 : bool bNotSupported = false;
246 6 : std::string osErrMsg;
247 3 : CPLErr eSavedClass = CE_None;
248 3 : int nSavedNo = CPLE_None;
249 : {
250 3 : CPLErrorHandlerPusher oQuiet(CPLQuietErrorHandler);
251 3 : CPLErr eErr = m_poArray->BuildOverviews(
252 : pszResampling, nOverviews, panOverviewList, pfnProgress,
253 3 : pProgressData, papszOptions);
254 3 : if (eErr == CE_None)
255 : {
256 1 : m_bOverviewsDiscovered = false;
257 1 : m_apoOverviews.clear();
258 1 : return CE_None;
259 : }
260 2 : nSavedNo = CPLGetLastErrorNo();
261 2 : eSavedClass = CPLGetLastErrorType();
262 2 : osErrMsg = CPLGetLastErrorMsg();
263 2 : bNotSupported = (nSavedNo == CPLE_NotSupported);
264 : }
265 2 : if (!bNotSupported)
266 : {
267 : // Re-emit the error that was suppressed by the quiet handler.
268 0 : CPLError(eSavedClass, nSavedNo, "%s", osErrMsg.c_str());
269 0 : return CE_Failure;
270 : }
271 : // Driver doesn't implement BuildOverviews - fall back to
272 : // default path (e.g. external .ovr file).
273 2 : CPLErrorReset();
274 :
275 2 : if (!HasRegisteredPAMOverviewFile())
276 : {
277 2 : const std::string &osFilename = m_poArray->GetFilename();
278 2 : if (osFilename.empty())
279 : {
280 0 : CPLError(CE_Failure, CPLE_AppDefined,
281 : "No filename associated with array %s",
282 0 : m_poArray->GetFullName().c_str());
283 0 : return CE_Failure;
284 : }
285 2 : auto poPAM = GDALPamMultiDim::GetPAM(m_poArray);
286 2 : if (!poPAM)
287 0 : poPAM = std::make_shared<GDALPamMultiDim>(osFilename);
288 8 : m_osOvrFilename = poPAM->GenerateOverviewFilename(
289 6 : m_poArray->GetFullName(), m_poArray->GetContext());
290 2 : if (m_osOvrFilename.empty())
291 0 : return CE_Failure;
292 2 : oOvManager.Initialize(this, ":::VIRTUAL:::");
293 : }
294 :
295 2 : return GDALDataset::IBuildOverviews(
296 : pszResampling, nOverviews, panOverviewList, nListBands, panBandList,
297 2 : pfnProgress, pProgressData, papszOptions);
298 : }
299 :
300 70 : void DiscoverOverviews()
301 : {
302 70 : if (!m_bOverviewsDiscovered)
303 : {
304 24 : m_bOverviewsDiscovered = true;
305 24 : if (const int nOverviews = m_poArray->GetOverviewCount())
306 : {
307 14 : if (auto poRootGroup = m_poArray->GetRootGroup())
308 : {
309 7 : const size_t nDims = m_poArray->GetDimensionCount();
310 14 : CPLStringList aosOptions(m_aosOptions);
311 7 : aosOptions.SetNameValue("LOAD_PAM", "NO");
312 18 : for (int iOvr = 0; iOvr < nOverviews; ++iOvr)
313 : {
314 22 : if (auto poOvrArray = m_poArray->GetOverview(iOvr))
315 : {
316 22 : if (poOvrArray->GetDimensionCount() == nDims &&
317 11 : poOvrArray->GetDataType() ==
318 11 : m_poArray->GetDataType())
319 : {
320 : auto poOvrDS =
321 11 : Create(poOvrArray, m_iXDim, m_iYDim,
322 22 : poRootGroup, aosOptions);
323 11 : if (poOvrDS)
324 : {
325 11 : m_apoOverviews.push_back(
326 22 : std::unique_ptr<
327 : GDALDataset,
328 : GDALDatasetUniquePtrReleaser>(
329 11 : poOvrDS.release()));
330 : }
331 : }
332 : }
333 : }
334 : }
335 : }
336 : }
337 70 : }
338 : };
339 :
340 952 : GDALDatasetFromArray::~GDALDatasetFromArray()
341 : {
342 476 : GDALDatasetFromArray::Close();
343 952 : }
344 :
345 : /************************************************************************/
346 : /* GDALRasterBandFromArray() */
347 : /************************************************************************/
348 :
349 565 : GDALRasterBandFromArray::GDALRasterBandFromArray(
350 : GDALDatasetFromArray *poDSIn, const std::vector<GUInt64> &anOtherDimCoord,
351 : const std::vector<std::vector<MetadataItem>> &aoBandParameterMetadataItems,
352 : const std::vector<BandImageryMetadata> &aoBandImageryMetadata,
353 565 : double dfDelay, time_t nStartTime, bool &bHasWarned)
354 : {
355 565 : const auto &poArray(poDSIn->m_poArray);
356 565 : const auto &dims(poArray->GetDimensions());
357 565 : const auto nDimCount(dims.size());
358 1130 : const auto blockSize(poArray->GetBlockSize());
359 :
360 546 : nBlockYSize = (nDimCount >= 2 && blockSize[poDSIn->m_iYDim])
361 1111 : ? static_cast<int>(std::min(static_cast<GUInt64>(INT_MAX),
362 387 : blockSize[poDSIn->m_iYDim]))
363 : : 1;
364 565 : nBlockXSize = blockSize[poDSIn->m_iXDim]
365 404 : ? static_cast<int>(std::min(static_cast<GUInt64>(INT_MAX),
366 404 : blockSize[poDSIn->m_iXDim]))
367 565 : : poDSIn->GetRasterXSize();
368 :
369 565 : eDataType = poArray->GetDataType().GetNumericDataType();
370 565 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
371 :
372 565 : if (nDTSize > 0)
373 : {
374 : // If the above computed block size exceeds INT_MAX or 1/100th of the
375 : // maximum allowed size for the block cache, divide its shape by two,
376 : // along the largest dimension. Only do that while there are at least
377 : // one dimension with 2 pixels.
378 12 : while (
379 1109 : (nBlockXSize >= 2 || nBlockYSize >= 2) &&
380 532 : (nBlockXSize > INT_MAX / nBlockYSize / nDTSize ||
381 529 : (nBlockXSize > GDALGetCacheMax64() / 100 / nBlockYSize / nDTSize)))
382 : {
383 12 : if (nBlockXSize > nBlockYSize)
384 12 : nBlockXSize /= 2;
385 : else
386 0 : nBlockYSize /= 2;
387 : }
388 : }
389 :
390 565 : eAccess = poDSIn->eAccess;
391 565 : m_anOffset.resize(nDimCount);
392 565 : m_anCount.resize(nDimCount, 1);
393 565 : m_anStride.resize(nDimCount);
394 1827 : for (size_t i = 0, j = 0; i < nDimCount; ++i)
395 : {
396 1262 : if (i != poDSIn->m_iXDim && !(nDimCount >= 2 && i == poDSIn->m_iYDim))
397 : {
398 302 : std::string dimName(dims[i]->GetName());
399 151 : GUInt64 nIndex = anOtherDimCoord[j];
400 : // Detect subset_{orig_dim_name}_{start}_{incr}_{size} names of
401 : // subsetted dimensions as generated by GetView()
402 151 : if (STARTS_WITH(dimName.c_str(), "subset_"))
403 : {
404 : CPLStringList aosTokens(
405 12 : CSLTokenizeString2(dimName.c_str(), "_", 0));
406 6 : if (aosTokens.size() == 5)
407 : {
408 6 : dimName = aosTokens[1];
409 18 : const auto nStartDim = static_cast<GUInt64>(CPLScanUIntBig(
410 6 : aosTokens[2], static_cast<int>(strlen(aosTokens[2]))));
411 6 : const auto nIncrDim = CPLAtoGIntBig(aosTokens[3]);
412 6 : nIndex = nIncrDim > 0 ? nStartDim + nIndex * nIncrDim
413 0 : : nStartDim - (nIndex * -nIncrDim);
414 : }
415 : }
416 151 : if (nDimCount != 3 || dimName != "Band")
417 : {
418 90 : SetMetadataItem(
419 : CPLSPrintf("DIM_%s_INDEX", dimName.c_str()),
420 : CPLSPrintf(CPL_FRMT_GUIB, static_cast<GUIntBig>(nIndex)));
421 : }
422 :
423 151 : auto indexingVar = dims[i]->GetIndexingVariable();
424 :
425 : // If the indexing variable is also listed in band parameter arrays,
426 : // then don't use our default formatting
427 151 : if (indexingVar)
428 : {
429 49 : for (const auto &oItem : aoBandParameterMetadataItems[j])
430 : {
431 14 : if (oItem.poArray->GetFullName() ==
432 14 : indexingVar->GetFullName())
433 : {
434 12 : indexingVar.reset();
435 12 : break;
436 : }
437 : }
438 : }
439 :
440 186 : if (indexingVar && indexingVar->GetDimensionCount() == 1 &&
441 35 : indexingVar->GetDimensions()[0]->GetSize() ==
442 35 : dims[i]->GetSize())
443 : {
444 35 : if (dfDelay >= 0 && time(nullptr) - nStartTime > dfDelay)
445 : {
446 0 : if (!bHasWarned)
447 : {
448 0 : CPLError(
449 : CE_Warning, CPLE_AppDefined,
450 : "Maximum delay to load band metadata from "
451 : "dimension indexing variables has expired. "
452 : "Increase the value of the "
453 : "LOAD_EXTRA_DIM_METADATA_DELAY "
454 : "option of GDALMDArray::AsClassicDataset() "
455 : "(also accessible as the "
456 : "GDAL_LOAD_EXTRA_DIM_METADATA_DELAY "
457 : "configuration option), "
458 : "or set it to 'unlimited' for unlimited delay. ");
459 0 : bHasWarned = true;
460 : }
461 : }
462 : else
463 : {
464 35 : size_t nCount = 1;
465 35 : const auto &dt(indexingVar->GetDataType());
466 70 : std::vector<GByte> abyTmp(dt.GetSize());
467 70 : if (indexingVar->Read(&(anOtherDimCoord[j]), &nCount,
468 35 : nullptr, nullptr, dt, &abyTmp[0]))
469 : {
470 35 : char *pszTmp = nullptr;
471 35 : GDALExtendedDataType::CopyValue(
472 35 : &abyTmp[0], dt, &pszTmp,
473 70 : GDALExtendedDataType::CreateString());
474 35 : dt.FreeDynamicMemory(abyTmp.data());
475 35 : if (pszTmp)
476 : {
477 35 : SetMetadataItem(
478 : CPLSPrintf("DIM_%s_VALUE", dimName.c_str()),
479 : pszTmp);
480 35 : CPLFree(pszTmp);
481 : }
482 :
483 35 : const auto &unit(indexingVar->GetUnit());
484 35 : if (!unit.empty())
485 : {
486 12 : SetMetadataItem(
487 : CPLSPrintf("DIM_%s_UNIT", dimName.c_str()),
488 : unit.c_str());
489 : }
490 : }
491 : }
492 : }
493 :
494 169 : for (const auto &oItem : aoBandParameterMetadataItems[j])
495 : {
496 36 : CPLString osVal;
497 :
498 18 : size_t nCount = 1;
499 18 : const auto &dt(oItem.poArray->GetDataType());
500 18 : if (oItem.bDefinitionUsesPctForG)
501 : {
502 : // There is one and only one %[x][.y]f|g in osDefinition
503 16 : std::vector<GByte> abyTmp(dt.GetSize());
504 16 : if (oItem.poArray->Read(&(anOtherDimCoord[j]), &nCount,
505 8 : nullptr, nullptr, dt, &abyTmp[0]))
506 : {
507 8 : double dfVal = 0;
508 8 : GDALExtendedDataType::CopyValue(
509 8 : &abyTmp[0], dt, &dfVal,
510 16 : GDALExtendedDataType::Create(GDT_Float64));
511 8 : osVal.Printf(oItem.osDefinition.c_str(), dfVal);
512 8 : dt.FreeDynamicMemory(abyTmp.data());
513 : }
514 : }
515 : else
516 : {
517 : // There should be zero or one %s in osDefinition
518 10 : char *pszValue = nullptr;
519 10 : if (dt.GetClass() == GEDTC_STRING)
520 : {
521 4 : CPL_IGNORE_RET_VAL(oItem.poArray->Read(
522 2 : &(anOtherDimCoord[j]), &nCount, nullptr, nullptr,
523 2 : dt, &pszValue));
524 : }
525 : else
526 : {
527 16 : std::vector<GByte> abyTmp(dt.GetSize());
528 16 : if (oItem.poArray->Read(&(anOtherDimCoord[j]), &nCount,
529 : nullptr, nullptr, dt,
530 8 : &abyTmp[0]))
531 : {
532 8 : GDALExtendedDataType::CopyValue(
533 8 : &abyTmp[0], dt, &pszValue,
534 16 : GDALExtendedDataType::CreateString());
535 : }
536 : }
537 :
538 10 : if (pszValue)
539 : {
540 10 : osVal.Printf(oItem.osDefinition.c_str(), pszValue);
541 10 : CPLFree(pszValue);
542 : }
543 : }
544 18 : if (!osVal.empty())
545 18 : SetMetadataItem(oItem.osName.c_str(), osVal);
546 : }
547 :
548 151 : if (aoBandImageryMetadata[j].poCentralWavelengthArray)
549 : {
550 : auto &poCentralWavelengthArray =
551 4 : aoBandImageryMetadata[j].poCentralWavelengthArray;
552 4 : size_t nCount = 1;
553 4 : const auto &dt(poCentralWavelengthArray->GetDataType());
554 8 : std::vector<GByte> abyTmp(dt.GetSize());
555 8 : if (poCentralWavelengthArray->Read(&(anOtherDimCoord[j]),
556 : &nCount, nullptr, nullptr,
557 4 : dt, &abyTmp[0]))
558 : {
559 4 : double dfVal = 0;
560 4 : GDALExtendedDataType::CopyValue(
561 4 : &abyTmp[0], dt, &dfVal,
562 8 : GDALExtendedDataType::Create(GDT_Float64));
563 4 : dt.FreeDynamicMemory(abyTmp.data());
564 4 : SetMetadataItem(
565 : "CENTRAL_WAVELENGTH_UM",
566 : CPLSPrintf(
567 4 : "%g", dfVal * aoBandImageryMetadata[j]
568 4 : .dfCentralWavelengthToMicrometer),
569 : "IMAGERY");
570 : }
571 : }
572 :
573 151 : if (aoBandImageryMetadata[j].poFWHMArray)
574 : {
575 2 : auto &poFWHMArray = aoBandImageryMetadata[j].poFWHMArray;
576 2 : size_t nCount = 1;
577 2 : const auto &dt(poFWHMArray->GetDataType());
578 4 : std::vector<GByte> abyTmp(dt.GetSize());
579 4 : if (poFWHMArray->Read(&(anOtherDimCoord[j]), &nCount, nullptr,
580 2 : nullptr, dt, &abyTmp[0]))
581 : {
582 2 : double dfVal = 0;
583 2 : GDALExtendedDataType::CopyValue(
584 2 : &abyTmp[0], dt, &dfVal,
585 4 : GDALExtendedDataType::Create(GDT_Float64));
586 2 : dt.FreeDynamicMemory(abyTmp.data());
587 2 : SetMetadataItem(
588 : "FWHM_UM",
589 2 : CPLSPrintf("%g", dfVal * aoBandImageryMetadata[j]
590 2 : .dfFWHMToMicrometer),
591 : "IMAGERY");
592 : }
593 : }
594 :
595 151 : m_anOffset[i] = anOtherDimCoord[j];
596 151 : j++;
597 : }
598 : }
599 565 : }
600 :
601 : /************************************************************************/
602 : /* GetNoDataValue() */
603 : /************************************************************************/
604 :
605 161 : double GDALRasterBandFromArray::GetNoDataValue(int *pbHasNoData)
606 : {
607 161 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
608 161 : const auto &poArray(l_poDS->m_poArray);
609 161 : bool bHasNodata = false;
610 161 : const auto res = poArray->GetNoDataValueAsDouble(&bHasNodata);
611 161 : if (pbHasNoData)
612 145 : *pbHasNoData = bHasNodata;
613 161 : return res;
614 : }
615 :
616 : /************************************************************************/
617 : /* GetNoDataValueAsInt64() */
618 : /************************************************************************/
619 :
620 1 : int64_t GDALRasterBandFromArray::GetNoDataValueAsInt64(int *pbHasNoData)
621 : {
622 1 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
623 1 : const auto &poArray(l_poDS->m_poArray);
624 1 : bool bHasNodata = false;
625 1 : const auto nodata = poArray->GetNoDataValueAsInt64(&bHasNodata);
626 1 : if (pbHasNoData)
627 1 : *pbHasNoData = bHasNodata;
628 1 : return nodata;
629 : }
630 :
631 : /************************************************************************/
632 : /* GetNoDataValueAsUInt64() */
633 : /************************************************************************/
634 :
635 1 : uint64_t GDALRasterBandFromArray::GetNoDataValueAsUInt64(int *pbHasNoData)
636 : {
637 1 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
638 1 : const auto &poArray(l_poDS->m_poArray);
639 1 : bool bHasNodata = false;
640 1 : const auto nodata = poArray->GetNoDataValueAsUInt64(&bHasNodata);
641 1 : if (pbHasNoData)
642 1 : *pbHasNoData = bHasNodata;
643 1 : return nodata;
644 : }
645 :
646 : /************************************************************************/
647 : /* GetOffset() */
648 : /************************************************************************/
649 :
650 42 : double GDALRasterBandFromArray::GetOffset(int *pbHasOffset)
651 : {
652 42 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
653 42 : const auto &poArray(l_poDS->m_poArray);
654 42 : bool bHasValue = false;
655 42 : double dfRes = poArray->GetOffset(&bHasValue);
656 42 : if (pbHasOffset)
657 23 : *pbHasOffset = bHasValue;
658 42 : return dfRes;
659 : }
660 :
661 : /************************************************************************/
662 : /* GetUnitType() */
663 : /************************************************************************/
664 :
665 52 : const char *GDALRasterBandFromArray::GetUnitType()
666 : {
667 52 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
668 52 : const auto &poArray(l_poDS->m_poArray);
669 52 : return poArray->GetUnit().c_str();
670 : }
671 :
672 : /************************************************************************/
673 : /* GetScale() */
674 : /************************************************************************/
675 :
676 40 : double GDALRasterBandFromArray::GetScale(int *pbHasScale)
677 : {
678 40 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
679 40 : const auto &poArray(l_poDS->m_poArray);
680 40 : bool bHasValue = false;
681 40 : double dfRes = poArray->GetScale(&bHasValue);
682 40 : if (pbHasScale)
683 21 : *pbHasScale = bHasValue;
684 40 : return dfRes;
685 : }
686 :
687 : /************************************************************************/
688 : /* IReadBlock() */
689 : /************************************************************************/
690 :
691 102 : CPLErr GDALRasterBandFromArray::IReadBlock(int nBlockXOff, int nBlockYOff,
692 : void *pImage)
693 : {
694 102 : const int nDTSize(GDALGetDataTypeSizeBytes(eDataType));
695 102 : const int nXOff = nBlockXOff * nBlockXSize;
696 102 : const int nYOff = nBlockYOff * nBlockYSize;
697 102 : const int nReqXSize = std::min(nRasterXSize - nXOff, nBlockXSize);
698 102 : const int nReqYSize = std::min(nRasterYSize - nYOff, nBlockYSize);
699 : GDALRasterIOExtraArg sExtraArg;
700 102 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
701 204 : return IRasterIO(GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, pImage,
702 : nReqXSize, nReqYSize, eDataType, nDTSize,
703 204 : static_cast<GSpacing>(nDTSize) * nBlockXSize, &sExtraArg);
704 : }
705 :
706 : /************************************************************************/
707 : /* IWriteBlock() */
708 : /************************************************************************/
709 :
710 2 : CPLErr GDALRasterBandFromArray::IWriteBlock(int nBlockXOff, int nBlockYOff,
711 : void *pImage)
712 : {
713 2 : const int nDTSize(GDALGetDataTypeSizeBytes(eDataType));
714 2 : const int nXOff = nBlockXOff * nBlockXSize;
715 2 : const int nYOff = nBlockYOff * nBlockYSize;
716 2 : const int nReqXSize = std::min(nRasterXSize - nXOff, nBlockXSize);
717 2 : const int nReqYSize = std::min(nRasterYSize - nYOff, nBlockYSize);
718 : GDALRasterIOExtraArg sExtraArg;
719 2 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
720 4 : return IRasterIO(GF_Write, nXOff, nYOff, nReqXSize, nReqYSize, pImage,
721 : nReqXSize, nReqYSize, eDataType, nDTSize,
722 4 : static_cast<GSpacing>(nDTSize) * nBlockXSize, &sExtraArg);
723 : }
724 :
725 : /************************************************************************/
726 : /* AdviseRead() */
727 : /************************************************************************/
728 :
729 45 : CPLErr GDALRasterBandFromArray::AdviseRead(int nXOff, int nYOff, int nXSize,
730 : int nYSize, int nBufXSize,
731 : int nBufYSize,
732 : GDALDataType /*eBufType*/,
733 : CSLConstList papszOptions)
734 : {
735 45 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
736 45 : int bStopProcessing = FALSE;
737 90 : const CPLErr eErr = l_poDS->ValidateRasterIOOrAdviseReadParameters(
738 : "AdviseRead()", &bStopProcessing, nXOff, nYOff, nXSize, nYSize,
739 45 : nBufXSize, nBufYSize, 1, &nBand);
740 45 : if (eErr != CE_None || bStopProcessing)
741 0 : return eErr;
742 :
743 45 : const auto &poArray(l_poDS->m_poArray);
744 90 : std::vector<GUInt64> anArrayStartIdx = m_anOffset;
745 45 : std::vector<size_t> anCount = m_anCount;
746 45 : anArrayStartIdx[l_poDS->m_iXDim] = nXOff;
747 45 : anCount[l_poDS->m_iXDim] = nXSize;
748 45 : if (poArray->GetDimensionCount() >= 2)
749 : {
750 42 : anArrayStartIdx[l_poDS->m_iYDim] = nYOff;
751 42 : anCount[l_poDS->m_iYDim] = nYSize;
752 : }
753 45 : return poArray->AdviseRead(anArrayStartIdx.data(), anCount.data(),
754 : papszOptions)
755 45 : ? CE_None
756 45 : : CE_Failure;
757 : }
758 :
759 : /************************************************************************/
760 : /* IRasterIO() */
761 : /************************************************************************/
762 :
763 525 : CPLErr GDALRasterBandFromArray::IRasterIO(GDALRWFlag eRWFlag, int nXOff,
764 : int nYOff, int nXSize, int nYSize,
765 : void *pData, int nBufXSize,
766 : int nBufYSize, GDALDataType eBufType,
767 : GSpacing nPixelSpaceBuf,
768 : GSpacing nLineSpaceBuf,
769 : GDALRasterIOExtraArg *psExtraArg)
770 : {
771 525 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
772 525 : const auto &poArray(l_poDS->m_poArray);
773 525 : const int nBufferDTSize(GDALGetDataTypeSizeBytes(eBufType));
774 : // If reading/writing at full resolution and with proper stride, go
775 : // directly to the array, but, for performance reasons,
776 : // only if exactly on chunk boundaries, otherwise go through the block cache.
777 525 : if (nXSize == nBufXSize && nYSize == nBufYSize && nBufferDTSize > 0 &&
778 525 : (nPixelSpaceBuf % nBufferDTSize) == 0 &&
779 525 : (nLineSpaceBuf % nBufferDTSize) == 0 && (nXOff % nBlockXSize) == 0 &&
780 523 : (nYOff % nBlockYSize) == 0 &&
781 523 : ((nXSize % nBlockXSize) == 0 || nXOff + nXSize == nRasterXSize) &&
782 520 : ((nYSize % nBlockYSize) == 0 || nYOff + nYSize == nRasterYSize))
783 : {
784 520 : m_anOffset[l_poDS->m_iXDim] = static_cast<GUInt64>(nXOff);
785 520 : m_anCount[l_poDS->m_iXDim] = static_cast<size_t>(nXSize);
786 1040 : m_anStride[l_poDS->m_iXDim] =
787 520 : static_cast<GPtrDiff_t>(nPixelSpaceBuf / nBufferDTSize);
788 520 : if (poArray->GetDimensionCount() >= 2)
789 : {
790 507 : m_anOffset[l_poDS->m_iYDim] = static_cast<GUInt64>(nYOff);
791 507 : m_anCount[l_poDS->m_iYDim] = static_cast<size_t>(nYSize);
792 507 : m_anStride[l_poDS->m_iYDim] =
793 507 : static_cast<GPtrDiff_t>(nLineSpaceBuf / nBufferDTSize);
794 : }
795 520 : if (eRWFlag == GF_Read)
796 : {
797 1002 : return poArray->Read(m_anOffset.data(), m_anCount.data(), nullptr,
798 501 : m_anStride.data(),
799 1002 : GDALExtendedDataType::Create(eBufType), pData)
800 501 : ? CE_None
801 501 : : CE_Failure;
802 : }
803 : else
804 : {
805 38 : return poArray->Write(m_anOffset.data(), m_anCount.data(), nullptr,
806 19 : m_anStride.data(),
807 38 : GDALExtendedDataType::Create(eBufType), pData)
808 19 : ? CE_None
809 19 : : CE_Failure;
810 : }
811 : }
812 : // For unaligned reads, give the array a chance to pre-populate its
813 : // internal chunk cache (e.g. Zarr v3 sharded batches I/O via
814 : // PreloadShardedBlocks). The block cache loop below then hits the
815 : // already-decompressed chunks instead of issuing individual reads.
816 : // Backends that don't override AdviseRead() return true (no-op).
817 5 : if (eRWFlag == GF_Read)
818 : {
819 4 : AdviseRead(nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, eBufType,
820 : nullptr);
821 : }
822 5 : return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
823 : pData, nBufXSize, nBufYSize, eBufType,
824 5 : nPixelSpaceBuf, nLineSpaceBuf, psExtraArg);
825 : }
826 :
827 : /************************************************************************/
828 : /* IsContiguousSequence() */
829 : /************************************************************************/
830 :
831 11 : static bool IsContiguousSequence(int nBandCount, const int *panBandMap)
832 : {
833 23 : for (int i = 1; i < nBandCount; ++i)
834 : {
835 12 : if (panBandMap[i] != panBandMap[i - 1] + 1)
836 0 : return false;
837 : }
838 11 : return true;
839 : }
840 :
841 : /************************************************************************/
842 : /* IRasterIO() */
843 : /************************************************************************/
844 :
845 53 : CPLErr GDALDatasetFromArray::IRasterIO(
846 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
847 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
848 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpaceBuf,
849 : GSpacing nLineSpaceBuf, GSpacing nBandSpaceBuf,
850 : GDALRasterIOExtraArg *psExtraArg)
851 : {
852 53 : const int nBufferDTSize(GDALGetDataTypeSizeBytes(eBufType));
853 : int nBlockXSize, nBlockYSize;
854 53 : papoBands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
855 :
856 : // If reading/writing at full resolution and with proper stride, go
857 : // directly to the array, but, for performance reasons,
858 : // only if exactly on chunk boundaries, or with a pixel interleaved dataset,
859 : // otherwise go through the block cache.
860 64 : if (m_poArray->GetDimensionCount() == 3 && nXSize == nBufXSize &&
861 11 : nYSize == nBufYSize && nBufferDTSize > 0 && nBandCount == nBands &&
862 11 : (nPixelSpaceBuf % nBufferDTSize) == 0 &&
863 11 : (nLineSpaceBuf % nBufferDTSize) == 0 &&
864 11 : (m_bPixelInterleaved ||
865 11 : ((nXOff % nBlockXSize) == 0 && (nYOff % nBlockYSize) == 0 &&
866 11 : ((nXSize % nBlockXSize) == 0 || nXOff + nXSize == nRasterXSize) &&
867 75 : ((nYSize % nBlockYSize) == 0 || nYOff + nYSize == nRasterYSize))) &&
868 11 : IsContiguousSequence(nBandCount, panBandMap))
869 : {
870 11 : m_anOffset[m_iXDim] = static_cast<GUInt64>(nXOff);
871 11 : m_anCount[m_iXDim] = static_cast<size_t>(nXSize);
872 22 : m_anStride[m_iXDim] =
873 11 : static_cast<GPtrDiff_t>(nPixelSpaceBuf / nBufferDTSize);
874 :
875 11 : m_anOffset[m_iYDim] = static_cast<GUInt64>(nYOff);
876 11 : m_anCount[m_iYDim] = static_cast<size_t>(nYSize);
877 22 : m_anStride[m_iYDim] =
878 11 : static_cast<GPtrDiff_t>(nLineSpaceBuf / nBufferDTSize);
879 :
880 11 : size_t iBandDim = 0;
881 11 : for (size_t i = 0; i < 3; ++i)
882 : {
883 11 : if (i != m_iXDim && i != m_iYDim)
884 : {
885 11 : iBandDim = i;
886 11 : break;
887 : }
888 : }
889 :
890 11 : m_anOffset[iBandDim] = panBandMap[0] - 1;
891 11 : m_anCount[iBandDim] = nBandCount;
892 22 : m_anStride[iBandDim] =
893 11 : static_cast<GPtrDiff_t>(nBandSpaceBuf / nBufferDTSize);
894 :
895 11 : if (eRWFlag == GF_Read)
896 : {
897 22 : return m_poArray->Read(m_anOffset.data(), m_anCount.data(), nullptr,
898 11 : m_anStride.data(),
899 22 : GDALExtendedDataType::Create(eBufType),
900 : pData)
901 11 : ? CE_None
902 11 : : CE_Failure;
903 : }
904 : else
905 : {
906 0 : return m_poArray->Write(m_anOffset.data(), m_anCount.data(),
907 0 : nullptr, m_anStride.data(),
908 0 : GDALExtendedDataType::Create(eBufType),
909 : pData)
910 0 : ? CE_None
911 0 : : CE_Failure;
912 : }
913 : }
914 :
915 42 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
916 : nBufXSize, nBufYSize, eBufType, nBandCount,
917 : panBandMap, nPixelSpaceBuf, nLineSpaceBuf,
918 42 : nBandSpaceBuf, psExtraArg);
919 : }
920 :
921 : /************************************************************************/
922 : /* GetColorInterpretation() */
923 : /************************************************************************/
924 :
925 83 : GDALColorInterp GDALRasterBandFromArray::GetColorInterpretation()
926 : {
927 83 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
928 83 : const auto &poArray(l_poDS->m_poArray);
929 249 : auto poAttr = poArray->GetAttribute("COLOR_INTERPRETATION");
930 83 : if (poAttr && poAttr->GetDataType().GetClass() == GEDTC_STRING)
931 : {
932 6 : bool bOK = false;
933 6 : GUInt64 nStartIndex = 0;
934 6 : if (poArray->GetDimensionCount() == 2 &&
935 0 : poAttr->GetDimensionCount() == 0)
936 : {
937 0 : bOK = true;
938 : }
939 6 : else if (poArray->GetDimensionCount() == 3)
940 : {
941 6 : uint64_t nExtraDimSamples = 1;
942 6 : const auto &apoDims = poArray->GetDimensions();
943 24 : for (size_t i = 0; i < apoDims.size(); ++i)
944 : {
945 18 : if (i != l_poDS->m_iXDim && i != l_poDS->m_iYDim)
946 6 : nExtraDimSamples *= apoDims[i]->GetSize();
947 : }
948 6 : if (poAttr->GetDimensionsSize() ==
949 12 : std::vector<GUInt64>{static_cast<GUInt64>(nExtraDimSamples)})
950 : {
951 6 : bOK = true;
952 : }
953 6 : nStartIndex = nBand - 1;
954 : }
955 6 : if (bOK)
956 : {
957 6 : const auto oStringDT = GDALExtendedDataType::CreateString();
958 6 : const size_t nCount = 1;
959 6 : const GInt64 arrayStep = 1;
960 6 : const GPtrDiff_t bufferStride = 1;
961 6 : char *pszValue = nullptr;
962 6 : poAttr->Read(&nStartIndex, &nCount, &arrayStep, &bufferStride,
963 6 : oStringDT, &pszValue);
964 6 : if (pszValue)
965 : {
966 : const auto eColorInterp =
967 6 : GDALGetColorInterpretationByName(pszValue);
968 6 : CPLFree(pszValue);
969 6 : return eColorInterp;
970 : }
971 : }
972 : }
973 77 : return GCI_Undefined;
974 : }
975 :
976 : /************************************************************************/
977 : /* GDALRasterBandFromArray::GetOverviewCount() */
978 : /************************************************************************/
979 :
980 41 : int GDALRasterBandFromArray::GetOverviewCount()
981 : {
982 41 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
983 41 : if (l_poDS->HasRegisteredPAMOverviewFile())
984 : {
985 6 : const int nPAMCount = GDALPamRasterBand::GetOverviewCount();
986 6 : if (nPAMCount)
987 4 : return nPAMCount;
988 : }
989 37 : l_poDS->DiscoverOverviews();
990 37 : return static_cast<int>(l_poDS->m_apoOverviews.size());
991 : }
992 :
993 : /************************************************************************/
994 : /* GDALRasterBandFromArray::GetOverview() */
995 : /************************************************************************/
996 :
997 34 : GDALRasterBand *GDALRasterBandFromArray::GetOverview(int idx)
998 : {
999 34 : auto l_poDS(cpl::down_cast<GDALDatasetFromArray *>(poDS));
1000 34 : if (l_poDS->HasRegisteredPAMOverviewFile())
1001 : {
1002 1 : const int nPAMCount = GDALPamRasterBand::GetOverviewCount();
1003 1 : if (nPAMCount)
1004 1 : return GDALPamRasterBand::GetOverview(idx);
1005 : }
1006 33 : l_poDS->DiscoverOverviews();
1007 33 : if (idx < 0 || static_cast<size_t>(idx) >= l_poDS->m_apoOverviews.size())
1008 : {
1009 8 : return nullptr;
1010 : }
1011 25 : return l_poDS->m_apoOverviews[idx]->GetRasterBand(nBand);
1012 : }
1013 :
1014 : /************************************************************************/
1015 : /* GDALDatasetFromArray::Create() */
1016 : /************************************************************************/
1017 :
1018 529 : std::unique_ptr<GDALDatasetFromArray> GDALDatasetFromArray::Create(
1019 : const std::shared_ptr<GDALMDArray> &array, size_t iXDim, size_t iYDim,
1020 : const std::shared_ptr<GDALGroup> &poRootGroup, CSLConstList papszOptions)
1021 :
1022 : {
1023 529 : const auto nDimCount(array->GetDimensionCount());
1024 529 : if (nDimCount == 0)
1025 : {
1026 1 : CPLError(CE_Failure, CPLE_NotSupported,
1027 : "Unsupported number of dimensions");
1028 1 : return nullptr;
1029 : }
1030 1055 : if (array->GetDataType().GetClass() != GEDTC_NUMERIC ||
1031 527 : array->GetDataType().GetNumericDataType() == GDT_Unknown)
1032 : {
1033 1 : CPLError(CE_Failure, CPLE_NotSupported,
1034 : "Only arrays with numeric data types "
1035 : "can be exposed as classic GDALDataset");
1036 1 : return nullptr;
1037 : }
1038 527 : if (iXDim >= nDimCount || iYDim >= nDimCount ||
1039 501 : (nDimCount >= 2 && iXDim == iYDim))
1040 : {
1041 8 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid iXDim and/or iYDim");
1042 8 : return nullptr;
1043 : }
1044 519 : GUInt64 nTotalBands = 1;
1045 519 : const auto &dims(array->GetDimensions());
1046 1636 : for (size_t i = 0; i < nDimCount; ++i)
1047 : {
1048 1118 : if (i != iXDim && !(nDimCount >= 2 && i == iYDim))
1049 : {
1050 101 : if (dims[i]->GetSize() > 65536 / nTotalBands)
1051 : {
1052 1 : CPLError(CE_Failure, CPLE_AppDefined,
1053 : "Too many bands. Operate on a sliced view");
1054 1 : return nullptr;
1055 : }
1056 100 : nTotalBands *= dims[i]->GetSize();
1057 : }
1058 : }
1059 :
1060 1036 : std::map<std::string, size_t> oMapArrayDimNameToExtraDimIdx;
1061 1036 : std::vector<size_t> oMapArrayExtraDimIdxToOriginalIdx;
1062 1635 : for (size_t i = 0, j = 0; i < nDimCount; ++i)
1063 : {
1064 1117 : if (i != iXDim && !(nDimCount >= 2 && i == iYDim))
1065 : {
1066 100 : oMapArrayDimNameToExtraDimIdx[dims[i]->GetName()] = j;
1067 100 : oMapArrayExtraDimIdxToOriginalIdx.push_back(i);
1068 100 : ++j;
1069 : }
1070 : }
1071 :
1072 518 : const size_t nNewDimCount = nDimCount >= 2 ? nDimCount - 2 : 0;
1073 :
1074 : const char *pszBandMetadata =
1075 518 : CSLFetchNameValue(papszOptions, "BAND_METADATA");
1076 : std::vector<std::vector<MetadataItem>> aoBandParameterMetadataItems(
1077 1036 : nNewDimCount);
1078 518 : if (pszBandMetadata)
1079 : {
1080 32 : if (!poRootGroup)
1081 : {
1082 1 : CPLError(CE_Failure, CPLE_AppDefined,
1083 : "Root group should be provided when BAND_METADATA is set");
1084 24 : return nullptr;
1085 : }
1086 31 : CPLJSONDocument oDoc;
1087 31 : if (!oDoc.LoadMemory(pszBandMetadata))
1088 : {
1089 1 : CPLError(CE_Failure, CPLE_AppDefined,
1090 : "Invalid JSON content for BAND_METADATA");
1091 1 : return nullptr;
1092 : }
1093 30 : auto oRoot = oDoc.GetRoot();
1094 30 : if (oRoot.GetType() != CPLJSONObject::Type::Array)
1095 : {
1096 1 : CPLError(CE_Failure, CPLE_AppDefined,
1097 : "Value of BAND_METADATA should be an array");
1098 1 : return nullptr;
1099 : }
1100 :
1101 29 : auto oArray = oRoot.ToArray();
1102 38 : for (int j = 0; j < oArray.Size(); ++j)
1103 : {
1104 30 : const auto oJsonItem = oArray[j];
1105 30 : MetadataItem oItem;
1106 30 : size_t iExtraDimIdx = 0;
1107 :
1108 60 : const auto osBandArrayFullname = oJsonItem.GetString("array");
1109 60 : const auto osBandAttributeName = oJsonItem.GetString("attribute");
1110 0 : std::shared_ptr<GDALMDArray> poArray;
1111 0 : std::shared_ptr<GDALAttribute> poAttribute;
1112 30 : if (osBandArrayFullname.empty() && osBandAttributeName.empty())
1113 : {
1114 1 : CPLError(CE_Failure, CPLE_AppDefined,
1115 : "BAND_METADATA[%d][\"array\"] or "
1116 : "BAND_METADATA[%d][\"attribute\"] is missing",
1117 : j, j);
1118 1 : return nullptr;
1119 : }
1120 48 : else if (!osBandArrayFullname.empty() &&
1121 19 : !osBandAttributeName.empty())
1122 : {
1123 1 : CPLError(
1124 : CE_Failure, CPLE_AppDefined,
1125 : "BAND_METADATA[%d][\"array\"] and "
1126 : "BAND_METADATA[%d][\"attribute\"] are mutually exclusive",
1127 : j, j);
1128 1 : return nullptr;
1129 : }
1130 28 : else if (!osBandArrayFullname.empty())
1131 : {
1132 : poArray =
1133 18 : poRootGroup->OpenMDArrayFromFullname(osBandArrayFullname);
1134 18 : if (!poArray)
1135 : {
1136 1 : CPLError(CE_Failure, CPLE_AppDefined,
1137 : "Array %s cannot be found",
1138 : osBandArrayFullname.c_str());
1139 3 : return nullptr;
1140 : }
1141 17 : if (poArray->GetDimensionCount() != 1)
1142 : {
1143 1 : CPLError(CE_Failure, CPLE_AppDefined,
1144 : "Array %s is not a 1D array",
1145 : osBandArrayFullname.c_str());
1146 1 : return nullptr;
1147 : }
1148 : const auto &osAuxArrayDimName =
1149 16 : poArray->GetDimensions()[0]->GetName();
1150 : auto oIter =
1151 16 : oMapArrayDimNameToExtraDimIdx.find(osAuxArrayDimName);
1152 16 : if (oIter == oMapArrayDimNameToExtraDimIdx.end())
1153 : {
1154 1 : CPLError(
1155 : CE_Failure, CPLE_AppDefined,
1156 : "Dimension %s of array %s is not a non-X/Y dimension "
1157 : "of array %s",
1158 : osAuxArrayDimName.c_str(), osBandArrayFullname.c_str(),
1159 1 : array->GetName().c_str());
1160 1 : return nullptr;
1161 : }
1162 15 : iExtraDimIdx = oIter->second;
1163 15 : CPLAssert(iExtraDimIdx < nNewDimCount);
1164 : }
1165 : else
1166 : {
1167 10 : CPLAssert(!osBandAttributeName.empty());
1168 10 : poAttribute = !osBandAttributeName.empty() &&
1169 10 : osBandAttributeName[0] == '/'
1170 24 : ? poRootGroup->OpenAttributeFromFullname(
1171 : osBandAttributeName)
1172 10 : : array->GetAttribute(osBandAttributeName);
1173 10 : if (!poAttribute)
1174 : {
1175 2 : CPLError(CE_Failure, CPLE_AppDefined,
1176 : "Attribute %s cannot be found",
1177 : osBandAttributeName.c_str());
1178 8 : return nullptr;
1179 : }
1180 8 : const auto aoAttrDims = poAttribute->GetDimensionsSize();
1181 8 : if (aoAttrDims.size() != 1)
1182 : {
1183 4 : CPLError(CE_Failure, CPLE_AppDefined,
1184 : "Attribute %s is not a 1D array",
1185 : osBandAttributeName.c_str());
1186 4 : return nullptr;
1187 : }
1188 4 : bool found = false;
1189 8 : for (const auto &iter : oMapArrayDimNameToExtraDimIdx)
1190 : {
1191 5 : if (dims[oMapArrayExtraDimIdxToOriginalIdx[iter.second]]
1192 5 : ->GetSize() == aoAttrDims[0])
1193 : {
1194 4 : if (found)
1195 : {
1196 2 : CPLError(CE_Failure, CPLE_AppDefined,
1197 : "Several dimensions of %s have the same "
1198 : "size as attribute %s. Cannot infer which "
1199 : "one to bind to!",
1200 1 : array->GetName().c_str(),
1201 : osBandAttributeName.c_str());
1202 1 : return nullptr;
1203 : }
1204 3 : found = true;
1205 3 : iExtraDimIdx = iter.second;
1206 : }
1207 : }
1208 3 : if (!found)
1209 : {
1210 2 : CPLError(
1211 : CE_Failure, CPLE_AppDefined,
1212 : "No dimension of %s has the same size as attribute %s",
1213 1 : array->GetName().c_str(), osBandAttributeName.c_str());
1214 1 : return nullptr;
1215 : }
1216 : }
1217 :
1218 17 : oItem.osName = oJsonItem.GetString("item_name");
1219 17 : if (oItem.osName.empty())
1220 : {
1221 1 : CPLError(CE_Failure, CPLE_AppDefined,
1222 : "BAND_METADATA[%d][\"item_name\"] is missing", j);
1223 1 : return nullptr;
1224 : }
1225 :
1226 32 : const auto osDefinition = oJsonItem.GetString("item_value", "%s");
1227 :
1228 : // Check correctness of definition
1229 16 : bool bFirstNumericFormatter = true;
1230 16 : std::string osModDefinition;
1231 16 : bool bDefinitionUsesPctForG = false;
1232 79 : for (size_t k = 0; k < osDefinition.size(); ++k)
1233 : {
1234 70 : if (osDefinition[k] == '%')
1235 : {
1236 15 : osModDefinition += osDefinition[k];
1237 15 : if (k + 1 == osDefinition.size())
1238 : {
1239 1 : CPLError(CE_Failure, CPLE_AppDefined,
1240 : "Value of "
1241 : "BAND_METADATA[%d][\"item_value\"] = "
1242 : "%s is invalid at offset %d",
1243 : j, osDefinition.c_str(), int(k));
1244 1 : return nullptr;
1245 : }
1246 14 : ++k;
1247 14 : if (osDefinition[k] == '%')
1248 : {
1249 1 : osModDefinition += osDefinition[k];
1250 1 : continue;
1251 : }
1252 13 : if (!bFirstNumericFormatter)
1253 : {
1254 1 : CPLError(CE_Failure, CPLE_AppDefined,
1255 : "Value of "
1256 : "BAND_METADATA[%d][\"item_value\"] = %s is "
1257 : "invalid at offset %d: %%[x][.y]f|g or %%s "
1258 : "formatters should be specified at most once",
1259 : j, osDefinition.c_str(), int(k));
1260 1 : return nullptr;
1261 : }
1262 12 : bFirstNumericFormatter = false;
1263 19 : for (; k < osDefinition.size(); ++k)
1264 : {
1265 19 : osModDefinition += osDefinition[k];
1266 38 : if (!((osDefinition[k] >= '0' &&
1267 16 : osDefinition[k] <= '9') ||
1268 15 : osDefinition[k] == '.'))
1269 12 : break;
1270 : }
1271 24 : if (k == osDefinition.size() ||
1272 12 : (osDefinition[k] != 'f' && osDefinition[k] != 'g' &&
1273 5 : osDefinition[k] != 's'))
1274 : {
1275 1 : CPLError(CE_Failure, CPLE_AppDefined,
1276 : "Value of "
1277 : "BAND_METADATA[%d][\"item_value\"] = "
1278 : "%s is invalid at offset %d: only "
1279 : "%%[x][.y]f|g or %%s formatters are accepted",
1280 : j, osDefinition.c_str(), int(k));
1281 1 : return nullptr;
1282 : }
1283 11 : bDefinitionUsesPctForG =
1284 11 : (osDefinition[k] == 'f' || osDefinition[k] == 'g');
1285 11 : if (bDefinitionUsesPctForG)
1286 : {
1287 12 : if (poArray &&
1288 12 : poArray->GetDataType().GetClass() != GEDTC_NUMERIC)
1289 : {
1290 1 : CPLError(CE_Failure, CPLE_AppDefined,
1291 : "Data type of %s array is not numeric",
1292 1 : poArray->GetName().c_str());
1293 1 : return nullptr;
1294 : }
1295 8 : else if (poAttribute &&
1296 2 : poAttribute->GetDataType().GetClass() !=
1297 6 : GEDTC_NUMERIC)
1298 : {
1299 0 : CPLError(CE_Failure, CPLE_AppDefined,
1300 : "Data type of %s attribute is not numeric",
1301 0 : poAttribute->GetFullName().c_str());
1302 0 : return nullptr;
1303 : }
1304 : }
1305 : }
1306 62 : else if (osDefinition[k] == '$' &&
1307 62 : k + 1 < osDefinition.size() &&
1308 7 : osDefinition[k + 1] == '{')
1309 : {
1310 7 : const auto nPos = osDefinition.find('}', k);
1311 7 : if (nPos == std::string::npos)
1312 : {
1313 1 : CPLError(CE_Failure, CPLE_AppDefined,
1314 : "Value of "
1315 : "BAND_METADATA[%d][\"item_value\"] = "
1316 : "%s is invalid at offset %d",
1317 : j, osDefinition.c_str(), int(k));
1318 3 : return nullptr;
1319 : }
1320 : const auto osAttrName =
1321 6 : osDefinition.substr(k + 2, nPos - (k + 2));
1322 0 : std::shared_ptr<GDALAttribute> poAttr;
1323 6 : if (poArray && !osAttrName.empty() && osAttrName[0] != '/')
1324 : {
1325 4 : poAttr = poArray->GetAttribute(osAttrName);
1326 4 : if (!poAttr)
1327 : {
1328 1 : CPLError(
1329 : CE_Failure, CPLE_AppDefined,
1330 : "Value of "
1331 : "BAND_METADATA[%d][\"item_value\"] = "
1332 : "%s is invalid: %s is not an attribute of %s",
1333 : j, osDefinition.c_str(), osAttrName.c_str(),
1334 1 : poArray->GetName().c_str());
1335 1 : return nullptr;
1336 : }
1337 : }
1338 : else
1339 : {
1340 : poAttr =
1341 2 : poRootGroup->OpenAttributeFromFullname(osAttrName);
1342 2 : if (!poAttr)
1343 : {
1344 1 : CPLError(CE_Failure, CPLE_AppDefined,
1345 : "Value of "
1346 : "BAND_METADATA[%d][\"item_value\"] = "
1347 : "%s is invalid: %s is not an attribute",
1348 : j, osDefinition.c_str(),
1349 : osAttrName.c_str());
1350 1 : return nullptr;
1351 : }
1352 : }
1353 4 : k = nPos;
1354 4 : const char *pszValue = poAttr->ReadAsString();
1355 4 : if (!pszValue)
1356 : {
1357 0 : CPLError(CE_Failure, CPLE_AppDefined,
1358 : "Cannot get value of attribute %s as a "
1359 : "string",
1360 : osAttrName.c_str());
1361 0 : return nullptr;
1362 : }
1363 4 : osModDefinition += pszValue;
1364 : }
1365 : else
1366 : {
1367 48 : osModDefinition += osDefinition[k];
1368 : }
1369 : }
1370 :
1371 9 : if (poArray)
1372 8 : oItem.poArray = std::move(poArray);
1373 : else
1374 1 : oItem.poArray = std::move(poAttribute);
1375 9 : oItem.osDefinition = std::move(osModDefinition);
1376 9 : oItem.bDefinitionUsesPctForG = bDefinitionUsesPctForG;
1377 :
1378 9 : aoBandParameterMetadataItems[iExtraDimIdx].emplace_back(
1379 9 : std::move(oItem));
1380 : }
1381 : }
1382 :
1383 988 : std::vector<BandImageryMetadata> aoBandImageryMetadata(nNewDimCount);
1384 : const char *pszBandImageryMetadata =
1385 494 : CSLFetchNameValue(papszOptions, "BAND_IMAGERY_METADATA");
1386 494 : if (pszBandImageryMetadata)
1387 : {
1388 20 : if (!poRootGroup)
1389 : {
1390 1 : CPLError(CE_Failure, CPLE_AppDefined,
1391 : "Root group should be provided when BAND_IMAGERY_METADATA "
1392 : "is set");
1393 17 : return nullptr;
1394 : }
1395 19 : CPLJSONDocument oDoc;
1396 19 : if (!oDoc.LoadMemory(pszBandImageryMetadata))
1397 : {
1398 1 : CPLError(CE_Failure, CPLE_AppDefined,
1399 : "Invalid JSON content for BAND_IMAGERY_METADATA");
1400 1 : return nullptr;
1401 : }
1402 18 : auto oRoot = oDoc.GetRoot();
1403 18 : if (oRoot.GetType() != CPLJSONObject::Type::Object)
1404 : {
1405 1 : CPLError(CE_Failure, CPLE_AppDefined,
1406 : "Value of BAND_IMAGERY_METADATA should be an object");
1407 1 : return nullptr;
1408 : }
1409 21 : for (const auto &oJsonItem : oRoot.GetChildren())
1410 : {
1411 38 : if (oJsonItem.GetName() == "CENTRAL_WAVELENGTH_UM" ||
1412 20 : oJsonItem.GetName() == "FWHM_UM")
1413 : {
1414 34 : const auto osBandArrayFullname = oJsonItem.GetString("array");
1415 : const auto osBandAttributeName =
1416 34 : oJsonItem.GetString("attribute");
1417 0 : std::shared_ptr<GDALMDArray> poArray;
1418 0 : std::shared_ptr<GDALAttribute> poAttribute;
1419 17 : size_t iExtraDimIdx = 0;
1420 17 : if (osBandArrayFullname.empty() && osBandAttributeName.empty())
1421 : {
1422 2 : CPLError(CE_Failure, CPLE_AppDefined,
1423 : "BAND_IMAGERY_METADATA[\"%s\"][\"array\"] or "
1424 : "BAND_IMAGERY_METADATA[\"%s\"][\"attribute\"] is "
1425 : "missing",
1426 2 : oJsonItem.GetName().c_str(),
1427 2 : oJsonItem.GetName().c_str());
1428 1 : return nullptr;
1429 : }
1430 25 : else if (!osBandArrayFullname.empty() &&
1431 9 : !osBandAttributeName.empty())
1432 : {
1433 2 : CPLError(CE_Failure, CPLE_AppDefined,
1434 : "BAND_IMAGERY_METADATA[\"%s\"][\"array\"] and "
1435 : "BAND_IMAGERY_METADATA[\"%s\"][\"attribute\"] are "
1436 : "mutually exclusive",
1437 2 : oJsonItem.GetName().c_str(),
1438 2 : oJsonItem.GetName().c_str());
1439 1 : return nullptr;
1440 : }
1441 15 : else if (!osBandArrayFullname.empty())
1442 : {
1443 16 : poArray = poRootGroup->OpenMDArrayFromFullname(
1444 8 : osBandArrayFullname);
1445 8 : if (!poArray)
1446 : {
1447 1 : CPLError(CE_Failure, CPLE_AppDefined,
1448 : "Array %s cannot be found",
1449 : osBandArrayFullname.c_str());
1450 3 : return nullptr;
1451 : }
1452 7 : if (poArray->GetDimensionCount() != 1)
1453 : {
1454 1 : CPLError(CE_Failure, CPLE_AppDefined,
1455 : "Array %s is not a 1D array",
1456 : osBandArrayFullname.c_str());
1457 1 : return nullptr;
1458 : }
1459 : const auto &osAuxArrayDimName =
1460 6 : poArray->GetDimensions()[0]->GetName();
1461 : auto oIter =
1462 6 : oMapArrayDimNameToExtraDimIdx.find(osAuxArrayDimName);
1463 6 : if (oIter == oMapArrayDimNameToExtraDimIdx.end())
1464 : {
1465 1 : CPLError(CE_Failure, CPLE_AppDefined,
1466 : "Dimension \"%s\" of array \"%s\" is not a "
1467 : "non-X/Y dimension of array \"%s\"",
1468 : osAuxArrayDimName.c_str(),
1469 : osBandArrayFullname.c_str(),
1470 1 : array->GetName().c_str());
1471 1 : return nullptr;
1472 : }
1473 5 : iExtraDimIdx = oIter->second;
1474 5 : CPLAssert(iExtraDimIdx < nNewDimCount);
1475 : }
1476 : else
1477 : {
1478 : poAttribute =
1479 7 : !osBandAttributeName.empty() &&
1480 7 : osBandAttributeName[0] == '/'
1481 16 : ? poRootGroup->OpenAttributeFromFullname(
1482 : osBandAttributeName)
1483 7 : : array->GetAttribute(osBandAttributeName);
1484 7 : if (!poAttribute)
1485 : {
1486 2 : CPLError(CE_Failure, CPLE_AppDefined,
1487 : "Attribute %s cannot be found",
1488 : osBandAttributeName.c_str());
1489 6 : return nullptr;
1490 : }
1491 5 : const auto aoAttrDims = poAttribute->GetDimensionsSize();
1492 5 : if (aoAttrDims.size() != 1)
1493 : {
1494 2 : CPLError(CE_Failure, CPLE_AppDefined,
1495 : "Attribute %s is not a 1D array",
1496 : osBandAttributeName.c_str());
1497 2 : return nullptr;
1498 : }
1499 3 : bool found = false;
1500 6 : for (const auto &iter : oMapArrayDimNameToExtraDimIdx)
1501 : {
1502 4 : if (dims[oMapArrayExtraDimIdxToOriginalIdx[iter.second]]
1503 4 : ->GetSize() == aoAttrDims[0])
1504 : {
1505 3 : if (found)
1506 : {
1507 2 : CPLError(CE_Failure, CPLE_AppDefined,
1508 : "Several dimensions of %s have the "
1509 : "same size as attribute %s. Cannot "
1510 : "infer which one to bind to!",
1511 1 : array->GetName().c_str(),
1512 : osBandAttributeName.c_str());
1513 1 : return nullptr;
1514 : }
1515 2 : found = true;
1516 2 : iExtraDimIdx = iter.second;
1517 : }
1518 : }
1519 2 : if (!found)
1520 : {
1521 2 : CPLError(CE_Failure, CPLE_AppDefined,
1522 : "No dimension of %s has the same size as "
1523 : "attribute %s",
1524 1 : array->GetName().c_str(),
1525 : osBandAttributeName.c_str());
1526 1 : return nullptr;
1527 : }
1528 : }
1529 :
1530 12 : std::string osUnit = oJsonItem.GetString("unit", "um");
1531 6 : if (STARTS_WITH(osUnit.c_str(), "${"))
1532 : {
1533 4 : if (osUnit.back() != '}')
1534 : {
1535 2 : CPLError(CE_Failure, CPLE_AppDefined,
1536 : "Value of "
1537 : "BAND_IMAGERY_METADATA[\"%s\"][\"unit\"] = "
1538 : "%s is invalid",
1539 2 : oJsonItem.GetName().c_str(), osUnit.c_str());
1540 2 : return nullptr;
1541 : }
1542 3 : const auto osAttrName = osUnit.substr(2, osUnit.size() - 3);
1543 0 : std::shared_ptr<GDALAttribute> poAttr;
1544 3 : if (poArray && !osAttrName.empty() && osAttrName[0] != '/')
1545 : {
1546 2 : poAttr = poArray->GetAttribute(osAttrName);
1547 2 : if (!poAttr)
1548 : {
1549 2 : CPLError(
1550 : CE_Failure, CPLE_AppDefined,
1551 : "Value of "
1552 : "BAND_IMAGERY_METADATA[\"%s\"][\"unit\"] = "
1553 : "%s is invalid: %s is not an attribute of %s",
1554 2 : oJsonItem.GetName().c_str(), osUnit.c_str(),
1555 : osAttrName.c_str(),
1556 : osBandArrayFullname.c_str());
1557 1 : return nullptr;
1558 : }
1559 : }
1560 : else
1561 : {
1562 : poAttr =
1563 1 : poRootGroup->OpenAttributeFromFullname(osAttrName);
1564 1 : if (!poAttr)
1565 : {
1566 0 : CPLError(
1567 : CE_Failure, CPLE_AppDefined,
1568 : "Value of "
1569 : "BAND_IMAGERY_METADATA[\"%s\"][\"unit\"] = "
1570 : "%s is invalid: %s is not an attribute",
1571 0 : oJsonItem.GetName().c_str(), osUnit.c_str(),
1572 : osAttrName.c_str());
1573 0 : return nullptr;
1574 : }
1575 : }
1576 :
1577 2 : const char *pszValue = poAttr->ReadAsString();
1578 2 : if (!pszValue)
1579 : {
1580 0 : CPLError(CE_Failure, CPLE_AppDefined,
1581 : "Cannot get value of attribute %s of %s as a "
1582 : "string",
1583 : osAttrName.c_str(),
1584 : osBandArrayFullname.c_str());
1585 0 : return nullptr;
1586 : }
1587 2 : osUnit = pszValue;
1588 : }
1589 4 : double dfConvToUM = 1.0;
1590 10 : if (osUnit == "nm" || osUnit == "nanometre" ||
1591 13 : osUnit == "nanometres" || osUnit == "nanometer" ||
1592 3 : osUnit == "nanometers")
1593 : {
1594 1 : dfConvToUM = 1e-3;
1595 : }
1596 5 : else if (!(osUnit == "um" || osUnit == "micrometre" ||
1597 2 : osUnit == "micrometres" || osUnit == "micrometer" ||
1598 1 : osUnit == "micrometers"))
1599 : {
1600 2 : CPLError(CE_Failure, CPLE_AppDefined,
1601 : "Unhandled value for "
1602 : "BAND_IMAGERY_METADATA[\"%s\"][\"unit\"] = %s",
1603 2 : oJsonItem.GetName().c_str(), osUnit.c_str());
1604 1 : return nullptr;
1605 : }
1606 :
1607 3 : BandImageryMetadata &item = aoBandImageryMetadata[iExtraDimIdx];
1608 :
1609 3 : std::shared_ptr<GDALAbstractMDArray> abstractArray;
1610 3 : if (poArray)
1611 2 : abstractArray = std::move(poArray);
1612 : else
1613 1 : abstractArray = std::move(poAttribute);
1614 3 : if (oJsonItem.GetName() == "CENTRAL_WAVELENGTH_UM")
1615 : {
1616 2 : item.poCentralWavelengthArray = std::move(abstractArray);
1617 2 : item.dfCentralWavelengthToMicrometer = dfConvToUM;
1618 : }
1619 : else
1620 : {
1621 1 : item.poFWHMArray = std::move(abstractArray);
1622 1 : item.dfFWHMToMicrometer = dfConvToUM;
1623 : }
1624 : }
1625 : else
1626 : {
1627 1 : CPLError(CE_Warning, CPLE_AppDefined,
1628 : "Ignored member \"%s\" in BAND_IMAGERY_METADATA",
1629 2 : oJsonItem.GetName().c_str());
1630 : }
1631 : }
1632 : }
1633 :
1634 458 : if ((nDimCount >= 2 &&
1635 954 : dims[iYDim]->GetSize() > static_cast<uint64_t>(INT_MAX)) ||
1636 477 : dims[iXDim]->GetSize() > static_cast<uint64_t>(INT_MAX))
1637 : {
1638 1 : CPLError(CE_Failure, CPLE_AppDefined,
1639 : "Array is too large to be exposed as a GDAL dataset");
1640 1 : return nullptr;
1641 : }
1642 :
1643 : auto poDS = std::make_unique<GDALDatasetFromArray>(
1644 952 : array, iXDim, iYDim, CPLStringList(papszOptions));
1645 :
1646 476 : poDS->eAccess = array->IsWritable() ? GA_Update : GA_ReadOnly;
1647 :
1648 476 : poDS->nRasterYSize =
1649 476 : nDimCount < 2 ? 1 : static_cast<int>(dims[iYDim]->GetSize());
1650 :
1651 476 : poDS->nRasterXSize = static_cast<int>(dims[iXDim]->GetSize());
1652 :
1653 952 : std::vector<GUInt64> anOtherDimCoord(nNewDimCount);
1654 952 : std::vector<GUInt64> anStackIters(nDimCount);
1655 952 : std::vector<size_t> anMapNewToOld(nNewDimCount);
1656 1466 : for (size_t i = 0, j = 0; i < nDimCount; ++i)
1657 : {
1658 990 : if (i != iXDim && !(nDimCount >= 2 && i == iYDim))
1659 : {
1660 57 : anMapNewToOld[j] = i;
1661 57 : j++;
1662 : }
1663 : }
1664 :
1665 476 : poDS->m_bHasGT = array->GuessGeoTransform(iXDim, iYDim, false, poDS->m_gt);
1666 :
1667 952 : const auto attrs(array->GetAttributes());
1668 650 : for (const auto &attr : attrs)
1669 : {
1670 174 : if (attr->GetName() == "spatial:registration")
1671 : {
1672 : // From https://github.com/zarr-conventions/spatial
1673 6 : const char *pszValue = attr->ReadAsString();
1674 6 : if (pszValue && strcmp(pszValue, "pixel") == 0)
1675 5 : poDS->m_oMDD.SetMetadataItem(GDALMD_AREA_OR_POINT,
1676 : GDALMD_AOP_AREA);
1677 1 : else if (pszValue && strcmp(pszValue, "node") == 0)
1678 1 : poDS->m_oMDD.SetMetadataItem(GDALMD_AREA_OR_POINT,
1679 : GDALMD_AOP_POINT);
1680 0 : else if (pszValue)
1681 0 : poDS->m_oMDD.SetMetadataItem(attr->GetName().c_str(), pszValue);
1682 : }
1683 168 : else if (attr->GetName() == "gdal:geotransform")
1684 : {
1685 : // From Zarr driver
1686 4 : const auto doubleArray = attr->ReadAsDoubleArray();
1687 2 : if (doubleArray.size() == 6)
1688 : {
1689 2 : poDS->m_bHasGT = true;
1690 2 : poDS->m_gt = GDALGeoTransform(doubleArray.data());
1691 : }
1692 : }
1693 166 : else if (attr->GetName() != "COLOR_INTERPRETATION")
1694 : {
1695 308 : auto stringArray = attr->ReadAsStringArray();
1696 308 : std::string val;
1697 154 : if (stringArray.size() > 1)
1698 : {
1699 61 : val += '{';
1700 : }
1701 595 : for (int i = 0; i < stringArray.size(); ++i)
1702 : {
1703 441 : if (i > 0)
1704 287 : val += ',';
1705 441 : val += stringArray[i];
1706 : }
1707 154 : if (stringArray.size() > 1)
1708 : {
1709 61 : val += '}';
1710 : }
1711 154 : poDS->m_oMDD.SetMetadataItem(attr->GetName().c_str(), val.c_str());
1712 : }
1713 : }
1714 :
1715 476 : const char *pszDelay = CSLFetchNameValueDef(
1716 : papszOptions, "LOAD_EXTRA_DIM_METADATA_DELAY",
1717 : CPLGetConfigOption("GDAL_LOAD_EXTRA_DIM_METADATA_DELAY", "5"));
1718 : const double dfDelay =
1719 476 : EQUAL(pszDelay, "unlimited") ? -1 : CPLAtof(pszDelay);
1720 476 : const auto nStartTime = time(nullptr);
1721 476 : bool bHasWarned = false;
1722 : // Instantiate bands by iterating over non-XY variables
1723 476 : size_t iDim = 0;
1724 476 : int nCurBand = 1;
1725 624 : lbl_next_depth:
1726 624 : if (iDim < nNewDimCount)
1727 : {
1728 59 : anStackIters[iDim] = dims[anMapNewToOld[iDim]]->GetSize();
1729 59 : anOtherDimCoord[iDim] = 0;
1730 : while (true)
1731 : {
1732 148 : ++iDim;
1733 148 : goto lbl_next_depth;
1734 148 : lbl_return_to_caller:
1735 148 : --iDim;
1736 148 : --anStackIters[iDim];
1737 148 : if (anStackIters[iDim] == 0)
1738 59 : break;
1739 89 : ++anOtherDimCoord[iDim];
1740 : }
1741 : }
1742 : else
1743 : {
1744 1130 : poDS->SetBand(nCurBand,
1745 565 : std::make_unique<GDALRasterBandFromArray>(
1746 565 : poDS.get(), anOtherDimCoord,
1747 : aoBandParameterMetadataItems, aoBandImageryMetadata,
1748 : dfDelay, nStartTime, bHasWarned));
1749 565 : ++nCurBand;
1750 : }
1751 624 : if (iDim > 0)
1752 148 : goto lbl_return_to_caller;
1753 :
1754 484 : if (nDimCount == 3 && iXDim <= 1 && iYDim <= 1 &&
1755 8 : poDS->GetRasterCount() > 1)
1756 : {
1757 8 : poDS->m_bPixelInterleaved = true;
1758 8 : poDS->m_oMDD.SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1759 : }
1760 :
1761 910 : if (!array->GetFilename().empty() &&
1762 434 : CPLTestBool(CSLFetchNameValueDef(papszOptions, "LOAD_PAM", "YES")))
1763 : {
1764 424 : poDS->SetPhysicalFilename(array->GetFilename().c_str());
1765 : std::string osDerivedDatasetName(
1766 : CPLSPrintf("AsClassicDataset(%d,%d) view of %s", int(iXDim),
1767 848 : int(iYDim), array->GetFullName().c_str()));
1768 424 : if (!array->GetContext().empty())
1769 : {
1770 2 : osDerivedDatasetName += " with context ";
1771 2 : osDerivedDatasetName += array->GetContext();
1772 : }
1773 424 : poDS->SetDerivedDatasetName(osDerivedDatasetName.c_str());
1774 424 : poDS->TryLoadXML();
1775 :
1776 2 : for (const auto &[pszKey, pszValue] :
1777 : cpl::IterateNameValue(static_cast<CSLConstList>(
1778 426 : poDS->GDALPamDataset::GetMetadata())))
1779 : {
1780 1 : poDS->m_oMDD.SetMetadataItem(pszKey, pszValue);
1781 : }
1782 : }
1783 :
1784 476 : return poDS;
1785 : }
1786 :
1787 : //! @endcond
1788 :
1789 : /************************************************************************/
1790 : /* AsClassicDataset() */
1791 : /************************************************************************/
1792 :
1793 : /** Return a view of this array as a "classic" GDALDataset (ie 2D)
1794 : *
1795 : * In the case of > 2D arrays, additional dimensions will be represented as
1796 : * raster bands.
1797 : *
1798 : * The "reverse" method is GDALRasterBand::AsMDArray().
1799 : *
1800 : * This is the same as the C function GDALMDArrayAsClassicDataset().
1801 : *
1802 : * @param iXDim Index of the dimension that will be used as the X/width axis.
1803 : * @param iYDim Index of the dimension that will be used as the Y/height axis.
1804 : * Ignored if the dimension count is 1.
1805 : * @param poRootGroup (Added in GDAL 3.8) Root group. Used with the BAND_METADATA
1806 : * and BAND_IMAGERY_METADATA option.
1807 : * @param papszOptions (Added in GDAL 3.8) Null-terminated list of options, or
1808 : * nullptr. Current supported options are:
1809 : * <ul>
1810 : * <li>BAND_METADATA: JSON serialized array defining which
1811 : * arrays of the poRootGroup, indexed by non-X and Y
1812 : * dimensions, should be mapped as band metadata items.
1813 : * Each array item should be an object with the
1814 : * following members:
1815 : * - "array": full name of a band parameter array.
1816 : * Such array must be a one
1817 : * dimensional array, and its dimension must be one of
1818 : * the dimensions of the array on which the method is
1819 : * called (excluding the X and Y dimensions).
1820 : * - "attribute": name relative to *this array or full
1821 : * name of a single dimension numeric array whose size
1822 : * must be one of the dimensions of *this array
1823 : * (excluding the X and Y dimensions).
1824 : * "array" and "attribute" are mutually exclusive.
1825 : * - "item_name": band metadata item name
1826 : * - "item_value": (optional) String, where "%[x][.y]f",
1827 : * "%[x][.y]g" or "%s" printf-like formatting can be
1828 : * used to format the corresponding value of the
1829 : * parameter array. The percentage character should be
1830 : * repeated: "%%"
1831 : * "${attribute_name}" can also be used to include the
1832 : * value of an attribute for "array" when set and if
1833 : * not starting with '/'. Otherwise if starting with
1834 : * '/', it is the full path to the attribute.
1835 : *
1836 : * If "item_value" is not provided, a default formatting
1837 : * of the value will be applied.
1838 : *
1839 : * Example:
1840 : * [
1841 : * {
1842 : * "array": "/sensor_band_parameters/wavelengths",
1843 : * "item_name": "WAVELENGTH",
1844 : * "item_value": "%.1f ${units}"
1845 : * },
1846 : * {
1847 : * "array": "/sensor_band_parameters/fwhm",
1848 : * "item_name": "FWHM"
1849 : * },
1850 : * {
1851 : * "array": "/sensor_band_parameters/fwhm",
1852 : * "item_name": "FWHM_UNIT",
1853 : * "item_value": "${units}"
1854 : * }
1855 : * ]
1856 : *
1857 : * Example for Planet Labs Tanager radiance products:
1858 : * [
1859 : * {
1860 : * "attribute": "center_wavelengths",
1861 : * "item_name": "WAVELENGTH",
1862 : * "item_value": "%.1f ${center_wavelengths_units}"
1863 : * }
1864 : * ]
1865 : *
1866 : * </li>
1867 : * <li>BAND_IMAGERY_METADATA: (GDAL >= 3.11)
1868 : * JSON serialized object defining which arrays of the
1869 : * poRootGroup, indexed by non-X and Y dimensions,
1870 : * should be mapped as band metadata items in the
1871 : * band IMAGERY domain.
1872 : * The object currently accepts 2 members:
1873 : * - "CENTRAL_WAVELENGTH_UM": Central Wavelength in
1874 : * micrometers.
1875 : * - "FWHM_UM": Full-width half-maximum
1876 : * in micrometers.
1877 : * The value of each member should be an object with the
1878 : * following members:
1879 : * - "array": full name of a band parameter array.
1880 : * Such array must be a one dimensional array, and its
1881 : * dimension must be one of the dimensions of the
1882 : * array on which the method is called
1883 : * (excluding the X and Y dimensions).
1884 : * - "attribute": name relative to *this array or full
1885 : * name of a single dimension numeric array whose size
1886 : * must be one of the dimensions of *this array
1887 : * (excluding the X and Y dimensions).
1888 : * "array" and "attribute" are mutually exclusive,
1889 : * and one of them is required.
1890 : * - "unit": (optional) unit of the values pointed in
1891 : * the above array.
1892 : * Can be a literal string or a string of the form
1893 : * "${attribute_name}" to point to an attribute for
1894 : * "array" when set and if no starting
1895 : * with '/'. Otherwise if starting with '/', it is
1896 : * the full path to the attribute.
1897 : * Accepted values are "um", "micrometer"
1898 : * (with UK vs US spelling, singular or plural), "nm",
1899 : * "nanometer" (with UK vs US spelling, singular or
1900 : * plural)
1901 : * If not provided, micrometer is assumed.
1902 : *
1903 : * Example for EMIT datasets:
1904 : * {
1905 : * "CENTRAL_WAVELENGTH_UM": {
1906 : * "array": "/sensor_band_parameters/wavelengths",
1907 : * "unit": "${units}"
1908 : * },
1909 : * "FWHM_UM": {
1910 : * "array": "/sensor_band_parameters/fwhm",
1911 : * "unit": "${units}"
1912 : * }
1913 : * }
1914 : *
1915 : * Example for Planet Labs Tanager radiance products:
1916 : * {
1917 : * "CENTRAL_WAVELENGTH_UM": {
1918 : * "attribute": "center_wavelengths",
1919 : * "unit": "${center_wavelengths_units}"
1920 : * },
1921 : * "FWHM_UM": {
1922 : * "attribute": "fwhm",
1923 : * "unit": "${fwhm_units}"
1924 : * }
1925 : * }
1926 : *
1927 : * </li>
1928 : * <li>LOAD_EXTRA_DIM_METADATA_DELAY: Maximum delay in
1929 : * seconds allowed to set the DIM_{dimname}_VALUE band
1930 : * metadata items from the indexing variable of the
1931 : * dimensions.
1932 : * Default value is 5. 'unlimited' can be used to mean
1933 : * unlimited delay. Can also be defined globally with
1934 : * the GDAL_LOAD_EXTRA_DIM_METADATA_DELAY configuration
1935 : * option.</li>
1936 : * </ul>
1937 : * @return a new GDALDataset that must be freed with GDALClose(), or nullptr
1938 : */
1939 : GDALDataset *
1940 518 : GDALMDArray::AsClassicDataset(size_t iXDim, size_t iYDim,
1941 : const std::shared_ptr<GDALGroup> &poRootGroup,
1942 : CSLConstList papszOptions) const
1943 : {
1944 1036 : auto self = std::dynamic_pointer_cast<GDALMDArray>(m_pSelf.lock());
1945 518 : if (!self)
1946 : {
1947 0 : CPLError(CE_Failure, CPLE_AppDefined,
1948 : "Driver implementation issue: m_pSelf not set !");
1949 0 : return nullptr;
1950 : }
1951 1036 : return GDALDatasetFromArray::Create(self, iXDim, iYDim, poRootGroup,
1952 : papszOptions)
1953 518 : .release();
1954 : }
|
