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