Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Zarr driver
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2021, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "zarr.h"
14 :
15 : #include "cpl_float.h"
16 : #include "cpl_mem_cache.h"
17 : #include "cpl_multiproc.h"
18 : #include "cpl_vsi_virtual.h"
19 : #include "ucs4_utf8.hpp"
20 :
21 : #include "netcdf_cf_constants.h" // for CF_UNITS, etc
22 :
23 : #include <algorithm>
24 : #include <cassert>
25 : #include <cmath>
26 : #include <cstdlib>
27 : #include <limits>
28 : #include <map>
29 : #include <set>
30 :
31 : #if defined(__clang__) || defined(_MSC_VER)
32 : #define COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT
33 : #endif
34 :
35 : namespace
36 : {
37 :
38 4 : inline std::vector<GByte> UTF8ToUCS4(const char *pszStr, bool needByteSwap)
39 : {
40 4 : const size_t nLen = strlen(pszStr);
41 : // Worst case if that we need 4 more bytes than the UTF-8 one
42 : // (when the content is pure ASCII)
43 4 : if (nLen > std::numeric_limits<size_t>::max() / sizeof(uint32_t))
44 0 : throw std::bad_alloc();
45 4 : std::vector<GByte> ret(nLen * sizeof(uint32_t));
46 4 : size_t outPos = 0;
47 27 : for (size_t i = 0; i < nLen; outPos += sizeof(uint32_t))
48 : {
49 23 : uint32_t ucs4 = 0;
50 23 : int consumed = FcUtf8ToUcs4(
51 23 : reinterpret_cast<const uint8_t *>(pszStr + i), &ucs4, nLen - i);
52 23 : if (consumed <= 0)
53 : {
54 0 : ret.resize(outPos);
55 : }
56 23 : if (needByteSwap)
57 : {
58 1 : CPL_SWAP32PTR(&ucs4);
59 : }
60 23 : memcpy(&ret[outPos], &ucs4, sizeof(uint32_t));
61 23 : i += consumed;
62 : }
63 4 : ret.resize(outPos);
64 4 : return ret;
65 : }
66 :
67 7 : inline char *UCS4ToUTF8(const uint8_t *ucs4Ptr, size_t nSize, bool needByteSwap)
68 : {
69 : // A UCS4 char can require up to 6 bytes in UTF8.
70 7 : if (nSize > (std::numeric_limits<size_t>::max() - 1) / 6 * 4)
71 0 : return nullptr;
72 7 : const size_t nOutSize = nSize / 4 * 6 + 1;
73 7 : char *ret = static_cast<char *>(VSI_MALLOC_VERBOSE(nOutSize));
74 7 : if (ret == nullptr)
75 0 : return nullptr;
76 7 : size_t outPos = 0;
77 32 : for (size_t i = 0; i + sizeof(uint32_t) - 1 < nSize; i += sizeof(uint32_t))
78 : {
79 : uint32_t ucs4;
80 25 : memcpy(&ucs4, ucs4Ptr + i, sizeof(uint32_t));
81 25 : if (needByteSwap)
82 : {
83 2 : CPL_SWAP32PTR(&ucs4);
84 : }
85 : int written =
86 25 : FcUcs4ToUtf8(ucs4, reinterpret_cast<uint8_t *>(ret + outPos));
87 25 : outPos += written;
88 : }
89 7 : ret[outPos] = 0;
90 7 : return ret;
91 : }
92 :
93 : } // namespace
94 :
95 : /************************************************************************/
96 : /* ZarrArray::ParseChunkSize() */
97 : /************************************************************************/
98 :
99 1781 : /* static */ bool ZarrArray::ParseChunkSize(const CPLJSONArray &oChunks,
100 : const GDALExtendedDataType &oType,
101 : std::vector<GUInt64> &anBlockSize)
102 : {
103 1781 : size_t nBlockSize = oType.GetSize();
104 4846 : for (const auto &item : oChunks)
105 : {
106 3068 : const auto nSize = static_cast<GUInt64>(item.ToLong());
107 3068 : if (nSize == 0)
108 : {
109 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid content for chunks");
110 3 : return false;
111 : }
112 3067 : if (nBlockSize > std::numeric_limits<size_t>::max() / nSize)
113 : {
114 2 : CPLError(CE_Failure, CPLE_AppDefined, "Too large chunks");
115 2 : return false;
116 : }
117 3065 : nBlockSize *= static_cast<size_t>(nSize);
118 3065 : anBlockSize.emplace_back(nSize);
119 : }
120 :
121 1778 : return true;
122 : }
123 :
124 : /************************************************************************/
125 : /* ZarrArray::ComputeBlockCount() */
126 : /************************************************************************/
127 :
128 2173 : /* static */ uint64_t ZarrArray::ComputeBlockCount(
129 : const std::string &osName,
130 : const std::vector<std::shared_ptr<GDALDimension>> &aoDims,
131 : const std::vector<GUInt64> &anBlockSize)
132 : {
133 2173 : uint64_t nTotalBlockCount = 1;
134 5811 : for (size_t i = 0; i < aoDims.size(); ++i)
135 : {
136 : const uint64_t nBlockThisDim =
137 3640 : cpl::div_round_up(aoDims[i]->GetSize(), anBlockSize[i]);
138 7280 : if (nBlockThisDim != 0 &&
139 : nTotalBlockCount >
140 3640 : std::numeric_limits<uint64_t>::max() / nBlockThisDim)
141 : {
142 2 : CPLError(
143 : CE_Failure, CPLE_NotSupported,
144 : "Array %s has more than 2^64 blocks. This is not supported.",
145 : osName.c_str());
146 2 : return 0;
147 : }
148 3638 : nTotalBlockCount *= nBlockThisDim;
149 : }
150 2171 : return nTotalBlockCount;
151 : }
152 :
153 : /************************************************************************/
154 : /* ComputeCountInnerBlockInOuter() */
155 : /************************************************************************/
156 :
157 : static std::vector<GUInt64>
158 2173 : ComputeCountInnerBlockInOuter(const std::vector<GUInt64> &anInnerBlockSize,
159 : const std::vector<GUInt64> &anOuterBlockSize)
160 : {
161 2173 : std::vector<GUInt64> ret;
162 2173 : CPLAssert(anInnerBlockSize.size() == anOuterBlockSize.size());
163 5813 : for (size_t i = 0; i < anInnerBlockSize.size(); ++i)
164 : {
165 : // All those assertions must be checked by the caller before
166 : // constructing the ZarrArray instance.
167 3640 : CPLAssert(anInnerBlockSize[i] > 0);
168 3640 : CPLAssert(anInnerBlockSize[i] <= anOuterBlockSize[i]);
169 3640 : CPLAssert((anOuterBlockSize[i] % anInnerBlockSize[i]) == 0);
170 3640 : ret.push_back(anOuterBlockSize[i] / anInnerBlockSize[i]);
171 : }
172 2173 : return ret;
173 : }
174 :
175 : /************************************************************************/
176 : /* ComputeInnerBlockSizeBytes() */
177 : /************************************************************************/
178 :
179 : static size_t
180 2171 : ComputeInnerBlockSizeBytes(const std::vector<DtypeElt> &aoDtypeElts,
181 : const std::vector<GUInt64> &anInnerBlockSize)
182 : {
183 : const size_t nSourceSize =
184 2171 : aoDtypeElts.back().nativeOffset + aoDtypeElts.back().nativeSize;
185 2171 : size_t nInnerBlockSizeBytes = nSourceSize;
186 5805 : for (const auto &nBlockSize : anInnerBlockSize)
187 : {
188 : // Given that ParseChunkSize() has checked that the outer block size
189 : // fits on size_t, and that m_anInnerBlockSize[i] <= m_anOuterBlockSize[i],
190 : // this cast is safe, and the multiplication cannot overflow.
191 3634 : nInnerBlockSizeBytes *= static_cast<size_t>(nBlockSize);
192 : }
193 2171 : return nInnerBlockSizeBytes;
194 : }
195 :
196 : /************************************************************************/
197 : /* ZarrArray::ZarrArray() */
198 : /************************************************************************/
199 :
200 2173 : ZarrArray::ZarrArray(
201 : const std::shared_ptr<ZarrSharedResource> &poSharedResource,
202 : const std::shared_ptr<ZarrGroupBase> &poParent, const std::string &osName,
203 : const std::vector<std::shared_ptr<GDALDimension>> &aoDims,
204 : const GDALExtendedDataType &oType, const std::vector<DtypeElt> &aoDtypeElts,
205 : const std::vector<GUInt64> &anOuterBlockSize,
206 0 : const std::vector<GUInt64> &anInnerBlockSize)
207 : :
208 : #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
209 : GDALAbstractMDArray(poParent->GetFullName(), osName),
210 : #endif
211 2173 : GDALPamMDArray(poParent->GetFullName(), osName,
212 : poSharedResource->GetPAM()),
213 : m_poSharedResource(poSharedResource), m_poParent(poParent),
214 : m_aoDims(aoDims), m_oType(oType), m_aoDtypeElts(aoDtypeElts),
215 : m_anOuterBlockSize(anOuterBlockSize),
216 : m_anInnerBlockSize(anInnerBlockSize),
217 : m_anCountInnerBlockInOuter(ComputeCountInnerBlockInOuter(
218 2173 : m_anInnerBlockSize, m_anOuterBlockSize)),
219 : m_nTotalInnerChunkCount(
220 2173 : ComputeBlockCount(osName, aoDims, m_anInnerBlockSize)),
221 : m_nInnerBlockSizeBytes(
222 2173 : m_nTotalInnerChunkCount > 0
223 2173 : ? ComputeInnerBlockSizeBytes(m_aoDtypeElts, m_anInnerBlockSize)
224 : : 0),
225 2173 : m_oAttrGroup(m_osFullName, /*bContainerIsGroup=*/false),
226 2173 : m_bUseOptimizedCodePaths(CPLTestBool(
227 10865 : CPLGetConfigOption("GDAL_ZARR_USE_OPTIMIZED_CODE_PATHS", "YES")))
228 : {
229 2173 : }
230 :
231 : /************************************************************************/
232 : /* ~ZarrArray() */
233 : /************************************************************************/
234 :
235 2173 : ZarrArray::~ZarrArray()
236 : {
237 2173 : if (m_pabyNoData)
238 : {
239 1083 : m_oType.FreeDynamicMemory(&m_pabyNoData[0]);
240 1083 : CPLFree(m_pabyNoData);
241 : }
242 :
243 2173 : DeallocateDecodedBlockData();
244 2173 : }
245 :
246 : /************************************************************************/
247 : /* ZarrArray::SerializeSpecialAttributes() */
248 : /************************************************************************/
249 :
250 392 : CPLJSONObject ZarrArray::SerializeSpecialAttributes()
251 : {
252 392 : m_bSRSModified = false;
253 392 : m_oAttrGroup.UnsetModified();
254 :
255 392 : auto oAttrs = m_oAttrGroup.Serialize();
256 :
257 : const bool bUseSpatialProjConventions =
258 392 : EQUAL(m_aosCreationOptions.FetchNameValueDef(
259 : "GEOREFERENCING_CONVENTION", "GDAL"),
260 : "SPATIAL_PROJ");
261 :
262 41 : const auto ExportToWkt2AndPROJJSON = [this](CPLJSONObject &oContainer,
263 : const char *pszWKT2AttrName,
264 82 : const char *pszPROJJSONAttrName)
265 : {
266 41 : const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
267 41 : char *pszWKT = nullptr;
268 41 : if (m_poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
269 : {
270 41 : oContainer.Set(pszWKT2AttrName, pszWKT);
271 : }
272 41 : CPLFree(pszWKT);
273 :
274 : {
275 82 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
276 41 : char *projjson = nullptr;
277 82 : if (m_poSRS->exportToPROJJSON(&projjson, nullptr) == OGRERR_NONE &&
278 41 : projjson != nullptr)
279 : {
280 82 : CPLJSONDocument oDocProjJSON;
281 41 : if (oDocProjJSON.LoadMemory(std::string(projjson)))
282 : {
283 41 : oContainer.Set(pszPROJJSONAttrName, oDocProjJSON.GetRoot());
284 : }
285 : }
286 41 : CPLFree(projjson);
287 : }
288 41 : };
289 :
290 784 : CPLJSONArray oZarrConventionsArray;
291 392 : if (bUseSpatialProjConventions)
292 : {
293 2 : if (m_poSRS)
294 : {
295 4 : CPLJSONObject oConventionProj;
296 2 : oConventionProj.Set(
297 : "schema_url",
298 : "https://raw.githubusercontent.com/zarr-experimental/geo-proj/"
299 : "refs/tags/v1/schema.json");
300 2 : oConventionProj.Set("spec_url",
301 : "https://github.com/zarr-experimental/geo-proj/"
302 : "blob/v1/README.md");
303 2 : oConventionProj.Set("uuid", "f17cb550-5864-4468-aeb7-f3180cfb622f");
304 2 : oConventionProj.Set("name", "proj:"); // ending colon intended
305 2 : oConventionProj.Set(
306 : "description",
307 : "Coordinate reference system information for geospatial data");
308 :
309 2 : oZarrConventionsArray.Add(oConventionProj);
310 :
311 2 : const char *pszAuthorityName = m_poSRS->GetAuthorityName(nullptr);
312 2 : const char *pszAuthorityCode = m_poSRS->GetAuthorityCode(nullptr);
313 2 : if (pszAuthorityName && pszAuthorityCode)
314 : {
315 1 : oAttrs.Set("proj:code", CPLSPrintf("%s:%s", pszAuthorityName,
316 : pszAuthorityCode));
317 : }
318 : else
319 : {
320 1 : ExportToWkt2AndPROJJSON(oAttrs, "proj:wkt2", "proj:projjson");
321 : }
322 : }
323 :
324 2 : if (GetDimensionCount() >= 2)
325 : {
326 2 : bool bAddSpatialProjConvention = false;
327 :
328 2 : double dfXOff = std::numeric_limits<double>::quiet_NaN();
329 2 : double dfXRes = std::numeric_limits<double>::quiet_NaN();
330 2 : double dfYOff = std::numeric_limits<double>::quiet_NaN();
331 2 : double dfYRes = std::numeric_limits<double>::quiet_NaN();
332 4 : std::string osDimXName;
333 4 : std::string osDimYName;
334 4 : std::string osDimZName;
335 2 : double dfWidth = 0, dfHeight = 0;
336 6 : for (const auto &poDim : GetDimensions())
337 : {
338 4 : if (poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_X)
339 : {
340 2 : osDimXName = poDim->GetName();
341 2 : dfWidth = static_cast<double>(poDim->GetSize());
342 4 : auto poVar = poDim->GetIndexingVariable();
343 2 : if (poVar && poVar->IsRegularlySpaced(dfXOff, dfXRes))
344 : {
345 1 : dfXOff -= dfXRes / 2;
346 : }
347 : }
348 2 : else if (poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_Y)
349 : {
350 2 : osDimYName = poDim->GetName();
351 2 : dfHeight = static_cast<double>(poDim->GetSize());
352 4 : auto poVar = poDim->GetIndexingVariable();
353 2 : if (poVar && poVar->IsRegularlySpaced(dfYOff, dfYRes))
354 : {
355 1 : dfYOff -= dfYRes / 2;
356 : }
357 : }
358 0 : else if (poDim->GetType() == GDAL_DIM_TYPE_VERTICAL)
359 : {
360 0 : osDimZName = poDim->GetName();
361 : }
362 : }
363 :
364 2 : GDALGeoTransform gt;
365 2 : if (!osDimXName.empty() && !osDimYName.empty())
366 : {
367 6 : const auto oGDALGeoTransform = oAttrs["gdal:geotransform"];
368 : const bool bHasGDALGeoTransform =
369 4 : (oGDALGeoTransform.GetType() ==
370 4 : CPLJSONObject::Type::Array &&
371 4 : oGDALGeoTransform.ToArray().size() == 6);
372 2 : if (bHasGDALGeoTransform)
373 : {
374 : const auto oGDALGeoTransformArray =
375 2 : oGDALGeoTransform.ToArray();
376 14 : for (int i = 0; i < 6; ++i)
377 : {
378 12 : gt[i] = oGDALGeoTransformArray[i].ToDouble();
379 : }
380 2 : bAddSpatialProjConvention = true;
381 : }
382 0 : else if (!std::isnan(dfXOff) && !std::isnan(dfXRes) &&
383 0 : !std::isnan(dfYOff) && !std::isnan(dfYRes))
384 : {
385 0 : gt.xorig = dfXOff;
386 0 : gt.xscale = dfXRes;
387 0 : gt.xrot = 0; // xrot
388 0 : gt.yorig = dfYOff;
389 0 : gt.yrot = 0; // yrot
390 0 : gt.yscale = dfYRes;
391 0 : bAddSpatialProjConvention = true;
392 : }
393 : }
394 :
395 2 : if (bAddSpatialProjConvention)
396 : {
397 : const auto osGDALMD_AREA_OR_POINT =
398 6 : oAttrs.GetString(GDALMD_AREA_OR_POINT);
399 2 : if (osGDALMD_AREA_OR_POINT == GDALMD_AOP_AREA)
400 : {
401 1 : oAttrs.Add("spatial:registration", "pixel");
402 1 : oAttrs.Delete(GDALMD_AREA_OR_POINT);
403 : }
404 1 : else if (osGDALMD_AREA_OR_POINT == GDALMD_AOP_POINT)
405 : {
406 1 : oAttrs.Add("spatial:registration", "node");
407 1 : oAttrs.Delete(GDALMD_AREA_OR_POINT);
408 :
409 : // Going from GDAL's corner convention to pixel center
410 1 : gt.xorig += 0.5 * gt.xscale + 0.5 * gt.xrot;
411 1 : gt.yorig += 0.5 * gt.yrot + 0.5 * gt.yscale;
412 1 : dfWidth -= 1.0;
413 1 : dfHeight -= 1.0;
414 : }
415 :
416 4 : CPLJSONArray oAttrSpatialTransform;
417 2 : oAttrSpatialTransform.Add(gt.xscale); // xres
418 2 : oAttrSpatialTransform.Add(gt.xrot); // xrot
419 2 : oAttrSpatialTransform.Add(gt.xorig); // xoff
420 2 : oAttrSpatialTransform.Add(gt.yrot); // yrot
421 2 : oAttrSpatialTransform.Add(gt.yscale); // yres
422 2 : oAttrSpatialTransform.Add(gt.yorig); // yoff
423 :
424 2 : oAttrs.Add("spatial:transform_type", "affine");
425 2 : oAttrs.Add("spatial:transform", oAttrSpatialTransform);
426 2 : oAttrs.Delete("gdal:geotransform");
427 :
428 : double dfX0, dfY0;
429 : double dfX1, dfY1;
430 : double dfX2, dfY2;
431 : double dfX3, dfY3;
432 2 : gt.Apply(0, 0, &dfX0, &dfY0);
433 2 : gt.Apply(dfWidth, 0, &dfX1, &dfY1);
434 2 : gt.Apply(0, dfHeight, &dfX2, &dfY2);
435 2 : gt.Apply(dfWidth, dfHeight, &dfX3, &dfY3);
436 : const double dfXMin =
437 2 : std::min(std::min(dfX0, dfX1), std::min(dfX2, dfX3));
438 : const double dfYMin =
439 2 : std::min(std::min(dfY0, dfY1), std::min(dfY2, dfY3));
440 : const double dfXMax =
441 2 : std::max(std::max(dfX0, dfX1), std::max(dfX2, dfX3));
442 : const double dfYMax =
443 2 : std::max(std::max(dfY0, dfY1), std::max(dfY2, dfY3));
444 :
445 4 : CPLJSONArray oAttrSpatialBBOX;
446 2 : oAttrSpatialBBOX.Add(dfXMin);
447 2 : oAttrSpatialBBOX.Add(dfYMin);
448 2 : oAttrSpatialBBOX.Add(dfXMax);
449 2 : oAttrSpatialBBOX.Add(dfYMax);
450 2 : oAttrs.Add("spatial:bbox", oAttrSpatialBBOX);
451 :
452 4 : CPLJSONArray aoSpatialDimensions;
453 2 : if (!osDimZName.empty())
454 0 : aoSpatialDimensions.Add(osDimZName);
455 2 : aoSpatialDimensions.Add(osDimYName);
456 2 : aoSpatialDimensions.Add(osDimXName);
457 2 : oAttrs.Add("spatial:dimensions", aoSpatialDimensions);
458 :
459 4 : CPLJSONObject oConventionSpatial;
460 2 : oConventionSpatial.Set(
461 : "schema_url",
462 : "https://raw.githubusercontent.com/zarr-conventions/"
463 : "spatial/refs/tags/v1/schema.json");
464 2 : oConventionSpatial.Set("spec_url",
465 : "https://github.com/zarr-conventions/"
466 : "spatial/blob/v1/README.md");
467 2 : oConventionSpatial.Set("uuid",
468 : "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4");
469 2 : oConventionSpatial.Set("name",
470 : "spatial:"); // ending colon intended
471 2 : oConventionSpatial.Set("description",
472 : "Spatial coordinate information");
473 :
474 2 : oZarrConventionsArray.Add(oConventionSpatial);
475 : }
476 : }
477 :
478 2 : if (oZarrConventionsArray.size() > 0)
479 : {
480 2 : oAttrs.Add("zarr_conventions", oZarrConventionsArray);
481 : }
482 : }
483 390 : else if (m_poSRS)
484 : {
485 : // GDAL convention
486 :
487 40 : CPLJSONObject oCRS;
488 40 : ExportToWkt2AndPROJJSON(oCRS, "wkt", "projjson");
489 :
490 40 : const char *pszAuthorityCode = m_poSRS->GetAuthorityCode(nullptr);
491 40 : const char *pszAuthorityName = m_poSRS->GetAuthorityName(nullptr);
492 40 : if (pszAuthorityCode && pszAuthorityName &&
493 7 : EQUAL(pszAuthorityName, "EPSG"))
494 : {
495 7 : oCRS.Add("url",
496 14 : std::string("http://www.opengis.net/def/crs/EPSG/0/") +
497 : pszAuthorityCode);
498 : }
499 :
500 40 : oAttrs.Add(CRS_ATTRIBUTE_NAME, oCRS);
501 : }
502 :
503 392 : if (m_osUnit.empty())
504 : {
505 383 : if (m_bUnitModified)
506 0 : oAttrs.Delete(CF_UNITS);
507 : }
508 : else
509 : {
510 9 : oAttrs.Set(CF_UNITS, m_osUnit);
511 : }
512 392 : m_bUnitModified = false;
513 :
514 392 : if (!m_bHasOffset)
515 : {
516 389 : oAttrs.Delete(CF_ADD_OFFSET);
517 : }
518 : else
519 : {
520 3 : oAttrs.Set(CF_ADD_OFFSET, m_dfOffset);
521 : }
522 392 : m_bOffsetModified = false;
523 :
524 392 : if (!m_bHasScale)
525 : {
526 389 : oAttrs.Delete(CF_SCALE_FACTOR);
527 : }
528 : else
529 : {
530 3 : oAttrs.Set(CF_SCALE_FACTOR, m_dfScale);
531 : }
532 392 : m_bScaleModified = false;
533 :
534 784 : return oAttrs;
535 : }
536 :
537 : /************************************************************************/
538 : /* FillBlockSize() */
539 : /************************************************************************/
540 :
541 : /* static */
542 450 : bool ZarrArray::FillBlockSize(
543 : const std::vector<std::shared_ptr<GDALDimension>> &aoDimensions,
544 : const GDALExtendedDataType &oDataType, std::vector<GUInt64> &anBlockSize,
545 : CSLConstList papszOptions)
546 : {
547 450 : const auto nDims = aoDimensions.size();
548 450 : anBlockSize.resize(nDims);
549 1118 : for (size_t i = 0; i < nDims; ++i)
550 668 : anBlockSize[i] = 1;
551 450 : if (nDims >= 2)
552 : {
553 440 : anBlockSize[nDims - 2] =
554 440 : std::min(std::max<GUInt64>(1, aoDimensions[nDims - 2]->GetSize()),
555 440 : static_cast<GUInt64>(256));
556 440 : anBlockSize[nDims - 1] =
557 440 : std::min(std::max<GUInt64>(1, aoDimensions[nDims - 1]->GetSize()),
558 660 : static_cast<GUInt64>(256));
559 : }
560 230 : else if (nDims == 1)
561 : {
562 183 : anBlockSize[0] = std::max<GUInt64>(1, aoDimensions[0]->GetSize());
563 : }
564 :
565 450 : const char *pszBlockSize = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
566 450 : if (pszBlockSize)
567 : {
568 : const auto aszTokens(
569 18 : CPLStringList(CSLTokenizeString2(pszBlockSize, ",", 0)));
570 18 : if (static_cast<size_t>(aszTokens.size()) != nDims)
571 : {
572 5 : CPLError(CE_Failure, CPLE_AppDefined,
573 : "Invalid number of values in BLOCKSIZE");
574 5 : return false;
575 : }
576 13 : size_t nBlockSize = oDataType.GetSize();
577 40 : for (size_t i = 0; i < nDims; ++i)
578 : {
579 28 : const auto v = static_cast<GUInt64>(CPLAtoGIntBig(aszTokens[i]));
580 28 : if (v > 0)
581 : {
582 28 : anBlockSize[i] = v;
583 : }
584 28 : if (anBlockSize[i] >
585 28 : std::numeric_limits<size_t>::max() / nBlockSize)
586 : {
587 1 : CPLError(CE_Failure, CPLE_AppDefined,
588 : "Too large values in BLOCKSIZE");
589 1 : return false;
590 : }
591 27 : nBlockSize *= static_cast<size_t>(anBlockSize[i]);
592 : }
593 : }
594 444 : return true;
595 : }
596 :
597 : /************************************************************************/
598 : /* DeallocateDecodedBlockData() */
599 : /************************************************************************/
600 :
601 3680 : void ZarrArray::DeallocateDecodedBlockData()
602 : {
603 3680 : if (!m_abyDecodedBlockData.empty())
604 : {
605 223 : const size_t nDTSize = m_oType.GetSize();
606 223 : GByte *pDst = &m_abyDecodedBlockData[0];
607 223 : const size_t nValues = m_abyDecodedBlockData.size() / nDTSize;
608 454 : for (const auto &elt : m_aoDtypeElts)
609 : {
610 231 : if (elt.nativeType == DtypeElt::NativeType::STRING_ASCII ||
611 230 : elt.nativeType == DtypeElt::NativeType::STRING_UNICODE)
612 : {
613 2 : for (size_t i = 0; i < nValues; i++, pDst += nDTSize)
614 : {
615 : char *ptr;
616 1 : char **pptr =
617 1 : reinterpret_cast<char **>(pDst + elt.gdalOffset);
618 1 : memcpy(&ptr, pptr, sizeof(ptr));
619 1 : VSIFree(ptr);
620 : }
621 : }
622 : }
623 : }
624 3680 : }
625 :
626 : /************************************************************************/
627 : /* EncodeElt() */
628 : /************************************************************************/
629 :
630 : /* Encode from GDAL raw type to Zarr native type */
631 : /*static*/
632 121 : void ZarrArray::EncodeElt(const std::vector<DtypeElt> &elts, const GByte *pSrc,
633 : GByte *pDst)
634 : {
635 243 : for (const auto &elt : elts)
636 : {
637 122 : if (elt.nativeType == DtypeElt::NativeType::STRING_UNICODE)
638 : {
639 0 : const char *pStr =
640 0 : *reinterpret_cast<const char *const *>(pSrc + elt.gdalOffset);
641 0 : if (pStr)
642 : {
643 : try
644 : {
645 0 : const auto ucs4 = UTF8ToUCS4(pStr, elt.needByteSwapping);
646 0 : const auto ucs4Len = ucs4.size();
647 0 : memcpy(pDst + elt.nativeOffset, ucs4.data(),
648 0 : std::min(ucs4Len, elt.nativeSize));
649 0 : if (ucs4Len > elt.nativeSize)
650 : {
651 0 : CPLError(CE_Warning, CPLE_AppDefined,
652 : "Too long string truncated");
653 : }
654 0 : else if (ucs4Len < elt.nativeSize)
655 : {
656 0 : memset(pDst + elt.nativeOffset + ucs4Len, 0,
657 0 : elt.nativeSize - ucs4Len);
658 : }
659 : }
660 0 : catch (const std::exception &)
661 : {
662 0 : memset(pDst + elt.nativeOffset, 0, elt.nativeSize);
663 : }
664 : }
665 : else
666 : {
667 0 : memset(pDst + elt.nativeOffset, 0, elt.nativeSize);
668 : }
669 : }
670 122 : else if (elt.needByteSwapping)
671 : {
672 120 : if (elt.nativeSize == 2)
673 : {
674 24 : if (elt.gdalTypeIsApproxOfNative)
675 : {
676 0 : CPLAssert(elt.nativeType == DtypeElt::NativeType::IEEEFP);
677 0 : CPLAssert(elt.gdalType.GetNumericDataType() == GDT_Float32);
678 0 : const uint32_t uint32Val =
679 0 : *reinterpret_cast<const uint32_t *>(pSrc +
680 0 : elt.gdalOffset);
681 0 : bool bHasWarned = false;
682 : uint16_t uint16Val =
683 0 : CPL_SWAP16(CPLFloatToHalf(uint32Val, bHasWarned));
684 0 : memcpy(pDst + elt.nativeOffset, &uint16Val,
685 : sizeof(uint16Val));
686 : }
687 : else
688 : {
689 24 : const uint16_t val =
690 24 : CPL_SWAP16(*reinterpret_cast<const uint16_t *>(
691 : pSrc + elt.gdalOffset));
692 24 : memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
693 : }
694 : }
695 96 : else if (elt.nativeSize == 4)
696 : {
697 36 : const uint32_t val = CPL_SWAP32(
698 : *reinterpret_cast<const uint32_t *>(pSrc + elt.gdalOffset));
699 36 : memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
700 : }
701 60 : else if (elt.nativeSize == 8)
702 : {
703 48 : if (elt.nativeType == DtypeElt::NativeType::COMPLEX_IEEEFP)
704 : {
705 12 : uint32_t val =
706 12 : CPL_SWAP32(*reinterpret_cast<const uint32_t *>(
707 : pSrc + elt.gdalOffset));
708 12 : memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
709 12 : val = CPL_SWAP32(*reinterpret_cast<const uint32_t *>(
710 : pSrc + elt.gdalOffset + 4));
711 12 : memcpy(pDst + elt.nativeOffset + 4, &val, sizeof(val));
712 : }
713 : else
714 : {
715 36 : const uint64_t val =
716 36 : CPL_SWAP64(*reinterpret_cast<const uint64_t *>(
717 : pSrc + elt.gdalOffset));
718 36 : memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
719 : }
720 : }
721 12 : else if (elt.nativeSize == 16)
722 : {
723 12 : uint64_t val = CPL_SWAP64(
724 : *reinterpret_cast<const uint64_t *>(pSrc + elt.gdalOffset));
725 12 : memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
726 12 : val = CPL_SWAP64(*reinterpret_cast<const uint64_t *>(
727 : pSrc + elt.gdalOffset + 8));
728 12 : memcpy(pDst + elt.nativeOffset + 8, &val, sizeof(val));
729 : }
730 : else
731 : {
732 0 : CPLAssert(false);
733 : }
734 : }
735 2 : else if (elt.gdalTypeIsApproxOfNative)
736 : {
737 0 : if (elt.nativeType == DtypeElt::NativeType::IEEEFP &&
738 0 : elt.nativeSize == 2)
739 : {
740 0 : CPLAssert(elt.gdalType.GetNumericDataType() == GDT_Float32);
741 0 : const uint32_t uint32Val =
742 0 : *reinterpret_cast<const uint32_t *>(pSrc + elt.gdalOffset);
743 0 : bool bHasWarned = false;
744 : const uint16_t uint16Val =
745 0 : CPLFloatToHalf(uint32Val, bHasWarned);
746 0 : memcpy(pDst + elt.nativeOffset, &uint16Val, sizeof(uint16Val));
747 : }
748 : else
749 : {
750 0 : CPLAssert(false);
751 : }
752 : }
753 2 : else if (elt.nativeType == DtypeElt::NativeType::STRING_ASCII)
754 : {
755 0 : const char *pStr =
756 0 : *reinterpret_cast<const char *const *>(pSrc + elt.gdalOffset);
757 0 : if (pStr)
758 : {
759 0 : const size_t nLen = strlen(pStr);
760 0 : memcpy(pDst + elt.nativeOffset, pStr,
761 0 : std::min(nLen, elt.nativeSize));
762 0 : if (nLen < elt.nativeSize)
763 0 : memset(pDst + elt.nativeOffset + nLen, 0,
764 0 : elt.nativeSize - nLen);
765 : }
766 : else
767 : {
768 0 : memset(pDst + elt.nativeOffset, 0, elt.nativeSize);
769 : }
770 : }
771 : else
772 : {
773 2 : CPLAssert(elt.nativeSize == elt.gdalSize);
774 2 : memcpy(pDst + elt.nativeOffset, pSrc + elt.gdalOffset,
775 2 : elt.nativeSize);
776 : }
777 : }
778 121 : }
779 :
780 : /************************************************************************/
781 : /* ZarrArray::SerializeNumericNoData() */
782 : /************************************************************************/
783 :
784 16 : void ZarrArray::SerializeNumericNoData(CPLJSONObject &oRoot) const
785 : {
786 16 : if (m_oType.GetNumericDataType() == GDT_Int64)
787 : {
788 0 : const auto nVal = GetNoDataValueAsInt64();
789 0 : oRoot.Add("fill_value", static_cast<GInt64>(nVal));
790 : }
791 16 : else if (m_oType.GetNumericDataType() == GDT_UInt64)
792 : {
793 0 : const auto nVal = GetNoDataValueAsUInt64();
794 0 : oRoot.Add("fill_value", static_cast<uint64_t>(nVal));
795 : }
796 : else
797 : {
798 16 : const double dfVal = GetNoDataValueAsDouble();
799 16 : if (std::isnan(dfVal))
800 2 : oRoot.Add("fill_value", "NaN");
801 14 : else if (dfVal == std::numeric_limits<double>::infinity())
802 2 : oRoot.Add("fill_value", "Infinity");
803 12 : else if (dfVal == -std::numeric_limits<double>::infinity())
804 2 : oRoot.Add("fill_value", "-Infinity");
805 10 : else if (GDALDataTypeIsInteger(m_oType.GetNumericDataType()))
806 8 : oRoot.Add("fill_value", static_cast<GInt64>(dfVal));
807 : else
808 2 : oRoot.Add("fill_value", dfVal);
809 : }
810 16 : }
811 :
812 : /************************************************************************/
813 : /* ZarrArray::GetSpatialRef() */
814 : /************************************************************************/
815 :
816 49 : std::shared_ptr<OGRSpatialReference> ZarrArray::GetSpatialRef() const
817 : {
818 49 : if (!CheckValidAndErrorOutIfNot())
819 0 : return nullptr;
820 :
821 49 : if (m_poSRS)
822 30 : return m_poSRS;
823 19 : return GDALPamMDArray::GetSpatialRef();
824 : }
825 :
826 : /************************************************************************/
827 : /* SetRawNoDataValue() */
828 : /************************************************************************/
829 :
830 18 : bool ZarrArray::SetRawNoDataValue(const void *pRawNoData)
831 : {
832 18 : if (!CheckValidAndErrorOutIfNot())
833 0 : return false;
834 :
835 18 : if (!m_bUpdatable)
836 : {
837 0 : CPLError(CE_Failure, CPLE_AppDefined, "Array opened in read-only mode");
838 0 : return false;
839 : }
840 18 : m_bDefinitionModified = true;
841 18 : RegisterNoDataValue(pRawNoData);
842 18 : return true;
843 : }
844 :
845 : /************************************************************************/
846 : /* RegisterNoDataValue() */
847 : /************************************************************************/
848 :
849 1083 : void ZarrArray::RegisterNoDataValue(const void *pNoData)
850 : {
851 1083 : if (m_pabyNoData)
852 : {
853 0 : m_oType.FreeDynamicMemory(&m_pabyNoData[0]);
854 : }
855 :
856 1083 : if (pNoData == nullptr)
857 : {
858 0 : CPLFree(m_pabyNoData);
859 0 : m_pabyNoData = nullptr;
860 : }
861 : else
862 : {
863 1083 : const auto nSize = m_oType.GetSize();
864 1083 : if (m_pabyNoData == nullptr)
865 : {
866 1083 : m_pabyNoData = static_cast<GByte *>(CPLMalloc(nSize));
867 : }
868 1083 : memset(m_pabyNoData, 0, nSize);
869 1083 : GDALExtendedDataType::CopyValue(pNoData, m_oType, m_pabyNoData,
870 1083 : m_oType);
871 : }
872 1083 : }
873 :
874 : /************************************************************************/
875 : /* DecodeSourceElt() */
876 : /************************************************************************/
877 :
878 : /* static */
879 1936 : void ZarrArray::DecodeSourceElt(const std::vector<DtypeElt> &elts,
880 : const GByte *pSrc, GByte *pDst)
881 : {
882 3904 : for (const auto &elt : elts)
883 : {
884 1968 : if (elt.nativeType == DtypeElt::NativeType::STRING_UNICODE)
885 : {
886 : char *ptr;
887 0 : char **pDstPtr = reinterpret_cast<char **>(pDst + elt.gdalOffset);
888 0 : memcpy(&ptr, pDstPtr, sizeof(ptr));
889 0 : VSIFree(ptr);
890 :
891 0 : char *pDstStr = UCS4ToUTF8(pSrc + elt.nativeOffset, elt.nativeSize,
892 0 : elt.needByteSwapping);
893 0 : memcpy(pDstPtr, &pDstStr, sizeof(pDstStr));
894 : }
895 1968 : else if (elt.needByteSwapping)
896 : {
897 1922 : if (elt.nativeSize == 2)
898 : {
899 : uint16_t val;
900 386 : memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
901 386 : if (elt.gdalTypeIsApproxOfNative)
902 : {
903 0 : CPLAssert(elt.nativeType == DtypeElt::NativeType::IEEEFP);
904 0 : CPLAssert(elt.gdalType.GetNumericDataType() == GDT_Float32);
905 0 : uint32_t uint32Val = CPLHalfToFloat(CPL_SWAP16(val));
906 0 : memcpy(pDst + elt.gdalOffset, &uint32Val,
907 : sizeof(uint32Val));
908 : }
909 : else
910 : {
911 386 : *reinterpret_cast<uint16_t *>(pDst + elt.gdalOffset) =
912 386 : CPL_SWAP16(val);
913 : }
914 : }
915 1536 : else if (elt.nativeSize == 4)
916 : {
917 : uint32_t val;
918 576 : memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
919 576 : *reinterpret_cast<uint32_t *>(pDst + elt.gdalOffset) =
920 576 : CPL_SWAP32(val);
921 : }
922 960 : else if (elt.nativeSize == 8)
923 : {
924 768 : if (elt.nativeType == DtypeElt::NativeType::COMPLEX_IEEEFP)
925 : {
926 : uint32_t val;
927 192 : memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
928 192 : *reinterpret_cast<uint32_t *>(pDst + elt.gdalOffset) =
929 192 : CPL_SWAP32(val);
930 192 : memcpy(&val, pSrc + elt.nativeOffset + 4, sizeof(val));
931 192 : *reinterpret_cast<uint32_t *>(pDst + elt.gdalOffset + 4) =
932 192 : CPL_SWAP32(val);
933 : }
934 : else
935 : {
936 : uint64_t val;
937 576 : memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
938 576 : *reinterpret_cast<uint64_t *>(pDst + elt.gdalOffset) =
939 576 : CPL_SWAP64(val);
940 : }
941 : }
942 192 : else if (elt.nativeSize == 16)
943 : {
944 : uint64_t val;
945 192 : memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
946 192 : *reinterpret_cast<uint64_t *>(pDst + elt.gdalOffset) =
947 192 : CPL_SWAP64(val);
948 192 : memcpy(&val, pSrc + elt.nativeOffset + 8, sizeof(val));
949 192 : *reinterpret_cast<uint64_t *>(pDst + elt.gdalOffset + 8) =
950 192 : CPL_SWAP64(val);
951 : }
952 : else
953 : {
954 0 : CPLAssert(false);
955 : }
956 : }
957 46 : else if (elt.gdalTypeIsApproxOfNative)
958 : {
959 0 : if (elt.nativeType == DtypeElt::NativeType::IEEEFP &&
960 0 : elt.nativeSize == 2)
961 : {
962 0 : CPLAssert(elt.gdalType.GetNumericDataType() == GDT_Float32);
963 : uint16_t uint16Val;
964 0 : memcpy(&uint16Val, pSrc + elt.nativeOffset, sizeof(uint16Val));
965 0 : uint32_t uint32Val = CPLHalfToFloat(uint16Val);
966 0 : memcpy(pDst + elt.gdalOffset, &uint32Val, sizeof(uint32Val));
967 : }
968 : else
969 : {
970 0 : CPLAssert(false);
971 : }
972 : }
973 46 : else if (elt.nativeType == DtypeElt::NativeType::STRING_ASCII)
974 : {
975 : char *ptr;
976 3 : char **pDstPtr = reinterpret_cast<char **>(pDst + elt.gdalOffset);
977 3 : memcpy(&ptr, pDstPtr, sizeof(ptr));
978 3 : VSIFree(ptr);
979 :
980 3 : char *pDstStr = static_cast<char *>(CPLMalloc(elt.nativeSize + 1));
981 3 : memcpy(pDstStr, pSrc + elt.nativeOffset, elt.nativeSize);
982 3 : pDstStr[elt.nativeSize] = 0;
983 3 : memcpy(pDstPtr, &pDstStr, sizeof(pDstStr));
984 : }
985 : else
986 : {
987 43 : CPLAssert(elt.nativeSize == elt.gdalSize);
988 43 : memcpy(pDst + elt.gdalOffset, pSrc + elt.nativeOffset,
989 43 : elt.nativeSize);
990 : }
991 : }
992 1936 : }
993 :
994 : /************************************************************************/
995 : /* ZarrArray::IAdviseReadCommon() */
996 : /************************************************************************/
997 :
998 13 : bool ZarrArray::IAdviseReadCommon(const GUInt64 *arrayStartIdx,
999 : const size_t *count,
1000 : CSLConstList papszOptions,
1001 : std::vector<uint64_t> &anIndicesCur,
1002 : int &nThreadsMax,
1003 : std::vector<uint64_t> &anReqBlocksIndices,
1004 : size_t &nReqBlocks) const
1005 : {
1006 13 : if (!CheckValidAndErrorOutIfNot())
1007 0 : return false;
1008 :
1009 13 : const size_t nDims = m_aoDims.size();
1010 13 : anIndicesCur.resize(nDims);
1011 26 : std::vector<uint64_t> anIndicesMin(nDims);
1012 26 : std::vector<uint64_t> anIndicesMax(nDims);
1013 :
1014 : // Compute min and max tile indices in each dimension, and the total
1015 : // number of tiles this represents.
1016 13 : nReqBlocks = 1;
1017 39 : for (size_t i = 0; i < nDims; ++i)
1018 : {
1019 26 : anIndicesMin[i] = arrayStartIdx[i] / m_anInnerBlockSize[i];
1020 52 : anIndicesMax[i] =
1021 26 : (arrayStartIdx[i] + count[i] - 1) / m_anInnerBlockSize[i];
1022 : // Overflow on number of tiles already checked in Create()
1023 26 : nReqBlocks *=
1024 26 : static_cast<size_t>(anIndicesMax[i] - anIndicesMin[i] + 1);
1025 : }
1026 :
1027 : // Find available cache size
1028 13 : const size_t nCacheSize = [papszOptions]()
1029 : {
1030 : size_t nCacheSizeTmp;
1031 : const char *pszCacheSize =
1032 13 : CSLFetchNameValue(papszOptions, "CACHE_SIZE");
1033 13 : if (pszCacheSize)
1034 : {
1035 4 : const auto nCacheSizeBig = CPLAtoGIntBig(pszCacheSize);
1036 8 : if (nCacheSizeBig < 0 || static_cast<uint64_t>(nCacheSizeBig) >
1037 4 : std::numeric_limits<size_t>::max() / 2)
1038 : {
1039 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "Too big CACHE_SIZE");
1040 0 : return std::numeric_limits<size_t>::max();
1041 : }
1042 4 : nCacheSizeTmp = static_cast<size_t>(nCacheSizeBig);
1043 : }
1044 : else
1045 : {
1046 : // Arbitrarily take half of remaining cache size
1047 9 : nCacheSizeTmp = static_cast<size_t>(std::min(
1048 18 : static_cast<uint64_t>(
1049 9 : (GDALGetCacheMax64() - GDALGetCacheUsed64()) / 2),
1050 18 : static_cast<uint64_t>(std::numeric_limits<size_t>::max() / 2)));
1051 9 : CPLDebug(ZARR_DEBUG_KEY, "Using implicit CACHE_SIZE=" CPL_FRMT_GUIB,
1052 : static_cast<GUIntBig>(nCacheSizeTmp));
1053 : }
1054 13 : return nCacheSizeTmp;
1055 13 : }();
1056 13 : if (nCacheSize == std::numeric_limits<size_t>::max())
1057 0 : return false;
1058 :
1059 : // Check that cache size is sufficient to hold all needed tiles.
1060 : // Also check that anReqBlocksIndices size computation won't overflow.
1061 13 : if (nReqBlocks > nCacheSize / std::max(m_nInnerBlockSizeBytes, nDims))
1062 : {
1063 4 : CPLError(CE_Failure, CPLE_OutOfMemory,
1064 : "CACHE_SIZE=" CPL_FRMT_GUIB " is not big enough to cache "
1065 : "all needed tiles. "
1066 : "At least " CPL_FRMT_GUIB " bytes would be needed",
1067 : static_cast<GUIntBig>(nCacheSize),
1068 : static_cast<GUIntBig>(
1069 4 : nReqBlocks * std::max(m_nInnerBlockSizeBytes, nDims)));
1070 4 : return false;
1071 : }
1072 :
1073 9 : const char *pszNumThreads = CSLFetchNameValueDef(
1074 : papszOptions, "NUM_THREADS",
1075 : CPLGetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"));
1076 9 : if (EQUAL(pszNumThreads, "ALL_CPUS"))
1077 9 : nThreadsMax = CPLGetNumCPUs();
1078 : else
1079 0 : nThreadsMax = std::max(1, atoi(pszNumThreads));
1080 9 : if (nThreadsMax > 1024)
1081 0 : nThreadsMax = 1024;
1082 9 : if (nThreadsMax <= 1)
1083 0 : return true;
1084 9 : CPLDebug(ZARR_DEBUG_KEY, "IAdviseRead(): Using up to %d threads",
1085 : nThreadsMax);
1086 :
1087 9 : m_oChunkCache.clear();
1088 :
1089 : // Overflow checked above
1090 : try
1091 : {
1092 9 : anReqBlocksIndices.resize(nDims * nReqBlocks);
1093 : }
1094 0 : catch (const std::bad_alloc &e)
1095 : {
1096 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
1097 0 : "Cannot allocate anReqBlocksIndices: %s", e.what());
1098 0 : return false;
1099 : }
1100 :
1101 9 : size_t dimIdx = 0;
1102 9 : size_t nBlockIter = 0;
1103 21639 : lbl_next_depth:
1104 21639 : if (dimIdx == nDims)
1105 : {
1106 21369 : if (nDims == 2)
1107 : {
1108 : // optimize in common case
1109 21369 : memcpy(&anReqBlocksIndices[nBlockIter * nDims], anIndicesCur.data(),
1110 : sizeof(uint64_t) * 2);
1111 : }
1112 0 : else if (nDims == 3)
1113 : {
1114 : // optimize in common case
1115 0 : memcpy(&anReqBlocksIndices[nBlockIter * nDims], anIndicesCur.data(),
1116 : sizeof(uint64_t) * 3);
1117 : }
1118 : else
1119 : {
1120 0 : memcpy(&anReqBlocksIndices[nBlockIter * nDims], anIndicesCur.data(),
1121 0 : sizeof(uint64_t) * nDims);
1122 : }
1123 21369 : nBlockIter++;
1124 : }
1125 : else
1126 : {
1127 : // This level of loop loops over blocks
1128 270 : anIndicesCur[dimIdx] = anIndicesMin[dimIdx];
1129 : while (true)
1130 : {
1131 21630 : dimIdx++;
1132 21630 : goto lbl_next_depth;
1133 21630 : lbl_return_to_caller:
1134 21630 : dimIdx--;
1135 21630 : if (anIndicesCur[dimIdx] == anIndicesMax[dimIdx])
1136 270 : break;
1137 21360 : ++anIndicesCur[dimIdx];
1138 : }
1139 : }
1140 21639 : if (dimIdx > 0)
1141 21630 : goto lbl_return_to_caller;
1142 9 : assert(nBlockIter == nReqBlocks);
1143 :
1144 9 : return true;
1145 : }
1146 :
1147 : /************************************************************************/
1148 : /* ZarrArray::IRead() */
1149 : /************************************************************************/
1150 :
1151 2285 : bool ZarrArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
1152 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
1153 : const GDALExtendedDataType &bufferDataType,
1154 : void *pDstBuffer) const
1155 : {
1156 2285 : if (!CheckValidAndErrorOutIfNot())
1157 0 : return false;
1158 :
1159 2285 : if (!AllocateWorkingBuffers())
1160 3 : return false;
1161 :
1162 : // Need to be kept in top-level scope
1163 4564 : std::vector<GUInt64> arrayStartIdxMod;
1164 4564 : std::vector<GInt64> arrayStepMod;
1165 4564 : std::vector<GPtrDiff_t> bufferStrideMod;
1166 :
1167 2282 : const size_t nDims = m_aoDims.size();
1168 2282 : bool negativeStep = false;
1169 6452 : for (size_t i = 0; i < nDims; ++i)
1170 : {
1171 4372 : if (arrayStep[i] < 0)
1172 : {
1173 202 : negativeStep = true;
1174 202 : break;
1175 : }
1176 : }
1177 :
1178 : // const auto eBufferDT = bufferDataType.GetNumericDataType();
1179 2282 : const auto nBufferDTSize = static_cast<int>(bufferDataType.GetSize());
1180 :
1181 : // Make sure that arrayStep[i] are positive for sake of simplicity
1182 2282 : if (negativeStep)
1183 : {
1184 : #if defined(__GNUC__)
1185 : #pragma GCC diagnostic push
1186 : #pragma GCC diagnostic ignored "-Wnull-dereference"
1187 : #endif
1188 202 : arrayStartIdxMod.resize(nDims);
1189 202 : arrayStepMod.resize(nDims);
1190 202 : bufferStrideMod.resize(nDims);
1191 : #if defined(__GNUC__)
1192 : #pragma GCC diagnostic pop
1193 : #endif
1194 606 : for (size_t i = 0; i < nDims; ++i)
1195 : {
1196 404 : if (arrayStep[i] < 0)
1197 : {
1198 808 : arrayStartIdxMod[i] =
1199 404 : arrayStartIdx[i] - (count[i] - 1) * (-arrayStep[i]);
1200 404 : arrayStepMod[i] = -arrayStep[i];
1201 404 : bufferStrideMod[i] = -bufferStride[i];
1202 404 : pDstBuffer =
1203 : static_cast<GByte *>(pDstBuffer) +
1204 404 : bufferStride[i] *
1205 404 : static_cast<GPtrDiff_t>(nBufferDTSize * (count[i] - 1));
1206 : }
1207 : else
1208 : {
1209 0 : arrayStartIdxMod[i] = arrayStartIdx[i];
1210 0 : arrayStepMod[i] = arrayStep[i];
1211 0 : bufferStrideMod[i] = bufferStride[i];
1212 : }
1213 : }
1214 202 : arrayStartIdx = arrayStartIdxMod.data();
1215 202 : arrayStep = arrayStepMod.data();
1216 202 : bufferStride = bufferStrideMod.data();
1217 : }
1218 :
1219 4564 : std::vector<uint64_t> indicesOuterLoop(nDims + 1);
1220 4564 : std::vector<GByte *> dstPtrStackOuterLoop(nDims + 1);
1221 :
1222 4564 : std::vector<uint64_t> indicesInnerLoop(nDims + 1);
1223 4564 : std::vector<GByte *> dstPtrStackInnerLoop(nDims + 1);
1224 :
1225 4564 : std::vector<GPtrDiff_t> dstBufferStrideBytes;
1226 6856 : for (size_t i = 0; i < nDims; ++i)
1227 : {
1228 4574 : dstBufferStrideBytes.push_back(bufferStride[i] *
1229 4574 : static_cast<GPtrDiff_t>(nBufferDTSize));
1230 : }
1231 2282 : dstBufferStrideBytes.push_back(0);
1232 :
1233 2282 : const auto nDTSize = m_oType.GetSize();
1234 :
1235 4564 : std::vector<uint64_t> blockIndices(nDims);
1236 : const size_t nSourceSize =
1237 2282 : m_aoDtypeElts.back().nativeOffset + m_aoDtypeElts.back().nativeSize;
1238 :
1239 4564 : std::vector<size_t> countInnerLoopInit(nDims + 1, 1);
1240 4564 : std::vector<size_t> countInnerLoop(nDims);
1241 :
1242 4546 : const bool bBothAreNumericDT = m_oType.GetClass() == GEDTC_NUMERIC &&
1243 2264 : bufferDataType.GetClass() == GEDTC_NUMERIC;
1244 : const bool bSameNumericDT =
1245 4546 : bBothAreNumericDT &&
1246 2264 : m_oType.GetNumericDataType() == bufferDataType.GetNumericDataType();
1247 2282 : const auto nSameDTSize = bSameNumericDT ? m_oType.GetSize() : 0;
1248 : const bool bSameCompoundAndNoDynamicMem =
1249 2285 : m_oType.GetClass() == GEDTC_COMPOUND && m_oType == bufferDataType &&
1250 3 : !m_oType.NeedsFreeDynamicMemory();
1251 4564 : std::vector<GByte> abyTargetNoData;
1252 2282 : bool bNoDataIsZero = false;
1253 :
1254 2282 : size_t dimIdx = 0;
1255 2282 : dstPtrStackOuterLoop[0] = static_cast<GByte *>(pDstBuffer);
1256 35010 : lbl_next_depth:
1257 35010 : if (dimIdx == nDims)
1258 : {
1259 28646 : size_t dimIdxSubLoop = 0;
1260 28646 : dstPtrStackInnerLoop[0] = dstPtrStackOuterLoop[nDims];
1261 28646 : bool bEmptyChunk = false;
1262 :
1263 28646 : const GByte *pabySrcBlock = m_abyDecodedBlockData.empty()
1264 27655 : ? m_abyRawBlockData.data()
1265 28646 : : m_abyDecodedBlockData.data();
1266 28646 : bool bMatchFoundInMapChunkIndexToCachedBlock = false;
1267 :
1268 : // Use cache built by IAdviseRead() if possible
1269 28646 : if (!m_oChunkCache.empty())
1270 : {
1271 24480 : const auto oIter = m_oChunkCache.find(blockIndices);
1272 24480 : if (oIter != m_oChunkCache.end())
1273 : {
1274 23366 : bMatchFoundInMapChunkIndexToCachedBlock = true;
1275 23366 : if (oIter->second.abyDecoded.empty())
1276 : {
1277 4 : bEmptyChunk = true;
1278 : }
1279 : else
1280 : {
1281 23362 : pabySrcBlock = oIter->second.abyDecoded.data();
1282 : }
1283 : }
1284 : else
1285 : {
1286 : #ifdef DEBUG
1287 2228 : std::string key;
1288 3342 : for (size_t j = 0; j < nDims; ++j)
1289 : {
1290 2228 : if (j)
1291 1114 : key += ',';
1292 2228 : key += std::to_string(blockIndices[j]);
1293 : }
1294 1114 : CPLDebugOnly(ZARR_DEBUG_KEY, "Cache miss for tile %s",
1295 : key.c_str());
1296 : #endif
1297 : }
1298 : }
1299 :
1300 28646 : if (!bMatchFoundInMapChunkIndexToCachedBlock)
1301 : {
1302 5280 : if (!blockIndices.empty() && blockIndices == m_anCachedBlockIndices)
1303 : {
1304 224 : if (!m_bCachedBlockValid)
1305 575 : return false;
1306 224 : bEmptyChunk = m_bCachedBlockEmpty;
1307 : }
1308 : else
1309 : {
1310 5056 : if (!FlushDirtyBlock())
1311 0 : return false;
1312 :
1313 5056 : m_anCachedBlockIndices = blockIndices;
1314 5056 : m_bCachedBlockValid =
1315 5056 : LoadBlockData(blockIndices.data(), bEmptyChunk);
1316 5056 : if (!m_bCachedBlockValid)
1317 : {
1318 575 : return false;
1319 : }
1320 4481 : m_bCachedBlockEmpty = bEmptyChunk;
1321 : }
1322 :
1323 9410 : pabySrcBlock = m_abyDecodedBlockData.empty()
1324 3714 : ? m_abyRawBlockData.data()
1325 991 : : m_abyDecodedBlockData.data();
1326 : }
1327 : const size_t nSrcDTSize =
1328 28071 : m_abyDecodedBlockData.empty() ? nSourceSize : nDTSize;
1329 :
1330 84330 : for (size_t i = 0; i < nDims; ++i)
1331 : {
1332 56259 : countInnerLoopInit[i] = 1;
1333 56259 : if (arrayStep[i] != 0)
1334 : {
1335 : const auto nextBlockIdx =
1336 55557 : std::min((1 + indicesOuterLoop[i] / m_anInnerBlockSize[i]) *
1337 55557 : m_anInnerBlockSize[i],
1338 111114 : arrayStartIdx[i] + count[i] * arrayStep[i]);
1339 55557 : countInnerLoopInit[i] = static_cast<size_t>(cpl::div_round_up(
1340 55557 : nextBlockIdx - indicesOuterLoop[i], arrayStep[i]));
1341 : }
1342 : }
1343 :
1344 28071 : if (bEmptyChunk && bBothAreNumericDT && abyTargetNoData.empty())
1345 : {
1346 605 : abyTargetNoData.resize(nBufferDTSize);
1347 605 : if (m_pabyNoData)
1348 : {
1349 157 : GDALExtendedDataType::CopyValue(
1350 157 : m_pabyNoData, m_oType, &abyTargetNoData[0], bufferDataType);
1351 157 : bNoDataIsZero = true;
1352 1311 : for (size_t i = 0; i < abyTargetNoData.size(); ++i)
1353 : {
1354 1154 : if (abyTargetNoData[i] != 0)
1355 335 : bNoDataIsZero = false;
1356 : }
1357 : }
1358 : else
1359 : {
1360 448 : bNoDataIsZero = true;
1361 448 : GByte zero = 0;
1362 448 : GDALCopyWords(&zero, GDT_UInt8, 0, &abyTargetNoData[0],
1363 : bufferDataType.GetNumericDataType(), 0, 1);
1364 : }
1365 : }
1366 :
1367 27466 : lbl_next_depth_inner_loop:
1368 483718 : if (nDims == 0 || dimIdxSubLoop == nDims - 1)
1369 : {
1370 455526 : indicesInnerLoop[dimIdxSubLoop] = indicesOuterLoop[dimIdxSubLoop];
1371 455526 : void *dst_ptr = dstPtrStackInnerLoop[dimIdxSubLoop];
1372 :
1373 453039 : if (m_bUseOptimizedCodePaths && bEmptyChunk && bBothAreNumericDT &&
1374 908565 : bNoDataIsZero &&
1375 1138 : nBufferDTSize == dstBufferStrideBytes[dimIdxSubLoop])
1376 : {
1377 1072 : memset(dst_ptr, 0,
1378 1072 : nBufferDTSize * countInnerLoopInit[dimIdxSubLoop]);
1379 1072 : goto end_inner_loop;
1380 : }
1381 451967 : else if (m_bUseOptimizedCodePaths && bEmptyChunk &&
1382 906802 : !abyTargetNoData.empty() && bBothAreNumericDT &&
1383 381 : dstBufferStrideBytes[dimIdxSubLoop] <
1384 381 : std::numeric_limits<int>::max())
1385 : {
1386 762 : GDALCopyWords64(
1387 381 : abyTargetNoData.data(), bufferDataType.GetNumericDataType(),
1388 : 0, dst_ptr, bufferDataType.GetNumericDataType(),
1389 381 : static_cast<int>(dstBufferStrideBytes[dimIdxSubLoop]),
1390 381 : static_cast<GPtrDiff_t>(countInnerLoopInit[dimIdxSubLoop]));
1391 381 : goto end_inner_loop;
1392 : }
1393 454073 : else if (bEmptyChunk)
1394 : {
1395 3875 : for (size_t i = 0; i < countInnerLoopInit[dimIdxSubLoop];
1396 2528 : ++i, dst_ptr = static_cast<uint8_t *>(dst_ptr) +
1397 2528 : dstBufferStrideBytes[dimIdxSubLoop])
1398 : {
1399 2528 : if (bNoDataIsZero)
1400 : {
1401 1930 : if (nBufferDTSize == 1)
1402 : {
1403 36 : *static_cast<uint8_t *>(dst_ptr) = 0;
1404 : }
1405 1894 : else if (nBufferDTSize == 2)
1406 : {
1407 48 : *static_cast<uint16_t *>(dst_ptr) = 0;
1408 : }
1409 1846 : else if (nBufferDTSize == 4)
1410 : {
1411 72 : *static_cast<uint32_t *>(dst_ptr) = 0;
1412 : }
1413 1774 : else if (nBufferDTSize == 8)
1414 : {
1415 1534 : *static_cast<uint64_t *>(dst_ptr) = 0;
1416 : }
1417 240 : else if (nBufferDTSize == 16)
1418 : {
1419 240 : static_cast<uint64_t *>(dst_ptr)[0] = 0;
1420 240 : static_cast<uint64_t *>(dst_ptr)[1] = 0;
1421 : }
1422 : else
1423 : {
1424 0 : CPLAssert(false);
1425 : }
1426 : }
1427 598 : else if (m_pabyNoData)
1428 : {
1429 597 : if (bBothAreNumericDT)
1430 : {
1431 588 : const void *src_ptr_v = abyTargetNoData.data();
1432 588 : if (nBufferDTSize == 1)
1433 24 : *static_cast<uint8_t *>(dst_ptr) =
1434 24 : *static_cast<const uint8_t *>(src_ptr_v);
1435 564 : else if (nBufferDTSize == 2)
1436 0 : *static_cast<uint16_t *>(dst_ptr) =
1437 0 : *static_cast<const uint16_t *>(src_ptr_v);
1438 564 : else if (nBufferDTSize == 4)
1439 60 : *static_cast<uint32_t *>(dst_ptr) =
1440 60 : *static_cast<const uint32_t *>(src_ptr_v);
1441 504 : else if (nBufferDTSize == 8)
1442 504 : *static_cast<uint64_t *>(dst_ptr) =
1443 504 : *static_cast<const uint64_t *>(src_ptr_v);
1444 0 : else if (nBufferDTSize == 16)
1445 : {
1446 0 : static_cast<uint64_t *>(dst_ptr)[0] =
1447 0 : static_cast<const uint64_t *>(src_ptr_v)[0];
1448 0 : static_cast<uint64_t *>(dst_ptr)[1] =
1449 0 : static_cast<const uint64_t *>(src_ptr_v)[1];
1450 : }
1451 : else
1452 : {
1453 0 : CPLAssert(false);
1454 : }
1455 : }
1456 : else
1457 : {
1458 9 : GDALExtendedDataType::CopyValue(
1459 9 : m_pabyNoData, m_oType, dst_ptr, bufferDataType);
1460 : }
1461 : }
1462 : else
1463 : {
1464 1 : memset(dst_ptr, 0, nBufferDTSize);
1465 : }
1466 : }
1467 :
1468 1347 : goto end_inner_loop;
1469 : }
1470 :
1471 452726 : size_t nOffset = 0;
1472 1370610 : for (size_t i = 0; i < nDims; i++)
1473 : {
1474 917882 : nOffset = static_cast<size_t>(
1475 917882 : nOffset * m_anInnerBlockSize[i] +
1476 917882 : (indicesInnerLoop[i] -
1477 917882 : blockIndices[i] * m_anInnerBlockSize[i]));
1478 : }
1479 452726 : const GByte *src_ptr = pabySrcBlock + nOffset * nSrcDTSize;
1480 452726 : const auto step = nDims == 0 ? 0 : arrayStep[dimIdxSubLoop];
1481 :
1482 451576 : if (m_bUseOptimizedCodePaths && bBothAreNumericDT &&
1483 451551 : step <= static_cast<GIntBig>(std::numeric_limits<int>::max() /
1484 904302 : nDTSize) &&
1485 451551 : dstBufferStrideBytes[dimIdxSubLoop] <=
1486 451551 : std::numeric_limits<int>::max())
1487 : {
1488 451551 : GDALCopyWords64(
1489 : src_ptr, m_oType.GetNumericDataType(),
1490 : static_cast<int>(step * nDTSize), dst_ptr,
1491 : bufferDataType.GetNumericDataType(),
1492 451551 : static_cast<int>(dstBufferStrideBytes[dimIdxSubLoop]),
1493 451551 : static_cast<GPtrDiff_t>(countInnerLoopInit[dimIdxSubLoop]));
1494 :
1495 451551 : goto end_inner_loop;
1496 : }
1497 :
1498 3393 : for (size_t i = 0; i < countInnerLoopInit[dimIdxSubLoop];
1499 2218 : ++i, src_ptr += step * nSrcDTSize,
1500 2218 : dst_ptr = static_cast<uint8_t *>(dst_ptr) +
1501 2218 : dstBufferStrideBytes[dimIdxSubLoop])
1502 : {
1503 2218 : if (bSameNumericDT)
1504 : {
1505 814 : const void *src_ptr_v = src_ptr;
1506 814 : if (nSameDTSize == 1)
1507 70 : *static_cast<uint8_t *>(dst_ptr) =
1508 70 : *static_cast<const uint8_t *>(src_ptr_v);
1509 744 : else if (nSameDTSize == 2)
1510 : {
1511 56 : *static_cast<uint16_t *>(dst_ptr) =
1512 56 : *static_cast<const uint16_t *>(src_ptr_v);
1513 : }
1514 688 : else if (nSameDTSize == 4)
1515 : {
1516 154 : *static_cast<uint32_t *>(dst_ptr) =
1517 154 : *static_cast<const uint32_t *>(src_ptr_v);
1518 : }
1519 534 : else if (nSameDTSize == 8)
1520 : {
1521 482 : *static_cast<uint64_t *>(dst_ptr) =
1522 482 : *static_cast<const uint64_t *>(src_ptr_v);
1523 : }
1524 52 : else if (nSameDTSize == 16)
1525 : {
1526 52 : static_cast<uint64_t *>(dst_ptr)[0] =
1527 52 : static_cast<const uint64_t *>(src_ptr_v)[0];
1528 52 : static_cast<uint64_t *>(dst_ptr)[1] =
1529 52 : static_cast<const uint64_t *>(src_ptr_v)[1];
1530 : }
1531 : else
1532 : {
1533 0 : CPLAssert(false);
1534 : }
1535 : }
1536 1404 : else if (bSameCompoundAndNoDynamicMem)
1537 : {
1538 4 : memcpy(dst_ptr, src_ptr, nDTSize);
1539 : }
1540 1400 : else if (m_oType.GetClass() == GEDTC_STRING)
1541 : {
1542 27 : if (m_aoDtypeElts.back().nativeType ==
1543 : DtypeElt::NativeType::STRING_UNICODE)
1544 : {
1545 : char *pDstStr =
1546 14 : UCS4ToUTF8(src_ptr, nSourceSize,
1547 7 : m_aoDtypeElts.back().needByteSwapping);
1548 7 : char **pDstPtr = static_cast<char **>(dst_ptr);
1549 7 : memcpy(pDstPtr, &pDstStr, sizeof(pDstStr));
1550 : }
1551 : else
1552 : {
1553 : char *pDstStr =
1554 20 : static_cast<char *>(CPLMalloc(nSourceSize + 1));
1555 20 : memcpy(pDstStr, src_ptr, nSourceSize);
1556 20 : pDstStr[nSourceSize] = 0;
1557 20 : char **pDstPtr = static_cast<char **>(dst_ptr);
1558 20 : memcpy(pDstPtr, &pDstStr, sizeof(char *));
1559 : }
1560 : }
1561 : else
1562 : {
1563 1373 : GDALExtendedDataType::CopyValue(src_ptr, m_oType, dst_ptr,
1564 : bufferDataType);
1565 : }
1566 1175 : }
1567 : }
1568 : else
1569 : {
1570 : // This level of loop loops over individual samples, within a
1571 : // block
1572 28192 : indicesInnerLoop[dimIdxSubLoop] = indicesOuterLoop[dimIdxSubLoop];
1573 28192 : countInnerLoop[dimIdxSubLoop] = countInnerLoopInit[dimIdxSubLoop];
1574 : while (true)
1575 : {
1576 455647 : dimIdxSubLoop++;
1577 455647 : dstPtrStackInnerLoop[dimIdxSubLoop] =
1578 455647 : dstPtrStackInnerLoop[dimIdxSubLoop - 1];
1579 455647 : goto lbl_next_depth_inner_loop;
1580 455647 : lbl_return_to_caller_inner_loop:
1581 455647 : dimIdxSubLoop--;
1582 455647 : --countInnerLoop[dimIdxSubLoop];
1583 455647 : if (countInnerLoop[dimIdxSubLoop] == 0)
1584 : {
1585 28192 : break;
1586 : }
1587 427455 : indicesInnerLoop[dimIdxSubLoop] += arrayStep[dimIdxSubLoop];
1588 427455 : dstPtrStackInnerLoop[dimIdxSubLoop] +=
1589 427455 : dstBufferStrideBytes[dimIdxSubLoop];
1590 : }
1591 : }
1592 483718 : end_inner_loop:
1593 483718 : if (dimIdxSubLoop > 0)
1594 455647 : goto lbl_return_to_caller_inner_loop;
1595 : }
1596 : else
1597 : {
1598 : // This level of loop loops over blocks
1599 6364 : indicesOuterLoop[dimIdx] = arrayStartIdx[dimIdx];
1600 6364 : blockIndices[dimIdx] =
1601 6364 : indicesOuterLoop[dimIdx] / m_anInnerBlockSize[dimIdx];
1602 : while (true)
1603 : {
1604 32728 : dimIdx++;
1605 32728 : dstPtrStackOuterLoop[dimIdx] = dstPtrStackOuterLoop[dimIdx - 1];
1606 32728 : goto lbl_next_depth;
1607 31582 : lbl_return_to_caller:
1608 31582 : dimIdx--;
1609 31582 : if (count[dimIdx] == 1 || arrayStep[dimIdx] == 0)
1610 : break;
1611 :
1612 : size_t nIncr;
1613 30488 : if (static_cast<GUInt64>(arrayStep[dimIdx]) <
1614 30488 : m_anInnerBlockSize[dimIdx])
1615 : {
1616 : // Compute index at next block boundary
1617 : auto newIdx =
1618 30075 : indicesOuterLoop[dimIdx] +
1619 30075 : (m_anInnerBlockSize[dimIdx] -
1620 30075 : (indicesOuterLoop[dimIdx] % m_anInnerBlockSize[dimIdx]));
1621 : // And round up compared to arrayStartIdx, arrayStep
1622 30075 : nIncr = static_cast<size_t>(cpl::div_round_up(
1623 30075 : newIdx - indicesOuterLoop[dimIdx], arrayStep[dimIdx]));
1624 : }
1625 : else
1626 : {
1627 413 : nIncr = 1;
1628 : }
1629 30488 : indicesOuterLoop[dimIdx] += nIncr * arrayStep[dimIdx];
1630 30488 : if (indicesOuterLoop[dimIdx] >
1631 30488 : arrayStartIdx[dimIdx] + (count[dimIdx] - 1) * arrayStep[dimIdx])
1632 4124 : break;
1633 26364 : dstPtrStackOuterLoop[dimIdx] +=
1634 26364 : bufferStride[dimIdx] *
1635 26364 : static_cast<GPtrDiff_t>(nIncr * nBufferDTSize);
1636 52728 : blockIndices[dimIdx] =
1637 26364 : indicesOuterLoop[dimIdx] / m_anInnerBlockSize[dimIdx];
1638 26364 : }
1639 : }
1640 33289 : if (dimIdx > 0)
1641 31582 : goto lbl_return_to_caller;
1642 :
1643 1707 : return true;
1644 : }
1645 :
1646 : /************************************************************************/
1647 : /* ZarrArray::IWrite() */
1648 : /************************************************************************/
1649 :
1650 498 : bool ZarrArray::IWrite(const GUInt64 *arrayStartIdx, const size_t *count,
1651 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
1652 : const GDALExtendedDataType &bufferDataType,
1653 : const void *pSrcBuffer)
1654 : {
1655 498 : if (!CheckValidAndErrorOutIfNot())
1656 0 : return false;
1657 :
1658 498 : if (!AllocateWorkingBuffers())
1659 0 : return false;
1660 :
1661 498 : m_oChunkCache.clear();
1662 :
1663 : // Need to be kept in top-level scope
1664 996 : std::vector<GUInt64> arrayStartIdxMod;
1665 996 : std::vector<GInt64> arrayStepMod;
1666 996 : std::vector<GPtrDiff_t> bufferStrideMod;
1667 :
1668 498 : const size_t nDims = m_aoDims.size();
1669 498 : bool negativeStep = false;
1670 498 : bool bWriteWholeBlockInit = true;
1671 1394 : for (size_t i = 0; i < nDims; ++i)
1672 : {
1673 896 : if (arrayStep[i] < 0)
1674 : {
1675 132 : negativeStep = true;
1676 132 : if (arrayStep[i] != -1 && count[i] > 1)
1677 0 : bWriteWholeBlockInit = false;
1678 : }
1679 764 : else if (arrayStep[i] != 1 && count[i] > 1)
1680 0 : bWriteWholeBlockInit = false;
1681 : }
1682 :
1683 498 : const auto nBufferDTSize = static_cast<int>(bufferDataType.GetSize());
1684 :
1685 : // Make sure that arrayStep[i] are positive for sake of simplicity
1686 498 : if (negativeStep)
1687 : {
1688 : #if defined(__GNUC__)
1689 : #pragma GCC diagnostic push
1690 : #pragma GCC diagnostic ignored "-Wnull-dereference"
1691 : #endif
1692 66 : arrayStartIdxMod.resize(nDims);
1693 66 : arrayStepMod.resize(nDims);
1694 66 : bufferStrideMod.resize(nDims);
1695 : #if defined(__GNUC__)
1696 : #pragma GCC diagnostic pop
1697 : #endif
1698 198 : for (size_t i = 0; i < nDims; ++i)
1699 : {
1700 132 : if (arrayStep[i] < 0)
1701 : {
1702 264 : arrayStartIdxMod[i] =
1703 132 : arrayStartIdx[i] - (count[i] - 1) * (-arrayStep[i]);
1704 132 : arrayStepMod[i] = -arrayStep[i];
1705 132 : bufferStrideMod[i] = -bufferStride[i];
1706 132 : pSrcBuffer =
1707 : static_cast<const GByte *>(pSrcBuffer) +
1708 132 : bufferStride[i] *
1709 132 : static_cast<GPtrDiff_t>(nBufferDTSize * (count[i] - 1));
1710 : }
1711 : else
1712 : {
1713 0 : arrayStartIdxMod[i] = arrayStartIdx[i];
1714 0 : arrayStepMod[i] = arrayStep[i];
1715 0 : bufferStrideMod[i] = bufferStride[i];
1716 : }
1717 : }
1718 66 : arrayStartIdx = arrayStartIdxMod.data();
1719 66 : arrayStep = arrayStepMod.data();
1720 66 : bufferStride = bufferStrideMod.data();
1721 : }
1722 :
1723 996 : std::vector<uint64_t> indicesOuterLoop(nDims + 1);
1724 996 : std::vector<const GByte *> srcPtrStackOuterLoop(nDims + 1);
1725 :
1726 996 : std::vector<size_t> offsetDstBuffer(nDims + 1);
1727 996 : std::vector<const GByte *> srcPtrStackInnerLoop(nDims + 1);
1728 :
1729 996 : std::vector<GPtrDiff_t> srcBufferStrideBytes;
1730 1394 : for (size_t i = 0; i < nDims; ++i)
1731 : {
1732 896 : srcBufferStrideBytes.push_back(bufferStride[i] *
1733 896 : static_cast<GPtrDiff_t>(nBufferDTSize));
1734 : }
1735 498 : srcBufferStrideBytes.push_back(0);
1736 :
1737 498 : const auto nDTSize = m_oType.GetSize();
1738 :
1739 996 : std::vector<uint64_t> blockIndices(nDims);
1740 : const size_t nNativeSize =
1741 498 : m_aoDtypeElts.back().nativeOffset + m_aoDtypeElts.back().nativeSize;
1742 :
1743 996 : std::vector<size_t> countInnerLoopInit(nDims + 1, 1);
1744 996 : std::vector<size_t> countInnerLoop(nDims);
1745 :
1746 992 : const bool bBothAreNumericDT = m_oType.GetClass() == GEDTC_NUMERIC &&
1747 494 : bufferDataType.GetClass() == GEDTC_NUMERIC;
1748 : const bool bSameNumericDT =
1749 992 : bBothAreNumericDT &&
1750 494 : m_oType.GetNumericDataType() == bufferDataType.GetNumericDataType();
1751 498 : const auto nSameDTSize = bSameNumericDT ? m_oType.GetSize() : 0;
1752 : const bool bSameCompoundAndNoDynamicMem =
1753 498 : m_oType.GetClass() == GEDTC_COMPOUND && m_oType == bufferDataType &&
1754 0 : !m_oType.NeedsFreeDynamicMemory();
1755 :
1756 498 : size_t dimIdx = 0;
1757 498 : size_t dimIdxForCopy = nDims == 0 ? 0 : nDims - 1;
1758 498 : if (nDims)
1759 : {
1760 566 : while (dimIdxForCopy > 0 && count[dimIdxForCopy] == 1)
1761 68 : --dimIdxForCopy;
1762 : }
1763 :
1764 498 : srcPtrStackOuterLoop[0] = static_cast<const GByte *>(pSrcBuffer);
1765 24654 : lbl_next_depth:
1766 24654 : if (dimIdx == nDims)
1767 : {
1768 22824 : bool bWriteWholeBlock = bWriteWholeBlockInit;
1769 22824 : bool bPartialBlock = false;
1770 68633 : for (size_t i = 0; i < nDims; ++i)
1771 : {
1772 45809 : countInnerLoopInit[i] = 1;
1773 45809 : if (arrayStep[i] != 0)
1774 : {
1775 : const auto nextBlockIdx =
1776 45526 : std::min((1 + indicesOuterLoop[i] / m_anOuterBlockSize[i]) *
1777 45526 : m_anOuterBlockSize[i],
1778 91052 : arrayStartIdx[i] + count[i] * arrayStep[i]);
1779 45526 : countInnerLoopInit[i] = static_cast<size_t>(cpl::div_round_up(
1780 45526 : nextBlockIdx - indicesOuterLoop[i], arrayStep[i]));
1781 : }
1782 45809 : if (bWriteWholeBlock)
1783 : {
1784 : const bool bWholePartialBlockThisDim =
1785 46701 : indicesOuterLoop[i] == 0 &&
1786 1840 : countInnerLoopInit[i] == m_aoDims[i]->GetSize();
1787 44861 : bWriteWholeBlock =
1788 44861 : (countInnerLoopInit[i] == m_anOuterBlockSize[i] ||
1789 : bWholePartialBlockThisDim);
1790 44861 : if (bWholePartialBlockThisDim)
1791 : {
1792 593 : bPartialBlock = true;
1793 : }
1794 : }
1795 : }
1796 :
1797 22824 : size_t dimIdxSubLoop = 0;
1798 22824 : srcPtrStackInnerLoop[0] = srcPtrStackOuterLoop[nDims];
1799 : const size_t nCacheDTSize =
1800 22824 : m_abyDecodedBlockData.empty() ? nNativeSize : nDTSize;
1801 22824 : auto &abyBlock = m_abyDecodedBlockData.empty() ? m_abyRawBlockData
1802 22824 : : m_abyDecodedBlockData;
1803 :
1804 22824 : if (!blockIndices.empty() && blockIndices == m_anCachedBlockIndices)
1805 : {
1806 4 : if (!m_bCachedBlockValid)
1807 0 : return false;
1808 : }
1809 : else
1810 : {
1811 22820 : if (!FlushDirtyBlock())
1812 0 : return false;
1813 :
1814 22820 : m_anCachedBlockIndices = blockIndices;
1815 22820 : m_bCachedBlockValid = true;
1816 :
1817 22820 : if (bWriteWholeBlock)
1818 : {
1819 21296 : if (bPartialBlock)
1820 : {
1821 292 : DeallocateDecodedBlockData();
1822 292 : memset(&abyBlock[0], 0, abyBlock.size());
1823 : }
1824 : }
1825 : else
1826 : {
1827 : // If we don't write the whole tile, we need to fetch a
1828 : // potentially existing one.
1829 1524 : bool bEmptyBlock = false;
1830 1524 : m_bCachedBlockValid =
1831 1524 : LoadBlockData(blockIndices.data(), bEmptyBlock);
1832 1524 : if (!m_bCachedBlockValid)
1833 : {
1834 0 : return false;
1835 : }
1836 :
1837 1524 : if (bEmptyBlock)
1838 : {
1839 1215 : DeallocateDecodedBlockData();
1840 :
1841 1215 : if (m_pabyNoData == nullptr)
1842 : {
1843 495 : memset(&abyBlock[0], 0, abyBlock.size());
1844 : }
1845 : else
1846 : {
1847 720 : const size_t nElts = abyBlock.size() / nCacheDTSize;
1848 720 : GByte *dstPtr = &abyBlock[0];
1849 720 : if (m_oType.GetClass() == GEDTC_NUMERIC)
1850 : {
1851 720 : GDALCopyWords64(
1852 720 : m_pabyNoData, m_oType.GetNumericDataType(), 0,
1853 : dstPtr, m_oType.GetNumericDataType(),
1854 720 : static_cast<int>(m_oType.GetSize()),
1855 : static_cast<GPtrDiff_t>(nElts));
1856 : }
1857 : else
1858 : {
1859 0 : for (size_t i = 0; i < nElts; ++i)
1860 : {
1861 0 : GDALExtendedDataType::CopyValue(
1862 0 : m_pabyNoData, m_oType, dstPtr, m_oType);
1863 0 : dstPtr += nCacheDTSize;
1864 : }
1865 : }
1866 : }
1867 : }
1868 : }
1869 : }
1870 22824 : m_bDirtyBlock = true;
1871 22824 : m_bCachedBlockEmpty = false;
1872 22824 : if (nDims)
1873 22824 : offsetDstBuffer[0] = static_cast<size_t>(
1874 22824 : indicesOuterLoop[0] - blockIndices[0] * m_anOuterBlockSize[0]);
1875 :
1876 22824 : GByte *pabyBlock = &abyBlock[0];
1877 :
1878 453451 : lbl_next_depth_inner_loop:
1879 453451 : if (dimIdxSubLoop == dimIdxForCopy)
1880 : {
1881 430691 : size_t nOffset = offsetDstBuffer[dimIdxSubLoop];
1882 430691 : GInt64 step = nDims == 0 ? 0 : arrayStep[dimIdxSubLoop];
1883 430916 : for (size_t i = dimIdxSubLoop + 1; i < nDims; ++i)
1884 : {
1885 225 : nOffset = static_cast<size_t>(
1886 225 : nOffset * m_anOuterBlockSize[i] +
1887 225 : (indicesOuterLoop[i] -
1888 225 : blockIndices[i] * m_anOuterBlockSize[i]));
1889 225 : step *= m_anOuterBlockSize[i];
1890 : }
1891 430691 : const void *src_ptr = srcPtrStackInnerLoop[dimIdxSubLoop];
1892 430691 : GByte *dst_ptr = pabyBlock + nOffset * nCacheDTSize;
1893 :
1894 430691 : if (m_bUseOptimizedCodePaths && bBothAreNumericDT)
1895 : {
1896 429961 : if (countInnerLoopInit[dimIdxSubLoop] == 1 && bSameNumericDT)
1897 : {
1898 86 : void *dst_ptr_v = dst_ptr;
1899 86 : if (nSameDTSize == 1)
1900 8 : *static_cast<uint8_t *>(dst_ptr_v) =
1901 8 : *static_cast<const uint8_t *>(src_ptr);
1902 78 : else if (nSameDTSize == 2)
1903 : {
1904 2 : *static_cast<uint16_t *>(dst_ptr_v) =
1905 2 : *static_cast<const uint16_t *>(src_ptr);
1906 : }
1907 76 : else if (nSameDTSize == 4)
1908 : {
1909 2 : *static_cast<uint32_t *>(dst_ptr_v) =
1910 2 : *static_cast<const uint32_t *>(src_ptr);
1911 : }
1912 74 : else if (nSameDTSize == 8)
1913 : {
1914 52 : *static_cast<uint64_t *>(dst_ptr_v) =
1915 52 : *static_cast<const uint64_t *>(src_ptr);
1916 : }
1917 22 : else if (nSameDTSize == 16)
1918 : {
1919 22 : static_cast<uint64_t *>(dst_ptr_v)[0] =
1920 22 : static_cast<const uint64_t *>(src_ptr)[0];
1921 22 : static_cast<uint64_t *>(dst_ptr_v)[1] =
1922 22 : static_cast<const uint64_t *>(src_ptr)[1];
1923 : }
1924 : else
1925 : {
1926 0 : CPLAssert(false);
1927 : }
1928 : }
1929 429875 : else if (step <=
1930 429875 : static_cast<GIntBig>(
1931 859750 : std::numeric_limits<int>::max() / nDTSize) &&
1932 429875 : srcBufferStrideBytes[dimIdxSubLoop] <=
1933 429875 : std::numeric_limits<int>::max())
1934 : {
1935 859750 : GDALCopyWords64(
1936 : src_ptr, bufferDataType.GetNumericDataType(),
1937 429875 : static_cast<int>(srcBufferStrideBytes[dimIdxSubLoop]),
1938 : dst_ptr, m_oType.GetNumericDataType(),
1939 : static_cast<int>(step * nDTSize),
1940 : static_cast<GPtrDiff_t>(
1941 429875 : countInnerLoopInit[dimIdxSubLoop]));
1942 : }
1943 429961 : goto end_inner_loop;
1944 : }
1945 :
1946 2188 : for (size_t i = 0; i < countInnerLoopInit[dimIdxSubLoop];
1947 1458 : ++i, dst_ptr += step * nCacheDTSize,
1948 1458 : src_ptr = static_cast<const uint8_t *>(src_ptr) +
1949 1458 : srcBufferStrideBytes[dimIdxSubLoop])
1950 : {
1951 1458 : if (bSameNumericDT)
1952 : {
1953 300 : void *dst_ptr_v = dst_ptr;
1954 300 : if (nSameDTSize == 1)
1955 0 : *static_cast<uint8_t *>(dst_ptr_v) =
1956 0 : *static_cast<const uint8_t *>(src_ptr);
1957 300 : else if (nSameDTSize == 2)
1958 : {
1959 0 : *static_cast<uint16_t *>(dst_ptr_v) =
1960 0 : *static_cast<const uint16_t *>(src_ptr);
1961 : }
1962 300 : else if (nSameDTSize == 4)
1963 : {
1964 0 : *static_cast<uint32_t *>(dst_ptr_v) =
1965 0 : *static_cast<const uint32_t *>(src_ptr);
1966 : }
1967 300 : else if (nSameDTSize == 8)
1968 : {
1969 220 : *static_cast<uint64_t *>(dst_ptr_v) =
1970 220 : *static_cast<const uint64_t *>(src_ptr);
1971 : }
1972 80 : else if (nSameDTSize == 16)
1973 : {
1974 80 : static_cast<uint64_t *>(dst_ptr_v)[0] =
1975 80 : static_cast<const uint64_t *>(src_ptr)[0];
1976 80 : static_cast<uint64_t *>(dst_ptr_v)[1] =
1977 80 : static_cast<const uint64_t *>(src_ptr)[1];
1978 : }
1979 : else
1980 : {
1981 0 : CPLAssert(false);
1982 : }
1983 : }
1984 1158 : else if (bSameCompoundAndNoDynamicMem)
1985 : {
1986 0 : memcpy(dst_ptr, src_ptr, nDTSize);
1987 : }
1988 1158 : else if (m_oType.GetClass() == GEDTC_STRING)
1989 : {
1990 6 : const char *pSrcStr =
1991 : *static_cast<const char *const *>(src_ptr);
1992 6 : if (pSrcStr)
1993 : {
1994 6 : const size_t nLen = strlen(pSrcStr);
1995 6 : if (m_aoDtypeElts.back().nativeType ==
1996 : DtypeElt::NativeType::STRING_UNICODE)
1997 : {
1998 : try
1999 : {
2000 : const auto ucs4 = UTF8ToUCS4(
2001 : pSrcStr,
2002 8 : m_aoDtypeElts.back().needByteSwapping);
2003 4 : const auto ucs4Len = ucs4.size();
2004 4 : memcpy(dst_ptr, ucs4.data(),
2005 4 : std::min(ucs4Len, nNativeSize));
2006 4 : if (ucs4Len > nNativeSize)
2007 : {
2008 1 : CPLError(CE_Warning, CPLE_AppDefined,
2009 : "Too long string truncated");
2010 : }
2011 3 : else if (ucs4Len < nNativeSize)
2012 : {
2013 1 : memset(dst_ptr + ucs4Len, 0,
2014 1 : nNativeSize - ucs4Len);
2015 : }
2016 : }
2017 0 : catch (const std::exception &)
2018 : {
2019 0 : memset(dst_ptr, 0, nNativeSize);
2020 : }
2021 : }
2022 : else
2023 : {
2024 2 : memcpy(dst_ptr, pSrcStr,
2025 2 : std::min(nLen, nNativeSize));
2026 2 : if (nLen < nNativeSize)
2027 1 : memset(dst_ptr + nLen, 0, nNativeSize - nLen);
2028 : }
2029 : }
2030 : else
2031 : {
2032 0 : memset(dst_ptr, 0, nNativeSize);
2033 : }
2034 : }
2035 : else
2036 : {
2037 1152 : if (m_oType.NeedsFreeDynamicMemory())
2038 0 : m_oType.FreeDynamicMemory(dst_ptr);
2039 1152 : GDALExtendedDataType::CopyValue(src_ptr, bufferDataType,
2040 1152 : dst_ptr, m_oType);
2041 : }
2042 : }
2043 : }
2044 : else
2045 : {
2046 : // This level of loop loops over individual samples, within a
2047 : // block
2048 22760 : countInnerLoop[dimIdxSubLoop] = countInnerLoopInit[dimIdxSubLoop];
2049 : while (true)
2050 : {
2051 430627 : dimIdxSubLoop++;
2052 430627 : srcPtrStackInnerLoop[dimIdxSubLoop] =
2053 430627 : srcPtrStackInnerLoop[dimIdxSubLoop - 1];
2054 861254 : offsetDstBuffer[dimIdxSubLoop] = static_cast<size_t>(
2055 1291880 : offsetDstBuffer[dimIdxSubLoop - 1] *
2056 430627 : m_anOuterBlockSize[dimIdxSubLoop] +
2057 430627 : (indicesOuterLoop[dimIdxSubLoop] -
2058 861254 : blockIndices[dimIdxSubLoop] *
2059 430627 : m_anOuterBlockSize[dimIdxSubLoop]));
2060 430627 : goto lbl_next_depth_inner_loop;
2061 430627 : lbl_return_to_caller_inner_loop:
2062 430627 : dimIdxSubLoop--;
2063 430627 : --countInnerLoop[dimIdxSubLoop];
2064 430627 : if (countInnerLoop[dimIdxSubLoop] == 0)
2065 : {
2066 22760 : break;
2067 : }
2068 407867 : srcPtrStackInnerLoop[dimIdxSubLoop] +=
2069 407867 : srcBufferStrideBytes[dimIdxSubLoop];
2070 407867 : offsetDstBuffer[dimIdxSubLoop] +=
2071 407867 : static_cast<size_t>(arrayStep[dimIdxSubLoop]);
2072 : }
2073 : }
2074 430691 : end_inner_loop:
2075 453451 : if (dimIdxSubLoop > 0)
2076 430627 : goto lbl_return_to_caller_inner_loop;
2077 : }
2078 : else
2079 : {
2080 : // This level of loop loops over blocks
2081 1830 : indicesOuterLoop[dimIdx] = arrayStartIdx[dimIdx];
2082 1830 : blockIndices[dimIdx] =
2083 1830 : indicesOuterLoop[dimIdx] / m_anOuterBlockSize[dimIdx];
2084 : while (true)
2085 : {
2086 24156 : dimIdx++;
2087 24156 : srcPtrStackOuterLoop[dimIdx] = srcPtrStackOuterLoop[dimIdx - 1];
2088 24156 : goto lbl_next_depth;
2089 24156 : lbl_return_to_caller:
2090 24156 : dimIdx--;
2091 24156 : if (count[dimIdx] == 1 || arrayStep[dimIdx] == 0)
2092 : break;
2093 :
2094 : size_t nIncr;
2095 23903 : if (static_cast<GUInt64>(arrayStep[dimIdx]) <
2096 23903 : m_anOuterBlockSize[dimIdx])
2097 : {
2098 : // Compute index at next block boundary
2099 : auto newIdx =
2100 23711 : indicesOuterLoop[dimIdx] +
2101 23711 : (m_anOuterBlockSize[dimIdx] -
2102 23711 : (indicesOuterLoop[dimIdx] % m_anOuterBlockSize[dimIdx]));
2103 : // And round up compared to arrayStartIdx, arrayStep
2104 23711 : nIncr = static_cast<size_t>(cpl::div_round_up(
2105 23711 : newIdx - indicesOuterLoop[dimIdx], arrayStep[dimIdx]));
2106 : }
2107 : else
2108 : {
2109 192 : nIncr = 1;
2110 : }
2111 23903 : indicesOuterLoop[dimIdx] += nIncr * arrayStep[dimIdx];
2112 23903 : if (indicesOuterLoop[dimIdx] >
2113 23903 : arrayStartIdx[dimIdx] + (count[dimIdx] - 1) * arrayStep[dimIdx])
2114 1577 : break;
2115 22326 : srcPtrStackOuterLoop[dimIdx] +=
2116 22326 : bufferStride[dimIdx] *
2117 22326 : static_cast<GPtrDiff_t>(nIncr * nBufferDTSize);
2118 44652 : blockIndices[dimIdx] =
2119 22326 : indicesOuterLoop[dimIdx] / m_anOuterBlockSize[dimIdx];
2120 22326 : }
2121 : }
2122 24654 : if (dimIdx > 0)
2123 24156 : goto lbl_return_to_caller;
2124 :
2125 498 : return true;
2126 : }
2127 :
2128 : /************************************************************************/
2129 : /* ZarrArray::IsEmptyBlock() */
2130 : /************************************************************************/
2131 :
2132 22824 : bool ZarrArray::IsEmptyBlock(const ZarrByteVectorQuickResize &abyBlock) const
2133 : {
2134 44463 : if (m_pabyNoData == nullptr || (m_oType.GetClass() == GEDTC_NUMERIC &&
2135 21639 : GetNoDataValueAsDouble() == 0.0))
2136 : {
2137 22516 : const size_t nBytes = abyBlock.size();
2138 22516 : size_t i = 0;
2139 25766 : for (; i + (sizeof(size_t) - 1) < nBytes; i += sizeof(size_t))
2140 : {
2141 25118 : if (*reinterpret_cast<const size_t *>(abyBlock.data() + i) != 0)
2142 : {
2143 21868 : return false;
2144 : }
2145 : }
2146 1583 : for (; i < nBytes; ++i)
2147 : {
2148 1000 : if (abyBlock[i] != 0)
2149 : {
2150 65 : return false;
2151 : }
2152 : }
2153 583 : return true;
2154 : }
2155 616 : else if (m_oType.GetClass() == GEDTC_NUMERIC &&
2156 308 : !GDALDataTypeIsComplex(m_oType.GetNumericDataType()))
2157 : {
2158 308 : const int nDTSize = static_cast<int>(m_oType.GetSize());
2159 308 : const size_t nElts = abyBlock.size() / nDTSize;
2160 308 : const auto eDT = m_oType.GetNumericDataType();
2161 616 : return GDALBufferHasOnlyNoData(
2162 308 : abyBlock.data(), GetNoDataValueAsDouble(),
2163 : nElts, // nWidth
2164 : 1, // nHeight
2165 : nElts, // nLineStride
2166 : 1, // nComponents
2167 : nDTSize * 8, // nBitsPerSample
2168 308 : GDALDataTypeIsInteger(eDT)
2169 112 : ? (GDALDataTypeIsSigned(eDT) ? GSF_SIGNED_INT
2170 : : GSF_UNSIGNED_INT)
2171 308 : : GSF_FLOATING_POINT);
2172 : }
2173 0 : return false;
2174 : }
2175 :
2176 : /************************************************************************/
2177 : /* ZarrArray::OpenBlockPresenceCache() */
2178 : /************************************************************************/
2179 :
2180 : std::shared_ptr<GDALMDArray>
2181 27967 : ZarrArray::OpenBlockPresenceCache(bool bCanCreate) const
2182 : {
2183 27967 : if (m_bHasTriedBlockCachePresenceArray)
2184 26871 : return m_poBlockCachePresenceArray;
2185 1096 : m_bHasTriedBlockCachePresenceArray = true;
2186 :
2187 1096 : if (m_nTotalInnerChunkCount == 1)
2188 391 : return nullptr;
2189 :
2190 1410 : std::string osCacheFilename;
2191 1410 : auto poRGCache = GetCacheRootGroup(bCanCreate, osCacheFilename);
2192 705 : if (!poRGCache)
2193 694 : return nullptr;
2194 :
2195 11 : const std::string osBlockPresenceArrayName(MassageName(GetFullName()) +
2196 22 : "_tile_presence");
2197 : auto poBlockPresenceArray =
2198 22 : poRGCache->OpenMDArray(osBlockPresenceArrayName);
2199 22 : const auto eByteDT = GDALExtendedDataType::Create(GDT_UInt8);
2200 11 : if (poBlockPresenceArray)
2201 : {
2202 8 : bool ok = true;
2203 8 : const auto &apoDimsCache = poBlockPresenceArray->GetDimensions();
2204 16 : if (poBlockPresenceArray->GetDataType() != eByteDT ||
2205 8 : apoDimsCache.size() != m_aoDims.size())
2206 : {
2207 0 : ok = false;
2208 : }
2209 : else
2210 : {
2211 24 : for (size_t i = 0; i < m_aoDims.size(); i++)
2212 : {
2213 16 : const auto nExpectedDimSize = cpl::div_round_up(
2214 16 : m_aoDims[i]->GetSize(), m_anInnerBlockSize[i]);
2215 16 : if (apoDimsCache[i]->GetSize() != nExpectedDimSize)
2216 : {
2217 0 : ok = false;
2218 0 : break;
2219 : }
2220 : }
2221 : }
2222 8 : if (!ok)
2223 : {
2224 0 : CPLError(CE_Failure, CPLE_NotSupported,
2225 : "Array %s in %s has not expected characteristics",
2226 : osBlockPresenceArrayName.c_str(), osCacheFilename.c_str());
2227 0 : return nullptr;
2228 : }
2229 :
2230 8 : if (!poBlockPresenceArray->GetAttribute("filling_status") &&
2231 0 : !bCanCreate)
2232 : {
2233 0 : CPLDebug(ZARR_DEBUG_KEY,
2234 : "Cache tile presence array for %s found, but filling not "
2235 : "finished",
2236 0 : GetFullName().c_str());
2237 0 : return nullptr;
2238 : }
2239 :
2240 8 : CPLDebug(ZARR_DEBUG_KEY, "Using cache tile presence for %s",
2241 8 : GetFullName().c_str());
2242 : }
2243 3 : else if (bCanCreate)
2244 : {
2245 3 : int idxDim = 0;
2246 3 : std::string osBlockSize;
2247 3 : std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
2248 9 : for (const auto &poDim : m_aoDims)
2249 : {
2250 6 : auto poNewDim = poRGCache->CreateDimension(
2251 12 : osBlockPresenceArrayName + '_' + std::to_string(idxDim),
2252 12 : std::string(), std::string(),
2253 : cpl::div_round_up(poDim->GetSize(),
2254 12 : m_anInnerBlockSize[idxDim]));
2255 6 : if (!poNewDim)
2256 0 : return nullptr;
2257 6 : apoNewDims.emplace_back(poNewDim);
2258 :
2259 6 : if (!osBlockSize.empty())
2260 3 : osBlockSize += ',';
2261 6 : constexpr GUInt64 BLOCKSIZE = 256;
2262 : osBlockSize +=
2263 6 : std::to_string(std::min(poNewDim->GetSize(), BLOCKSIZE));
2264 :
2265 6 : idxDim++;
2266 : }
2267 :
2268 3 : CPLStringList aosOptionsBlockPresence;
2269 3 : aosOptionsBlockPresence.SetNameValue("BLOCKSIZE", osBlockSize.c_str());
2270 : poBlockPresenceArray =
2271 9 : poRGCache->CreateMDArray(osBlockPresenceArrayName, apoNewDims,
2272 6 : eByteDT, aosOptionsBlockPresence.List());
2273 3 : if (!poBlockPresenceArray)
2274 : {
2275 0 : CPLError(CE_Failure, CPLE_NotSupported, "Cannot create %s in %s",
2276 : osBlockPresenceArrayName.c_str(), osCacheFilename.c_str());
2277 0 : return nullptr;
2278 : }
2279 3 : poBlockPresenceArray->SetNoDataValue(0);
2280 : }
2281 : else
2282 : {
2283 0 : return nullptr;
2284 : }
2285 :
2286 11 : m_poBlockCachePresenceArray = poBlockPresenceArray;
2287 :
2288 11 : return poBlockPresenceArray;
2289 : }
2290 :
2291 : /************************************************************************/
2292 : /* ZarrArray::BlockCachePresence() */
2293 : /************************************************************************/
2294 :
2295 5 : bool ZarrArray::BlockCachePresence()
2296 : {
2297 5 : if (m_nTotalInnerChunkCount == 1)
2298 0 : return true;
2299 :
2300 10 : const std::string osDirectoryName = GetDataDirectory();
2301 :
2302 : auto psDir = std::unique_ptr<VSIDIR, decltype(&VSICloseDir)>(
2303 10 : VSIOpenDir(osDirectoryName.c_str(), -1, nullptr), VSICloseDir);
2304 5 : if (!psDir)
2305 0 : return false;
2306 :
2307 10 : auto poBlockPresenceArray = OpenBlockPresenceCache(true);
2308 5 : if (!poBlockPresenceArray)
2309 : {
2310 0 : return false;
2311 : }
2312 :
2313 5 : if (poBlockPresenceArray->GetAttribute("filling_status"))
2314 : {
2315 2 : CPLDebug(ZARR_DEBUG_KEY,
2316 : "BlockCachePresence(): %s already filled. Nothing to do",
2317 2 : poBlockPresenceArray->GetName().c_str());
2318 2 : return true;
2319 : }
2320 :
2321 3 : const auto nDims = m_aoDims.size();
2322 6 : std::vector<GUInt64> anInnerBlockIdx(nDims);
2323 6 : std::vector<GUInt64> anInnerBlockCounter(nDims);
2324 6 : const std::vector<size_t> anCount(nDims, 1);
2325 6 : const std::vector<GInt64> anArrayStep(nDims, 0);
2326 6 : const std::vector<GPtrDiff_t> anBufferStride(nDims, 0);
2327 3 : const auto &apoDimsCache = poBlockPresenceArray->GetDimensions();
2328 6 : const auto eByteDT = GDALExtendedDataType::Create(GDT_UInt8);
2329 :
2330 3 : CPLDebug(ZARR_DEBUG_KEY,
2331 : "BlockCachePresence(): Iterating over %s to find which tiles are "
2332 : "present...",
2333 : osDirectoryName.c_str());
2334 3 : uint64_t nCounter = 0;
2335 : const char chSrcFilenameDirSeparator =
2336 3 : VSIGetDirectorySeparator(osDirectoryName.c_str())[0];
2337 30 : while (const VSIDIREntry *psEntry = VSIGetNextDirEntry(psDir.get()))
2338 : {
2339 27 : if (!VSI_ISDIR(psEntry->nMode))
2340 : {
2341 : const CPLStringList aosTokens = GetChunkIndicesFromFilename(
2342 0 : CPLString(psEntry->pszName)
2343 20 : .replaceAll(chSrcFilenameDirSeparator, '/')
2344 40 : .c_str());
2345 20 : if (aosTokens.size() == static_cast<int>(nDims))
2346 : {
2347 : // Get tile indices from filename
2348 14 : bool unexpectedIndex = false;
2349 14 : uint64_t nInnerChunksInOuter = 1;
2350 42 : for (int i = 0; i < aosTokens.size(); ++i)
2351 : {
2352 28 : if (CPLGetValueType(aosTokens[i]) != CPL_VALUE_INTEGER)
2353 : {
2354 2 : unexpectedIndex = true;
2355 : }
2356 56 : anInnerBlockIdx[i] =
2357 28 : static_cast<GUInt64>(CPLAtoGIntBig(aosTokens[i]));
2358 : const auto nInnerChunkCounterThisDim =
2359 28 : m_anCountInnerBlockInOuter[i];
2360 28 : nInnerChunksInOuter *= nInnerChunkCounterThisDim;
2361 28 : if (anInnerBlockIdx[i] >=
2362 28 : apoDimsCache[i]->GetSize() / nInnerChunkCounterThisDim)
2363 : {
2364 6 : unexpectedIndex = true;
2365 : }
2366 28 : anInnerBlockIdx[i] *= nInnerChunkCounterThisDim;
2367 : }
2368 14 : if (unexpectedIndex)
2369 : {
2370 6 : continue;
2371 : }
2372 :
2373 8 : std::fill(anInnerBlockCounter.begin(),
2374 8 : anInnerBlockCounter.end(), 0);
2375 :
2376 28 : for (uint64_t iInnerChunk = 0;
2377 28 : iInnerChunk < nInnerChunksInOuter; ++iInnerChunk)
2378 : {
2379 20 : if (iInnerChunk > 0)
2380 : {
2381 : // Update chunk coordinates
2382 12 : size_t iDim = m_anInnerBlockSize.size() - 1;
2383 : const auto nInnerChunkCounterThisDim =
2384 12 : m_anCountInnerBlockInOuter[iDim];
2385 :
2386 12 : ++anInnerBlockIdx[iDim];
2387 12 : ++anInnerBlockCounter[iDim];
2388 :
2389 16 : while (anInnerBlockCounter[iDim] ==
2390 : nInnerChunkCounterThisDim)
2391 : {
2392 4 : anInnerBlockIdx[iDim] -= nInnerChunkCounterThisDim;
2393 4 : anInnerBlockCounter[iDim] = 0;
2394 4 : --iDim;
2395 :
2396 4 : ++anInnerBlockIdx[iDim];
2397 4 : ++anInnerBlockCounter[iDim];
2398 : }
2399 : }
2400 :
2401 20 : nCounter++;
2402 20 : if ((nCounter % 1000) == 0)
2403 : {
2404 0 : CPLDebug(
2405 : ZARR_DEBUG_KEY,
2406 : "BlockCachePresence(): Listing in progress "
2407 : "(last examined %s, at least %.02f %% completed)",
2408 0 : psEntry->pszName,
2409 0 : 100.0 * double(nCounter) /
2410 0 : double(m_nTotalInnerChunkCount));
2411 : }
2412 20 : constexpr GByte byOne = 1;
2413 : // CPLDebugOnly(ZARR_DEBUG_KEY, "Marking %s has present",
2414 : // psEntry->pszName);
2415 40 : if (!poBlockPresenceArray->Write(
2416 20 : anInnerBlockIdx.data(), anCount.data(),
2417 : anArrayStep.data(), anBufferStride.data(), eByteDT,
2418 : &byOne))
2419 : {
2420 0 : return false;
2421 : }
2422 : }
2423 : }
2424 : }
2425 27 : }
2426 3 : CPLDebug(ZARR_DEBUG_KEY, "BlockCachePresence(): finished");
2427 :
2428 : // Write filling_status attribute
2429 3 : auto poAttr = poBlockPresenceArray->CreateAttribute(
2430 6 : "filling_status", {}, GDALExtendedDataType::CreateString(), nullptr);
2431 3 : if (poAttr)
2432 : {
2433 3 : if (nCounter == 0)
2434 0 : poAttr->Write("no_tile_present");
2435 3 : else if (nCounter == m_nTotalInnerChunkCount)
2436 0 : poAttr->Write("all_tiles_present");
2437 : else
2438 3 : poAttr->Write("some_tiles_missing");
2439 : }
2440 :
2441 : // Force closing
2442 3 : m_poBlockCachePresenceArray = nullptr;
2443 3 : m_bHasTriedBlockCachePresenceArray = false;
2444 :
2445 3 : return true;
2446 : }
2447 :
2448 : /************************************************************************/
2449 : /* ZarrArray::CreateAttribute() */
2450 : /************************************************************************/
2451 :
2452 120 : std::shared_ptr<GDALAttribute> ZarrArray::CreateAttribute(
2453 : const std::string &osName, const std::vector<GUInt64> &anDimensions,
2454 : const GDALExtendedDataType &oDataType, CSLConstList papszOptions)
2455 : {
2456 120 : if (!CheckValidAndErrorOutIfNot())
2457 0 : return nullptr;
2458 :
2459 120 : if (!m_bUpdatable)
2460 : {
2461 2 : CPLError(CE_Failure, CPLE_NotSupported,
2462 : "Dataset not open in update mode");
2463 2 : return nullptr;
2464 : }
2465 118 : if (anDimensions.size() >= 2)
2466 : {
2467 2 : CPLError(CE_Failure, CPLE_NotSupported,
2468 : "Cannot create attributes of dimension >= 2");
2469 2 : return nullptr;
2470 : }
2471 : return m_oAttrGroup.CreateAttribute(osName, anDimensions, oDataType,
2472 116 : papszOptions);
2473 : }
2474 :
2475 : /************************************************************************/
2476 : /* ZarrGroupBase::DeleteAttribute() */
2477 : /************************************************************************/
2478 :
2479 18 : bool ZarrArray::DeleteAttribute(const std::string &osName, CSLConstList)
2480 : {
2481 18 : if (!CheckValidAndErrorOutIfNot())
2482 0 : return false;
2483 :
2484 18 : if (!m_bUpdatable)
2485 : {
2486 6 : CPLError(CE_Failure, CPLE_NotSupported,
2487 : "Dataset not open in update mode");
2488 6 : return false;
2489 : }
2490 :
2491 12 : return m_oAttrGroup.DeleteAttribute(osName);
2492 : }
2493 :
2494 : /************************************************************************/
2495 : /* ZarrArray::SetSpatialRef() */
2496 : /************************************************************************/
2497 :
2498 44 : bool ZarrArray::SetSpatialRef(const OGRSpatialReference *poSRS)
2499 : {
2500 44 : if (!CheckValidAndErrorOutIfNot())
2501 0 : return false;
2502 :
2503 44 : if (!m_bUpdatable)
2504 : {
2505 2 : return GDALPamMDArray::SetSpatialRef(poSRS);
2506 : }
2507 42 : m_poSRS.reset();
2508 42 : if (poSRS)
2509 42 : m_poSRS.reset(poSRS->Clone());
2510 42 : m_bSRSModified = true;
2511 42 : return true;
2512 : }
2513 :
2514 : /************************************************************************/
2515 : /* ZarrArray::SetUnit() */
2516 : /************************************************************************/
2517 :
2518 9 : bool ZarrArray::SetUnit(const std::string &osUnit)
2519 : {
2520 9 : if (!CheckValidAndErrorOutIfNot())
2521 0 : return false;
2522 :
2523 9 : if (!m_bUpdatable)
2524 : {
2525 0 : CPLError(CE_Failure, CPLE_NotSupported,
2526 : "Dataset not open in update mode");
2527 0 : return false;
2528 : }
2529 9 : m_osUnit = osUnit;
2530 9 : m_bUnitModified = true;
2531 9 : return true;
2532 : }
2533 :
2534 : /************************************************************************/
2535 : /* ZarrArray::GetOffset() */
2536 : /************************************************************************/
2537 :
2538 87 : double ZarrArray::GetOffset(bool *pbHasOffset,
2539 : GDALDataType *peStorageType) const
2540 : {
2541 87 : if (pbHasOffset)
2542 87 : *pbHasOffset = m_bHasOffset;
2543 87 : if (peStorageType)
2544 0 : *peStorageType = GDT_Unknown;
2545 87 : return m_dfOffset;
2546 : }
2547 :
2548 : /************************************************************************/
2549 : /* ZarrArray::GetScale() */
2550 : /************************************************************************/
2551 :
2552 85 : double ZarrArray::GetScale(bool *pbHasScale, GDALDataType *peStorageType) const
2553 : {
2554 85 : if (pbHasScale)
2555 85 : *pbHasScale = m_bHasScale;
2556 85 : if (peStorageType)
2557 0 : *peStorageType = GDT_Unknown;
2558 85 : return m_dfScale;
2559 : }
2560 :
2561 : /************************************************************************/
2562 : /* ZarrArray::SetOffset() */
2563 : /************************************************************************/
2564 :
2565 3 : bool ZarrArray::SetOffset(double dfOffset, GDALDataType /* eStorageType */)
2566 : {
2567 3 : if (!CheckValidAndErrorOutIfNot())
2568 0 : return false;
2569 :
2570 3 : m_dfOffset = dfOffset;
2571 3 : m_bHasOffset = true;
2572 3 : m_bOffsetModified = true;
2573 3 : return true;
2574 : }
2575 :
2576 : /************************************************************************/
2577 : /* ZarrArray::SetScale() */
2578 : /************************************************************************/
2579 :
2580 3 : bool ZarrArray::SetScale(double dfScale, GDALDataType /* eStorageType */)
2581 : {
2582 3 : if (!CheckValidAndErrorOutIfNot())
2583 0 : return false;
2584 :
2585 3 : m_dfScale = dfScale;
2586 3 : m_bHasScale = true;
2587 3 : m_bScaleModified = true;
2588 3 : return true;
2589 : }
2590 :
2591 : /************************************************************************/
2592 : /* GetDimensionTypeDirection() */
2593 : /************************************************************************/
2594 :
2595 : /* static */
2596 226 : void ZarrArray::GetDimensionTypeDirection(CPLJSONObject &oAttributes,
2597 : std::string &osType,
2598 : std::string &osDirection)
2599 : {
2600 452 : std::string osUnit;
2601 678 : const auto unit = oAttributes[CF_UNITS];
2602 226 : if (unit.GetType() == CPLJSONObject::Type::String)
2603 : {
2604 52 : osUnit = unit.ToString();
2605 : }
2606 :
2607 678 : const auto oStdName = oAttributes[CF_STD_NAME];
2608 226 : if (oStdName.GetType() == CPLJSONObject::Type::String)
2609 : {
2610 156 : const auto osStdName = oStdName.ToString();
2611 52 : if (osStdName == CF_PROJ_X_COORD || osStdName == CF_LONGITUDE_STD_NAME)
2612 : {
2613 26 : osType = GDAL_DIM_TYPE_HORIZONTAL_X;
2614 26 : oAttributes.Delete(CF_STD_NAME);
2615 26 : if (osUnit == CF_DEGREES_EAST)
2616 : {
2617 24 : osDirection = "EAST";
2618 : }
2619 : }
2620 50 : else if (osStdName == CF_PROJ_Y_COORD ||
2621 24 : osStdName == CF_LATITUDE_STD_NAME)
2622 : {
2623 26 : osType = GDAL_DIM_TYPE_HORIZONTAL_Y;
2624 26 : oAttributes.Delete(CF_STD_NAME);
2625 26 : if (osUnit == CF_DEGREES_NORTH)
2626 : {
2627 24 : osDirection = "NORTH";
2628 : }
2629 : }
2630 0 : else if (osStdName == "time")
2631 : {
2632 0 : osType = GDAL_DIM_TYPE_TEMPORAL;
2633 0 : oAttributes.Delete(CF_STD_NAME);
2634 : }
2635 : }
2636 :
2637 678 : const auto osAxis = oAttributes[CF_AXIS].ToString();
2638 226 : if (osAxis == "Z")
2639 : {
2640 0 : osType = GDAL_DIM_TYPE_VERTICAL;
2641 0 : const auto osPositive = oAttributes["positive"].ToString();
2642 0 : if (osPositive == "up")
2643 : {
2644 0 : osDirection = "UP";
2645 0 : oAttributes.Delete("positive");
2646 : }
2647 0 : else if (osPositive == "down")
2648 : {
2649 0 : osDirection = "DOWN";
2650 0 : oAttributes.Delete("positive");
2651 : }
2652 0 : oAttributes.Delete(CF_AXIS);
2653 : }
2654 226 : }
2655 :
2656 : /************************************************************************/
2657 : /* GetCoordinateVariables() */
2658 : /************************************************************************/
2659 :
2660 : std::vector<std::shared_ptr<GDALMDArray>>
2661 2 : ZarrArray::GetCoordinateVariables() const
2662 : {
2663 2 : if (!CheckValidAndErrorOutIfNot())
2664 0 : return {};
2665 :
2666 4 : std::vector<std::shared_ptr<GDALMDArray>> ret;
2667 6 : const auto poCoordinates = GetAttribute("coordinates");
2668 1 : if (poCoordinates &&
2669 3 : poCoordinates->GetDataType().GetClass() == GEDTC_STRING &&
2670 1 : poCoordinates->GetDimensionCount() == 0)
2671 : {
2672 1 : const char *pszCoordinates = poCoordinates->ReadAsString();
2673 1 : if (pszCoordinates)
2674 : {
2675 2 : auto poGroup = m_poGroupWeak.lock();
2676 1 : if (!poGroup)
2677 : {
2678 0 : CPLError(CE_Failure, CPLE_AppDefined,
2679 : "Cannot access coordinate variables of %s has "
2680 : "belonging group has gone out of scope",
2681 0 : GetName().c_str());
2682 : }
2683 : else
2684 : {
2685 : const CPLStringList aosNames(
2686 2 : CSLTokenizeString2(pszCoordinates, " ", 0));
2687 3 : for (int i = 0; i < aosNames.size(); i++)
2688 : {
2689 6 : auto poCoordinateVar = poGroup->OpenMDArray(aosNames[i]);
2690 2 : if (poCoordinateVar)
2691 : {
2692 2 : ret.emplace_back(poCoordinateVar);
2693 : }
2694 : else
2695 : {
2696 0 : CPLError(CE_Warning, CPLE_AppDefined,
2697 : "Cannot find variable corresponding to "
2698 : "coordinate %s",
2699 : aosNames[i]);
2700 : }
2701 : }
2702 : }
2703 : }
2704 : }
2705 :
2706 2 : return ret;
2707 : }
2708 :
2709 : /************************************************************************/
2710 : /* Resize() */
2711 : /************************************************************************/
2712 :
2713 16 : bool ZarrArray::Resize(const std::vector<GUInt64> &anNewDimSizes,
2714 : CSLConstList /* papszOptions */)
2715 : {
2716 16 : if (!CheckValidAndErrorOutIfNot())
2717 0 : return false;
2718 :
2719 16 : if (!IsWritable())
2720 : {
2721 3 : CPLError(CE_Failure, CPLE_AppDefined,
2722 : "Resize() not supported on read-only file");
2723 3 : return false;
2724 : }
2725 :
2726 13 : const auto nDimCount = GetDimensionCount();
2727 13 : if (anNewDimSizes.size() != nDimCount)
2728 : {
2729 0 : CPLError(CE_Failure, CPLE_IllegalArg,
2730 : "Not expected number of values in anNewDimSizes.");
2731 0 : return false;
2732 : }
2733 :
2734 13 : auto &dims = GetDimensions();
2735 26 : std::vector<size_t> anGrownDimIdx;
2736 26 : std::map<GDALDimension *, GUInt64> oMapDimToSize;
2737 27 : for (size_t i = 0; i < nDimCount; ++i)
2738 : {
2739 21 : auto oIter = oMapDimToSize.find(dims[i].get());
2740 21 : if (oIter != oMapDimToSize.end() && oIter->second != anNewDimSizes[i])
2741 : {
2742 2 : CPLError(CE_Failure, CPLE_AppDefined,
2743 : "Cannot resize a dimension referenced several times "
2744 : "to different sizes");
2745 7 : return false;
2746 : }
2747 19 : if (anNewDimSizes[i] != dims[i]->GetSize())
2748 : {
2749 14 : if (anNewDimSizes[i] < dims[i]->GetSize())
2750 : {
2751 5 : CPLError(CE_Failure, CPLE_NotSupported,
2752 : "Resize() does not support shrinking the array.");
2753 5 : return false;
2754 : }
2755 :
2756 9 : oMapDimToSize[dims[i].get()] = anNewDimSizes[i];
2757 9 : anGrownDimIdx.push_back(i);
2758 : }
2759 : else
2760 : {
2761 5 : oMapDimToSize[dims[i].get()] = dims[i]->GetSize();
2762 : }
2763 : }
2764 6 : if (!anGrownDimIdx.empty())
2765 : {
2766 6 : m_bDefinitionModified = true;
2767 13 : for (size_t dimIdx : anGrownDimIdx)
2768 : {
2769 14 : auto dim = std::dynamic_pointer_cast<ZarrDimension>(dims[dimIdx]);
2770 7 : if (dim)
2771 : {
2772 7 : dim->SetSize(anNewDimSizes[dimIdx]);
2773 7 : if (dim->GetName() != dim->GetFullName())
2774 : {
2775 : // This is not a local dimension
2776 7 : m_poSharedResource->UpdateDimensionSize(dim);
2777 : }
2778 : }
2779 : else
2780 : {
2781 0 : CPLAssert(false);
2782 : }
2783 : }
2784 : }
2785 6 : return true;
2786 : }
2787 :
2788 : /************************************************************************/
2789 : /* NotifyChildrenOfRenaming() */
2790 : /************************************************************************/
2791 :
2792 15 : void ZarrArray::NotifyChildrenOfRenaming()
2793 : {
2794 15 : m_oAttrGroup.ParentRenamed(m_osFullName);
2795 15 : }
2796 :
2797 : /************************************************************************/
2798 : /* ParentRenamed() */
2799 : /************************************************************************/
2800 :
2801 9 : void ZarrArray::ParentRenamed(const std::string &osNewParentFullName)
2802 : {
2803 9 : GDALMDArray::ParentRenamed(osNewParentFullName);
2804 :
2805 9 : auto poParent = m_poGroupWeak.lock();
2806 : // The parent necessarily exist, since it notified us
2807 9 : CPLAssert(poParent);
2808 :
2809 27 : m_osFilename = CPLFormFilenameSafe(
2810 18 : CPLFormFilenameSafe(poParent->GetDirectoryName().c_str(),
2811 9 : m_osName.c_str(), nullptr)
2812 : .c_str(),
2813 9 : CPLGetFilename(m_osFilename.c_str()), nullptr);
2814 9 : }
2815 :
2816 : /************************************************************************/
2817 : /* Rename() */
2818 : /************************************************************************/
2819 :
2820 21 : bool ZarrArray::Rename(const std::string &osNewName)
2821 : {
2822 21 : if (!CheckValidAndErrorOutIfNot())
2823 0 : return false;
2824 :
2825 21 : if (!m_bUpdatable)
2826 : {
2827 6 : CPLError(CE_Failure, CPLE_NotSupported,
2828 : "Dataset not open in update mode");
2829 6 : return false;
2830 : }
2831 15 : if (!ZarrGroupBase::IsValidObjectName(osNewName))
2832 : {
2833 3 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid array name");
2834 3 : return false;
2835 : }
2836 :
2837 24 : auto poParent = m_poGroupWeak.lock();
2838 12 : if (poParent)
2839 : {
2840 12 : if (!poParent->CheckArrayOrGroupWithSameNameDoesNotExist(osNewName))
2841 6 : return false;
2842 : }
2843 :
2844 : const std::string osRootDirectoryName(
2845 12 : CPLGetDirnameSafe(CPLGetDirnameSafe(m_osFilename.c_str()).c_str()));
2846 : const std::string osOldDirectoryName = CPLFormFilenameSafe(
2847 12 : osRootDirectoryName.c_str(), m_osName.c_str(), nullptr);
2848 : const std::string osNewDirectoryName = CPLFormFilenameSafe(
2849 12 : osRootDirectoryName.c_str(), osNewName.c_str(), nullptr);
2850 :
2851 6 : if (VSIRename(osOldDirectoryName.c_str(), osNewDirectoryName.c_str()) != 0)
2852 : {
2853 0 : CPLError(CE_Failure, CPLE_AppDefined, "Renaming of %s to %s failed",
2854 : osOldDirectoryName.c_str(), osNewDirectoryName.c_str());
2855 0 : return false;
2856 : }
2857 :
2858 6 : m_poSharedResource->RenameZMetadataRecursive(osOldDirectoryName,
2859 : osNewDirectoryName);
2860 :
2861 : m_osFilename =
2862 12 : CPLFormFilenameSafe(osNewDirectoryName.c_str(),
2863 6 : CPLGetFilename(m_osFilename.c_str()), nullptr);
2864 :
2865 6 : if (poParent)
2866 : {
2867 6 : poParent->NotifyArrayRenamed(m_osName, osNewName);
2868 : }
2869 :
2870 6 : BaseRename(osNewName);
2871 :
2872 6 : return true;
2873 : }
2874 :
2875 : /************************************************************************/
2876 : /* NotifyChildrenOfDeletion() */
2877 : /************************************************************************/
2878 :
2879 8 : void ZarrArray::NotifyChildrenOfDeletion()
2880 : {
2881 8 : m_oAttrGroup.ParentDeleted();
2882 8 : }
2883 :
2884 : /************************************************************************/
2885 : /* ParseProjCRS() */
2886 : /************************************************************************/
2887 :
2888 1754 : static void ParseProjCRS(const ZarrAttributeGroup *poAttrGroup,
2889 : CPLJSONObject &oAttributes, bool bFoundProjUUID,
2890 : std::shared_ptr<OGRSpatialReference> &poSRS)
2891 : {
2892 : const auto poAttrProjCode =
2893 3515 : bFoundProjUUID ? poAttrGroup->GetAttribute("proj:code") : nullptr;
2894 : const char *pszProjCode =
2895 1754 : poAttrProjCode ? poAttrProjCode->ReadAsString() : nullptr;
2896 1754 : if (pszProjCode)
2897 : {
2898 2 : poSRS = std::make_shared<OGRSpatialReference>();
2899 2 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2900 2 : if (poSRS->SetFromUserInput(
2901 : pszProjCode,
2902 2 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
2903 : OGRERR_NONE)
2904 : {
2905 0 : poSRS.reset();
2906 : }
2907 : else
2908 : {
2909 2 : oAttributes.Delete("proj:code");
2910 : }
2911 : }
2912 : else
2913 : {
2914 : // EOP Sentinel Zarr Samples Service only
2915 5256 : const auto poAttrProjEPSG = poAttrGroup->GetAttribute("proj:epsg");
2916 1752 : if (poAttrProjEPSG)
2917 : {
2918 1 : poSRS = std::make_shared<OGRSpatialReference>();
2919 1 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2920 1 : if (poSRS->importFromEPSG(poAttrProjEPSG->ReadAsInt()) !=
2921 : OGRERR_NONE)
2922 : {
2923 0 : poSRS.reset();
2924 : }
2925 : else
2926 : {
2927 1 : oAttributes.Delete("proj:epsg");
2928 : }
2929 : }
2930 : else
2931 : {
2932 : // Both EOPF Sentinel Zarr Samples Service and geo-proj convention
2933 5253 : const auto poAttrProjWKT2 = poAttrGroup->GetAttribute("proj:wkt2");
2934 : const char *pszProjWKT2 =
2935 1751 : poAttrProjWKT2 ? poAttrProjWKT2->ReadAsString() : nullptr;
2936 1751 : if (pszProjWKT2)
2937 : {
2938 3 : poSRS = std::make_shared<OGRSpatialReference>();
2939 3 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2940 3 : if (poSRS->importFromWkt(pszProjWKT2) != OGRERR_NONE)
2941 : {
2942 0 : poSRS.reset();
2943 : }
2944 : else
2945 : {
2946 3 : oAttributes.Delete("proj:wkt2");
2947 : }
2948 : }
2949 1748 : else if (bFoundProjUUID)
2950 : {
2951 : // geo-proj convention
2952 : const auto poAttrProjPROJJSON =
2953 9 : poAttrGroup->GetAttribute("proj:projjson");
2954 : const char *pszProjPROJJSON =
2955 3 : poAttrProjPROJJSON ? poAttrProjPROJJSON->ReadAsString()
2956 3 : : nullptr;
2957 3 : if (pszProjPROJJSON)
2958 : {
2959 1 : poSRS = std::make_shared<OGRSpatialReference>();
2960 1 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2961 1 : if (poSRS->SetFromUserInput(
2962 : pszProjPROJJSON,
2963 : OGRSpatialReference::
2964 1 : SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
2965 : OGRERR_NONE)
2966 : {
2967 0 : poSRS.reset();
2968 : }
2969 : else
2970 : {
2971 1 : oAttributes.Delete("proj:projjson");
2972 : }
2973 : }
2974 : }
2975 : }
2976 : }
2977 1754 : }
2978 :
2979 : /************************************************************************/
2980 : /* ParseSpatialConventions() */
2981 : /************************************************************************/
2982 :
2983 116 : static void ParseSpatialConventions(
2984 : const std::shared_ptr<ZarrSharedResource> &poSharedResource,
2985 : const ZarrAttributeGroup *poAttrGroup, CPLJSONObject &oAttributes,
2986 : std::shared_ptr<OGRSpatialReference> &poSRS, bool &bAxisAssigned,
2987 : const std::vector<std::shared_ptr<GDALDimension>> &apoDims)
2988 : {
2989 : // From https://github.com/zarr-conventions/spatial
2990 : const auto poAttrSpatialDimensions =
2991 232 : poAttrGroup->GetAttribute("spatial:dimensions");
2992 116 : if (!poAttrSpatialDimensions)
2993 110 : return;
2994 :
2995 : const auto aosSpatialDimensions =
2996 6 : poAttrSpatialDimensions->ReadAsStringArray();
2997 6 : if (aosSpatialDimensions.size() < 2)
2998 0 : return;
2999 :
3000 6 : int iDimNameY = 0;
3001 6 : int iDimNameX = 0;
3002 : const char *pszNameY =
3003 6 : aosSpatialDimensions[aosSpatialDimensions.size() - 2];
3004 : const char *pszNameX =
3005 6 : aosSpatialDimensions[aosSpatialDimensions.size() - 1];
3006 6 : int iDim = 1;
3007 18 : for (const auto &poDim : apoDims)
3008 : {
3009 12 : if (poDim->GetName() == pszNameX)
3010 6 : iDimNameX = iDim;
3011 6 : else if (poDim->GetName() == pszNameY)
3012 6 : iDimNameY = iDim;
3013 12 : ++iDim;
3014 : }
3015 6 : if (iDimNameX == 0)
3016 : {
3017 0 : CPLError(CE_Warning, CPLE_AppDefined,
3018 : "spatial:dimensions[%d] = %s is a unknown "
3019 : "Zarr dimension",
3020 0 : static_cast<int>(aosSpatialDimensions.size() - 1), pszNameX);
3021 : }
3022 6 : if (iDimNameY == 0)
3023 : {
3024 0 : CPLError(CE_Warning, CPLE_AppDefined,
3025 : "spatial_dimensions[%d] = %s is a unknown "
3026 : "Zarr dimension",
3027 0 : static_cast<int>(aosSpatialDimensions.size() - 2), pszNameY);
3028 : }
3029 :
3030 6 : if (iDimNameX > 0 && iDimNameY > 0)
3031 : {
3032 6 : oAttributes.Delete("spatial:dimensions");
3033 :
3034 6 : if (!bAxisAssigned && poSRS)
3035 : {
3036 5 : const auto &oMapping = poSRS->GetDataAxisToSRSAxisMapping();
3037 15 : if (oMapping == std::vector<int>{2, 1} ||
3038 10 : oMapping == std::vector<int>{2, 1, 3})
3039 0 : poSRS->SetDataAxisToSRSAxisMapping({iDimNameY, iDimNameX});
3040 10 : else if (oMapping == std::vector<int>{1, 2} ||
3041 5 : oMapping == std::vector<int>{1, 2, 3})
3042 5 : poSRS->SetDataAxisToSRSAxisMapping({iDimNameX, iDimNameY});
3043 :
3044 5 : bAxisAssigned = true;
3045 : }
3046 : }
3047 :
3048 : const auto poAttrSpatialRegistration =
3049 18 : poAttrGroup->GetAttribute("spatial:registration");
3050 6 : bool bIsNodeRegistration = false;
3051 6 : if (!poAttrSpatialRegistration)
3052 3 : oAttributes.Set("spatial:registration", "pixel"); // default value
3053 : else
3054 : {
3055 : const char *pszSpatialRegistration =
3056 3 : poAttrSpatialRegistration->ReadAsString();
3057 3 : if (pszSpatialRegistration &&
3058 3 : strcmp(pszSpatialRegistration, "node") == 0)
3059 1 : bIsNodeRegistration = true;
3060 : }
3061 :
3062 : const auto poAttrSpatialTransform =
3063 18 : poAttrGroup->GetAttribute("spatial:transform");
3064 : const auto poAttrSpatialTransformType =
3065 18 : poAttrGroup->GetAttribute("spatial:transform_type");
3066 : const char *pszAttrSpatialTransformType =
3067 6 : poAttrSpatialTransformType ? poAttrSpatialTransformType->ReadAsString()
3068 6 : : nullptr;
3069 :
3070 8 : if (poAttrSpatialTransform &&
3071 2 : (!pszAttrSpatialTransformType ||
3072 8 : strcmp(pszAttrSpatialTransformType, "affine") == 0))
3073 : {
3074 12 : auto adfSpatialTransform = poAttrSpatialTransform->ReadAsDoubleArray();
3075 6 : if (adfSpatialTransform.size() == 6)
3076 : {
3077 6 : oAttributes.Delete("spatial:transform");
3078 6 : oAttributes.Delete("spatial:transform_type");
3079 :
3080 : // If we have rotation/shear coefficients, expose a gdal:geotransform
3081 : // attributes
3082 6 : if (adfSpatialTransform[1] != 0 || adfSpatialTransform[3] != 0)
3083 : {
3084 2 : if (bIsNodeRegistration)
3085 : {
3086 : // From pixel center convention to GDAL's corner convention
3087 1 : adfSpatialTransform[2] -= 0.5 * adfSpatialTransform[0] +
3088 1 : 0.5 * adfSpatialTransform[1];
3089 2 : adfSpatialTransform[5] -= 0.5 * adfSpatialTransform[3] +
3090 1 : 0.5 * adfSpatialTransform[4];
3091 : }
3092 :
3093 2 : CPLJSONArray oGeoTransform;
3094 : // Reorder coefficients to GDAL convention
3095 14 : for (int idx : {2, 0, 1, 5, 3, 4})
3096 12 : oGeoTransform.Add(adfSpatialTransform[idx]);
3097 2 : oAttributes["gdal:geotransform"] = oGeoTransform;
3098 : }
3099 : else
3100 : {
3101 4 : auto &poDimX = apoDims[iDimNameX - 1];
3102 4 : auto &poDimY = apoDims[iDimNameY - 1];
3103 4 : if (!dynamic_cast<GDALMDArrayRegularlySpaced *>(
3104 16 : poDimX->GetIndexingVariable().get()) &&
3105 4 : !dynamic_cast<GDALMDArrayRegularlySpaced *>(
3106 9 : poDimY->GetIndexingVariable().get()))
3107 : {
3108 : auto poIndexingVarX = GDALMDArrayRegularlySpaced::Create(
3109 4 : std::string(), poDimX->GetName(), poDimX,
3110 12 : adfSpatialTransform[2] + adfSpatialTransform[0] / 2,
3111 12 : adfSpatialTransform[0], 0);
3112 4 : poDimX->SetIndexingVariable(poIndexingVarX);
3113 :
3114 : // Make the shared resource hold a strong
3115 : // reference on the indexing variable,
3116 : // so that it remains available to anyone
3117 : // querying the dimension for it.
3118 4 : poSharedResource->RegisterIndexingVariable(
3119 : poDimX->GetFullName(), poIndexingVarX);
3120 :
3121 : auto poIndexingVarY = GDALMDArrayRegularlySpaced::Create(
3122 4 : std::string(), poDimY->GetName(), poDimY,
3123 12 : adfSpatialTransform[5] + adfSpatialTransform[4] / 2,
3124 8 : adfSpatialTransform[4], 0);
3125 4 : poDimY->SetIndexingVariable(poIndexingVarY);
3126 4 : poSharedResource->RegisterIndexingVariable(
3127 : poDimY->GetFullName(), poIndexingVarY);
3128 : }
3129 : }
3130 : }
3131 : else
3132 : {
3133 0 : CPLError(CE_Warning, CPLE_AppDefined,
3134 : "spatial:transform[] contains an "
3135 : "unexpected number of values: %d",
3136 0 : static_cast<int>(adfSpatialTransform.size()));
3137 : }
3138 : }
3139 : }
3140 :
3141 : /************************************************************************/
3142 : /* DetectSRSFromEOPFSampleServiceMetadata() */
3143 : /************************************************************************/
3144 :
3145 : /* This function is derived from ExtractCoordinateMetadata() of
3146 : * https://github.com/EOPF-Sample-Service/GDAL-ZARR-EOPF/blob/main/src/eopf_metadata.cpp
3147 : * released under the MIT license and
3148 : * Copyright (c) 2024 Yuvraj Adagale and contributors
3149 : *
3150 : * Note: it does not handle defaulting to EPSG:4326 as it is not clear to me
3151 : * (E. Rouault) why this would be needed, at least for Sentinel2 L1C or L2
3152 : * products
3153 : */
3154 13 : static void DetectSRSFromEOPFSampleServiceMetadata(
3155 : const std::string &osRootDirectoryName,
3156 : const std::shared_ptr<ZarrGroupBase> &poGroup,
3157 : std::shared_ptr<OGRSpatialReference> &poSRS)
3158 : {
3159 13 : const CPLJSONObject obj = poGroup->GetAttributeGroup().Serialize();
3160 :
3161 : // -----------------------------------
3162 : // STEP 1: Extract spatial reference information
3163 : // -----------------------------------
3164 :
3165 : // Find EPSG code directly or in STAC properties
3166 13 : int nEPSGCode = 0;
3167 13 : const CPLJSONObject &stacDiscovery = obj.GetObj("stac_discovery");
3168 13 : if (stacDiscovery.IsValid())
3169 : {
3170 14 : const CPLJSONObject &properties = stacDiscovery.GetObj("properties");
3171 7 : if (properties.IsValid())
3172 : {
3173 : // Try to get proj:epsg as a number first, then as a string
3174 6 : nEPSGCode = properties.GetInteger("proj:epsg", 0);
3175 6 : if (nEPSGCode <= 0)
3176 : {
3177 : nEPSGCode =
3178 4 : std::atoi(properties.GetString("proj:epsg", "").c_str());
3179 : }
3180 6 : if (nEPSGCode > 0)
3181 : {
3182 2 : CPLDebugOnly(ZARR_DEBUG_KEY,
3183 : "Found proj:epsg in STAC properties: %d",
3184 : nEPSGCode);
3185 : }
3186 : }
3187 : }
3188 :
3189 : // If not found in STAC, try top level
3190 13 : if (nEPSGCode <= 0)
3191 : {
3192 11 : nEPSGCode = obj.GetInteger("proj:epsg", obj.GetInteger("epsg", 0));
3193 11 : if (nEPSGCode <= 0)
3194 : {
3195 9 : nEPSGCode = std::atoi(
3196 18 : obj.GetString("proj:epsg", obj.GetString("epsg", "")).c_str());
3197 : }
3198 11 : if (nEPSGCode > 0)
3199 : {
3200 2 : CPLDebugOnly(ZARR_DEBUG_KEY, "Found proj:epsg at top level: %d",
3201 : nEPSGCode);
3202 : }
3203 : }
3204 :
3205 : // If still not found, simple search in common locations
3206 13 : if (nEPSGCode <= 0)
3207 : {
3208 17 : for (const auto &child : obj.GetChildren())
3209 : {
3210 8 : if (child.GetType() == CPLJSONObject::Type::Object)
3211 : {
3212 : nEPSGCode =
3213 7 : child.GetInteger("proj:epsg", child.GetInteger("epsg", 0));
3214 7 : if (nEPSGCode <= 0)
3215 : {
3216 5 : nEPSGCode = std::atoi(
3217 : child
3218 10 : .GetString("proj:epsg", child.GetString("epsg", ""))
3219 : .c_str());
3220 : }
3221 7 : if (nEPSGCode > 0)
3222 : {
3223 2 : CPLDebugOnly(ZARR_DEBUG_KEY,
3224 : "Found proj:epsg in child %s: %d",
3225 : child.GetName().c_str(), nEPSGCode);
3226 2 : break;
3227 : }
3228 : }
3229 : }
3230 : }
3231 :
3232 : // Enhanced search for STAC discovery metadata with better structure parsing
3233 13 : if (nEPSGCode <= 0 && stacDiscovery.IsValid())
3234 : {
3235 : // Try to get the full STAC item
3236 10 : const CPLJSONObject &geometry = stacDiscovery.GetObj("geometry");
3237 :
3238 5 : if (geometry.IsValid())
3239 : {
3240 2 : const CPLJSONObject &geomCrs = geometry.GetObj("crs");
3241 1 : if (geomCrs.IsValid())
3242 : {
3243 2 : const CPLJSONObject &geomProps = geomCrs.GetObj("properties");
3244 1 : if (geomProps.IsValid())
3245 : {
3246 1 : const int nGeomEpsg = geomProps.GetInteger("code", 0);
3247 1 : if (nGeomEpsg != 0)
3248 : {
3249 1 : nEPSGCode = nGeomEpsg;
3250 1 : CPLDebugOnly(ZARR_DEBUG_KEY,
3251 : "Found CRS code in STAC geometry: %d",
3252 : nEPSGCode);
3253 : }
3254 : }
3255 : }
3256 : }
3257 : }
3258 :
3259 : // Try to infer CRS from Sentinel-2 tile naming convention
3260 13 : if (nEPSGCode <= 0)
3261 : {
3262 : // Look for Sentinel-2 tile ID pattern in dataset name or metadata
3263 12 : std::string tileName;
3264 :
3265 : // First, try to extract from dataset name if it contains T##XXX pattern
3266 12 : std::string dsNameStr(CPLGetFilename(osRootDirectoryName.c_str()));
3267 6 : const size_t tilePos = dsNameStr.find("_T");
3268 6 : if (tilePos != std::string::npos && tilePos + 6 < dsNameStr.length())
3269 : {
3270 1 : tileName = dsNameStr.substr(tilePos + 1, 6); // Extract T##XXX
3271 1 : CPLDebugOnly(ZARR_DEBUG_KEY,
3272 : "Extracted tile name from dataset name: %s",
3273 : tileName.c_str());
3274 : }
3275 :
3276 : // Also check in STAC discovery metadata
3277 6 : if (tileName.empty() && stacDiscovery.IsValid())
3278 : {
3279 : const CPLJSONObject &properties =
3280 8 : stacDiscovery.GetObj("properties");
3281 4 : if (properties.IsValid())
3282 : {
3283 8 : tileName = properties.GetString(
3284 : "s2:mgrs_tile",
3285 8 : properties.GetString("mgrs_tile",
3286 12 : properties.GetString("tile_id", "")));
3287 4 : if (!tileName.empty())
3288 : {
3289 3 : CPLDebug("EOPFZARR",
3290 : "Found tile name in STAC properties: %s",
3291 : tileName.c_str());
3292 : }
3293 : }
3294 : }
3295 :
3296 : // Parse tile name to get EPSG code (T##XXX -> UTM Zone ## North/South)
3297 6 : if (!tileName.empty() && tileName.length() >= 3 && tileName[0] == 'T')
3298 : {
3299 : // Extract zone number (characters 1-2)
3300 8 : const std::string zoneStr = tileName.substr(1, 2);
3301 4 : const int zone = std::atoi(zoneStr.c_str());
3302 :
3303 4 : if (zone >= 1 && zone <= 60)
3304 : {
3305 : // Determine hemisphere from the third character
3306 : const char hemisphere =
3307 4 : tileName.length() > 3 ? tileName[3] : 'N';
3308 :
3309 : // For Sentinel-2, assume Northern hemisphere unless explicitly Southern
3310 : // Most Sentinel-2 data is Northern hemisphere
3311 : // Cf https://en.wikipedia.org/wiki/Military_Grid_Reference_System#Grid_zone_designation
3312 4 : const bool isNorth = (hemisphere >= 'N' && hemisphere <= 'X');
3313 :
3314 4 : nEPSGCode = isNorth ? (32600 + zone) : (32700 + zone);
3315 4 : CPLDebugOnly(ZARR_DEBUG_KEY,
3316 : "Inferred EPSG %d from Sentinel-2 tile %s (zone "
3317 : "%d, %s hemisphere)",
3318 : nEPSGCode, tileName.c_str(), zone,
3319 : isNorth ? "North" : "South");
3320 : }
3321 : }
3322 : }
3323 :
3324 13 : if (nEPSGCode > 0)
3325 : {
3326 11 : poSRS = std::make_shared<OGRSpatialReference>();
3327 11 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3328 11 : if (poSRS->importFromEPSG(nEPSGCode) == OGRERR_NONE)
3329 : {
3330 11 : return;
3331 : }
3332 0 : poSRS.reset();
3333 : }
3334 :
3335 : // Look for WKT
3336 4 : std::string wkt = obj.GetString("spatial_ref", "");
3337 2 : if (wkt.empty() && stacDiscovery.IsValid())
3338 : {
3339 2 : const CPLJSONObject &properties = stacDiscovery.GetObj("properties");
3340 1 : if (properties.IsValid())
3341 : {
3342 1 : wkt = properties.GetString("spatial_ref", "");
3343 : }
3344 : }
3345 2 : if (!wkt.empty())
3346 : {
3347 2 : poSRS = std::make_shared<OGRSpatialReference>();
3348 2 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3349 2 : if (poSRS->importFromWkt(wkt.c_str()) == OGRERR_NONE)
3350 : {
3351 2 : return;
3352 : }
3353 0 : poSRS.reset();
3354 : }
3355 : }
3356 :
3357 : /************************************************************************/
3358 : /* SetAttributes() */
3359 : /************************************************************************/
3360 :
3361 1730 : void ZarrArray::SetAttributes(const std::shared_ptr<ZarrGroupBase> &poGroup,
3362 : CPLJSONObject &oAttributes)
3363 : {
3364 5190 : const auto crs = oAttributes[CRS_ATTRIBUTE_NAME];
3365 1730 : std::shared_ptr<OGRSpatialReference> poSRS;
3366 1730 : if (crs.GetType() == CPLJSONObject::Type::Object)
3367 : {
3368 48 : for (const char *key : {"url", "wkt", "projjson"})
3369 : {
3370 96 : const auto item = crs[key];
3371 48 : if (item.IsValid())
3372 : {
3373 29 : poSRS = std::make_shared<OGRSpatialReference>();
3374 29 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3375 58 : if (poSRS->SetFromUserInput(
3376 58 : item.ToString().c_str(),
3377 : OGRSpatialReference::
3378 29 : SET_FROM_USER_INPUT_LIMITATIONS_get()) ==
3379 : OGRERR_NONE)
3380 : {
3381 29 : oAttributes.Delete(CRS_ATTRIBUTE_NAME);
3382 29 : break;
3383 : }
3384 0 : poSRS.reset();
3385 : }
3386 : }
3387 : }
3388 : else
3389 : {
3390 : // Check if SRS is using CF-1 conventions
3391 5103 : const auto gridMapping = oAttributes["grid_mapping"];
3392 1701 : if (gridMapping.GetType() == CPLJSONObject::Type::String)
3393 : {
3394 : const auto gridMappingArray =
3395 6 : poGroup->OpenMDArray(gridMapping.ToString());
3396 2 : if (gridMappingArray)
3397 : {
3398 2 : poSRS = std::make_shared<OGRSpatialReference>();
3399 2 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3400 4 : CPLStringList aosKeyValues;
3401 22 : for (const auto &poAttr : gridMappingArray->GetAttributes())
3402 : {
3403 20 : if (poAttr->GetDataType().GetClass() == GEDTC_STRING)
3404 : {
3405 4 : aosKeyValues.SetNameValue(poAttr->GetName().c_str(),
3406 8 : poAttr->ReadAsString());
3407 : }
3408 16 : else if (poAttr->GetDataType().GetClass() == GEDTC_NUMERIC)
3409 : {
3410 32 : std::string osVal;
3411 32 : for (double val : poAttr->ReadAsDoubleArray())
3412 : {
3413 16 : if (!osVal.empty())
3414 0 : osVal += ',';
3415 16 : osVal += CPLSPrintf("%.17g", val);
3416 : }
3417 16 : aosKeyValues.SetNameValue(poAttr->GetName().c_str(),
3418 32 : osVal.c_str());
3419 : }
3420 : }
3421 2 : if (poSRS->importFromCF1(aosKeyValues.List(), nullptr) !=
3422 : OGRERR_NONE)
3423 : {
3424 0 : poSRS.reset();
3425 : }
3426 : }
3427 : }
3428 : }
3429 :
3430 : // For EOPF Sentinel Zarr Samples Service datasets, read attributes from
3431 : // the STAC Proj extension attributes to get the CRS.
3432 : // There is also partly intersection with https://github.com/zarr-conventions/geo-proj
3433 : // For zarr-conventions/geo-proj and zarr-conventions/spatial, try first
3434 : // at array level, then parent and finally grandparent.
3435 :
3436 3460 : auto poRootGroup = std::dynamic_pointer_cast<ZarrGroupBase>(GetRootGroup());
3437 :
3438 3460 : std::vector<const ZarrAttributeGroup *> apoAttrGroup;
3439 :
3440 1730 : ZarrAttributeGroup oThisAttrGroup(std::string(),
3441 3460 : /* bContainerIsGroup = */ false);
3442 1730 : oThisAttrGroup.Init(oAttributes, /* bUpdatable=*/false);
3443 :
3444 1730 : apoAttrGroup.push_back(&oThisAttrGroup);
3445 1730 : if (GetDimensionCount() >= 2)
3446 : {
3447 1245 : apoAttrGroup.push_back(&(poGroup->GetAttributeGroup()));
3448 : // Only use root group to detect conventions
3449 1245 : if (poRootGroup)
3450 1245 : apoAttrGroup.push_back(&(poRootGroup->GetAttributeGroup()));
3451 : }
3452 :
3453 : // Look for declaration of geo-proj and spatial conventions
3454 1730 : bool bFoundSpatialUUID = false;
3455 1730 : bool bFoundProjUUID = false;
3456 5814 : for (const ZarrAttributeGroup *poAttrGroup : apoAttrGroup)
3457 : {
3458 : const auto poAttrZarrConventions =
3459 8426 : poAttrGroup->GetAttribute("zarr_conventions");
3460 4213 : if (poAttrZarrConventions)
3461 : {
3462 : const char *pszZarrConventions =
3463 129 : poAttrZarrConventions->ReadAsString();
3464 129 : if (pszZarrConventions)
3465 : {
3466 258 : CPLJSONDocument oDoc;
3467 129 : if (oDoc.LoadMemory(pszZarrConventions))
3468 : {
3469 258 : const auto oZarrConventions = oDoc.GetRoot();
3470 129 : if (oZarrConventions.GetType() ==
3471 : CPLJSONObject::Type::Array)
3472 : {
3473 : const auto oZarrConventionsArray =
3474 129 : oZarrConventions.ToArray();
3475 :
3476 : const auto hasSpatialUUIDLambda =
3477 186 : [](const CPLJSONObject &obj)
3478 : {
3479 186 : constexpr const char *SPATIAL_UUID =
3480 : "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4";
3481 186 : return obj.GetString("uuid") == SPATIAL_UUID;
3482 : };
3483 : bFoundSpatialUUID =
3484 258 : std::find_if(oZarrConventionsArray.begin(),
3485 258 : oZarrConventionsArray.end(),
3486 129 : hasSpatialUUIDLambda) !=
3487 258 : oZarrConventionsArray.end();
3488 :
3489 : const auto hasProjUUIDLambda =
3490 181 : [](const CPLJSONObject &obj)
3491 : {
3492 181 : constexpr const char *PROJ_UUID =
3493 : "f17cb550-5864-4468-aeb7-f3180cfb622f";
3494 181 : return obj.GetString("uuid") == PROJ_UUID;
3495 : };
3496 : bFoundProjUUID =
3497 258 : std::find_if(oZarrConventionsArray.begin(),
3498 258 : oZarrConventionsArray.end(),
3499 129 : hasProjUUIDLambda) !=
3500 258 : oZarrConventionsArray.end();
3501 : }
3502 : }
3503 : }
3504 129 : break;
3505 : }
3506 : }
3507 :
3508 : // If there is neither spatial nor geo-proj, just consider the current array
3509 : // for EOPF Sentinel Zarr Samples Service datasets
3510 1730 : if (!bFoundSpatialUUID && !bFoundProjUUID)
3511 1672 : apoAttrGroup.resize(1);
3512 58 : else if (apoAttrGroup.size() == 3)
3513 58 : apoAttrGroup.resize(2);
3514 :
3515 1730 : bool bAxisAssigned = false;
3516 3518 : for (const ZarrAttributeGroup *poAttrGroup : apoAttrGroup)
3517 : {
3518 1788 : if (!poSRS)
3519 : {
3520 1754 : ParseProjCRS(poAttrGroup, oAttributes, bFoundProjUUID, poSRS);
3521 : }
3522 :
3523 1788 : if (GetDimensionCount() >= 2 && bFoundSpatialUUID)
3524 : {
3525 116 : const bool bAxisAssignedBefore = bAxisAssigned;
3526 116 : ParseSpatialConventions(m_poSharedResource, poAttrGroup,
3527 : oAttributes, poSRS, bAxisAssigned,
3528 116 : m_aoDims);
3529 116 : if (bAxisAssigned && !bAxisAssignedBefore)
3530 5 : SetSRS(poSRS);
3531 :
3532 : // Note: we ignore EOPF Sentinel Zarr Samples Service "proj:transform"
3533 : // attribute, as we don't need to
3534 : // use it since the x and y dimensions are already associated with a
3535 : // 1-dimensional array with the values.
3536 : }
3537 : }
3538 :
3539 1730 : if (!poSRS && poRootGroup && oAttributes.GetObj("_eopf_attrs").IsValid())
3540 : {
3541 13 : DetectSRSFromEOPFSampleServiceMetadata(
3542 13 : m_poSharedResource->GetRootDirectoryName(), poRootGroup, poSRS);
3543 : }
3544 :
3545 1730 : if (poSRS && !bAxisAssigned)
3546 : {
3547 46 : int iDimX = 0;
3548 46 : int iDimY = 0;
3549 46 : int iCount = 1;
3550 138 : for (const auto &poDim : GetDimensions())
3551 : {
3552 92 : if (poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_X)
3553 2 : iDimX = iCount;
3554 90 : else if (poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_Y)
3555 2 : iDimY = iCount;
3556 92 : iCount++;
3557 : }
3558 46 : if ((iDimX == 0 || iDimY == 0) && GetDimensionCount() >= 2)
3559 : {
3560 44 : iDimX = static_cast<int>(GetDimensionCount());
3561 44 : iDimY = iDimX - 1;
3562 : }
3563 46 : if (iDimX > 0 && iDimY > 0)
3564 : {
3565 46 : const auto &oMapping = poSRS->GetDataAxisToSRSAxisMapping();
3566 131 : if (oMapping == std::vector<int>{2, 1} ||
3567 85 : oMapping == std::vector<int>{2, 1, 3})
3568 7 : poSRS->SetDataAxisToSRSAxisMapping({iDimY, iDimX});
3569 79 : else if (oMapping == std::vector<int>{1, 2} ||
3570 40 : oMapping == std::vector<int>{1, 2, 3})
3571 39 : poSRS->SetDataAxisToSRSAxisMapping({iDimX, iDimY});
3572 : }
3573 :
3574 46 : SetSRS(poSRS);
3575 : }
3576 :
3577 5190 : const auto unit = oAttributes[CF_UNITS];
3578 1730 : if (unit.GetType() == CPLJSONObject::Type::String)
3579 : {
3580 165 : std::string osUnit = unit.ToString();
3581 55 : oAttributes.Delete(CF_UNITS);
3582 55 : RegisterUnit(osUnit);
3583 : }
3584 :
3585 5190 : const auto offset = oAttributes[CF_ADD_OFFSET];
3586 1730 : const auto offsetType = offset.GetType();
3587 1730 : if (offsetType == CPLJSONObject::Type::Integer ||
3588 1730 : offsetType == CPLJSONObject::Type::Long ||
3589 : offsetType == CPLJSONObject::Type::Double)
3590 : {
3591 3 : double dfOffset = offset.ToDouble();
3592 3 : oAttributes.Delete(CF_ADD_OFFSET);
3593 3 : RegisterOffset(dfOffset);
3594 : }
3595 :
3596 5190 : const auto scale = oAttributes[CF_SCALE_FACTOR];
3597 1730 : const auto scaleType = scale.GetType();
3598 1730 : if (scaleType == CPLJSONObject::Type::Integer ||
3599 1730 : scaleType == CPLJSONObject::Type::Long ||
3600 : scaleType == CPLJSONObject::Type::Double)
3601 : {
3602 3 : double dfScale = scale.ToDouble();
3603 3 : oAttributes.Delete(CF_SCALE_FACTOR);
3604 3 : RegisterScale(dfScale);
3605 : }
3606 :
3607 1730 : m_oAttrGroup.Init(oAttributes, m_bUpdatable);
3608 1730 : }
3609 :
3610 : /************************************************************************/
3611 : /* SetStatistics() */
3612 : /************************************************************************/
3613 :
3614 1 : bool ZarrArray::SetStatistics(bool bApproxStats, double dfMin, double dfMax,
3615 : double dfMean, double dfStdDev,
3616 : GUInt64 nValidCount, CSLConstList papszOptions)
3617 : {
3618 2 : if (!bApproxStats && m_bUpdatable &&
3619 1 : CPLTestBool(
3620 : CSLFetchNameValueDef(papszOptions, "UPDATE_METADATA", "NO")))
3621 : {
3622 3 : auto poAttr = GetAttribute("actual_range");
3623 1 : if (!poAttr)
3624 : {
3625 : poAttr =
3626 1 : CreateAttribute("actual_range", {2}, GetDataType(), nullptr);
3627 : }
3628 1 : if (poAttr)
3629 : {
3630 2 : std::vector<GUInt64> startIdx = {0};
3631 2 : std::vector<size_t> count = {2};
3632 1 : std::vector<double> values = {dfMin, dfMax};
3633 2 : poAttr->Write(startIdx.data(), count.data(), nullptr, nullptr,
3634 2 : GDALExtendedDataType::Create(GDT_Float64),
3635 1 : values.data(), nullptr, 0);
3636 : }
3637 : }
3638 1 : return GDALPamMDArray::SetStatistics(bApproxStats, dfMin, dfMax, dfMean,
3639 1 : dfStdDev, nValidCount, papszOptions);
3640 : }
3641 :
3642 : /************************************************************************/
3643 : /* ZarrArray::IsBlockMissingFromCacheInfo() */
3644 : /************************************************************************/
3645 :
3646 27962 : bool ZarrArray::IsBlockMissingFromCacheInfo(const std::string &osFilename,
3647 : const uint64_t *blockIndices) const
3648 : {
3649 27962 : CPL_IGNORE_RET_VAL(osFilename);
3650 55924 : auto poBlockPresenceArray = OpenBlockPresenceCache(false);
3651 27962 : if (poBlockPresenceArray)
3652 : {
3653 36 : std::vector<GUInt64> anBlockIdx(m_aoDims.size());
3654 36 : const std::vector<size_t> anCount(m_aoDims.size(), 1);
3655 36 : const std::vector<GInt64> anArrayStep(m_aoDims.size(), 0);
3656 36 : const std::vector<GPtrDiff_t> anBufferStride(m_aoDims.size(), 0);
3657 36 : const auto eByteDT = GDALExtendedDataType::Create(GDT_UInt8);
3658 108 : for (size_t i = 0; i < m_aoDims.size(); ++i)
3659 : {
3660 72 : anBlockIdx[i] = static_cast<GUInt64>(blockIndices[i]);
3661 : }
3662 36 : GByte byValue = 0;
3663 72 : if (poBlockPresenceArray->Read(
3664 36 : anBlockIdx.data(), anCount.data(), anArrayStep.data(),
3665 72 : anBufferStride.data(), eByteDT, &byValue) &&
3666 36 : byValue == 0)
3667 : {
3668 26 : CPLDebugOnly(ZARR_DEBUG_KEY, "Block %s missing (=nodata)",
3669 : osFilename.c_str());
3670 26 : return true;
3671 : }
3672 : }
3673 27936 : return false;
3674 : }
3675 :
3676 : /************************************************************************/
3677 : /* ZarrArray::GetRawBlockInfo() */
3678 : /************************************************************************/
3679 :
3680 16 : bool ZarrArray::GetRawBlockInfo(const uint64_t *panBlockCoordinates,
3681 : GDALMDArrayRawBlockInfo &info) const
3682 : {
3683 16 : info.clear();
3684 33 : for (size_t i = 0; i < m_anInnerBlockSize.size(); ++i)
3685 : {
3686 20 : const auto nBlockSize = m_anInnerBlockSize[i];
3687 : const auto nBlockCount =
3688 20 : cpl::div_round_up(m_aoDims[i]->GetSize(), nBlockSize);
3689 20 : if (panBlockCoordinates[i] >= nBlockCount)
3690 : {
3691 3 : CPLError(CE_Failure, CPLE_AppDefined,
3692 : "GetRawBlockInfo() failed: array %s: "
3693 : "invalid block coordinate (%u) for dimension %u",
3694 3 : GetName().c_str(),
3695 3 : static_cast<unsigned>(panBlockCoordinates[i]),
3696 : static_cast<unsigned>(i));
3697 3 : return false;
3698 : }
3699 : }
3700 :
3701 26 : std::vector<uint64_t> anOuterBlockIndices;
3702 30 : for (size_t i = 0; i < m_anCountInnerBlockInOuter.size(); ++i)
3703 : {
3704 17 : anOuterBlockIndices.push_back(panBlockCoordinates[i] /
3705 17 : m_anCountInnerBlockInOuter[i]);
3706 : }
3707 :
3708 26 : std::string osFilename = BuildChunkFilename(anOuterBlockIndices.data());
3709 :
3710 : // For network file systems, get the streaming version of the filename,
3711 : // as we don't need arbitrary seeking in the file
3712 13 : osFilename = VSIFileManager::GetHandler(osFilename.c_str())
3713 13 : ->GetStreamingFilename(osFilename);
3714 : {
3715 13 : std::lock_guard<std::mutex> oLock(m_oMutex);
3716 13 : if (IsBlockMissingFromCacheInfo(osFilename, panBlockCoordinates))
3717 0 : return true;
3718 : }
3719 :
3720 13 : VSILFILE *fp = nullptr;
3721 : // This is the number of files returned in a S3 directory listing operation
3722 13 : const char *const apszOpenOptions[] = {"IGNORE_FILENAME_RESTRICTIONS=YES",
3723 : nullptr};
3724 13 : const auto nErrorBefore = CPLGetErrorCounter();
3725 : {
3726 : // Avoid issuing ReadDir() when a lot of files are expected
3727 : CPLConfigOptionSetter optionSetter("GDAL_DISABLE_READDIR_ON_OPEN",
3728 13 : "YES", true);
3729 13 : fp = VSIFOpenEx2L(osFilename.c_str(), "rb", 0, apszOpenOptions);
3730 : }
3731 13 : if (fp == nullptr)
3732 : {
3733 1 : if (nErrorBefore != CPLGetErrorCounter())
3734 : {
3735 0 : return false;
3736 : }
3737 : else
3738 : {
3739 : // Missing files are OK and indicate nodata_value
3740 1 : return true;
3741 : }
3742 : }
3743 12 : VSIFSeekL(fp, 0, SEEK_END);
3744 12 : const auto nFileSize = VSIFTellL(fp);
3745 12 : VSIFCloseL(fp);
3746 :
3747 : // For Kerchunk files, get information on the actual location
3748 : const CPLStringList aosMetadata(
3749 12 : VSIGetFileMetadata(osFilename.c_str(), "CHUNK_INFO", nullptr));
3750 12 : if (!aosMetadata.empty())
3751 : {
3752 2 : const char *pszFilename = aosMetadata.FetchNameValue("FILENAME");
3753 2 : if (pszFilename)
3754 1 : info.pszFilename = CPLStrdup(pszFilename);
3755 2 : info.nOffset = std::strtoull(
3756 : aosMetadata.FetchNameValueDef("OFFSET", "0"), nullptr, 10);
3757 2 : info.nSize = std::strtoull(aosMetadata.FetchNameValueDef("SIZE", "0"),
3758 : nullptr, 10);
3759 2 : const char *pszBase64 = aosMetadata.FetchNameValue("BASE64");
3760 2 : if (pszBase64)
3761 : {
3762 1 : const size_t nSizeBase64 = strlen(pszBase64) + 1;
3763 1 : info.pabyInlineData = static_cast<GByte *>(CPLMalloc(nSizeBase64));
3764 1 : memcpy(info.pabyInlineData, pszBase64, nSizeBase64);
3765 : const int nDecodedSize =
3766 1 : CPLBase64DecodeInPlace(info.pabyInlineData);
3767 1 : CPLAssert(static_cast<size_t>(nDecodedSize) ==
3768 : static_cast<size_t>(info.nSize));
3769 1 : CPL_IGNORE_RET_VAL(nDecodedSize);
3770 : }
3771 : }
3772 : else
3773 : {
3774 10 : info.pszFilename = CPLStrdup(osFilename.c_str());
3775 10 : info.nOffset = 0;
3776 10 : info.nSize = nFileSize;
3777 : }
3778 :
3779 12 : info.papszInfo = CSLDuplicate(GetRawBlockInfoInfo().List());
3780 :
3781 12 : return true;
3782 : }
3783 :
3784 : /************************************************************************/
3785 : /* ZarrArray::GetParentGroup() */
3786 : /************************************************************************/
3787 :
3788 43 : std::shared_ptr<ZarrGroupBase> ZarrArray::GetParentGroup() const
3789 : {
3790 43 : std::shared_ptr<ZarrGroupBase> poGroup = m_poParent.lock();
3791 43 : if (!poGroup)
3792 : {
3793 8 : if (auto poRootGroup = m_poSharedResource->GetRootGroup())
3794 : {
3795 4 : const auto nPos = m_osFullName.rfind('/');
3796 4 : if (nPos != 0 && nPos != std::string::npos)
3797 : {
3798 8 : poGroup = std::dynamic_pointer_cast<ZarrGroupBase>(
3799 12 : poRootGroup->OpenGroupFromFullname(
3800 12 : m_osFullName.substr(0, nPos)));
3801 : }
3802 : }
3803 : }
3804 43 : return poGroup;
3805 : }
3806 :
3807 : /************************************************************************/
3808 : /* ZarrArray::IsRegularlySpaced() */
3809 : /************************************************************************/
3810 :
3811 : // Process-level LRU cache for coordinate array regularity results.
3812 : // Keyed by root directory + array full name.
3813 : // Avoids redundant HTTP reads for immutable cloud-hosted coordinate arrays.
3814 : // Thread-safe: lru11::Cache with std::mutex handles locking internally.
3815 :
3816 : struct CoordCacheEntry
3817 : {
3818 : bool bIsRegular;
3819 : double dfStart;
3820 : double dfIncrement;
3821 : };
3822 :
3823 : static lru11::Cache<std::string, CoordCacheEntry, std::mutex> g_oCoordCache{
3824 : 128};
3825 :
3826 2 : void ZarrClearCoordinateCache()
3827 : {
3828 2 : g_oCoordCache.clear();
3829 2 : }
3830 :
3831 86 : bool ZarrArray::IsRegularlySpaced(double &dfStart, double &dfIncrement) const
3832 : {
3833 : // Only cache 1D coordinate arrays (the ones that trigger HTTP reads)
3834 86 : if (GetDimensionCount() != 1)
3835 0 : return GDALMDArray::IsRegularlySpaced(dfStart, dfIncrement);
3836 :
3837 86 : const std::string &osKey = GetFilename();
3838 :
3839 : CoordCacheEntry entry;
3840 86 : if (g_oCoordCache.tryGet(osKey, entry))
3841 : {
3842 2 : CPLDebugOnly("ZARR", "IsRegularlySpaced cache hit for %s",
3843 : osKey.c_str());
3844 2 : dfStart = entry.dfStart;
3845 2 : dfIncrement = entry.dfIncrement;
3846 2 : return entry.bIsRegular;
3847 : }
3848 :
3849 : // Cache miss: perform the full coordinate read
3850 84 : const bool bResult = GDALMDArray::IsRegularlySpaced(dfStart, dfIncrement);
3851 :
3852 84 : g_oCoordCache.insert(osKey, {bResult, dfStart, dfIncrement});
3853 :
3854 84 : CPLDebugOnly("ZARR", "IsRegularlySpaced cached for %s: %s", osKey.c_str(),
3855 : bResult ? "regular" : "irregular");
3856 :
3857 84 : return bResult;
3858 : }
|
