Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: gdalmultidim_array_view.cpp
4 : * Project: GDAL Core
5 : * Purpose: GDALMDArray::GetView() implementation
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_multidim.h"
16 : #include "ogr_spatialref.h"
17 :
18 : #include <limits>
19 :
20 : //! @cond Doxygen_Suppress
21 :
22 : /************************************************************************/
23 : /* GDALSlicedMDArray */
24 : /************************************************************************/
25 :
26 : class GDALSlicedMDArray final : public GDALPamMDArray
27 : {
28 : private:
29 : std::shared_ptr<GDALMDArray> m_poParent{};
30 : std::vector<std::shared_ptr<GDALDimension>> m_dims{};
31 : std::vector<size_t> m_mapDimIdxToParentDimIdx{}; // of size m_dims.size()
32 : std::vector<std::shared_ptr<GDALMDArray>> m_apoNewIndexingVariables{};
33 : std::vector<Range>
34 : m_parentRanges{}; // of size m_poParent->GetDimensionCount()
35 :
36 : mutable std::vector<GUInt64> m_parentStart;
37 : mutable std::vector<size_t> m_parentCount;
38 : mutable std::vector<GInt64> m_parentStep;
39 : mutable std::vector<GPtrDiff_t> m_parentStride;
40 :
41 : void PrepareParentArrays(const GUInt64 *arrayStartIdx, const size_t *count,
42 : const GInt64 *arrayStep,
43 : const GPtrDiff_t *bufferStride) const;
44 :
45 : protected:
46 862 : explicit GDALSlicedMDArray(
47 : const std::shared_ptr<GDALMDArray> &poParent,
48 : const std::string &viewExpr,
49 : std::vector<std::shared_ptr<GDALDimension>> &&dims,
50 : std::vector<size_t> &&mapDimIdxToParentDimIdx,
51 : std::vector<std::shared_ptr<GDALMDArray>> &&apoNewIndexingVariables,
52 : std::vector<Range> &&parentRanges)
53 2586 : : GDALAbstractMDArray(std::string(), "Sliced view of " +
54 2586 : poParent->GetFullName() +
55 1724 : " (" + viewExpr + ")"),
56 1724 : GDALPamMDArray(std::string(),
57 1724 : "Sliced view of " + poParent->GetFullName() + " (" +
58 1724 : viewExpr + ")",
59 1724 : GDALPamMultiDim::GetPAM(poParent),
60 : poParent->GetContext()),
61 1724 : m_poParent(std::move(poParent)), m_dims(std::move(dims)),
62 862 : m_mapDimIdxToParentDimIdx(std::move(mapDimIdxToParentDimIdx)),
63 862 : m_apoNewIndexingVariables(std::move(apoNewIndexingVariables)),
64 862 : m_parentRanges(std::move(parentRanges)),
65 862 : m_parentStart(m_poParent->GetDimensionCount()),
66 862 : m_parentCount(m_poParent->GetDimensionCount(), 1),
67 862 : m_parentStep(m_poParent->GetDimensionCount()),
68 6896 : m_parentStride(m_poParent->GetDimensionCount())
69 : {
70 862 : }
71 :
72 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
73 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
74 : const GDALExtendedDataType &bufferDataType,
75 : void *pDstBuffer) const override;
76 :
77 : bool IWrite(const GUInt64 *arrayStartIdx, const size_t *count,
78 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
79 : const GDALExtendedDataType &bufferDataType,
80 : const void *pSrcBuffer) override;
81 :
82 : bool IAdviseRead(const GUInt64 *arrayStartIdx, const size_t *count,
83 : CSLConstList papszOptions) const override;
84 :
85 : public:
86 : static std::shared_ptr<GDALSlicedMDArray>
87 862 : Create(const std::shared_ptr<GDALMDArray> &poParent,
88 : const std::string &viewExpr,
89 : std::vector<std::shared_ptr<GDALDimension>> &&dims,
90 : std::vector<size_t> &&mapDimIdxToParentDimIdx,
91 : std::vector<std::shared_ptr<GDALMDArray>> &&apoNewIndexingVariables,
92 : std::vector<Range> &&parentRanges)
93 : {
94 862 : CPLAssert(dims.size() == mapDimIdxToParentDimIdx.size());
95 862 : CPLAssert(parentRanges.size() == poParent->GetDimensionCount());
96 :
97 : auto newAr(std::shared_ptr<GDALSlicedMDArray>(new GDALSlicedMDArray(
98 862 : poParent, viewExpr, std::move(dims),
99 862 : std::move(mapDimIdxToParentDimIdx),
100 862 : std::move(apoNewIndexingVariables), std::move(parentRanges))));
101 862 : newAr->SetSelf(newAr);
102 862 : return newAr;
103 : }
104 :
105 228 : bool IsWritable() const override
106 : {
107 228 : return m_poParent->IsWritable();
108 : }
109 :
110 1614 : const std::string &GetFilename() const override
111 : {
112 1614 : return m_poParent->GetFilename();
113 : }
114 :
115 : const std::vector<std::shared_ptr<GDALDimension>> &
116 6215 : GetDimensions() const override
117 : {
118 6215 : return m_dims;
119 : }
120 :
121 2445 : const GDALExtendedDataType &GetDataType() const override
122 : {
123 2445 : return m_poParent->GetDataType();
124 : }
125 :
126 4 : const std::string &GetUnit() const override
127 : {
128 4 : return m_poParent->GetUnit();
129 : }
130 :
131 : // bool SetUnit(const std::string& osUnit) override { return
132 : // m_poParent->SetUnit(osUnit); }
133 :
134 2 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
135 : {
136 4 : auto poSrcSRS = m_poParent->GetSpatialRef();
137 2 : if (!poSrcSRS)
138 1 : return nullptr;
139 2 : auto srcMapping = poSrcSRS->GetDataAxisToSRSAxisMapping();
140 2 : std::vector<int> dstMapping;
141 3 : for (int srcAxis : srcMapping)
142 : {
143 2 : bool bFound = false;
144 3 : for (size_t i = 0; i < m_mapDimIdxToParentDimIdx.size(); i++)
145 : {
146 3 : if (static_cast<int>(m_mapDimIdxToParentDimIdx[i]) ==
147 3 : srcAxis - 1)
148 : {
149 2 : dstMapping.push_back(static_cast<int>(i) + 1);
150 2 : bFound = true;
151 2 : break;
152 : }
153 : }
154 2 : if (!bFound)
155 : {
156 0 : dstMapping.push_back(0);
157 : }
158 : }
159 2 : auto poClone(std::shared_ptr<OGRSpatialReference>(poSrcSRS->Clone()));
160 1 : poClone->SetDataAxisToSRSAxisMapping(dstMapping);
161 1 : return poClone;
162 : }
163 :
164 104 : const void *GetRawNoDataValue() const override
165 : {
166 104 : return m_poParent->GetRawNoDataValue();
167 : }
168 :
169 : // bool SetRawNoDataValue(const void* pRawNoData) override { return
170 : // m_poParent->SetRawNoDataValue(pRawNoData); }
171 :
172 4 : double GetOffset(bool *pbHasOffset,
173 : GDALDataType *peStorageType) const override
174 : {
175 4 : return m_poParent->GetOffset(pbHasOffset, peStorageType);
176 : }
177 :
178 4 : double GetScale(bool *pbHasScale,
179 : GDALDataType *peStorageType) const override
180 : {
181 4 : return m_poParent->GetScale(pbHasScale, peStorageType);
182 : }
183 :
184 : // bool SetOffset(double dfOffset) override { return
185 : // m_poParent->SetOffset(dfOffset); }
186 :
187 : // bool SetScale(double dfScale) override { return
188 : // m_poParent->SetScale(dfScale); }
189 :
190 580 : std::vector<GUInt64> GetBlockSize() const override
191 : {
192 580 : std::vector<GUInt64> ret(GetDimensionCount());
193 1160 : const auto parentBlockSize(m_poParent->GetBlockSize());
194 1564 : for (size_t i = 0; i < m_mapDimIdxToParentDimIdx.size(); ++i)
195 : {
196 984 : const auto iOldAxis = m_mapDimIdxToParentDimIdx[i];
197 984 : if (iOldAxis != static_cast<size_t>(-1))
198 : {
199 984 : ret[i] = parentBlockSize[iOldAxis];
200 : }
201 : }
202 1160 : return ret;
203 : }
204 :
205 : std::shared_ptr<GDALAttribute>
206 14 : GetAttribute(const std::string &osName) const override
207 : {
208 14 : return m_poParent->GetAttribute(osName);
209 : }
210 :
211 : std::vector<std::shared_ptr<GDALAttribute>>
212 37 : GetAttributes(CSLConstList papszOptions = nullptr) const override
213 : {
214 37 : return m_poParent->GetAttributes(papszOptions);
215 : }
216 : };
217 :
218 : /************************************************************************/
219 : /* PrepareParentArrays() */
220 : /************************************************************************/
221 :
222 763 : void GDALSlicedMDArray::PrepareParentArrays(
223 : const GUInt64 *arrayStartIdx, const size_t *count, const GInt64 *arrayStep,
224 : const GPtrDiff_t *bufferStride) const
225 : {
226 763 : const size_t nParentDimCount = m_parentRanges.size();
227 2148 : for (size_t i = 0; i < nParentDimCount; i++)
228 : {
229 : // For dimensions in parent that have no existence in sliced array
230 1385 : m_parentStart[i] = m_parentRanges[i].m_nStartIdx;
231 : }
232 :
233 1919 : for (size_t i = 0; i < m_dims.size(); i++)
234 : {
235 1156 : const auto iParent = m_mapDimIdxToParentDimIdx[i];
236 1156 : if (iParent != static_cast<size_t>(-1))
237 : {
238 1154 : m_parentStart[iParent] =
239 1154 : m_parentRanges[iParent].m_nIncr >= 0
240 1154 : ? m_parentRanges[iParent].m_nStartIdx +
241 852 : arrayStartIdx[i] * m_parentRanges[iParent].m_nIncr
242 302 : : m_parentRanges[iParent].m_nStartIdx -
243 604 : arrayStartIdx[i] *
244 302 : static_cast<GUInt64>(
245 302 : -m_parentRanges[iParent].m_nIncr);
246 1154 : m_parentCount[iParent] = count[i];
247 1154 : if (arrayStep)
248 : {
249 1151 : m_parentStep[iParent] =
250 1151 : count[i] == 1 ? 1 :
251 : // other checks should have ensured this does
252 : // not overflow
253 947 : arrayStep[i] * m_parentRanges[iParent].m_nIncr;
254 : }
255 1154 : if (bufferStride)
256 : {
257 1151 : m_parentStride[iParent] = bufferStride[i];
258 : }
259 : }
260 : }
261 763 : }
262 :
263 : /************************************************************************/
264 : /* IRead() */
265 : /************************************************************************/
266 :
267 732 : bool GDALSlicedMDArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
268 : const GInt64 *arrayStep,
269 : const GPtrDiff_t *bufferStride,
270 : const GDALExtendedDataType &bufferDataType,
271 : void *pDstBuffer) const
272 : {
273 732 : PrepareParentArrays(arrayStartIdx, count, arrayStep, bufferStride);
274 1464 : return m_poParent->Read(m_parentStart.data(), m_parentCount.data(),
275 732 : m_parentStep.data(), m_parentStride.data(),
276 732 : bufferDataType, pDstBuffer);
277 : }
278 :
279 : /************************************************************************/
280 : /* IWrite() */
281 : /************************************************************************/
282 :
283 29 : bool GDALSlicedMDArray::IWrite(const GUInt64 *arrayStartIdx,
284 : const size_t *count, const GInt64 *arrayStep,
285 : const GPtrDiff_t *bufferStride,
286 : const GDALExtendedDataType &bufferDataType,
287 : const void *pSrcBuffer)
288 : {
289 29 : PrepareParentArrays(arrayStartIdx, count, arrayStep, bufferStride);
290 58 : return m_poParent->Write(m_parentStart.data(), m_parentCount.data(),
291 29 : m_parentStep.data(), m_parentStride.data(),
292 29 : bufferDataType, pSrcBuffer);
293 : }
294 :
295 : /************************************************************************/
296 : /* IAdviseRead() */
297 : /************************************************************************/
298 :
299 2 : bool GDALSlicedMDArray::IAdviseRead(const GUInt64 *arrayStartIdx,
300 : const size_t *count,
301 : CSLConstList papszOptions) const
302 : {
303 2 : PrepareParentArrays(arrayStartIdx, count, nullptr, nullptr);
304 2 : return m_poParent->AdviseRead(m_parentStart.data(), m_parentCount.data(),
305 2 : papszOptions);
306 : }
307 :
308 : /************************************************************************/
309 : /* CreateSlicedArray() */
310 : /************************************************************************/
311 :
312 : static std::shared_ptr<GDALMDArray>
313 724 : CreateSlicedArray(const std::shared_ptr<GDALMDArray> &self,
314 : const std::string &viewExpr, const std::string &activeSlice,
315 : bool bRenameDimensions,
316 : std::vector<GDALMDArray::ViewSpec> &viewSpecs)
317 : {
318 724 : const auto &srcDims(self->GetDimensions());
319 724 : if (srcDims.empty())
320 : {
321 2 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot slice a 0-d array");
322 2 : return nullptr;
323 : }
324 :
325 1444 : CPLStringList aosTokens(CSLTokenizeString2(activeSlice.c_str(), ",", 0));
326 722 : const auto nTokens = static_cast<size_t>(aosTokens.size());
327 :
328 1444 : std::vector<std::shared_ptr<GDALDimension>> newDims;
329 1444 : std::vector<size_t> mapDimIdxToParentDimIdx;
330 1444 : std::vector<GDALSlicedMDArray::Range> parentRanges;
331 722 : newDims.reserve(nTokens);
332 722 : mapDimIdxToParentDimIdx.reserve(nTokens);
333 722 : parentRanges.reserve(nTokens);
334 :
335 722 : bool bGotEllipsis = false;
336 722 : size_t nCurSrcDim = 0;
337 1444 : std::vector<std::shared_ptr<GDALMDArray>> apoNewIndexingVariables;
338 2151 : for (size_t i = 0; i < nTokens; i++)
339 : {
340 1446 : const char *pszIdxSpec = aosTokens[i];
341 1446 : if (EQUAL(pszIdxSpec, "..."))
342 : {
343 129 : if (bGotEllipsis)
344 : {
345 2 : CPLError(CE_Failure, CPLE_AppDefined,
346 : "Only one single ellipsis is supported");
347 2 : return nullptr;
348 : }
349 127 : bGotEllipsis = true;
350 127 : const auto nSubstitutionCount = srcDims.size() - (nTokens - 1);
351 263 : for (size_t j = 0; j < nSubstitutionCount; j++, nCurSrcDim++)
352 : {
353 136 : parentRanges.emplace_back(0, 1);
354 136 : newDims.push_back(srcDims[nCurSrcDim]);
355 136 : mapDimIdxToParentDimIdx.push_back(nCurSrcDim);
356 : }
357 127 : continue;
358 : }
359 1317 : else if (EQUAL(pszIdxSpec, "newaxis") ||
360 1314 : EQUAL(pszIdxSpec, "np.newaxis"))
361 : {
362 3 : newDims.push_back(std::make_shared<GDALDimension>(
363 6 : std::string(), "newaxis", std::string(), std::string(), 1));
364 3 : mapDimIdxToParentDimIdx.push_back(static_cast<size_t>(-1));
365 3 : continue;
366 : }
367 1314 : else if (CPLGetValueType(pszIdxSpec) == CPL_VALUE_INTEGER)
368 : {
369 342 : if (nCurSrcDim >= srcDims.size())
370 : {
371 2 : CPLError(CE_Failure, CPLE_AppDefined, "Too many values in %s",
372 : activeSlice.c_str());
373 7 : return nullptr;
374 : }
375 :
376 340 : auto nVal = CPLAtoGIntBig(pszIdxSpec);
377 340 : GUInt64 nDimSize = srcDims[nCurSrcDim]->GetSize();
378 340 : if ((nVal >= 0 && static_cast<GUInt64>(nVal) >= nDimSize) ||
379 336 : (nVal < 0 && nDimSize < static_cast<GUInt64>(-nVal)))
380 : {
381 5 : CPLError(CE_Failure, CPLE_AppDefined,
382 : "Index " CPL_FRMT_GIB " is out of bounds", nVal);
383 5 : return nullptr;
384 : }
385 335 : if (nVal < 0)
386 0 : nVal += nDimSize;
387 335 : parentRanges.emplace_back(nVal, 0);
388 : }
389 : else
390 : {
391 972 : if (nCurSrcDim >= srcDims.size())
392 : {
393 1 : CPLError(CE_Failure, CPLE_AppDefined, "Too many values in %s",
394 : activeSlice.c_str());
395 8 : return nullptr;
396 : }
397 :
398 : CPLStringList aosRangeTokens(
399 971 : CSLTokenizeString2(pszIdxSpec, ":", CSLT_ALLOWEMPTYTOKENS));
400 971 : int nRangeTokens = aosRangeTokens.size();
401 971 : if (nRangeTokens > 3)
402 : {
403 1 : CPLError(CE_Failure, CPLE_AppDefined, "Too many : in %s",
404 : pszIdxSpec);
405 1 : return nullptr;
406 : }
407 970 : if (nRangeTokens <= 1)
408 : {
409 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value %s",
410 : pszIdxSpec);
411 1 : return nullptr;
412 : }
413 969 : const char *pszStart = aosRangeTokens[0];
414 969 : const char *pszEnd = aosRangeTokens[1];
415 969 : const char *pszInc = (nRangeTokens == 3) ? aosRangeTokens[2] : "";
416 969 : GDALSlicedMDArray::Range range;
417 969 : const GUInt64 nDimSize(srcDims[nCurSrcDim]->GetSize());
418 969 : range.m_nIncr = EQUAL(pszInc, "") ? 1 : CPLAtoGIntBig(pszInc);
419 1937 : if (range.m_nIncr == 0 ||
420 968 : range.m_nIncr == std::numeric_limits<GInt64>::min())
421 : {
422 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid increment");
423 1 : return nullptr;
424 : }
425 968 : auto startIdx(CPLAtoGIntBig(pszStart));
426 968 : if (startIdx < 0)
427 : {
428 0 : if (nDimSize < static_cast<GUInt64>(-startIdx))
429 0 : startIdx = 0;
430 : else
431 0 : startIdx = nDimSize + startIdx;
432 : }
433 968 : const bool bPosIncr = range.m_nIncr > 0;
434 968 : range.m_nStartIdx = startIdx;
435 1936 : range.m_nStartIdx = EQUAL(pszStart, "")
436 968 : ? (bPosIncr ? 0 : nDimSize - 1)
437 : : range.m_nStartIdx;
438 968 : if (range.m_nStartIdx >= nDimSize - 1)
439 286 : range.m_nStartIdx = nDimSize - 1;
440 968 : auto endIdx(CPLAtoGIntBig(pszEnd));
441 968 : if (endIdx < 0)
442 : {
443 1 : const auto positiveEndIdx = static_cast<GUInt64>(-endIdx);
444 1 : if (nDimSize < positiveEndIdx)
445 0 : endIdx = 0;
446 : else
447 1 : endIdx = nDimSize - positiveEndIdx;
448 : }
449 968 : GUInt64 nEndIdx = endIdx;
450 968 : nEndIdx = EQUAL(pszEnd, "") ? (!bPosIncr ? 0 : nDimSize) : nEndIdx;
451 968 : if (pszStart[0] || pszEnd[0])
452 : {
453 636 : if ((bPosIncr && range.m_nStartIdx >= nEndIdx) ||
454 633 : (!bPosIncr && range.m_nStartIdx <= nEndIdx))
455 : {
456 4 : CPLError(CE_Failure, CPLE_AppDefined,
457 : "Output dimension of size 0 is not allowed");
458 4 : return nullptr;
459 : }
460 : }
461 964 : int inc = (EQUAL(pszEnd, "") && !bPosIncr) ? 1 : 0;
462 964 : const auto nAbsIncr = std::abs(range.m_nIncr);
463 : const GUInt64 newSize =
464 332 : (pszStart[0] == 0 && pszEnd[0] == 0 &&
465 332 : range.m_nStartIdx == nEndIdx)
466 1928 : ? 1
467 : : bPosIncr
468 963 : ? cpl::div_round_up(nEndIdx - range.m_nStartIdx, nAbsIncr)
469 128 : : cpl::div_round_up(inc + range.m_nStartIdx - nEndIdx,
470 964 : nAbsIncr);
471 964 : const auto &poSrcDim = srcDims[nCurSrcDim];
472 1544 : if (range.m_nStartIdx == 0 && range.m_nIncr == 1 &&
473 580 : newSize == poSrcDim->GetSize())
474 : {
475 211 : newDims.push_back(poSrcDim);
476 : }
477 : else
478 : {
479 1506 : std::string osNewDimName(poSrcDim->GetName());
480 753 : if (bRenameDimensions)
481 : {
482 : osNewDimName =
483 1410 : "subset_" + poSrcDim->GetName() +
484 : CPLSPrintf("_" CPL_FRMT_GUIB "_" CPL_FRMT_GIB
485 : "_" CPL_FRMT_GUIB,
486 705 : static_cast<GUIntBig>(range.m_nStartIdx),
487 705 : static_cast<GIntBig>(range.m_nIncr),
488 705 : static_cast<GUIntBig>(newSize));
489 : }
490 : auto poNewDim = std::make_shared<GDALDimensionWeakIndexingVar>(
491 1506 : std::string(), osNewDimName, poSrcDim->GetType(),
492 753 : range.m_nIncr > 0 ? poSrcDim->GetDirection()
493 : : std::string(),
494 1506 : newSize);
495 753 : auto poSrcIndexingVar = poSrcDim->GetIndexingVariable();
496 910 : if (poSrcIndexingVar &&
497 910 : poSrcIndexingVar->GetDimensionCount() == 1 &&
498 157 : poSrcIndexingVar->GetDimensions()[0] == poSrcDim)
499 : {
500 : std::vector<std::shared_ptr<GDALDimension>>
501 628 : indexingVarNewDims{poNewDim};
502 314 : std::vector<size_t> indexingVarMapDimIdxToParentDimIdx{0};
503 : std::vector<std::shared_ptr<GDALMDArray>>
504 314 : indexingVarNewIndexingVar;
505 : std::vector<GDALSlicedMDArray::Range>
506 314 : indexingVarParentRanges{range};
507 : auto poNewIndexingVar = GDALSlicedMDArray::Create(
508 : poSrcIndexingVar, pszIdxSpec,
509 157 : std::move(indexingVarNewDims),
510 157 : std::move(indexingVarMapDimIdxToParentDimIdx),
511 157 : std::move(indexingVarNewIndexingVar),
512 471 : std::move(indexingVarParentRanges));
513 157 : poNewDim->SetIndexingVariable(poNewIndexingVar);
514 157 : apoNewIndexingVariables.push_back(
515 157 : std::move(poNewIndexingVar));
516 : }
517 753 : newDims.push_back(std::move(poNewDim));
518 : }
519 964 : mapDimIdxToParentDimIdx.push_back(nCurSrcDim);
520 964 : parentRanges.emplace_back(range);
521 : }
522 :
523 1299 : nCurSrcDim++;
524 : }
525 778 : for (; nCurSrcDim < srcDims.size(); nCurSrcDim++)
526 : {
527 73 : parentRanges.emplace_back(0, 1);
528 73 : newDims.push_back(srcDims[nCurSrcDim]);
529 73 : mapDimIdxToParentDimIdx.push_back(nCurSrcDim);
530 : }
531 :
532 705 : GDALMDArray::ViewSpec viewSpec;
533 705 : viewSpec.m_mapDimIdxToParentDimIdx = mapDimIdxToParentDimIdx;
534 705 : viewSpec.m_parentRanges = parentRanges;
535 705 : viewSpecs.emplace_back(std::move(viewSpec));
536 :
537 1410 : return GDALSlicedMDArray::Create(
538 705 : self, viewExpr, std::move(newDims), std::move(mapDimIdxToParentDimIdx),
539 1410 : std::move(apoNewIndexingVariables), std::move(parentRanges));
540 : }
541 :
542 : /************************************************************************/
543 : /* GDALExtractFieldMDArray */
544 : /************************************************************************/
545 :
546 : class GDALExtractFieldMDArray final : public GDALPamMDArray
547 : {
548 : private:
549 : std::shared_ptr<GDALMDArray> m_poParent{};
550 : GDALExtendedDataType m_dt;
551 : std::string m_srcCompName;
552 : mutable std::vector<GByte> m_pabyNoData{};
553 :
554 : protected:
555 465 : GDALExtractFieldMDArray(const std::shared_ptr<GDALMDArray> &poParent,
556 : const std::string &fieldName,
557 : const std::unique_ptr<GDALEDTComponent> &srcComp)
558 1860 : : GDALAbstractMDArray(std::string(), "Extract field " + fieldName +
559 930 : " of " +
560 465 : poParent->GetFullName()),
561 : GDALPamMDArray(
562 930 : std::string(),
563 930 : "Extract field " + fieldName + " of " + poParent->GetFullName(),
564 930 : GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
565 : m_poParent(poParent), m_dt(srcComp->GetType()),
566 2325 : m_srcCompName(srcComp->GetName())
567 : {
568 465 : m_pabyNoData.resize(m_dt.GetSize());
569 465 : }
570 :
571 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
572 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
573 : const GDALExtendedDataType &bufferDataType,
574 : void *pDstBuffer) const override;
575 :
576 1 : bool IAdviseRead(const GUInt64 *arrayStartIdx, const size_t *count,
577 : CSLConstList papszOptions) const override
578 : {
579 1 : return m_poParent->AdviseRead(arrayStartIdx, count, papszOptions);
580 : }
581 :
582 : public:
583 : static std::shared_ptr<GDALExtractFieldMDArray>
584 465 : Create(const std::shared_ptr<GDALMDArray> &poParent,
585 : const std::string &fieldName,
586 : const std::unique_ptr<GDALEDTComponent> &srcComp)
587 : {
588 : auto newAr(std::shared_ptr<GDALExtractFieldMDArray>(
589 465 : new GDALExtractFieldMDArray(poParent, fieldName, srcComp)));
590 465 : newAr->SetSelf(newAr);
591 465 : return newAr;
592 : }
593 :
594 930 : ~GDALExtractFieldMDArray() override
595 465 : {
596 465 : m_dt.FreeDynamicMemory(&m_pabyNoData[0]);
597 930 : }
598 :
599 200 : bool IsWritable() const override
600 : {
601 200 : return m_poParent->IsWritable();
602 : }
603 :
604 1116 : const std::string &GetFilename() const override
605 : {
606 1116 : return m_poParent->GetFilename();
607 : }
608 :
609 : const std::vector<std::shared_ptr<GDALDimension>> &
610 1583 : GetDimensions() const override
611 : {
612 1583 : return m_poParent->GetDimensions();
613 : }
614 :
615 1137 : const GDALExtendedDataType &GetDataType() const override
616 : {
617 1137 : return m_dt;
618 : }
619 :
620 2 : const std::string &GetUnit() const override
621 : {
622 2 : return m_poParent->GetUnit();
623 : }
624 :
625 2 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
626 : {
627 2 : return m_poParent->GetSpatialRef();
628 : }
629 :
630 99 : const void *GetRawNoDataValue() const override
631 : {
632 99 : const void *parentNoData = m_poParent->GetRawNoDataValue();
633 99 : if (parentNoData == nullptr)
634 6 : return nullptr;
635 :
636 93 : m_dt.FreeDynamicMemory(&m_pabyNoData[0]);
637 93 : memset(&m_pabyNoData[0], 0, m_dt.GetSize());
638 :
639 186 : std::vector<std::unique_ptr<GDALEDTComponent>> comps;
640 186 : comps.emplace_back(std::unique_ptr<GDALEDTComponent>(
641 186 : new GDALEDTComponent(m_srcCompName, 0, m_dt)));
642 93 : auto tmpDT(GDALExtendedDataType::Create(std::string(), m_dt.GetSize(),
643 279 : std::move(comps)));
644 :
645 93 : GDALExtendedDataType::CopyValue(parentNoData, m_poParent->GetDataType(),
646 93 : &m_pabyNoData[0], tmpDT);
647 :
648 93 : return &m_pabyNoData[0];
649 : }
650 :
651 2 : double GetOffset(bool *pbHasOffset,
652 : GDALDataType *peStorageType) const override
653 : {
654 2 : return m_poParent->GetOffset(pbHasOffset, peStorageType);
655 : }
656 :
657 2 : double GetScale(bool *pbHasScale,
658 : GDALDataType *peStorageType) const override
659 : {
660 2 : return m_poParent->GetScale(pbHasScale, peStorageType);
661 : }
662 :
663 201 : std::vector<GUInt64> GetBlockSize() const override
664 : {
665 201 : return m_poParent->GetBlockSize();
666 : }
667 : };
668 :
669 : /************************************************************************/
670 : /* IRead() */
671 : /************************************************************************/
672 :
673 414 : bool GDALExtractFieldMDArray::IRead(const GUInt64 *arrayStartIdx,
674 : const size_t *count,
675 : const GInt64 *arrayStep,
676 : const GPtrDiff_t *bufferStride,
677 : const GDALExtendedDataType &bufferDataType,
678 : void *pDstBuffer) const
679 : {
680 828 : std::vector<std::unique_ptr<GDALEDTComponent>> comps;
681 828 : comps.emplace_back(std::unique_ptr<GDALEDTComponent>(
682 828 : new GDALEDTComponent(m_srcCompName, 0, bufferDataType)));
683 : auto tmpDT(GDALExtendedDataType::Create(
684 828 : std::string(), bufferDataType.GetSize(), std::move(comps)));
685 :
686 414 : return m_poParent->Read(arrayStartIdx, count, arrayStep, bufferStride,
687 828 : tmpDT, pDstBuffer);
688 : }
689 :
690 : /************************************************************************/
691 : /* CreateFieldNameExtractArray() */
692 : /************************************************************************/
693 :
694 : static std::shared_ptr<GDALMDArray>
695 466 : CreateFieldNameExtractArray(const std::shared_ptr<GDALMDArray> &self,
696 : const std::string &fieldName)
697 : {
698 466 : CPLAssert(self->GetDataType().GetClass() == GEDTC_COMPOUND);
699 466 : const std::unique_ptr<GDALEDTComponent> *srcComp = nullptr;
700 983 : for (const auto &comp : self->GetDataType().GetComponents())
701 : {
702 982 : if (comp->GetName() == fieldName)
703 : {
704 465 : srcComp = ∁
705 465 : break;
706 : }
707 : }
708 466 : if (srcComp == nullptr)
709 : {
710 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find field %s",
711 : fieldName.c_str());
712 1 : return nullptr;
713 : }
714 465 : return GDALExtractFieldMDArray::Create(self, fieldName, *srcComp);
715 : }
716 :
717 : //! @endcond
718 :
719 : /************************************************************************/
720 : /* GetView() */
721 : /************************************************************************/
722 :
723 : // clang-format off
724 : /** Return a view of the array using slicing or field access.
725 : *
726 : * The slice expression uses the same syntax as NumPy basic slicing and
727 : * indexing. See
728 : * https://www.numpy.org/devdocs/reference/arrays.indexing.html#basic-slicing-and-indexing
729 : * Or it can use field access by name. See
730 : * https://www.numpy.org/devdocs/reference/arrays.indexing.html#field-access
731 : *
732 : * Multiple [] bracket elements can be concatenated, with a slice expression
733 : * or field name inside each.
734 : *
735 : * For basic slicing and indexing, inside each [] bracket element, a list of
736 : * indexes that apply to successive source dimensions, can be specified, using
737 : * integer indexing (e.g. 1), range indexing (start:stop:step), ellipsis (...)
738 : * or newaxis, using a comma separator.
739 : *
740 : * Examples with a 2-dimensional array whose content is [[0,1,2,3],[4,5,6,7]].
741 : * <ul>
742 : * <li>GetView("[1][2]"): returns a 0-dimensional/scalar array with the value
743 : * at index 1 in the first dimension, and index 2 in the second dimension
744 : * from the source array. That is 5</li>
745 : * <li>GetView("[1]")->GetView("[2]"): same as above. Above is actually
746 : * implemented internally doing this intermediate slicing approach.</li>
747 : * <li>GetView("[1,2]"): same as above, but a bit more performant.</li>
748 : * <li>GetView("[1]"): returns a 1-dimensional array, sliced at index 1 in the
749 : * first dimension. That is [4,5,6,7].</li>
750 : * <li>GetView("[:,2]"): returns a 1-dimensional array, sliced at index 2 in the
751 : * second dimension. That is [2,6].</li>
752 : * <li>GetView("[:,2:3:]"): returns a 2-dimensional array, sliced at index 2 in
753 : * the second dimension. That is [[2],[6]].</li>
754 : * <li>GetView("[::,2]"): Same as
755 : * above.</li> <li>GetView("[...,2]"): same as above, in that case, since the
756 : * ellipsis only expands to one dimension here.</li>
757 : * <li>GetView("[:,::2]"):
758 : * returns a 2-dimensional array, with even-indexed elements of the second
759 : * dimension. That is [[0,2],[4,6]].</li>
760 : * <li>GetView("[:,1::2]"): returns a
761 : * 2-dimensional array, with odd-indexed elements of the second dimension. That
762 : * is [[1,3],[5,7]].</li>
763 : * <li>GetView("[:,1:3:]"): returns a 2-dimensional
764 : * array, with elements of the second dimension with index in the range [1,3[.
765 : * That is [[1,2],[5,6]].</li>
766 : * <li>GetView("[::-1,:]"): returns a 2-dimensional
767 : * array, with the values in first dimension reversed. That is
768 : * [[4,5,6,7],[0,1,2,3]].</li>
769 : * <li>GetView("[newaxis,...]"): returns a
770 : * 3-dimensional array, with an additional dimension of size 1 put at the
771 : * beginning. That is [[[0,1,2,3],[4,5,6,7]]].</li>
772 : * </ul>
773 : *
774 : * One difference with NumPy behavior is that ranges that would result in
775 : * zero elements are not allowed (dimensions of size 0 not being allowed in the
776 : * GDAL multidimensional model).
777 : *
778 : * For field access, the syntax to use is ["field_name"] or ['field_name'].
779 : * Multiple field specification is not supported currently.
780 : *
781 : * Both type of access can be combined, e.g. GetView("[1]['field_name']")
782 : *
783 : * \note When using the GDAL Python bindings, natural Python syntax can be
784 : * used. That is ar[0,::,1]["foo"] will be internally translated to
785 : * ar.GetView("[0,::,1]['foo']")
786 : * \note When using the C++ API and integer indexing only, you may use the
787 : * at(idx0, idx1, ...) method.
788 : *
789 : * The returned array holds a reference to the original one, and thus is
790 : * a view of it (not a copy). If the content of the original array changes,
791 : * the content of the view array too. When using basic slicing and indexing,
792 : * the view can be written if the underlying array is writable.
793 : *
794 : * This is the same as the C function GDALMDArrayGetView()
795 : *
796 : * @param viewExpr Expression expressing basic slicing and indexing, or field
797 : * access.
798 : * @return a new array, that holds a reference to the original one, and thus is
799 : * a view of it (not a copy), or nullptr in case of error.
800 : */
801 : // clang-format on
802 :
803 : std::shared_ptr<GDALMDArray>
804 1124 : GDALMDArray::GetView(const std::string &viewExpr) const
805 : {
806 2248 : std::vector<ViewSpec> viewSpecs;
807 2248 : return GetView(viewExpr, true, viewSpecs);
808 : }
809 :
810 : //! @cond Doxygen_Suppress
811 : std::shared_ptr<GDALMDArray>
812 1196 : GDALMDArray::GetView(const std::string &viewExpr, bool bRenameDimensions,
813 : std::vector<ViewSpec> &viewSpecs) const
814 : {
815 2392 : auto self = std::dynamic_pointer_cast<GDALMDArray>(m_pSelf.lock());
816 1196 : if (!self)
817 : {
818 1 : CPLError(CE_Failure, CPLE_AppDefined,
819 : "Driver implementation issue: m_pSelf not set !");
820 1 : return nullptr;
821 : }
822 1195 : std::string curExpr(viewExpr);
823 : while (true)
824 : {
825 1198 : if (curExpr.empty() || curExpr[0] != '[')
826 : {
827 2 : CPLError(CE_Failure, CPLE_AppDefined,
828 : "Slice string should start with ['");
829 1195 : return nullptr;
830 : }
831 :
832 1196 : std::string fieldName;
833 : size_t endExpr;
834 1196 : if (curExpr.size() > 2 && (curExpr[1] == '"' || curExpr[1] == '\''))
835 : {
836 470 : if (self->GetDataType().GetClass() != GEDTC_COMPOUND)
837 : {
838 2 : CPLError(CE_Failure, CPLE_AppDefined,
839 : "Field access not allowed on non-compound data type");
840 2 : return nullptr;
841 : }
842 468 : size_t idx = 2;
843 5119 : for (; idx < curExpr.size(); idx++)
844 : {
845 5118 : const char ch = curExpr[idx];
846 5118 : if (ch == curExpr[1])
847 467 : break;
848 4651 : if (ch == '\\' && idx + 1 < curExpr.size())
849 : {
850 1 : fieldName += curExpr[idx + 1];
851 1 : idx++;
852 : }
853 : else
854 : {
855 4650 : fieldName += ch;
856 : }
857 : }
858 468 : if (idx + 1 >= curExpr.size() || curExpr[idx + 1] != ']')
859 : {
860 2 : CPLError(CE_Failure, CPLE_AppDefined,
861 : "Invalid field access specification");
862 2 : return nullptr;
863 : }
864 466 : endExpr = idx + 1;
865 : }
866 : else
867 : {
868 726 : endExpr = curExpr.find(']');
869 : }
870 1192 : if (endExpr == std::string::npos)
871 : {
872 1 : CPLError(CE_Failure, CPLE_AppDefined, "Missing ]'");
873 1 : return nullptr;
874 : }
875 1191 : if (endExpr == 1)
876 : {
877 1 : CPLError(CE_Failure, CPLE_AppDefined, "[] not allowed");
878 1 : return nullptr;
879 : }
880 1190 : std::string activeSlice(curExpr.substr(1, endExpr - 1));
881 :
882 1190 : if (!fieldName.empty())
883 : {
884 932 : ViewSpec viewSpec;
885 466 : viewSpec.m_osFieldName = fieldName;
886 466 : viewSpecs.emplace_back(std::move(viewSpec));
887 : }
888 :
889 1190 : auto newArray = !fieldName.empty()
890 : ? CreateFieldNameExtractArray(self, fieldName)
891 : : CreateSlicedArray(self, viewExpr, activeSlice,
892 1190 : bRenameDimensions, viewSpecs);
893 :
894 1190 : if (endExpr == curExpr.size() - 1)
895 : {
896 1187 : return newArray;
897 : }
898 3 : self = std::move(newArray);
899 3 : curExpr = curExpr.substr(endExpr + 1);
900 3 : }
901 : }
902 :
903 : //! @endcond
904 :
905 : std::shared_ptr<GDALMDArray>
906 19 : GDALMDArray::GetView(const std::vector<GUInt64> &indices) const
907 : {
908 19 : std::string osExpr("[");
909 19 : bool bFirst = true;
910 45 : for (const auto &idx : indices)
911 : {
912 26 : if (!bFirst)
913 7 : osExpr += ',';
914 26 : bFirst = false;
915 26 : osExpr += CPLSPrintf(CPL_FRMT_GUIB, static_cast<GUIntBig>(idx));
916 : }
917 57 : return GetView(osExpr + ']');
918 : }
919 :
920 : /************************************************************************/
921 : /* operator[] */
922 : /************************************************************************/
923 :
924 : /** Return a view of the array using field access
925 : *
926 : * Equivalent of GetView("['fieldName']")
927 : *
928 : * \note When operating on a shared_ptr, use (*array)["fieldName"] syntax.
929 : */
930 : std::shared_ptr<GDALMDArray>
931 2 : GDALMDArray::operator[](const std::string &fieldName) const
932 : {
933 2 : return GetView(CPLSPrintf("['%s']", CPLString(fieldName)
934 4 : .replaceAll('\\', "\\\\")
935 4 : .replaceAll('\'', "\\\'")
936 6 : .c_str()));
937 : }
|
