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