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