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