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