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