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