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