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