Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: gdalmultidim_array.cpp
4 : * Project: GDAL Core
5 : * Purpose: Implementation of core GDALMDArray functionality
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 "cpl_float.h"
15 : #include "gdal_multidim.h"
16 : #include "gdalmultidim_array_unscaled.h"
17 : #include "gdal_pam.h"
18 : #include "memmultidim.h"
19 :
20 : #include <algorithm>
21 : #include <cmath>
22 : #include <ctype.h> // isalnum
23 : #include <limits>
24 :
25 : #if defined(__clang__) || defined(_MSC_VER)
26 : #define COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT
27 : #endif
28 :
29 : /************************************************************************/
30 : /* SetUnit() */
31 : /************************************************************************/
32 :
33 : /** Set the variable unit.
34 : *
35 : * Values should conform as much as possible with those allowed by
36 : * the NetCDF CF conventions:
37 : * http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#units
38 : * but others might be returned.
39 : *
40 : * Few examples are "meter", "degrees", "second", ...
41 : * Empty value means unknown.
42 : *
43 : * This is the same as the C function GDALMDArraySetUnit()
44 : *
45 : * @note Driver implementation: optionally implemented.
46 : *
47 : * @param osUnit unit name.
48 : * @return true in case of success.
49 : */
50 0 : bool GDALMDArray::SetUnit(CPL_UNUSED const std::string &osUnit)
51 : {
52 0 : CPLError(CE_Failure, CPLE_NotSupported, "SetUnit() not implemented");
53 0 : return false;
54 : }
55 :
56 : /************************************************************************/
57 : /* GetUnit() */
58 : /************************************************************************/
59 :
60 : /** Return the array unit.
61 : *
62 : * Values should conform as much as possible with those allowed by
63 : * the NetCDF CF conventions:
64 : * http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#units
65 : * but others might be returned.
66 : *
67 : * Few examples are "meter", "degrees", "second", ...
68 : * Empty value means unknown.
69 : *
70 : * This is the same as the C function GDALMDArrayGetUnit()
71 : */
72 9 : const std::string &GDALMDArray::GetUnit() const
73 : {
74 9 : static const std::string emptyString;
75 9 : return emptyString;
76 : }
77 :
78 : /************************************************************************/
79 : /* SetSpatialRef() */
80 : /************************************************************************/
81 :
82 : /** Assign a spatial reference system object to the array.
83 : *
84 : * This is the same as the C function GDALMDArraySetSpatialRef().
85 : */
86 0 : bool GDALMDArray::SetSpatialRef(CPL_UNUSED const OGRSpatialReference *poSRS)
87 : {
88 0 : CPLError(CE_Failure, CPLE_NotSupported, "SetSpatialRef() not implemented");
89 0 : return false;
90 : }
91 :
92 : /************************************************************************/
93 : /* GetSpatialRef() */
94 : /************************************************************************/
95 :
96 : /** Return the spatial reference system object associated with the array.
97 : *
98 : * This is the same as the C function GDALMDArrayGetSpatialRef().
99 : */
100 8 : std::shared_ptr<OGRSpatialReference> GDALMDArray::GetSpatialRef() const
101 : {
102 8 : return nullptr;
103 : }
104 :
105 : /************************************************************************/
106 : /* GetRawNoDataValue() */
107 : /************************************************************************/
108 :
109 : /** Return the nodata value as a "raw" value.
110 : *
111 : * The value returned might be nullptr in case of no nodata value. When
112 : * a nodata value is registered, a non-nullptr will be returned whose size in
113 : * bytes is GetDataType().GetSize().
114 : *
115 : * The returned value should not be modified or freed. It is valid until
116 : * the array is destroyed, or the next call to GetRawNoDataValue() or
117 : * SetRawNoDataValue(), or any similar methods.
118 : *
119 : * @note Driver implementation: this method shall be implemented if nodata
120 : * is supported.
121 : *
122 : * This is the same as the C function GDALMDArrayGetRawNoDataValue().
123 : *
124 : * @return nullptr or a pointer to GetDataType().GetSize() bytes.
125 : */
126 9 : const void *GDALMDArray::GetRawNoDataValue() const
127 : {
128 9 : return nullptr;
129 : }
130 :
131 : /************************************************************************/
132 : /* GetNoDataValueAsDouble() */
133 : /************************************************************************/
134 :
135 : /** Return the nodata value as a double.
136 : *
137 : * This is the same as the C function GDALMDArrayGetNoDataValueAsDouble().
138 : *
139 : * @param pbHasNoData Pointer to a output boolean that will be set to true if
140 : * a nodata value exists and can be converted to double. Might be nullptr.
141 : *
142 : * @return the nodata value as a double. A 0.0 value might also indicate the
143 : * absence of a nodata value or an error in the conversion (*pbHasNoData will be
144 : * set to false then).
145 : */
146 23219 : double GDALMDArray::GetNoDataValueAsDouble(bool *pbHasNoData) const
147 : {
148 23219 : const void *pNoData = GetRawNoDataValue();
149 23219 : double dfNoData = 0.0;
150 23219 : const auto &eDT = GetDataType();
151 23219 : const bool ok = pNoData != nullptr && eDT.GetClass() == GEDTC_NUMERIC;
152 23219 : if (ok)
153 : {
154 22876 : GDALCopyWords64(pNoData, eDT.GetNumericDataType(), 0, &dfNoData,
155 : GDT_Float64, 0, 1);
156 : }
157 23219 : if (pbHasNoData)
158 525 : *pbHasNoData = ok;
159 23219 : return dfNoData;
160 : }
161 :
162 : /************************************************************************/
163 : /* GetNoDataValueAsInt64() */
164 : /************************************************************************/
165 :
166 : /** Return the nodata value as a Int64.
167 : *
168 : * @param pbHasNoData Pointer to a output boolean that will be set to true if
169 : * a nodata value exists and can be converted to Int64. Might be nullptr.
170 : *
171 : * This is the same as the C function GDALMDArrayGetNoDataValueAsInt64().
172 : *
173 : * @return the nodata value as a Int64
174 : *
175 : * @since GDAL 3.5
176 : */
177 14 : int64_t GDALMDArray::GetNoDataValueAsInt64(bool *pbHasNoData) const
178 : {
179 14 : const void *pNoData = GetRawNoDataValue();
180 14 : int64_t nNoData = GDAL_PAM_DEFAULT_NODATA_VALUE_INT64;
181 14 : const auto &eDT = GetDataType();
182 14 : const bool ok = pNoData != nullptr && eDT.GetClass() == GEDTC_NUMERIC;
183 14 : if (ok)
184 : {
185 10 : GDALCopyWords64(pNoData, eDT.GetNumericDataType(), 0, &nNoData,
186 : GDT_Int64, 0, 1);
187 : }
188 14 : if (pbHasNoData)
189 14 : *pbHasNoData = ok;
190 14 : return nNoData;
191 : }
192 :
193 : /************************************************************************/
194 : /* GetNoDataValueAsUInt64() */
195 : /************************************************************************/
196 :
197 : /** Return the nodata value as a UInt64.
198 : *
199 : * This is the same as the C function GDALMDArrayGetNoDataValueAsUInt64().
200 :
201 : * @param pbHasNoData Pointer to a output boolean that will be set to true if
202 : * a nodata value exists and can be converted to UInt64. Might be nullptr.
203 : *
204 : * @return the nodata value as a UInt64
205 : *
206 : * @since GDAL 3.5
207 : */
208 8 : uint64_t GDALMDArray::GetNoDataValueAsUInt64(bool *pbHasNoData) const
209 : {
210 8 : const void *pNoData = GetRawNoDataValue();
211 8 : uint64_t nNoData = GDAL_PAM_DEFAULT_NODATA_VALUE_UINT64;
212 8 : const auto &eDT = GetDataType();
213 8 : const bool ok = pNoData != nullptr && eDT.GetClass() == GEDTC_NUMERIC;
214 8 : if (ok)
215 : {
216 6 : GDALCopyWords64(pNoData, eDT.GetNumericDataType(), 0, &nNoData,
217 : GDT_UInt64, 0, 1);
218 : }
219 8 : if (pbHasNoData)
220 8 : *pbHasNoData = ok;
221 8 : return nNoData;
222 : }
223 :
224 : /************************************************************************/
225 : /* SetRawNoDataValue() */
226 : /************************************************************************/
227 :
228 : /** Set the nodata value as a "raw" value.
229 : *
230 : * The value passed might be nullptr in case of no nodata value. When
231 : * a nodata value is registered, a non-nullptr whose size in
232 : * bytes is GetDataType().GetSize() must be passed.
233 : *
234 : * This is the same as the C function GDALMDArraySetRawNoDataValue().
235 : *
236 : * @note Driver implementation: this method shall be implemented if setting
237 : nodata
238 : * is supported.
239 :
240 : * @return true in case of success.
241 : */
242 0 : bool GDALMDArray::SetRawNoDataValue(CPL_UNUSED const void *pRawNoData)
243 : {
244 0 : CPLError(CE_Failure, CPLE_NotSupported,
245 : "SetRawNoDataValue() not implemented");
246 0 : return false;
247 : }
248 :
249 : /************************************************************************/
250 : /* SetNoDataValue() */
251 : /************************************************************************/
252 :
253 : /** Set the nodata value as a double.
254 : *
255 : * If the natural data type of the attribute/array is not double, type
256 : * conversion will occur to the type returned by GetDataType().
257 : *
258 : * This is the same as the C function GDALMDArraySetNoDataValueAsDouble().
259 : *
260 : * @return true in case of success.
261 : */
262 64 : bool GDALMDArray::SetNoDataValue(double dfNoData)
263 : {
264 64 : void *pRawNoData = CPLMalloc(GetDataType().GetSize());
265 64 : bool bRet = false;
266 64 : if (GDALExtendedDataType::CopyValue(
267 128 : &dfNoData, GDALExtendedDataType::Create(GDT_Float64), pRawNoData,
268 64 : GetDataType()))
269 : {
270 64 : bRet = SetRawNoDataValue(pRawNoData);
271 : }
272 64 : CPLFree(pRawNoData);
273 64 : return bRet;
274 : }
275 :
276 : /************************************************************************/
277 : /* SetNoDataValue() */
278 : /************************************************************************/
279 :
280 : /** Set the nodata value as a Int64.
281 : *
282 : * If the natural data type of the attribute/array is not Int64, type conversion
283 : * will occur to the type returned by GetDataType().
284 : *
285 : * This is the same as the C function GDALMDArraySetNoDataValueAsInt64().
286 : *
287 : * @return true in case of success.
288 : *
289 : * @since GDAL 3.5
290 : */
291 6 : bool GDALMDArray::SetNoDataValue(int64_t nNoData)
292 : {
293 6 : void *pRawNoData = CPLMalloc(GetDataType().GetSize());
294 6 : bool bRet = false;
295 6 : if (GDALExtendedDataType::CopyValue(&nNoData,
296 12 : GDALExtendedDataType::Create(GDT_Int64),
297 6 : pRawNoData, GetDataType()))
298 : {
299 6 : bRet = SetRawNoDataValue(pRawNoData);
300 : }
301 6 : CPLFree(pRawNoData);
302 6 : return bRet;
303 : }
304 :
305 : /************************************************************************/
306 : /* SetNoDataValue() */
307 : /************************************************************************/
308 :
309 : /** Set the nodata value as a Int64.
310 : *
311 : * If the natural data type of the attribute/array is not Int64, type conversion
312 : * will occur to the type returned by GetDataType().
313 : *
314 : * This is the same as the C function GDALMDArraySetNoDataValueAsUInt64().
315 : *
316 : * @return true in case of success.
317 : *
318 : * @since GDAL 3.5
319 : */
320 1 : bool GDALMDArray::SetNoDataValue(uint64_t nNoData)
321 : {
322 1 : void *pRawNoData = CPLMalloc(GetDataType().GetSize());
323 1 : bool bRet = false;
324 1 : if (GDALExtendedDataType::CopyValue(
325 2 : &nNoData, GDALExtendedDataType::Create(GDT_UInt64), pRawNoData,
326 1 : GetDataType()))
327 : {
328 1 : bRet = SetRawNoDataValue(pRawNoData);
329 : }
330 1 : CPLFree(pRawNoData);
331 1 : return bRet;
332 : }
333 :
334 : /************************************************************************/
335 : /* Resize() */
336 : /************************************************************************/
337 :
338 : /** Resize an array to new dimensions.
339 : *
340 : * Not all drivers may allow this operation, and with restrictions (e.g.
341 : * for netCDF, this is limited to growing of "unlimited" dimensions)
342 : *
343 : * Resizing a dimension used in other arrays will cause those other arrays
344 : * to be resized.
345 : *
346 : * This is the same as the C function GDALMDArrayResize().
347 : *
348 : * @param anNewDimSizes Array of GetDimensionCount() values containing the
349 : * new size of each indexing dimension.
350 : * @param papszOptions Options. (Driver specific)
351 : * @return true in case of success.
352 : * @since GDAL 3.7
353 : */
354 0 : bool GDALMDArray::Resize(CPL_UNUSED const std::vector<GUInt64> &anNewDimSizes,
355 : CPL_UNUSED CSLConstList papszOptions)
356 : {
357 0 : CPLError(CE_Failure, CPLE_NotSupported,
358 : "Resize() is not supported for this array");
359 0 : return false;
360 : }
361 :
362 : /************************************************************************/
363 : /* SetScale() */
364 : /************************************************************************/
365 :
366 : /** Set the scale value to apply to raw values.
367 : *
368 : * unscaled_value = raw_value * GetScale() + GetOffset()
369 : *
370 : * This is the same as the C function GDALMDArraySetScale() /
371 : * GDALMDArraySetScaleEx().
372 : *
373 : * @note Driver implementation: this method shall be implemented if setting
374 : * scale is supported.
375 : *
376 : * @param dfScale scale
377 : * @param eStorageType Data type to which create the potential attribute that
378 : * will store the scale. Added in GDAL 3.3 If let to its GDT_Unknown value, the
379 : * implementation will decide automatically the data type. Note that changing
380 : * the data type after initial setting might not be supported.
381 : * @return true in case of success.
382 : */
383 0 : bool GDALMDArray::SetScale(CPL_UNUSED double dfScale,
384 : CPL_UNUSED GDALDataType eStorageType)
385 : {
386 0 : CPLError(CE_Failure, CPLE_NotSupported, "SetScale() not implemented");
387 0 : return false;
388 : }
389 :
390 : /************************************************************************/
391 : /* SetOffset) */
392 : /************************************************************************/
393 :
394 : /** Set the offset value to apply to raw values.
395 : *
396 : * unscaled_value = raw_value * GetScale() + GetOffset()
397 : *
398 : * This is the same as the C function GDALMDArraySetOffset() /
399 : * GDALMDArraySetOffsetEx().
400 : *
401 : * @note Driver implementation: this method shall be implemented if setting
402 : * offset is supported.
403 : *
404 : * @param dfOffset Offset
405 : * @param eStorageType Data type to which create the potential attribute that
406 : * will store the offset. Added in GDAL 3.3 If let to its GDT_Unknown value, the
407 : * implementation will decide automatically the data type. Note that changing
408 : * the data type after initial setting might not be supported.
409 : * @return true in case of success.
410 : */
411 0 : bool GDALMDArray::SetOffset(CPL_UNUSED double dfOffset,
412 : CPL_UNUSED GDALDataType eStorageType)
413 : {
414 0 : CPLError(CE_Failure, CPLE_NotSupported, "SetOffset() not implemented");
415 0 : return false;
416 : }
417 :
418 : /************************************************************************/
419 : /* GetScale() */
420 : /************************************************************************/
421 :
422 : /** Get the scale value to apply to raw values.
423 : *
424 : * unscaled_value = raw_value * GetScale() + GetOffset()
425 : *
426 : * This is the same as the C function GDALMDArrayGetScale().
427 : *
428 : * @note Driver implementation: this method shall be implemented if getting
429 : * scale is supported.
430 : *
431 : * @param pbHasScale Pointer to a output boolean that will be set to true if
432 : * a scale value exists. Might be nullptr.
433 : * @param peStorageType Pointer to a output GDALDataType that will be set to
434 : * the storage type of the scale value, when known/relevant. Otherwise will be
435 : * set to GDT_Unknown. Might be nullptr. Since GDAL 3.3
436 : *
437 : * @return the scale value. A 1.0 value might also indicate the
438 : * absence of a scale value.
439 : */
440 26 : double GDALMDArray::GetScale(CPL_UNUSED bool *pbHasScale,
441 : CPL_UNUSED GDALDataType *peStorageType) const
442 : {
443 26 : if (pbHasScale)
444 26 : *pbHasScale = false;
445 26 : return 1.0;
446 : }
447 :
448 : /************************************************************************/
449 : /* GetOffset() */
450 : /************************************************************************/
451 :
452 : /** Get the offset value to apply to raw values.
453 : *
454 : * unscaled_value = raw_value * GetScale() + GetOffset()
455 : *
456 : * This is the same as the C function GDALMDArrayGetOffset().
457 : *
458 : * @note Driver implementation: this method shall be implemented if getting
459 : * offset is supported.
460 : *
461 : * @param pbHasOffset Pointer to a output boolean that will be set to true if
462 : * a offset value exists. Might be nullptr.
463 : * @param peStorageType Pointer to a output GDALDataType that will be set to
464 : * the storage type of the offset value, when known/relevant. Otherwise will be
465 : * set to GDT_Unknown. Might be nullptr. Since GDAL 3.3
466 : *
467 : * @return the offset value. A 0.0 value might also indicate the
468 : * absence of a offset value.
469 : */
470 26 : double GDALMDArray::GetOffset(CPL_UNUSED bool *pbHasOffset,
471 : CPL_UNUSED GDALDataType *peStorageType) const
472 : {
473 26 : if (pbHasOffset)
474 26 : *pbHasOffset = false;
475 26 : return 0.0;
476 : }
477 :
478 : /************************************************************************/
479 : /* GDALMDArray() */
480 : /************************************************************************/
481 :
482 : //! @cond Doxygen_Suppress
483 10176 : GDALMDArray::GDALMDArray(CPL_UNUSED const std::string &osParentName,
484 : CPL_UNUSED const std::string &osName,
485 0 : const std::string &osContext)
486 : :
487 : #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
488 : GDALAbstractMDArray(osParentName, osName),
489 : #endif
490 10176 : m_osContext(osContext)
491 : {
492 10176 : }
493 :
494 : //! @endcond
495 :
496 : /************************************************************************/
497 : /* GetTotalCopyCost() */
498 : /************************************************************************/
499 :
500 : /** Return a total "cost" to copy the array.
501 : *
502 : * Used as a parameter for CopyFrom()
503 : */
504 73 : GUInt64 GDALMDArray::GetTotalCopyCost() const
505 : {
506 146 : return COPY_COST + GetAttributes().size() * GDALAttribute::COPY_COST +
507 146 : GetTotalElementsCount() * GetDataType().GetSize();
508 : }
509 :
510 : /************************************************************************/
511 : /* CopyFromAllExceptValues() */
512 : /************************************************************************/
513 :
514 : //! @cond Doxygen_Suppress
515 :
516 215 : bool GDALMDArray::CopyFromAllExceptValues(const GDALMDArray *poSrcArray,
517 : bool bStrict, GUInt64 &nCurCost,
518 : const GUInt64 nTotalCost,
519 : GDALProgressFunc pfnProgress,
520 : void *pProgressData)
521 : {
522 : // Nodata setting must be one of the first things done for TileDB
523 215 : const void *pNoData = poSrcArray->GetRawNoDataValue();
524 215 : if (pNoData && poSrcArray->GetDataType() == GetDataType())
525 : {
526 13 : SetRawNoDataValue(pNoData);
527 : }
528 :
529 215 : const bool bThisIsUnscaledArray =
530 215 : dynamic_cast<GDALMDArrayUnscaled *>(this) != nullptr;
531 430 : auto attrs = poSrcArray->GetAttributes();
532 298 : for (const auto &attr : attrs)
533 : {
534 83 : const auto &osAttrName = attr->GetName();
535 83 : if (bThisIsUnscaledArray)
536 : {
537 6 : if (osAttrName == "missing_value" || osAttrName == "_FillValue" ||
538 7 : osAttrName == "valid_min" || osAttrName == "valid_max" ||
539 1 : osAttrName == "valid_range")
540 : {
541 1 : continue;
542 : }
543 : }
544 :
545 82 : auto dstAttr = CreateAttribute(osAttrName, attr->GetDimensionsSize(),
546 164 : attr->GetDataType());
547 82 : if (!dstAttr)
548 : {
549 0 : if (bStrict)
550 0 : return false;
551 0 : continue;
552 : }
553 82 : auto raw = attr->ReadAsRaw();
554 82 : if (!dstAttr->Write(raw.data(), raw.size()) && bStrict)
555 0 : return false;
556 : }
557 215 : if (!attrs.empty())
558 : {
559 47 : nCurCost += attrs.size() * GDALAttribute::COPY_COST;
560 81 : if (pfnProgress &&
561 34 : !pfnProgress(double(nCurCost) / nTotalCost, "", pProgressData))
562 0 : return false;
563 : }
564 :
565 215 : auto srcSRS = poSrcArray->GetSpatialRef();
566 215 : if (srcSRS)
567 : {
568 22 : SetSpatialRef(srcSRS.get());
569 : }
570 :
571 215 : const std::string &osUnit(poSrcArray->GetUnit());
572 215 : if (!osUnit.empty())
573 : {
574 24 : SetUnit(osUnit);
575 : }
576 :
577 215 : bool bGotValue = false;
578 215 : GDALDataType eOffsetStorageType = GDT_Unknown;
579 : const double dfOffset =
580 215 : poSrcArray->GetOffset(&bGotValue, &eOffsetStorageType);
581 215 : if (bGotValue)
582 : {
583 3 : SetOffset(dfOffset, eOffsetStorageType);
584 : }
585 :
586 215 : bGotValue = false;
587 215 : GDALDataType eScaleStorageType = GDT_Unknown;
588 215 : const double dfScale = poSrcArray->GetScale(&bGotValue, &eScaleStorageType);
589 215 : if (bGotValue)
590 : {
591 3 : SetScale(dfScale, eScaleStorageType);
592 : }
593 :
594 215 : return true;
595 : }
596 :
597 : //! @endcond
598 :
599 : /************************************************************************/
600 : /* CopyFrom() */
601 : /************************************************************************/
602 :
603 : /** Copy the content of an array into a new (generally empty) array.
604 : *
605 : * @param poSrcDS Source dataset. Might be nullptr (but for correct behavior
606 : * of some output drivers this is not recommended)
607 : * @param poSrcArray Source array. Should NOT be nullptr.
608 : * @param bStrict Whether to enable strict mode. In strict mode, any error will
609 : * stop the copy. In relaxed mode, the copy will be attempted to
610 : * be pursued.
611 : * @param nCurCost Should be provided as a variable initially set to 0.
612 : * @param nTotalCost Total cost from GetTotalCopyCost().
613 : * @param pfnProgress Progress callback, or nullptr.
614 : * @param pProgressData Progress user data, or nullptr.
615 : *
616 : * @return true in case of success (or partial success if bStrict == false).
617 : */
618 68 : bool GDALMDArray::CopyFrom(CPL_UNUSED GDALDataset *poSrcDS,
619 : const GDALMDArray *poSrcArray, bool bStrict,
620 : GUInt64 &nCurCost, const GUInt64 nTotalCost,
621 : GDALProgressFunc pfnProgress, void *pProgressData)
622 : {
623 68 : if (pfnProgress == nullptr)
624 4 : pfnProgress = GDALDummyProgress;
625 :
626 68 : nCurCost += GDALMDArray::COPY_COST;
627 :
628 68 : if (!CopyFromAllExceptValues(poSrcArray, bStrict, nCurCost, nTotalCost,
629 : pfnProgress, pProgressData))
630 : {
631 0 : return false;
632 : }
633 :
634 68 : const auto &dims = poSrcArray->GetDimensions();
635 68 : const auto nDTSize = poSrcArray->GetDataType().GetSize();
636 68 : if (dims.empty())
637 : {
638 2 : std::vector<GByte> abyTmp(nDTSize);
639 2 : if (!(poSrcArray->Read(nullptr, nullptr, nullptr, nullptr,
640 2 : GetDataType(), &abyTmp[0]) &&
641 2 : Write(nullptr, nullptr, nullptr, nullptr, GetDataType(),
642 4 : &abyTmp[0])) &&
643 : bStrict)
644 : {
645 0 : return false;
646 : }
647 2 : nCurCost += GetTotalElementsCount() * GetDataType().GetSize();
648 2 : if (!pfnProgress(double(nCurCost) / nTotalCost, "", pProgressData))
649 0 : return false;
650 : }
651 : else
652 : {
653 66 : std::vector<GUInt64> arrayStartIdx(dims.size());
654 66 : std::vector<GUInt64> count(dims.size());
655 172 : for (size_t i = 0; i < dims.size(); i++)
656 : {
657 106 : count[i] = static_cast<size_t>(dims[i]->GetSize());
658 : }
659 :
660 : struct CopyFunc
661 : {
662 : GDALMDArray *poDstArray = nullptr;
663 : std::vector<GByte> abyTmp{};
664 : GDALProgressFunc pfnProgress = nullptr;
665 : void *pProgressData = nullptr;
666 : GUInt64 nCurCost = 0;
667 : GUInt64 nTotalCost = 0;
668 : GUInt64 nTotalBytesThisArray = 0;
669 : bool bStop = false;
670 :
671 84 : static bool f(GDALAbstractMDArray *l_poSrcArray,
672 : const GUInt64 *chunkArrayStartIdx,
673 : const size_t *chunkCount, GUInt64 iCurChunk,
674 : GUInt64 nChunkCount, void *pUserData)
675 : {
676 84 : const auto &dt(l_poSrcArray->GetDataType());
677 84 : auto data = static_cast<CopyFunc *>(pUserData);
678 84 : auto poDstArray = data->poDstArray;
679 84 : if (!l_poSrcArray->Read(chunkArrayStartIdx, chunkCount, nullptr,
680 84 : nullptr, dt, &data->abyTmp[0]))
681 : {
682 1 : return false;
683 : }
684 : bool bRet =
685 83 : poDstArray->Write(chunkArrayStartIdx, chunkCount, nullptr,
686 83 : nullptr, dt, &data->abyTmp[0]);
687 83 : if (dt.NeedsFreeDynamicMemory())
688 : {
689 5 : const auto l_nDTSize = dt.GetSize();
690 5 : GByte *ptr = &data->abyTmp[0];
691 5 : const size_t l_nDims(l_poSrcArray->GetDimensionCount());
692 5 : size_t nEltCount = 1;
693 10 : for (size_t i = 0; i < l_nDims; ++i)
694 : {
695 5 : nEltCount *= chunkCount[i];
696 : }
697 22 : for (size_t i = 0; i < nEltCount; i++)
698 : {
699 17 : dt.FreeDynamicMemory(ptr);
700 17 : ptr += l_nDTSize;
701 : }
702 : }
703 83 : if (!bRet)
704 : {
705 0 : return false;
706 : }
707 :
708 83 : double dfCurCost =
709 83 : double(data->nCurCost) + double(iCurChunk) / nChunkCount *
710 83 : data->nTotalBytesThisArray;
711 83 : if (!data->pfnProgress(dfCurCost / data->nTotalCost, "",
712 : data->pProgressData))
713 : {
714 0 : data->bStop = true;
715 0 : return false;
716 : }
717 :
718 83 : return true;
719 : }
720 : };
721 :
722 66 : CopyFunc copyFunc;
723 66 : copyFunc.poDstArray = this;
724 66 : copyFunc.nCurCost = nCurCost;
725 66 : copyFunc.nTotalCost = nTotalCost;
726 66 : copyFunc.nTotalBytesThisArray = GetTotalElementsCount() * nDTSize;
727 66 : copyFunc.pfnProgress = pfnProgress;
728 66 : copyFunc.pProgressData = pProgressData;
729 : const char *pszSwathSize =
730 66 : CPLGetConfigOption("GDAL_SWATH_SIZE", nullptr);
731 : const size_t nMaxChunkSize =
732 : pszSwathSize
733 66 : ? static_cast<size_t>(
734 1 : std::min(GIntBig(std::numeric_limits<size_t>::max() / 2),
735 1 : CPLAtoGIntBig(pszSwathSize)))
736 : : static_cast<size_t>(
737 65 : std::min(GIntBig(std::numeric_limits<size_t>::max() / 2),
738 65 : GDALGetCacheMax64() / 4));
739 66 : const auto anChunkSizes(GetProcessingChunkSize(nMaxChunkSize));
740 66 : size_t nRealChunkSize = nDTSize;
741 172 : for (const auto &nChunkSize : anChunkSizes)
742 : {
743 106 : nRealChunkSize *= nChunkSize;
744 : }
745 : try
746 : {
747 66 : copyFunc.abyTmp.resize(nRealChunkSize);
748 : }
749 0 : catch (const std::exception &)
750 : {
751 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
752 : "Cannot allocate temporary buffer");
753 0 : nCurCost += copyFunc.nTotalBytesThisArray;
754 0 : return false;
755 : }
756 197 : if (copyFunc.nTotalBytesThisArray != 0 &&
757 65 : !const_cast<GDALMDArray *>(poSrcArray)
758 65 : ->ProcessPerChunk(arrayStartIdx.data(), count.data(),
759 : anChunkSizes.data(), CopyFunc::f,
760 132 : ©Func) &&
761 1 : (bStrict || copyFunc.bStop))
762 : {
763 0 : nCurCost += copyFunc.nTotalBytesThisArray;
764 0 : return false;
765 : }
766 66 : nCurCost += copyFunc.nTotalBytesThisArray;
767 : }
768 :
769 68 : return true;
770 : }
771 :
772 : /************************************************************************/
773 : /* GetStructuralInfo() */
774 : /************************************************************************/
775 :
776 : /** Return structural information on the array.
777 : *
778 : * This may be the compression, etc..
779 : *
780 : * The return value should not be freed and is valid until GDALMDArray is
781 : * released or this function called again.
782 : *
783 : * This is the same as the C function GDALMDArrayGetStructuralInfo().
784 : */
785 99 : CSLConstList GDALMDArray::GetStructuralInfo() const
786 : {
787 99 : return nullptr;
788 : }
789 :
790 : /************************************************************************/
791 : /* AdviseRead() */
792 : /************************************************************************/
793 :
794 : /** Advise driver of upcoming read requests.
795 : *
796 : * Some GDAL drivers operate more efficiently if they know in advance what
797 : * set of upcoming read requests will be made. The AdviseRead() method allows
798 : * an application to notify the driver of the region of interest.
799 : *
800 : * Many drivers just ignore the AdviseRead() call, but it can dramatically
801 : * accelerate access via some drivers. One such case is when reading through
802 : * a DAP dataset with the netCDF driver (a in-memory cache array is then created
803 : * with the region of interest defined by AdviseRead())
804 : *
805 : * This is the same as the C function GDALMDArrayAdviseRead().
806 : *
807 : * @param arrayStartIdx Values representing the starting index to read
808 : * in each dimension (in [0, aoDims[i].GetSize()-1] range).
809 : * Array of GetDimensionCount() values.
810 : * Can be nullptr as a synonymous for [0 for i in
811 : * range(GetDimensionCount() ]
812 : *
813 : * @param count Values representing the number of values to extract in
814 : * each dimension.
815 : * Array of GetDimensionCount() values.
816 : * Can be nullptr as a synonymous for
817 : * [ aoDims[i].GetSize() - arrayStartIdx[i] for i in
818 : * range(GetDimensionCount() ]
819 : *
820 : * @param papszOptions Driver specific options, or nullptr. Consult driver
821 : * documentation.
822 : *
823 : * @return true in case of success (ignoring the advice is a success)
824 : *
825 : * @since GDAL 3.2
826 : */
827 68 : bool GDALMDArray::AdviseRead(const GUInt64 *arrayStartIdx, const size_t *count,
828 : CSLConstList papszOptions) const
829 : {
830 68 : const auto nDimCount = GetDimensionCount();
831 68 : if (nDimCount == 0)
832 2 : return true;
833 :
834 132 : std::vector<GUInt64> tmp_arrayStartIdx;
835 66 : if (arrayStartIdx == nullptr)
836 : {
837 0 : tmp_arrayStartIdx.resize(nDimCount);
838 0 : arrayStartIdx = tmp_arrayStartIdx.data();
839 : }
840 :
841 132 : std::vector<size_t> tmp_count;
842 66 : if (count == nullptr)
843 : {
844 0 : tmp_count.resize(nDimCount);
845 0 : const auto &dims = GetDimensions();
846 0 : for (size_t i = 0; i < nDimCount; i++)
847 : {
848 0 : const GUInt64 nSize = dims[i]->GetSize() - arrayStartIdx[i];
849 : #if SIZEOF_VOIDP < 8
850 : if (nSize != static_cast<size_t>(nSize))
851 : {
852 : CPLError(CE_Failure, CPLE_AppDefined, "Integer overflow");
853 : return false;
854 : }
855 : #endif
856 0 : tmp_count[i] = static_cast<size_t>(nSize);
857 : }
858 0 : count = tmp_count.data();
859 : }
860 :
861 132 : std::vector<GInt64> tmp_arrayStep;
862 132 : std::vector<GPtrDiff_t> tmp_bufferStride;
863 66 : const GInt64 *arrayStep = nullptr;
864 66 : const GPtrDiff_t *bufferStride = nullptr;
865 66 : if (!CheckReadWriteParams(arrayStartIdx, count, arrayStep, bufferStride,
866 132 : GDALExtendedDataType::Create(GDT_Unknown),
867 : nullptr, nullptr, 0, tmp_arrayStep,
868 : tmp_bufferStride))
869 : {
870 2 : return false;
871 : }
872 :
873 64 : return IAdviseRead(arrayStartIdx, count, papszOptions);
874 : }
875 :
876 : /************************************************************************/
877 : /* IAdviseRead() */
878 : /************************************************************************/
879 :
880 : //! @cond Doxygen_Suppress
881 19 : bool GDALMDArray::IAdviseRead(const GUInt64 *, const size_t *,
882 : CSLConstList /* papszOptions*/) const
883 : {
884 19 : return true;
885 : }
886 :
887 : //! @endcond
888 :
889 : /************************************************************************/
890 : /* MassageName() */
891 : /************************************************************************/
892 :
893 : //! @cond Doxygen_Suppress
894 49 : /*static*/ std::string GDALMDArray::MassageName(const std::string &inputName)
895 : {
896 49 : std::string ret;
897 717 : for (const char ch : inputName)
898 : {
899 668 : if (!isalnum(static_cast<unsigned char>(ch)))
900 156 : ret += '_';
901 : else
902 512 : ret += ch;
903 : }
904 49 : return ret;
905 : }
906 :
907 : //! @endcond
908 :
909 : /************************************************************************/
910 : /* GetCacheRootGroup() */
911 : /************************************************************************/
912 :
913 : //! @cond Doxygen_Suppress
914 : std::shared_ptr<GDALGroup>
915 3569 : GDALMDArray::GetCacheRootGroup(bool bCanCreate,
916 : std::string &osCacheFilenameOut) const
917 : {
918 3569 : const auto &osFilename = GetFilename();
919 3569 : if (osFilename.empty())
920 : {
921 1 : CPLError(CE_Failure, CPLE_AppDefined,
922 : "Cannot cache an array with an empty filename");
923 1 : return nullptr;
924 : }
925 :
926 3568 : osCacheFilenameOut = osFilename + ".gmac";
927 3568 : if (STARTS_WITH(osFilename.c_str(), "/vsicurl/http"))
928 : {
929 2 : const auto nPosQuestionMark = osFilename.find('?');
930 2 : if (nPosQuestionMark != std::string::npos)
931 : {
932 : osCacheFilenameOut =
933 0 : osFilename.substr(0, nPosQuestionMark)
934 0 : .append(".gmac")
935 0 : .append(osFilename.substr(nPosQuestionMark));
936 : }
937 : }
938 3568 : const char *pszProxy = PamGetProxy(osCacheFilenameOut.c_str());
939 3568 : if (pszProxy != nullptr)
940 7 : osCacheFilenameOut = pszProxy;
941 :
942 : // .gmac sidecars are local-only; skip stat for non-local filesystems.
943 7116 : if (!bCanCreate && pszProxy == nullptr &&
944 3548 : !VSIIsLocal(osCacheFilenameOut.c_str()))
945 : {
946 2 : return nullptr;
947 : }
948 :
949 3566 : std::unique_ptr<GDALDataset> poDS;
950 : VSIStatBufL sStat;
951 3566 : if (VSIStatL(osCacheFilenameOut.c_str(), &sStat) == 0)
952 : {
953 42 : poDS.reset(GDALDataset::Open(osCacheFilenameOut.c_str(),
954 : GDAL_OF_MULTIDIM_RASTER | GDAL_OF_UPDATE,
955 : nullptr, nullptr, nullptr));
956 : }
957 3566 : if (poDS)
958 : {
959 42 : CPLDebug("GDAL", "Opening cache %s", osCacheFilenameOut.c_str());
960 42 : return poDS->GetRootGroup();
961 : }
962 :
963 3524 : if (bCanCreate)
964 : {
965 7 : const char *pszDrvName = "netCDF";
966 7 : GDALDriver *poDrv = GetGDALDriverManager()->GetDriverByName(pszDrvName);
967 7 : if (poDrv == nullptr)
968 : {
969 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver %s",
970 : pszDrvName);
971 0 : return nullptr;
972 : }
973 : {
974 14 : CPLErrorHandlerPusher oHandlerPusher(CPLQuietErrorHandler);
975 14 : CPLErrorStateBackuper oErrorStateBackuper;
976 7 : poDS.reset(poDrv->CreateMultiDimensional(osCacheFilenameOut.c_str(),
977 : nullptr, nullptr));
978 : }
979 7 : if (!poDS)
980 : {
981 1 : pszProxy = PamAllocateProxy(osCacheFilenameOut.c_str());
982 1 : if (pszProxy)
983 : {
984 1 : osCacheFilenameOut = pszProxy;
985 1 : poDS.reset(poDrv->CreateMultiDimensional(
986 : osCacheFilenameOut.c_str(), nullptr, nullptr));
987 : }
988 : }
989 7 : if (poDS)
990 : {
991 7 : CPLDebug("GDAL", "Creating cache %s", osCacheFilenameOut.c_str());
992 7 : return poDS->GetRootGroup();
993 : }
994 : else
995 : {
996 0 : CPLError(CE_Failure, CPLE_AppDefined,
997 : "Cannot create %s. Set the GDAL_PAM_PROXY_DIR "
998 : "configuration option to write the cache in "
999 : "another directory",
1000 : osCacheFilenameOut.c_str());
1001 : }
1002 : }
1003 :
1004 3517 : return nullptr;
1005 : }
1006 :
1007 : //! @endcond
1008 :
1009 : /************************************************************************/
1010 : /* Cache() */
1011 : /************************************************************************/
1012 :
1013 : /** Cache the content of the array into an auxiliary filename.
1014 : *
1015 : * The main purpose of this method is to be able to cache views that are
1016 : * expensive to compute, such as transposed arrays.
1017 : *
1018 : * The array will be stored in a file whose name is the one of
1019 : * GetFilename(), with an extra .gmac extension (stands for GDAL
1020 : * Multidimensional Array Cache). The cache is a netCDF dataset.
1021 : *
1022 : * If the .gmac file cannot be written next to the dataset, the
1023 : * GDAL_PAM_PROXY_DIR will be used, if set, to write the cache file into that
1024 : * directory.
1025 : *
1026 : * The GDALMDArray::Read() method will automatically use the cache when it
1027 : * exists. There is no timestamp checks between the source array and the cached
1028 : * array. If the source arrays changes, the cache must be manually deleted.
1029 : *
1030 : * This is the same as the C function GDALMDArrayCache()
1031 : *
1032 : * @note Driver implementation: optionally implemented.
1033 : *
1034 : * @param papszOptions List of options, null terminated, or NULL. Currently
1035 : * the only option supported is BLOCKSIZE=bs0,bs1,...,bsN
1036 : * to specify the block size of the cached array.
1037 : * @return true in case of success.
1038 : */
1039 7 : bool GDALMDArray::Cache(CSLConstList papszOptions) const
1040 : {
1041 14 : std::string osCacheFilename;
1042 14 : auto poRG = GetCacheRootGroup(true, osCacheFilename);
1043 7 : if (!poRG)
1044 1 : return false;
1045 :
1046 12 : const std::string osCachedArrayName(MassageName(GetFullName()));
1047 6 : if (poRG->OpenMDArray(osCachedArrayName))
1048 : {
1049 2 : CPLError(CE_Failure, CPLE_NotSupported,
1050 : "An array with same name %s already exists in %s",
1051 : osCachedArrayName.c_str(), osCacheFilename.c_str());
1052 2 : return false;
1053 : }
1054 :
1055 8 : CPLStringList aosOptions;
1056 4 : aosOptions.SetNameValue("COMPRESS", "DEFLATE");
1057 4 : const auto &aoDims = GetDimensions();
1058 8 : std::vector<std::shared_ptr<GDALDimension>> aoNewDims;
1059 4 : if (!aoDims.empty())
1060 : {
1061 : std::string osBlockSize(
1062 4 : CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", ""));
1063 4 : if (osBlockSize.empty())
1064 : {
1065 6 : const auto anBlockSize = GetBlockSize();
1066 3 : int idxDim = 0;
1067 10 : for (auto nBlockSize : anBlockSize)
1068 : {
1069 7 : if (idxDim > 0)
1070 4 : osBlockSize += ',';
1071 7 : if (nBlockSize == 0)
1072 7 : nBlockSize = 256;
1073 7 : nBlockSize = std::min(nBlockSize, aoDims[idxDim]->GetSize());
1074 : osBlockSize +=
1075 7 : std::to_string(static_cast<uint64_t>(nBlockSize));
1076 7 : idxDim++;
1077 : }
1078 : }
1079 4 : aosOptions.SetNameValue("BLOCKSIZE", osBlockSize.c_str());
1080 :
1081 4 : int idxDim = 0;
1082 13 : for (const auto &poDim : aoDims)
1083 : {
1084 9 : auto poNewDim = poRG->CreateDimension(
1085 18 : osCachedArrayName + '_' + std::to_string(idxDim),
1086 18 : poDim->GetType(), poDim->GetDirection(), poDim->GetSize());
1087 9 : if (!poNewDim)
1088 0 : return false;
1089 9 : aoNewDims.emplace_back(poNewDim);
1090 9 : idxDim++;
1091 : }
1092 : }
1093 :
1094 4 : auto poCachedArray = poRG->CreateMDArray(osCachedArrayName, aoNewDims,
1095 8 : GetDataType(), aosOptions.List());
1096 4 : if (!poCachedArray)
1097 : {
1098 0 : CPLError(CE_Failure, CPLE_NotSupported, "Cannot create %s in %s",
1099 : osCachedArrayName.c_str(), osCacheFilename.c_str());
1100 0 : return false;
1101 : }
1102 :
1103 4 : GUInt64 nCost = 0;
1104 8 : return poCachedArray->CopyFrom(nullptr, this,
1105 : false, // strict
1106 4 : nCost, GetTotalCopyCost(), nullptr, nullptr);
1107 : }
1108 :
1109 : /************************************************************************/
1110 : /* Read() */
1111 : /************************************************************************/
1112 :
1113 6457 : bool GDALMDArray::Read(const GUInt64 *arrayStartIdx, const size_t *count,
1114 : const GInt64 *arrayStep, // step in elements
1115 : const GPtrDiff_t *bufferStride, // stride in elements
1116 : const GDALExtendedDataType &bufferDataType,
1117 : void *pDstBuffer, const void *pDstBufferAllocStart,
1118 : size_t nDstBufferAllocSize) const
1119 : {
1120 6457 : if (!m_bHasTriedCachedArray)
1121 : {
1122 3490 : m_bHasTriedCachedArray = true;
1123 3490 : if (IsCacheable())
1124 : {
1125 3490 : const auto &osFilename = GetFilename();
1126 6252 : if (!osFilename.empty() &&
1127 6252 : !EQUAL(CPLGetExtensionSafe(osFilename.c_str()).c_str(), "gmac"))
1128 : {
1129 5488 : std::string osCacheFilename;
1130 5488 : auto poRG = GetCacheRootGroup(false, osCacheFilename);
1131 2744 : if (poRG)
1132 : {
1133 : const std::string osCachedArrayName(
1134 44 : MassageName(GetFullName()));
1135 22 : m_poCachedArray = poRG->OpenMDArray(osCachedArrayName);
1136 22 : if (m_poCachedArray)
1137 : {
1138 6 : const auto &dims = GetDimensions();
1139 : const auto &cachedDims =
1140 6 : m_poCachedArray->GetDimensions();
1141 6 : const size_t nDims = dims.size();
1142 : bool ok =
1143 12 : m_poCachedArray->GetDataType() == GetDataType() &&
1144 6 : cachedDims.size() == nDims;
1145 19 : for (size_t i = 0; ok && i < nDims; ++i)
1146 : {
1147 13 : ok = dims[i]->GetSize() == cachedDims[i]->GetSize();
1148 : }
1149 6 : if (ok)
1150 : {
1151 6 : CPLDebug("GDAL", "Cached array for %s found in %s",
1152 : osCachedArrayName.c_str(),
1153 : osCacheFilename.c_str());
1154 : }
1155 : else
1156 : {
1157 0 : CPLError(CE_Warning, CPLE_AppDefined,
1158 : "Cached array %s in %s has incompatible "
1159 : "characteristics with current array.",
1160 : osCachedArrayName.c_str(),
1161 : osCacheFilename.c_str());
1162 0 : m_poCachedArray.reset();
1163 : }
1164 : }
1165 : }
1166 : }
1167 : }
1168 : }
1169 :
1170 6457 : const auto array = m_poCachedArray ? m_poCachedArray.get() : this;
1171 6457 : if (!array->GetDataType().CanConvertTo(bufferDataType))
1172 : {
1173 0 : CPLError(CE_Failure, CPLE_AppDefined,
1174 : "Array data type is not convertible to buffer data type");
1175 0 : return false;
1176 : }
1177 :
1178 12914 : std::vector<GInt64> tmp_arrayStep;
1179 12914 : std::vector<GPtrDiff_t> tmp_bufferStride;
1180 6457 : if (!array->CheckReadWriteParams(arrayStartIdx, count, arrayStep,
1181 : bufferStride, bufferDataType, pDstBuffer,
1182 : pDstBufferAllocStart, nDstBufferAllocSize,
1183 : tmp_arrayStep, tmp_bufferStride))
1184 : {
1185 9 : return false;
1186 : }
1187 :
1188 6448 : return array->IRead(arrayStartIdx, count, arrayStep, bufferStride,
1189 6448 : bufferDataType, pDstBuffer);
1190 : }
1191 :
1192 : /************************************************************************/
1193 : /* GetRootGroup() */
1194 : /************************************************************************/
1195 :
1196 : /** Return the root group to which this arrays belongs too.
1197 : *
1198 : * Note that arrays may be free standing and some drivers may not implement
1199 : * this method, hence nullptr may be returned.
1200 : *
1201 : * It is used internally by the GetResampled() method to detect if GLT
1202 : * orthorectification is available.
1203 : *
1204 : * @return the root group, or nullptr.
1205 : * @since GDAL 3.8
1206 : */
1207 0 : std::shared_ptr<GDALGroup> GDALMDArray::GetRootGroup() const
1208 : {
1209 0 : return nullptr;
1210 : }
1211 :
1212 : //! @cond Doxygen_Suppress
1213 :
1214 : /************************************************************************/
1215 : /* IsStepOneContiguousRowMajorOrderedSameDataType() */
1216 : /************************************************************************/
1217 :
1218 : // Returns true if at all following conditions are met:
1219 : // arrayStep[] == 1, bufferDataType == GetDataType() and bufferStride[]
1220 : // defines a row-major ordered contiguous buffer.
1221 78 : bool GDALMDArray::IsStepOneContiguousRowMajorOrderedSameDataType(
1222 : const size_t *count, const GInt64 *arrayStep,
1223 : const GPtrDiff_t *bufferStride,
1224 : const GDALExtendedDataType &bufferDataType) const
1225 : {
1226 78 : if (bufferDataType != GetDataType())
1227 5 : return false;
1228 73 : size_t nExpectedStride = 1;
1229 166 : for (size_t i = GetDimensionCount(); i > 0;)
1230 : {
1231 96 : --i;
1232 96 : if (arrayStep[i] != 1 || bufferStride[i] < 0 ||
1233 94 : static_cast<size_t>(bufferStride[i]) != nExpectedStride)
1234 : {
1235 3 : return false;
1236 : }
1237 93 : nExpectedStride *= count[i];
1238 : }
1239 70 : return true;
1240 : }
1241 :
1242 : /************************************************************************/
1243 : /* ReadUsingContiguousIRead() */
1244 : /************************************************************************/
1245 :
1246 : // Used for example by the TileDB driver when requesting it with
1247 : // arrayStep[] != 1, bufferDataType != GetDataType() or bufferStride[]
1248 : // not defining a row-major ordered contiguous buffer.
1249 : // Should only be called when at least one of the above conditions are true,
1250 : // which can be tested with IsStepOneContiguousRowMajorOrderedSameDataType()
1251 : // returning none.
1252 : // This method will call IRead() again with arrayStep[] == 1,
1253 : // bufferDataType == GetDataType() and bufferStride[] defining a row-major
1254 : // ordered contiguous buffer, on a temporary buffer. And it will rearrange the
1255 : // content of that temporary buffer onto pDstBuffer.
1256 7 : bool GDALMDArray::ReadUsingContiguousIRead(
1257 : const GUInt64 *arrayStartIdx, const size_t *count, const GInt64 *arrayStep,
1258 : const GPtrDiff_t *bufferStride, const GDALExtendedDataType &bufferDataType,
1259 : void *pDstBuffer) const
1260 : {
1261 7 : const size_t nDims(GetDimensionCount());
1262 14 : std::vector<GUInt64> anTmpStartIdx(nDims);
1263 14 : std::vector<size_t> anTmpCount(nDims);
1264 7 : const auto &oType = GetDataType();
1265 7 : size_t nMemArraySize = oType.GetSize();
1266 14 : std::vector<GPtrDiff_t> anTmpStride(nDims);
1267 7 : GPtrDiff_t nStride = 1;
1268 18 : for (size_t i = nDims; i > 0;)
1269 : {
1270 11 : --i;
1271 11 : if (arrayStep[i] > 0)
1272 9 : anTmpStartIdx[i] = arrayStartIdx[i];
1273 : else
1274 2 : anTmpStartIdx[i] =
1275 2 : arrayStartIdx[i] - (count[i] - 1) * (-arrayStep[i]);
1276 : const uint64_t nCount =
1277 11 : (count[i] - 1) * static_cast<uint64_t>(std::abs(arrayStep[i])) + 1;
1278 11 : if (nCount > std::numeric_limits<size_t>::max() / nMemArraySize)
1279 : {
1280 0 : CPLError(CE_Failure, CPLE_AppDefined,
1281 : "Read() failed due to too large memory requirement");
1282 0 : return false;
1283 : }
1284 11 : anTmpCount[i] = static_cast<size_t>(nCount);
1285 11 : nMemArraySize *= anTmpCount[i];
1286 11 : anTmpStride[i] = nStride;
1287 11 : nStride *= anTmpCount[i];
1288 : }
1289 : std::unique_ptr<void, decltype(&VSIFree)> pTmpBuffer(
1290 14 : VSI_MALLOC_VERBOSE(nMemArraySize), VSIFree);
1291 7 : if (!pTmpBuffer)
1292 0 : return false;
1293 7 : if (!IRead(anTmpStartIdx.data(), anTmpCount.data(),
1294 14 : std::vector<GInt64>(nDims, 1).data(), // steps
1295 7 : anTmpStride.data(), oType, pTmpBuffer.get()))
1296 : {
1297 0 : return false;
1298 : }
1299 14 : std::vector<std::shared_ptr<GDALDimension>> apoTmpDims(nDims);
1300 18 : for (size_t i = 0; i < nDims; ++i)
1301 : {
1302 11 : if (arrayStep[i] > 0)
1303 9 : anTmpStartIdx[i] = 0;
1304 : else
1305 2 : anTmpStartIdx[i] = anTmpCount[i] - 1;
1306 22 : apoTmpDims[i] = std::make_shared<GDALDimension>(
1307 22 : std::string(), std::string(), std::string(), std::string(),
1308 22 : anTmpCount[i]);
1309 : }
1310 : auto poMEMArray =
1311 14 : MEMMDArray::Create(std::string(), std::string(), apoTmpDims, oType);
1312 21 : return poMEMArray->Init(static_cast<GByte *>(pTmpBuffer.get())) &&
1313 7 : poMEMArray->Read(anTmpStartIdx.data(), count, arrayStep,
1314 7 : bufferStride, bufferDataType, pDstBuffer);
1315 : }
1316 :
1317 : //! @endcond
1318 :
1319 : /************************************************************************/
1320 : /* GuessGeoTransform() */
1321 : /************************************************************************/
1322 :
1323 : /** Returns whether 2 specified dimensions form a geotransform
1324 : *
1325 : * @param nDimX Index of the X axis.
1326 : * @param nDimY Index of the Y axis.
1327 : * @param bPixelIsPoint Whether the geotransform should be returned
1328 : * with the pixel-is-point (pixel-center) convention
1329 : * (bPixelIsPoint = true), or with the pixel-is-area
1330 : * (top left corner convention)
1331 : * (bPixelIsPoint = false)
1332 : * @param[out] gt Computed geotransform
1333 : * @return true if a geotransform could be computed.
1334 : */
1335 500 : bool GDALMDArray::GuessGeoTransform(size_t nDimX, size_t nDimY,
1336 : bool bPixelIsPoint,
1337 : GDALGeoTransform >) const
1338 : {
1339 500 : const auto &dims(GetDimensions());
1340 1000 : auto poVarX = dims[nDimX]->GetIndexingVariable();
1341 1000 : auto poVarY = dims[nDimY]->GetIndexingVariable();
1342 500 : double dfXStart = 0.0;
1343 500 : double dfXSpacing = 0.0;
1344 500 : double dfYStart = 0.0;
1345 500 : double dfYSpacing = 0.0;
1346 1118 : if (poVarX && poVarX->GetDimensionCount() == 1 &&
1347 618 : poVarX->GetDimensions()[0]->GetSize() == dims[nDimX]->GetSize() &&
1348 903 : poVarY && poVarY->GetDimensionCount() == 1 &&
1349 297 : poVarY->GetDimensions()[0]->GetSize() == dims[nDimY]->GetSize() &&
1350 1101 : poVarX->IsRegularlySpaced(dfXStart, dfXSpacing) &&
1351 292 : poVarY->IsRegularlySpaced(dfYStart, dfYSpacing))
1352 : {
1353 290 : gt.xorig = dfXStart - (bPixelIsPoint ? 0 : dfXSpacing / 2);
1354 290 : gt.xscale = dfXSpacing;
1355 290 : gt.xrot = 0;
1356 290 : gt.yorig = dfYStart - (bPixelIsPoint ? 0 : dfYSpacing / 2);
1357 290 : gt.yrot = 0;
1358 290 : gt.yscale = dfYSpacing;
1359 290 : return true;
1360 : }
1361 210 : return false;
1362 : }
1363 :
1364 : /** Returns whether 2 specified dimensions form a geotransform
1365 : *
1366 : * @param nDimX Index of the X axis.
1367 : * @param nDimY Index of the Y axis.
1368 : * @param bPixelIsPoint Whether the geotransform should be returned
1369 : * with the pixel-is-point (pixel-center) convention
1370 : * (bPixelIsPoint = true), or with the pixel-is-area
1371 : * (top left corner convention)
1372 : * (bPixelIsPoint = false)
1373 : * @param[out] adfGeoTransform Computed geotransform
1374 : * @return true if a geotransform could be computed.
1375 : */
1376 0 : bool GDALMDArray::GuessGeoTransform(size_t nDimX, size_t nDimY,
1377 : bool bPixelIsPoint,
1378 : double adfGeoTransform[6]) const
1379 : {
1380 0 : GDALGeoTransform *gt =
1381 : reinterpret_cast<GDALGeoTransform *>(adfGeoTransform);
1382 0 : return GuessGeoTransform(nDimX, nDimY, bPixelIsPoint, *gt);
1383 : }
1384 :
1385 : /************************************************************************/
1386 : /* GetStatistics() */
1387 : /************************************************************************/
1388 :
1389 : /**
1390 : * \brief Fetch statistics.
1391 : *
1392 : * Returns the minimum, maximum, mean and standard deviation of all
1393 : * pixel values in this array.
1394 : *
1395 : * If bForce is FALSE results will only be returned if it can be done
1396 : * quickly (i.e. without scanning the data). If bForce is FALSE and
1397 : * results cannot be returned efficiently, the method will return CE_Warning
1398 : * but no warning will have been issued. This is a non-standard use of
1399 : * the CE_Warning return value to indicate "nothing done".
1400 : *
1401 : * When cached statistics are not available, and bForce is TRUE,
1402 : * ComputeStatistics() is called.
1403 : *
1404 : * Note that file formats using PAM (Persistent Auxiliary Metadata) services
1405 : * will generally cache statistics in the .aux.xml file allowing fast fetch
1406 : * after the first request.
1407 : *
1408 : * Cached statistics can be cleared with GDALDataset::ClearStatistics().
1409 : *
1410 : * This method is the same as the C function GDALMDArrayGetStatistics().
1411 : *
1412 : * @param bApproxOK Currently ignored. In the future, should be set to true
1413 : * if statistics on the whole array are wished, or to false if a subset of it
1414 : * may be used.
1415 : *
1416 : * @param bForce If false statistics will only be returned if it can
1417 : * be done without rescanning the image.
1418 : *
1419 : * @param pdfMin Location into which to load image minimum (may be NULL).
1420 : *
1421 : * @param pdfMax Location into which to load image maximum (may be NULL).-
1422 : *
1423 : * @param pdfMean Location into which to load image mean (may be NULL).
1424 : *
1425 : * @param pdfStdDev Location into which to load image standard deviation
1426 : * (may be NULL).
1427 : *
1428 : * @param pnValidCount Number of samples whose value is different from the
1429 : * nodata value. (may be NULL)
1430 : *
1431 : * @param pfnProgress a function to call to report progress, or NULL.
1432 : *
1433 : * @param pProgressData application data to pass to the progress function.
1434 : *
1435 : * @return CE_None on success, CE_Warning if no values returned,
1436 : * CE_Failure if an error occurs.
1437 : *
1438 : * @since GDAL 3.2
1439 : */
1440 :
1441 10 : CPLErr GDALMDArray::GetStatistics(bool bApproxOK, bool bForce, double *pdfMin,
1442 : double *pdfMax, double *pdfMean,
1443 : double *pdfStdDev, GUInt64 *pnValidCount,
1444 : GDALProgressFunc pfnProgress,
1445 : void *pProgressData)
1446 : {
1447 10 : if (!bForce)
1448 1 : return CE_Warning;
1449 :
1450 18 : return ComputeStatistics(bApproxOK, pdfMin, pdfMax, pdfMean, pdfStdDev,
1451 9 : pnValidCount, pfnProgress, pProgressData, nullptr)
1452 9 : ? CE_None
1453 9 : : CE_Failure;
1454 : }
1455 :
1456 : /************************************************************************/
1457 : /* ComputeStatistics() */
1458 : /************************************************************************/
1459 :
1460 : /**
1461 : * \brief Compute statistics.
1462 : *
1463 : * Returns the minimum, maximum, mean and standard deviation of all
1464 : * pixel values in this array.
1465 : *
1466 : * Pixels taken into account in statistics are those whose mask value
1467 : * (as determined by GetMask()) is non-zero.
1468 : *
1469 : * Once computed, the statistics will generally be "set" back on the
1470 : * owing dataset.
1471 : *
1472 : * Cached statistics can be cleared with GDALDataset::ClearStatistics().
1473 : *
1474 : * This method is the same as the C functions GDALMDArrayComputeStatistics().
1475 : * and GDALMDArrayComputeStatisticsEx().
1476 : *
1477 : * @param bApproxOK Currently ignored. In the future, should be set to true
1478 : * if statistics on the whole array are wished, or to false if a subset of it
1479 : * may be used.
1480 : *
1481 : * @param pdfMin Location into which to load image minimum (may be NULL).
1482 : *
1483 : * @param pdfMax Location into which to load image maximum (may be NULL).-
1484 : *
1485 : * @param pdfMean Location into which to load image mean (may be NULL).
1486 : *
1487 : * @param pdfStdDev Location into which to load image standard deviation
1488 : * (may be NULL).
1489 : *
1490 : * @param pnValidCount Number of samples whose value is different from the
1491 : * nodata value. (may be NULL)
1492 : *
1493 : * @param pfnProgress a function to call to report progress, or NULL.
1494 : *
1495 : * @param pProgressData application data to pass to the progress function.
1496 : *
1497 : * @param papszOptions NULL-terminated list of options, of NULL. Added in 3.8.
1498 : * Options are driver specific. For now the netCDF and Zarr
1499 : * drivers recognize UPDATE_METADATA=YES, whose effect is
1500 : * to add or update the actual_range attribute with the
1501 : * computed min/max, only if done on the full array, in non
1502 : * approximate mode, and the dataset is opened in update
1503 : * mode.
1504 : *
1505 : * @return true on success
1506 : *
1507 : * @since GDAL 3.2
1508 : */
1509 :
1510 13 : bool GDALMDArray::ComputeStatistics(bool bApproxOK, double *pdfMin,
1511 : double *pdfMax, double *pdfMean,
1512 : double *pdfStdDev, GUInt64 *pnValidCount,
1513 : GDALProgressFunc pfnProgress,
1514 : void *pProgressData,
1515 : CSLConstList papszOptions)
1516 : {
1517 : struct StatsPerChunkType
1518 : {
1519 : const GDALMDArray *array = nullptr;
1520 : std::shared_ptr<GDALMDArray> poMask{};
1521 : double dfMin = cpl::NumericLimits<double>::max();
1522 : double dfMax = -cpl::NumericLimits<double>::max();
1523 : double dfMean = 0.0;
1524 : double dfM2 = 0.0;
1525 : GUInt64 nValidCount = 0;
1526 : std::vector<GByte> abyData{};
1527 : std::vector<double> adfData{};
1528 : std::vector<GByte> abyMaskData{};
1529 : GDALProgressFunc pfnProgress = nullptr;
1530 : void *pProgressData = nullptr;
1531 : };
1532 :
1533 13 : const auto PerChunkFunc = [](GDALAbstractMDArray *,
1534 : const GUInt64 *chunkArrayStartIdx,
1535 : const size_t *chunkCount, GUInt64 iCurChunk,
1536 : GUInt64 nChunkCount, void *pUserData)
1537 : {
1538 13 : StatsPerChunkType *data = static_cast<StatsPerChunkType *>(pUserData);
1539 13 : const GDALMDArray *array = data->array;
1540 13 : const GDALMDArray *poMask = data->poMask.get();
1541 13 : const size_t nDims = array->GetDimensionCount();
1542 13 : size_t nVals = 1;
1543 34 : for (size_t i = 0; i < nDims; i++)
1544 21 : nVals *= chunkCount[i];
1545 :
1546 : // Get mask
1547 13 : data->abyMaskData.resize(nVals);
1548 13 : if (!(poMask->Read(chunkArrayStartIdx, chunkCount, nullptr, nullptr,
1549 13 : poMask->GetDataType(), &data->abyMaskData[0])))
1550 : {
1551 0 : return false;
1552 : }
1553 :
1554 : // Get data
1555 13 : const auto &oType = array->GetDataType();
1556 13 : if (oType.GetNumericDataType() == GDT_Float64)
1557 : {
1558 6 : data->adfData.resize(nVals);
1559 6 : if (!array->Read(chunkArrayStartIdx, chunkCount, nullptr, nullptr,
1560 6 : oType, &data->adfData[0]))
1561 : {
1562 0 : return false;
1563 : }
1564 : }
1565 : else
1566 : {
1567 7 : data->abyData.resize(nVals * oType.GetSize());
1568 7 : if (!array->Read(chunkArrayStartIdx, chunkCount, nullptr, nullptr,
1569 7 : oType, &data->abyData[0]))
1570 : {
1571 0 : return false;
1572 : }
1573 7 : data->adfData.resize(nVals);
1574 7 : GDALCopyWords64(&data->abyData[0], oType.GetNumericDataType(),
1575 7 : static_cast<int>(oType.GetSize()),
1576 7 : &data->adfData[0], GDT_Float64,
1577 : static_cast<int>(sizeof(double)),
1578 : static_cast<GPtrDiff_t>(nVals));
1579 : }
1580 912 : for (size_t i = 0; i < nVals; i++)
1581 : {
1582 899 : if (data->abyMaskData[i])
1583 : {
1584 894 : const double dfValue = data->adfData[i];
1585 894 : data->dfMin = std::min(data->dfMin, dfValue);
1586 894 : data->dfMax = std::max(data->dfMax, dfValue);
1587 894 : data->nValidCount++;
1588 894 : const double dfDelta = dfValue - data->dfMean;
1589 894 : data->dfMean += dfDelta / data->nValidCount;
1590 894 : data->dfM2 += dfDelta * (dfValue - data->dfMean);
1591 : }
1592 : }
1593 13 : if (data->pfnProgress &&
1594 0 : !data->pfnProgress(static_cast<double>(iCurChunk + 1) / nChunkCount,
1595 : "", data->pProgressData))
1596 : {
1597 0 : return false;
1598 : }
1599 13 : return true;
1600 : };
1601 :
1602 13 : const auto &oType = GetDataType();
1603 26 : if (oType.GetClass() != GEDTC_NUMERIC ||
1604 13 : GDALDataTypeIsComplex(oType.GetNumericDataType()))
1605 : {
1606 0 : CPLError(
1607 : CE_Failure, CPLE_NotSupported,
1608 : "Statistics can only be computed on non-complex numeric data type");
1609 0 : return false;
1610 : }
1611 :
1612 13 : const size_t nDims = GetDimensionCount();
1613 26 : std::vector<GUInt64> arrayStartIdx(nDims);
1614 26 : std::vector<GUInt64> count(nDims);
1615 13 : const auto &poDims = GetDimensions();
1616 34 : for (size_t i = 0; i < nDims; i++)
1617 : {
1618 21 : count[i] = poDims[i]->GetSize();
1619 : }
1620 13 : const char *pszSwathSize = CPLGetConfigOption("GDAL_SWATH_SIZE", nullptr);
1621 : const size_t nMaxChunkSize =
1622 : pszSwathSize
1623 13 : ? static_cast<size_t>(
1624 0 : std::min(GIntBig(std::numeric_limits<size_t>::max() / 2),
1625 0 : CPLAtoGIntBig(pszSwathSize)))
1626 : : static_cast<size_t>(
1627 13 : std::min(GIntBig(std::numeric_limits<size_t>::max() / 2),
1628 13 : GDALGetCacheMax64() / 4));
1629 26 : StatsPerChunkType sData;
1630 13 : sData.array = this;
1631 13 : sData.poMask = GetMask(nullptr);
1632 13 : if (sData.poMask == nullptr)
1633 : {
1634 0 : return false;
1635 : }
1636 13 : sData.pfnProgress = pfnProgress;
1637 13 : sData.pProgressData = pProgressData;
1638 13 : if (!ProcessPerChunk(arrayStartIdx.data(), count.data(),
1639 26 : GetProcessingChunkSize(nMaxChunkSize).data(),
1640 13 : PerChunkFunc, &sData))
1641 : {
1642 0 : return false;
1643 : }
1644 :
1645 13 : if (pdfMin)
1646 13 : *pdfMin = sData.dfMin;
1647 :
1648 13 : if (pdfMax)
1649 13 : *pdfMax = sData.dfMax;
1650 :
1651 13 : if (pdfMean)
1652 11 : *pdfMean = sData.dfMean;
1653 :
1654 : const double dfStdDev =
1655 13 : sData.nValidCount > 0 ? sqrt(sData.dfM2 / sData.nValidCount) : 0.0;
1656 13 : if (pdfStdDev)
1657 11 : *pdfStdDev = dfStdDev;
1658 :
1659 13 : if (pnValidCount)
1660 11 : *pnValidCount = sData.nValidCount;
1661 :
1662 13 : SetStatistics(bApproxOK, sData.dfMin, sData.dfMax, sData.dfMean, dfStdDev,
1663 13 : sData.nValidCount, papszOptions);
1664 :
1665 13 : return true;
1666 : }
1667 :
1668 : /************************************************************************/
1669 : /* SetStatistics() */
1670 : /************************************************************************/
1671 : //! @cond Doxygen_Suppress
1672 5 : bool GDALMDArray::SetStatistics(bool /* bApproxStats */, double /* dfMin */,
1673 : double /* dfMax */, double /* dfMean */,
1674 : double /* dfStdDev */,
1675 : GUInt64 /* nValidCount */,
1676 : CSLConstList /* papszOptions */)
1677 : {
1678 5 : CPLDebug("GDAL", "Cannot save statistics on a non-PAM MDArray");
1679 5 : return false;
1680 : }
1681 :
1682 : //! @endcond
1683 :
1684 : /************************************************************************/
1685 : /* ClearStatistics() */
1686 : /************************************************************************/
1687 :
1688 : /**
1689 : * \brief Clear statistics.
1690 : *
1691 : * @since GDAL 3.4
1692 : */
1693 0 : void GDALMDArray::ClearStatistics()
1694 : {
1695 0 : }
1696 :
1697 : /************************************************************************/
1698 : /* GetCoordinateVariables() */
1699 : /************************************************************************/
1700 :
1701 : /**
1702 : * \brief Return coordinate variables.
1703 : *
1704 : * Coordinate variables are an alternate way of indexing an array that can
1705 : * be sometimes used. For example, an array collected through remote sensing
1706 : * might be indexed by (scanline, pixel). But there can be
1707 : * a longitude and latitude arrays alongside that are also both indexed by
1708 : * (scanline, pixel), and are referenced from operational arrays for
1709 : * reprojection purposes.
1710 : *
1711 : * For netCDF, this will return the arrays referenced by the "coordinates"
1712 : * attribute.
1713 : *
1714 : * This method is the same as the C function
1715 : * GDALMDArrayGetCoordinateVariables().
1716 : *
1717 : * @return a vector of arrays
1718 : *
1719 : * @since GDAL 3.4
1720 : */
1721 :
1722 : std::vector<std::shared_ptr<GDALMDArray>>
1723 13 : GDALMDArray::GetCoordinateVariables() const
1724 : {
1725 13 : return {};
1726 : }
|
