Line data Source code
1 : /******************************************************************************
2 : * Name: gdalmultidim_subsetdimension.cpp
3 : * Project: GDAL Core
4 : * Purpose: GDALGroup::SubsetDimensionFromSelection() implementation
5 : * Author: Even Rouault <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdal_priv.h"
14 : #include "gdal_pam.h"
15 :
16 : #include <algorithm>
17 :
18 : /************************************************************************/
19 : /* GetParentName() */
20 : /************************************************************************/
21 :
22 42 : static std::string GetParentName(const std::string &osPath)
23 : {
24 42 : if (osPath == "/" || osPath.rfind('/') == 0)
25 36 : return "/";
26 6 : return osPath.substr(0, osPath.rfind('/'));
27 : }
28 :
29 : /************************************************************************/
30 : /* GDALSubsetGroupSharedResources */
31 : /************************************************************************/
32 :
33 : struct GDALSubsetGroupSharedResources
34 : {
35 : std::shared_ptr<GDALGroup> m_poRootGroup{}; // may be nullptr
36 : std::string m_osDimFullName{};
37 : std::vector<int> m_anMapNewDimToOldDim{};
38 : std::string m_osSelection{};
39 : std::shared_ptr<GDALDimension> m_poNewDim{};
40 : std::shared_ptr<GDALMDArray> m_poNewIndexingVar{};
41 : };
42 :
43 : /************************************************************************/
44 : /* CreateContext() */
45 : /************************************************************************/
46 :
47 : static std::string
48 14 : CreateContext(const std::string &osParentContext,
49 : const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared)
50 : {
51 14 : std::string osRet(osParentContext);
52 14 : if (!osRet.empty())
53 0 : osRet += ". ";
54 14 : osRet += "Selection ";
55 14 : osRet += poShared->m_osSelection;
56 14 : return osRet;
57 : }
58 :
59 : /************************************************************************/
60 : /* GDALSubsetGroup */
61 : /************************************************************************/
62 :
63 : class GDALSubsetGroup final : public GDALGroup
64 : {
65 : std::shared_ptr<GDALGroup> m_poParent{};
66 : std::shared_ptr<GDALSubsetGroupSharedResources> m_poShared{};
67 :
68 9 : GDALSubsetGroup(
69 : const std::shared_ptr<GDALGroup> &poParent,
70 : const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared)
71 18 : : GDALGroup(GetParentName(poParent->GetFullName()), poParent->GetName(),
72 18 : CreateContext(poParent->GetContext(), poShared)),
73 36 : m_poParent(std::move(poParent)), m_poShared(std::move(poShared))
74 : {
75 9 : }
76 :
77 : public:
78 : static std::shared_ptr<GDALGroup>
79 9 : Create(const std::shared_ptr<GDALGroup> &poParent,
80 : const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared)
81 : {
82 : auto poGroup = std::shared_ptr<GDALSubsetGroup>(
83 18 : new GDALSubsetGroup(poParent, poShared));
84 9 : poGroup->SetSelf(poGroup);
85 18 : return poGroup;
86 : }
87 :
88 : std::vector<std::string>
89 2 : GetMDArrayNames(CSLConstList papszOptions = nullptr) const override
90 : {
91 2 : return m_poParent->GetMDArrayNames(papszOptions);
92 : }
93 :
94 : std::shared_ptr<GDALMDArray>
95 : OpenMDArray(const std::string &osName,
96 : CSLConstList papszOptions = nullptr) const override;
97 :
98 : std::vector<std::string>
99 1 : GetGroupNames(CSLConstList papszOptions = nullptr) const override
100 : {
101 1 : return m_poParent->GetGroupNames(papszOptions);
102 : }
103 :
104 : std::shared_ptr<GDALGroup>
105 : OpenGroup(const std::string &osName,
106 : CSLConstList papszOptions = nullptr) const override;
107 :
108 : std::vector<std::shared_ptr<GDALDimension>>
109 : GetDimensions(CSLConstList papszOptions = nullptr) const override;
110 :
111 : std::shared_ptr<GDALAttribute>
112 2 : GetAttribute(const std::string &osName) const override
113 : {
114 2 : return m_poParent->GetAttribute(osName);
115 : }
116 :
117 : std::vector<std::shared_ptr<GDALAttribute>>
118 1 : GetAttributes(CSLConstList papszOptions = nullptr) const override
119 : {
120 1 : return m_poParent->GetAttributes(papszOptions);
121 : }
122 : };
123 :
124 : /************************************************************************/
125 : /* GDALSubsetArray */
126 : /************************************************************************/
127 :
128 : class GDALSubsetArray final : public GDALPamMDArray
129 : {
130 : private:
131 : std::shared_ptr<GDALMDArray> m_poParent{};
132 : std::shared_ptr<GDALSubsetGroupSharedResources> m_poShared{};
133 : std::vector<std::shared_ptr<GDALDimension>> m_apoDims{};
134 : std::vector<bool> m_abPatchedDim{};
135 : bool m_bPatchedDimIsFirst = false;
136 :
137 : protected:
138 14 : GDALSubsetArray(
139 : const std::shared_ptr<GDALMDArray> &poParent,
140 : const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared,
141 : const std::string &osContext)
142 28 : : GDALAbstractMDArray(GetParentName(poParent->GetFullName()),
143 14 : poParent->GetName()),
144 28 : GDALPamMDArray(GetParentName(poParent->GetFullName()),
145 28 : poParent->GetName(), GDALPamMultiDim::GetPAM(poParent),
146 : osContext),
147 84 : m_poParent(std::move(poParent)), m_poShared(std::move(poShared))
148 : {
149 14 : m_apoDims = m_poParent->GetDimensions();
150 39 : for (size_t i = 0; i < m_apoDims.size(); ++i)
151 : {
152 25 : auto &poDim = m_apoDims[i];
153 25 : if (poDim->GetFullName() == m_poShared->m_osDimFullName)
154 : {
155 15 : m_bPatchedDimIsFirst = (i == 0);
156 15 : poDim = m_poShared->m_poNewDim;
157 15 : m_abPatchedDim.push_back(true);
158 : }
159 : else
160 : {
161 10 : m_abPatchedDim.push_back(false);
162 : }
163 : }
164 14 : }
165 :
166 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
167 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
168 : const GDALExtendedDataType &bufferDataType,
169 : void *pDstBuffer) const override;
170 :
171 : public:
172 : static std::shared_ptr<GDALSubsetArray>
173 14 : Create(const std::shared_ptr<GDALMDArray> &poParent,
174 : const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared,
175 : const std::string &osContext)
176 : {
177 : auto newAr(std::shared_ptr<GDALSubsetArray>(
178 14 : new GDALSubsetArray(poParent, poShared, osContext)));
179 14 : newAr->SetSelf(newAr);
180 14 : return newAr;
181 : }
182 :
183 2 : bool IsWritable() const override
184 : {
185 2 : return false;
186 : }
187 :
188 19 : const std::string &GetFilename() const override
189 : {
190 19 : return m_poParent->GetFilename();
191 : }
192 :
193 : const std::vector<std::shared_ptr<GDALDimension>> &
194 109 : GetDimensions() const override
195 : {
196 109 : return m_apoDims;
197 : }
198 :
199 62 : const GDALExtendedDataType &GetDataType() const override
200 : {
201 62 : return m_poParent->GetDataType();
202 : }
203 :
204 3 : const std::string &GetUnit() const override
205 : {
206 3 : return m_poParent->GetUnit();
207 : }
208 :
209 1 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
210 : {
211 1 : return m_poParent->GetSpatialRef();
212 : }
213 :
214 8 : const void *GetRawNoDataValue() const override
215 : {
216 8 : return m_poParent->GetRawNoDataValue();
217 : }
218 :
219 4 : std::vector<GUInt64> GetBlockSize() const override
220 : {
221 4 : std::vector<GUInt64> ret(m_poParent->GetBlockSize());
222 13 : for (size_t i = 0; i < m_apoDims.size(); ++i)
223 : {
224 9 : if (m_abPatchedDim[i])
225 4 : ret[1] = 1;
226 : }
227 4 : return ret;
228 : }
229 :
230 : std::shared_ptr<GDALAttribute>
231 8 : GetAttribute(const std::string &osName) const override
232 : {
233 8 : return m_poParent->GetAttribute(osName);
234 : }
235 :
236 : std::vector<std::shared_ptr<GDALAttribute>>
237 3 : GetAttributes(CSLConstList papszOptions = nullptr) const override
238 : {
239 3 : return m_poParent->GetAttributes(papszOptions);
240 : }
241 :
242 1 : std::shared_ptr<GDALGroup> GetRootGroup() const override
243 : {
244 1 : if (m_poShared->m_poRootGroup)
245 : {
246 2 : return GDALSubsetGroup::Create(m_poShared->m_poRootGroup,
247 1 : m_poShared);
248 : }
249 0 : return nullptr;
250 : }
251 : };
252 :
253 : /************************************************************************/
254 : /* OpenMDArray() */
255 : /************************************************************************/
256 :
257 : std::shared_ptr<GDALMDArray>
258 12 : GDALSubsetGroup::OpenMDArray(const std::string &osName,
259 : CSLConstList papszOptions) const
260 : {
261 24 : auto poArray = m_poParent->OpenMDArray(osName, papszOptions);
262 12 : if (poArray)
263 : {
264 22 : for (const auto &poDim : poArray->GetDimensions())
265 : {
266 20 : if (poDim->GetFullName() == m_poShared->m_osDimFullName)
267 : {
268 18 : return GDALSubsetArray::Create(poArray, m_poShared,
269 9 : GetContext());
270 : }
271 : }
272 : }
273 3 : return poArray;
274 : }
275 :
276 : /************************************************************************/
277 : /* OpenGroup() */
278 : /************************************************************************/
279 :
280 : std::shared_ptr<GDALGroup>
281 5 : GDALSubsetGroup::OpenGroup(const std::string &osName,
282 : CSLConstList papszOptions) const
283 : {
284 5 : auto poSubGroup = m_poParent->OpenGroup(osName, papszOptions);
285 5 : if (poSubGroup)
286 : {
287 3 : poSubGroup = GDALSubsetGroup::Create(poSubGroup, m_poShared);
288 : }
289 5 : return poSubGroup;
290 : }
291 :
292 : /************************************************************************/
293 : /* GetDimensions() */
294 : /************************************************************************/
295 :
296 : std::vector<std::shared_ptr<GDALDimension>>
297 1 : GDALSubsetGroup::GetDimensions(CSLConstList papszOptions) const
298 : {
299 1 : auto apoDims = m_poParent->GetDimensions(papszOptions);
300 5 : for (auto &poDim : apoDims)
301 : {
302 4 : if (poDim->GetFullName() == m_poShared->m_osDimFullName)
303 : {
304 1 : poDim = m_poShared->m_poNewDim;
305 : }
306 : }
307 1 : return apoDims;
308 : }
309 :
310 : /************************************************************************/
311 : /* IRead() */
312 : /************************************************************************/
313 :
314 28 : bool GDALSubsetArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
315 : const GInt64 *arrayStep,
316 : const GPtrDiff_t *bufferStride,
317 : const GDALExtendedDataType &bufferDataType,
318 : void *pDstBuffer) const
319 : {
320 28 : const auto nDims = m_apoDims.size();
321 56 : std::vector<GUInt64> newArrayStartIdx(nDims);
322 : // the +1 in nDims + 1 is to make happy -Werror=null-dereference when
323 : // doing newCount[0] = 1 and newArrayStep[0] = 1
324 56 : std::vector<size_t> newCount(nDims + 1, 1);
325 56 : std::vector<GInt64> newArrayStep(nDims + 1, 1);
326 28 : const size_t nBufferDTSize = bufferDataType.GetSize();
327 :
328 28 : if (m_bPatchedDimIsFirst)
329 : {
330 : // Optimized case when the only patched dimension is the first one.
331 3 : std::copy_n(arrayStartIdx, nDims, newArrayStartIdx.data());
332 3 : std::copy_n(count, nDims, newCount.data());
333 3 : std::copy_n(arrayStep, nDims, newArrayStep.data());
334 3 : GUInt64 arrayIdx = arrayStartIdx[0];
335 3 : GByte *pabyDstBuffer = static_cast<GByte *>(pDstBuffer);
336 : #if defined(__GNUC__)
337 : #pragma GCC diagnostic push
338 : #pragma GCC diagnostic ignored "-Wnull-dereference"
339 : #endif
340 3 : newCount[0] = 1;
341 3 : newArrayStep[0] = 1;
342 : #if defined(__GNUC__)
343 : #pragma GCC diagnostic pop
344 : #endif
345 9 : for (size_t i = 0; i < count[0]; ++i)
346 : {
347 6 : if (i > 0)
348 : {
349 3 : if (arrayStep[0] > 0)
350 2 : arrayIdx += arrayStep[0];
351 : else
352 1 : arrayIdx -= static_cast<GUInt64>(-arrayStep[0]);
353 3 : pabyDstBuffer += bufferStride[0] * nBufferDTSize;
354 : }
355 6 : newArrayStartIdx[0] =
356 6 : m_poShared->m_anMapNewDimToOldDim[static_cast<int>(arrayIdx)];
357 12 : if (!m_poParent->Read(newArrayStartIdx.data(), newCount.data(),
358 6 : newArrayStep.data(), bufferStride,
359 : bufferDataType, pabyDstBuffer))
360 : {
361 0 : return false;
362 : }
363 : }
364 3 : return true;
365 : }
366 :
367 : // Slow/unoptimized case
368 50 : std::vector<size_t> anStackIter(nDims);
369 50 : std::vector<GUInt64> anStackArrayIdx(nDims);
370 50 : std::vector<GByte *> pabyDstBufferStack(nDims + 1);
371 : #if defined(__GNUC__)
372 : #pragma GCC diagnostic push
373 : #pragma GCC diagnostic ignored "-Wnull-dereference"
374 : #endif
375 25 : pabyDstBufferStack[0] = static_cast<GByte *>(pDstBuffer);
376 : #if defined(__GNUC__)
377 : #pragma GCC diagnostic pop
378 : #endif
379 25 : size_t iDim = 0;
380 183 : lbl_next_depth:
381 183 : if (iDim == nDims)
382 : {
383 160 : if (!m_poParent->Read(newArrayStartIdx.data(), newCount.data(),
384 80 : newArrayStep.data(), bufferStride, bufferDataType,
385 80 : pabyDstBufferStack[iDim]))
386 : {
387 0 : return false;
388 : }
389 : }
390 : else
391 : {
392 103 : anStackIter[iDim] = 0;
393 103 : anStackArrayIdx[iDim] = arrayStartIdx[iDim];
394 : while (true)
395 : {
396 158 : if (m_abPatchedDim[iDim])
397 : {
398 84 : newArrayStartIdx[iDim] =
399 84 : m_poShared->m_anMapNewDimToOldDim[static_cast<int>(
400 84 : anStackArrayIdx[iDim])];
401 : }
402 : else
403 : {
404 74 : newArrayStartIdx[iDim] = anStackArrayIdx[iDim];
405 : }
406 158 : ++iDim;
407 158 : pabyDstBufferStack[iDim] = pabyDstBufferStack[iDim - 1];
408 158 : goto lbl_next_depth;
409 158 : lbl_return_to_caller_in_loop:
410 158 : --iDim;
411 158 : ++anStackIter[iDim];
412 158 : if (anStackIter[iDim] == count[iDim])
413 103 : break;
414 55 : if (arrayStep[iDim] > 0)
415 54 : anStackArrayIdx[iDim] += arrayStep[iDim];
416 : else
417 1 : anStackArrayIdx[iDim] -= -arrayStep[iDim];
418 55 : pabyDstBufferStack[iDim] += bufferStride[iDim] * nBufferDTSize;
419 : }
420 : }
421 183 : if (iDim > 0)
422 158 : goto lbl_return_to_caller_in_loop;
423 :
424 25 : return true;
425 : }
426 :
427 : /************************************************************************/
428 : /* SubsetDimensionFromSelection() */
429 : /************************************************************************/
430 :
431 : /** Return a virtual group whose one dimension has been subset according to a
432 : * selection.
433 : *
434 : * The selection criterion is currently restricted to the form
435 : * "/path/to/array=numeric_value" (no spaces around equal)
436 : *
437 : * This is similar to XArray indexing by name and label on a XArray Dataset
438 : * using the sel() method.
439 : * Cf https://docs.xarray.dev/en/latest/user-guide/indexing.html#quick-overview
440 : *
441 : * For example on a EMIT L2A product
442 : * (https://github.com/nasa/EMIT-Data-Resources/blob/main/python/tutorials/Exploring_EMIT_L2A_Reflectance.ipynb),
443 : * this can be used to keep only valid bands with
444 : * SubsetDimensionFromSelection("/sensor_band_parameters/good_wavelengths=1")
445 : *
446 : * This is the same as the C function GDALGroupSubsetDimensionFromSelection().
447 : *
448 : * @param osSelection Selection criterion.
449 : * @return a virtual group, or nullptr in case of error
450 : * @since 3.8
451 : */
452 : std::shared_ptr<GDALGroup>
453 14 : GDALGroup::SubsetDimensionFromSelection(const std::string &osSelection) const
454 : {
455 28 : auto self = std::dynamic_pointer_cast<GDALGroup>(m_pSelf.lock());
456 14 : if (!self)
457 : {
458 0 : CPLError(CE_Failure, CPLE_AppDefined,
459 : "Driver implementation issue: m_pSelf not set !");
460 0 : return nullptr;
461 : }
462 :
463 14 : const auto nEqualPos = osSelection.find('=');
464 14 : if (nEqualPos == std::string::npos)
465 : {
466 2 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for selection");
467 2 : return nullptr;
468 : }
469 24 : const auto osArrayName = osSelection.substr(0, nEqualPos);
470 24 : const auto osValue = osSelection.substr(nEqualPos + 1);
471 13 : if (CPLGetValueType(osValue.c_str()) != CPL_VALUE_INTEGER &&
472 1 : CPLGetValueType(osValue.c_str()) != CPL_VALUE_REAL)
473 : {
474 1 : CPLError(CE_Failure, CPLE_AppDefined,
475 : "Non-numeric value in selection criterion");
476 1 : return nullptr;
477 : }
478 22 : auto poArray = OpenMDArrayFromFullname(osArrayName);
479 11 : if (!poArray)
480 : {
481 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find array %s",
482 : osArrayName.c_str());
483 1 : return nullptr;
484 : }
485 10 : if (poArray->GetDimensionCount() != 1)
486 : {
487 1 : CPLError(CE_Failure, CPLE_AppDefined,
488 : "Array %s is not single dimensional", osArrayName.c_str());
489 1 : return nullptr;
490 : }
491 9 : if (poArray->GetDataType().GetClass() != GEDTC_NUMERIC)
492 : {
493 1 : CPLError(CE_Failure, CPLE_AppDefined, "Array %s is not of numeric type",
494 : osArrayName.c_str());
495 1 : return nullptr;
496 : }
497 :
498 8 : const auto nElts = poArray->GetTotalElementsCount();
499 8 : if (nElts > 10 * 1024 * 1024)
500 : {
501 1 : CPLError(CE_Failure, CPLE_AppDefined, "Too many values in %s",
502 : osArrayName.c_str());
503 1 : return nullptr;
504 : }
505 14 : std::vector<double> values;
506 : try
507 : {
508 7 : values.resize(static_cast<size_t>(nElts));
509 : }
510 0 : catch (const std::bad_alloc &e)
511 : {
512 0 : CPLError(CE_Failure, CPLE_AppDefined, "Out of memory: %s", e.what());
513 0 : return nullptr;
514 : }
515 7 : const GUInt64 startIdx[1] = {0};
516 7 : const size_t count[1] = {values.size()};
517 14 : if (!poArray->Read(startIdx, count, nullptr, nullptr,
518 14 : GDALExtendedDataType::Create(GDT_Float64), &values[0],
519 7 : values.data(), values.size() * sizeof(values[0])))
520 : {
521 0 : return nullptr;
522 : }
523 7 : const double dfSelectionValue = CPLAtof(osValue.c_str());
524 14 : std::vector<int> anMapNewDimToOldDim;
525 77 : for (int i = 0; i < static_cast<int>(nElts); ++i)
526 : {
527 70 : if (values[i] == dfSelectionValue)
528 8 : anMapNewDimToOldDim.push_back(i);
529 : }
530 7 : if (anMapNewDimToOldDim.empty())
531 : {
532 1 : CPLError(CE_Failure, CPLE_AppDefined, "No value in %s matching %f",
533 : osArrayName.c_str(), dfSelectionValue);
534 1 : return nullptr;
535 : }
536 6 : if (anMapNewDimToOldDim.size() == nElts)
537 : {
538 1 : return self;
539 : }
540 :
541 10 : auto poDim = poArray->GetDimensions()[0];
542 10 : auto poShared = std::make_shared<GDALSubsetGroupSharedResources>();
543 5 : if (GetFullName() == "/")
544 5 : poShared->m_poRootGroup = self;
545 5 : poShared->m_osSelection = osSelection;
546 5 : poShared->m_osDimFullName = poArray->GetDimensions()[0]->GetFullName();
547 5 : poShared->m_anMapNewDimToOldDim = std::move(anMapNewDimToOldDim);
548 :
549 : // Create a modified dimension of reduced size
550 : auto poNewDim = std::make_shared<GDALDimensionWeakIndexingVar>(
551 5 : GetParentName(poDim->GetFullName()), poDim->GetName(), poDim->GetType(),
552 15 : poDim->GetDirection(), poShared->m_anMapNewDimToOldDim.size());
553 5 : poShared->m_poNewDim = poNewDim;
554 :
555 10 : auto poIndexingVar = poDim->GetIndexingVariable();
556 5 : if (poIndexingVar)
557 : {
558 : // poNewIndexingVar must be created with a different GDALSubsetGroupSharedResources
559 : // instance than poShared, to avoid cross reference, that would result in
560 : // objects not being freed !
561 : auto poSpecificShared =
562 10 : std::make_shared<GDALSubsetGroupSharedResources>();
563 5 : poSpecificShared->m_poRootGroup = poShared->m_poRootGroup;
564 5 : poSpecificShared->m_osSelection = osSelection;
565 5 : poSpecificShared->m_osDimFullName =
566 10 : poArray->GetDimensions()[0]->GetFullName();
567 5 : poSpecificShared->m_anMapNewDimToOldDim =
568 5 : poShared->m_anMapNewDimToOldDim;
569 5 : poSpecificShared->m_poNewDim = poNewDim;
570 : auto poNewIndexingVar =
571 : GDALSubsetArray::Create(poIndexingVar, poSpecificShared,
572 10 : CreateContext(GetContext(), poShared));
573 5 : poNewDim->SetIndexingVariable(poNewIndexingVar);
574 5 : poShared->m_poNewIndexingVar = poNewIndexingVar;
575 : }
576 :
577 5 : return GDALSubsetGroup::Create(self, poShared);
578 : }
|
