Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Hierarchical Data Format Release 5 (HDF5)
4 : * Purpose: Read BAG datasets.
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2009, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2009-2018, Even Rouault <even dot rouault at spatialys dot com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "hdf5dataset.h"
16 : #include "hdf5drivercore.h"
17 : #include "gh5_convenience.h"
18 :
19 : #include "cpl_mem_cache.h"
20 : #include "cpl_string.h"
21 : #include "cpl_time.h"
22 : #include "gdal_alg.h"
23 : #include "gdal_frmts.h"
24 : #include "gdal_pam.h"
25 : #include "gdal_priv.h"
26 : #include "gdal_rat.h"
27 : #include "iso19115_srs.h"
28 : #include "ogr_core.h"
29 : #include "ogr_spatialref.h"
30 : #include "ogrsf_frmts.h"
31 : #include "rat.h"
32 :
33 : #ifdef EMBED_RESOURCE_FILES
34 : #include "embedded_resources.h"
35 : #endif
36 :
37 : #include <cassert>
38 : #include <cmath>
39 : #include <algorithm>
40 : #include <limits>
41 : #include <map>
42 : #include <utility>
43 : #include <set>
44 :
45 : struct BAGRefinementGrid
46 : {
47 : unsigned nIndex = 0;
48 : unsigned nWidth = 0;
49 : unsigned nHeight = 0;
50 : float fResX = 0.0f;
51 : float fResY = 0.0f;
52 : float fSWX = 0.0f; // offset from (bottom left corner of) the south west
53 : // low resolution grid, in center-pixel convention
54 : float fSWY = 0.0f; // offset from (bottom left corner of) the south west
55 : // low resolution grid, in center-pixel convention
56 : };
57 :
58 : constexpr float fDEFAULT_NODATA = 1000000.0f;
59 :
60 : /************************************************************************/
61 : /* ==================================================================== */
62 : /* BAGDataset */
63 : /* ==================================================================== */
64 : /************************************************************************/
65 : class BAGDataset final : public GDALPamDataset
66 : {
67 : friend class BAGRasterBand;
68 : friend class BAGSuperGridBand;
69 : friend class BAGBaseBand;
70 : friend class BAGResampledBand;
71 : friend class BAGGeorefMDSuperGridBand;
72 : friend class BAGInterpolatedBand;
73 :
74 : bool m_bReportVertCRS = true;
75 :
76 : enum class Population
77 : {
78 : MAX,
79 : MIN,
80 : MEAN,
81 : COUNT
82 : };
83 : Population m_ePopulation = Population::MAX;
84 : bool m_bMask = false;
85 :
86 : bool m_bIsChild = false;
87 : std::vector<std::unique_ptr<BAGDataset>> m_apoOverviewDS{};
88 :
89 : std::shared_ptr<GDAL::HDF5SharedResources> m_poSharedResources{};
90 : std::shared_ptr<GDALGroup> m_poRootGroup{};
91 :
92 : std::unique_ptr<OGRLayer> m_poTrackingListLayer{};
93 :
94 : OGRSpatialReference m_oSRS{};
95 : GDALGeoTransform m_gt{};
96 :
97 : int m_nLowResWidth = 0;
98 : int m_nLowResHeight = 0;
99 :
100 : double m_dfLowResMinX = 0.0;
101 : double m_dfLowResMinY = 0.0;
102 : double m_dfLowResMaxX = 0.0;
103 : double m_dfLowResMaxY = 0.0;
104 :
105 : void LoadMetadata();
106 :
107 : char *pszXMLMetadata = nullptr;
108 : char *apszMDList[2]{};
109 :
110 : int m_nChunkXSizeVarresMD = 0;
111 : int m_nChunkYSizeVarresMD = 0;
112 : void GetVarresMetadataChunkSizes(int &nChunkXSize, int &nChunkYSize);
113 :
114 : unsigned m_nChunkSizeVarresRefinement = 0;
115 : void GetVarresRefinementChunkSize(unsigned &nChunkSize);
116 :
117 : bool ReadVarresMetadataValue(int y, int x, hid_t memspace,
118 : BAGRefinementGrid *rgrid, int height,
119 : int width);
120 :
121 : bool LookForRefinementGrids(CSLConstList l_papszOpenOptions, int nY,
122 : int nX);
123 :
124 : hid_t m_hVarresMetadata = -1;
125 : hid_t m_hVarresMetadataDataType = -1;
126 : hid_t m_hVarresMetadataDataspace = -1;
127 : hid_t m_hVarresMetadataNative = -1;
128 : std::map<int, BAGRefinementGrid> m_oMapRefinemendGrids{};
129 :
130 : CPLStringList m_aosSubdatasets{};
131 :
132 : hid_t m_hVarresRefinements = -1;
133 : hid_t m_hVarresRefinementsDataType = -1;
134 : hid_t m_hVarresRefinementsDataspace = -1;
135 : hid_t m_hVarresRefinementsNative = -1;
136 : unsigned m_nRefinementsSize = 0;
137 :
138 : unsigned m_nSuperGridRefinementStartIndex = 0;
139 :
140 : lru11::Cache<unsigned, std::vector<float>> m_oCacheRefinementValues{};
141 : const float *GetRefinementValues(unsigned nRefinementIndex);
142 :
143 : bool GetMeanSupergridsResolution(double &dfResX, double &dfResY);
144 :
145 : double m_dfResFilterMin = 0;
146 : double m_dfResFilterMax = std::numeric_limits<double>::infinity();
147 :
148 : void InitOverviewDS(BAGDataset *poParentDS, int nXSize, int nYSize);
149 :
150 : bool m_bMetadataWritten = false;
151 : CPLStringList m_aosCreationOptions{};
152 : bool WriteMetadataIfNeeded();
153 :
154 : bool OpenRaster(GDALOpenInfo *poOpenInfo, const CPLString &osFilename,
155 : bool bOpenSuperGrid, int nX, int nY, bool bIsSubdataset,
156 : const CPLString &osGeorefMetadataLayer,
157 : CPLString &outOsSubDsName);
158 : bool OpenVector();
159 :
160 154 : inline hid_t GetHDF5Handle()
161 : {
162 154 : return m_poSharedResources->m_hHDF5;
163 : }
164 :
165 : CPL_DISALLOW_COPY_ASSIGN(BAGDataset)
166 :
167 : public:
168 : BAGDataset();
169 : BAGDataset(BAGDataset *poParentDS, int nOvrFactor);
170 : BAGDataset(BAGDataset *poParentDS, int nXSize, int nYSize);
171 : ~BAGDataset() override;
172 :
173 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
174 : const OGRSpatialReference *GetSpatialRef() const override;
175 :
176 : CPLErr SetGeoTransform(const GDALGeoTransform >) override;
177 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
178 :
179 : char **GetMetadataDomainList() override;
180 : char **GetMetadata(const char *pszDomain = "") override;
181 :
182 139 : int GetLayerCount() const override
183 : {
184 139 : return m_poTrackingListLayer ? 1 : 0;
185 : }
186 :
187 : const OGRLayer *GetLayer(int idx) const override;
188 :
189 : static GDALDataset *Open(GDALOpenInfo *);
190 : static GDALDataset *OpenForCreate(GDALOpenInfo *, int nXSizeIn,
191 : int nYSizeIn, int nBandsIn,
192 : CSLConstList papszCreationOptions);
193 : static GDALDataset *CreateCopy(const char *pszFilename,
194 : GDALDataset *poSrcDS, int bStrict,
195 : char **papszOptions,
196 : GDALProgressFunc pfnProgress,
197 : void *pProgressData);
198 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
199 : int nBandsIn, GDALDataType eType,
200 : char **papszOptions);
201 :
202 : OGRErr ParseWKTFromXML(const char *pszISOXML);
203 : };
204 :
205 : /************************************************************************/
206 : /* BAGCreator */
207 : /************************************************************************/
208 :
209 : class BAGCreator
210 : {
211 : hid_t m_hdf5 = -1;
212 : hid_t m_bagRoot = -1;
213 :
214 : bool CreateBase(const char *pszFilename, char **papszOptions);
215 : bool CreateTrackingListDataset();
216 : bool CreateElevationOrUncertainty(
217 : GDALDataset *poSrcDS, int nBand, const char *pszDSName,
218 : const char *pszMaxAttrName, const char *pszMinAttrName,
219 : char **papszOptions, GDALProgressFunc pfnProgress, void *pProgressData);
220 : bool Close();
221 :
222 : public:
223 74 : BAGCreator() = default;
224 : ~BAGCreator();
225 :
226 : static bool SubstituteVariables(CPLXMLNode *psNode, char **papszDict);
227 : static CPLString GenerateMetadata(int nXSize, int nYSize,
228 : const GDALGeoTransform >,
229 : const OGRSpatialReference *poSRS,
230 : char **papszOptions);
231 : static bool CreateAndWriteMetadata(hid_t hdf5,
232 : const CPLString &osXMLMetadata);
233 :
234 : bool Create(const char *pszFilename, GDALDataset *poSrcDS,
235 : char **papszOptions, GDALProgressFunc pfnProgress,
236 : void *pProgressData);
237 :
238 : bool Create(const char *pszFilename, int nBands, GDALDataType eType,
239 : char **papszOptions);
240 : };
241 :
242 : /************************************************************************/
243 : /* ==================================================================== */
244 : /* BAGRasterBand */
245 : /* ==================================================================== */
246 : /************************************************************************/
247 : class BAGRasterBand final : public GDALPamRasterBand
248 : {
249 : friend class BAGDataset;
250 :
251 : hid_t m_hDatasetID = -1;
252 : hid_t m_hNative = -1;
253 : hid_t m_hDataspace = -1;
254 :
255 : bool m_bMinMaxSet = false;
256 : double m_dfMinimum = std::numeric_limits<double>::max();
257 : double m_dfMaximum = -std::numeric_limits<double>::max();
258 :
259 : bool m_bHasNoData = false;
260 : float m_fNoDataValue = std::numeric_limits<float>::quiet_NaN();
261 :
262 : bool CreateDatasetIfNeeded();
263 : void FinalizeDataset();
264 :
265 : public:
266 : BAGRasterBand(BAGDataset *, int);
267 : ~BAGRasterBand() override;
268 :
269 : bool Initialize(hid_t hDataset, const char *pszName);
270 :
271 : CPLErr IReadBlock(int, int, void *) override;
272 : CPLErr IWriteBlock(int, int, void *) override;
273 : double GetNoDataValue(int *) override;
274 : CPLErr SetNoDataValue(double dfNoData) override;
275 :
276 : double GetMinimum(int *pbSuccess = nullptr) override;
277 : double GetMaximum(int *pbSuccess = nullptr) override;
278 : };
279 :
280 : /************************************************************************/
281 : /* ==================================================================== */
282 : /* BAGBaseBand */
283 : /* ==================================================================== */
284 : /************************************************************************/
285 :
286 : class BAGBaseBand CPL_NON_FINAL : public GDALRasterBand
287 : {
288 : protected:
289 114 : BAGBaseBand() = default;
290 :
291 : bool m_bHasNoData = false;
292 : float m_fNoDataValue = std::numeric_limits<float>::quiet_NaN();
293 :
294 : public:
295 : double GetNoDataValue(int *) override;
296 :
297 : int GetOverviewCount() override;
298 : GDALRasterBand *GetOverview(int) override;
299 : };
300 :
301 : /************************************************************************/
302 : /* GetNoDataValue() */
303 : /************************************************************************/
304 :
305 59 : double BAGBaseBand::GetNoDataValue(int *pbSuccess)
306 :
307 : {
308 59 : if (pbSuccess)
309 43 : *pbSuccess = m_bHasNoData;
310 59 : if (m_bHasNoData)
311 58 : return m_fNoDataValue;
312 :
313 1 : return GDALRasterBand::GetNoDataValue(pbSuccess);
314 : }
315 :
316 : /************************************************************************/
317 : /* GetOverviewCount() */
318 : /************************************************************************/
319 :
320 53 : int BAGBaseBand::GetOverviewCount()
321 : {
322 53 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
323 53 : return static_cast<int>(poGDS->m_apoOverviewDS.size());
324 : }
325 :
326 : /************************************************************************/
327 : /* GetOverview() */
328 : /************************************************************************/
329 :
330 23 : GDALRasterBand *BAGBaseBand::GetOverview(int i)
331 :
332 : {
333 23 : if (i < 0 || i >= GetOverviewCount())
334 2 : return nullptr;
335 21 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
336 21 : return poGDS->m_apoOverviewDS[i]->GetRasterBand(nBand);
337 : }
338 :
339 : /************************************************************************/
340 : /* ==================================================================== */
341 : /* BAGSuperGridBand */
342 : /* ==================================================================== */
343 : /************************************************************************/
344 :
345 28 : class BAGSuperGridBand final : public BAGBaseBand
346 : {
347 : friend class BAGDataset;
348 :
349 : public:
350 : BAGSuperGridBand(BAGDataset *, int, bool bHasNoData, float fNoDataValue);
351 : ~BAGSuperGridBand() override;
352 :
353 : CPLErr IReadBlock(int, int, void *) override;
354 : };
355 :
356 : /************************************************************************/
357 : /* ==================================================================== */
358 : /* BAGResampledBand */
359 : /* ==================================================================== */
360 : /************************************************************************/
361 :
362 180 : class BAGResampledBand final : public BAGBaseBand
363 : {
364 : friend class BAGDataset;
365 :
366 : bool m_bMinMaxSet = false;
367 : double m_dfMinimum = 0.0;
368 : double m_dfMaximum = 0.0;
369 :
370 : public:
371 : BAGResampledBand(BAGDataset *, int nBandIn, bool bHasNoData,
372 : float fNoDataValue, bool bInitializeMinMax);
373 : ~BAGResampledBand() override;
374 :
375 : void InitializeMinMax();
376 :
377 : CPLErr IReadBlock(int, int, void *) override;
378 :
379 : double GetMinimum(int *pbSuccess = nullptr) override;
380 : double GetMaximum(int *pbSuccess = nullptr) override;
381 : };
382 :
383 : /************************************************************************/
384 : /* ==================================================================== */
385 : /* BAGInterpolatedBand */
386 : /* ==================================================================== */
387 : /************************************************************************/
388 :
389 20 : class BAGInterpolatedBand final : public BAGBaseBand
390 : {
391 : friend class BAGDataset;
392 :
393 : bool m_bMinMaxSet = false;
394 : double m_dfMinimum = 0.0;
395 : double m_dfMaximum = 0.0;
396 :
397 : void LoadClosestRefinedNodes(
398 : double dfX, double dfY, int iXRefinedGrid, int iYRefinedGrid,
399 : const std::vector<BAGRefinementGrid> &rgrids, int nLowResMinIdxX,
400 : int nLowResMinIdxY, int nCountLowResX, int nCountLowResY,
401 : double dfLowResMinX, double dfLowResMinY, double dfLowResResX,
402 : double dfLowResResY, std::vector<double> &adfX,
403 : std::vector<double> &adfY, std::vector<float> &afDepth,
404 : std::vector<float> &afUncrt);
405 :
406 : public:
407 : BAGInterpolatedBand(BAGDataset *, int nBandIn, bool bHasNoData,
408 : float fNoDataValue, bool bInitializeMinMax);
409 : ~BAGInterpolatedBand() override;
410 :
411 : void InitializeMinMax();
412 :
413 : CPLErr IReadBlock(int, int, void *) override;
414 :
415 : double GetMinimum(int *pbSuccess = nullptr) override;
416 : double GetMaximum(int *pbSuccess = nullptr) override;
417 : };
418 :
419 : /************************************************************************/
420 : /* BAGRasterBand() */
421 : /************************************************************************/
422 154 : BAGRasterBand::BAGRasterBand(BAGDataset *poDSIn, int nBandIn)
423 : {
424 154 : poDS = poDSIn;
425 154 : nBand = nBandIn;
426 154 : }
427 :
428 : /************************************************************************/
429 : /* ~BAGRasterBand() */
430 : /************************************************************************/
431 :
432 308 : BAGRasterBand::~BAGRasterBand()
433 : {
434 : HDF5_GLOBAL_LOCK();
435 :
436 154 : if (eAccess == GA_Update)
437 : {
438 8 : CreateDatasetIfNeeded();
439 8 : FinalizeDataset();
440 : }
441 :
442 154 : if (m_hDataspace > 0)
443 153 : H5Sclose(m_hDataspace);
444 :
445 154 : if (m_hNative > 0)
446 153 : H5Tclose(m_hNative);
447 :
448 154 : if (m_hDatasetID > 0)
449 153 : H5Dclose(m_hDatasetID);
450 308 : }
451 :
452 : /************************************************************************/
453 : /* Initialize() */
454 : /************************************************************************/
455 :
456 145 : bool BAGRasterBand::Initialize(hid_t hDatasetIDIn, const char *pszName)
457 :
458 : {
459 145 : GDALRasterBand::SetDescription(pszName);
460 :
461 145 : m_hDatasetID = hDatasetIDIn;
462 :
463 145 : const hid_t datatype = H5Dget_type(m_hDatasetID);
464 145 : m_hDataspace = H5Dget_space(m_hDatasetID);
465 145 : const int n_dims = H5Sget_simple_extent_ndims(m_hDataspace);
466 145 : m_hNative = H5Tget_native_type(datatype, H5T_DIR_ASCEND);
467 :
468 145 : eDataType = GH5_GetDataType(m_hNative);
469 :
470 145 : if (n_dims == 2)
471 : {
472 145 : hsize_t dims[2] = {static_cast<hsize_t>(0), static_cast<hsize_t>(0)};
473 145 : hsize_t maxdims[2] = {static_cast<hsize_t>(0), static_cast<hsize_t>(0)};
474 :
475 145 : H5Sget_simple_extent_dims(m_hDataspace, dims, maxdims);
476 :
477 145 : if (dims[0] > INT_MAX || dims[1] > INT_MAX)
478 : {
479 1 : CPLError(CE_Failure, CPLE_NotSupported,
480 : "At least one dimension size exceeds INT_MAX !");
481 1 : return false;
482 : }
483 144 : nRasterXSize = static_cast<int>(dims[1]);
484 144 : nRasterYSize = static_cast<int>(dims[0]);
485 : }
486 : else
487 : {
488 0 : CPLError(CE_Failure, CPLE_AppDefined, "Dataset not of rank 2.");
489 0 : return false;
490 : }
491 :
492 144 : nBlockXSize = nRasterXSize;
493 144 : nBlockYSize = 1;
494 :
495 : // Check for chunksize, and use it as blocksize for optimized reading.
496 144 : const hid_t listid = H5Dget_create_plist(hDatasetIDIn);
497 144 : if (listid > 0)
498 : {
499 144 : if (H5Pget_layout(listid) == H5D_CHUNKED)
500 : {
501 28 : hsize_t panChunkDims[3] = {static_cast<hsize_t>(0),
502 : static_cast<hsize_t>(0),
503 : static_cast<hsize_t>(0)};
504 28 : const int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
505 28 : nBlockXSize = static_cast<int>(panChunkDims[nDimSize - 1]);
506 28 : nBlockYSize = static_cast<int>(panChunkDims[nDimSize - 2]);
507 : }
508 :
509 144 : H5D_fill_value_t fillType = H5D_FILL_VALUE_UNDEFINED;
510 288 : if (H5Pfill_value_defined(listid, &fillType) >= 0 &&
511 144 : fillType == H5D_FILL_VALUE_USER_DEFINED)
512 : {
513 141 : float fNoDataValue = 0.0f;
514 141 : if (H5Pget_fill_value(listid, H5T_NATIVE_FLOAT, &fNoDataValue) >= 0)
515 : {
516 141 : m_bHasNoData = true;
517 141 : m_fNoDataValue = fNoDataValue;
518 : }
519 : }
520 :
521 144 : int nfilters = H5Pget_nfilters(listid);
522 :
523 144 : char name[120] = {};
524 144 : size_t cd_nelmts = 20;
525 144 : unsigned int cd_values[20] = {};
526 144 : unsigned int flags = 0;
527 172 : for (int i = 0; i < nfilters; i++)
528 : {
529 28 : const H5Z_filter_t filter = H5Pget_filter(
530 : listid, i, &flags, &cd_nelmts, cd_values, sizeof(name), name);
531 28 : if (filter == H5Z_FILTER_DEFLATE)
532 28 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "DEFLATE",
533 : "IMAGE_STRUCTURE");
534 0 : else if (filter == H5Z_FILTER_NBIT)
535 0 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "NBIT",
536 : "IMAGE_STRUCTURE");
537 0 : else if (filter == H5Z_FILTER_SCALEOFFSET)
538 0 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "SCALEOFFSET",
539 : "IMAGE_STRUCTURE");
540 0 : else if (filter == H5Z_FILTER_SZIP)
541 0 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "SZIP",
542 : "IMAGE_STRUCTURE");
543 : }
544 :
545 144 : H5Pclose(listid);
546 : }
547 :
548 : // Load min/max information.
549 387 : if (EQUAL(pszName, "elevation") &&
550 99 : GH5_FetchAttribute(hDatasetIDIn, "Maximum Elevation Value",
551 243 : m_dfMaximum) &&
552 87 : GH5_FetchAttribute(hDatasetIDIn, "Minimum Elevation Value",
553 87 : m_dfMinimum))
554 : {
555 87 : m_bMinMaxSet = true;
556 : }
557 157 : else if (EQUAL(pszName, "uncertainty") &&
558 43 : GH5_FetchAttribute(hDatasetIDIn, "Maximum Uncertainty Value",
559 100 : m_dfMaximum) &&
560 38 : GH5_FetchAttribute(hDatasetIDIn, "Minimum Uncertainty Value",
561 38 : m_dfMinimum))
562 : {
563 : // Some products where uncertainty band is completely set to nodata
564 : // wrongly declare minimum and maximum to 0.0.
565 38 : if (m_dfMinimum != 0.0 || m_dfMaximum != 0.0)
566 36 : m_bMinMaxSet = true;
567 : }
568 40 : else if (EQUAL(pszName, "nominal_elevation") &&
569 21 : GH5_FetchAttribute(hDatasetIDIn, "max_value", m_dfMaximum) &&
570 2 : GH5_FetchAttribute(hDatasetIDIn, "min_value", m_dfMinimum))
571 : {
572 2 : m_bMinMaxSet = true;
573 : }
574 :
575 144 : return true;
576 : }
577 :
578 : /************************************************************************/
579 : /* CreateDatasetIfNeeded() */
580 : /************************************************************************/
581 :
582 283 : bool BAGRasterBand::CreateDatasetIfNeeded()
583 : {
584 283 : if (m_hDatasetID > 0 || eAccess == GA_ReadOnly)
585 275 : return true;
586 :
587 8 : hsize_t dims[2] = {static_cast<hsize_t>(nRasterYSize),
588 8 : static_cast<hsize_t>(nRasterXSize)};
589 :
590 8 : m_hDataspace = H5_CHECK(H5Screate_simple(2, dims, nullptr));
591 8 : if (m_hDataspace < 0)
592 0 : return false;
593 :
594 8 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
595 8 : bool bDeflate = EQUAL(
596 : poGDS->m_aosCreationOptions.FetchNameValueDef("COMPRESS", "DEFLATE"),
597 : "DEFLATE");
598 : int nCompressionLevel =
599 8 : atoi(poGDS->m_aosCreationOptions.FetchNameValueDef("ZLEVEL", "6"));
600 :
601 8 : bool ret = false;
602 8 : hid_t hDataType = -1;
603 8 : hid_t hParams = -1;
604 : do
605 : {
606 8 : hDataType = H5_CHECK(H5Tcopy(H5T_NATIVE_FLOAT));
607 8 : if (hDataType < 0)
608 0 : break;
609 :
610 8 : if (H5_CHECK(H5Tset_order(hDataType, H5T_ORDER_LE)) < 0)
611 0 : break;
612 :
613 8 : hParams = H5_CHECK(H5Pcreate(H5P_DATASET_CREATE));
614 8 : if (hParams < 0)
615 0 : break;
616 :
617 8 : if (H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) < 0)
618 0 : break;
619 :
620 8 : if (H5_CHECK(H5Pset_fill_value(hParams, hDataType, &m_fNoDataValue)) <
621 : 0)
622 0 : break;
623 :
624 8 : if (H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) < 0)
625 0 : break;
626 8 : hsize_t chunk_size[2] = {static_cast<hsize_t>(nBlockYSize),
627 8 : static_cast<hsize_t>(nBlockXSize)};
628 8 : if (H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) < 0)
629 0 : break;
630 :
631 8 : if (bDeflate)
632 : {
633 8 : if (H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) < 0)
634 0 : break;
635 : }
636 :
637 8 : m_hDatasetID = H5_CHECK(H5Dcreate(poGDS->GetHDF5Handle(),
638 : nBand == 1 ? "/BAG_root/elevation"
639 : : "/BAG_root/uncertainty",
640 : hDataType, m_hDataspace, hParams));
641 8 : if (m_hDatasetID < 0)
642 0 : break;
643 :
644 8 : ret = true;
645 : } while (false);
646 :
647 8 : if (hParams >= 0)
648 8 : H5_CHECK(H5Pclose(hParams));
649 8 : if (hDataType > 0)
650 8 : H5_CHECK(H5Tclose(hDataType));
651 :
652 8 : m_hNative = H5_CHECK(H5Tcopy(H5T_NATIVE_FLOAT));
653 :
654 8 : return ret;
655 : }
656 :
657 : /************************************************************************/
658 : /* FinalizeDataset() */
659 : /************************************************************************/
660 :
661 8 : void BAGRasterBand::FinalizeDataset()
662 : {
663 8 : if (m_dfMinimum > m_dfMaximum)
664 3 : return;
665 :
666 5 : const char *pszMaxAttrName =
667 5 : nBand == 1 ? "Maximum Elevation Value" : "Maximum Uncertainty Value";
668 5 : const char *pszMinAttrName =
669 5 : nBand == 1 ? "Minimum Elevation Value" : "Minimum Uncertainty Value";
670 :
671 5 : if (!GH5_CreateAttribute(m_hDatasetID, pszMaxAttrName, m_hNative))
672 0 : return;
673 :
674 5 : if (!GH5_CreateAttribute(m_hDatasetID, pszMinAttrName, m_hNative))
675 0 : return;
676 :
677 5 : if (!GH5_WriteAttribute(m_hDatasetID, pszMaxAttrName, m_dfMaximum))
678 0 : return;
679 :
680 5 : if (!GH5_WriteAttribute(m_hDatasetID, pszMinAttrName, m_dfMinimum))
681 0 : return;
682 : }
683 :
684 : /************************************************************************/
685 : /* GetMinimum() */
686 : /************************************************************************/
687 :
688 8 : double BAGRasterBand::GetMinimum(int *pbSuccess)
689 :
690 : {
691 8 : if (m_bMinMaxSet)
692 : {
693 8 : if (pbSuccess)
694 8 : *pbSuccess = TRUE;
695 8 : return m_dfMinimum;
696 : }
697 :
698 0 : return GDALRasterBand::GetMinimum(pbSuccess);
699 : }
700 :
701 : /************************************************************************/
702 : /* GetMaximum() */
703 : /************************************************************************/
704 :
705 8 : double BAGRasterBand::GetMaximum(int *pbSuccess)
706 :
707 : {
708 8 : if (m_bMinMaxSet)
709 : {
710 8 : if (pbSuccess)
711 8 : *pbSuccess = TRUE;
712 8 : return m_dfMaximum;
713 : }
714 :
715 0 : return GDALRasterBand::GetMaximum(pbSuccess);
716 : }
717 :
718 : /************************************************************************/
719 : /* GetNoDataValue() */
720 : /************************************************************************/
721 107 : double BAGRasterBand::GetNoDataValue(int *pbSuccess)
722 :
723 : {
724 107 : if (pbSuccess)
725 97 : *pbSuccess = m_bHasNoData;
726 107 : if (m_bHasNoData)
727 107 : return m_fNoDataValue;
728 :
729 0 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
730 : }
731 :
732 : /************************************************************************/
733 : /* SetNoDataValue() */
734 : /************************************************************************/
735 :
736 4 : CPLErr BAGRasterBand::SetNoDataValue(double dfNoData)
737 : {
738 4 : if (eAccess == GA_ReadOnly)
739 0 : return GDALPamRasterBand::SetNoDataValue(dfNoData);
740 :
741 4 : if (m_hDatasetID > 0)
742 : {
743 0 : CPLError(CE_Failure, CPLE_NotSupported,
744 : "Setting the nodata value after grid values have been written "
745 : "is not supported");
746 0 : return CE_Failure;
747 : }
748 :
749 4 : m_bHasNoData = true;
750 4 : m_fNoDataValue = static_cast<float>(dfNoData);
751 4 : return CE_None;
752 : }
753 :
754 : /************************************************************************/
755 : /* IReadBlock() */
756 : /************************************************************************/
757 252 : CPLErr BAGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
758 : {
759 : HDF5_GLOBAL_LOCK();
760 :
761 252 : if (!CreateDatasetIfNeeded())
762 0 : return CE_Failure;
763 :
764 252 : const int nXOff = nBlockXOff * nBlockXSize;
765 : H5OFFSET_TYPE offset[2] = {
766 252 : static_cast<H5OFFSET_TYPE>(
767 252 : std::max(0, nRasterYSize - (nBlockYOff + 1) * nBlockYSize)),
768 252 : static_cast<H5OFFSET_TYPE>(nXOff)};
769 :
770 252 : const int nSizeOfData = static_cast<int>(H5Tget_size(m_hNative));
771 252 : memset(pImage, 0,
772 252 : static_cast<size_t>(nBlockXSize) * nBlockYSize * nSizeOfData);
773 :
774 : // Blocksize may not be a multiple of imagesize.
775 252 : hsize_t count[3] = {
776 252 : std::min(static_cast<hsize_t>(nBlockYSize), GetYSize() - offset[0]),
777 252 : std::min(static_cast<hsize_t>(nBlockXSize), GetXSize() - offset[1]),
778 252 : static_cast<hsize_t>(0)};
779 :
780 252 : if (nRasterYSize - (nBlockYOff + 1) * nBlockYSize < 0)
781 : {
782 0 : count[0] +=
783 0 : (nRasterYSize - static_cast<hsize_t>(nBlockYOff + 1) * nBlockYSize);
784 : }
785 :
786 : // Select block from file space.
787 : {
788 252 : const herr_t status = H5Sselect_hyperslab(
789 : m_hDataspace, H5S_SELECT_SET, offset, nullptr, count, nullptr);
790 252 : if (status < 0)
791 0 : return CE_Failure;
792 : }
793 :
794 : // Create memory space to receive the data.
795 252 : hsize_t col_dims[2] = {static_cast<hsize_t>(nBlockYSize),
796 252 : static_cast<hsize_t>(nBlockXSize)};
797 252 : const int rank = 2;
798 252 : const hid_t memspace = H5Screate_simple(rank, col_dims, nullptr);
799 252 : H5OFFSET_TYPE mem_offset[2] = {0, 0};
800 252 : const herr_t status = H5Sselect_hyperslab(
801 : memspace, H5S_SELECT_SET, mem_offset, nullptr, count, nullptr);
802 252 : if (status < 0)
803 : {
804 0 : H5Sclose(memspace);
805 0 : return CE_Failure;
806 : }
807 :
808 252 : const herr_t status_read = H5Dread(m_hDatasetID, m_hNative, memspace,
809 : m_hDataspace, H5P_DEFAULT, pImage);
810 :
811 252 : H5Sclose(memspace);
812 :
813 : // Y flip the data.
814 252 : const int nLinesToFlip = static_cast<int>(count[0]);
815 252 : const int nLineSize = nSizeOfData * nBlockXSize;
816 252 : GByte *const pabyTemp = static_cast<GByte *>(CPLMalloc(nLineSize));
817 252 : GByte *const pbyImage = static_cast<GByte *>(pImage);
818 :
819 316 : for (int iY = 0; iY < nLinesToFlip / 2; iY++)
820 : {
821 64 : memcpy(pabyTemp, pbyImage + iY * nLineSize, nLineSize);
822 64 : memcpy(pbyImage + iY * nLineSize,
823 64 : pbyImage + (nLinesToFlip - iY - 1) * nLineSize, nLineSize);
824 64 : memcpy(pbyImage + (nLinesToFlip - iY - 1) * nLineSize, pabyTemp,
825 : nLineSize);
826 : }
827 :
828 252 : CPLFree(pabyTemp);
829 :
830 : // Return success or failure.
831 252 : if (status_read < 0)
832 : {
833 0 : CPLError(CE_Failure, CPLE_AppDefined, "H5Dread() failed for block.");
834 0 : return CE_Failure;
835 : }
836 :
837 252 : return CE_None;
838 : }
839 :
840 : /************************************************************************/
841 : /* IWriteBlock() */
842 : /************************************************************************/
843 :
844 15 : CPLErr BAGRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)
845 : {
846 : HDF5_GLOBAL_LOCK();
847 :
848 15 : if (!CreateDatasetIfNeeded())
849 0 : return CE_Failure;
850 :
851 15 : const int nXOff = nBlockXOff * nBlockXSize;
852 15 : H5OFFSET_TYPE offset[3] = {
853 15 : static_cast<H5OFFSET_TYPE>(
854 15 : std::max(0, nRasterYSize - (nBlockYOff + 1) * nBlockYSize)),
855 15 : static_cast<H5OFFSET_TYPE>(nXOff)};
856 :
857 : // Blocksize may not be a multiple of imagesize.
858 15 : hsize_t count[3] = {
859 15 : std::min(static_cast<hsize_t>(nBlockYSize), GetYSize() - offset[0]),
860 15 : std::min(static_cast<hsize_t>(nBlockXSize), GetXSize() - offset[1])};
861 :
862 15 : if (nRasterYSize - (nBlockYOff + 1) * nBlockYSize < 0)
863 : {
864 0 : count[0] +=
865 0 : (nRasterYSize - static_cast<hsize_t>(nBlockYOff + 1) * nBlockYSize);
866 : }
867 :
868 : // Select block from file space.
869 : {
870 15 : const herr_t status = H5Sselect_hyperslab(
871 : m_hDataspace, H5S_SELECT_SET, offset, nullptr, count, nullptr);
872 15 : if (status < 0)
873 0 : return CE_Failure;
874 : }
875 :
876 : // Create memory space to receive the data.
877 15 : hsize_t col_dims[2] = {static_cast<hsize_t>(nBlockYSize),
878 15 : static_cast<hsize_t>(nBlockXSize)};
879 15 : const int rank = 2;
880 15 : const hid_t memspace = H5Screate_simple(rank, col_dims, nullptr);
881 15 : H5OFFSET_TYPE mem_offset[2] = {0, 0};
882 15 : const herr_t status = H5Sselect_hyperslab(
883 : memspace, H5S_SELECT_SET, mem_offset, nullptr, count, nullptr);
884 15 : if (status < 0)
885 : {
886 0 : H5Sclose(memspace);
887 0 : return CE_Failure;
888 : }
889 :
890 : // Y flip the data.
891 15 : const int nLinesToFlip = static_cast<int>(count[0]);
892 15 : const int nSizeOfData = static_cast<int>(H5Tget_size(m_hNative));
893 15 : const int nLineSize = nSizeOfData * nBlockXSize;
894 : GByte *const pabyTemp = static_cast<GByte *>(
895 15 : CPLMalloc(static_cast<size_t>(nLineSize) * nLinesToFlip));
896 15 : GByte *const pbyImage = static_cast<GByte *>(pImage);
897 :
898 65 : for (int iY = 0; iY < nLinesToFlip; iY++)
899 : {
900 50 : memcpy(pabyTemp + iY * nLineSize,
901 50 : pbyImage + (nLinesToFlip - iY - 1) * nLineSize, nLineSize);
902 504 : for (int iX = 0; iX < static_cast<int>(count[1]); iX++)
903 : {
904 : float f;
905 454 : GDALCopyWords(pabyTemp + iY * nLineSize + iX * nSizeOfData,
906 : eDataType, 0, &f, GDT_Float32, 0, 1);
907 454 : if (!m_bHasNoData || m_fNoDataValue != f)
908 : {
909 452 : m_dfMinimum = std::min(m_dfMinimum, static_cast<double>(f));
910 452 : m_dfMaximum = std::max(m_dfMaximum, static_cast<double>(f));
911 : }
912 : }
913 : }
914 :
915 15 : const herr_t status_write = H5Dwrite(m_hDatasetID, m_hNative, memspace,
916 : m_hDataspace, H5P_DEFAULT, pabyTemp);
917 :
918 15 : H5Sclose(memspace);
919 :
920 15 : CPLFree(pabyTemp);
921 :
922 : // Return success or failure.
923 15 : if (status_write < 0)
924 : {
925 0 : CPLError(CE_Failure, CPLE_AppDefined, "H5Dwrite() failed for block.");
926 0 : return CE_Failure;
927 : }
928 :
929 15 : return CE_None;
930 : }
931 :
932 : /************************************************************************/
933 : /* BAGSuperGridBand() */
934 : /************************************************************************/
935 :
936 14 : BAGSuperGridBand::BAGSuperGridBand(BAGDataset *poDSIn, int nBandIn,
937 14 : bool bHasNoData, float fNoDataValue)
938 : {
939 14 : poDS = poDSIn;
940 14 : nBand = nBandIn;
941 14 : nRasterXSize = poDS->GetRasterXSize();
942 14 : nRasterYSize = poDS->GetRasterYSize();
943 14 : nBlockXSize = nRasterXSize;
944 14 : nBlockYSize = 1;
945 14 : eDataType = GDT_Float32;
946 14 : GDALRasterBand::SetDescription(nBand == 1 ? "elevation" : "uncertainty");
947 14 : m_bHasNoData = bHasNoData;
948 14 : m_fNoDataValue = fNoDataValue;
949 14 : }
950 :
951 : /************************************************************************/
952 : /* ~BAGSuperGridBand() */
953 : /************************************************************************/
954 :
955 : BAGSuperGridBand::~BAGSuperGridBand() = default;
956 :
957 : /************************************************************************/
958 : /* IReadBlock() */
959 : /************************************************************************/
960 :
961 8 : CPLErr BAGSuperGridBand::IReadBlock(int, int nBlockYOff, void *pImage)
962 : {
963 : HDF5_GLOBAL_LOCK();
964 :
965 8 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
966 8 : H5OFFSET_TYPE offset[2] = {
967 : static_cast<H5OFFSET_TYPE>(0),
968 : static_cast<H5OFFSET_TYPE>(
969 8 : poGDS->m_nSuperGridRefinementStartIndex +
970 8 : static_cast<H5OFFSET_TYPE>(nRasterYSize - 1 - nBlockYOff) *
971 8 : nBlockXSize)};
972 8 : hsize_t count[2] = {1, static_cast<hsize_t>(nBlockXSize)};
973 : {
974 8 : herr_t status = H5Sselect_hyperslab(
975 : poGDS->m_hVarresRefinementsDataspace, H5S_SELECT_SET, offset,
976 : nullptr, count, nullptr);
977 8 : if (status < 0)
978 0 : return CE_Failure;
979 : }
980 :
981 : // Create memory space to receive the data.
982 8 : const hid_t memspace = H5Screate_simple(2, count, nullptr);
983 8 : H5OFFSET_TYPE mem_offset[2] = {0, 0};
984 : {
985 8 : const herr_t status = H5Sselect_hyperslab(
986 : memspace, H5S_SELECT_SET, mem_offset, nullptr, count, nullptr);
987 8 : if (status < 0)
988 : {
989 0 : H5Sclose(memspace);
990 0 : return CE_Failure;
991 : }
992 : }
993 :
994 8 : float *afBuffer = new float[2 * nBlockXSize];
995 : {
996 8 : const herr_t status = H5Dread(
997 : poGDS->m_hVarresRefinements, poGDS->m_hVarresRefinementsNative,
998 : memspace, poGDS->m_hVarresRefinementsDataspace, H5P_DEFAULT,
999 : afBuffer);
1000 8 : if (status < 0)
1001 : {
1002 0 : H5Sclose(memspace);
1003 0 : delete[] afBuffer;
1004 0 : return CE_Failure;
1005 : }
1006 : }
1007 :
1008 8 : GDALCopyWords(afBuffer + nBand - 1, GDT_Float32, 2 * sizeof(float), pImage,
1009 : GDT_Float32, sizeof(float), nBlockXSize);
1010 :
1011 8 : H5Sclose(memspace);
1012 8 : delete[] afBuffer;
1013 8 : return CE_None;
1014 : }
1015 :
1016 : /************************************************************************/
1017 : /* BAGResampledBand() */
1018 : /************************************************************************/
1019 :
1020 90 : BAGResampledBand::BAGResampledBand(BAGDataset *poDSIn, int nBandIn,
1021 : bool bHasNoData, float fNoDataValue,
1022 90 : bool bInitializeMinMax)
1023 : {
1024 90 : poDS = poDSIn;
1025 90 : nBand = nBandIn;
1026 90 : nRasterXSize = poDS->GetRasterXSize();
1027 90 : nRasterYSize = poDS->GetRasterYSize();
1028 : // Mostly for autotest purposes
1029 : const int nBlockSize =
1030 90 : std::max(1, atoi(CPLGetConfigOption("GDAL_BAG_BLOCK_SIZE", "256")));
1031 90 : nBlockXSize = std::min(nBlockSize, poDS->GetRasterXSize());
1032 90 : nBlockYSize = std::min(nBlockSize, poDS->GetRasterYSize());
1033 90 : if (poDSIn->m_bMask)
1034 : {
1035 1 : eDataType = GDT_UInt8;
1036 : }
1037 89 : else if (poDSIn->m_ePopulation == BAGDataset::Population::COUNT)
1038 : {
1039 1 : eDataType = GDT_UInt32;
1040 1 : GDALRasterBand::SetDescription("count");
1041 : }
1042 : else
1043 : {
1044 88 : m_bHasNoData = true;
1045 88 : m_fNoDataValue = bHasNoData ? fNoDataValue : fDEFAULT_NODATA;
1046 88 : eDataType = GDT_Float32;
1047 88 : GDALRasterBand::SetDescription(nBand == 1 ? "elevation"
1048 : : "uncertainty");
1049 : }
1050 90 : if (bInitializeMinMax)
1051 : {
1052 36 : InitializeMinMax();
1053 : }
1054 90 : }
1055 :
1056 : /************************************************************************/
1057 : /* ~BAGResampledBand() */
1058 : /************************************************************************/
1059 :
1060 : BAGResampledBand::~BAGResampledBand() = default;
1061 :
1062 : /************************************************************************/
1063 : /* InitializeMinMax() */
1064 : /************************************************************************/
1065 :
1066 36 : void BAGResampledBand::InitializeMinMax()
1067 : {
1068 36 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
1069 90 : if (nBand == 1 &&
1070 18 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "max_depth",
1071 54 : m_dfMaximum) &&
1072 18 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "min_depth",
1073 18 : m_dfMinimum))
1074 : {
1075 18 : m_bMinMaxSet = true;
1076 : }
1077 54 : else if (nBand == 2 &&
1078 18 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "max_uncrt",
1079 36 : m_dfMaximum) &&
1080 18 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "min_uncrt",
1081 18 : m_dfMinimum))
1082 : {
1083 18 : m_bMinMaxSet = true;
1084 : }
1085 36 : }
1086 :
1087 : /************************************************************************/
1088 : /* GetMinimum() */
1089 : /************************************************************************/
1090 :
1091 2 : double BAGResampledBand::GetMinimum(int *pbSuccess)
1092 :
1093 : {
1094 2 : if (m_bMinMaxSet)
1095 : {
1096 2 : if (pbSuccess)
1097 2 : *pbSuccess = TRUE;
1098 2 : return m_dfMinimum;
1099 : }
1100 :
1101 0 : return GDALRasterBand::GetMinimum(pbSuccess);
1102 : }
1103 :
1104 : /************************************************************************/
1105 : /* GetMaximum() */
1106 : /************************************************************************/
1107 :
1108 2 : double BAGResampledBand::GetMaximum(int *pbSuccess)
1109 :
1110 : {
1111 2 : if (m_bMinMaxSet)
1112 : {
1113 2 : if (pbSuccess)
1114 2 : *pbSuccess = TRUE;
1115 2 : return m_dfMaximum;
1116 : }
1117 :
1118 0 : return GDALRasterBand::GetMaximum(pbSuccess);
1119 : }
1120 :
1121 : /************************************************************************/
1122 : /* IReadBlock() */
1123 : /************************************************************************/
1124 :
1125 289 : CPLErr BAGResampledBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1126 : void *pImage)
1127 : {
1128 : HDF5_GLOBAL_LOCK();
1129 :
1130 289 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
1131 : #ifdef DEBUG_VERBOSE
1132 : CPLDebug(
1133 : "BAG",
1134 : "IReadBlock: nRasterXSize=%d, nBlockXOff=%d, nBlockYOff=%d, nBand=%d",
1135 : nRasterXSize, nBlockXOff, nBlockYOff, nBand);
1136 : #endif
1137 :
1138 289 : const float fNoDataValue = m_fNoDataValue;
1139 289 : const float fNoSuperGridValue = m_fNoDataValue;
1140 289 : float *depthsPtr = nullptr;
1141 289 : float *uncrtPtr = nullptr;
1142 :
1143 289 : GDALRasterBlock *poBlock = nullptr;
1144 289 : if (poGDS->nBands == 2)
1145 : {
1146 287 : if (nBand == 1)
1147 : {
1148 287 : depthsPtr = static_cast<float *>(pImage);
1149 287 : poBlock = poGDS->GetRasterBand(2)->GetLockedBlockRef(
1150 287 : nBlockXOff, nBlockYOff, TRUE);
1151 287 : if (poBlock)
1152 : {
1153 287 : uncrtPtr = static_cast<float *>(poBlock->GetDataRef());
1154 : }
1155 : else
1156 : {
1157 0 : return CE_Failure;
1158 : }
1159 : }
1160 : else
1161 : {
1162 0 : uncrtPtr = static_cast<float *>(pImage);
1163 0 : poBlock = poGDS->GetRasterBand(1)->GetLockedBlockRef(
1164 0 : nBlockXOff, nBlockYOff, TRUE);
1165 0 : if (poBlock)
1166 : {
1167 0 : depthsPtr = static_cast<float *>(poBlock->GetDataRef());
1168 : }
1169 : else
1170 : {
1171 0 : return CE_Failure;
1172 : }
1173 : }
1174 : }
1175 :
1176 289 : if (depthsPtr)
1177 : {
1178 287 : GDALCopyWords(&fNoSuperGridValue, GDT_Float32, 0, depthsPtr,
1179 : GDT_Float32, static_cast<int>(sizeof(float)),
1180 287 : nBlockXSize * nBlockYSize);
1181 : }
1182 :
1183 289 : if (uncrtPtr)
1184 : {
1185 287 : GDALCopyWords(&fNoSuperGridValue, GDT_Float32, 0, uncrtPtr, GDT_Float32,
1186 : static_cast<int>(sizeof(float)),
1187 287 : nBlockXSize * nBlockYSize);
1188 : }
1189 :
1190 578 : std::vector<int> counts;
1191 289 : if (poGDS->m_bMask)
1192 : {
1193 1 : CPLAssert(pImage); // to make CLang Static Analyzer happy
1194 1 : memset(pImage, 0, static_cast<size_t>(nBlockXSize) * nBlockYSize);
1195 : }
1196 288 : else if (poGDS->m_ePopulation == BAGDataset::Population::MEAN)
1197 : {
1198 1 : counts.resize(static_cast<size_t>(nBlockXSize) * nBlockYSize);
1199 : }
1200 287 : else if (poGDS->m_ePopulation == BAGDataset::Population::COUNT)
1201 : {
1202 1 : CPLAssert(pImage); // to make CLang Static Analyzer happy
1203 1 : memset(pImage, 0,
1204 1 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
1205 1 : GDALGetDataTypeSizeBytes(eDataType));
1206 : }
1207 :
1208 : const int nReqCountX =
1209 289 : std::min(nBlockXSize, nRasterXSize - nBlockXOff * nBlockXSize);
1210 : const int nReqCountY =
1211 289 : std::min(nBlockYSize, nRasterYSize - nBlockYOff * nBlockYSize);
1212 : // Compute extent of block in georeferenced coordinates
1213 : double dfBlockMinX =
1214 289 : poGDS->m_gt[0] + nBlockXOff * nBlockXSize * poGDS->m_gt[1];
1215 289 : double dfBlockMaxX = dfBlockMinX + nReqCountX * poGDS->m_gt[1];
1216 : double dfBlockMaxY =
1217 289 : poGDS->m_gt[3] + nBlockYOff * nBlockYSize * poGDS->m_gt[5];
1218 289 : double dfBlockMinY = dfBlockMaxY + nReqCountY * poGDS->m_gt[5];
1219 :
1220 : // Compute min/max indices of intersecting supergrids (origin bottom-left)
1221 289 : const double dfLowResResX =
1222 289 : (poGDS->m_dfLowResMaxX - poGDS->m_dfLowResMinX) / poGDS->m_nLowResWidth;
1223 289 : const double dfLowResResY =
1224 289 : (poGDS->m_dfLowResMaxY - poGDS->m_dfLowResMinY) /
1225 289 : poGDS->m_nLowResHeight;
1226 : int nLowResMinIdxX =
1227 578 : std::max(0, static_cast<int>((dfBlockMinX - poGDS->m_dfLowResMinX) /
1228 289 : dfLowResResX));
1229 : int nLowResMinIdxY =
1230 578 : std::max(0, static_cast<int>((dfBlockMinY - poGDS->m_dfLowResMinY) /
1231 289 : dfLowResResY));
1232 : int nLowResMaxIdxX = std::min(
1233 578 : poGDS->m_nLowResWidth - 1,
1234 289 : static_cast<int>((dfBlockMaxX - poGDS->m_dfLowResMinX) / dfLowResResX));
1235 : int nLowResMaxIdxY = std::min(
1236 578 : poGDS->m_nLowResHeight - 1,
1237 289 : static_cast<int>((dfBlockMaxY - poGDS->m_dfLowResMinY) / dfLowResResY));
1238 :
1239 : // Create memory space to receive the data.
1240 289 : const int nCountLowResX = nLowResMaxIdxX - nLowResMinIdxX + 1;
1241 289 : const int nCountLowResY = nLowResMaxIdxY - nLowResMinIdxY + 1;
1242 289 : hsize_t countVarresMD[2] = {static_cast<hsize_t>(nCountLowResY),
1243 289 : static_cast<hsize_t>(nCountLowResX)};
1244 289 : const hid_t memspaceVarresMD = H5Screate_simple(2, countVarresMD, nullptr);
1245 289 : H5OFFSET_TYPE mem_offset[2] = {static_cast<H5OFFSET_TYPE>(0),
1246 : static_cast<H5OFFSET_TYPE>(0)};
1247 289 : if (H5Sselect_hyperslab(memspaceVarresMD, H5S_SELECT_SET, mem_offset,
1248 289 : nullptr, countVarresMD, nullptr) < 0)
1249 : {
1250 0 : H5Sclose(memspaceVarresMD);
1251 0 : if (poBlock != nullptr)
1252 : {
1253 0 : poBlock->DropLock();
1254 0 : poBlock = nullptr;
1255 : }
1256 0 : return CE_Failure;
1257 : }
1258 :
1259 289 : std::vector<BAGRefinementGrid> rgrids(static_cast<size_t>(nCountLowResY) *
1260 578 : nCountLowResX);
1261 289 : if (!(poGDS->ReadVarresMetadataValue(nLowResMinIdxY, nLowResMinIdxX,
1262 : memspaceVarresMD, rgrids.data(),
1263 : nCountLowResY, nCountLowResX)))
1264 : {
1265 0 : H5Sclose(memspaceVarresMD);
1266 0 : if (poBlock != nullptr)
1267 : {
1268 0 : poBlock->DropLock();
1269 0 : poBlock = nullptr;
1270 : }
1271 0 : return CE_Failure;
1272 : }
1273 :
1274 289 : H5Sclose(memspaceVarresMD);
1275 :
1276 289 : CPLErr eErr = CE_None;
1277 722 : for (int y = nLowResMinIdxY; y <= nLowResMaxIdxY; y++)
1278 : {
1279 1336 : for (int x = nLowResMinIdxX; x <= nLowResMaxIdxX; x++)
1280 : {
1281 903 : const auto &rgrid = rgrids[(y - nLowResMinIdxY) * nCountLowResX +
1282 903 : (x - nLowResMinIdxX)];
1283 903 : if (rgrid.nWidth == 0)
1284 : {
1285 80 : continue;
1286 : }
1287 903 : const float gridRes = std::max(rgrid.fResX, rgrid.fResY);
1288 903 : if (!(gridRes > poGDS->m_dfResFilterMin &&
1289 851 : gridRes <= poGDS->m_dfResFilterMax))
1290 : {
1291 80 : continue;
1292 : }
1293 :
1294 : // Super grid bounding box with pixel-center convention
1295 823 : const double dfMinX =
1296 823 : poGDS->m_dfLowResMinX + x * dfLowResResX + rgrid.fSWX;
1297 823 : const double dfMaxX =
1298 823 : dfMinX + (rgrid.nWidth - 1) * static_cast<double>(rgrid.fResX);
1299 823 : const double dfMinY =
1300 823 : poGDS->m_dfLowResMinY + y * dfLowResResY + rgrid.fSWY;
1301 823 : const double dfMaxY =
1302 823 : dfMinY + (rgrid.nHeight - 1) * static_cast<double>(rgrid.fResY);
1303 :
1304 : // Intersection of super grid with block
1305 823 : const double dfInterMinX = std::max(dfBlockMinX, dfMinX);
1306 823 : const double dfInterMinY = std::max(dfBlockMinY, dfMinY);
1307 823 : const double dfInterMaxX = std::min(dfBlockMaxX, dfMaxX);
1308 823 : const double dfInterMaxY = std::min(dfBlockMaxY, dfMaxY);
1309 :
1310 : // Min/max indices in the super grid
1311 : const int nMinSrcX = std::max(
1312 823 : 0, static_cast<int>((dfInterMinX - dfMinX) / rgrid.fResX));
1313 : const int nMinSrcY = std::max(
1314 823 : 0, static_cast<int>((dfInterMinY - dfMinY) / rgrid.fResY));
1315 : // Need to use ceil due to numerical imprecision
1316 : const int nMaxSrcX =
1317 1646 : std::min(static_cast<int>(rgrid.nWidth) - 1,
1318 1646 : static_cast<int>(
1319 823 : std::ceil((dfInterMaxX - dfMinX) / rgrid.fResX)));
1320 : const int nMaxSrcY =
1321 1646 : std::min(static_cast<int>(rgrid.nHeight) - 1,
1322 1646 : static_cast<int>(
1323 823 : std::ceil((dfInterMaxY - dfMinY) / rgrid.fResY)));
1324 : #ifdef DEBUG_VERBOSE
1325 : CPLDebug(
1326 : "BAG",
1327 : "y = %d, x = %d, minx = %d, miny = %d, maxx = %d, maxy = %d", y,
1328 : x, nMinSrcX, nMinSrcY, nMaxSrcX, nMaxSrcY);
1329 : #endif
1330 823 : const double dfCstX = (dfMinX - dfBlockMinX) / poGDS->m_gt[1];
1331 823 : const double dfMulX = rgrid.fResX / poGDS->m_gt[1];
1332 :
1333 3709 : for (int super_y = nMinSrcY; super_y <= nMaxSrcY; super_y++)
1334 : {
1335 2886 : const double dfSrcY =
1336 2886 : dfMinY + super_y * static_cast<double>(rgrid.fResY);
1337 : const int nTargetY = static_cast<int>(
1338 2886 : std::floor((dfBlockMaxY - dfSrcY) / -poGDS->m_gt[5]));
1339 2886 : if (!(nTargetY >= 0 && nTargetY < nReqCountY))
1340 : {
1341 752 : continue;
1342 : }
1343 :
1344 2134 : const unsigned nTargetIdxBase = nTargetY * nBlockXSize;
1345 2134 : const unsigned nRefinementIndexase =
1346 2134 : rgrid.nIndex + super_y * rgrid.nWidth;
1347 :
1348 11648 : for (int super_x = nMinSrcX; super_x <= nMaxSrcX; super_x++)
1349 : {
1350 : /*
1351 : const double dfSrcX = dfMinX + super_x * rgrid.fResX;
1352 : const int nTargetX = static_cast<int>(std::floor(
1353 : (dfSrcX - dfBlockMinX) / poGDS->m_gt[1]));
1354 : */
1355 9514 : const int nTargetX =
1356 9514 : static_cast<int>(std::floor(dfCstX + super_x * dfMulX));
1357 9514 : if (!(nTargetX >= 0 && nTargetX < nReqCountX))
1358 : {
1359 1260 : continue;
1360 : }
1361 :
1362 8254 : const unsigned nTargetIdx = nTargetIdxBase + nTargetX;
1363 8254 : if (poGDS->m_bMask)
1364 : {
1365 502 : static_cast<GByte *>(pImage)[nTargetIdx] = 255;
1366 502 : continue;
1367 : }
1368 :
1369 7752 : if (poGDS->m_ePopulation == BAGDataset::Population::COUNT)
1370 : {
1371 556 : static_cast<GUInt32 *>(pImage)[nTargetIdx]++;
1372 556 : continue;
1373 : }
1374 :
1375 7196 : CPLAssert(depthsPtr);
1376 7196 : CPLAssert(uncrtPtr);
1377 :
1378 7196 : const unsigned nRefinementIndex =
1379 7196 : nRefinementIndexase + super_x;
1380 : const auto *pafRefValues =
1381 7196 : poGDS->GetRefinementValues(nRefinementIndex);
1382 7196 : if (!pafRefValues)
1383 : {
1384 0 : eErr = CE_Failure;
1385 0 : goto end;
1386 : }
1387 :
1388 7196 : float depth = pafRefValues[0];
1389 7196 : if (depth == fNoDataValue)
1390 : {
1391 0 : if (depthsPtr[nTargetIdx] == fNoSuperGridValue)
1392 : {
1393 0 : depthsPtr[nTargetIdx] = fNoDataValue;
1394 : }
1395 0 : continue;
1396 : }
1397 :
1398 7196 : if (poGDS->m_ePopulation == BAGDataset::Population::MEAN)
1399 : {
1400 :
1401 556 : if (counts[nTargetIdx] == 0)
1402 : {
1403 133 : depthsPtr[nTargetIdx] = depth;
1404 : }
1405 : else
1406 : {
1407 423 : depthsPtr[nTargetIdx] += depth;
1408 : }
1409 556 : counts[nTargetIdx]++;
1410 :
1411 556 : auto uncrt = pafRefValues[1];
1412 556 : auto &target_uncrt_ptr = uncrtPtr[nTargetIdx];
1413 556 : if (uncrt > target_uncrt_ptr ||
1414 408 : target_uncrt_ptr == fNoDataValue)
1415 : {
1416 281 : target_uncrt_ptr = uncrt;
1417 : }
1418 : }
1419 6640 : else if ((poGDS->m_ePopulation ==
1420 6084 : BAGDataset::Population::MAX &&
1421 6084 : depth > depthsPtr[nTargetIdx]) ||
1422 5371 : (poGDS->m_ePopulation ==
1423 556 : BAGDataset::Population::MIN &&
1424 556 : depth < depthsPtr[nTargetIdx]) ||
1425 5145 : depthsPtr[nTargetIdx] == fNoDataValue ||
1426 2172 : depthsPtr[nTargetIdx] == fNoSuperGridValue)
1427 : {
1428 4468 : depthsPtr[nTargetIdx] = depth;
1429 4468 : uncrtPtr[nTargetIdx] = pafRefValues[1];
1430 : }
1431 : }
1432 : }
1433 : }
1434 : }
1435 :
1436 289 : if (poGDS->m_ePopulation == BAGDataset::Population::MEAN && depthsPtr)
1437 : {
1438 151 : for (int i = 0; i < nBlockXSize * nBlockYSize; i++)
1439 : {
1440 150 : if (counts[i])
1441 : {
1442 133 : depthsPtr[i] /= counts[i];
1443 : }
1444 : }
1445 : }
1446 :
1447 288 : end:
1448 289 : if (poBlock != nullptr)
1449 : {
1450 287 : poBlock->DropLock();
1451 287 : poBlock = nullptr;
1452 : }
1453 :
1454 289 : return eErr;
1455 : }
1456 :
1457 : /************************************************************************/
1458 : /* BAGInterpolatedBand() */
1459 : /************************************************************************/
1460 :
1461 10 : BAGInterpolatedBand::BAGInterpolatedBand(BAGDataset *poDSIn, int nBandIn,
1462 : bool bHasNoData, float fNoDataValue,
1463 10 : bool bInitializeMinMax)
1464 : {
1465 10 : poDS = poDSIn;
1466 10 : nBand = nBandIn;
1467 10 : nRasterXSize = poDS->GetRasterXSize();
1468 10 : nRasterYSize = poDS->GetRasterYSize();
1469 : // Mostly for autotest purposes
1470 : const int nBlockSize =
1471 10 : std::max(1, atoi(CPLGetConfigOption("GDAL_BAG_BLOCK_SIZE", "256")));
1472 10 : nBlockXSize = std::min(nBlockSize, poDS->GetRasterXSize());
1473 10 : nBlockYSize = std::min(nBlockSize, poDS->GetRasterYSize());
1474 :
1475 10 : m_bHasNoData = true;
1476 10 : m_fNoDataValue = bHasNoData ? fNoDataValue : fDEFAULT_NODATA;
1477 10 : eDataType = GDT_Float32;
1478 10 : GDALRasterBand::SetDescription(nBand == 1 ? "elevation" : "uncertainty");
1479 10 : if (bInitializeMinMax)
1480 : {
1481 2 : InitializeMinMax();
1482 : }
1483 10 : }
1484 :
1485 : /************************************************************************/
1486 : /* ~BAGInterpolatedBand() */
1487 : /************************************************************************/
1488 :
1489 : BAGInterpolatedBand::~BAGInterpolatedBand() = default;
1490 :
1491 : /************************************************************************/
1492 : /* InitializeMinMax() */
1493 : /************************************************************************/
1494 :
1495 2 : void BAGInterpolatedBand::InitializeMinMax()
1496 : {
1497 2 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
1498 5 : if (nBand == 1 &&
1499 1 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "max_depth",
1500 3 : m_dfMaximum) &&
1501 0 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "min_depth",
1502 0 : m_dfMinimum))
1503 : {
1504 0 : m_bMinMaxSet = true;
1505 : }
1506 5 : else if (nBand == 2 &&
1507 1 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "max_uncrt",
1508 3 : m_dfMaximum) &&
1509 0 : GH5_FetchAttribute(poGDS->m_hVarresRefinements, "min_uncrt",
1510 0 : m_dfMinimum))
1511 : {
1512 0 : m_bMinMaxSet = true;
1513 : }
1514 2 : }
1515 :
1516 : /************************************************************************/
1517 : /* GetMinimum() */
1518 : /************************************************************************/
1519 :
1520 2 : double BAGInterpolatedBand::GetMinimum(int *pbSuccess)
1521 :
1522 : {
1523 2 : if (m_bMinMaxSet)
1524 : {
1525 0 : if (pbSuccess)
1526 0 : *pbSuccess = TRUE;
1527 0 : return m_dfMinimum;
1528 : }
1529 :
1530 2 : return GDALRasterBand::GetMinimum(pbSuccess);
1531 : }
1532 :
1533 : /************************************************************************/
1534 : /* GetMaximum() */
1535 : /************************************************************************/
1536 :
1537 2 : double BAGInterpolatedBand::GetMaximum(int *pbSuccess)
1538 :
1539 : {
1540 2 : if (m_bMinMaxSet)
1541 : {
1542 0 : if (pbSuccess)
1543 0 : *pbSuccess = TRUE;
1544 0 : return m_dfMaximum;
1545 : }
1546 :
1547 2 : return GDALRasterBand::GetMaximum(pbSuccess);
1548 : }
1549 :
1550 : /************************************************************************/
1551 : /* BarycentricInterpolation() */
1552 : /************************************************************************/
1553 :
1554 203 : static bool BarycentricInterpolation(double dfX, double dfY,
1555 : const double adfX[3], const double adfY[3],
1556 : double &dfCoord0, double &dfCoord1,
1557 : double &dfCoord2)
1558 : {
1559 : /* See https://en.wikipedia.org/wiki/Barycentric_coordinate_system */
1560 203 : const double dfDenom = (adfY[1] - adfY[2]) * (adfX[0] - adfX[2]) +
1561 203 : (adfX[2] - adfX[1]) * (adfY[0] - adfY[2]);
1562 203 : if (fabs(dfDenom) < 1e-5)
1563 : {
1564 : // Degenerate triangle
1565 0 : return false;
1566 : }
1567 203 : const double dfTerm0X = adfY[1] - adfY[2];
1568 203 : const double dfTerm0Y = adfX[2] - adfX[1];
1569 203 : const double dfTerm1X = adfY[2] - adfY[0];
1570 203 : const double dfTerm1Y = adfX[0] - adfX[2];
1571 203 : const double dfDiffX = dfX - adfX[2];
1572 203 : const double dfDiffY = dfY - adfY[2];
1573 203 : dfCoord0 = (dfTerm0X * dfDiffX + dfTerm0Y * dfDiffY) / dfDenom;
1574 203 : dfCoord1 = (dfTerm1X * dfDiffX + dfTerm1Y * dfDiffY) / dfDenom;
1575 203 : dfCoord2 = 1 - dfCoord0 - dfCoord1;
1576 203 : return (dfCoord0 >= 0 && dfCoord0 <= 1 && dfCoord1 >= 0 && dfCoord1 <= 1 &&
1577 406 : dfCoord2 >= 0 && dfCoord2 <= 1);
1578 : }
1579 :
1580 : /************************************************************************/
1581 : /* SQ() */
1582 : /************************************************************************/
1583 :
1584 205823 : static inline double SQ(double v)
1585 : {
1586 205823 : return v * v;
1587 : }
1588 :
1589 : /************************************************************************/
1590 : /* IReadBlock() */
1591 : /************************************************************************/
1592 :
1593 5 : CPLErr BAGInterpolatedBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1594 : void *pImage)
1595 : {
1596 : HDF5_GLOBAL_LOCK();
1597 :
1598 5 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
1599 : #ifdef DEBUG_VERBOSE
1600 : CPLDebug(
1601 : "BAG",
1602 : "IReadBlock: nRasterXSize=%d, nBlockXOff=%d, nBlockYOff=%d, nBand=%d",
1603 : nRasterXSize, nBlockXOff, nBlockYOff, nBand);
1604 : #endif
1605 :
1606 5 : float *depthsPtr = nullptr;
1607 5 : float *uncrtPtr = nullptr;
1608 :
1609 5 : GDALRasterBlock *poBlock = nullptr;
1610 5 : if (nBand == 1)
1611 : {
1612 5 : depthsPtr = static_cast<float *>(pImage);
1613 5 : poBlock = poGDS->GetRasterBand(2)->GetLockedBlockRef(nBlockXOff,
1614 5 : nBlockYOff, TRUE);
1615 5 : if (poBlock)
1616 : {
1617 5 : uncrtPtr = static_cast<float *>(poBlock->GetDataRef());
1618 : }
1619 : }
1620 : else
1621 : {
1622 0 : uncrtPtr = static_cast<float *>(pImage);
1623 0 : poBlock = poGDS->GetRasterBand(1)->GetLockedBlockRef(nBlockXOff,
1624 0 : nBlockYOff, TRUE);
1625 0 : if (poBlock)
1626 : {
1627 0 : depthsPtr = static_cast<float *>(poBlock->GetDataRef());
1628 : }
1629 : }
1630 :
1631 5 : if (!depthsPtr || !uncrtPtr)
1632 0 : return CE_Failure;
1633 :
1634 5 : GDALCopyWords(&m_fNoDataValue, GDT_Float32, 0, depthsPtr, GDT_Float32,
1635 5 : static_cast<int>(sizeof(float)), nBlockXSize * nBlockYSize);
1636 :
1637 5 : GDALCopyWords(&m_fNoDataValue, GDT_Float32, 0, uncrtPtr, GDT_Float32,
1638 5 : static_cast<int>(sizeof(float)), nBlockXSize * nBlockYSize);
1639 :
1640 : const int nReqCountX =
1641 5 : std::min(nBlockXSize, nRasterXSize - nBlockXOff * nBlockXSize);
1642 : const int nReqCountY =
1643 5 : std::min(nBlockYSize, nRasterYSize - nBlockYOff * nBlockYSize);
1644 : // Compute extent of block in georeferenced coordinates
1645 : const double dfBlockMinX =
1646 5 : poGDS->m_gt[0] + nBlockXOff * nBlockXSize * poGDS->m_gt[1];
1647 5 : const double dfBlockMaxX = dfBlockMinX + nReqCountX * poGDS->m_gt[1];
1648 : const double dfBlockMaxY =
1649 5 : poGDS->m_gt[3] + nBlockYOff * nBlockYSize * poGDS->m_gt[5];
1650 5 : const double dfBlockMinY = dfBlockMaxY + nReqCountY * poGDS->m_gt[5];
1651 :
1652 : // Compute min/max indices of intersecting supergrids (origin bottom-left)
1653 : // We add a margin of (dfLowResResX, dfLowResResY) to be able to
1654 : // properly interpolate if the edges of a GDAL block align with the edge
1655 : // of a low-res cell
1656 5 : const double dfLowResResX =
1657 5 : (poGDS->m_dfLowResMaxX - poGDS->m_dfLowResMinX) / poGDS->m_nLowResWidth;
1658 5 : const double dfLowResResY =
1659 5 : (poGDS->m_dfLowResMaxY - poGDS->m_dfLowResMinY) /
1660 5 : poGDS->m_nLowResHeight;
1661 : const int nLowResMinIdxX = std::max(
1662 10 : 0,
1663 10 : static_cast<int>((dfBlockMinX - poGDS->m_dfLowResMinX - dfLowResResX) /
1664 5 : dfLowResResX));
1665 : const int nLowResMinIdxY = std::max(
1666 10 : 0,
1667 10 : static_cast<int>((dfBlockMinY - poGDS->m_dfLowResMinY - dfLowResResY) /
1668 5 : dfLowResResY));
1669 : const int nLowResMaxIdxX = std::min(
1670 10 : poGDS->m_nLowResWidth - 1,
1671 10 : static_cast<int>((dfBlockMaxX - poGDS->m_dfLowResMinX + dfLowResResX) /
1672 5 : dfLowResResX));
1673 : const int nLowResMaxIdxY = std::min(
1674 10 : poGDS->m_nLowResHeight - 1,
1675 10 : static_cast<int>((dfBlockMaxY - poGDS->m_dfLowResMinY + dfLowResResY) /
1676 5 : dfLowResResY));
1677 5 : if (nLowResMinIdxX >= poGDS->m_nLowResWidth ||
1678 5 : nLowResMinIdxY >= poGDS->m_nLowResHeight || nLowResMaxIdxX < 0 ||
1679 : nLowResMaxIdxY < 0)
1680 : {
1681 0 : return CE_None;
1682 : }
1683 :
1684 : // Create memory space to receive the data.
1685 5 : const int nCountLowResX = nLowResMaxIdxX - nLowResMinIdxX + 1;
1686 5 : const int nCountLowResY = nLowResMaxIdxY - nLowResMinIdxY + 1;
1687 5 : hsize_t countVarresMD[2] = {static_cast<hsize_t>(nCountLowResY),
1688 5 : static_cast<hsize_t>(nCountLowResX)};
1689 5 : const hid_t memspaceVarresMD = H5Screate_simple(2, countVarresMD, nullptr);
1690 5 : H5OFFSET_TYPE mem_offset[2] = {static_cast<H5OFFSET_TYPE>(0),
1691 : static_cast<H5OFFSET_TYPE>(0)};
1692 5 : if (H5Sselect_hyperslab(memspaceVarresMD, H5S_SELECT_SET, mem_offset,
1693 5 : nullptr, countVarresMD, nullptr) < 0)
1694 : {
1695 0 : H5Sclose(memspaceVarresMD);
1696 0 : if (poBlock != nullptr)
1697 : {
1698 0 : poBlock->DropLock();
1699 0 : poBlock = nullptr;
1700 : }
1701 0 : return CE_Failure;
1702 : }
1703 :
1704 5 : std::vector<BAGRefinementGrid> rgrids(static_cast<size_t>(nCountLowResY) *
1705 10 : nCountLowResX);
1706 5 : if (!(poGDS->ReadVarresMetadataValue(nLowResMinIdxY, nLowResMinIdxX,
1707 : memspaceVarresMD, rgrids.data(),
1708 : nCountLowResY, nCountLowResX)))
1709 : {
1710 0 : H5Sclose(memspaceVarresMD);
1711 0 : if (poBlock != nullptr)
1712 : {
1713 0 : poBlock->DropLock();
1714 0 : poBlock = nullptr;
1715 : }
1716 0 : return CE_Failure;
1717 : }
1718 :
1719 5 : H5Sclose(memspaceVarresMD);
1720 :
1721 5 : CPLErr eErr = CE_None;
1722 5 : const double dfLowResMinX = poGDS->m_dfLowResMinX;
1723 5 : const double dfLowResMinY = poGDS->m_dfLowResMinY;
1724 :
1725 : // georeferenced (X,Y) coordinates of the source nodes from the refinement
1726 : //grids
1727 10 : std::vector<double> adfX, adfY;
1728 : // Depth and uncertainty values from the source nodes from the refinement
1729 : //grids
1730 10 : std::vector<float> afDepth, afUncrt;
1731 : // Work variable to sort adfX, adfY, afDepth, afUncrt w.r.t to their
1732 : // distance with the (dfX, dfY) point to be interpolated
1733 5 : std::vector<int> anIndices;
1734 :
1735 : // Maximum distance of candidate source nodes to the point to be
1736 : // interpolated
1737 5 : const double dfMaxDistance = 0.5 * std::max(dfLowResResX, dfLowResResY);
1738 :
1739 253 : for (int y = 0; y < nReqCountY; ++y)
1740 : {
1741 : // Y georeference ordinate of the center of the cell to interpolate
1742 248 : const double dfY = dfBlockMaxY + (y + 0.5) * poGDS->m_gt[5];
1743 : // Y index of the corresponding refinement grid
1744 248 : const int iYRefinedGrid =
1745 248 : static_cast<int>(floor((dfY - dfLowResMinY) / dfLowResResY));
1746 248 : if (iYRefinedGrid < nLowResMinIdxY || iYRefinedGrid > nLowResMaxIdxY)
1747 0 : continue;
1748 30928 : for (int x = 0; x < nReqCountX; ++x)
1749 : {
1750 : // X georeference ordinate of the center of the cell to interpolate
1751 30680 : const double dfX = dfBlockMinX + (x + 0.5) * poGDS->m_gt[1];
1752 : // X index of the corresponding refinement grid
1753 30680 : const int iXRefinedGrid =
1754 30680 : static_cast<int>((dfX - dfLowResMinX) / dfLowResResX);
1755 30680 : if (iXRefinedGrid < nLowResMinIdxX ||
1756 : iXRefinedGrid > nLowResMaxIdxX)
1757 0 : continue;
1758 :
1759 : // Correspond refinement grid
1760 : const auto &rgrid =
1761 30680 : rgrids[(iYRefinedGrid - nLowResMinIdxY) * nCountLowResX +
1762 30680 : (iXRefinedGrid - nLowResMinIdxX)];
1763 30680 : if (rgrid.nWidth == 0)
1764 : {
1765 0 : continue;
1766 : }
1767 30680 : const float gridRes = std::max(rgrid.fResX, rgrid.fResY);
1768 30680 : if (!(gridRes > poGDS->m_dfResFilterMin &&
1769 30680 : gridRes <= poGDS->m_dfResFilterMax))
1770 : {
1771 0 : continue;
1772 : }
1773 :
1774 : // (dfMinRefinedX, dfMinRefinedY) is the georeferenced coordinate
1775 : // of the bottom-left corner of the refinement grid
1776 30680 : const double dfMinRefinedX =
1777 30680 : dfLowResMinX + iXRefinedGrid * dfLowResResX + rgrid.fSWX;
1778 30680 : const double dfMinRefinedY =
1779 30680 : dfLowResMinY + iYRefinedGrid * dfLowResResY + rgrid.fSWY;
1780 :
1781 : // (iXInRefinedGrid, iYInRefinedGrid) is the index of the cell
1782 : // within the refinement grid into which (dfX, dfY) falls into.
1783 30680 : const int iXInRefinedGrid =
1784 30680 : static_cast<int>(floor((dfX - dfMinRefinedX) / rgrid.fResX));
1785 30680 : const int iYInRefinedGrid =
1786 30680 : static_cast<int>(floor((dfY - dfMinRefinedY) / rgrid.fResY));
1787 :
1788 30680 : if (iXInRefinedGrid >= 0 &&
1789 30680 : iXInRefinedGrid < static_cast<int>(rgrid.nWidth) - 1 &&
1790 2088 : iYInRefinedGrid >= 0 &&
1791 2088 : iYInRefinedGrid < static_cast<int>(rgrid.nHeight) - 1)
1792 : {
1793 : // The point to interpolate is fully within a single refinement
1794 : // grid
1795 126 : const unsigned nRefinementIndexBase =
1796 126 : rgrid.nIndex + iYInRefinedGrid * rgrid.nWidth +
1797 126 : iXInRefinedGrid;
1798 : float d[2][2];
1799 : float u[2][2];
1800 126 : int nCountNoData = 0;
1801 :
1802 : // Load the depth and uncertainty values of the 4 nodes of
1803 : // the refinement grid surrounding the center of the target
1804 : // cell.
1805 378 : for (int j = 0; j < 2; ++j)
1806 : {
1807 756 : for (int i = 0; i < 2; ++i)
1808 : {
1809 504 : const unsigned nRefinementIndex =
1810 504 : nRefinementIndexBase + j * rgrid.nWidth + i;
1811 : const auto pafRefValues =
1812 504 : poGDS->GetRefinementValues(nRefinementIndex);
1813 504 : if (!pafRefValues)
1814 : {
1815 0 : eErr = CE_Failure;
1816 0 : goto end;
1817 : }
1818 504 : d[j][i] = pafRefValues[0];
1819 504 : u[j][i] = pafRefValues[1];
1820 504 : if (d[j][i] == m_fNoDataValue)
1821 121 : ++nCountNoData;
1822 : }
1823 : }
1824 :
1825 : // Compute the relative distance of the point to be
1826 : // interpolated compared to the closest bottom-left most node
1827 : // of the refinement grid.
1828 : // (alphaX,alphaY)=(0,0): point to be interpolated matches the
1829 : // closest bottom-left most node
1830 : // (alphaX,alphaY)=(1,1): point to be interpolated matches the
1831 : // closest top-right most node
1832 : const double alphaX =
1833 126 : fmod(dfX - dfMinRefinedX, rgrid.fResX) / rgrid.fResX;
1834 : const double alphaY =
1835 126 : fmod(dfY - dfMinRefinedY, rgrid.fResY) / rgrid.fResY;
1836 126 : if (nCountNoData == 0)
1837 : {
1838 : // If the 4 nodes of the supergrid around the point
1839 : // to be interpolated are valid, do bilinear interpolation
1840 : #define BILINEAR_INTERP(var) \
1841 : ((1 - alphaY) * ((1 - alphaX) * var[0][0] + alphaX * var[0][1]) + \
1842 : alphaY * ((1 - alphaX) * var[1][0] + alphaX * var[1][1]))
1843 :
1844 5 : depthsPtr[y * nBlockXSize + x] =
1845 5 : static_cast<float>(BILINEAR_INTERP(d));
1846 5 : uncrtPtr[y * nBlockXSize + x] =
1847 5 : static_cast<float>(BILINEAR_INTERP(u));
1848 : }
1849 121 : else if (nCountNoData == 1)
1850 : {
1851 : // If only one of the 4 nodes is at nodata, determine if the
1852 : // point to be interpolated is within the triangle formed
1853 : // by the remaining 3 valid nodes.
1854 : // If so, do barycentric interpolation
1855 121 : adfX.resize(3);
1856 121 : adfY.resize(3);
1857 121 : afDepth.resize(3);
1858 121 : afUncrt.resize(3);
1859 :
1860 121 : int idx = 0;
1861 363 : for (int j = 0; j < 2; ++j)
1862 : {
1863 726 : for (int i = 0; i < 2; ++i)
1864 : {
1865 484 : if (d[j][i] != m_fNoDataValue)
1866 : {
1867 363 : CPLAssert(idx < 3);
1868 363 : adfX[idx] = i;
1869 363 : adfY[idx] = j;
1870 363 : afDepth[idx] = d[j][i];
1871 363 : afUncrt[idx] = u[j][i];
1872 363 : ++idx;
1873 : }
1874 : }
1875 : }
1876 121 : CPLAssert(idx == 3);
1877 : double dfCoord0;
1878 : double dfCoord1;
1879 : double dfCoord2;
1880 121 : if (BarycentricInterpolation(alphaX, alphaY, adfX.data(),
1881 121 : adfY.data(), dfCoord0,
1882 : dfCoord1, dfCoord2))
1883 : {
1884 : // Inside triangle
1885 98 : depthsPtr[y * nBlockXSize + x] = static_cast<float>(
1886 98 : dfCoord0 * afDepth[0] + dfCoord1 * afDepth[1] +
1887 98 : dfCoord2 * afDepth[2]);
1888 98 : uncrtPtr[y * nBlockXSize + x] = static_cast<float>(
1889 98 : dfCoord0 * afUncrt[0] + dfCoord1 * afUncrt[1] +
1890 98 : dfCoord2 * afUncrt[2]);
1891 : }
1892 126 : }
1893 : // else: 2 or more nodes invalid. Target point is set at nodata
1894 : }
1895 : else
1896 : {
1897 : // Point to interpolate is on an edge or corner of the
1898 : // refinement grid
1899 30554 : adfX.clear();
1900 30554 : adfY.clear();
1901 30554 : afDepth.clear();
1902 30554 : afUncrt.clear();
1903 :
1904 : const auto LoadValues =
1905 92778 : [this, dfX, dfY, &rgrids, nLowResMinIdxX, nLowResMinIdxY,
1906 : nCountLowResX, nCountLowResY, dfLowResMinX, dfLowResMinY,
1907 : dfLowResResX, dfLowResResY, &adfX, &adfY, &afDepth,
1908 92778 : &afUncrt](int iX, int iY)
1909 : {
1910 92778 : LoadClosestRefinedNodes(
1911 : dfX, dfY, iX, iY, rgrids, nLowResMinIdxX,
1912 : nLowResMinIdxY, nCountLowResX, nCountLowResY,
1913 : dfLowResMinX, dfLowResMinY, dfLowResResX, dfLowResResY,
1914 : adfX, adfY, afDepth, afUncrt);
1915 123332 : };
1916 :
1917 : // Load values of the closest point to the point to be
1918 : // interpolated in the current refinement grid
1919 30554 : LoadValues(iXRefinedGrid, iYRefinedGrid);
1920 :
1921 30554 : const bool bAtLeft = iXInRefinedGrid < 0 && iXRefinedGrid > 0;
1922 30554 : const bool bAtRight =
1923 59146 : iXInRefinedGrid >= static_cast<int>(rgrid.nWidth) - 1 &&
1924 28592 : iXRefinedGrid + 1 < poGDS->m_nLowResWidth;
1925 30554 : const bool bAtBottom = iYInRefinedGrid < 0 && iYRefinedGrid > 0;
1926 30554 : const bool bAtTop =
1927 59396 : iYInRefinedGrid >= static_cast<int>(rgrid.nHeight) - 1 &&
1928 28842 : iYRefinedGrid + 1 < poGDS->m_nLowResHeight;
1929 30554 : const int nXShift = bAtLeft ? -1 : bAtRight ? 1 : 0;
1930 30554 : const int nYShift = bAtBottom ? -1 : bAtTop ? 1 : 0;
1931 :
1932 : // Load values of the closest point to the point to be
1933 : // interpolated in the surrounding refinement grids.
1934 30554 : if (nXShift)
1935 : {
1936 23824 : LoadValues(iXRefinedGrid + nXShift, iYRefinedGrid);
1937 23824 : if (nYShift)
1938 : {
1939 16778 : LoadValues(iXRefinedGrid + nXShift,
1940 : iYRefinedGrid + nYShift);
1941 : }
1942 : }
1943 30554 : if (nYShift)
1944 : {
1945 21622 : LoadValues(iXRefinedGrid, iYRefinedGrid + nYShift);
1946 : }
1947 :
1948 : // Filter out candidate source points that are away from
1949 : // target point of more than dfMaxDistance
1950 30554 : size_t j = 0;
1951 129888 : for (size_t i = 0; i < adfX.size(); ++i)
1952 : {
1953 99334 : if (SQ(adfX[i] - dfX) + SQ(adfY[i] - dfY) <=
1954 99334 : dfMaxDistance * dfMaxDistance)
1955 : {
1956 28244 : adfX[j] = adfX[i];
1957 28244 : adfY[j] = adfY[i];
1958 28244 : afDepth[j] = afDepth[i];
1959 28244 : afUncrt[j] = afUncrt[i];
1960 28244 : ++j;
1961 : }
1962 : }
1963 30554 : adfX.resize(j);
1964 30554 : adfY.resize(j);
1965 30554 : afDepth.resize(j);
1966 30554 : afUncrt.resize(j);
1967 :
1968 : // Now interpolate the target point from the source values
1969 30554 : bool bTryIDW = false;
1970 30554 : if (adfX.size() >= 3)
1971 : {
1972 : // If there are at least 3 source nodes, sort them
1973 : // by increasing distance w.r.t the point to be interpolated
1974 319 : anIndices.clear();
1975 1571 : for (size_t i = 0; i < adfX.size(); ++i)
1976 1252 : anIndices.push_back(static_cast<int>(i));
1977 : // Sort nodes by increasing distance w.r.t (dfX, dfY)
1978 319 : std::sort(
1979 : anIndices.begin(), anIndices.end(),
1980 14240 : [&adfX, &adfY, dfX, dfY](int i1, int i2)
1981 : {
1982 1780 : return SQ(adfX[i1] - dfX) + SQ(adfY[i1] - dfY) <
1983 1780 : SQ(adfX[i2] - dfX) + SQ(adfY[i2] - dfY);
1984 : });
1985 : double adfXSorted[3];
1986 : double adfYSorted[3];
1987 : float afDepthSorted[3];
1988 : float afUncrtSorted[3];
1989 1276 : for (int i = 0; i < 3; ++i)
1990 : {
1991 957 : adfXSorted[i] = adfX[anIndices[i]];
1992 957 : adfYSorted[i] = adfY[anIndices[i]];
1993 957 : afDepthSorted[i] = afDepth[anIndices[i]];
1994 957 : afUncrtSorted[i] = afUncrt[anIndices[i]];
1995 : }
1996 : double dfCoord0;
1997 : double dfCoord1;
1998 : double dfCoord2;
1999 : // Perform barycentric interpolation with those 3
2000 : // points if they are all valid, and if the point to
2001 : // be interpolated falls into it.
2002 319 : if (afDepthSorted[0] != m_fNoDataValue &&
2003 230 : afDepthSorted[1] != m_fNoDataValue &&
2004 174 : afDepthSorted[2] != m_fNoDataValue)
2005 : {
2006 82 : if (BarycentricInterpolation(dfX, dfY, adfXSorted,
2007 : adfYSorted, dfCoord0,
2008 : dfCoord1, dfCoord2))
2009 : {
2010 : // Inside triangle
2011 77 : depthsPtr[y * nBlockXSize + x] =
2012 77 : static_cast<float>(dfCoord0 * afDepthSorted[0] +
2013 77 : dfCoord1 * afDepthSorted[1] +
2014 77 : dfCoord2 * afDepthSorted[2]);
2015 77 : uncrtPtr[y * nBlockXSize + x] =
2016 77 : static_cast<float>(dfCoord0 * afUncrtSorted[0] +
2017 77 : dfCoord1 * afUncrtSorted[1] +
2018 77 : dfCoord2 * afUncrtSorted[2]);
2019 : }
2020 : else
2021 : {
2022 : // Attempt inverse distance weighting in the
2023 : // cases where the point to be interpolated doesn't
2024 : // fall within the triangle formed by the 3
2025 : // closes points.
2026 5 : bTryIDW = true;
2027 : }
2028 : }
2029 : }
2030 30554 : if (bTryIDW)
2031 : {
2032 : // Do inverse distance weighting a a fallback.
2033 :
2034 5 : int nCountValid = 0;
2035 5 : double dfTotalDepth = 0;
2036 5 : double dfTotalUncrt = 0;
2037 5 : double dfTotalWeight = 0;
2038 : // Epsilon value to add to weights to avoid potential
2039 : // divergence to infinity if a source node is too close
2040 : // to the target point
2041 : const double EPS =
2042 5 : SQ(std::min(poGDS->m_gt[1], -poGDS->m_gt[5]) / 10);
2043 20 : for (size_t i = 0; i < adfX.size(); ++i)
2044 : {
2045 15 : if (afDepth[i] != m_fNoDataValue)
2046 : {
2047 15 : nCountValid++;
2048 : double dfSqrDistance =
2049 15 : SQ(adfX[i] - dfX) + SQ(adfY[i] - dfY) + EPS;
2050 15 : double dfWeight = 1. / dfSqrDistance;
2051 15 : dfTotalDepth += dfWeight * afDepth[i];
2052 15 : dfTotalUncrt += dfWeight * afUncrt[i];
2053 15 : dfTotalWeight += dfWeight;
2054 : }
2055 : }
2056 5 : if (nCountValid >= 3)
2057 : {
2058 5 : depthsPtr[y * nBlockXSize + x] =
2059 5 : static_cast<float>(dfTotalDepth / dfTotalWeight);
2060 5 : uncrtPtr[y * nBlockXSize + x] =
2061 5 : static_cast<float>(dfTotalUncrt / dfTotalWeight);
2062 : }
2063 : }
2064 : }
2065 : }
2066 : }
2067 5 : end:
2068 5 : if (poBlock != nullptr)
2069 : {
2070 5 : poBlock->DropLock();
2071 5 : poBlock = nullptr;
2072 : }
2073 :
2074 5 : return eErr;
2075 : }
2076 :
2077 : /************************************************************************/
2078 : /* LoadClosestRefinedNodes() */
2079 : /************************************************************************/
2080 :
2081 92778 : void BAGInterpolatedBand::LoadClosestRefinedNodes(
2082 : double dfX, double dfY, int iXRefinedGrid, int iYRefinedGrid,
2083 : const std::vector<BAGRefinementGrid> &rgrids, int nLowResMinIdxX,
2084 : int nLowResMinIdxY, int nCountLowResX, int nCountLowResY,
2085 : double dfLowResMinX, double dfLowResMinY, double dfLowResResX,
2086 : double dfLowResResY, std::vector<double> &adfX, std::vector<double> &adfY,
2087 : std::vector<float> &afDepth, std::vector<float> &afUncrt)
2088 : {
2089 92778 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
2090 92778 : CPLAssert(iXRefinedGrid >= nLowResMinIdxX);
2091 92778 : CPLAssert(iXRefinedGrid < nLowResMinIdxX + nCountLowResX);
2092 92778 : CPLAssert(iYRefinedGrid >= nLowResMinIdxY);
2093 92778 : CPLAssert(iYRefinedGrid < nLowResMinIdxY + nCountLowResY);
2094 92778 : CPL_IGNORE_RET_VAL(nCountLowResY);
2095 : const auto &rgrid =
2096 92778 : rgrids[(iYRefinedGrid - nLowResMinIdxY) * nCountLowResX +
2097 92778 : (iXRefinedGrid - nLowResMinIdxX)];
2098 92778 : if (rgrid.nWidth == 0)
2099 : {
2100 0 : return;
2101 : }
2102 92778 : const float gridRes = std::max(rgrid.fResX, rgrid.fResY);
2103 92778 : if (!(gridRes > poGDS->m_dfResFilterMin &&
2104 92778 : gridRes <= poGDS->m_dfResFilterMax))
2105 : {
2106 0 : return;
2107 : }
2108 :
2109 92778 : const double dfMinRefinedX =
2110 92778 : dfLowResMinX + iXRefinedGrid * dfLowResResX + rgrid.fSWX;
2111 92778 : const double dfMinRefinedY =
2112 92778 : dfLowResMinY + iYRefinedGrid * dfLowResResY + rgrid.fSWY;
2113 92778 : const int iXInRefinedGrid =
2114 92778 : static_cast<int>(floor((dfX - dfMinRefinedX) / rgrid.fResX));
2115 92778 : const int iYInRefinedGrid =
2116 92778 : static_cast<int>(floor((dfY - dfMinRefinedY) / rgrid.fResY));
2117 :
2118 99334 : const auto LoadValues = [poGDS, dfMinRefinedX, dfMinRefinedY, &rgrid, &adfX,
2119 : &adfY, &afDepth,
2120 496670 : &afUncrt](int iXAdjusted, int iYAdjusted)
2121 : {
2122 99334 : const unsigned nRefinementIndex =
2123 99334 : rgrid.nIndex + iYAdjusted * rgrid.nWidth + iXAdjusted;
2124 99334 : const auto pafRefValues = poGDS->GetRefinementValues(nRefinementIndex);
2125 99334 : if (pafRefValues)
2126 : {
2127 0 : adfX.push_back(dfMinRefinedX +
2128 99334 : iXAdjusted * static_cast<double>(rgrid.fResX));
2129 198668 : adfY.push_back(dfMinRefinedY +
2130 99334 : iYAdjusted * static_cast<double>(rgrid.fResY));
2131 99334 : afDepth.push_back(pafRefValues[0]);
2132 99334 : afUncrt.push_back(pafRefValues[1]);
2133 : }
2134 99334 : };
2135 :
2136 : const int iXAdjusted = std::max(
2137 92778 : 0, std::min(iXInRefinedGrid, static_cast<int>(rgrid.nWidth) - 1));
2138 : const int iYAdjusted = std::max(
2139 92778 : 0, std::min(iYInRefinedGrid, static_cast<int>(rgrid.nHeight) - 1));
2140 92778 : LoadValues(iXAdjusted, iYAdjusted);
2141 92778 : if (iYInRefinedGrid >= 0 &&
2142 54378 : iYInRefinedGrid < static_cast<int>(rgrid.nHeight) - 1)
2143 3145 : LoadValues(iXAdjusted, iYInRefinedGrid + 1);
2144 92778 : if (iXInRefinedGrid >= 0 &&
2145 52176 : iXInRefinedGrid < static_cast<int>(rgrid.nWidth) - 1)
2146 3411 : LoadValues(iXInRefinedGrid + 1, iYAdjusted);
2147 : }
2148 :
2149 : /************************************************************************/
2150 : /* ==================================================================== */
2151 : /* BAGGeorefMDBandBase */
2152 : /* ==================================================================== */
2153 : /************************************************************************/
2154 :
2155 : class BAGGeorefMDBandBase CPL_NON_FINAL : public GDALPamRasterBand
2156 : {
2157 : protected:
2158 : std::shared_ptr<GDALMDArray> m_poKeys;
2159 : std::unique_ptr<GDALRasterBand> m_poElevBand;
2160 : std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
2161 :
2162 6 : BAGGeorefMDBandBase(const std::shared_ptr<GDALMDArray> &poValues,
2163 : const std::shared_ptr<GDALMDArray> &poKeys,
2164 : GDALRasterBand *poElevBand)
2165 6 : : m_poKeys(poKeys), m_poElevBand(poElevBand),
2166 6 : m_poRAT(HDF5CreateRAT(poValues, false))
2167 : {
2168 6 : }
2169 :
2170 : CPLErr IReadBlockFromElevBand(int nBlockXOff, int nBlockYOff, void *pImage);
2171 :
2172 : public:
2173 1 : GDALRasterAttributeTable *GetDefaultRAT() override
2174 : {
2175 1 : return m_poRAT.get();
2176 : }
2177 :
2178 : double GetNoDataValue(int *pbSuccess) override;
2179 : };
2180 :
2181 : /************************************************************************/
2182 : /* GetNoDataValue() */
2183 : /************************************************************************/
2184 :
2185 1 : double BAGGeorefMDBandBase::GetNoDataValue(int *pbSuccess)
2186 : {
2187 1 : if (pbSuccess)
2188 1 : *pbSuccess = TRUE;
2189 1 : return 0;
2190 : }
2191 :
2192 : /************************************************************************/
2193 : /* IReadBlockFromElevBand() */
2194 : /************************************************************************/
2195 :
2196 8 : CPLErr BAGGeorefMDBandBase::IReadBlockFromElevBand(int nBlockXOff,
2197 : int nBlockYOff, void *pImage)
2198 : {
2199 16 : std::vector<float> afData(static_cast<size_t>(nBlockXSize) * nBlockYSize);
2200 8 : const int nXOff = nBlockXOff * nBlockXSize;
2201 8 : const int nReqXSize = std::min(nBlockXSize, nRasterXSize - nXOff);
2202 8 : const int nYOff = nBlockYOff * nBlockYSize;
2203 8 : const int nReqYSize = std::min(nBlockYSize, nRasterYSize - nYOff);
2204 16 : if (m_poElevBand->RasterIO(GF_Read, nXOff, nYOff, nReqXSize, nReqYSize,
2205 8 : &afData[0], nReqXSize, nReqYSize, GDT_Float32, 4,
2206 16 : nBlockXSize * 4, nullptr) != CE_None)
2207 : {
2208 0 : return CE_Failure;
2209 : }
2210 8 : int bHasNoData = FALSE;
2211 : const float fNoDataValue =
2212 8 : static_cast<float>(m_poElevBand->GetNoDataValue(&bHasNoData));
2213 8 : GByte *const pbyImage = static_cast<GByte *>(pImage);
2214 16 : for (int y = 0; y < nReqYSize; y++)
2215 : {
2216 40 : for (int x = 0; x < nReqXSize; x++)
2217 : {
2218 32 : pbyImage[y * nBlockXSize + x] =
2219 32 : (afData[y * nBlockXSize + x] == fNoDataValue ||
2220 29 : std::isnan(afData[y * nBlockXSize + x]))
2221 61 : ? 0
2222 : : 1;
2223 : }
2224 : }
2225 :
2226 8 : return CE_None;
2227 : }
2228 :
2229 : /************************************************************************/
2230 : /* ==================================================================== */
2231 : /* BAGGeorefMDBand */
2232 : /* ==================================================================== */
2233 : /************************************************************************/
2234 :
2235 : class BAGGeorefMDBand final : public BAGGeorefMDBandBase
2236 : {
2237 : public:
2238 : BAGGeorefMDBand(const std::shared_ptr<GDALMDArray> &poValues,
2239 : const std::shared_ptr<GDALMDArray> &poKeys,
2240 : GDALRasterBand *poElevBand);
2241 :
2242 : CPLErr IReadBlock(int, int, void *) override;
2243 : };
2244 :
2245 : /************************************************************************/
2246 : /* BAGGeorefMDBand() */
2247 : /************************************************************************/
2248 :
2249 2 : BAGGeorefMDBand::BAGGeorefMDBand(const std::shared_ptr<GDALMDArray> &poValues,
2250 : const std::shared_ptr<GDALMDArray> &poKeys,
2251 2 : GDALRasterBand *poElevBand)
2252 2 : : BAGGeorefMDBandBase(poValues, poKeys, poElevBand)
2253 : {
2254 2 : nRasterXSize = poElevBand->GetXSize();
2255 2 : nRasterYSize = poElevBand->GetYSize();
2256 2 : if (poKeys)
2257 : {
2258 2 : auto blockSize = poKeys->GetBlockSize();
2259 1 : CPLAssert(blockSize.size() == 2);
2260 1 : nBlockYSize = static_cast<int>(blockSize[0]);
2261 1 : nBlockXSize = static_cast<int>(blockSize[1]);
2262 1 : eDataType = poKeys->GetDataType().GetNumericDataType();
2263 1 : if (nBlockXSize == 0 || nBlockYSize == 0)
2264 : {
2265 0 : nBlockXSize = nRasterXSize;
2266 0 : nBlockYSize = 1;
2267 : }
2268 : }
2269 : else
2270 : {
2271 1 : eDataType = GDT_UInt8;
2272 1 : m_poElevBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
2273 : }
2274 :
2275 : // For testing purposes
2276 : const char *pszBlockXSize =
2277 2 : CPLGetConfigOption("BAG_GEOREF_MD_BLOCKXSIZE", nullptr);
2278 2 : if (pszBlockXSize)
2279 0 : nBlockXSize = atoi(pszBlockXSize);
2280 : const char *pszBlockYSize =
2281 2 : CPLGetConfigOption("BAG_GEOREF_MD_BLOCKYSIZE", nullptr);
2282 2 : if (pszBlockYSize)
2283 0 : nBlockYSize = atoi(pszBlockYSize);
2284 2 : }
2285 :
2286 : /************************************************************************/
2287 : /* IReadBlock() */
2288 : /************************************************************************/
2289 :
2290 8 : CPLErr BAGGeorefMDBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
2291 : {
2292 : HDF5_GLOBAL_LOCK();
2293 :
2294 8 : if (m_poKeys)
2295 : {
2296 : const GUInt64 arrayStartIdx[2] = {
2297 4 : static_cast<GUInt64>(
2298 4 : std::max(0, nRasterYSize - (nBlockYOff + 1) * nBlockYSize)),
2299 4 : static_cast<GUInt64>(nBlockXOff) * nBlockXSize};
2300 : size_t count[2] = {
2301 8 : std::min(static_cast<size_t>(nBlockYSize),
2302 4 : static_cast<size_t>(GetYSize() - arrayStartIdx[0])),
2303 8 : std::min(static_cast<size_t>(nBlockXSize),
2304 4 : static_cast<size_t>(GetXSize() - arrayStartIdx[1]))};
2305 4 : if (nRasterYSize - (nBlockYOff + 1) * nBlockYSize < 0)
2306 : {
2307 0 : count[0] += (nRasterYSize - (nBlockYOff + 1) * nBlockYSize);
2308 : }
2309 4 : const GInt64 arrayStep[2] = {1, 1};
2310 4 : const GPtrDiff_t bufferStride[2] = {nBlockXSize, 1};
2311 :
2312 8 : if (!m_poKeys->Read(arrayStartIdx, count, arrayStep, bufferStride,
2313 4 : m_poKeys->GetDataType(), pImage))
2314 : {
2315 0 : return CE_Failure;
2316 : }
2317 :
2318 : // Y flip the data.
2319 4 : const int nLinesToFlip = static_cast<int>(count[0]);
2320 4 : if (nLinesToFlip > 1)
2321 : {
2322 : const int nLineSize =
2323 4 : GDALGetDataTypeSizeBytes(eDataType) * nBlockXSize;
2324 4 : GByte *const pabyTemp = static_cast<GByte *>(CPLMalloc(nLineSize));
2325 4 : GByte *const pbyImage = static_cast<GByte *>(pImage);
2326 :
2327 8 : for (int iY = 0; iY < nLinesToFlip / 2; iY++)
2328 : {
2329 4 : memcpy(pabyTemp, pbyImage + iY * nLineSize, nLineSize);
2330 4 : memcpy(pbyImage + iY * nLineSize,
2331 4 : pbyImage + (nLinesToFlip - iY - 1) * nLineSize,
2332 : nLineSize);
2333 4 : memcpy(pbyImage + (nLinesToFlip - iY - 1) * nLineSize, pabyTemp,
2334 : nLineSize);
2335 : }
2336 :
2337 4 : CPLFree(pabyTemp);
2338 : }
2339 4 : return CE_None;
2340 : }
2341 : else
2342 : {
2343 4 : return IReadBlockFromElevBand(nBlockXOff, nBlockYOff, pImage);
2344 : }
2345 : }
2346 :
2347 : /************************************************************************/
2348 : /* ==================================================================== */
2349 : /* BAGGeorefMDSuperGridBand */
2350 : /* ==================================================================== */
2351 : /************************************************************************/
2352 :
2353 : class BAGGeorefMDSuperGridBand final : public BAGGeorefMDBandBase
2354 : {
2355 : public:
2356 : BAGGeorefMDSuperGridBand(const std::shared_ptr<GDALMDArray> &poValues,
2357 : const std::shared_ptr<GDALMDArray> &poKeys,
2358 : GDALRasterBand *poElevBand);
2359 :
2360 : CPLErr IReadBlock(int, int, void *) override;
2361 : };
2362 :
2363 : /************************************************************************/
2364 : /* BAGGeorefMDSuperGridBand() */
2365 : /************************************************************************/
2366 :
2367 4 : BAGGeorefMDSuperGridBand::BAGGeorefMDSuperGridBand(
2368 : const std::shared_ptr<GDALMDArray> &poValues,
2369 4 : const std::shared_ptr<GDALMDArray> &poKeys, GDALRasterBand *poElevBand)
2370 4 : : BAGGeorefMDBandBase(poValues, poKeys, poElevBand)
2371 : {
2372 4 : nRasterXSize = poElevBand->GetXSize();
2373 4 : nRasterYSize = poElevBand->GetYSize();
2374 4 : if (poKeys)
2375 : {
2376 2 : nBlockYSize = 1;
2377 2 : nBlockXSize = nRasterXSize;
2378 2 : eDataType = poKeys->GetDataType().GetNumericDataType();
2379 : }
2380 : else
2381 : {
2382 2 : eDataType = GDT_UInt8;
2383 2 : m_poElevBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
2384 : }
2385 4 : }
2386 :
2387 : /************************************************************************/
2388 : /* IReadBlock() */
2389 : /************************************************************************/
2390 :
2391 8 : CPLErr BAGGeorefMDSuperGridBand::IReadBlock(int nBlockXOff, int nBlockYOff,
2392 : void *pImage)
2393 : {
2394 8 : BAGDataset *poGDS = cpl::down_cast<BAGDataset *>(poDS);
2395 8 : if (m_poKeys)
2396 : {
2397 4 : const GUInt64 arrayStartIdx[2] = {
2398 4 : 0, poGDS->m_nSuperGridRefinementStartIndex +
2399 4 : static_cast<GUInt64>(nRasterYSize - 1 - nBlockYOff) *
2400 4 : nBlockXSize};
2401 4 : size_t count[2] = {1, static_cast<size_t>(nBlockXSize)};
2402 4 : const GInt64 arrayStep[2] = {1, 1};
2403 4 : const GPtrDiff_t bufferStride[2] = {nBlockXSize, 1};
2404 :
2405 8 : if (!m_poKeys->Read(arrayStartIdx, count, arrayStep, bufferStride,
2406 4 : m_poKeys->GetDataType(), pImage))
2407 : {
2408 0 : return CE_Failure;
2409 : }
2410 4 : return CE_None;
2411 : }
2412 : else
2413 : {
2414 4 : return IReadBlockFromElevBand(nBlockXOff, nBlockYOff, pImage);
2415 : }
2416 : }
2417 :
2418 : /************************************************************************/
2419 : /* ==================================================================== */
2420 : /* BAGDataset */
2421 : /* ==================================================================== */
2422 : /************************************************************************/
2423 :
2424 : /************************************************************************/
2425 : /* BAGDataset() */
2426 : /************************************************************************/
2427 :
2428 89 : BAGDataset::BAGDataset()
2429 : {
2430 89 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2431 89 : }
2432 :
2433 28 : BAGDataset::BAGDataset(BAGDataset *poParentDS, int nOvrFactor)
2434 : {
2435 28 : const int nXSize = poParentDS->nRasterXSize / nOvrFactor;
2436 28 : const int nYSize = poParentDS->nRasterYSize / nOvrFactor;
2437 28 : InitOverviewDS(poParentDS, nXSize, nYSize);
2438 28 : }
2439 :
2440 17 : BAGDataset::BAGDataset(BAGDataset *poParentDS, int nXSize, int nYSize)
2441 : {
2442 17 : InitOverviewDS(poParentDS, nXSize, nYSize);
2443 17 : }
2444 :
2445 45 : void BAGDataset::InitOverviewDS(BAGDataset *poParentDS, int nXSize, int nYSize)
2446 : {
2447 45 : m_ePopulation = poParentDS->m_ePopulation;
2448 45 : m_bMask = poParentDS->m_bMask;
2449 45 : m_bIsChild = true;
2450 : // m_apoOverviewDS
2451 45 : m_poSharedResources = poParentDS->m_poSharedResources;
2452 45 : m_poRootGroup = poParentDS->m_poRootGroup;
2453 45 : m_oSRS = poParentDS->m_oSRS;
2454 45 : nRasterXSize = nXSize;
2455 45 : nRasterYSize = nYSize;
2456 45 : m_gt[0] = poParentDS->m_gt[0];
2457 45 : m_gt[1] = poParentDS->m_gt[1] * poParentDS->nRasterXSize / nRasterXSize;
2458 45 : m_gt[2] = poParentDS->m_gt[2];
2459 45 : m_gt[3] = poParentDS->m_gt[3];
2460 45 : m_gt[4] = poParentDS->m_gt[4];
2461 45 : m_gt[5] = poParentDS->m_gt[5] * poParentDS->nRasterYSize / nRasterYSize;
2462 45 : m_nLowResWidth = poParentDS->m_nLowResWidth;
2463 45 : m_nLowResHeight = poParentDS->m_nLowResHeight;
2464 45 : m_dfLowResMinX = poParentDS->m_dfLowResMinX;
2465 45 : m_dfLowResMinY = poParentDS->m_dfLowResMinY;
2466 45 : m_dfLowResMaxX = poParentDS->m_dfLowResMaxX;
2467 45 : m_dfLowResMaxY = poParentDS->m_dfLowResMaxY;
2468 : // char *pszXMLMetadata = nullptr;
2469 : // char *apszMDList[2]{};
2470 45 : m_nChunkXSizeVarresMD = poParentDS->m_nChunkXSizeVarresMD;
2471 45 : m_nChunkYSizeVarresMD = poParentDS->m_nChunkYSizeVarresMD;
2472 45 : m_nChunkSizeVarresRefinement = poParentDS->m_nChunkSizeVarresRefinement;
2473 :
2474 45 : m_hVarresMetadata = poParentDS->m_hVarresMetadata;
2475 45 : m_hVarresMetadataDataType = poParentDS->m_hVarresMetadataDataType;
2476 45 : m_hVarresMetadataDataspace = poParentDS->m_hVarresMetadataDataspace;
2477 45 : m_hVarresMetadataNative = poParentDS->m_hVarresMetadataNative;
2478 : // m_oMapRefinemendGrids;
2479 :
2480 : // m_aosSubdatasets;
2481 :
2482 45 : m_hVarresRefinements = poParentDS->m_hVarresRefinements;
2483 45 : m_hVarresRefinementsDataType = poParentDS->m_hVarresRefinementsDataType;
2484 45 : m_hVarresRefinementsDataspace = poParentDS->m_hVarresRefinementsDataspace;
2485 45 : m_hVarresRefinementsNative = poParentDS->m_hVarresRefinementsNative;
2486 45 : m_nRefinementsSize = poParentDS->m_nRefinementsSize;
2487 :
2488 45 : m_nSuperGridRefinementStartIndex =
2489 45 : poParentDS->m_nSuperGridRefinementStartIndex;
2490 45 : m_dfResFilterMin = poParentDS->m_dfResFilterMin;
2491 45 : m_dfResFilterMax = poParentDS->m_dfResFilterMax;
2492 :
2493 45 : if (poParentDS->GetRasterCount() > 1)
2494 : {
2495 45 : GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
2496 : }
2497 45 : }
2498 :
2499 : /************************************************************************/
2500 : /* ~BAGDataset() */
2501 : /************************************************************************/
2502 268 : BAGDataset::~BAGDataset()
2503 : {
2504 134 : if (eAccess == GA_Update && nBands == 1)
2505 : {
2506 2 : auto poFirstBand = cpl::down_cast<BAGRasterBand *>(GetRasterBand(1));
2507 2 : auto poBand = new BAGRasterBand(this, 2);
2508 2 : poBand->nBlockXSize = poFirstBand->nBlockXSize;
2509 2 : poBand->nBlockYSize = poFirstBand->nBlockYSize;
2510 2 : poBand->eDataType = GDT_Float32;
2511 2 : poBand->m_bHasNoData = true;
2512 2 : poBand->m_fNoDataValue = poFirstBand->m_fNoDataValue;
2513 2 : SetBand(2, poBand);
2514 : }
2515 :
2516 134 : if (eAccess == GA_Update)
2517 : {
2518 12 : for (int i = 0; i < nBands; i++)
2519 : {
2520 : cpl::down_cast<BAGRasterBand *>(GetRasterBand(i + 1))
2521 8 : ->CreateDatasetIfNeeded();
2522 : }
2523 : }
2524 :
2525 134 : FlushCache(true);
2526 :
2527 134 : m_apoOverviewDS.clear();
2528 134 : if (!m_bIsChild)
2529 : {
2530 89 : if (m_hVarresMetadataDataType >= 0)
2531 63 : H5Tclose(m_hVarresMetadataDataType);
2532 :
2533 89 : if (m_hVarresMetadataDataspace >= 0)
2534 63 : H5Sclose(m_hVarresMetadataDataspace);
2535 :
2536 89 : if (m_hVarresMetadataNative >= 0)
2537 63 : H5Tclose(m_hVarresMetadataNative);
2538 :
2539 89 : if (m_hVarresMetadata >= 0)
2540 63 : H5Dclose(m_hVarresMetadata);
2541 :
2542 89 : if (m_hVarresRefinementsDataType >= 0)
2543 63 : H5Tclose(m_hVarresRefinementsDataType);
2544 :
2545 89 : if (m_hVarresRefinementsDataspace >= 0)
2546 63 : H5Sclose(m_hVarresRefinementsDataspace);
2547 :
2548 89 : if (m_hVarresRefinementsNative >= 0)
2549 63 : H5Tclose(m_hVarresRefinementsNative);
2550 :
2551 89 : if (m_hVarresRefinements >= 0)
2552 63 : H5Dclose(m_hVarresRefinements);
2553 :
2554 89 : CPLFree(pszXMLMetadata);
2555 : }
2556 268 : }
2557 :
2558 : /************************************************************************/
2559 : /* GH5DopenNoWarning() */
2560 : /************************************************************************/
2561 :
2562 83 : static hid_t GH5DopenNoWarning(hid_t hHDF5, const char *pszDatasetName)
2563 : {
2564 : hid_t hDataset;
2565 : #ifdef HAVE_GCC_WARNING_ZERO_AS_NULL_POINTER_CONSTANT
2566 : #pragma GCC diagnostic push
2567 : #pragma GCC diagnostic ignored "-Wzero-as-null-pointer-constant"
2568 : #endif
2569 83 : H5E_BEGIN_TRY
2570 : {
2571 83 : hDataset = H5Dopen(hHDF5, pszDatasetName);
2572 : }
2573 83 : H5E_END_TRY;
2574 :
2575 : #ifdef HAVE_GCC_WARNING_ZERO_AS_NULL_POINTER_CONSTANT
2576 : #pragma GCC diagnostic pop
2577 : #endif
2578 83 : return hDataset;
2579 : }
2580 :
2581 : /************************************************************************/
2582 : /* Open() */
2583 : /************************************************************************/
2584 :
2585 98 : GDALDataset *BAGDataset::Open(GDALOpenInfo *poOpenInfo)
2586 :
2587 : {
2588 : // Confirm that this appears to be a BAG file.
2589 98 : if (!BAGDatasetIdentify(poOpenInfo))
2590 12 : return nullptr;
2591 :
2592 : HDF5_GLOBAL_LOCK();
2593 :
2594 86 : if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
2595 : {
2596 0 : return HDF5Dataset::OpenMultiDim(poOpenInfo);
2597 : }
2598 :
2599 : // Confirm the requested access is supported.
2600 86 : if (poOpenInfo->eAccess == GA_Update)
2601 : {
2602 0 : ReportUpdateNotSupportedByDriver("BAG");
2603 0 : return nullptr;
2604 : }
2605 :
2606 86 : bool bOpenSuperGrid = false;
2607 86 : int nX = -1;
2608 86 : int nY = -1;
2609 172 : CPLString osFilename(poOpenInfo->pszFilename);
2610 172 : CPLString osGeorefMetadataLayer;
2611 86 : bool bIsSubdataset = false;
2612 86 : if (STARTS_WITH(poOpenInfo->pszFilename, "BAG:"))
2613 : {
2614 16 : bIsSubdataset = true;
2615 : char **papszTokens =
2616 16 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
2617 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES);
2618 :
2619 17 : if (CSLCount(papszTokens) == 3 &&
2620 1 : EQUAL(papszTokens[2], "bathymetry_coverage"))
2621 : {
2622 1 : osFilename = papszTokens[1];
2623 : }
2624 18 : else if (CSLCount(papszTokens) == 4 &&
2625 3 : EQUAL(papszTokens[2], "georef_metadata"))
2626 : {
2627 3 : osFilename = papszTokens[1];
2628 3 : osGeorefMetadataLayer = papszTokens[3];
2629 : }
2630 16 : else if (CSLCount(papszTokens) == 6 &&
2631 4 : EQUAL(papszTokens[2], "georef_metadata"))
2632 : {
2633 4 : osFilename = papszTokens[1];
2634 4 : osGeorefMetadataLayer = papszTokens[3];
2635 4 : bOpenSuperGrid = true;
2636 4 : nY = atoi(papszTokens[4]);
2637 4 : nX = atoi(papszTokens[5]);
2638 : }
2639 : else
2640 : {
2641 8 : if (CSLCount(papszTokens) != 5)
2642 : {
2643 0 : CSLDestroy(papszTokens);
2644 0 : return nullptr;
2645 : }
2646 8 : bOpenSuperGrid = true;
2647 8 : osFilename = papszTokens[1];
2648 8 : nY = atoi(papszTokens[3]);
2649 8 : nX = atoi(papszTokens[4]);
2650 : }
2651 16 : if (bOpenSuperGrid)
2652 : {
2653 :
2654 12 : if (CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINX") !=
2655 11 : nullptr ||
2656 11 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINY") !=
2657 11 : nullptr ||
2658 11 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXX") !=
2659 11 : nullptr ||
2660 11 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXY") !=
2661 23 : nullptr ||
2662 11 : CSLFetchNameValue(poOpenInfo->papszOpenOptions,
2663 : "SUPERGRIDS_INDICES") != nullptr)
2664 : {
2665 1 : CPLError(
2666 : CE_Warning, CPLE_AppDefined,
2667 : "Open options MINX/MINY/MAXX/MAXY/SUPERGRIDS_INDICES are "
2668 : "ignored when opening a supergrid");
2669 : }
2670 : }
2671 16 : CSLDestroy(papszTokens);
2672 : }
2673 :
2674 : // Open the file as an HDF5 file.
2675 86 : hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
2676 86 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
2677 86 : hid_t hHDF5 = H5Fopen(
2678 : osFilename,
2679 86 : /*poOpenInfo->eAccess == GA_Update ? H5F_ACC_RDWR : */ H5F_ACC_RDONLY,
2680 : fapl);
2681 86 : H5Pclose(fapl);
2682 86 : if (hHDF5 < 0)
2683 1 : return nullptr;
2684 :
2685 : // Confirm it is a BAG dataset by checking for the
2686 : // BAG_root/Bag Version attribute.
2687 85 : const hid_t hBagRoot = H5Gopen(hHDF5, "/BAG_root");
2688 : const hid_t hVersion =
2689 85 : hBagRoot >= 0 ? H5Aopen_name(hBagRoot, "Bag Version") : -1;
2690 :
2691 85 : if (hVersion < 0)
2692 : {
2693 0 : if (hBagRoot >= 0)
2694 0 : H5Gclose(hBagRoot);
2695 0 : H5Fclose(hHDF5);
2696 0 : return nullptr;
2697 : }
2698 85 : H5Aclose(hVersion);
2699 :
2700 170 : auto poSharedResources = GDAL::HDF5SharedResources::Create(osFilename);
2701 85 : poSharedResources->m_hHDF5 = hHDF5;
2702 :
2703 170 : auto poRootGroup = HDF5Dataset::OpenGroup(poSharedResources);
2704 85 : if (poRootGroup == nullptr)
2705 0 : return nullptr;
2706 :
2707 : // Create a corresponding dataset.
2708 85 : BAGDataset *const poDS = new BAGDataset();
2709 :
2710 85 : poDS->eAccess = poOpenInfo->eAccess;
2711 85 : poDS->m_poRootGroup = std::move(poRootGroup);
2712 85 : poDS->m_poSharedResources = std::move(poSharedResources);
2713 :
2714 : // Extract version as metadata.
2715 170 : CPLString osVersion;
2716 :
2717 85 : if (GH5_FetchAttribute(hBagRoot, "Bag Version", osVersion))
2718 85 : poDS->GDALDataset::SetMetadataItem("BagVersion", osVersion);
2719 :
2720 85 : H5Gclose(hBagRoot);
2721 :
2722 170 : CPLString osSubDsName;
2723 85 : if (poOpenInfo->nOpenFlags & GDAL_OF_RASTER)
2724 : {
2725 84 : if (poDS->OpenRaster(poOpenInfo, osFilename, bOpenSuperGrid, nX, nY,
2726 : bIsSubdataset, osGeorefMetadataLayer, osSubDsName))
2727 : {
2728 76 : if (!osSubDsName.empty())
2729 : {
2730 2 : delete poDS;
2731 4 : GDALOpenInfo oOpenInfo(osSubDsName, GA_ReadOnly);
2732 2 : oOpenInfo.nOpenFlags = poOpenInfo->nOpenFlags;
2733 2 : return Open(&oOpenInfo);
2734 : }
2735 : }
2736 : else
2737 : {
2738 8 : delete poDS;
2739 8 : return nullptr;
2740 : }
2741 : }
2742 :
2743 75 : if (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR)
2744 : {
2745 47 : if (!poDS->OpenVector() &&
2746 1 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0)
2747 : {
2748 0 : delete poDS;
2749 0 : return nullptr;
2750 : }
2751 : }
2752 :
2753 75 : return poDS;
2754 : }
2755 :
2756 : /************************************************************************/
2757 : /* OpenRaster() */
2758 : /************************************************************************/
2759 :
2760 84 : bool BAGDataset::OpenRaster(GDALOpenInfo *poOpenInfo,
2761 : const CPLString &osFilename, bool bOpenSuperGrid,
2762 : int nX, int nY, bool bIsSubdataset,
2763 : const CPLString &osGeorefMetadataLayer,
2764 : CPLString &outOsSubDsName)
2765 : {
2766 : const char *pszMode =
2767 84 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "MODE", "AUTO");
2768 84 : const bool bLowResGrid = EQUAL(pszMode, "LOW_RES_GRID");
2769 85 : if (bLowResGrid &&
2770 1 : (CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINX") != nullptr ||
2771 0 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINY") != nullptr ||
2772 0 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXX") != nullptr ||
2773 0 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXY") != nullptr ||
2774 0 : CSLFetchNameValue(poOpenInfo->papszOpenOptions,
2775 : "SUPERGRIDS_INDICES") != nullptr))
2776 : {
2777 1 : CPLError(CE_Warning, CPLE_AppDefined,
2778 : "Open options MINX/MINY/MAXX/MAXY/SUPERGRIDS_INDICES are "
2779 : "ignored when opening the low resolution grid");
2780 : }
2781 :
2782 : const bool bListSubDS =
2783 167 : !bLowResGrid &&
2784 83 : (EQUAL(pszMode, "LIST_SUPERGRIDS") ||
2785 78 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINX") != nullptr ||
2786 76 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINY") != nullptr ||
2787 75 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXX") != nullptr ||
2788 75 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXY") != nullptr ||
2789 75 : CSLFetchNameValue(poOpenInfo->papszOpenOptions,
2790 84 : "SUPERGRIDS_INDICES") != nullptr);
2791 84 : const bool bResampledGrid = EQUAL(pszMode, "RESAMPLED_GRID");
2792 84 : const bool bInterpolate = EQUAL(pszMode, "INTERPOLATED");
2793 :
2794 : const char *pszNoDataValue =
2795 84 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "NODATA_VALUE");
2796 84 : bool bHasNoData = pszNoDataValue != nullptr;
2797 : float fNoDataValue =
2798 84 : bHasNoData ? static_cast<float>(CPLAtof(pszNoDataValue)) : 0.0f;
2799 :
2800 : // Fetch the elevation dataset and attach as a band.
2801 84 : int nNextBand = 1;
2802 :
2803 84 : BAGRasterBand *poElevBand = new BAGRasterBand(this, nNextBand);
2804 :
2805 : {
2806 : const hid_t hElevation =
2807 84 : H5Dopen(GetHDF5Handle(), "/BAG_root/elevation");
2808 84 : if (hElevation < 0 || !poElevBand->Initialize(hElevation, "elevation"))
2809 : {
2810 2 : delete poElevBand;
2811 2 : return false;
2812 : }
2813 : }
2814 :
2815 82 : m_nLowResWidth = poElevBand->nRasterXSize;
2816 82 : m_nLowResHeight = poElevBand->nRasterYSize;
2817 :
2818 82 : if (bOpenSuperGrid || bListSubDS || bResampledGrid || bInterpolate)
2819 : {
2820 53 : if (!bHasNoData)
2821 : {
2822 53 : int nHasNoData = FALSE;
2823 53 : double dfNoData = poElevBand->GetNoDataValue(&nHasNoData);
2824 53 : if (nHasNoData)
2825 : {
2826 53 : bHasNoData = true;
2827 53 : fNoDataValue = static_cast<float>(dfNoData);
2828 : }
2829 : }
2830 53 : delete poElevBand;
2831 53 : nRasterXSize = 0;
2832 53 : nRasterYSize = 0;
2833 : }
2834 29 : else if (!osGeorefMetadataLayer.empty())
2835 : {
2836 3 : auto poGeoref_metadataLayer = m_poRootGroup->OpenGroupFromFullname(
2837 3 : "/BAG_root/georef_metadata/" + osGeorefMetadataLayer, nullptr);
2838 3 : if (poGeoref_metadataLayer == nullptr)
2839 : {
2840 1 : CPLError(CE_Failure, CPLE_AppDefined,
2841 : "Cannot find georef_metadata layer %s",
2842 : osGeorefMetadataLayer.c_str());
2843 1 : delete poElevBand;
2844 1 : return false;
2845 : }
2846 :
2847 4 : auto poKeys = poGeoref_metadataLayer->OpenMDArray("keys");
2848 2 : if (poKeys != nullptr)
2849 : {
2850 1 : auto poDims = poKeys->GetDimensions();
2851 1 : if (poDims.size() != 2 ||
2852 1 : poDims[0]->GetSize() !=
2853 3 : static_cast<size_t>(poElevBand->nRasterYSize) ||
2854 1 : poDims[1]->GetSize() !=
2855 1 : static_cast<size_t>(poElevBand->nRasterXSize))
2856 : {
2857 0 : CPLError(CE_Failure, CPLE_AppDefined,
2858 : "Wrong dimensions for %s/keys",
2859 : osGeorefMetadataLayer.c_str());
2860 0 : delete poElevBand;
2861 0 : return false;
2862 : }
2863 2 : if (poKeys->GetDataType().GetClass() != GEDTC_NUMERIC ||
2864 1 : !GDALDataTypeIsInteger(
2865 1 : poKeys->GetDataType().GetNumericDataType()))
2866 : {
2867 0 : CPLError(CE_Failure, CPLE_AppDefined,
2868 : "Only integer data type supported for %s/keys",
2869 : osGeorefMetadataLayer.c_str());
2870 0 : delete poElevBand;
2871 0 : return false;
2872 : }
2873 : }
2874 :
2875 4 : auto poValues = poGeoref_metadataLayer->OpenMDArray("values");
2876 2 : if (poValues == nullptr)
2877 : {
2878 0 : CPLError(CE_Failure, CPLE_AppDefined,
2879 : "Cannot find array values of georef_metadata layer %s",
2880 : osGeorefMetadataLayer.c_str());
2881 0 : delete poElevBand;
2882 0 : return false;
2883 : }
2884 2 : if (poValues->GetDimensionCount() != 1)
2885 : {
2886 0 : CPLError(CE_Failure, CPLE_AppDefined,
2887 : "Wrong dimensions for %s/values",
2888 : osGeorefMetadataLayer.c_str());
2889 0 : delete poElevBand;
2890 0 : return false;
2891 : }
2892 2 : if (poValues->GetDataType().GetClass() != GEDTC_COMPOUND)
2893 : {
2894 0 : CPLError(CE_Failure, CPLE_AppDefined,
2895 : "Only compound data type supported for %s/values",
2896 : osGeorefMetadataLayer.c_str());
2897 0 : delete poElevBand;
2898 0 : return false;
2899 : }
2900 :
2901 2 : nRasterXSize = poElevBand->nRasterXSize;
2902 2 : nRasterYSize = poElevBand->nRasterYSize;
2903 2 : SetBand(1, new BAGGeorefMDBand(poValues, poKeys, poElevBand));
2904 : }
2905 : else
2906 : {
2907 26 : nRasterXSize = poElevBand->nRasterXSize;
2908 26 : nRasterYSize = poElevBand->nRasterYSize;
2909 :
2910 26 : SetBand(nNextBand++, poElevBand);
2911 :
2912 : // Try to do the same for the uncertainty band.
2913 : const hid_t hUncertainty =
2914 26 : H5Dopen(GetHDF5Handle(), "/BAG_root/uncertainty");
2915 26 : BAGRasterBand *poUBand = new BAGRasterBand(this, nNextBand);
2916 :
2917 52 : if (hUncertainty >= 0 &&
2918 26 : poUBand->Initialize(hUncertainty, "uncertainty"))
2919 : {
2920 26 : SetBand(nNextBand++, poUBand);
2921 : }
2922 : else
2923 : {
2924 0 : delete poUBand;
2925 : }
2926 :
2927 : // Load other root datasets (such as nominal_elevation)
2928 78 : auto poBAG_root = m_poRootGroup->OpenGroup("BAG_root", nullptr);
2929 26 : if (poBAG_root)
2930 : {
2931 52 : const auto arrayNames = poBAG_root->GetMDArrayNames(nullptr);
2932 153 : for (const auto &arrayName : arrayNames)
2933 : {
2934 127 : if (arrayName != "elevation" && arrayName != "uncertainty")
2935 : {
2936 150 : auto poArray = poBAG_root->OpenMDArray(arrayName, nullptr);
2937 150 : if (poArray && poArray->GetDimensions().size() == 2 &&
2938 18 : poArray->GetDimensions()[0]->GetSize() ==
2939 18 : static_cast<unsigned>(nRasterYSize) &&
2940 10 : poArray->GetDimensions()[1]->GetSize() ==
2941 160 : static_cast<unsigned>(nRasterXSize) &&
2942 10 : poArray->GetDataType().GetClass() == GEDTC_NUMERIC)
2943 : {
2944 2 : hid_t hBandId = GH5DopenNoWarning(
2945 : GetHDF5Handle(),
2946 4 : ("/BAG_root/" + arrayName).c_str());
2947 : BAGRasterBand *const poBand =
2948 2 : new BAGRasterBand(this, nNextBand);
2949 4 : if (hBandId >= 0 &&
2950 2 : poBand->Initialize(hBandId, arrayName.c_str()))
2951 : {
2952 2 : SetBand(nNextBand++, poBand);
2953 : }
2954 : else
2955 : {
2956 0 : delete poBand;
2957 : }
2958 : }
2959 : }
2960 : }
2961 : }
2962 : }
2963 :
2964 81 : SetDescription(poOpenInfo->pszFilename);
2965 :
2966 81 : m_bReportVertCRS = CPLTestBool(CSLFetchNameValueDef(
2967 81 : poOpenInfo->papszOpenOptions, "REPORT_VERTCRS", "YES"));
2968 :
2969 : // Load the XML metadata.
2970 81 : LoadMetadata();
2971 :
2972 81 : if (bResampledGrid)
2973 : {
2974 25 : m_bMask = CPLTestBool(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
2975 : "SUPERGRIDS_MASK", "NO"));
2976 : }
2977 :
2978 81 : if (!m_bMask)
2979 : {
2980 80 : GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
2981 : }
2982 :
2983 : // Load for refinements grids for variable resolution datasets
2984 81 : bool bHasRefinementGrids = false;
2985 81 : if (bOpenSuperGrid || bListSubDS || bResampledGrid || bInterpolate)
2986 : {
2987 : bHasRefinementGrids =
2988 53 : LookForRefinementGrids(poOpenInfo->papszOpenOptions, nY, nX);
2989 55 : if (!bOpenSuperGrid && m_aosSubdatasets.size() == 2 &&
2990 2 : EQUAL(pszMode, "AUTO"))
2991 : {
2992 : outOsSubDsName =
2993 2 : CSLFetchNameValueDef(m_aosSubdatasets, "SUBDATASET_1_NAME", "");
2994 2 : return true;
2995 : }
2996 : }
2997 : else
2998 : {
2999 28 : if (LookForRefinementGrids(poOpenInfo->papszOpenOptions, 0, 0))
3000 : {
3001 10 : GDALDataset::SetMetadataItem("HAS_SUPERGRIDS", "TRUE");
3002 : }
3003 28 : m_aosSubdatasets.Clear();
3004 : }
3005 :
3006 79 : if (!bIsSubdataset && osGeorefMetadataLayer.empty())
3007 : {
3008 65 : auto poGeoref_metadata = m_poRootGroup->OpenGroupFromFullname(
3009 195 : "/BAG_root/georef_metadata", nullptr);
3010 65 : if (poGeoref_metadata)
3011 : {
3012 4 : const auto groupNames = poGeoref_metadata->GetGroupNames(nullptr);
3013 2 : if (!groupNames.empty())
3014 : {
3015 2 : if (m_aosSubdatasets.empty())
3016 : {
3017 1 : const int nIdx = 1;
3018 : m_aosSubdatasets.AddNameValue(
3019 : CPLSPrintf("SUBDATASET_%d_NAME", nIdx),
3020 : CPLSPrintf("BAG:\"%s\":bathymetry_coverage",
3021 1 : GetDescription()));
3022 : m_aosSubdatasets.AddNameValue(
3023 : CPLSPrintf("SUBDATASET_%d_DESC", nIdx),
3024 1 : "Bathymetry gridded data");
3025 : }
3026 6 : for (const auto &groupName : groupNames)
3027 : {
3028 4 : const int nIdx = m_aosSubdatasets.size() / 2 + 1;
3029 : m_aosSubdatasets.AddNameValue(
3030 : CPLSPrintf("SUBDATASET_%d_NAME", nIdx),
3031 : CPLSPrintf("BAG:\"%s\":georef_metadata:%s",
3032 4 : GetDescription(), groupName.c_str()));
3033 : m_aosSubdatasets.AddNameValue(
3034 : CPLSPrintf("SUBDATASET_%d_DESC", nIdx),
3035 : CPLSPrintf("Georeferenced metadata %s",
3036 4 : groupName.c_str()));
3037 : }
3038 : }
3039 : }
3040 : }
3041 :
3042 79 : double dfMinResX = 0.0;
3043 79 : double dfMinResY = 0.0;
3044 79 : double dfMaxResX = 0.0;
3045 79 : double dfMaxResY = 0.0;
3046 79 : if (m_hVarresMetadata >= 0)
3047 : {
3048 61 : if (!GH5_FetchAttribute(m_hVarresMetadata, "min_resolution_x",
3049 122 : dfMinResX) ||
3050 61 : !GH5_FetchAttribute(m_hVarresMetadata, "min_resolution_y",
3051 : dfMinResY))
3052 : {
3053 0 : CPLError(CE_Failure, CPLE_AppDefined,
3054 : "Cannot get min_resolution_x and/or min_resolution_y");
3055 0 : return false;
3056 : }
3057 :
3058 61 : if (!GH5_FetchAttribute(m_hVarresMetadata, "max_resolution_x",
3059 122 : dfMaxResX) ||
3060 61 : !GH5_FetchAttribute(m_hVarresMetadata, "max_resolution_y",
3061 : dfMaxResY))
3062 : {
3063 0 : CPLError(CE_Failure, CPLE_AppDefined,
3064 : "Cannot get max_resolution_x and/or max_resolution_y");
3065 0 : return false;
3066 : }
3067 :
3068 61 : if (!bOpenSuperGrid && !bResampledGrid && !bInterpolate)
3069 : {
3070 24 : GDALDataset::SetMetadataItem("MIN_RESOLUTION_X",
3071 : CPLSPrintf("%f", dfMinResX));
3072 24 : GDALDataset::SetMetadataItem("MIN_RESOLUTION_Y",
3073 : CPLSPrintf("%f", dfMinResY));
3074 24 : GDALDataset::SetMetadataItem("MAX_RESOLUTION_X",
3075 : CPLSPrintf("%f", dfMaxResX));
3076 24 : GDALDataset::SetMetadataItem("MAX_RESOLUTION_Y",
3077 : CPLSPrintf("%f", dfMaxResY));
3078 : }
3079 : }
3080 :
3081 79 : if (bResampledGrid || bInterpolate)
3082 : {
3083 26 : if (!bHasRefinementGrids)
3084 : {
3085 0 : CPLError(CE_Failure, CPLE_AppDefined,
3086 : "No supergrids available. "
3087 : "%s mode not available",
3088 : pszMode);
3089 3 : return false;
3090 : }
3091 :
3092 52 : const char *pszValuePopStrategy = CSLFetchNameValueDef(
3093 26 : poOpenInfo->papszOpenOptions, "VALUE_POPULATION", "MAX");
3094 26 : if (EQUAL(pszValuePopStrategy, "MIN"))
3095 : {
3096 1 : m_ePopulation = BAGDataset::Population::MIN;
3097 : }
3098 25 : else if (EQUAL(pszValuePopStrategy, "MEAN"))
3099 : {
3100 1 : m_ePopulation = BAGDataset::Population::MEAN;
3101 : }
3102 24 : else if (EQUAL(pszValuePopStrategy, "MAX"))
3103 : {
3104 23 : m_ePopulation = BAGDataset::Population::MAX;
3105 : }
3106 : else
3107 : {
3108 1 : m_ePopulation = BAGDataset::Population::COUNT;
3109 1 : bHasNoData = false;
3110 1 : fNoDataValue = 0;
3111 : }
3112 :
3113 : const char *pszResX =
3114 26 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "RESX");
3115 : const char *pszResY =
3116 26 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "RESY");
3117 52 : const char *pszResStrategy = CSLFetchNameValueDef(
3118 26 : poOpenInfo->papszOpenOptions, "RES_STRATEGY", "AUTO");
3119 26 : double dfDefaultResX = 0.0;
3120 26 : double dfDefaultResY = 0.0;
3121 :
3122 : const char *pszResFilterMin =
3123 26 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "RES_FILTER_MIN");
3124 : const char *pszResFilterMax =
3125 26 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "RES_FILTER_MAX");
3126 :
3127 26 : double dfResFilterMin = 0;
3128 26 : if (pszResFilterMin != nullptr)
3129 : {
3130 6 : dfResFilterMin = CPLAtof(pszResFilterMin);
3131 6 : const double dfMaxRes = std::min(dfMaxResX, dfMaxResY);
3132 6 : if (dfResFilterMin >= dfMaxRes)
3133 : {
3134 1 : CPLError(CE_Failure, CPLE_AppDefined,
3135 : "Cannot specified RES_FILTER_MIN >= %g", dfMaxRes);
3136 1 : return false;
3137 : }
3138 5 : GDALDataset::SetMetadataItem("RES_FILTER_MIN",
3139 : CPLSPrintf("%g", dfResFilterMin));
3140 : }
3141 :
3142 25 : double dfResFilterMax = std::numeric_limits<double>::infinity();
3143 25 : if (pszResFilterMax != nullptr)
3144 : {
3145 6 : dfResFilterMax = CPLAtof(pszResFilterMax);
3146 6 : const double dfMinRes = std::min(dfMinResX, dfMinResY);
3147 6 : if (dfResFilterMax < dfMinRes)
3148 : {
3149 2 : CPLError(CE_Failure, CPLE_AppDefined,
3150 : "Cannot specified RES_FILTER_MAX < %g", dfMinRes);
3151 2 : return false;
3152 : }
3153 4 : GDALDataset::SetMetadataItem("RES_FILTER_MAX",
3154 : CPLSPrintf("%g", dfResFilterMax));
3155 : }
3156 :
3157 23 : if (dfResFilterMin >= dfResFilterMax)
3158 : {
3159 0 : CPLError(CE_Failure, CPLE_AppDefined,
3160 : "Cannot specified RES_FILTER_MIN >= RES_FILTER_MAX");
3161 0 : return false;
3162 : }
3163 :
3164 23 : if (EQUAL(pszResStrategy, "AUTO") &&
3165 12 : (pszResFilterMin || pszResFilterMax))
3166 : {
3167 5 : if (pszResFilterMax)
3168 : {
3169 4 : dfDefaultResX = dfResFilterMax;
3170 4 : dfDefaultResY = dfResFilterMax;
3171 : }
3172 : else
3173 : {
3174 1 : dfDefaultResX = dfMaxResX;
3175 1 : dfDefaultResY = dfMaxResY;
3176 : }
3177 : }
3178 18 : else if (EQUAL(pszResStrategy, "AUTO") || EQUAL(pszResStrategy, "MIN"))
3179 : {
3180 12 : dfDefaultResX = dfMinResX;
3181 12 : dfDefaultResY = dfMinResY;
3182 : }
3183 6 : else if (EQUAL(pszResStrategy, "MAX"))
3184 : {
3185 1 : dfDefaultResX = dfMaxResX;
3186 1 : dfDefaultResY = dfMaxResY;
3187 : }
3188 5 : else if (EQUAL(pszResStrategy, "MEAN"))
3189 : {
3190 5 : if (!GetMeanSupergridsResolution(dfDefaultResX, dfDefaultResY))
3191 : {
3192 0 : return false;
3193 : }
3194 : }
3195 :
3196 : const char *pszMinX =
3197 23 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINX");
3198 : const char *pszMinY =
3199 23 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MINY");
3200 : const char *pszMaxX =
3201 23 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXX");
3202 : const char *pszMaxY =
3203 23 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAXY");
3204 :
3205 23 : double dfMinX = m_dfLowResMinX;
3206 23 : double dfMinY = m_dfLowResMinY;
3207 23 : double dfMaxX = m_dfLowResMaxX;
3208 23 : double dfMaxY = m_dfLowResMaxY;
3209 23 : double dfResX = dfDefaultResX;
3210 23 : double dfResY = dfDefaultResY;
3211 23 : if (pszMinX)
3212 1 : dfMinX = CPLAtof(pszMinX);
3213 23 : if (pszMinY)
3214 1 : dfMinY = CPLAtof(pszMinY);
3215 23 : if (pszMaxX)
3216 1 : dfMaxX = CPLAtof(pszMaxX);
3217 23 : if (pszMaxY)
3218 1 : dfMaxY = CPLAtof(pszMaxY);
3219 23 : if (pszResX)
3220 1 : dfResX = CPLAtof(pszResX);
3221 23 : if (pszResY)
3222 1 : dfResY = CPLAtof(pszResY);
3223 :
3224 23 : if (dfResX <= 0.0 || dfResY <= 0.0)
3225 : {
3226 0 : CPLError(CE_Failure, CPLE_NotSupported,
3227 : "Invalid resolution: %f x %f", dfResX, dfResY);
3228 0 : return false;
3229 : }
3230 23 : const double dfRasterXSize = (dfMaxX - dfMinX) / dfResX;
3231 23 : const double dfRasterYSize = (dfMaxY - dfMinY) / dfResY;
3232 23 : if (dfRasterXSize <= 1 || dfRasterYSize <= 1 ||
3233 23 : dfRasterXSize > INT_MAX || dfRasterYSize > INT_MAX)
3234 : {
3235 0 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid raster dimension");
3236 0 : return false;
3237 : }
3238 23 : nRasterXSize = static_cast<int>(dfRasterXSize + 0.5);
3239 23 : nRasterYSize = static_cast<int>(dfRasterYSize + 0.5);
3240 23 : m_gt[0] = dfMinX;
3241 23 : m_gt[1] = dfResX;
3242 23 : m_gt[3] = dfMaxY;
3243 23 : m_gt[5] = -dfResY;
3244 23 : if (pszMaxY == nullptr || pszMinY != nullptr)
3245 : {
3246 : // if the constraint is not given by MAXY, we may need to tweak
3247 : // m_gt[3] / maxy, so that we get the requested MINY
3248 : // value
3249 22 : m_gt[3] += dfMinY - (dfMaxY - nRasterYSize * dfResY);
3250 : }
3251 :
3252 23 : const double dfMinRes = std::min(dfMinResX, dfMinResY);
3253 23 : if (dfResFilterMin > dfMinRes)
3254 : {
3255 3 : m_dfResFilterMin = dfResFilterMin;
3256 : }
3257 23 : m_dfResFilterMax = dfResFilterMax;
3258 :
3259 : // Use min/max BAG refinement metadata items only if the
3260 : // GDAL dataset bounding box is equal or larger to the BAG dataset
3261 23 : const bool bInitializeMinMax =
3262 22 : (!m_bMask && m_ePopulation != BAGDataset::Population::COUNT &&
3263 21 : dfMinX <= m_dfLowResMinX && dfMinY <= m_dfLowResMinY &&
3264 45 : dfMaxX >= m_dfLowResMaxX && dfMaxY >= m_dfLowResMaxY);
3265 :
3266 23 : if (bInterpolate)
3267 : {
3268 : // Depth
3269 1 : SetBand(1,
3270 : new BAGInterpolatedBand(this, 1, bHasNoData, fNoDataValue,
3271 1 : bInitializeMinMax));
3272 :
3273 : // Uncertainty
3274 1 : SetBand(2,
3275 : new BAGInterpolatedBand(this, 2, bHasNoData, fNoDataValue,
3276 1 : bInitializeMinMax));
3277 : }
3278 22 : else if (m_bMask || m_ePopulation == BAGDataset::Population::COUNT)
3279 : {
3280 2 : SetBand(1, new BAGResampledBand(this, 1, false, 0.0f, false));
3281 : }
3282 : else
3283 : {
3284 20 : SetBand(1, new BAGResampledBand(this, 1, bHasNoData, fNoDataValue,
3285 20 : bInitializeMinMax));
3286 :
3287 20 : SetBand(2, new BAGResampledBand(this, 2, bHasNoData, fNoDataValue,
3288 20 : bInitializeMinMax));
3289 : }
3290 :
3291 23 : if (GetRasterCount() > 1)
3292 : {
3293 21 : GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
3294 : "IMAGE_STRUCTURE");
3295 : }
3296 :
3297 42 : const bool bCanUseLowResAsOvr = nRasterXSize > m_nLowResWidth &&
3298 42 : nRasterYSize > m_nLowResHeight &&
3299 19 : GetRasterCount() > 1;
3300 : const int nMinOvrSize =
3301 23 : bCanUseLowResAsOvr ? 1 + std::min(m_nLowResWidth, m_nLowResHeight)
3302 23 : : 256;
3303 51 : for (int nOvrFactor = 2; nRasterXSize / nOvrFactor >= nMinOvrSize &&
3304 32 : nRasterYSize / nOvrFactor >= nMinOvrSize;
3305 28 : nOvrFactor *= 2)
3306 : {
3307 56 : auto poOvrDS = std::make_unique<BAGDataset>(this, nOvrFactor);
3308 :
3309 84 : for (int i = 1; i <= GetRasterCount(); i++)
3310 : {
3311 56 : if (bInterpolate)
3312 : {
3313 16 : poOvrDS->SetBand(
3314 8 : i, new BAGInterpolatedBand(poOvrDS.get(), i, bHasNoData,
3315 8 : fNoDataValue, false));
3316 : }
3317 : else
3318 : {
3319 96 : poOvrDS->SetBand(
3320 48 : i, new BAGResampledBand(poOvrDS.get(), i, bHasNoData,
3321 48 : fNoDataValue, false));
3322 : }
3323 : }
3324 :
3325 28 : m_apoOverviewDS.emplace_back(std::move(poOvrDS));
3326 : }
3327 :
3328 : // Use the low resolution grid as the last overview level
3329 23 : if (bCanUseLowResAsOvr)
3330 : {
3331 17 : auto poOvrDS = std::make_unique<BAGDataset>(this, m_nLowResWidth,
3332 17 : m_nLowResHeight);
3333 :
3334 17 : poElevBand = new BAGRasterBand(poOvrDS.get(), 1);
3335 : const hid_t hElevation =
3336 17 : H5Dopen(GetHDF5Handle(), "/BAG_root/elevation");
3337 34 : if (hElevation < 0 ||
3338 17 : !poElevBand->Initialize(hElevation, "elevation"))
3339 : {
3340 0 : delete poElevBand;
3341 0 : return false;
3342 : }
3343 17 : poOvrDS->SetBand(1, poElevBand);
3344 :
3345 : const hid_t hUncertainty =
3346 17 : H5Dopen(GetHDF5Handle(), "/BAG_root/uncertainty");
3347 17 : BAGRasterBand *poUBand = new BAGRasterBand(poOvrDS.get(), 2);
3348 34 : if (hUncertainty < 0 ||
3349 17 : !poUBand->Initialize(hUncertainty, "uncertainty"))
3350 : {
3351 0 : delete poUBand;
3352 0 : return false;
3353 : }
3354 17 : poOvrDS->SetBand(2, poUBand);
3355 :
3356 17 : m_apoOverviewDS.emplace_back(std::move(poOvrDS));
3357 23 : }
3358 : }
3359 53 : else if (bOpenSuperGrid)
3360 : {
3361 11 : auto oIter = m_oMapRefinemendGrids.find(nY * m_nLowResWidth + nX);
3362 11 : if (oIter == m_oMapRefinemendGrids.end())
3363 : {
3364 2 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid subdataset");
3365 2 : return false;
3366 : }
3367 :
3368 9 : m_aosSubdatasets.Clear();
3369 9 : const auto &pSuperGrid = oIter->second;
3370 9 : nRasterXSize = static_cast<int>(pSuperGrid.nWidth);
3371 9 : nRasterYSize = static_cast<int>(pSuperGrid.nHeight);
3372 :
3373 : // Convert from pixel-center convention to corner-pixel convention
3374 : const double dfMinX =
3375 9 : m_gt[0] + nX * m_gt[1] + pSuperGrid.fSWX - pSuperGrid.fResX / 2;
3376 9 : const double dfMinY = m_gt[3] + m_nLowResHeight * m_gt[5] +
3377 9 : nY * -m_gt[5] + pSuperGrid.fSWY -
3378 9 : pSuperGrid.fResY / 2;
3379 9 : const double dfMaxY =
3380 9 : dfMinY + pSuperGrid.nHeight * static_cast<double>(pSuperGrid.fResY);
3381 :
3382 9 : m_gt[0] = dfMinX;
3383 9 : m_gt[1] = pSuperGrid.fResX;
3384 9 : m_gt[3] = dfMaxY;
3385 9 : m_gt[5] = -pSuperGrid.fResY;
3386 9 : m_nSuperGridRefinementStartIndex = pSuperGrid.nIndex;
3387 :
3388 9 : if (!osGeorefMetadataLayer.empty())
3389 : {
3390 4 : auto poGeoref_metadataLayer = m_poRootGroup->OpenGroupFromFullname(
3391 4 : "/BAG_root/georef_metadata/" + osGeorefMetadataLayer, nullptr);
3392 4 : if (poGeoref_metadataLayer == nullptr)
3393 : {
3394 0 : CPLError(CE_Failure, CPLE_AppDefined,
3395 : "Cannot find georef_metadata layer %s",
3396 : osGeorefMetadataLayer.c_str());
3397 0 : return false;
3398 : }
3399 :
3400 8 : auto poKeys = poGeoref_metadataLayer->OpenMDArray("varres_keys");
3401 4 : if (poKeys != nullptr)
3402 : {
3403 2 : auto poDims = poKeys->GetDimensions();
3404 4 : if (poDims.size() != 2 || poDims[0]->GetSize() != 1 ||
3405 2 : poDims[1]->GetSize() != m_nRefinementsSize)
3406 : {
3407 0 : CPLError(CE_Failure, CPLE_AppDefined,
3408 : "Wrong dimensions for %s/varres_keys",
3409 : osGeorefMetadataLayer.c_str());
3410 0 : return false;
3411 : }
3412 4 : if (poKeys->GetDataType().GetClass() != GEDTC_NUMERIC ||
3413 2 : !GDALDataTypeIsInteger(
3414 2 : poKeys->GetDataType().GetNumericDataType()))
3415 : {
3416 0 : CPLError(
3417 : CE_Failure, CPLE_AppDefined,
3418 : "Only integer data type supported for %s/varres_keys",
3419 : osGeorefMetadataLayer.c_str());
3420 0 : return false;
3421 : }
3422 : }
3423 :
3424 8 : auto poValues = poGeoref_metadataLayer->OpenMDArray("values");
3425 4 : if (poValues == nullptr)
3426 : {
3427 0 : CPLError(CE_Failure, CPLE_AppDefined,
3428 : "Cannot find array values of georef_metadata layer %s",
3429 : osGeorefMetadataLayer.c_str());
3430 0 : return false;
3431 : }
3432 4 : if (poValues->GetDimensionCount() != 1)
3433 : {
3434 0 : CPLError(CE_Failure, CPLE_AppDefined,
3435 : "Wrong dimensions for %s/values",
3436 : osGeorefMetadataLayer.c_str());
3437 0 : return false;
3438 : }
3439 4 : if (poValues->GetDataType().GetClass() != GEDTC_COMPOUND)
3440 : {
3441 0 : CPLError(CE_Failure, CPLE_AppDefined,
3442 : "Only compound data type supported for %s/values",
3443 : osGeorefMetadataLayer.c_str());
3444 0 : return false;
3445 : }
3446 4 : SetBand(1, new BAGGeorefMDSuperGridBand(
3447 : poValues, poKeys,
3448 4 : new BAGSuperGridBand(this, 1, bHasNoData,
3449 4 : fNoDataValue)));
3450 : }
3451 : else
3452 : {
3453 15 : for (int i = 0; i < 2; i++)
3454 : {
3455 10 : SetBand(i + 1, new BAGSuperGridBand(this, i + 1, bHasNoData,
3456 10 : fNoDataValue));
3457 : }
3458 :
3459 5 : GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
3460 : "IMAGE_STRUCTURE");
3461 : }
3462 :
3463 9 : SetPhysicalFilename(osFilename);
3464 :
3465 9 : m_oMapRefinemendGrids.clear();
3466 : }
3467 :
3468 : // Setup/check for pam .aux.xml.
3469 74 : TryLoadXML();
3470 :
3471 : // Setup overviews.
3472 74 : oOvManager.Initialize(this, poOpenInfo->pszFilename);
3473 :
3474 74 : return true;
3475 : }
3476 :
3477 : /************************************************************************/
3478 : /* GetLayer() */
3479 : /************************************************************************/
3480 :
3481 47 : const OGRLayer *BAGDataset::GetLayer(int idx) const
3482 : {
3483 47 : if (idx != 0)
3484 1 : return nullptr;
3485 46 : return m_poTrackingListLayer.get();
3486 : }
3487 :
3488 : /************************************************************************/
3489 : /* BAGTrackingListLayer */
3490 : /************************************************************************/
3491 :
3492 : class BAGTrackingListLayer final
3493 : : public OGRLayer,
3494 : public OGRGetNextFeatureThroughRaw<BAGTrackingListLayer>
3495 : {
3496 : std::shared_ptr<GDALMDArray> m_poArray{};
3497 : OGRFeatureDefn *m_poFeatureDefn = nullptr;
3498 : int m_nIdx = 0;
3499 :
3500 : OGRFeature *GetNextRawFeature();
3501 :
3502 : CPL_DISALLOW_COPY_ASSIGN(BAGTrackingListLayer)
3503 :
3504 : public:
3505 : explicit BAGTrackingListLayer(const std::shared_ptr<GDALMDArray> &poArray);
3506 : ~BAGTrackingListLayer() override;
3507 :
3508 0 : const OGRFeatureDefn *GetLayerDefn() const override
3509 : {
3510 0 : return m_poFeatureDefn;
3511 : }
3512 :
3513 : void ResetReading() override;
3514 5 : DEFINE_GET_NEXT_FEATURE_THROUGH_RAW(BAGTrackingListLayer)
3515 :
3516 0 : int TestCapability(const char *) const override
3517 : {
3518 0 : return false;
3519 : }
3520 : };
3521 :
3522 : /************************************************************************/
3523 : /* BAGTrackingListLayer() */
3524 : /************************************************************************/
3525 :
3526 45 : BAGTrackingListLayer::BAGTrackingListLayer(
3527 45 : const std::shared_ptr<GDALMDArray> &poArray)
3528 45 : : m_poArray(poArray)
3529 : {
3530 45 : m_poFeatureDefn = new OGRFeatureDefn("tracking_list");
3531 45 : SetDescription(m_poFeatureDefn->GetName());
3532 45 : m_poFeatureDefn->Reference();
3533 45 : m_poFeatureDefn->SetGeomType(wkbNone);
3534 :
3535 45 : const auto &poComponents = poArray->GetDataType().GetComponents();
3536 315 : for (const auto &poComponent : poComponents)
3537 : {
3538 270 : if (poComponent->GetType().GetClass() == GEDTC_NUMERIC)
3539 : {
3540 : OGRFieldType eType;
3541 270 : if (GDALDataTypeIsInteger(
3542 540 : poComponent->GetType().GetNumericDataType()))
3543 180 : eType = OFTInteger;
3544 : else
3545 90 : eType = OFTReal;
3546 540 : OGRFieldDefn oFieldDefn(poComponent->GetName().c_str(), eType);
3547 270 : m_poFeatureDefn->AddFieldDefn(&oFieldDefn);
3548 : }
3549 : }
3550 45 : }
3551 :
3552 : /************************************************************************/
3553 : /* ~BAGTrackingListLayer() */
3554 : /************************************************************************/
3555 :
3556 90 : BAGTrackingListLayer::~BAGTrackingListLayer()
3557 : {
3558 45 : m_poFeatureDefn->Release();
3559 90 : }
3560 :
3561 : /************************************************************************/
3562 : /* ResetReading() */
3563 : /************************************************************************/
3564 :
3565 3 : void BAGTrackingListLayer::ResetReading()
3566 : {
3567 3 : m_nIdx = 0;
3568 3 : }
3569 :
3570 : /************************************************************************/
3571 : /* GetNextRawFeature() */
3572 : /************************************************************************/
3573 :
3574 5 : OGRFeature *BAGTrackingListLayer::GetNextRawFeature()
3575 : {
3576 5 : if (static_cast<GUInt64>(m_nIdx) >=
3577 5 : m_poArray->GetDimensions()[0]->GetSize())
3578 1 : return nullptr;
3579 :
3580 4 : const auto &oDataType = m_poArray->GetDataType();
3581 4 : std::vector<GByte> abyRow(oDataType.GetSize());
3582 :
3583 4 : const GUInt64 arrayStartIdx = static_cast<GUInt64>(m_nIdx);
3584 4 : const size_t count = 1;
3585 4 : const GInt64 arrayStep = 0;
3586 4 : const GPtrDiff_t bufferStride = 0;
3587 8 : m_poArray->Read(&arrayStartIdx, &count, &arrayStep, &bufferStride,
3588 4 : oDataType, &abyRow[0]);
3589 4 : int iCol = 0;
3590 4 : auto poFeature = new OGRFeature(m_poFeatureDefn);
3591 4 : poFeature->SetFID(m_nIdx);
3592 4 : m_nIdx++;
3593 :
3594 4 : const auto &poComponents = oDataType.GetComponents();
3595 28 : for (const auto &poComponent : poComponents)
3596 : {
3597 24 : if (poComponent->GetType().GetClass() == GEDTC_NUMERIC)
3598 : {
3599 24 : if (GDALDataTypeIsInteger(
3600 48 : poComponent->GetType().GetNumericDataType()))
3601 : {
3602 16 : int nValue = 0;
3603 16 : GDALCopyWords(&abyRow[poComponent->GetOffset()],
3604 16 : poComponent->GetType().GetNumericDataType(), 0,
3605 : &nValue, GDT_Int32, 0, 1);
3606 16 : poFeature->SetField(iCol, nValue);
3607 : }
3608 : else
3609 : {
3610 8 : double dfValue = 0;
3611 8 : GDALCopyWords(&abyRow[poComponent->GetOffset()],
3612 8 : poComponent->GetType().GetNumericDataType(), 0,
3613 : &dfValue, GDT_Float64, 0, 1);
3614 8 : poFeature->SetField(iCol, dfValue);
3615 : }
3616 24 : iCol++;
3617 : }
3618 : }
3619 :
3620 4 : return poFeature;
3621 : }
3622 :
3623 : /************************************************************************/
3624 : /* OpenVector() */
3625 : /************************************************************************/
3626 :
3627 46 : bool BAGDataset::OpenVector()
3628 : {
3629 : auto poTrackingList =
3630 138 : m_poRootGroup->OpenMDArrayFromFullname("/BAG_root/tracking_list");
3631 46 : if (!poTrackingList)
3632 1 : return false;
3633 45 : if (poTrackingList->GetDimensions().size() != 1)
3634 0 : return false;
3635 45 : if (poTrackingList->GetDataType().GetClass() != GEDTC_COMPOUND)
3636 0 : return false;
3637 :
3638 45 : m_poTrackingListLayer.reset(new BAGTrackingListLayer(poTrackingList));
3639 45 : return true;
3640 : }
3641 :
3642 : /************************************************************************/
3643 : /* OpenForCreate() */
3644 : /************************************************************************/
3645 :
3646 4 : GDALDataset *BAGDataset::OpenForCreate(GDALOpenInfo *poOpenInfo, int nXSizeIn,
3647 : int nYSizeIn, int nBandsIn,
3648 : CSLConstList papszCreationOptions)
3649 : {
3650 8 : CPLString osFilename(poOpenInfo->pszFilename);
3651 :
3652 : // Open the file as an HDF5 file.
3653 4 : hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
3654 4 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
3655 4 : hid_t hHDF5 = H5Fopen(osFilename, H5F_ACC_RDWR, fapl);
3656 4 : H5Pclose(fapl);
3657 4 : if (hHDF5 < 0)
3658 0 : return nullptr;
3659 :
3660 8 : auto poSharedResources = GDAL::HDF5SharedResources::Create(osFilename);
3661 4 : poSharedResources->m_hHDF5 = hHDF5;
3662 :
3663 8 : auto poRootGroup = HDF5Dataset::OpenGroup(poSharedResources);
3664 4 : if (poRootGroup == nullptr)
3665 0 : return nullptr;
3666 :
3667 : // Create a corresponding dataset.
3668 4 : BAGDataset *const poDS = new BAGDataset();
3669 :
3670 4 : poDS->eAccess = poOpenInfo->eAccess;
3671 4 : poDS->m_poRootGroup = std::move(poRootGroup);
3672 4 : poDS->m_poSharedResources = std::move(poSharedResources);
3673 4 : poDS->m_aosCreationOptions = papszCreationOptions;
3674 :
3675 4 : poDS->nRasterXSize = nXSizeIn;
3676 4 : poDS->nRasterYSize = nYSizeIn;
3677 :
3678 4 : const int nBlockSize = std::min(
3679 8 : 4096,
3680 4 : atoi(CSLFetchNameValueDef(papszCreationOptions, "BLOCK_SIZE", "100")));
3681 4 : const int nBlockXSize = std::min(poDS->nRasterXSize, nBlockSize);
3682 4 : const int nBlockYSize = std::min(poDS->nRasterYSize, nBlockSize);
3683 :
3684 10 : for (int i = 0; i < nBandsIn; i++)
3685 : {
3686 6 : auto poBand = new BAGRasterBand(poDS, i + 1);
3687 6 : poBand->nBlockXSize = nBlockXSize;
3688 6 : poBand->nBlockYSize = nBlockYSize;
3689 6 : poBand->eDataType = GDT_Float32;
3690 6 : poBand->m_bHasNoData = true;
3691 6 : poBand->m_fNoDataValue = fDEFAULT_NODATA;
3692 6 : poBand->GDALRasterBand::SetDescription(i == 0 ? "elevation"
3693 : : "uncertainty");
3694 6 : poDS->SetBand(i + 1, poBand);
3695 : }
3696 :
3697 4 : poDS->SetDescription(poOpenInfo->pszFilename);
3698 :
3699 4 : poDS->m_bReportVertCRS = CPLTestBool(CSLFetchNameValueDef(
3700 4 : poOpenInfo->papszOpenOptions, "REPORT_VERTCRS", "YES"));
3701 :
3702 4 : poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
3703 :
3704 : // Setup/check for pam .aux.xml.
3705 4 : poDS->TryLoadXML();
3706 :
3707 : // Setup overviews.
3708 4 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
3709 :
3710 4 : return poDS;
3711 : }
3712 :
3713 : /************************************************************************/
3714 : /* GetMeanSupergridsResolution() */
3715 : /************************************************************************/
3716 :
3717 5 : bool BAGDataset::GetMeanSupergridsResolution(double &dfResX, double &dfResY)
3718 : {
3719 5 : const int nChunkXSize = m_nChunkXSizeVarresMD;
3720 5 : const int nChunkYSize = m_nChunkYSizeVarresMD;
3721 :
3722 5 : dfResX = 0.0;
3723 5 : dfResY = 0.0;
3724 5 : int nValidSuperGrids = 0;
3725 5 : std::vector<BAGRefinementGrid> rgrids(static_cast<size_t>(nChunkXSize) *
3726 10 : nChunkYSize);
3727 5 : const int county = DIV_ROUND_UP(m_nLowResHeight, nChunkYSize);
3728 5 : const int countx = DIV_ROUND_UP(m_nLowResWidth, nChunkXSize);
3729 10 : for (int y = 0; y < county; y++)
3730 : {
3731 : const int nReqCountY =
3732 5 : std::min(nChunkYSize, m_nLowResHeight - y * nChunkYSize);
3733 10 : for (int x = 0; x < countx; x++)
3734 : {
3735 : const int nReqCountX =
3736 5 : std::min(nChunkXSize, m_nLowResWidth - x * nChunkXSize);
3737 :
3738 : // Create memory space to receive the data.
3739 5 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
3740 5 : static_cast<hsize_t>(nReqCountX)};
3741 5 : const hid_t memspace = H5Screate_simple(2, count, nullptr);
3742 5 : H5OFFSET_TYPE mem_offset[2] = {static_cast<H5OFFSET_TYPE>(0),
3743 : static_cast<H5OFFSET_TYPE>(0)};
3744 5 : if (H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
3745 5 : nullptr, count, nullptr) < 0)
3746 : {
3747 0 : H5Sclose(memspace);
3748 0 : return false;
3749 : }
3750 :
3751 5 : if (ReadVarresMetadataValue(y * nChunkYSize, x * nChunkXSize,
3752 : memspace, rgrids.data(), nReqCountY,
3753 : nReqCountX))
3754 : {
3755 125 : for (int i = 0; i < nReqCountX * nReqCountY; i++)
3756 : {
3757 120 : if (rgrids[i].nWidth > 0)
3758 : {
3759 120 : dfResX += rgrids[i].fResX;
3760 120 : dfResY += rgrids[i].fResY;
3761 120 : nValidSuperGrids++;
3762 : }
3763 : }
3764 : }
3765 5 : H5Sclose(memspace);
3766 : }
3767 : }
3768 :
3769 5 : if (nValidSuperGrids == 0)
3770 : {
3771 0 : CPLError(CE_Failure, CPLE_AppDefined, "No valid supergrids");
3772 0 : return false;
3773 : }
3774 :
3775 5 : dfResX /= nValidSuperGrids;
3776 5 : dfResY /= nValidSuperGrids;
3777 5 : return true;
3778 : }
3779 :
3780 : /************************************************************************/
3781 : /* GetVarresMetadataChunkSizes() */
3782 : /************************************************************************/
3783 :
3784 63 : void BAGDataset::GetVarresMetadataChunkSizes(int &nChunkXSize, int &nChunkYSize)
3785 : {
3786 63 : const hid_t listid = H5Dget_create_plist(m_hVarresMetadata);
3787 63 : nChunkXSize = m_nLowResWidth;
3788 63 : nChunkYSize = std::max(
3789 63 : 1, std::min(10 * 1024 * 1024 / m_nLowResWidth, m_nLowResHeight));
3790 63 : if (listid > 0)
3791 : {
3792 63 : if (H5Pget_layout(listid) == H5D_CHUNKED)
3793 : {
3794 0 : hsize_t panChunkDims[2] = {0, 0};
3795 0 : const int nDimSize = H5Pget_chunk(listid, 2, panChunkDims);
3796 0 : CPL_IGNORE_RET_VAL(nDimSize);
3797 0 : CPLAssert(nDimSize == 2);
3798 0 : nChunkXSize = static_cast<int>(panChunkDims[1]);
3799 0 : nChunkYSize = static_cast<int>(panChunkDims[0]);
3800 : }
3801 :
3802 63 : H5Pclose(listid);
3803 : }
3804 63 : }
3805 :
3806 : /************************************************************************/
3807 : /* GetVarresRefinementChunkSize() */
3808 : /************************************************************************/
3809 :
3810 63 : void BAGDataset::GetVarresRefinementChunkSize(unsigned &nChunkSize)
3811 : {
3812 63 : const hid_t listid = H5Dget_create_plist(m_hVarresRefinements);
3813 63 : nChunkSize = 1024;
3814 63 : if (listid > 0)
3815 : {
3816 63 : if (H5Pget_layout(listid) == H5D_CHUNKED)
3817 : {
3818 53 : hsize_t panChunkDims[2] = {0, 0};
3819 53 : const int nDimSize = H5Pget_chunk(listid, 2, panChunkDims);
3820 53 : CPL_IGNORE_RET_VAL(nDimSize);
3821 53 : CPLAssert(nDimSize == 2);
3822 53 : nChunkSize = static_cast<unsigned>(panChunkDims[1]);
3823 : }
3824 :
3825 63 : H5Pclose(listid);
3826 : }
3827 63 : }
3828 :
3829 : /************************************************************************/
3830 : /* GetRefinementValues() */
3831 : /************************************************************************/
3832 :
3833 107034 : const float *BAGDataset::GetRefinementValues(unsigned nRefinementIndex)
3834 : {
3835 107034 : unsigned nStartIndex = (nRefinementIndex / m_nChunkSizeVarresRefinement) *
3836 107034 : m_nChunkSizeVarresRefinement;
3837 107034 : const auto vPtr = m_oCacheRefinementValues.getPtr(nStartIndex);
3838 107034 : if (vPtr)
3839 107013 : return vPtr->data() + 2 * (nRefinementIndex - nStartIndex);
3840 :
3841 : const unsigned nCachedRefinementCount = std::min(
3842 21 : m_nChunkSizeVarresRefinement, m_nRefinementsSize - nStartIndex);
3843 42 : std::vector<float> values(2 * nCachedRefinementCount);
3844 :
3845 21 : hsize_t countVarresRefinements[2] = {
3846 21 : static_cast<hsize_t>(1), static_cast<hsize_t>(nCachedRefinementCount)};
3847 : const hid_t memspaceVarresRefinements =
3848 21 : H5Screate_simple(2, countVarresRefinements, nullptr);
3849 21 : H5OFFSET_TYPE mem_offset[2] = {static_cast<H5OFFSET_TYPE>(0),
3850 : static_cast<H5OFFSET_TYPE>(0)};
3851 21 : if (H5Sselect_hyperslab(memspaceVarresRefinements, H5S_SELECT_SET,
3852 : mem_offset, nullptr, countVarresRefinements,
3853 21 : nullptr) < 0)
3854 : {
3855 0 : H5Sclose(memspaceVarresRefinements);
3856 0 : return nullptr;
3857 : }
3858 :
3859 21 : H5OFFSET_TYPE offsetRefinement[2] = {
3860 21 : static_cast<H5OFFSET_TYPE>(0), static_cast<H5OFFSET_TYPE>(nStartIndex)};
3861 21 : if (H5Sselect_hyperslab(m_hVarresRefinementsDataspace, H5S_SELECT_SET,
3862 : offsetRefinement, nullptr, countVarresRefinements,
3863 21 : nullptr) < 0)
3864 : {
3865 0 : H5Sclose(memspaceVarresRefinements);
3866 0 : return nullptr;
3867 : }
3868 21 : if (H5Dread(m_hVarresRefinements, m_hVarresRefinementsNative,
3869 : memspaceVarresRefinements, m_hVarresRefinementsDataspace,
3870 42 : H5P_DEFAULT, values.data()) < 0)
3871 : {
3872 0 : H5Sclose(memspaceVarresRefinements);
3873 0 : return nullptr;
3874 : }
3875 21 : H5Sclose(memspaceVarresRefinements);
3876 : const auto &vRef =
3877 21 : m_oCacheRefinementValues.insert(nStartIndex, std::move(values));
3878 21 : return vRef.data() + 2 * (nRefinementIndex - nStartIndex);
3879 : }
3880 :
3881 : /************************************************************************/
3882 : /* ReadVarresMetadataValue() */
3883 : /************************************************************************/
3884 :
3885 336 : bool BAGDataset::ReadVarresMetadataValue(int y, int x, hid_t memspace,
3886 : BAGRefinementGrid *rgrid, int height,
3887 : int width)
3888 : {
3889 336 : constexpr int metadata_elt_size = 3 * 4 + 4 * 4; // 3 uint and 4 float
3890 336 : std::vector<char> buffer(static_cast<size_t>(metadata_elt_size) * height *
3891 672 : width);
3892 :
3893 336 : hsize_t count[2] = {static_cast<hsize_t>(height),
3894 336 : static_cast<hsize_t>(width)};
3895 336 : H5OFFSET_TYPE offset[2] = {static_cast<H5OFFSET_TYPE>(y),
3896 336 : static_cast<H5OFFSET_TYPE>(x)};
3897 336 : if (H5Sselect_hyperslab(m_hVarresMetadataDataspace, H5S_SELECT_SET, offset,
3898 336 : nullptr, count, nullptr) < 0)
3899 : {
3900 0 : CPLError(CE_Failure, CPLE_AppDefined,
3901 : "ReadVarresMetadataValue(): H5Sselect_hyperslab() failed");
3902 0 : return false;
3903 : }
3904 :
3905 336 : if (H5Dread(m_hVarresMetadata, m_hVarresMetadataNative, memspace,
3906 672 : m_hVarresMetadataDataspace, H5P_DEFAULT, buffer.data()) < 0)
3907 : {
3908 0 : CPLError(CE_Failure, CPLE_AppDefined,
3909 : "ReadVarresMetadataValue(): H5Dread() failed");
3910 0 : return false;
3911 : }
3912 :
3913 2319 : for (int i = 0; i < width * height; i++)
3914 : {
3915 1983 : const char *src_ptr = buffer.data() + metadata_elt_size * i;
3916 1983 : memcpy(&rgrid[i].nIndex, src_ptr, 4);
3917 1983 : memcpy(&rgrid[i].nWidth, src_ptr + 4, 4);
3918 1983 : memcpy(&rgrid[i].nHeight, src_ptr + 8, 4);
3919 1983 : memcpy(&rgrid[i].fResX, src_ptr + 12, 4);
3920 1983 : memcpy(&rgrid[i].fResY, src_ptr + 16, 4);
3921 1983 : memcpy(&rgrid[i].fSWX, src_ptr + 20, 4);
3922 1983 : memcpy(&rgrid[i].fSWY, src_ptr + 24, 4);
3923 : }
3924 336 : return true;
3925 : }
3926 :
3927 : /************************************************************************/
3928 : /* LookForRefinementGrids() */
3929 : /************************************************************************/
3930 :
3931 81 : bool BAGDataset::LookForRefinementGrids(CSLConstList l_papszOpenOptions,
3932 : int nYSubDS, int nXSubDS)
3933 :
3934 : {
3935 81 : m_hVarresMetadata = GH5DopenNoWarning(m_poSharedResources->m_hHDF5,
3936 : "/BAG_root/varres_metadata");
3937 81 : if (m_hVarresMetadata < 0)
3938 18 : return false;
3939 63 : m_hVarresRefinements =
3940 63 : H5Dopen(m_poSharedResources->m_hHDF5, "/BAG_root/varres_refinements");
3941 63 : if (m_hVarresRefinements < 0)
3942 0 : return false;
3943 :
3944 63 : m_hVarresMetadataDataType = H5Dget_type(m_hVarresMetadata);
3945 63 : if (H5Tget_class(m_hVarresMetadataDataType) != H5T_COMPOUND)
3946 : {
3947 0 : CPLError(CE_Failure, CPLE_NotSupported,
3948 : "m_hVarresMetadataDataType is not compound");
3949 0 : return false;
3950 : }
3951 :
3952 : const struct
3953 : {
3954 : const char *pszName;
3955 : hid_t eType;
3956 63 : } asMetadataFields[] = {
3957 63 : {"index", H5T_NATIVE_UINT}, {"dimensions_x", H5T_NATIVE_UINT},
3958 126 : {"dimensions_y", H5T_NATIVE_UINT}, {"resolution_x", H5T_NATIVE_FLOAT},
3959 126 : {"resolution_y", H5T_NATIVE_FLOAT}, {"sw_corner_x", H5T_NATIVE_FLOAT},
3960 63 : {"sw_corner_y", H5T_NATIVE_FLOAT},
3961 63 : };
3962 :
3963 63 : if (H5Tget_nmembers(m_hVarresMetadataDataType) !=
3964 : CPL_ARRAYSIZE(asMetadataFields))
3965 : {
3966 0 : CPLError(CE_Failure, CPLE_NotSupported,
3967 : "m_hVarresMetadataDataType has not %u members",
3968 : static_cast<unsigned>(CPL_ARRAYSIZE(asMetadataFields)));
3969 0 : return false;
3970 : }
3971 :
3972 504 : for (unsigned i = 0; i < CPL_ARRAYSIZE(asMetadataFields); i++)
3973 : {
3974 441 : char *pszName = H5Tget_member_name(m_hVarresMetadataDataType, i);
3975 441 : if (strcmp(pszName, asMetadataFields[i].pszName) != 0)
3976 : {
3977 0 : CPLError(CE_Failure, CPLE_NotSupported,
3978 : "asMetadataFields[%u].pszName = %s instead of %s",
3979 : static_cast<unsigned>(i), pszName,
3980 0 : asMetadataFields[i].pszName);
3981 0 : H5free_memory(pszName);
3982 0 : return false;
3983 : }
3984 441 : H5free_memory(pszName);
3985 441 : const hid_t type = H5Tget_member_type(m_hVarresMetadataDataType, i);
3986 441 : const hid_t hNativeType = H5Tget_native_type(type, H5T_DIR_DEFAULT);
3987 441 : bool bTypeOK = H5Tequal(asMetadataFields[i].eType, hNativeType) > 0;
3988 441 : H5Tclose(hNativeType);
3989 441 : H5Tclose(type);
3990 441 : if (!bTypeOK)
3991 : {
3992 0 : CPLError(CE_Failure, CPLE_NotSupported,
3993 : "asMetadataFields[%u].eType is not of expected type", i);
3994 0 : return false;
3995 : }
3996 : }
3997 :
3998 63 : m_hVarresMetadataDataspace = H5Dget_space(m_hVarresMetadata);
3999 63 : if (H5Sget_simple_extent_ndims(m_hVarresMetadataDataspace) != 2)
4000 : {
4001 0 : CPLDebug("BAG",
4002 : "H5Sget_simple_extent_ndims(m_hVarresMetadataDataspace) != 2");
4003 0 : return false;
4004 : }
4005 :
4006 : {
4007 63 : hsize_t dims[2] = {static_cast<hsize_t>(0), static_cast<hsize_t>(0)};
4008 63 : hsize_t maxdims[2] = {static_cast<hsize_t>(0), static_cast<hsize_t>(0)};
4009 :
4010 63 : H5Sget_simple_extent_dims(m_hVarresMetadataDataspace, dims, maxdims);
4011 63 : if (dims[0] != static_cast<hsize_t>(m_nLowResHeight) ||
4012 63 : dims[1] != static_cast<hsize_t>(m_nLowResWidth))
4013 : {
4014 0 : CPLDebug("BAG", "Unexpected dimension for m_hVarresMetadata");
4015 0 : return false;
4016 : }
4017 : }
4018 :
4019 : // We could potentially go beyond but we'd need to make sure that
4020 : // m_oMapRefinemendGrids is indexed by a int64_t
4021 63 : if (m_nLowResWidth > std::numeric_limits<int>::max() / m_nLowResHeight)
4022 : {
4023 0 : CPLError(CE_Failure, CPLE_NotSupported, "Too many refinement grids");
4024 0 : return false;
4025 : }
4026 :
4027 63 : m_hVarresMetadataNative =
4028 63 : H5Tget_native_type(m_hVarresMetadataDataType, H5T_DIR_ASCEND);
4029 :
4030 63 : m_hVarresRefinementsDataType = H5Dget_type(m_hVarresRefinements);
4031 63 : if (H5Tget_class(m_hVarresRefinementsDataType) != H5T_COMPOUND)
4032 : {
4033 0 : CPLError(CE_Failure, CPLE_NotSupported,
4034 : "m_hVarresRefinementsDataType is not compound");
4035 0 : return false;
4036 : }
4037 :
4038 : const struct
4039 : {
4040 : const char *pszName;
4041 : hid_t eType;
4042 63 : } asRefinementsFields[] = {
4043 63 : {"depth", H5T_NATIVE_FLOAT},
4044 63 : {"depth_uncrt", H5T_NATIVE_FLOAT},
4045 63 : };
4046 :
4047 63 : if (H5Tget_nmembers(m_hVarresRefinementsDataType) !=
4048 : CPL_ARRAYSIZE(asRefinementsFields))
4049 : {
4050 0 : CPLError(CE_Failure, CPLE_NotSupported,
4051 : "m_hVarresRefinementsDataType has not %u members",
4052 : static_cast<unsigned>(CPL_ARRAYSIZE(asRefinementsFields)));
4053 0 : return false;
4054 : }
4055 :
4056 189 : for (unsigned i = 0; i < CPL_ARRAYSIZE(asRefinementsFields); i++)
4057 : {
4058 126 : char *pszName = H5Tget_member_name(m_hVarresRefinementsDataType, i);
4059 126 : if (strcmp(pszName, asRefinementsFields[i].pszName) != 0)
4060 : {
4061 0 : CPLError(CE_Failure, CPLE_NotSupported,
4062 : "asRefinementsFields[%u].pszName = %s instead of %s",
4063 : static_cast<unsigned>(i), pszName,
4064 0 : asRefinementsFields[i].pszName);
4065 0 : H5free_memory(pszName);
4066 0 : return false;
4067 : }
4068 126 : H5free_memory(pszName);
4069 126 : const hid_t type = H5Tget_member_type(m_hVarresRefinementsDataType, i);
4070 126 : const hid_t hNativeType = H5Tget_native_type(type, H5T_DIR_DEFAULT);
4071 126 : bool bTypeOK = H5Tequal(asRefinementsFields[i].eType, hNativeType) > 0;
4072 126 : H5Tclose(hNativeType);
4073 126 : H5Tclose(type);
4074 126 : if (!bTypeOK)
4075 : {
4076 0 : CPLError(CE_Failure, CPLE_NotSupported,
4077 : "asRefinementsFields[%u].eType is not of expected type",
4078 : i);
4079 0 : return false;
4080 : }
4081 : }
4082 :
4083 63 : m_hVarresRefinementsDataspace = H5Dget_space(m_hVarresRefinements);
4084 63 : if (H5Sget_simple_extent_ndims(m_hVarresRefinementsDataspace) != 2)
4085 : {
4086 0 : CPLDebug(
4087 : "BAG",
4088 : "H5Sget_simple_extent_ndims(m_hVarresRefinementsDataspace) != 2");
4089 0 : return false;
4090 : }
4091 :
4092 63 : m_hVarresRefinementsNative =
4093 63 : H5Tget_native_type(m_hVarresRefinementsDataType, H5T_DIR_ASCEND);
4094 :
4095 : hsize_t nRefinementsSize;
4096 : {
4097 63 : hsize_t dims[2] = {static_cast<hsize_t>(0), static_cast<hsize_t>(0)};
4098 63 : hsize_t maxdims[2] = {static_cast<hsize_t>(0), static_cast<hsize_t>(0)};
4099 :
4100 63 : H5Sget_simple_extent_dims(m_hVarresRefinementsDataspace, dims, maxdims);
4101 63 : if (dims[0] != 1)
4102 : {
4103 0 : CPLDebug("BAG", "Unexpected dimension for m_hVarresRefinements");
4104 0 : return false;
4105 : }
4106 63 : nRefinementsSize = dims[1];
4107 63 : m_nRefinementsSize = static_cast<unsigned>(nRefinementsSize);
4108 63 : CPLDebug("BAG", "m_nRefinementsSize = %u", m_nRefinementsSize);
4109 : }
4110 :
4111 63 : GetVarresMetadataChunkSizes(m_nChunkXSizeVarresMD, m_nChunkYSizeVarresMD);
4112 63 : CPLDebug("BAG", "m_nChunkXSizeVarresMD = %d, m_nChunkYSizeVarresMD = %d",
4113 : m_nChunkXSizeVarresMD, m_nChunkYSizeVarresMD);
4114 63 : GetVarresRefinementChunkSize(m_nChunkSizeVarresRefinement);
4115 63 : CPLDebug("BAG", "m_nChunkSizeVarresRefinement = %u",
4116 : m_nChunkSizeVarresRefinement);
4117 :
4118 63 : const char *pszMode = CSLFetchNameValueDef(l_papszOpenOptions, "MODE", "");
4119 63 : if (EQUAL(pszMode, "RESAMPLED_GRID") || EQUAL(pszMode, "INTERPOLATED"))
4120 : {
4121 26 : return true;
4122 : }
4123 :
4124 : const char *pszSUPERGRIDS =
4125 37 : CSLFetchNameValue(l_papszOpenOptions, "SUPERGRIDS_INDICES");
4126 :
4127 : struct yx
4128 : {
4129 : int y;
4130 : int x;
4131 :
4132 14 : yx(int yin, int xin) : y(yin), x(xin)
4133 : {
4134 14 : }
4135 :
4136 30 : bool operator<(const yx &other) const
4137 : {
4138 30 : return y < other.y || (y == other.y && x < other.x);
4139 : }
4140 : };
4141 :
4142 74 : std::set<yx> oSupergrids;
4143 37 : int nMinX = 0;
4144 37 : int nMinY = 0;
4145 37 : int nMaxX = m_nLowResWidth - 1;
4146 37 : int nMaxY = m_nLowResHeight - 1;
4147 37 : if (nYSubDS >= 0 && nXSubDS >= 0)
4148 : {
4149 21 : nMinX = nXSubDS;
4150 21 : nMaxX = nXSubDS;
4151 21 : nMinY = nYSubDS;
4152 21 : nMaxY = nYSubDS;
4153 : }
4154 16 : else if (pszSUPERGRIDS)
4155 : {
4156 11 : char chExpectedChar = '(';
4157 11 : bool bNextIsY = false;
4158 11 : bool bNextIsX = false;
4159 11 : bool bHasY = false;
4160 11 : bool bHasX = false;
4161 11 : int nY = 0;
4162 11 : int i = 0;
4163 46 : for (; pszSUPERGRIDS[i]; i++)
4164 : {
4165 39 : if (chExpectedChar && pszSUPERGRIDS[i] != chExpectedChar)
4166 : {
4167 1 : CPLError(
4168 : CE_Warning, CPLE_AppDefined,
4169 : "Invalid formatting for SUPERGRIDS_INDICES at index %d. "
4170 : "Expecting %c, got %c",
4171 1 : i, chExpectedChar, pszSUPERGRIDS[i]);
4172 1 : break;
4173 : }
4174 38 : else if (chExpectedChar == '(')
4175 : {
4176 10 : chExpectedChar = 0;
4177 10 : bNextIsY = true;
4178 : }
4179 28 : else if (chExpectedChar == ',')
4180 : {
4181 2 : chExpectedChar = '(';
4182 : }
4183 : else
4184 : {
4185 26 : CPLAssert(chExpectedChar == 0);
4186 26 : if (bNextIsY && pszSUPERGRIDS[i] >= '0' &&
4187 9 : pszSUPERGRIDS[i] <= '9')
4188 : {
4189 8 : nY = atoi(pszSUPERGRIDS + i);
4190 8 : bNextIsY = false;
4191 8 : bHasY = true;
4192 : }
4193 18 : else if (bNextIsX && pszSUPERGRIDS[i] >= '0' &&
4194 5 : pszSUPERGRIDS[i] <= '9')
4195 : {
4196 4 : int nX = atoi(pszSUPERGRIDS + i);
4197 4 : bNextIsX = false;
4198 4 : oSupergrids.insert(yx(nY, nX));
4199 4 : bHasX = true;
4200 : }
4201 14 : else if ((bHasX || bHasY) && pszSUPERGRIDS[i] >= '0' &&
4202 1 : pszSUPERGRIDS[i] <= '9')
4203 : {
4204 : // ok
4205 : }
4206 14 : else if (bHasY && pszSUPERGRIDS[i] == ',')
4207 : {
4208 7 : bNextIsX = true;
4209 : }
4210 7 : else if (bHasX && bHasY && pszSUPERGRIDS[i] == ')')
4211 : {
4212 4 : chExpectedChar = ',';
4213 4 : bHasX = false;
4214 4 : bHasY = false;
4215 : }
4216 : else
4217 : {
4218 3 : CPLError(CE_Warning, CPLE_AppDefined,
4219 : "Invalid formatting for SUPERGRIDS_INDICES at "
4220 : "index %d. "
4221 : "Got %c",
4222 3 : i, pszSUPERGRIDS[i]);
4223 3 : break;
4224 : }
4225 : }
4226 : }
4227 11 : if (pszSUPERGRIDS[i] == 0 && chExpectedChar != ',')
4228 : {
4229 5 : CPLError(CE_Warning, CPLE_AppDefined,
4230 : "Invalid formatting for SUPERGRIDS_INDICES at index %d.",
4231 : i);
4232 : }
4233 :
4234 11 : bool bFirst = true;
4235 15 : for (const auto &yxPair : oSupergrids)
4236 : {
4237 4 : if (bFirst)
4238 : {
4239 3 : nMinX = yxPair.x;
4240 3 : nMaxX = yxPair.x;
4241 3 : nMinY = yxPair.y;
4242 3 : nMaxY = yxPair.y;
4243 3 : bFirst = false;
4244 : }
4245 : else
4246 : {
4247 1 : nMinX = std::min(nMinX, yxPair.x);
4248 1 : nMaxX = std::max(nMaxX, yxPair.x);
4249 1 : nMinY = std::min(nMinY, yxPair.y);
4250 1 : nMaxY = std::max(nMaxY, yxPair.y);
4251 : }
4252 : }
4253 : }
4254 37 : const char *pszMinX = CSLFetchNameValue(l_papszOpenOptions, "MINX");
4255 37 : const char *pszMinY = CSLFetchNameValue(l_papszOpenOptions, "MINY");
4256 37 : const char *pszMaxX = CSLFetchNameValue(l_papszOpenOptions, "MAXX");
4257 37 : const char *pszMaxY = CSLFetchNameValue(l_papszOpenOptions, "MAXY");
4258 37 : const int nCountBBoxElts = (pszMinX ? 1 : 0) + (pszMinY ? 1 : 0) +
4259 37 : (pszMaxX ? 1 : 0) + (pszMaxY ? 1 : 0);
4260 37 : const bool bHasBoundingBoxFilter =
4261 37 : !(nYSubDS >= 0 && nXSubDS >= 0) && (nCountBBoxElts == 4);
4262 37 : double dfFilterMinX = 0.0;
4263 37 : double dfFilterMinY = 0.0;
4264 37 : double dfFilterMaxX = 0.0;
4265 37 : double dfFilterMaxY = 0.0;
4266 37 : if (nYSubDS >= 0 && nXSubDS >= 0)
4267 : {
4268 : // do nothing
4269 : }
4270 16 : else if (bHasBoundingBoxFilter)
4271 : {
4272 1 : dfFilterMinX = CPLAtof(pszMinX);
4273 1 : dfFilterMinY = CPLAtof(pszMinY);
4274 1 : dfFilterMaxX = CPLAtof(pszMaxX);
4275 1 : dfFilterMaxY = CPLAtof(pszMaxY);
4276 :
4277 1 : nMinX = std::max(nMinX,
4278 1 : static_cast<int>((dfFilterMinX - m_gt[0]) / m_gt[1]));
4279 1 : nMaxX = std::min(nMaxX,
4280 1 : static_cast<int>((dfFilterMaxX - m_gt[0]) / m_gt[1]));
4281 :
4282 1 : nMinY = std::max(
4283 2 : nMinY, static_cast<int>(
4284 1 : (dfFilterMinY - (m_gt[3] + m_nLowResHeight * m_gt[5])) /
4285 1 : -m_gt[5]));
4286 1 : nMaxY = std::min(
4287 1 : nMaxY, static_cast<int>(
4288 1 : (dfFilterMaxY - (m_gt[3] + m_nLowResHeight * m_gt[5])) /
4289 1 : -m_gt[5]));
4290 : }
4291 15 : else if (nCountBBoxElts > 0)
4292 : {
4293 1 : CPLError(CE_Warning, CPLE_AppDefined,
4294 : "Bounding box filter ignored since only part of "
4295 : "MINX, MINY, MAXX and MAXY has been specified");
4296 : }
4297 :
4298 37 : const double dfResFilterMin = CPLAtof(
4299 : CSLFetchNameValueDef(l_papszOpenOptions, "RES_FILTER_MIN", "0"));
4300 37 : const double dfResFilterMax = CPLAtof(
4301 : CSLFetchNameValueDef(l_papszOpenOptions, "RES_FILTER_MAX", "inf"));
4302 :
4303 74 : std::vector<std::string> georefMDLayerNames;
4304 37 : auto poGeoref_metadata = m_poRootGroup->OpenGroupFromFullname(
4305 111 : "/BAG_root/georef_metadata", nullptr);
4306 37 : if (poGeoref_metadata)
4307 : {
4308 18 : const auto groupNames = poGeoref_metadata->GetGroupNames(nullptr);
4309 27 : for (const auto &groupName : groupNames)
4310 : {
4311 18 : georefMDLayerNames.push_back(groupName);
4312 : }
4313 : }
4314 :
4315 : const int nMaxSizeMap =
4316 37 : atoi(CPLGetConfigOption("GDAL_BAG_MAX_SIZE_VARRES_MAP", "50000000"));
4317 37 : const int nChunkXSize = m_nChunkXSizeVarresMD;
4318 37 : const int nChunkYSize = m_nChunkYSizeVarresMD;
4319 37 : std::vector<BAGRefinementGrid> rgrids(static_cast<size_t>(nChunkXSize) *
4320 74 : nChunkYSize);
4321 37 : bool bOK = true;
4322 74 : for (int blockY = nMinY / nChunkYSize; bOK && blockY <= nMaxY / nChunkYSize;
4323 : blockY++)
4324 : {
4325 : int nReqCountY =
4326 37 : std::min(nChunkYSize, m_nLowResHeight - blockY * nChunkYSize);
4327 74 : for (int blockX = nMinX / nChunkXSize;
4328 74 : bOK && blockX <= nMaxX / nChunkXSize; blockX++)
4329 : {
4330 : int nReqCountX =
4331 37 : std::min(nChunkXSize, m_nLowResWidth - blockX * nChunkXSize);
4332 :
4333 : // Create memory space to receive the data.
4334 37 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
4335 37 : static_cast<hsize_t>(nReqCountX)};
4336 37 : const hid_t memspace = H5Screate_simple(2, count, nullptr);
4337 37 : H5OFFSET_TYPE mem_offset[2] = {static_cast<H5OFFSET_TYPE>(0),
4338 : static_cast<H5OFFSET_TYPE>(0)};
4339 37 : if (H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
4340 37 : nullptr, count, nullptr) < 0)
4341 : {
4342 0 : H5Sclose(memspace);
4343 0 : bOK = false;
4344 0 : break;
4345 : }
4346 :
4347 74 : if (!ReadVarresMetadataValue(blockY * nChunkYSize,
4348 37 : blockX * nChunkXSize, memspace,
4349 : rgrids.data(), nReqCountY, nReqCountX))
4350 : {
4351 0 : bOK = false;
4352 0 : H5Sclose(memspace);
4353 0 : break;
4354 : }
4355 37 : H5Sclose(memspace);
4356 :
4357 113 : for (int yInBlock = std::max(0, nMinY - blockY * nChunkYSize);
4358 226 : bOK && yInBlock <= std::min(nReqCountY - 1,
4359 226 : nMaxY - blockY * nChunkYSize);
4360 : yInBlock++)
4361 : {
4362 413 : for (int xInBlock = std::max(0, nMinX - blockX * nChunkXSize);
4363 826 : bOK && xInBlock <= std::min(nReqCountX - 1,
4364 826 : nMaxX - blockX * nChunkXSize);
4365 : xInBlock++)
4366 : {
4367 337 : const int y = yInBlock + blockY * nChunkYSize;
4368 337 : const int x = xInBlock + blockX * nChunkXSize;
4369 : const auto &rgrid =
4370 337 : rgrids[yInBlock * nReqCountX + xInBlock];
4371 337 : if (rgrid.nWidth > 0)
4372 : {
4373 337 : if (rgrid.fResX <= 0 || rgrid.fResY <= 0)
4374 : {
4375 0 : CPLError(CE_Failure, CPLE_NotSupported,
4376 : "Incorrect resolution for supergrid "
4377 : "(%d, %d).",
4378 : static_cast<int>(y), static_cast<int>(x));
4379 0 : bOK = false;
4380 0 : break;
4381 : }
4382 337 : if (rgrid.nIndex + static_cast<GUInt64>(rgrid.nWidth) *
4383 337 : rgrid.nHeight >
4384 : nRefinementsSize)
4385 : {
4386 0 : CPLError(
4387 : CE_Failure, CPLE_NotSupported,
4388 : "Incorrect index / dimensions for supergrid "
4389 : "(%d, %d).",
4390 : static_cast<int>(y), static_cast<int>(x));
4391 0 : bOK = false;
4392 0 : break;
4393 : }
4394 :
4395 337 : if (rgrid.fSWX < 0.0f || rgrid.fSWY < 0.0f ||
4396 : // 0.1 is to deal with numeric imprecisions
4397 337 : rgrid.fSWX +
4398 337 : (rgrid.nWidth - 1 - 0.1) * rgrid.fResX >
4399 1011 : m_gt[1] ||
4400 337 : rgrid.fSWY +
4401 337 : (rgrid.nHeight - 1 - 0.1) * rgrid.fResY >
4402 337 : -m_gt[5])
4403 : {
4404 0 : CPLError(
4405 : CE_Failure, CPLE_NotSupported,
4406 : "Incorrect bounds for supergrid "
4407 : "(%d, %d): %f, %f, %f, %f.",
4408 : static_cast<int>(y), static_cast<int>(x),
4409 0 : rgrid.fSWX, rgrid.fSWY,
4410 0 : rgrid.fSWX + (rgrid.nWidth - 1) * rgrid.fResX,
4411 0 : rgrid.fSWY + (rgrid.nHeight - 1) * rgrid.fResY);
4412 0 : bOK = false;
4413 0 : break;
4414 : }
4415 :
4416 : const float gridRes =
4417 337 : std::max(rgrid.fResX, rgrid.fResY);
4418 337 : if (gridRes < dfResFilterMin ||
4419 337 : gridRes >= dfResFilterMax)
4420 : {
4421 12 : continue;
4422 : }
4423 :
4424 325 : if (static_cast<int>(m_oMapRefinemendGrids.size()) ==
4425 : nMaxSizeMap)
4426 : {
4427 0 : CPLError(
4428 : CE_Failure, CPLE_AppDefined,
4429 : "Size of map of refinement grids has reached "
4430 : "%d entries. "
4431 : "Set the GDAL_BAG_MAX_SIZE_VARRES_MAP "
4432 : "configuration option "
4433 : "to an higher value if you want to allow more",
4434 : nMaxSizeMap);
4435 0 : return false;
4436 : }
4437 :
4438 : try
4439 : {
4440 325 : m_oMapRefinemendGrids[y * m_nLowResWidth + x] =
4441 : rgrid;
4442 : }
4443 0 : catch (const std::exception &e)
4444 : {
4445 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
4446 : "Out of memory adding entries to map of "
4447 : "refinement "
4448 : "grids: %s",
4449 0 : e.what());
4450 0 : return false;
4451 : }
4452 :
4453 325 : const double dfMinX = m_gt[0] + x * m_gt[1] +
4454 325 : rgrid.fSWX - rgrid.fResX / 2;
4455 325 : const double dfMaxX =
4456 : dfMinX +
4457 325 : rgrid.nWidth * static_cast<double>(rgrid.fResX);
4458 : const double dfMinY =
4459 325 : m_gt[3] + m_nLowResHeight * m_gt[5] + y * -m_gt[5] +
4460 325 : rgrid.fSWY - rgrid.fResY / 2;
4461 325 : const double dfMaxY =
4462 : dfMinY +
4463 325 : static_cast<double>(rgrid.nHeight) * rgrid.fResY;
4464 :
4465 335 : if ((oSupergrids.empty() ||
4466 10 : oSupergrids.find(yx(static_cast<int>(y),
4467 10 : static_cast<int>(x))) !=
4468 650 : oSupergrids.end()) &&
4469 319 : (!bHasBoundingBoxFilter ||
4470 16 : (dfMinX >= dfFilterMinX &&
4471 12 : dfMinY >= dfFilterMinY &&
4472 9 : dfMaxX <= dfFilterMaxX &&
4473 : dfMaxY <= dfFilterMaxY)))
4474 : {
4475 : {
4476 : const int nIdx =
4477 305 : m_aosSubdatasets.size() / 2 + 1;
4478 : m_aosSubdatasets.AddNameValue(
4479 : CPLSPrintf("SUBDATASET_%d_NAME", nIdx),
4480 : CPLSPrintf("BAG:\"%s\":supergrid:%d:%d",
4481 305 : GetDescription(),
4482 : static_cast<int>(y),
4483 305 : static_cast<int>(x)));
4484 : m_aosSubdatasets.AddNameValue(
4485 : CPLSPrintf("SUBDATASET_%d_DESC", nIdx),
4486 : CPLSPrintf(
4487 : "Supergrid (y=%d, x=%d) from "
4488 : "(x=%f,y=%f) to "
4489 : "(x=%f,y=%f), resolution (x=%f,y=%f)",
4490 : static_cast<int>(y),
4491 : static_cast<int>(x), dfMinX, dfMinY,
4492 305 : dfMaxX, dfMaxY, rgrid.fResX,
4493 305 : rgrid.fResY));
4494 : }
4495 :
4496 369 : for (const auto &groupName : georefMDLayerNames)
4497 : {
4498 : const int nIdx =
4499 64 : m_aosSubdatasets.size() / 2 + 1;
4500 : m_aosSubdatasets.AddNameValue(
4501 : CPLSPrintf("SUBDATASET_%d_NAME", nIdx),
4502 : CPLSPrintf(
4503 : "BAG:\"%s\":georef_metadata:%s:%d:%d",
4504 64 : GetDescription(), groupName.c_str(),
4505 : static_cast<int>(y),
4506 64 : static_cast<int>(x)));
4507 : m_aosSubdatasets.AddNameValue(
4508 : CPLSPrintf("SUBDATASET_%d_DESC", nIdx),
4509 : CPLSPrintf("Georeferenced metadata %s of "
4510 : "supergrid (y=%d, x=%d)",
4511 : groupName.c_str(),
4512 : static_cast<int>(y),
4513 64 : static_cast<int>(x)));
4514 : }
4515 : }
4516 : }
4517 : }
4518 : }
4519 : }
4520 : }
4521 :
4522 37 : if (!bOK || m_oMapRefinemendGrids.empty())
4523 : {
4524 2 : m_aosSubdatasets.Clear();
4525 2 : m_oMapRefinemendGrids.clear();
4526 2 : return false;
4527 : }
4528 :
4529 35 : CPLDebug("BAG", "variable resolution extensions detected");
4530 35 : return true;
4531 : }
4532 :
4533 : /************************************************************************/
4534 : /* LoadMetadata() */
4535 : /************************************************************************/
4536 :
4537 81 : void BAGDataset::LoadMetadata()
4538 :
4539 : {
4540 : // Load the metadata from the file.
4541 : const hid_t hMDDS =
4542 81 : H5Dopen(m_poSharedResources->m_hHDF5, "/BAG_root/metadata");
4543 81 : const hid_t datatype = H5Dget_type(hMDDS);
4544 81 : const hid_t dataspace = H5Dget_space(hMDDS);
4545 81 : const hid_t native = H5Tget_native_type(datatype, H5T_DIR_ASCEND);
4546 :
4547 81 : const int n_dims = H5Sget_simple_extent_ndims(dataspace);
4548 81 : hsize_t dims[1] = {static_cast<hsize_t>(0)};
4549 81 : hsize_t maxdims[1] = {static_cast<hsize_t>(0)};
4550 :
4551 80 : if (n_dims == 1 && H5Tget_class(native) == H5T_STRING &&
4552 161 : !H5Tis_variable_str(native) && H5Tget_size(native) == 1)
4553 : {
4554 80 : H5Sget_simple_extent_dims(dataspace, dims, maxdims);
4555 :
4556 80 : pszXMLMetadata =
4557 80 : static_cast<char *>(CPLCalloc(static_cast<int>(dims[0] + 1), 1));
4558 :
4559 80 : H5Dread(hMDDS, native, H5S_ALL, dataspace, H5P_DEFAULT, pszXMLMetadata);
4560 : }
4561 :
4562 81 : H5Tclose(native);
4563 81 : H5Sclose(dataspace);
4564 81 : H5Tclose(datatype);
4565 81 : H5Dclose(hMDDS);
4566 :
4567 81 : if (pszXMLMetadata == nullptr || pszXMLMetadata[0] == 0)
4568 1 : return;
4569 :
4570 : // Try to get the geotransform.
4571 80 : CPLXMLNode *psRoot = CPLParseXMLString(pszXMLMetadata);
4572 :
4573 80 : if (psRoot == nullptr)
4574 0 : return;
4575 :
4576 80 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
4577 :
4578 80 : CPLXMLNode *const psGeo = CPLSearchXMLNode(psRoot, "=MD_Georectified");
4579 :
4580 80 : if (psGeo != nullptr)
4581 : {
4582 160 : CPLString osResHeight, osResWidth;
4583 722 : for (const auto *psIter = psGeo->psChild; psIter;
4584 642 : psIter = psIter->psNext)
4585 : {
4586 642 : if (strcmp(psIter->pszValue, "axisDimensionProperties") == 0)
4587 : {
4588 : // since BAG format 1.5 version
4589 160 : const char *pszDim = CPLGetXMLValue(
4590 : psIter,
4591 : "MD_Dimension.dimensionName.MD_DimensionNameTypeCode",
4592 : nullptr);
4593 160 : const char *pszRes = nullptr;
4594 160 : if (pszDim)
4595 : {
4596 152 : pszRes = CPLGetXMLValue(
4597 : psIter, "MD_Dimension.resolution.Measure", nullptr);
4598 : }
4599 : else
4600 : {
4601 : // prior to BAG format 1.5 version
4602 8 : pszDim = CPLGetXMLValue(
4603 : psIter, "MD_Dimension.dimensionName", nullptr);
4604 8 : pszRes = CPLGetXMLValue(
4605 : psIter, "MD_Dimension.resolution.Measure.value",
4606 : nullptr);
4607 : }
4608 :
4609 160 : if (pszDim && EQUAL(pszDim, "row") && pszRes)
4610 : {
4611 80 : osResHeight = pszRes;
4612 : }
4613 80 : else if (pszDim && EQUAL(pszDim, "column") && pszRes)
4614 : {
4615 80 : osResWidth = pszRes;
4616 : }
4617 : }
4618 : }
4619 :
4620 80 : char **papszCornerTokens = CSLTokenizeStringComplex(
4621 : CPLGetXMLValue(psGeo, "cornerPoints.Point.coordinates", ""), " ,",
4622 : FALSE, FALSE);
4623 :
4624 80 : if (CSLCount(papszCornerTokens) == 4)
4625 : {
4626 80 : const double dfLLX = CPLAtof(papszCornerTokens[0]);
4627 80 : const double dfLLY = CPLAtof(papszCornerTokens[1]);
4628 80 : const double dfURX = CPLAtof(papszCornerTokens[2]);
4629 80 : const double dfURY = CPLAtof(papszCornerTokens[3]);
4630 :
4631 80 : double dfResWidth = CPLAtof(osResWidth);
4632 80 : double dfResHeight = CPLAtof(osResHeight);
4633 80 : if (dfResWidth > 0 && dfResHeight > 0)
4634 : {
4635 80 : if (fabs((dfURX - dfLLX) / dfResWidth - m_nLowResWidth) <
4636 0 : 1e-2 &&
4637 0 : fabs((dfURY - dfLLY) / dfResHeight - m_nLowResHeight) <
4638 : 1e-2)
4639 : {
4640 : // Found with
4641 : // https://data.ngdc.noaa.gov/platforms/ocean/nos/coast/H12001-H14000/H12525/BAG/H12525_MB_4m_MLLW_1of2.bag
4642 : // to address issue
4643 : // https://github.com/OSGeo/gdal/issues/1643
4644 0 : CPLError(CE_Warning, CPLE_AppDefined,
4645 : "cornerPoints not consistent with resolution "
4646 : "given in metadata");
4647 : }
4648 80 : else if (fabs((dfURX - dfLLX) / dfResWidth -
4649 80 : (m_nLowResWidth - 1)) < 1e-2 &&
4650 78 : fabs((dfURY - dfLLY) / dfResHeight -
4651 78 : (m_nLowResHeight - 1)) < 1e-2)
4652 : {
4653 : // pixel center convention. OK
4654 : }
4655 : else
4656 : {
4657 2 : CPLDebug("BAG", "cornerPoints not consistent with "
4658 : "resolution given in metadata");
4659 2 : CPLDebug("BAG",
4660 : "Metadata horizontal resolution: %f. "
4661 : "Computed resolution: %f. "
4662 : "Computed width: %f vs %d",
4663 2 : dfResWidth, (dfURX - dfLLX) / (m_nLowResWidth - 1),
4664 2 : (dfURX - dfLLX) / dfResWidth, m_nLowResWidth);
4665 2 : CPLDebug("BAG",
4666 : "Metadata vertical resolution: %f. "
4667 : "Computed resolution: %f. "
4668 : "Computed height: %f vs %d",
4669 : dfResHeight,
4670 2 : (dfURY - dfLLY) / (m_nLowResHeight - 1),
4671 2 : (dfURY - dfLLY) / dfResHeight, m_nLowResHeight);
4672 2 : CPLError(CE_Warning, CPLE_AppDefined,
4673 : "cornerPoints not consistent with resolution "
4674 : "given in metadata");
4675 : }
4676 : }
4677 :
4678 80 : m_gt[0] = dfLLX;
4679 80 : m_gt[1] = dfResWidth;
4680 80 : m_gt[3] = dfLLY + dfResHeight * (m_nLowResHeight - 1);
4681 80 : m_gt[5] = dfResHeight * (-1);
4682 :
4683 : // shift to pixel corner convention
4684 80 : m_gt[0] -= m_gt[1] * 0.5;
4685 80 : m_gt[3] -= m_gt[5] * 0.5;
4686 :
4687 80 : m_dfLowResMinX = m_gt[0];
4688 80 : m_dfLowResMaxX = m_dfLowResMinX + m_nLowResWidth * m_gt[1];
4689 80 : m_dfLowResMaxY = m_gt[3];
4690 80 : m_dfLowResMinY = m_dfLowResMaxY + m_nLowResHeight * m_gt[5];
4691 : }
4692 80 : CSLDestroy(papszCornerTokens);
4693 : }
4694 :
4695 : // Try to get the coordinate system.
4696 80 : if (OGR_SRS_ImportFromISO19115(&m_oSRS, pszXMLMetadata) != OGRERR_NONE)
4697 : {
4698 78 : ParseWKTFromXML(pszXMLMetadata);
4699 : }
4700 :
4701 : // Fetch acquisition date.
4702 80 : CPLXMLNode *const psDateTime = CPLSearchXMLNode(psRoot, "=dateTime");
4703 80 : if (psDateTime != nullptr)
4704 : {
4705 : const char *pszDateTimeValue =
4706 80 : psDateTime->psChild && psDateTime->psChild->eType == CXT_Element
4707 156 : ? CPLGetXMLValue(psDateTime->psChild, nullptr, nullptr)
4708 4 : : CPLGetXMLValue(psDateTime, nullptr, nullptr);
4709 80 : if (pszDateTimeValue)
4710 80 : GDALDataset::SetMetadataItem("BAG_DATETIME", pszDateTimeValue);
4711 : }
4712 :
4713 80 : CPLDestroyXMLNode(psRoot);
4714 : }
4715 :
4716 : /************************************************************************/
4717 : /* ParseWKTFromXML() */
4718 : /************************************************************************/
4719 78 : OGRErr BAGDataset::ParseWKTFromXML(const char *pszISOXML)
4720 : {
4721 78 : CPLXMLNode *const psRoot = CPLParseXMLString(pszISOXML);
4722 :
4723 78 : if (psRoot == nullptr)
4724 0 : return OGRERR_FAILURE;
4725 :
4726 78 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
4727 :
4728 78 : CPLXMLNode *psRSI = CPLSearchXMLNode(psRoot, "=referenceSystemInfo");
4729 78 : if (psRSI == nullptr)
4730 : {
4731 0 : CPLError(CE_Failure, CPLE_AppDefined,
4732 : "Unable to find <referenceSystemInfo> in metadata.");
4733 0 : CPLDestroyXMLNode(psRoot);
4734 0 : return OGRERR_FAILURE;
4735 : }
4736 :
4737 : const char *pszSRCodeString =
4738 78 : CPLGetXMLValue(psRSI,
4739 : "MD_ReferenceSystem.referenceSystemIdentifier."
4740 : "RS_Identifier.code.CharacterString",
4741 : nullptr);
4742 78 : if (pszSRCodeString == nullptr)
4743 : {
4744 2 : CPLDebug("BAG",
4745 : "Unable to find /MI_Metadata/referenceSystemInfo[1]/"
4746 : "MD_ReferenceSystem[1]/referenceSystemIdentifier[1]/"
4747 : "RS_Identifier[1]/code[1]/CharacterString[1] in metadata.");
4748 2 : CPLDestroyXMLNode(psRoot);
4749 2 : return OGRERR_FAILURE;
4750 : }
4751 :
4752 : const char *pszSRCodeSpace =
4753 76 : CPLGetXMLValue(psRSI,
4754 : "MD_ReferenceSystem.referenceSystemIdentifier."
4755 : "RS_Identifier.codeSpace.CharacterString",
4756 : "");
4757 76 : if (!EQUAL(pszSRCodeSpace, "WKT"))
4758 : {
4759 0 : CPLError(CE_Failure, CPLE_AppDefined,
4760 : "Spatial reference string is not in WKT.");
4761 0 : CPLDestroyXMLNode(psRoot);
4762 0 : return OGRERR_FAILURE;
4763 : }
4764 :
4765 76 : if (m_oSRS.importFromWkt(pszSRCodeString) != OGRERR_NONE)
4766 : {
4767 0 : CPLError(CE_Failure, CPLE_AppDefined,
4768 : "Failed parsing WKT string \"%s\".", pszSRCodeString);
4769 0 : CPLDestroyXMLNode(psRoot);
4770 0 : return OGRERR_FAILURE;
4771 : }
4772 :
4773 76 : psRSI = CPLSearchXMLNode(psRSI->psNext, "=referenceSystemInfo");
4774 76 : if (psRSI == nullptr)
4775 : {
4776 0 : CPLError(CE_Failure, CPLE_AppDefined,
4777 : "Unable to find second instance of <referenceSystemInfo> "
4778 : "in metadata.");
4779 0 : CPLDestroyXMLNode(psRoot);
4780 0 : return OGRERR_NONE;
4781 : }
4782 :
4783 : pszSRCodeString =
4784 76 : CPLGetXMLValue(psRSI,
4785 : "MD_ReferenceSystem.referenceSystemIdentifier."
4786 : "RS_Identifier.code.CharacterString",
4787 : nullptr);
4788 76 : if (pszSRCodeString == nullptr)
4789 : {
4790 0 : CPLError(CE_Failure, CPLE_AppDefined,
4791 : "Unable to find /MI_Metadata/referenceSystemInfo[2]/"
4792 : "MD_ReferenceSystem[1]/referenceSystemIdentifier[1]/"
4793 : "RS_Identifier[1]/code[1]/CharacterString[1] in metadata.");
4794 0 : CPLDestroyXMLNode(psRoot);
4795 0 : return OGRERR_NONE;
4796 : }
4797 :
4798 : pszSRCodeSpace =
4799 76 : CPLGetXMLValue(psRSI,
4800 : "MD_ReferenceSystem.referenceSystemIdentifier."
4801 : "RS_Identifier.codeSpace.CharacterString",
4802 : "");
4803 76 : if (!EQUAL(pszSRCodeSpace, "WKT"))
4804 : {
4805 0 : CPLError(CE_Failure, CPLE_AppDefined,
4806 : "Spatial reference string is not in WKT.");
4807 0 : CPLDestroyXMLNode(psRoot);
4808 0 : return OGRERR_NONE;
4809 : }
4810 :
4811 76 : if (m_bReportVertCRS && (STARTS_WITH_CI(pszSRCodeString, "VERTCS") ||
4812 73 : STARTS_WITH_CI(pszSRCodeString, "VERT_CS")))
4813 : {
4814 146 : OGR_SRSNode oVertCRSRootNode;
4815 73 : const char *pszInput = pszSRCodeString;
4816 73 : if (oVertCRSRootNode.importFromWkt(&pszInput) == OGRERR_NONE)
4817 : {
4818 73 : if (oVertCRSRootNode.GetNode("UNIT") == nullptr)
4819 : {
4820 : // UNIT is required
4821 68 : auto poUnits = new OGR_SRSNode("UNIT");
4822 68 : poUnits->AddChild(new OGR_SRSNode("metre"));
4823 68 : poUnits->AddChild(new OGR_SRSNode("1.0"));
4824 68 : oVertCRSRootNode.AddChild(poUnits);
4825 : }
4826 73 : if (oVertCRSRootNode.GetNode("AXIS") == nullptr)
4827 : {
4828 : // If AXIS is missing, add an explicit Depth AXIS
4829 68 : auto poAxis = new OGR_SRSNode("AXIS");
4830 68 : poAxis->AddChild(new OGR_SRSNode("Depth"));
4831 68 : poAxis->AddChild(new OGR_SRSNode("DOWN"));
4832 68 : oVertCRSRootNode.AddChild(poAxis);
4833 : }
4834 :
4835 73 : char *pszVertCRSWKT = nullptr;
4836 73 : oVertCRSRootNode.exportToWkt(&pszVertCRSWKT);
4837 :
4838 146 : OGRSpatialReference oVertCRS;
4839 73 : if (oVertCRS.importFromWkt(pszVertCRSWKT) == OGRERR_NONE)
4840 : {
4841 73 : if (EQUAL(oVertCRS.GetName(), "MLLW"))
4842 : {
4843 60 : oVertCRS.importFromEPSG(5866);
4844 : }
4845 :
4846 146 : OGRSpatialReference oCompoundCRS;
4847 73 : oCompoundCRS.SetCompoundCS(
4848 146 : (CPLString(m_oSRS.GetName()) + " + " + oVertCRS.GetName())
4849 : .c_str(),
4850 73 : &m_oSRS, &oVertCRS);
4851 73 : oCompoundCRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4852 :
4853 73 : m_oSRS = std::move(oCompoundCRS);
4854 : }
4855 :
4856 73 : CPLFree(pszVertCRSWKT);
4857 : }
4858 : }
4859 :
4860 76 : CPLDestroyXMLNode(psRoot);
4861 :
4862 76 : return OGRERR_NONE;
4863 : }
4864 :
4865 : /************************************************************************/
4866 : /* GetGeoTransform() */
4867 : /************************************************************************/
4868 :
4869 78 : CPLErr BAGDataset::GetGeoTransform(GDALGeoTransform >) const
4870 :
4871 : {
4872 78 : if (m_gt[0] != 0.0 || m_gt[3] != 0.0)
4873 : {
4874 78 : gt = m_gt;
4875 78 : return CE_None;
4876 : }
4877 :
4878 0 : return GDALPamDataset::GetGeoTransform(gt);
4879 : }
4880 :
4881 : /************************************************************************/
4882 : /* GetSpatialRef() */
4883 : /************************************************************************/
4884 :
4885 21 : const OGRSpatialReference *BAGDataset::GetSpatialRef() const
4886 : {
4887 21 : if (!m_oSRS.IsEmpty())
4888 21 : return &m_oSRS;
4889 0 : return GDALPamDataset::GetSpatialRef();
4890 : }
4891 :
4892 : /************************************************************************/
4893 : /* SetGeoTransform() */
4894 : /************************************************************************/
4895 :
4896 3 : CPLErr BAGDataset::SetGeoTransform(const GDALGeoTransform >)
4897 : {
4898 3 : if (eAccess == GA_ReadOnly)
4899 0 : return GDALPamDataset::SetGeoTransform(gt);
4900 :
4901 3 : if (gt[2] != 0 || gt[4] != 0)
4902 : {
4903 0 : CPLError(CE_Failure, CPLE_NotSupported,
4904 : "BAG driver requires a non-rotated geotransform");
4905 0 : return CE_Failure;
4906 : }
4907 3 : m_gt = gt;
4908 3 : return WriteMetadataIfNeeded() ? CE_None : CE_Failure;
4909 : }
4910 :
4911 : /************************************************************************/
4912 : /* SetSpatialRef() */
4913 : /************************************************************************/
4914 :
4915 3 : CPLErr BAGDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
4916 : {
4917 3 : if (eAccess == GA_ReadOnly)
4918 0 : return GDALPamDataset::SetSpatialRef(poSRS);
4919 :
4920 3 : if (poSRS == nullptr || poSRS->IsEmpty())
4921 : {
4922 0 : CPLError(CE_Failure, CPLE_NotSupported,
4923 : "BAG driver requires a valid SRS");
4924 0 : return CE_Failure;
4925 : }
4926 :
4927 3 : m_oSRS = *poSRS;
4928 3 : return WriteMetadataIfNeeded() ? CE_None : CE_Failure;
4929 : }
4930 :
4931 : /************************************************************************/
4932 : /* WriteMetadataIfNeeded() */
4933 : /************************************************************************/
4934 :
4935 6 : bool BAGDataset::WriteMetadataIfNeeded()
4936 : {
4937 6 : if (m_bMetadataWritten)
4938 : {
4939 0 : return true;
4940 : }
4941 6 : if (m_gt == GDALGeoTransform() || m_oSRS.IsEmpty())
4942 : {
4943 3 : return true;
4944 : }
4945 3 : m_bMetadataWritten = true;
4946 :
4947 : CPLString osXMLMetadata = BAGCreator::GenerateMetadata(
4948 3 : nRasterXSize, nRasterYSize, m_gt, m_oSRS.IsEmpty() ? nullptr : &m_oSRS,
4949 9 : m_aosCreationOptions.List());
4950 3 : if (osXMLMetadata.empty())
4951 : {
4952 0 : return false;
4953 : }
4954 :
4955 3 : if (!BAGCreator::CreateAndWriteMetadata(m_poSharedResources->m_hHDF5,
4956 : osXMLMetadata))
4957 : {
4958 0 : return false;
4959 : }
4960 :
4961 3 : return true;
4962 : }
4963 :
4964 : /************************************************************************/
4965 : /* GetMetadataDomainList() */
4966 : /************************************************************************/
4967 :
4968 0 : char **BAGDataset::GetMetadataDomainList()
4969 : {
4970 0 : return BuildMetadataDomainList(GDALPamDataset::GetMetadataDomainList(),
4971 0 : TRUE, "xml:BAG", nullptr);
4972 : }
4973 :
4974 : /************************************************************************/
4975 : /* GetMetadata() */
4976 : /************************************************************************/
4977 :
4978 67 : char **BAGDataset::GetMetadata(const char *pszDomain)
4979 :
4980 : {
4981 67 : if (pszDomain != nullptr && EQUAL(pszDomain, "xml:BAG"))
4982 : {
4983 5 : apszMDList[0] = pszXMLMetadata;
4984 5 : apszMDList[1] = nullptr;
4985 :
4986 5 : return apszMDList;
4987 : }
4988 :
4989 62 : if (pszDomain != nullptr && EQUAL(pszDomain, "SUBDATASETS"))
4990 : {
4991 16 : return m_aosSubdatasets.List();
4992 : }
4993 :
4994 46 : return GDALPamDataset::GetMetadata(pszDomain);
4995 : }
4996 :
4997 : /************************************************************************/
4998 : /* BAGDatasetDriverUnload() */
4999 : /************************************************************************/
5000 :
5001 9 : static void BAGDatasetDriverUnload(GDALDriver *)
5002 : {
5003 9 : HDF5UnloadFileDriver();
5004 9 : }
5005 :
5006 : /************************************************************************/
5007 : /* ~BAGCreator()() */
5008 : /************************************************************************/
5009 :
5010 148 : BAGCreator::~BAGCreator()
5011 : {
5012 74 : Close();
5013 74 : }
5014 :
5015 : /************************************************************************/
5016 : /* Close() */
5017 : /************************************************************************/
5018 :
5019 94 : bool BAGCreator::Close()
5020 : {
5021 94 : bool ret = true;
5022 94 : if (m_bagRoot >= 0)
5023 : {
5024 20 : ret = (H5_CHECK(H5Gclose(m_bagRoot)) >= 0) && ret;
5025 20 : m_bagRoot = -1;
5026 : }
5027 94 : if (m_hdf5 >= 0)
5028 : {
5029 20 : ret = (H5_CHECK(H5Fclose(m_hdf5)) >= 0) && ret;
5030 20 : m_hdf5 = -1;
5031 : }
5032 94 : return ret;
5033 : }
5034 :
5035 : /************************************************************************/
5036 : /* SubstituteVariables() */
5037 : /************************************************************************/
5038 :
5039 4893 : bool BAGCreator::SubstituteVariables(CPLXMLNode *psNode, char **papszDict)
5040 : {
5041 4893 : if (psNode->eType == CXT_Text && psNode->pszValue &&
5042 1617 : strstr(psNode->pszValue, "${"))
5043 : {
5044 735 : CPLString osVal(psNode->pszValue);
5045 735 : size_t nPos = 0;
5046 : while (true)
5047 : {
5048 1470 : nPos = osVal.find("${", nPos);
5049 1470 : if (nPos == std::string::npos)
5050 : {
5051 735 : break;
5052 : }
5053 735 : CPLString osKeyName;
5054 735 : bool bHasDefaultValue = false;
5055 735 : CPLString osDefaultValue;
5056 735 : size_t nAfterKeyName = 0;
5057 14616 : for (size_t i = nPos + 2; i < osVal.size(); i++)
5058 : {
5059 14616 : if (osVal[i] == ':')
5060 : {
5061 378 : osKeyName = osVal.substr(nPos + 2, i - (nPos + 2));
5062 : }
5063 14238 : else if (osVal[i] == '}')
5064 : {
5065 735 : if (osKeyName.empty())
5066 : {
5067 357 : osKeyName = osVal.substr(nPos + 2, i - (nPos + 2));
5068 : }
5069 : else
5070 : {
5071 378 : bHasDefaultValue = true;
5072 : size_t nStartDefaultVal =
5073 378 : nPos + 2 + osKeyName.size() + 1;
5074 756 : osDefaultValue = osVal.substr(nStartDefaultVal,
5075 378 : i - nStartDefaultVal);
5076 : }
5077 735 : nAfterKeyName = i + 1;
5078 735 : break;
5079 : }
5080 : }
5081 735 : if (nAfterKeyName == 0)
5082 : {
5083 0 : CPLError(CE_Failure, CPLE_AppDefined,
5084 : "Invalid variable name in template");
5085 0 : return false;
5086 : }
5087 :
5088 735 : bool bSubstFound = false;
5089 9511 : for (char **papszIter = papszDict;
5090 9511 : !bSubstFound && papszIter && *papszIter; papszIter++)
5091 : {
5092 8776 : if (STARTS_WITH_CI(*papszIter, "VAR_"))
5093 : {
5094 8706 : char *pszKey = nullptr;
5095 : const char *pszValue =
5096 8706 : CPLParseNameValue(*papszIter, &pszKey);
5097 8706 : if (pszKey && pszValue)
5098 : {
5099 8706 : const char *pszVarName = pszKey + strlen("VAR_");
5100 8706 : if (EQUAL(pszVarName, osKeyName))
5101 : {
5102 362 : bSubstFound = true;
5103 724 : osVal = osVal.substr(0, nPos) + pszValue +
5104 1086 : osVal.substr(nAfterKeyName);
5105 : }
5106 8706 : CPLFree(pszKey);
5107 : }
5108 : }
5109 : }
5110 735 : if (!bSubstFound)
5111 : {
5112 373 : if (bHasDefaultValue)
5113 : {
5114 746 : osVal = osVal.substr(0, nPos) + osDefaultValue +
5115 1119 : osVal.substr(nAfterKeyName);
5116 : }
5117 : else
5118 : {
5119 0 : CPLError(CE_Warning, CPLE_AppDefined,
5120 : "%s could not be substituted", osKeyName.c_str());
5121 0 : return false;
5122 : }
5123 : }
5124 735 : }
5125 :
5126 735 : if (!osVal.empty() && osVal[0] == '<' && osVal.back() == '>')
5127 : {
5128 2 : CPLXMLNode *psSubNode = CPLParseXMLString(osVal);
5129 2 : if (psSubNode)
5130 : {
5131 2 : CPLFree(psNode->pszValue);
5132 2 : psNode->eType = psSubNode->eType;
5133 2 : psNode->pszValue = psSubNode->pszValue;
5134 2 : psNode->psChild = psSubNode->psChild;
5135 2 : psSubNode->pszValue = nullptr;
5136 2 : psSubNode->psChild = nullptr;
5137 2 : CPLDestroyXMLNode(psSubNode);
5138 : }
5139 : else
5140 : {
5141 0 : CPLFree(psNode->pszValue);
5142 0 : psNode->pszValue = CPLStrdup(osVal);
5143 : }
5144 : }
5145 : else
5146 : {
5147 733 : CPLFree(psNode->pszValue);
5148 733 : psNode->pszValue = CPLStrdup(osVal);
5149 : }
5150 : }
5151 :
5152 9765 : for (CPLXMLNode *psIter = psNode->psChild; psIter; psIter = psIter->psNext)
5153 : {
5154 4872 : if (!SubstituteVariables(psIter, papszDict))
5155 0 : return false;
5156 : }
5157 4893 : return true;
5158 : }
5159 :
5160 : /************************************************************************/
5161 : /* GenerateMetadata() */
5162 : /************************************************************************/
5163 :
5164 21 : CPLString BAGCreator::GenerateMetadata(int nXSize, int nYSize,
5165 : const GDALGeoTransform >,
5166 : const OGRSpatialReference *poSRS,
5167 : char **papszOptions)
5168 : {
5169 : CPLXMLNode *psRoot;
5170 : CPLString osTemplateFilename =
5171 42 : CSLFetchNameValueDef(papszOptions, "TEMPLATE", "");
5172 21 : if (!osTemplateFilename.empty())
5173 : {
5174 0 : psRoot = CPLParseXMLFile(osTemplateFilename);
5175 : }
5176 : else
5177 : {
5178 : #ifndef USE_ONLY_EMBEDDED_RESOURCE_FILES
5179 : #ifdef EMBED_RESOURCE_FILES
5180 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
5181 : #endif
5182 : const char *pszDefaultTemplateFilename =
5183 21 : CPLFindFile("gdal", "bag_template.xml");
5184 21 : if (pszDefaultTemplateFilename == nullptr)
5185 : #endif
5186 : {
5187 : #ifdef EMBED_RESOURCE_FILES
5188 : static const bool bOnce [[maybe_unused]] = []()
5189 : {
5190 : CPLDebug("BAG", "Using embedded bag_template.xml");
5191 : return true;
5192 : }();
5193 : psRoot = CPLParseXMLString(BAGGetEmbeddedTemplateFile());
5194 : #else
5195 0 : CPLError(CE_Failure, CPLE_AppDefined,
5196 : "Cannot find bag_template.xml and TEMPLATE "
5197 : "creation option not specified");
5198 0 : return CPLString();
5199 : #endif
5200 : }
5201 : #ifndef USE_ONLY_EMBEDDED_RESOURCE_FILES
5202 : else
5203 : {
5204 21 : psRoot = CPLParseXMLFile(pszDefaultTemplateFilename);
5205 : }
5206 : #endif
5207 : }
5208 21 : if (psRoot == nullptr)
5209 0 : return CPLString();
5210 42 : CPLXMLTreeCloser oCloser(psRoot);
5211 21 : CPL_IGNORE_RET_VAL(oCloser);
5212 :
5213 21 : CPLXMLNode *psMain = psRoot;
5214 42 : for (; psMain; psMain = psMain->psNext)
5215 : {
5216 42 : if (psMain->eType == CXT_Element && !STARTS_WITH(psMain->pszValue, "?"))
5217 : {
5218 21 : break;
5219 : }
5220 : }
5221 21 : if (psMain == nullptr)
5222 : {
5223 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find main XML node");
5224 0 : return CPLString();
5225 : }
5226 :
5227 42 : CPLStringList osOptions(papszOptions, FALSE);
5228 21 : if (osOptions.FetchNameValue("VAR_PROCESS_STEP_DESCRIPTION") == nullptr)
5229 : {
5230 : osOptions.SetNameValue("VAR_PROCESS_STEP_DESCRIPTION",
5231 : CPLSPrintf("Generated by GDAL %s",
5232 21 : GDALVersionInfo("RELEASE_NAME")));
5233 : }
5234 21 : osOptions.SetNameValue("VAR_HEIGHT", CPLSPrintf("%d", nYSize));
5235 21 : osOptions.SetNameValue("VAR_WIDTH", CPLSPrintf("%d", nXSize));
5236 :
5237 : struct tm brokenDown;
5238 21 : CPLUnixTimeToYMDHMS(time(nullptr), &brokenDown);
5239 21 : if (osOptions.FetchNameValue("VAR_DATE") == nullptr)
5240 : {
5241 : osOptions.SetNameValue(
5242 21 : "VAR_DATE", CPLSPrintf("%04d-%02d-%02d", brokenDown.tm_year + 1900,
5243 21 : brokenDown.tm_mon + 1, brokenDown.tm_mday));
5244 : }
5245 21 : if (osOptions.FetchNameValue("VAR_DATETIME") == nullptr)
5246 : {
5247 : osOptions.SetNameValue(
5248 : "VAR_DATETIME",
5249 : CPLSPrintf("%04d-%02d-%02dT%02d:%02d:%02d",
5250 21 : brokenDown.tm_year + 1900, brokenDown.tm_mon + 1,
5251 : brokenDown.tm_mday, brokenDown.tm_hour,
5252 21 : brokenDown.tm_min, brokenDown.tm_sec));
5253 : }
5254 :
5255 21 : osOptions.SetNameValue("VAR_RESX", CPLSPrintf("%.17g", gt[1]));
5256 21 : osOptions.SetNameValue("VAR_RESY", CPLSPrintf("%.17g", fabs(gt[5])));
5257 : osOptions.SetNameValue("VAR_RES",
5258 21 : CPLSPrintf("%.17g", std::max(gt[1], fabs(gt[5]))));
5259 :
5260 21 : char *pszProjection = nullptr;
5261 21 : if (poSRS)
5262 21 : poSRS->exportToWkt(&pszProjection);
5263 21 : if (pszProjection == nullptr || EQUAL(pszProjection, ""))
5264 : {
5265 0 : CPLError(CE_Failure, CPLE_NotSupported,
5266 : "BAG driver requires a source dataset with a projection");
5267 : }
5268 21 : osOptions.SetNameValue("VAR_HORIZ_WKT", pszProjection);
5269 21 : CPLFree(pszProjection);
5270 :
5271 42 : OGRSpatialReference oSRS;
5272 21 : if (poSRS)
5273 21 : oSRS = *poSRS;
5274 21 : if (oSRS.IsCompound())
5275 : {
5276 2 : auto node = oSRS.GetRoot();
5277 2 : if (node && node->GetChildCount() == 3)
5278 : {
5279 2 : char *pszHorizWKT = nullptr;
5280 2 : node->GetChild(1)->exportToWkt(&pszHorizWKT);
5281 2 : char *pszVertWKT = nullptr;
5282 2 : node->GetChild(2)->exportToWkt(&pszVertWKT);
5283 :
5284 2 : oSRS.StripVertical();
5285 :
5286 2 : osOptions.SetNameValue("VAR_HORIZ_WKT", pszHorizWKT);
5287 2 : if (osOptions.FetchNameValue("VAR_VERT_WKT") == nullptr)
5288 : {
5289 2 : osOptions.SetNameValue("VAR_VERT_WKT", pszVertWKT);
5290 : }
5291 2 : CPLFree(pszHorizWKT);
5292 2 : CPLFree(pszVertWKT);
5293 : }
5294 : }
5295 :
5296 21 : const char *pszUnits = "m";
5297 21 : if (oSRS.IsProjected())
5298 : {
5299 6 : oSRS.GetLinearUnits(&pszUnits);
5300 6 : if (EQUAL(pszUnits, "metre"))
5301 6 : pszUnits = "m";
5302 : }
5303 : else
5304 : {
5305 15 : pszUnits = "deg";
5306 : }
5307 21 : osOptions.SetNameValue("VAR_RES_UNIT", pszUnits);
5308 :
5309 : // get bounds as pixel center
5310 21 : double dfMinX = gt[0] + gt[1] / 2;
5311 21 : double dfMaxX = dfMinX + (nXSize - 1) * gt[1];
5312 21 : double dfMaxY = gt[3] + gt[5] / 2;
5313 21 : double dfMinY = dfMaxY + (nYSize - 1) * gt[5];
5314 :
5315 21 : if (gt[5] > 0)
5316 : {
5317 1 : std::swap(dfMinY, dfMaxY);
5318 : }
5319 : osOptions.SetNameValue(
5320 : "VAR_CORNER_POINTS",
5321 21 : CPLSPrintf("%.17g,%.17g %.17g,%.17g", dfMinX, dfMinY, dfMaxX, dfMaxY));
5322 :
5323 21 : double adfCornerX[4] = {dfMinX, dfMinX, dfMaxX, dfMaxX};
5324 21 : double adfCornerY[4] = {dfMinY, dfMaxY, dfMaxY, dfMinY};
5325 42 : OGRSpatialReference oSRS_WGS84;
5326 21 : oSRS_WGS84.SetFromUserInput("WGS84");
5327 21 : oSRS_WGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
5328 : OGRCoordinateTransformation *poCT =
5329 21 : OGRCreateCoordinateTransformation(&oSRS, &oSRS_WGS84);
5330 21 : if (!poCT)
5331 0 : return CPLString();
5332 21 : if (!poCT->Transform(4, adfCornerX, adfCornerY))
5333 : {
5334 0 : CPLError(CE_Failure, CPLE_AppDefined,
5335 : "Cannot compute raster extent in geodetic coordinates");
5336 0 : delete poCT;
5337 0 : return CPLString();
5338 : }
5339 21 : delete poCT;
5340 : double dfWest = std::min(std::min(adfCornerX[0], adfCornerX[1]),
5341 21 : std::min(adfCornerX[2], adfCornerX[3]));
5342 : double dfSouth = std::min(std::min(adfCornerY[0], adfCornerY[1]),
5343 21 : std::min(adfCornerY[2], adfCornerY[3]));
5344 : double dfEast = std::max(std::max(adfCornerX[0], adfCornerX[1]),
5345 21 : std::max(adfCornerX[2], adfCornerX[3]));
5346 : double dfNorth = std::max(std::max(adfCornerY[0], adfCornerY[1]),
5347 21 : std::max(adfCornerY[2], adfCornerY[3]));
5348 21 : osOptions.SetNameValue("VAR_WEST_LONGITUDE", CPLSPrintf("%.17g", dfWest));
5349 21 : osOptions.SetNameValue("VAR_SOUTH_LATITUDE", CPLSPrintf("%.17g", dfSouth));
5350 21 : osOptions.SetNameValue("VAR_EAST_LONGITUDE", CPLSPrintf("%.17g", dfEast));
5351 21 : osOptions.SetNameValue("VAR_NORTH_LATITUDE", CPLSPrintf("%.17g", dfNorth));
5352 :
5353 21 : if (!SubstituteVariables(psMain, osOptions.List()))
5354 : {
5355 0 : return CPLString();
5356 : }
5357 :
5358 21 : char *pszXML = CPLSerializeXMLTree(psRoot);
5359 42 : CPLString osXML(pszXML);
5360 21 : CPLFree(pszXML);
5361 21 : return osXML;
5362 : }
5363 :
5364 : /************************************************************************/
5365 : /* CreateAndWriteMetadata() */
5366 : /************************************************************************/
5367 :
5368 19 : bool BAGCreator::CreateAndWriteMetadata(hid_t hdf5,
5369 : const CPLString &osXMLMetadata)
5370 : {
5371 19 : hsize_t dim_init[1] = {1 + osXMLMetadata.size()};
5372 19 : hsize_t dim_max[1] = {H5S_UNLIMITED};
5373 :
5374 19 : hid_t hDataSpace = H5_CHECK(H5Screate_simple(1, dim_init, dim_max));
5375 19 : if (hDataSpace < 0)
5376 0 : return false;
5377 :
5378 19 : hid_t hParams = -1;
5379 19 : hid_t hDataType = -1;
5380 19 : hid_t hDatasetID = -1;
5381 19 : hid_t hFileSpace = -1;
5382 19 : bool ret = false;
5383 : do
5384 : {
5385 19 : hParams = H5_CHECK(H5Pcreate(H5P_DATASET_CREATE));
5386 19 : if (hParams < 0)
5387 0 : break;
5388 :
5389 19 : hsize_t chunk_dims[1] = {1024};
5390 19 : if (H5_CHECK(H5Pset_chunk(hParams, 1, chunk_dims)) < 0)
5391 0 : break;
5392 :
5393 19 : hDataType = H5_CHECK(H5Tcopy(H5T_C_S1));
5394 19 : if (hDataType < 0)
5395 0 : break;
5396 :
5397 19 : hDatasetID = H5_CHECK(H5Dcreate(hdf5, "/BAG_root/metadata", hDataType,
5398 : hDataSpace, hParams));
5399 19 : if (hDatasetID < 0)
5400 0 : break;
5401 :
5402 19 : if (H5_CHECK(H5Dextend(hDatasetID, dim_init)) < 0)
5403 0 : break;
5404 :
5405 19 : hFileSpace = H5_CHECK(H5Dget_space(hDatasetID));
5406 19 : if (hFileSpace < 0)
5407 0 : break;
5408 :
5409 19 : H5OFFSET_TYPE offset[1] = {0};
5410 19 : if (H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
5411 19 : nullptr, dim_init, nullptr)) < 0)
5412 : {
5413 0 : break;
5414 : }
5415 :
5416 19 : if (H5_CHECK(H5Dwrite(hDatasetID, hDataType, hDataSpace, hFileSpace,
5417 19 : H5P_DEFAULT, osXMLMetadata.data())) < 0)
5418 : {
5419 0 : break;
5420 : }
5421 :
5422 19 : ret = true;
5423 : } while (0);
5424 :
5425 19 : if (hParams >= 0)
5426 19 : H5_CHECK(H5Pclose(hParams));
5427 19 : if (hDataType >= 0)
5428 19 : H5_CHECK(H5Tclose(hDataType));
5429 19 : if (hFileSpace >= 0)
5430 19 : H5_CHECK(H5Sclose(hFileSpace));
5431 19 : if (hDatasetID >= 0)
5432 19 : H5_CHECK(H5Dclose(hDatasetID));
5433 19 : H5_CHECK(H5Sclose(hDataSpace));
5434 :
5435 19 : return ret;
5436 : }
5437 :
5438 : /************************************************************************/
5439 : /* CreateTrackingListDataset() */
5440 : /************************************************************************/
5441 :
5442 20 : bool BAGCreator::CreateTrackingListDataset()
5443 : {
5444 : typedef struct
5445 : {
5446 : uint32_t row;
5447 : uint32_t col;
5448 : float depth;
5449 : float uncertainty;
5450 : uint8_t track_code;
5451 : uint16_t list_series;
5452 : } TrackingListItem;
5453 :
5454 20 : hsize_t dim_init[1] = {0};
5455 20 : hsize_t dim_max[1] = {H5S_UNLIMITED};
5456 :
5457 20 : hid_t hDataSpace = H5_CHECK(H5Screate_simple(1, dim_init, dim_max));
5458 20 : if (hDataSpace < 0)
5459 0 : return false;
5460 :
5461 20 : hid_t hParams = -1;
5462 20 : hid_t hDataType = -1;
5463 20 : hid_t hDatasetID = -1;
5464 20 : bool ret = false;
5465 : do
5466 : {
5467 20 : hParams = H5_CHECK(H5Pcreate(H5P_DATASET_CREATE));
5468 20 : if (hParams < 0)
5469 0 : break;
5470 :
5471 20 : hsize_t chunk_dims[1] = {10};
5472 20 : if (H5_CHECK(H5Pset_chunk(hParams, 1, chunk_dims)) < 0)
5473 0 : break;
5474 :
5475 : TrackingListItem unusedItem;
5476 : (void)unusedItem.row;
5477 : (void)unusedItem.col;
5478 : (void)unusedItem.depth;
5479 : (void)unusedItem.uncertainty;
5480 : (void)unusedItem.track_code;
5481 : (void)unusedItem.list_series;
5482 20 : hDataType = H5_CHECK(H5Tcreate(H5T_COMPOUND, sizeof(unusedItem)));
5483 20 : if (hDataType < 0)
5484 0 : break;
5485 :
5486 20 : if (H5Tinsert(hDataType, "row", HOFFSET(TrackingListItem, row),
5487 20 : H5T_NATIVE_UINT) < 0 ||
5488 20 : H5Tinsert(hDataType, "col", HOFFSET(TrackingListItem, col),
5489 20 : H5T_NATIVE_UINT) < 0 ||
5490 20 : H5Tinsert(hDataType, "depth", HOFFSET(TrackingListItem, depth),
5491 20 : H5T_NATIVE_FLOAT) < 0 ||
5492 20 : H5Tinsert(hDataType, "uncertainty",
5493 : HOFFSET(TrackingListItem, uncertainty),
5494 20 : H5T_NATIVE_FLOAT) < 0 ||
5495 20 : H5Tinsert(hDataType, "track_code",
5496 : HOFFSET(TrackingListItem, track_code),
5497 60 : H5T_NATIVE_UCHAR) < 0 ||
5498 20 : H5Tinsert(hDataType, "list_series",
5499 : HOFFSET(TrackingListItem, list_series),
5500 20 : H5T_NATIVE_SHORT) < 0)
5501 : {
5502 0 : break;
5503 : }
5504 :
5505 20 : hDatasetID = H5_CHECK(H5Dcreate(m_hdf5, "/BAG_root/tracking_list",
5506 : hDataType, hDataSpace, hParams));
5507 20 : if (hDatasetID < 0)
5508 0 : break;
5509 :
5510 20 : if (H5_CHECK(H5Dextend(hDatasetID, dim_init)) < 0)
5511 0 : break;
5512 :
5513 20 : if (!GH5_CreateAttribute(hDatasetID, "Tracking List Length",
5514 20 : H5T_NATIVE_UINT))
5515 0 : break;
5516 :
5517 20 : if (!GH5_WriteAttribute(hDatasetID, "Tracking List Length", 0U))
5518 0 : break;
5519 :
5520 20 : ret = true;
5521 : } while (0);
5522 :
5523 20 : if (hParams >= 0)
5524 20 : H5_CHECK(H5Pclose(hParams));
5525 20 : if (hDataType >= 0)
5526 20 : H5_CHECK(H5Tclose(hDataType));
5527 20 : if (hDatasetID >= 0)
5528 20 : H5_CHECK(H5Dclose(hDatasetID));
5529 20 : H5_CHECK(H5Sclose(hDataSpace));
5530 :
5531 20 : return ret;
5532 : }
5533 :
5534 : /************************************************************************/
5535 : /* CreateElevationOrUncertainty() */
5536 : /************************************************************************/
5537 :
5538 32 : bool BAGCreator::CreateElevationOrUncertainty(
5539 : GDALDataset *poSrcDS, int nBand, const char *pszDSName,
5540 : const char *pszMaxAttrName, const char *pszMinAttrName, char **papszOptions,
5541 : GDALProgressFunc pfnProgress, void *pProgressData)
5542 : {
5543 32 : const int nYSize = poSrcDS->GetRasterYSize();
5544 32 : const int nXSize = poSrcDS->GetRasterXSize();
5545 :
5546 32 : GDALGeoTransform gt;
5547 32 : poSrcDS->GetGeoTransform(gt);
5548 :
5549 32 : hsize_t dims[2] = {static_cast<hsize_t>(nYSize),
5550 32 : static_cast<hsize_t>(nXSize)};
5551 :
5552 32 : hid_t hDataSpace = H5_CHECK(H5Screate_simple(2, dims, nullptr));
5553 32 : if (hDataSpace < 0)
5554 0 : return false;
5555 :
5556 32 : hid_t hParams = -1;
5557 32 : hid_t hDataType = -1;
5558 32 : hid_t hDatasetID = -1;
5559 32 : hid_t hFileSpace = -1;
5560 32 : bool bDeflate = EQUAL(
5561 : CSLFetchNameValueDef(papszOptions, "COMPRESS", "DEFLATE"), "DEFLATE");
5562 : int nCompressionLevel =
5563 32 : atoi(CSLFetchNameValueDef(papszOptions, "ZLEVEL", "6"));
5564 32 : const int nBlockSize = std::min(
5565 32 : 4096, atoi(CSLFetchNameValueDef(papszOptions, "BLOCK_SIZE", "100")));
5566 32 : const int nBlockXSize = std::min(nXSize, nBlockSize);
5567 32 : const int nBlockYSize = std::min(nYSize, nBlockSize);
5568 32 : bool ret = false;
5569 32 : const float fNoDataValue = fDEFAULT_NODATA;
5570 : do
5571 : {
5572 32 : hDataType = H5_CHECK(H5Tcopy(H5T_NATIVE_FLOAT));
5573 32 : if (hDataType < 0)
5574 0 : break;
5575 :
5576 32 : if (H5_CHECK(H5Tset_order(hDataType, H5T_ORDER_LE)) < 0)
5577 0 : break;
5578 :
5579 32 : hParams = H5_CHECK(H5Pcreate(H5P_DATASET_CREATE));
5580 32 : if (hParams < 0)
5581 0 : break;
5582 :
5583 32 : if (H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) < 0)
5584 0 : break;
5585 :
5586 32 : if (H5_CHECK(H5Pset_fill_value(hParams, hDataType, &fNoDataValue)) < 0)
5587 0 : break;
5588 :
5589 32 : if (H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) < 0)
5590 0 : break;
5591 32 : hsize_t chunk_size[2] = {static_cast<hsize_t>(nBlockYSize),
5592 32 : static_cast<hsize_t>(nBlockXSize)};
5593 32 : if (H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) < 0)
5594 0 : break;
5595 :
5596 32 : if (bDeflate)
5597 : {
5598 32 : if (H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) < 0)
5599 0 : break;
5600 : }
5601 :
5602 32 : hDatasetID = H5_CHECK(
5603 : H5Dcreate(m_hdf5, pszDSName, hDataType, hDataSpace, hParams));
5604 32 : if (hDatasetID < 0)
5605 0 : break;
5606 :
5607 32 : if (!GH5_CreateAttribute(hDatasetID, pszMaxAttrName, hDataType))
5608 0 : break;
5609 :
5610 32 : if (!GH5_CreateAttribute(hDatasetID, pszMinAttrName, hDataType))
5611 0 : break;
5612 :
5613 32 : hFileSpace = H5_CHECK(H5Dget_space(hDatasetID));
5614 32 : if (hFileSpace < 0)
5615 0 : break;
5616 :
5617 32 : const int nYBlocks =
5618 32 : static_cast<int>(DIV_ROUND_UP(nYSize, nBlockYSize));
5619 32 : const int nXBlocks =
5620 32 : static_cast<int>(DIV_ROUND_UP(nXSize, nBlockXSize));
5621 32 : std::vector<float> afValues(static_cast<size_t>(nBlockYSize) *
5622 32 : nBlockXSize);
5623 32 : ret = true;
5624 32 : const bool bReverseY = gt[5] < 0;
5625 :
5626 32 : float fMin = std::numeric_limits<float>::infinity();
5627 32 : float fMax = -std::numeric_limits<float>::infinity();
5628 :
5629 32 : if (nBand == 1 || poSrcDS->GetRasterCount() == 2)
5630 : {
5631 18 : int bHasNoData = FALSE;
5632 : const double dfSrcNoData =
5633 18 : poSrcDS->GetRasterBand(nBand)->GetNoDataValue(&bHasNoData);
5634 18 : const float fSrcNoData = static_cast<float>(dfSrcNoData);
5635 :
5636 38 : for (int iY = 0; ret && iY < nYBlocks; iY++)
5637 : {
5638 : const int nSrcYOff =
5639 20 : bReverseY ? std::max(0, nYSize - (iY + 1) * nBlockYSize)
5640 1 : : iY * nBlockYSize;
5641 : const int nReqCountY =
5642 20 : std::min(nBlockYSize, nYSize - iY * nBlockYSize);
5643 48 : for (int iX = 0; iX < nXBlocks; iX++)
5644 : {
5645 : const int nReqCountX =
5646 28 : std::min(nBlockXSize, nXSize - iX * nBlockXSize);
5647 :
5648 55 : if (poSrcDS->GetRasterBand(nBand)->RasterIO(
5649 28 : GF_Read, iX * nBlockXSize, nSrcYOff, nReqCountX,
5650 : nReqCountY,
5651 27 : bReverseY ? afValues.data() +
5652 27 : (nReqCountY - 1) * nReqCountX
5653 1 : : afValues.data(),
5654 : nReqCountX, nReqCountY, GDT_Float32, 0,
5655 27 : bReverseY ? -4 * nReqCountX : 0,
5656 28 : nullptr) != CE_None)
5657 : {
5658 0 : ret = false;
5659 0 : break;
5660 : }
5661 :
5662 2376 : for (int i = 0; i < nReqCountY * nReqCountX; i++)
5663 : {
5664 2348 : const float fVal = afValues[i];
5665 4696 : if ((bHasNoData && fVal == fSrcNoData) ||
5666 2348 : std::isnan(fVal))
5667 : {
5668 0 : afValues[i] = fNoDataValue;
5669 : }
5670 : else
5671 : {
5672 2348 : fMin = std::min(fMin, fVal);
5673 2348 : fMax = std::max(fMax, fVal);
5674 : }
5675 : }
5676 :
5677 : H5OFFSET_TYPE offset[2] = {
5678 28 : static_cast<H5OFFSET_TYPE>(iY) *
5679 28 : static_cast<H5OFFSET_TYPE>(nBlockYSize),
5680 28 : static_cast<H5OFFSET_TYPE>(iX) *
5681 28 : static_cast<H5OFFSET_TYPE>(nBlockXSize)};
5682 28 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
5683 28 : static_cast<hsize_t>(nReqCountX)};
5684 28 : if (H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET,
5685 : offset, nullptr, count,
5686 28 : nullptr)) < 0)
5687 : {
5688 0 : ret = false;
5689 0 : break;
5690 : }
5691 :
5692 28 : hid_t hMemSpace = H5Screate_simple(2, count, nullptr);
5693 28 : if (hMemSpace < 0)
5694 0 : break;
5695 :
5696 28 : if (H5_CHECK(H5Dwrite(hDatasetID, H5T_NATIVE_FLOAT,
5697 : hMemSpace, hFileSpace, H5P_DEFAULT,
5698 28 : afValues.data())) < 0)
5699 : {
5700 0 : H5Sclose(hMemSpace);
5701 0 : ret = false;
5702 0 : break;
5703 : }
5704 :
5705 28 : H5Sclose(hMemSpace);
5706 :
5707 28 : if (!pfnProgress(
5708 28 : (static_cast<double>(iY) * nXBlocks + iX + 1) /
5709 28 : (static_cast<double>(nXBlocks) * nYBlocks),
5710 : "", pProgressData))
5711 : {
5712 0 : ret = false;
5713 0 : break;
5714 : }
5715 : }
5716 : }
5717 : }
5718 32 : if (!ret)
5719 0 : break;
5720 :
5721 32 : if (fMin > fMax)
5722 14 : fMin = fMax = fNoDataValue;
5723 :
5724 32 : if (!GH5_WriteAttribute(hDatasetID, pszMaxAttrName, fMax))
5725 0 : break;
5726 :
5727 32 : if (!GH5_WriteAttribute(hDatasetID, pszMinAttrName, fMin))
5728 0 : break;
5729 :
5730 32 : ret = true;
5731 : } while (0);
5732 :
5733 32 : if (hParams >= 0)
5734 32 : H5_CHECK(H5Pclose(hParams));
5735 32 : if (hDataType >= 0)
5736 32 : H5_CHECK(H5Tclose(hDataType));
5737 32 : if (hFileSpace >= 0)
5738 32 : H5_CHECK(H5Sclose(hFileSpace));
5739 32 : if (hDatasetID >= 0)
5740 32 : H5_CHECK(H5Dclose(hDatasetID));
5741 32 : H5_CHECK(H5Sclose(hDataSpace));
5742 :
5743 32 : return ret;
5744 : }
5745 :
5746 : /************************************************************************/
5747 : /* CreateBase() */
5748 : /************************************************************************/
5749 :
5750 22 : bool BAGCreator::CreateBase(const char *pszFilename, char **papszOptions)
5751 : {
5752 22 : hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
5753 22 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
5754 22 : m_hdf5 = H5Fcreate(pszFilename, H5F_ACC_TRUNC, H5P_DEFAULT, fapl);
5755 22 : H5Pclose(fapl);
5756 22 : if (m_hdf5 < 0)
5757 : {
5758 2 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create file");
5759 2 : return false;
5760 : }
5761 :
5762 20 : m_bagRoot = H5_CHECK(H5Gcreate(m_hdf5, "/BAG_root", 0));
5763 20 : if (m_bagRoot < 0)
5764 : {
5765 0 : return false;
5766 : }
5767 :
5768 : const char *pszVersion =
5769 20 : CSLFetchNameValueDef(papszOptions, "BAG_VERSION", "1.6.2");
5770 20 : constexpr unsigned knVersionLength = 32;
5771 20 : char szVersion[knVersionLength] = {};
5772 20 : snprintf(szVersion, sizeof(szVersion), "%s", pszVersion);
5773 20 : if (!GH5_CreateAttribute(m_bagRoot, "Bag Version", H5T_C_S1,
5774 40 : knVersionLength) ||
5775 20 : !GH5_WriteAttribute(m_bagRoot, "Bag Version", szVersion))
5776 : {
5777 0 : return false;
5778 : }
5779 :
5780 20 : if (!CreateTrackingListDataset())
5781 : {
5782 0 : return false;
5783 : }
5784 20 : return true;
5785 : }
5786 :
5787 : /************************************************************************/
5788 : /* Create() */
5789 : /************************************************************************/
5790 :
5791 22 : bool BAGCreator::Create(const char *pszFilename, GDALDataset *poSrcDS,
5792 : char **papszOptions, GDALProgressFunc pfnProgress,
5793 : void *pProgressData)
5794 : {
5795 22 : const int nBands = poSrcDS->GetRasterCount();
5796 22 : if (nBands != 1 && nBands != 2)
5797 : {
5798 4 : CPLError(CE_Failure, CPLE_NotSupported,
5799 : "BAG driver doesn't support %d bands. Must be 1 or 2.",
5800 : nBands);
5801 4 : return false;
5802 : }
5803 18 : GDALGeoTransform gt;
5804 18 : if (poSrcDS->GetGeoTransform(gt) != CE_None)
5805 : {
5806 0 : CPLError(CE_Failure, CPLE_NotSupported,
5807 : "BAG driver requires a source dataset with a geotransform");
5808 0 : return false;
5809 : }
5810 18 : if (gt[2] != 0 || gt[4] != 0)
5811 : {
5812 0 : CPLError(CE_Failure, CPLE_NotSupported,
5813 : "BAG driver requires a source dataset with a non-rotated "
5814 : "geotransform");
5815 0 : return false;
5816 : }
5817 :
5818 : CPLString osXMLMetadata =
5819 : GenerateMetadata(poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
5820 36 : gt, poSrcDS->GetSpatialRef(), papszOptions);
5821 18 : if (osXMLMetadata.empty())
5822 : {
5823 0 : return false;
5824 : }
5825 :
5826 18 : if (!CreateBase(pszFilename, papszOptions))
5827 : {
5828 2 : return false;
5829 : }
5830 :
5831 16 : if (!CreateAndWriteMetadata(m_hdf5, osXMLMetadata))
5832 : {
5833 0 : return false;
5834 : }
5835 :
5836 16 : void *pScaled = GDALCreateScaledProgress(0, 1. / poSrcDS->GetRasterCount(),
5837 : pfnProgress, pProgressData);
5838 : bool bRet;
5839 16 : bRet = CreateElevationOrUncertainty(
5840 : poSrcDS, 1, "/BAG_root/elevation", "Maximum Elevation Value",
5841 : "Minimum Elevation Value", papszOptions, GDALScaledProgress, pScaled);
5842 16 : GDALDestroyScaledProgress(pScaled);
5843 16 : if (!bRet)
5844 : {
5845 0 : return false;
5846 : }
5847 :
5848 16 : pScaled = GDALCreateScaledProgress(1. / poSrcDS->GetRasterCount(), 1.0,
5849 : pfnProgress, pProgressData);
5850 16 : bRet = CreateElevationOrUncertainty(
5851 : poSrcDS, 2, "/BAG_root/uncertainty", "Maximum Uncertainty Value",
5852 : "Minimum Uncertainty Value", papszOptions, GDALScaledProgress, pScaled);
5853 16 : GDALDestroyScaledProgress(pScaled);
5854 16 : if (!bRet)
5855 : {
5856 0 : return false;
5857 : }
5858 :
5859 16 : return Close();
5860 : }
5861 :
5862 : /************************************************************************/
5863 : /* Create() */
5864 : /************************************************************************/
5865 :
5866 52 : bool BAGCreator::Create(const char *pszFilename, int nBands, GDALDataType eType,
5867 : char **papszOptions)
5868 : {
5869 52 : if (nBands != 1 && nBands != 2)
5870 : {
5871 34 : CPLError(CE_Failure, CPLE_NotSupported,
5872 : "BAG driver doesn't support %d bands. Must be 1 or 2.",
5873 : nBands);
5874 34 : return false;
5875 : }
5876 18 : if (eType != GDT_Float32)
5877 : {
5878 14 : CPLError(CE_Failure, CPLE_NotSupported,
5879 : "BAG driver only supports Float32");
5880 14 : return false;
5881 : }
5882 :
5883 4 : if (!CreateBase(pszFilename, papszOptions))
5884 : {
5885 0 : return false;
5886 : }
5887 :
5888 4 : return Close();
5889 : }
5890 :
5891 : /************************************************************************/
5892 : /* CreateCopy() */
5893 : /************************************************************************/
5894 :
5895 22 : GDALDataset *BAGDataset::CreateCopy(const char *pszFilename,
5896 : GDALDataset *poSrcDS, int /* bStrict */,
5897 : char **papszOptions,
5898 : GDALProgressFunc pfnProgress,
5899 : void *pProgressData)
5900 :
5901 : {
5902 22 : if (!BAGCreator().Create(pszFilename, poSrcDS, papszOptions, pfnProgress,
5903 : pProgressData))
5904 : {
5905 6 : return nullptr;
5906 : }
5907 :
5908 32 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
5909 16 : oOpenInfo.nOpenFlags = GDAL_OF_RASTER;
5910 16 : return Open(&oOpenInfo);
5911 : }
5912 :
5913 : /************************************************************************/
5914 : /* Create() */
5915 : /************************************************************************/
5916 :
5917 52 : GDALDataset *BAGDataset::Create(const char *pszFilename, int nXSize, int nYSize,
5918 : int nBandsIn, GDALDataType eType,
5919 : char **papszOptions)
5920 : {
5921 52 : if (!BAGCreator().Create(pszFilename, nBandsIn, eType, papszOptions))
5922 : {
5923 48 : return nullptr;
5924 : }
5925 :
5926 8 : GDALOpenInfo oOpenInfo(pszFilename, GA_Update);
5927 4 : oOpenInfo.nOpenFlags = GDAL_OF_RASTER;
5928 4 : return OpenForCreate(&oOpenInfo, nXSize, nYSize, nBandsIn, papszOptions);
5929 : }
5930 :
5931 : /************************************************************************/
5932 : /* GDALRegister_BAG() */
5933 : /************************************************************************/
5934 14 : void GDALRegister_BAG()
5935 :
5936 : {
5937 14 : if (!GDAL_CHECK_VERSION("BAG"))
5938 0 : return;
5939 :
5940 14 : if (GDALGetDriverByName(BAG_DRIVER_NAME) != nullptr)
5941 0 : return;
5942 :
5943 14 : GDALDriver *poDriver = new GDALDriver();
5944 :
5945 14 : BAGDriverSetCommonMetadata(poDriver);
5946 :
5947 14 : poDriver->pfnOpen = BAGDataset::Open;
5948 14 : poDriver->pfnUnloadDriver = BAGDatasetDriverUnload;
5949 14 : poDriver->pfnCreateCopy = BAGDataset::CreateCopy;
5950 14 : poDriver->pfnCreate = BAGDataset::Create;
5951 :
5952 14 : GetGDALDriverManager()->RegisterDriver(poDriver);
5953 : }
|
