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