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