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