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