Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Hierarchical Data Format Release 5 (HDF5)
4 : * Purpose: Read S102 bathymetric datasets.
5 : * Author: Even Rouault <even dot rouault at spatialys dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2023-2025, Even Rouault <even dot rouault at spatialys dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_port.h"
14 : #include "cpl_vsi.h"
15 :
16 : #include "hdf5dataset.h"
17 : #include "hdf5drivercore.h"
18 : #include "gh5_convenience.h"
19 : #include "rat.h"
20 : #include "s100.h"
21 :
22 : #include "gdal_frmts.h"
23 : #include "gdal_priv.h"
24 : #include "gdal_proxy.h"
25 : #include "gdal_rat.h"
26 :
27 : #include <algorithm>
28 : #include <array>
29 : #include <cmath>
30 : #include <limits>
31 : #include <map>
32 : #include <set>
33 : #include <string_view>
34 :
35 : /************************************************************************/
36 : /* S102Dataset */
37 : /************************************************************************/
38 :
39 98 : class S102Dataset final : public S100BaseDataset
40 : {
41 : bool OpenQuality(GDALOpenInfo *poOpenInfo,
42 : const std::shared_ptr<GDALGroup> &poRootGroup);
43 :
44 : public:
45 49 : explicit S102Dataset(const std::string &osFilename)
46 49 : : S100BaseDataset(osFilename)
47 : {
48 49 : }
49 :
50 : ~S102Dataset() override;
51 :
52 : static GDALDataset *Open(GDALOpenInfo *);
53 : static GDALDataset *CreateCopy(const char *pszFilename,
54 : GDALDataset *poSrcDS, int bStrict,
55 : char **papszOptions,
56 : GDALProgressFunc pfnProgress,
57 : void *pProgressData);
58 : };
59 :
60 : S102Dataset::~S102Dataset() = default;
61 :
62 : /************************************************************************/
63 : /* S102RasterBand */
64 : /************************************************************************/
65 :
66 : class S102RasterBand final : public GDALProxyRasterBand
67 : {
68 : friend class S102Dataset;
69 : std::unique_ptr<GDALDataset> m_poDS{};
70 : GDALRasterBand *m_poUnderlyingBand = nullptr;
71 : double m_dfMinimum = std::numeric_limits<double>::quiet_NaN();
72 : double m_dfMaximum = std::numeric_limits<double>::quiet_NaN();
73 :
74 : CPL_DISALLOW_COPY_ASSIGN(S102RasterBand)
75 :
76 : public:
77 55 : explicit S102RasterBand(std::unique_ptr<GDALDataset> &&poDSIn)
78 110 : : m_poDS(std::move(poDSIn)),
79 55 : m_poUnderlyingBand(m_poDS->GetRasterBand(1))
80 : {
81 55 : eDataType = m_poUnderlyingBand->GetRasterDataType();
82 55 : m_poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
83 55 : }
84 :
85 : GDALRasterBand *
86 : RefUnderlyingRasterBand(bool /*bForceOpen*/ = true) const override;
87 :
88 10 : double GetMinimum(int *pbSuccess = nullptr) override
89 : {
90 10 : if (pbSuccess)
91 10 : *pbSuccess = !std::isnan(m_dfMinimum);
92 10 : return m_dfMinimum;
93 : }
94 :
95 10 : double GetMaximum(int *pbSuccess = nullptr) override
96 : {
97 10 : if (pbSuccess)
98 10 : *pbSuccess = !std::isnan(m_dfMaximum);
99 10 : return m_dfMaximum;
100 : }
101 :
102 4 : const char *GetUnitType() override
103 : {
104 4 : return "metre";
105 : }
106 : };
107 :
108 : GDALRasterBand *
109 49 : S102RasterBand::RefUnderlyingRasterBand(bool /*bForceOpen*/) const
110 : {
111 49 : return m_poUnderlyingBand;
112 : }
113 :
114 : /************************************************************************/
115 : /* S102GeoreferencedMetadataRasterBand */
116 : /************************************************************************/
117 :
118 : class S102GeoreferencedMetadataRasterBand final : public GDALProxyRasterBand
119 : {
120 : friend class S102Dataset;
121 :
122 : std::unique_ptr<GDALDataset> m_poDS{};
123 : GDALRasterBand *m_poUnderlyingBand = nullptr;
124 : std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
125 :
126 : CPL_DISALLOW_COPY_ASSIGN(S102GeoreferencedMetadataRasterBand)
127 :
128 : public:
129 7 : explicit S102GeoreferencedMetadataRasterBand(
130 : std::unique_ptr<GDALDataset> &&poDSIn,
131 : std::unique_ptr<GDALRasterAttributeTable> &&poRAT)
132 14 : : m_poDS(std::move(poDSIn)),
133 7 : m_poUnderlyingBand(m_poDS->GetRasterBand(1)),
134 14 : m_poRAT(std::move(poRAT))
135 : {
136 7 : eDataType = m_poUnderlyingBand->GetRasterDataType();
137 7 : m_poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
138 7 : }
139 :
140 : GDALRasterBand *
141 : RefUnderlyingRasterBand(bool /*bForceOpen*/ = true) const override;
142 :
143 6 : GDALRasterAttributeTable *GetDefaultRAT() override
144 : {
145 6 : return m_poRAT.get();
146 : }
147 : };
148 :
149 43 : GDALRasterBand *S102GeoreferencedMetadataRasterBand::RefUnderlyingRasterBand(
150 : bool /*bForceOpen*/) const
151 : {
152 43 : return m_poUnderlyingBand;
153 : }
154 :
155 : /************************************************************************/
156 : /* Open() */
157 : /************************************************************************/
158 :
159 54 : GDALDataset *S102Dataset::Open(GDALOpenInfo *poOpenInfo)
160 :
161 : {
162 : // Confirm that this appears to be a S102 file.
163 54 : if (!S102DatasetIdentify(poOpenInfo))
164 1 : return nullptr;
165 :
166 : HDF5_GLOBAL_LOCK();
167 :
168 53 : if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
169 : {
170 2 : return HDF5Dataset::OpenMultiDim(poOpenInfo);
171 : }
172 :
173 : // Confirm the requested access is supported.
174 51 : if (poOpenInfo->eAccess == GA_Update)
175 : {
176 0 : ReportUpdateNotSupportedByDriver("S102");
177 0 : return nullptr;
178 : }
179 :
180 102 : std::string osFilename(poOpenInfo->pszFilename);
181 51 : bool bIsSubdataset = false;
182 51 : bool bIsQuality = false;
183 102 : std::string osBathymetryCoverageName = "BathymetryCoverage.01";
184 51 : if (STARTS_WITH(poOpenInfo->pszFilename, "S102:"))
185 : {
186 : const CPLStringList aosTokens(
187 19 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
188 19 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
189 :
190 19 : if (aosTokens.size() == 2)
191 : {
192 1 : osFilename = aosTokens[1];
193 : }
194 18 : else if (aosTokens.size() == 3)
195 : {
196 18 : bIsSubdataset = true;
197 18 : osFilename = aosTokens[1];
198 18 : if (EQUAL(aosTokens[2], "BathymetryCoverage"))
199 : {
200 : // Default dataset
201 : }
202 16 : else if (STARTS_WITH(aosTokens[2], "BathymetryCoverage"))
203 : {
204 5 : osBathymetryCoverageName = aosTokens[2];
205 : }
206 18 : else if (EQUAL(aosTokens[2], "QualityOfSurvey") || // < v3
207 7 : EQUAL(aosTokens[2], "QualityOfBathymetryCoverage")) // v3
208 : {
209 9 : bIsQuality = true;
210 : }
211 : else
212 : {
213 2 : CPLError(CE_Failure, CPLE_NotSupported,
214 : "Unsupported subdataset component: '%s'. Expected "
215 : "'QualityOfSurvey'",
216 : aosTokens[2]);
217 2 : return nullptr;
218 : }
219 : }
220 : else
221 : {
222 0 : return nullptr;
223 : }
224 : }
225 :
226 98 : auto poDS = std::make_unique<S102Dataset>(osFilename);
227 49 : if (!poDS->Init())
228 0 : return nullptr;
229 :
230 49 : const auto &poRootGroup = poDS->m_poRootGroup;
231 :
232 147 : auto poBathymetryCoverage = poRootGroup->OpenGroup("BathymetryCoverage");
233 49 : if (!poBathymetryCoverage)
234 : {
235 0 : CPLError(CE_Failure, CPLE_AppDefined,
236 : "S102: Cannot find /BathymetryCoverage group");
237 0 : return nullptr;
238 : }
239 :
240 49 : if (!bIsSubdataset)
241 : {
242 : auto poNumInstances =
243 66 : poBathymetryCoverage->GetAttribute("numInstances");
244 52 : if (poNumInstances &&
245 52 : poNumInstances->GetDataType().GetClass() == GEDTC_NUMERIC)
246 : {
247 19 : const int nNumInstances = poNumInstances->ReadAsInt();
248 19 : if (nNumInstances != 1)
249 : {
250 4 : CPLStringList aosSubDSList;
251 2 : int iSubDS = 0;
252 4 : for (const std::string &coverageName :
253 10 : poBathymetryCoverage->GetGroupNames())
254 : {
255 : auto poCoverage =
256 8 : poBathymetryCoverage->OpenGroup(coverageName);
257 4 : if (poCoverage)
258 : {
259 8 : GDALMajorObject mo;
260 : // Read first vertical datum from root group and let the
261 : // coverage override it.
262 4 : S100ReadVerticalDatum(&mo, poRootGroup.get());
263 4 : S100ReadVerticalDatum(&mo, poCoverage.get());
264 4 : ++iSubDS;
265 : aosSubDSList.SetNameValue(
266 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
267 : CPLSPrintf("S102:\"%s\":%s", osFilename.c_str(),
268 4 : coverageName.c_str()));
269 8 : std::string verticalDatum;
270 : const char *pszValue =
271 4 : mo.GetMetadataItem(S100_VERTICAL_DATUM_NAME);
272 4 : if (pszValue)
273 : {
274 4 : verticalDatum = ", vertical datum ";
275 4 : verticalDatum += pszValue;
276 : pszValue =
277 4 : mo.GetMetadataItem(S100_VERTICAL_DATUM_ABBREV);
278 4 : if (pszValue)
279 : {
280 4 : verticalDatum += " (";
281 4 : verticalDatum += pszValue;
282 4 : verticalDatum += ')';
283 : }
284 : }
285 : aosSubDSList.SetNameValue(
286 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
287 : CPLSPrintf(
288 : "Bathymetric gridded data, instance %s%s",
289 4 : coverageName.c_str(), verticalDatum.c_str()));
290 : }
291 : }
292 : auto poGroupQuality =
293 6 : poRootGroup->OpenGroup("QualityOfBathymetryCoverage");
294 2 : if (poGroupQuality)
295 : {
296 : auto poQualityOfBathymetryCoverage01 =
297 1 : poGroupQuality->OpenGroup(
298 3 : "QualityOfBathymetryCoverage.01");
299 1 : if (poQualityOfBathymetryCoverage01)
300 : {
301 1 : ++iSubDS;
302 : aosSubDSList.SetNameValue(
303 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
304 : CPLSPrintf(
305 : "S102:\"%s\":QualityOfBathymetryCoverage",
306 1 : osFilename.c_str()));
307 : aosSubDSList.SetNameValue(
308 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
309 : "Georeferenced metadata "
310 1 : "QualityOfBathymetryCoverage");
311 : }
312 : }
313 :
314 2 : poDS->GDALDataset::SetMetadata(aosSubDSList.List(),
315 : "SUBDATASETS");
316 :
317 : // Setup/check for pam .aux.xml.
318 2 : poDS->SetDescription(osFilename.c_str());
319 2 : poDS->TryLoadXML();
320 :
321 : // Setup overviews.
322 2 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
323 :
324 2 : return poDS.release();
325 : }
326 : }
327 : }
328 :
329 47 : if (bIsQuality)
330 : {
331 9 : if (!poDS->OpenQuality(poOpenInfo, poRootGroup))
332 2 : return nullptr;
333 :
334 : // Setup/check for pam .aux.xml.
335 7 : poDS->SetDescription(osFilename.c_str());
336 7 : poDS->TryLoadXML();
337 :
338 : // Setup overviews.
339 7 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
340 :
341 7 : return poDS.release();
342 : }
343 :
344 : auto poBathymetryCoverageInstance =
345 76 : poBathymetryCoverage->OpenGroup(osBathymetryCoverageName);
346 38 : if (!poBathymetryCoverageInstance)
347 : {
348 1 : CPLError(CE_Failure, CPLE_AppDefined,
349 : "S102: Cannot find %s group in BathymetryCoverage group",
350 : osBathymetryCoverageName.c_str());
351 1 : return nullptr;
352 : }
353 :
354 37 : if (auto poStartSequence =
355 74 : poBathymetryCoverage->GetAttribute("startSequence"))
356 : {
357 0 : const char *pszStartSequence = poStartSequence->ReadAsString();
358 0 : if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
359 : {
360 : // Shouldn't happen given this is imposed by the spec.
361 : // Cf 4.2.1.1.1.12 "startSequence" of Ed 3.0 spec, page 13
362 0 : CPLError(CE_Failure, CPLE_AppDefined,
363 : "startSequence (=%s) != 0,0 is not supported",
364 : pszStartSequence);
365 0 : return nullptr;
366 : }
367 : }
368 :
369 : // Potentially override vertical datum
370 37 : S100ReadVerticalDatum(poDS.get(), poBathymetryCoverageInstance.get());
371 :
372 37 : const bool bNorthUp = CPLTestBool(
373 37 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
374 :
375 : // Compute geotransform
376 37 : poDS->m_bHasGT = S100GetGeoTransform(poBathymetryCoverageInstance.get(),
377 37 : poDS->m_gt, bNorthUp);
378 :
379 111 : auto poGroup001 = poBathymetryCoverageInstance->OpenGroup("Group_001");
380 37 : if (!poGroup001)
381 : {
382 0 : CPLError(CE_Failure, CPLE_AppDefined,
383 : "S102: Cannot find "
384 : "/BathymetryCoverage/BathymetryCoverage.01/Group_001");
385 0 : return nullptr;
386 : }
387 111 : auto poValuesArray = poGroup001->OpenMDArray("values");
388 37 : if (!poValuesArray || poValuesArray->GetDimensionCount() != 2)
389 : {
390 0 : CPLError(CE_Failure, CPLE_AppDefined,
391 : "S102: Cannot find "
392 : "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values");
393 0 : return nullptr;
394 : }
395 37 : const auto &oType = poValuesArray->GetDataType();
396 37 : if (oType.GetClass() != GEDTC_COMPOUND)
397 : {
398 0 : CPLError(CE_Failure, CPLE_AppDefined,
399 : "S102: Wrong type for "
400 : "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values");
401 0 : return nullptr;
402 : }
403 37 : const auto &oComponents = oType.GetComponents();
404 37 : if (oComponents.size() == 0 || oComponents[0]->GetName() != "depth")
405 : {
406 0 : CPLError(CE_Failure, CPLE_AppDefined,
407 : "S102: Wrong type for "
408 : "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values");
409 0 : return nullptr;
410 : }
411 :
412 37 : if (bNorthUp)
413 36 : poValuesArray = poValuesArray->GetView("[::-1,...]");
414 :
415 111 : auto poDepth = poValuesArray->GetView("[\"depth\"]");
416 :
417 : // Mandatory in v2.2. Since v3.0, EPSG:6498 is the only allowed value
418 37 : bool bCSIsElevation = false;
419 111 : auto poVerticalCS = poRootGroup->GetAttribute("verticalCS");
420 37 : if (poVerticalCS && poVerticalCS->GetDataType().GetClass() == GEDTC_NUMERIC)
421 : {
422 27 : const auto nVal = poVerticalCS->ReadAsInt();
423 27 : if (nVal == 6498) // Depth metre
424 : {
425 : // nothing to do
426 : }
427 0 : else if (nVal == 6499) // Height metre
428 : {
429 0 : bCSIsElevation = true;
430 : }
431 : else
432 : {
433 0 : CPLError(CE_Warning, CPLE_NotSupported, "Unsupported verticalCS=%d",
434 : nVal);
435 : }
436 : }
437 :
438 : const bool bUseElevation =
439 37 : EQUAL(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
440 : "DEPTH_OR_ELEVATION", "DEPTH"),
441 : "ELEVATION");
442 73 : const bool bInvertDepth = (bUseElevation && !bCSIsElevation) ||
443 36 : (!bUseElevation && bCSIsElevation);
444 37 : const double dfDepthNoData = poDepth->GetNoDataValueAsDouble();
445 75 : auto poDepthDS = [&poDepth, bInvertDepth, dfDepthNoData]()
446 : {
447 37 : if (bInvertDepth)
448 : {
449 1 : auto poInverted = poDepth->GetUnscaled(-1, 0, dfDepthNoData);
450 : return std::unique_ptr<GDALDataset>(
451 1 : poInverted->AsClassicDataset(1, 0));
452 : }
453 : else
454 : {
455 : return std::unique_ptr<GDALDataset>(
456 36 : poDepth->AsClassicDataset(1, 0));
457 : }
458 74 : }();
459 :
460 37 : poDS->nRasterXSize = poDepthDS->GetRasterXSize();
461 37 : poDS->nRasterYSize = poDepthDS->GetRasterYSize();
462 :
463 : // Create depth (or elevation) band
464 37 : auto poDepthBand = new S102RasterBand(std::move(poDepthDS));
465 37 : poDepthBand->SetDescription(bUseElevation ? "elevation" : "depth");
466 :
467 111 : auto poMinimumDepth = poGroup001->GetAttribute("minimumDepth");
468 74 : if (poMinimumDepth &&
469 74 : poMinimumDepth->GetDataType().GetClass() == GEDTC_NUMERIC)
470 : {
471 37 : const double dfVal = poMinimumDepth->ReadAsDouble();
472 37 : if (dfVal != dfDepthNoData)
473 : {
474 33 : if (bInvertDepth)
475 1 : poDepthBand->m_dfMaximum = -dfVal;
476 : else
477 32 : poDepthBand->m_dfMinimum = dfVal;
478 : }
479 : }
480 :
481 111 : auto poMaximumDepth = poGroup001->GetAttribute("maximumDepth");
482 74 : if (poMaximumDepth &&
483 74 : poMaximumDepth->GetDataType().GetClass() == GEDTC_NUMERIC)
484 : {
485 37 : const double dfVal = poMaximumDepth->ReadAsDouble();
486 37 : if (dfVal != dfDepthNoData)
487 : {
488 37 : if (bInvertDepth)
489 1 : poDepthBand->m_dfMinimum = -dfVal;
490 : else
491 36 : poDepthBand->m_dfMaximum = dfVal;
492 : }
493 : }
494 :
495 37 : poDS->SetBand(1, poDepthBand);
496 :
497 : const bool bHasUncertainty =
498 37 : oComponents.size() >= 2 && oComponents[1]->GetName() == "uncertainty";
499 37 : if (bHasUncertainty)
500 : {
501 : // Create uncertainty band
502 54 : auto poUncertainty = poValuesArray->GetView("[\"uncertainty\"]");
503 : const double dfUncertaintyNoData =
504 18 : poUncertainty->GetNoDataValueAsDouble();
505 : auto poUncertaintyDS =
506 36 : std::unique_ptr<GDALDataset>(poUncertainty->AsClassicDataset(1, 0));
507 :
508 18 : auto poUncertaintyBand = new S102RasterBand(std::move(poUncertaintyDS));
509 18 : poUncertaintyBand->SetDescription("uncertainty");
510 :
511 : auto poMinimumUncertainty =
512 54 : poGroup001->GetAttribute("minimumUncertainty");
513 36 : if (poMinimumUncertainty &&
514 36 : poMinimumUncertainty->GetDataType().GetClass() == GEDTC_NUMERIC)
515 : {
516 18 : const double dfVal = poMinimumUncertainty->ReadAsDouble();
517 18 : if (dfVal != dfUncertaintyNoData)
518 : {
519 18 : poUncertaintyBand->m_dfMinimum = dfVal;
520 : }
521 : }
522 :
523 : auto poMaximumUncertainty =
524 54 : poGroup001->GetAttribute("maximumUncertainty");
525 36 : if (poMaximumUncertainty &&
526 36 : poMaximumUncertainty->GetDataType().GetClass() == GEDTC_NUMERIC)
527 : {
528 18 : const double dfVal = poMaximumUncertainty->ReadAsDouble();
529 18 : if (dfVal != dfUncertaintyNoData)
530 : {
531 18 : poUncertaintyBand->m_dfMaximum = dfVal;
532 : }
533 : }
534 :
535 18 : poDS->SetBand(2, poUncertaintyBand);
536 : }
537 :
538 37 : poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
539 :
540 111 : auto poGroupQuality = poRootGroup->OpenGroup("QualityOfSurvey");
541 37 : const bool bIsNamedQualityOfSurvey = poGroupQuality != nullptr;
542 37 : if (!bIsNamedQualityOfSurvey)
543 : {
544 : // S102 v3 now uses QualityOfBathymetryCoverage instead of QualityOfSurvey
545 35 : poGroupQuality = poRootGroup->OpenGroup("QualityOfBathymetryCoverage");
546 : }
547 37 : if (!bIsSubdataset && poGroupQuality)
548 : {
549 6 : const char *pszNameOfQualityGroup = bIsNamedQualityOfSurvey
550 6 : ? "QualityOfSurvey"
551 : : "QualityOfBathymetryCoverage";
552 6 : auto poGroupQuality01 = poGroupQuality->OpenGroup(
553 18 : CPLSPrintf("%s.01", pszNameOfQualityGroup));
554 6 : if (poGroupQuality01)
555 : {
556 6 : poDS->GDALDataset::SetMetadataItem(
557 : "SUBDATASET_1_NAME",
558 : CPLSPrintf("S102:\"%s\":BathymetryCoverage",
559 : osFilename.c_str()),
560 : "SUBDATASETS");
561 6 : poDS->GDALDataset::SetMetadataItem(
562 : "SUBDATASET_1_DESC", "Bathymetric gridded data", "SUBDATASETS");
563 :
564 6 : poDS->GDALDataset::SetMetadataItem(
565 : "SUBDATASET_2_NAME",
566 : CPLSPrintf("S102:\"%s\":%s", osFilename.c_str(),
567 : pszNameOfQualityGroup),
568 : "SUBDATASETS");
569 6 : poDS->GDALDataset::SetMetadataItem(
570 : "SUBDATASET_2_DESC",
571 : CPLSPrintf("Georeferenced metadata %s", pszNameOfQualityGroup),
572 : "SUBDATASETS");
573 : }
574 : }
575 :
576 : // Setup/check for pam .aux.xml.
577 37 : poDS->SetDescription(osFilename.c_str());
578 37 : poDS->TryLoadXML();
579 :
580 : // Setup overviews.
581 37 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
582 :
583 37 : return poDS.release();
584 : }
585 :
586 : /************************************************************************/
587 : /* OpenQuality() */
588 : /************************************************************************/
589 :
590 9 : bool S102Dataset::OpenQuality(GDALOpenInfo *poOpenInfo,
591 : const std::shared_ptr<GDALGroup> &poRootGroup)
592 : {
593 9 : const bool bNorthUp = CPLTestBool(
594 9 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
595 :
596 9 : const char *pszNameOfQualityGroup = "QualityOfSurvey";
597 27 : auto poGroupQuality = poRootGroup->OpenGroup(pszNameOfQualityGroup);
598 9 : if (!poGroupQuality)
599 : {
600 7 : pszNameOfQualityGroup = "QualityOfBathymetryCoverage";
601 7 : poGroupQuality = poRootGroup->OpenGroup(pszNameOfQualityGroup);
602 7 : if (!poGroupQuality)
603 : {
604 2 : CPLError(CE_Failure, CPLE_AppDefined,
605 : "Cannot find group /QualityOfSurvey or "
606 : "/QualityOfBathymetryCoverage");
607 2 : return false;
608 : }
609 : }
610 :
611 : const std::string osQuality01Name =
612 21 : std::string(pszNameOfQualityGroup).append(".01");
613 14 : const std::string osQuality01FullName = std::string("/")
614 7 : .append(pszNameOfQualityGroup)
615 7 : .append("/")
616 14 : .append(osQuality01Name);
617 14 : auto poGroupQuality01 = poGroupQuality->OpenGroup(osQuality01Name);
618 7 : if (!poGroupQuality01)
619 : {
620 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find group %s",
621 : osQuality01FullName.c_str());
622 0 : return false;
623 : }
624 :
625 14 : if (auto poStartSequence = poGroupQuality01->GetAttribute("startSequence"))
626 : {
627 3 : const char *pszStartSequence = poStartSequence->ReadAsString();
628 3 : if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
629 : {
630 0 : CPLError(CE_Failure, CPLE_AppDefined,
631 : "startSequence (=%s) != 0,0 is not supported",
632 : pszStartSequence);
633 0 : return false;
634 : }
635 : }
636 :
637 : // Compute geotransform
638 7 : m_bHasGT = S100GetGeoTransform(poGroupQuality01.get(), m_gt, bNorthUp);
639 :
640 21 : auto poGroup001 = poGroupQuality01->OpenGroup("Group_001");
641 7 : if (!poGroup001)
642 : {
643 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find group %s/Group_001",
644 : osQuality01FullName.c_str());
645 0 : return false;
646 : }
647 :
648 21 : auto poValuesArray = poGroup001->OpenMDArray("values");
649 7 : if (!poValuesArray)
650 : {
651 0 : CPLError(CE_Failure, CPLE_AppDefined,
652 : "Cannot find array "
653 : "%s/Group_001/values",
654 : osQuality01FullName.c_str());
655 0 : return false;
656 : }
657 :
658 : {
659 7 : const auto &oType = poValuesArray->GetDataType();
660 12 : if (oType.GetClass() == GEDTC_NUMERIC &&
661 5 : oType.GetNumericDataType() == GDT_UInt32)
662 : {
663 : // ok
664 : }
665 4 : else if (oType.GetClass() == GEDTC_COMPOUND &&
666 4 : oType.GetComponents().size() == 1 &&
667 2 : oType.GetComponents()[0]->GetType().GetClass() ==
668 4 : GEDTC_NUMERIC &&
669 2 : oType.GetComponents()[0]->GetType().GetNumericDataType() ==
670 : GDT_UInt32)
671 : {
672 : // seen in a S102 v3 product (102DE00CA22_UNC_MD.H5), although
673 : // I believe this is non-conformant.
674 :
675 : // Escape potentials single-quote and double-quote with back-slash
676 2 : CPLString osEscapedCompName(oType.GetComponents()[0]->GetName());
677 4 : osEscapedCompName.replaceAll("\\", "\\\\")
678 4 : .replaceAll("'", "\\'")
679 2 : .replaceAll("\"", "\\\"");
680 :
681 : // Gets a view with that single component extracted.
682 4 : poValuesArray = poValuesArray->GetView(
683 4 : std::string("['").append(osEscapedCompName).append("']"));
684 2 : if (!poValuesArray)
685 0 : return false;
686 : }
687 : else
688 : {
689 0 : CPLError(CE_Failure, CPLE_NotSupported,
690 : "Unsupported data type for %s",
691 0 : poValuesArray->GetFullName().c_str());
692 0 : return false;
693 : }
694 : }
695 :
696 7 : if (poValuesArray->GetDimensionCount() != 2)
697 : {
698 0 : CPLError(CE_Failure, CPLE_NotSupported,
699 : "Unsupported number of dimensions for %s",
700 0 : poValuesArray->GetFullName().c_str());
701 0 : return false;
702 : }
703 :
704 : auto poFeatureAttributeTable =
705 21 : poGroupQuality->OpenMDArray("featureAttributeTable");
706 7 : if (!poFeatureAttributeTable)
707 : {
708 0 : CPLError(CE_Failure, CPLE_AppDefined,
709 : "Cannot find array /%s/featureAttributeTable",
710 : pszNameOfQualityGroup);
711 0 : return false;
712 : }
713 :
714 : {
715 7 : const auto &oType = poFeatureAttributeTable->GetDataType();
716 7 : if (oType.GetClass() != GEDTC_COMPOUND)
717 : {
718 0 : CPLError(CE_Failure, CPLE_NotSupported,
719 : "Unsupported data type for %s",
720 0 : poFeatureAttributeTable->GetFullName().c_str());
721 0 : return false;
722 : }
723 :
724 7 : const auto &poComponents = oType.GetComponents();
725 7 : if (poComponents.size() >= 1 && poComponents[0]->GetName() != "id")
726 : {
727 0 : CPLError(CE_Failure, CPLE_AppDefined,
728 : "Missing 'id' component in %s",
729 0 : poFeatureAttributeTable->GetFullName().c_str());
730 0 : return false;
731 : }
732 : }
733 :
734 7 : if (bNorthUp)
735 5 : poValuesArray = poValuesArray->GetView("[::-1,...]");
736 :
737 : auto poDS =
738 14 : std::unique_ptr<GDALDataset>(poValuesArray->AsClassicDataset(1, 0));
739 7 : if (!poDS)
740 0 : return false;
741 :
742 7 : nRasterXSize = poDS->GetRasterXSize();
743 7 : nRasterYSize = poDS->GetRasterYSize();
744 :
745 : auto poRAT =
746 14 : HDF5CreateRAT(poFeatureAttributeTable, /* bFirstColIsMinMax = */ true);
747 : auto poBand = std::make_unique<S102GeoreferencedMetadataRasterBand>(
748 7 : std::move(poDS), std::move(poRAT));
749 7 : SetBand(1, poBand.release());
750 :
751 7 : return true;
752 : }
753 :
754 : /************************************************************************/
755 : /* S102Creator */
756 : /************************************************************************/
757 :
758 : class S102Creator final : public S100BaseWriter
759 : {
760 : public:
761 56 : S102Creator(const char *pszDestFilename, GDALDataset *poSrcDS,
762 : CSLConstList papszOptions)
763 56 : : S100BaseWriter(pszDestFilename, poSrcDS, papszOptions)
764 : {
765 56 : }
766 :
767 : ~S102Creator() override;
768 :
769 : bool Create(GDALProgressFunc pfnProgress, void *pProgressData);
770 :
771 : // From the S102 spec
772 : static constexpr float NODATA = 1000000.0f;
773 : static constexpr const char *FEATURE_TYPE = "BathymetryCoverage";
774 : static constexpr const char *QUALITY_FEATURE_TYPE =
775 : "QualityOfBathymetryCoverage";
776 :
777 : protected:
778 73 : bool Close() override
779 : {
780 73 : return BaseClose();
781 : }
782 :
783 : private:
784 : bool WriteFeatureGroupAttributes(bool isQuality);
785 : bool CopyValues(GDALProgressFunc pfnProgress, void *pProgressData);
786 : bool CopyQualityValues(GDALDataset *poQualityDS,
787 : const std::set<int> &oSetRATId,
788 : GDALProgressFunc pfnProgress, void *pProgressData);
789 : bool WriteFeatureAttributeTable(const GDALRasterAttributeTable *poRAT);
790 : bool CreateGroupF(bool hasQualityOfBathymetryCoverage);
791 : };
792 :
793 : /************************************************************************/
794 : /* S102Creator::~S102Creator() */
795 : /************************************************************************/
796 :
797 56 : S102Creator::~S102Creator()
798 : {
799 56 : S102Creator::Close();
800 56 : }
801 :
802 : /************************************************************************/
803 : /* S102Creator::Create() */
804 : /************************************************************************/
805 :
806 : // S102 v3.0 Table 10-8 - Elements of featureAttributeTable compound datatype
807 : static const struct
808 : {
809 : const char *pszName;
810 : const char *pszType;
811 : } gasFeatureAttributeTableMembers[] = {
812 : {"id", "uint32"},
813 : {"dataAssessment", "uint8"},
814 : {"featuresDetected.leastDepthOfDetectedFeaturesMeasured", "boolean"},
815 : {"featuresDetected.significantFeaturesDetected", "boolean"},
816 : {"featuresDetected.sizeOfFeaturesDetected", "float32"},
817 : {"featureSizeVar", "float32"},
818 : {"fullSeafloorCoverageAchieved", "boolean"},
819 : {"bathyCoverage", "boolean"},
820 : {"zoneOfConfidence.horizontalPositionUncertainty.uncertaintyFixed",
821 : "float32"},
822 : {"zoneOfConfidence.horizontalPositionUncertainty.uncertaintyVariableFactor",
823 : "float32"},
824 : {"surveyDateRange.dateStart", "date"},
825 : {"surveyDateRange.dateEnd", "date"},
826 : {"sourceSurveyID", "string"},
827 : {"surveyAuthority", "string"},
828 : {"typeOfBathymetricEstimationUncertainty", "enumeration"},
829 : };
830 :
831 56 : bool S102Creator::Create(GDALProgressFunc pfnProgress, void *pProgressData)
832 : {
833 56 : if (m_poSrcDS->GetRasterCount() != 1 && m_poSrcDS->GetRasterCount() != 2)
834 : {
835 4 : CPLError(CE_Failure, CPLE_NotSupported,
836 : "Source dataset must have one or two bands");
837 4 : return false;
838 : }
839 :
840 52 : if (!BaseChecks("S102", /* crsMustBeEPSG = */ true,
841 : /* verticalDatumRequired = */ true))
842 21 : return false;
843 :
844 : const bool bAppendSubdataset =
845 31 : CPLTestBool(m_aosOptions.FetchNameValueDef("APPEND_SUBDATASET", "NO"));
846 :
847 31 : std::unique_ptr<GDALDataset> poQualityDS;
848 : const char *pszQualityDataset =
849 31 : m_aosOptions.FetchNameValue("QUALITY_DATASET");
850 31 : const GDALRasterAttributeTable *poRAT = nullptr;
851 31 : if (!pszQualityDataset && !bAppendSubdataset)
852 : {
853 : const char *pszSubDSName =
854 13 : m_poSrcDS->GetMetadataItem("SUBDATASET_2_NAME", "SUBDATASETS");
855 13 : if (pszSubDSName &&
856 13 : cpl::starts_with(std::string_view(pszSubDSName), "S102:") &&
857 13 : cpl::ends_with(std::string_view(pszSubDSName),
858 : ":QualityOfBathymetryCoverage"))
859 : {
860 0 : pszQualityDataset = pszSubDSName;
861 : }
862 : }
863 :
864 62 : std::set<int> oSetRATId;
865 31 : if (pszQualityDataset)
866 : {
867 15 : if (bAppendSubdataset)
868 : {
869 0 : CPLError(CE_Failure, CPLE_NotSupported,
870 : "Quality dataset can only be set on initial creation");
871 12 : return false;
872 : }
873 15 : poQualityDS.reset(GDALDataset::Open(
874 : pszQualityDataset, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
875 : nullptr, nullptr));
876 15 : if (!poQualityDS)
877 1 : return false;
878 :
879 14 : if (poQualityDS->GetRasterCount() != 1)
880 : {
881 1 : CPLError(CE_Failure, CPLE_AppDefined,
882 : "%s does not have a single band.", pszQualityDataset);
883 1 : return false;
884 : }
885 13 : if (!GDALDataTypeIsInteger(
886 : poQualityDS->GetRasterBand(1)->GetRasterDataType()))
887 : {
888 1 : CPLError(CE_Failure, CPLE_AppDefined,
889 : "%s band is not of an integer data type.",
890 : pszQualityDataset);
891 1 : return false;
892 : }
893 23 : if (poQualityDS->GetRasterXSize() != m_poSrcDS->GetRasterXSize() ||
894 11 : poQualityDS->GetRasterYSize() != m_poSrcDS->GetRasterYSize())
895 : {
896 2 : CPLError(CE_Failure, CPLE_AppDefined,
897 : "%s does not have the same dimensions as %s.",
898 2 : pszQualityDataset, m_poSrcDS->GetDescription());
899 2 : return false;
900 : }
901 :
902 10 : const auto poQualityDS_SRS = poQualityDS->GetSpatialRef();
903 10 : if (!poQualityDS_SRS || !poQualityDS_SRS->IsSame(m_poSRS))
904 : {
905 1 : CPLError(CE_Failure, CPLE_AppDefined,
906 : "%s does not have the same CRS as %s.", pszQualityDataset,
907 1 : m_poSrcDS->GetDescription());
908 1 : return false;
909 : }
910 :
911 9 : GDALGeoTransform gt;
912 9 : if (poQualityDS->GetGeoTransform(gt) != CE_None || gt != m_gt)
913 : {
914 1 : CPLError(CE_Failure, CPLE_AppDefined,
915 : "%s does not have the same geotransform as %s.",
916 1 : pszQualityDataset, m_poSrcDS->GetDescription());
917 1 : return false;
918 : }
919 :
920 8 : poRAT = poQualityDS->GetRasterBand(1)->GetDefaultRAT();
921 8 : if (!poRAT)
922 : {
923 1 : CPLError(CE_Failure, CPLE_AppDefined,
924 : "%s does not have a raster attribute table.",
925 1 : poQualityDS->GetDescription());
926 1 : return false;
927 : }
928 :
929 7 : const int nRATColumnCount = poRAT->GetColumnCount();
930 7 : std::set<std::string_view> setKnownColumnNames;
931 112 : for (const auto &entry : gasFeatureAttributeTableMembers)
932 105 : setKnownColumnNames.insert(entry.pszName);
933 7 : int iRATIdField = -1;
934 24 : for (int i = 0; i < nRATColumnCount; ++i)
935 : {
936 17 : const char *pszColName = poRAT->GetNameOfCol(i);
937 17 : if (strcmp(pszColName, "id") == 0)
938 : {
939 6 : iRATIdField = i;
940 : }
941 11 : else if (!cpl::contains(setKnownColumnNames, pszColName))
942 : {
943 5 : CPLError(CE_Warning, CPLE_AppDefined,
944 : "'%s' is not a valid S102 feature attribute table "
945 : "column name.",
946 : pszColName);
947 : }
948 : }
949 7 : if (iRATIdField < 0)
950 : {
951 1 : CPLError(
952 : CE_Failure, CPLE_AppDefined,
953 : "Input raster attribute table lacks an integer 'id' field");
954 1 : return false;
955 : }
956 6 : const int nRATRowCount = poRAT->GetRowCount();
957 11 : for (int i = 0; i < nRATRowCount; ++i)
958 : {
959 8 : const int nID = poRAT->GetValueAsInt(i, iRATIdField);
960 8 : if (nID == 0)
961 : {
962 1 : CPLError(CE_Failure, CPLE_AppDefined,
963 : "id=0 is not allowed in input raster attribute table");
964 3 : return false;
965 : }
966 7 : else if (nID < 0)
967 : {
968 1 : CPLError(CE_Failure, CPLE_AppDefined,
969 : "Negative id is not allowed in input raster attribute "
970 : "table");
971 1 : return false;
972 : }
973 6 : else if (cpl::contains(oSetRATId, nID))
974 : {
975 1 : CPLError(
976 : CE_Failure, CPLE_AppDefined,
977 : "Several rows of input raster attribute table have id=%d",
978 : nID);
979 1 : return false;
980 : }
981 5 : oSetRATId.insert(nID);
982 : }
983 : }
984 :
985 19 : if (!((m_nVerticalDatum >= 1 && m_nVerticalDatum <= 30) ||
986 0 : m_nVerticalDatum == 44))
987 : {
988 0 : CPLError(CE_Warning, CPLE_AppDefined,
989 : "VERTICAL_DATUM=%d value is a valid S100 value but not "
990 : "allowed in S102. Valid values are [1, 30] or 44",
991 : m_nVerticalDatum);
992 : }
993 :
994 19 : if (!(m_nEPSGCode == 4326 || m_nEPSGCode == 5041 || m_nEPSGCode == 5042 ||
995 19 : (m_nEPSGCode >= 32601 && m_nEPSGCode <= 32660) ||
996 1 : (m_nEPSGCode >= 32701 && m_nEPSGCode <= 32760)))
997 : {
998 1 : CPLError(CE_Warning, CPLE_NotSupported,
999 : "The EPSG code of the CRS is %d. "
1000 : "Only EPSG codes 4326, 5041, 5042, [32601, 32660], "
1001 : "[32701, 32760] are officially supported. "
1002 : "The dataset may not be recognized by other software",
1003 : m_nEPSGCode);
1004 : }
1005 :
1006 19 : if (bAppendSubdataset)
1007 : {
1008 6 : GDALOpenInfo oOpenInfo(m_osDestFilename.c_str(), GA_ReadOnly);
1009 : auto poOriDS =
1010 6 : std::unique_ptr<GDALDataset>(S102Dataset::Open(&oOpenInfo));
1011 3 : if (!poOriDS)
1012 : {
1013 1 : CPLError(CE_Failure, CPLE_AppDefined,
1014 : "%s is not a valid existing S102 dataset",
1015 : m_osDestFilename.c_str());
1016 1 : return false;
1017 : }
1018 2 : const auto poOriSRS = poOriDS->GetSpatialRef();
1019 2 : if (!poOriSRS)
1020 : {
1021 : // shouldn't happen
1022 0 : return false;
1023 : }
1024 2 : if (!poOriSRS->IsSame(m_poSRS))
1025 : {
1026 1 : CPLError(CE_Failure, CPLE_AppDefined,
1027 : "CRS of %s is not the same as the one of %s",
1028 1 : m_osDestFilename.c_str(), m_poSrcDS->GetDescription());
1029 1 : return false;
1030 : }
1031 1 : poOriDS.reset();
1032 :
1033 1 : OGREnvelope sExtent;
1034 1 : if (m_poSrcDS->GetExtentWGS84LongLat(&sExtent) != OGRERR_NONE)
1035 : {
1036 0 : CPLError(CE_Failure, CPLE_AppDefined,
1037 : "Cannot get dataset extent in WGS84 longitude/latitude");
1038 0 : return false;
1039 : }
1040 :
1041 1 : bool ret = OpenFileUpdateMode();
1042 1 : if (ret)
1043 : {
1044 1 : m_featureGroup.reset(
1045 : H5_CHECK(H5Gopen(m_hdf5, "BathymetryCoverage")));
1046 : }
1047 :
1048 1 : ret = ret && m_featureGroup;
1049 1 : double dfNumInstances = 0;
1050 1 : ret = ret && GH5_FetchAttribute(m_featureGroup, "numInstances",
1051 : dfNumInstances, true);
1052 1 : if (ret && !(dfNumInstances >= 1 && dfNumInstances <= 99 &&
1053 1 : std::round(dfNumInstances) == dfNumInstances))
1054 : {
1055 0 : CPLError(CE_Failure, CPLE_AppDefined,
1056 : "Invalid value for numInstances");
1057 0 : ret = false;
1058 : }
1059 1 : else if (ret && dfNumInstances == 99)
1060 : {
1061 0 : CPLError(CE_Failure, CPLE_AppDefined,
1062 : "Too many existing feature instances");
1063 0 : ret = false;
1064 : }
1065 : else
1066 : {
1067 1 : double dfMainVerticalDatum = 0;
1068 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "verticalDatum",
1069 : dfMainVerticalDatum, true);
1070 :
1071 1 : const int newNumInstances = static_cast<int>(dfNumInstances) + 1;
1072 1 : ret = ret && GH5_WriteAttribute(m_featureGroup, "numInstances",
1073 : newNumInstances);
1074 1 : ret = ret && CreateFeatureInstanceGroup(CPLSPrintf(
1075 : "BathymetryCoverage.%02d", newNumInstances));
1076 1 : ret = ret && WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1077 1 : if (dfMainVerticalDatum != m_nVerticalDatum)
1078 : {
1079 2 : ret = ret &&
1080 1 : GH5_CreateAttribute(m_featureInstanceGroup,
1081 : "verticalDatumReference",
1082 2 : H5T_STD_U8LE) &&
1083 : // s100VerticalDatum
1084 1 : GH5_WriteAttribute(m_featureInstanceGroup,
1085 : "verticalDatumReference", 1);
1086 1 : ret =
1087 2 : ret && WriteVerticalDatum(m_featureInstanceGroup,
1088 1 : H5T_STD_U16LE, m_nVerticalDatum);
1089 : }
1090 :
1091 1 : ret = ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U8LE, 1);
1092 :
1093 1 : ret = ret && CreateValuesGroup("Group_001");
1094 1 : ret = ret && WriteVarLengthStringValue(m_valuesGroup, "timePoint",
1095 : "00010101T000000Z");
1096 1 : ret = ret && CopyValues(pfnProgress, pProgressData);
1097 : }
1098 :
1099 : // Update global bounding box
1100 1 : OGREnvelope sExistingExtent;
1101 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "westBoundLongitude",
1102 : sExistingExtent.MinX, true);
1103 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "southBoundLatitude",
1104 : sExistingExtent.MinY, true);
1105 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "eastBoundLongitude",
1106 : sExistingExtent.MaxX, true);
1107 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "northBoundLatitude",
1108 : sExistingExtent.MaxY, true);
1109 :
1110 1 : sExtent.Merge(sExistingExtent);
1111 2 : ret = ret &&
1112 1 : GH5_WriteAttribute(m_hdf5, "westBoundLongitude", sExtent.MinX);
1113 2 : ret = ret &&
1114 1 : GH5_WriteAttribute(m_hdf5, "southBoundLatitude", sExtent.MinY);
1115 2 : ret = ret &&
1116 1 : GH5_WriteAttribute(m_hdf5, "eastBoundLongitude", sExtent.MaxX);
1117 2 : ret = ret &&
1118 1 : GH5_WriteAttribute(m_hdf5, "northBoundLatitude", sExtent.MaxY);
1119 :
1120 1 : return Close() && ret;
1121 : }
1122 : else
1123 : {
1124 16 : bool ret = CreateFile();
1125 16 : ret = ret && WriteProductSpecification("INT.IHO.S-102.3.0.0");
1126 16 : ret = ret && WriteIssueDate();
1127 16 : ret = ret && WriteIssueTime(/* bAutogenerateFromCurrent = */ false);
1128 16 : ret = ret && WriteHorizontalCRS();
1129 16 : ret = ret && WriteTopLevelBoundingBox();
1130 16 : ret = ret && WriteVerticalCS(6498); // Depth, metre, down
1131 16 : ret = ret && WriteVerticalCoordinateBase(2); // verticalDatum
1132 : // s100VerticalDatum
1133 16 : ret = ret && WriteVerticalDatumReference(m_hdf5, 1);
1134 16 : ret =
1135 16 : ret && WriteVerticalDatum(m_hdf5, H5T_STD_U16LE, m_nVerticalDatum);
1136 :
1137 : // BathymetryCoverage
1138 16 : ret = ret && CreateFeatureGroup(FEATURE_TYPE);
1139 16 : ret = ret && WriteFeatureGroupAttributes(/* isQuality = */ false);
1140 16 : ret = ret && WriteAxisNames(m_featureGroup);
1141 :
1142 16 : ret = ret && CreateFeatureInstanceGroup("BathymetryCoverage.01");
1143 16 : ret = ret && WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1144 16 : ret = ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U8LE, 1);
1145 :
1146 16 : ret = ret && CreateValuesGroup("Group_001");
1147 :
1148 16 : ret = ret && WriteVarLengthStringValue(m_valuesGroup, "timePoint",
1149 : "00010101T000000Z");
1150 :
1151 : const double dfIntermediatePct =
1152 16 : m_poSrcDS->GetRasterCount() /
1153 16 : (m_poSrcDS->GetRasterCount() + (poQualityDS ? 1.0 : 0.0));
1154 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>
1155 : pScaledProgressData(GDALCreateScaledProgress(0.0, dfIntermediatePct,
1156 : pfnProgress,
1157 : pProgressData),
1158 16 : GDALDestroyScaledProgress);
1159 16 : ret = ret && CopyValues(GDALScaledProgress, pScaledProgressData.get());
1160 :
1161 16 : if (poQualityDS)
1162 : {
1163 : // QualityOfBathymetryCoverage group
1164 3 : ret = ret && CreateFeatureGroup(QUALITY_FEATURE_TYPE);
1165 3 : ret = ret && WriteFeatureGroupAttributes(/* isQuality = */ true);
1166 3 : ret = ret && WriteAxisNames(m_featureGroup);
1167 3 : ret = ret && WriteFeatureAttributeTable(poRAT);
1168 :
1169 6 : ret = ret &&
1170 3 : CreateFeatureInstanceGroup("QualityOfBathymetryCoverage.01");
1171 3 : ret = ret && WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1172 3 : ret = ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U8LE, 1);
1173 :
1174 3 : ret = ret && CreateValuesGroup("Group_001");
1175 3 : pScaledProgressData.reset(GDALCreateScaledProgress(
1176 : dfIntermediatePct, 1.0, pfnProgress, pProgressData));
1177 3 : ret = ret && CopyQualityValues(poQualityDS.get(), oSetRATId,
1178 : GDALScaledProgress,
1179 : pScaledProgressData.get());
1180 : }
1181 :
1182 16 : ret = ret && CreateGroupF(poQualityDS != nullptr);
1183 :
1184 16 : return Close() && ret;
1185 : }
1186 : }
1187 :
1188 : /************************************************************************/
1189 : /* S102Creator::WriteFeatureGroupAttributes() */
1190 : /************************************************************************/
1191 :
1192 18 : bool S102Creator::WriteFeatureGroupAttributes(bool isQuality)
1193 : {
1194 18 : CPLAssert(m_featureGroup);
1195 :
1196 18 : bool ret = WriteCommonPointRule(m_featureGroup, 2); // low
1197 18 : if (isQuality)
1198 : {
1199 : // Feature oriented Regular Grid
1200 3 : ret = ret && WriteDataCodingFormat(m_featureGroup, 9);
1201 : }
1202 : else
1203 : {
1204 15 : ret = ret && WriteDataCodingFormat(m_featureGroup, 2); // Regular grid
1205 : }
1206 18 : ret = ret && WriteDataOffsetCode(m_featureGroup, 5); // Center of cell
1207 18 : ret = ret && WriteDimension(m_featureGroup, 2);
1208 : const char *pszHorizontalPositionUncertainty =
1209 18 : m_aosOptions.FetchNameValue("HORIZONTAL_POSITION_UNCERTAINTY");
1210 18 : ret =
1211 36 : ret &&
1212 18 : WriteHorizontalPositionUncertainty(
1213 : m_featureGroup,
1214 0 : pszHorizontalPositionUncertainty &&
1215 0 : pszHorizontalPositionUncertainty[0]
1216 0 : ? static_cast<float>(CPLAtof(pszHorizontalPositionUncertainty))
1217 : : -1.0f);
1218 : const char *pszVerticalUncertainty =
1219 18 : m_aosOptions.FetchNameValue("VERTICAL_UNCERTAINTY");
1220 18 : ret = ret && WriteVerticalUncertainty(
1221 : m_featureGroup,
1222 0 : pszVerticalUncertainty && pszVerticalUncertainty[0]
1223 0 : ? static_cast<float>(CPLAtof(pszVerticalUncertainty))
1224 : : -1.0f);
1225 18 : ret = ret && WriteInterpolationType(m_featureGroup, 1); // Nearest neighbor
1226 18 : ret = ret && WriteNumInstances(m_featureGroup, H5T_STD_U8LE, 1);
1227 36 : ret = ret && WriteSequencingRuleScanDirection(m_featureGroup,
1228 18 : m_poSRS->IsProjected()
1229 : ? "Easting, Northing"
1230 : : "Longitude, Latitude");
1231 18 : ret = ret && WriteSequencingRuleType(m_featureGroup, 1); // Linear
1232 18 : return ret;
1233 : }
1234 :
1235 : /************************************************************************/
1236 : /* S102Creator::WriteFeatureAttributeTable() */
1237 : /************************************************************************/
1238 :
1239 3 : bool S102Creator::WriteFeatureAttributeTable(
1240 : const GDALRasterAttributeTable *poRAT)
1241 : {
1242 3 : CPLAssert(m_featureGroup);
1243 :
1244 6 : std::map<std::string_view, const char *> mapKnownColumns;
1245 48 : for (const auto &entry : gasFeatureAttributeTableMembers)
1246 45 : mapKnownColumns[entry.pszName] = entry.pszType;
1247 :
1248 3 : const int nColCount = poRAT->GetColumnCount();
1249 :
1250 3 : size_t nCompoundSize = 0;
1251 3 : size_t nMEMCompoundSize = 0;
1252 16 : for (int i = 0; i < nColCount; ++i)
1253 : {
1254 13 : const char *pszColName = poRAT->GetNameOfCol(i);
1255 13 : const auto iter = mapKnownColumns.find(pszColName);
1256 13 : size_t nMemberSize = sizeof(char *);
1257 13 : if (iter != mapKnownColumns.end())
1258 : {
1259 9 : const char *pszType = iter->second;
1260 9 : if (strcmp(pszType, "uint8") == 0 ||
1261 8 : strcmp(pszType, "boolean") == 0 ||
1262 7 : strcmp(pszType, "enumeration") == 0)
1263 : {
1264 3 : nMemberSize = sizeof(uint8_t);
1265 : }
1266 6 : else if (strcmp(pszType, "uint32") == 0)
1267 : {
1268 3 : nMemberSize = sizeof(uint32_t);
1269 : }
1270 3 : else if (strcmp(pszType, "float32") == 0)
1271 : {
1272 1 : nMemberSize = sizeof(float);
1273 : }
1274 2 : else if (strcmp(pszType, "string") == 0 ||
1275 1 : strcmp(pszType, "date") == 0)
1276 : {
1277 2 : nMemberSize = sizeof(char *);
1278 : }
1279 : else
1280 : {
1281 0 : CPLAssert(false);
1282 : }
1283 : }
1284 : else
1285 : {
1286 4 : GDALRATFieldType eType = poRAT->GetTypeOfCol(i);
1287 4 : switch (eType)
1288 : {
1289 1 : case GFT_Integer:
1290 1 : nMemberSize = sizeof(int32_t);
1291 1 : break;
1292 1 : case GFT_Real:
1293 1 : nMemberSize = sizeof(double);
1294 1 : break;
1295 1 : case GFT_Boolean:
1296 1 : nMemberSize = sizeof(uint8_t);
1297 1 : break;
1298 1 : case GFT_String:
1299 : case GFT_DateTime:
1300 : case GFT_WKBGeometry:
1301 1 : nMemberSize = sizeof(char *);
1302 1 : break;
1303 : }
1304 : }
1305 13 : nCompoundSize += nMemberSize;
1306 13 : if ((nMEMCompoundSize % nMemberSize) != 0)
1307 2 : nMEMCompoundSize += nMemberSize - (nMEMCompoundSize % nMemberSize);
1308 13 : nMEMCompoundSize += nMemberSize;
1309 : }
1310 :
1311 : GH5_HIDTypeHolder hDataType(
1312 6 : H5_CHECK(H5Tcreate(H5T_COMPOUND, nCompoundSize)));
1313 : GH5_HIDTypeHolder hDataTypeMEM(
1314 6 : H5_CHECK(H5Tcreate(H5T_COMPOUND, nMEMCompoundSize)));
1315 6 : GH5_HIDTypeHolder hVarLengthType(H5_CHECK(H5Tcopy(H5T_C_S1)));
1316 6 : bool bRet = hDataType && hDataTypeMEM && hVarLengthType &&
1317 9 : H5_CHECK(H5Tset_size(hVarLengthType, H5T_VARIABLE)) >= 0 &&
1318 3 : H5_CHECK(H5Tset_strpad(hVarLengthType, H5T_STR_NULLTERM)) >= 0;
1319 :
1320 6 : GH5_HIDTypeHolder hEnumType;
1321 6 : std::vector<const char *> apszTypes;
1322 :
1323 3 : size_t nOffset = 0;
1324 3 : size_t nMEMOffset = 0;
1325 6 : std::vector<size_t> anMEMOffsets;
1326 16 : for (int i = 0; i < nColCount && bRet; ++i)
1327 : {
1328 13 : const char *pszColName = poRAT->GetNameOfCol(i);
1329 13 : const auto iter = mapKnownColumns.find(pszColName);
1330 13 : hid_t hMemberType = hVarLengthType.get();
1331 13 : hid_t hMemberNativeType = hVarLengthType.get();
1332 13 : if (iter != mapKnownColumns.end())
1333 : {
1334 9 : const char *pszType = iter->second;
1335 9 : if (strcmp(pszType, "uint8") == 0 ||
1336 8 : strcmp(pszType, "boolean") == 0)
1337 : {
1338 2 : hMemberType = H5T_STD_U8LE;
1339 2 : hMemberNativeType = H5T_NATIVE_UCHAR;
1340 : }
1341 7 : else if (strcmp(pszType, "uint32") == 0)
1342 : {
1343 3 : hMemberType = H5T_STD_U32LE;
1344 3 : hMemberNativeType = H5T_NATIVE_UINT;
1345 : }
1346 4 : else if (strcmp(pszType, "float32") == 0)
1347 : {
1348 1 : hMemberType = H5T_IEEE_F32LE;
1349 1 : hMemberNativeType = H5T_NATIVE_FLOAT;
1350 : }
1351 3 : else if (strcmp(pszType, "string") == 0 ||
1352 2 : strcmp(pszType, "date") == 0)
1353 : {
1354 2 : hMemberType = hVarLengthType.get();
1355 2 : hMemberNativeType = hVarLengthType.get();
1356 : }
1357 1 : else if (strcmp(pszType, "enumeration") == 0 &&
1358 1 : strcmp(pszColName,
1359 : "typeOfBathymetricEstimationUncertainty") == 0)
1360 : {
1361 1 : hEnumType.reset(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1362 1 : bRet = hEnumType;
1363 1 : if (bRet)
1364 : {
1365 : uint8_t val;
1366 1 : val = 1;
1367 2 : bRet = bRet &&
1368 1 : H5_CHECK(H5Tenum_insert(
1369 : hEnumType, "rawStandardDeviation", &val)) >= 0;
1370 1 : val = 2;
1371 2 : bRet = bRet &&
1372 1 : H5_CHECK(H5Tenum_insert(
1373 : hEnumType, "cUBEStandardDeviation", &val)) >= 0;
1374 1 : val = 3;
1375 2 : bRet = bRet &&
1376 1 : H5_CHECK(H5Tenum_insert(
1377 : hEnumType, "productUncertainty", &val)) >= 0;
1378 1 : val = 4;
1379 1 : bRet = bRet && H5_CHECK(H5Tenum_insert(
1380 : hEnumType, "historicalStandardDeviation",
1381 : &val)) >= 0;
1382 :
1383 1 : hMemberType = hEnumType.get();
1384 1 : hMemberNativeType = hEnumType.get();
1385 1 : }
1386 : }
1387 : else
1388 : {
1389 0 : CPLAssert(false);
1390 : }
1391 9 : apszTypes.push_back(pszType);
1392 : }
1393 : else
1394 : {
1395 4 : GDALRATFieldType eType = poRAT->GetTypeOfCol(i);
1396 4 : switch (eType)
1397 : {
1398 1 : case GFT_Integer:
1399 1 : hMemberType = H5T_STD_I32LE;
1400 1 : hMemberNativeType = H5T_NATIVE_INT;
1401 1 : apszTypes.push_back("int32");
1402 1 : break;
1403 1 : case GFT_Real:
1404 1 : hMemberType = H5T_IEEE_F64LE;
1405 1 : hMemberNativeType = H5T_NATIVE_DOUBLE;
1406 1 : apszTypes.push_back("float64");
1407 1 : break;
1408 1 : case GFT_Boolean:
1409 1 : hMemberType = H5T_STD_U8LE;
1410 1 : hMemberNativeType = H5T_NATIVE_UCHAR;
1411 1 : apszTypes.push_back("boolean");
1412 1 : break;
1413 1 : case GFT_String:
1414 : case GFT_DateTime:
1415 : case GFT_WKBGeometry:
1416 1 : apszTypes.push_back("string");
1417 1 : break;
1418 : }
1419 : }
1420 :
1421 13 : CPLAssert(H5Tget_size(hMemberType) == H5Tget_size(hMemberNativeType));
1422 :
1423 13 : bRet = bRet && H5_CHECK(H5Tinsert(hDataType, pszColName, nOffset,
1424 : hMemberType)) >= 0;
1425 :
1426 13 : const size_t nMemberSize = H5Tget_size(hMemberType);
1427 13 : if ((nMEMOffset % nMemberSize) != 0)
1428 2 : nMEMOffset += nMemberSize - (nMEMOffset % nMemberSize);
1429 13 : anMEMOffsets.push_back(nMEMOffset);
1430 13 : bRet = bRet && H5_CHECK(H5Tinsert(hDataTypeMEM, pszColName, nMEMOffset,
1431 : hMemberNativeType)) >= 0;
1432 13 : nOffset += nMemberSize;
1433 13 : nMEMOffset += nMemberSize;
1434 : }
1435 3 : CPLAssert(nOffset == nCompoundSize);
1436 3 : CPLAssert(nMEMOffset == nMEMCompoundSize);
1437 :
1438 3 : CPLAssert(apszTypes.size() == static_cast<size_t>(nColCount));
1439 :
1440 3 : const int nRowCount = poRAT->GetRowCount();
1441 3 : hsize_t dims[] = {static_cast<hsize_t>(nRowCount)};
1442 6 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(1, dims, nullptr)));
1443 3 : bRet = bRet && hDataSpace;
1444 6 : GH5_HIDDatasetHolder hDatasetID;
1445 6 : GH5_HIDSpaceHolder hFileSpace;
1446 6 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1447 3 : bRet = bRet && hParams;
1448 3 : if (bRet)
1449 : {
1450 3 : bRet = H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0;
1451 3 : hsize_t chunk_size[] = {static_cast<hsize_t>(1)};
1452 3 : bRet = bRet && H5_CHECK(H5Pset_chunk(hParams, 1, chunk_size)) >= 0;
1453 3 : hDatasetID.reset(
1454 : H5_CHECK(H5Dcreate(m_featureGroup, "featureAttributeTable",
1455 : hDataType, hDataSpace, hParams)));
1456 3 : bRet = bRet && hDatasetID;
1457 : }
1458 3 : if (bRet)
1459 : {
1460 3 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
1461 3 : bRet = hFileSpace;
1462 : }
1463 :
1464 3 : hsize_t count[] = {1};
1465 6 : GH5_HIDSpaceHolder hMemSpace(H5_CHECK(H5Screate_simple(1, count, nullptr)));
1466 3 : bRet = bRet && hMemSpace;
1467 :
1468 6 : std::vector<GByte> abyBuffer(nMEMCompoundSize);
1469 3 : std::vector<std::string> asBuffers(nColCount);
1470 7 : for (int iRow = 0; iRow < nRowCount && bRet; ++iRow)
1471 : {
1472 24 : for (int iCol = 0; iCol < nColCount && bRet; ++iCol)
1473 : {
1474 20 : const char *const pszType = apszTypes[iCol];
1475 20 : GByte *const pabyDst = abyBuffer.data() + anMEMOffsets[iCol];
1476 20 : if (strcmp(pszType, "uint8") == 0 ||
1477 18 : strcmp(pszType, "boolean") == 0 ||
1478 15 : strcmp(pszType, "enumeration") == 0)
1479 : {
1480 : const uint8_t nVal =
1481 7 : static_cast<uint8_t>(poRAT->GetValueAsInt(iRow, iCol));
1482 7 : *pabyDst = nVal;
1483 : }
1484 13 : else if (strcmp(pszType, "int32") == 0 ||
1485 12 : strcmp(pszType, "uint32") == 0)
1486 : {
1487 5 : const int nVal = poRAT->GetValueAsInt(iRow, iCol);
1488 5 : memcpy(pabyDst, &nVal, sizeof(nVal));
1489 : }
1490 8 : else if (strcmp(pszType, "float32") == 0)
1491 : {
1492 : const float fVal =
1493 2 : static_cast<float>(poRAT->GetValueAsDouble(iRow, iCol));
1494 2 : memcpy(pabyDst, &fVal, sizeof(fVal));
1495 : }
1496 6 : else if (strcmp(pszType, "float64") == 0)
1497 : {
1498 1 : const double dfVal = poRAT->GetValueAsDouble(iRow, iCol);
1499 1 : memcpy(pabyDst, &dfVal, sizeof(dfVal));
1500 : }
1501 5 : else if (strcmp(pszType, "string") == 0)
1502 : {
1503 3 : asBuffers[iCol] = poRAT->GetValueAsString(iRow, iCol);
1504 3 : const char *pszStr = asBuffers[iCol].c_str();
1505 3 : memcpy(pabyDst, &pszStr, sizeof(pszStr));
1506 : }
1507 2 : else if (strcmp(pszType, "date") == 0)
1508 : {
1509 2 : asBuffers[iCol] = poRAT->GetValueAsString(iRow, iCol);
1510 2 : if (asBuffers[iCol].size() != 8)
1511 : {
1512 : OGRField sField;
1513 1 : if (OGRParseDate(asBuffers[iCol].c_str(), &sField, 0))
1514 : {
1515 2 : asBuffers[iCol] = CPLString().Printf(
1516 1 : "%04d%02d%02d", sField.Date.Year, sField.Date.Month,
1517 1 : sField.Date.Day);
1518 : }
1519 : }
1520 2 : const char *pszStr = asBuffers[iCol].c_str();
1521 2 : memcpy(pabyDst, &pszStr, sizeof(pszStr));
1522 : }
1523 : else
1524 : {
1525 0 : CPLAssert(false);
1526 : }
1527 : }
1528 :
1529 4 : H5OFFSET_TYPE offset[] = {static_cast<H5OFFSET_TYPE>(iRow)};
1530 4 : bRet =
1531 4 : bRet &&
1532 4 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1533 8 : nullptr, count, nullptr)) >= 0 &&
1534 4 : H5_CHECK(H5Dwrite(hDatasetID, hDataTypeMEM, hMemSpace, hFileSpace,
1535 : H5P_DEFAULT, abyBuffer.data())) >= 0;
1536 : }
1537 :
1538 6 : return bRet;
1539 : }
1540 :
1541 : /************************************************************************/
1542 : /* S102Creator::CreateGroupF() */
1543 : /************************************************************************/
1544 :
1545 : // Per S-102 v3.0 spec
1546 : #define MIN_DEPTH_VALUE -14
1547 : #define MAX_DEPTH_VALUE 11050
1548 :
1549 : #define STRINGIFY(x) #x
1550 : #define XSTRINGIFY(x) STRINGIFY(x)
1551 :
1552 15 : bool S102Creator::CreateGroupF(bool hasQualityOfBathymetryCoverage)
1553 : {
1554 15 : bool ret = S100BaseWriter::CreateGroupF();
1555 :
1556 15 : CPLStringList aosFeatureCodes;
1557 15 : aosFeatureCodes.push_back(FEATURE_TYPE);
1558 15 : if (hasQualityOfBathymetryCoverage)
1559 3 : aosFeatureCodes.push_back(QUALITY_FEATURE_TYPE);
1560 30 : ret = ret && WriteOneDimensionalVarLengthStringArray(
1561 15 : m_GroupF, "featureCode", aosFeatureCodes.List());
1562 :
1563 : {
1564 : std::vector<std::array<const char *, GROUP_F_DATASET_FIELD_COUNT>> rows{
1565 : {"depth", "depth", "metres", "1000000", "H5T_FLOAT",
1566 : XSTRINGIFY(MIN_DEPTH_VALUE), XSTRINGIFY(MAX_DEPTH_VALUE),
1567 : "closedInterval"},
1568 : {"uncertainty", "uncertainty", "metres", "1000000", "H5T_FLOAT",
1569 15 : "0", "", "geSemiInterval"}};
1570 15 : rows.resize(m_poSrcDS->GetRasterCount());
1571 15 : ret = ret && WriteGroupFDataset(FEATURE_TYPE, rows);
1572 : }
1573 : {
1574 : std::vector<std::array<const char *, GROUP_F_DATASET_FIELD_COUNT>> rows{
1575 15 : {"iD", "ID", "", "0", "H5T_INTEGER", "1", "", "geSemiInterval"}};
1576 15 : ret = ret && WriteGroupFDataset(QUALITY_FEATURE_TYPE, rows);
1577 : }
1578 :
1579 30 : return ret;
1580 : }
1581 :
1582 : /************************************************************************/
1583 : /* S102Creator::CopyValues() */
1584 : /************************************************************************/
1585 :
1586 16 : bool S102Creator::CopyValues(GDALProgressFunc pfnProgress, void *pProgressData)
1587 : {
1588 16 : CPLAssert(m_valuesGroup.get() >= 0);
1589 :
1590 16 : const int nYSize = m_poSrcDS->GetRasterYSize();
1591 16 : const int nXSize = m_poSrcDS->GetRasterXSize();
1592 :
1593 16 : hsize_t dims[] = {static_cast<hsize_t>(nYSize),
1594 16 : static_cast<hsize_t>(nXSize)};
1595 :
1596 32 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(2, dims, nullptr)));
1597 16 : bool bRet = hDataSpace;
1598 :
1599 : const bool bDeflate =
1600 16 : EQUAL(m_aosOptions.FetchNameValueDef("COMPRESS", "DEFLATE"), "DEFLATE");
1601 : const int nCompressionLevel =
1602 16 : atoi(m_aosOptions.FetchNameValueDef("ZLEVEL", "6"));
1603 : const int nBlockSize =
1604 16 : std::min(4096, std::max(100, atoi(m_aosOptions.FetchNameValueDef(
1605 16 : "BLOCK_SIZE", "100"))));
1606 16 : const int nBlockXSize = std::min(nXSize, nBlockSize);
1607 16 : const int nBlockYSize = std::min(nYSize, nBlockSize);
1608 16 : const float fNoDataValue = NODATA;
1609 16 : const int nComponents = m_poSrcDS->GetRasterCount();
1610 :
1611 : GH5_HIDTypeHolder hDataType(
1612 32 : H5_CHECK(H5Tcreate(H5T_COMPOUND, nComponents * sizeof(float))));
1613 16 : bRet = bRet && hDataType &&
1614 35 : H5_CHECK(H5Tinsert(hDataType, "depth", 0, H5T_IEEE_F32LE)) >= 0 &&
1615 3 : (nComponents == 1 ||
1616 3 : H5_CHECK(H5Tinsert(hDataType, "uncertainty", sizeof(float),
1617 : H5T_IEEE_F32LE)) >= 0);
1618 :
1619 16 : hsize_t chunk_size[] = {static_cast<hsize_t>(nBlockYSize),
1620 16 : static_cast<hsize_t>(nBlockXSize)};
1621 :
1622 16 : const float afFillValue[] = {fNoDataValue, fNoDataValue};
1623 32 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1624 16 : bRet = bRet && hParams &&
1625 16 : H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) >= 0 &&
1626 16 : H5_CHECK(H5Pset_fill_value(hParams, hDataType, afFillValue)) >= 0 &&
1627 48 : H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0 &&
1628 16 : H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) >= 0;
1629 :
1630 16 : if (bRet && bDeflate)
1631 : {
1632 15 : bRet = H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) >= 0;
1633 : }
1634 :
1635 32 : GH5_HIDDatasetHolder hDatasetID;
1636 16 : if (bRet)
1637 : {
1638 16 : hDatasetID.reset(H5_CHECK(H5Dcreate(m_valuesGroup, "values", hDataType,
1639 : hDataSpace, hParams)));
1640 16 : bRet = hDatasetID;
1641 : }
1642 :
1643 32 : GH5_HIDSpaceHolder hFileSpace;
1644 16 : if (bRet)
1645 : {
1646 16 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
1647 16 : bRet = hFileSpace;
1648 : }
1649 :
1650 16 : const int nYBlocks = static_cast<int>(DIV_ROUND_UP(nYSize, nBlockYSize));
1651 16 : const int nXBlocks = static_cast<int>(DIV_ROUND_UP(nXSize, nBlockXSize));
1652 16 : std::vector<float> afValues(static_cast<size_t>(nBlockYSize) * nBlockXSize *
1653 16 : nComponents);
1654 16 : const bool bReverseY = m_gt[5] < 0;
1655 :
1656 16 : float fMinDepth = std::numeric_limits<float>::infinity();
1657 16 : float fMaxDepth = -std::numeric_limits<float>::infinity();
1658 16 : float fMinUncertainty = std::numeric_limits<float>::infinity();
1659 16 : float fMaxUncertainty = -std::numeric_limits<float>::infinity();
1660 :
1661 16 : int bHasNoDataBand1 = FALSE;
1662 : const char *pszFirstBandDesc =
1663 16 : m_poSrcDS->GetRasterBand(1)->GetDescription();
1664 16 : const float fMulFactor =
1665 16 : EQUAL(pszFirstBandDesc, "elevation") ? -1.0f : 1.0f;
1666 16 : if (fMulFactor < 0.0f)
1667 : {
1668 1 : CPLError(CE_Warning, CPLE_AppDefined,
1669 : "Automatically convert from elevation to depth by negating "
1670 : "elevation values");
1671 : }
1672 : const double dfSrcNoDataBand1 =
1673 16 : m_poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoDataBand1);
1674 16 : const float fSrcNoDataBand1 = static_cast<float>(dfSrcNoDataBand1);
1675 16 : int bHasNoDataBand2 = FALSE;
1676 : const double dfSrcNoDataBand2 =
1677 : nComponents == 2
1678 16 : ? m_poSrcDS->GetRasterBand(2)->GetNoDataValue(&bHasNoDataBand2)
1679 16 : : 0.0;
1680 16 : const float fSrcNoDataBand2 = static_cast<float>(dfSrcNoDataBand2);
1681 :
1682 43 : for (int iY = 0; iY < nYBlocks && bRet; iY++)
1683 : {
1684 : const int nSrcYOff = bReverseY
1685 27 : ? std::max(0, nYSize - (iY + 1) * nBlockYSize)
1686 19 : : iY * nBlockYSize;
1687 27 : const int nReqCountY = std::min(nBlockYSize, nYSize - iY * nBlockYSize);
1688 186 : for (int iX = 0; iX < nXBlocks && bRet; iX++)
1689 : {
1690 : const int nReqCountX =
1691 159 : std::min(nBlockXSize, nXSize - iX * nBlockXSize);
1692 :
1693 159 : bRet =
1694 159 : m_poSrcDS->RasterIO(
1695 159 : GF_Read, iX * nBlockXSize, nSrcYOff, nReqCountX, nReqCountY,
1696 8 : bReverseY ? afValues.data() +
1697 8 : (nReqCountY - 1) * nReqCountX * nComponents
1698 151 : : afValues.data(),
1699 : nReqCountX, nReqCountY, GDT_Float32, nComponents, nullptr,
1700 159 : static_cast<int>(sizeof(float)) * nComponents,
1701 : bReverseY ? -static_cast<GPtrDiff_t>(sizeof(float)) *
1702 8 : nComponents * nReqCountX
1703 : : 0,
1704 : sizeof(float), nullptr) == CE_None;
1705 :
1706 159 : if (bRet)
1707 : {
1708 1441810 : for (int i = 0; i < nReqCountY * nReqCountX; i++)
1709 : {
1710 : {
1711 1441650 : float fVal = afValues[i * nComponents];
1712 2883290 : if ((bHasNoDataBand1 && fVal == fSrcNoDataBand1) ||
1713 1441650 : std::isnan(fVal))
1714 : {
1715 8 : afValues[i * nComponents] = fNoDataValue;
1716 : }
1717 : else
1718 : {
1719 1441640 : fVal *= fMulFactor;
1720 1441640 : afValues[i * nComponents] = fVal;
1721 1441640 : fMinDepth = std::min(fMinDepth, fVal);
1722 1441640 : fMaxDepth = std::max(fMaxDepth, fVal);
1723 : }
1724 : }
1725 1441650 : if (nComponents == 2)
1726 : {
1727 1440010 : const float fVal = afValues[i * nComponents + 1];
1728 2880020 : if ((bHasNoDataBand2 && fVal == fSrcNoDataBand2) ||
1729 1440010 : std::isnan(fVal))
1730 : {
1731 1 : afValues[i * nComponents + 1] = fNoDataValue;
1732 : }
1733 : else
1734 : {
1735 1440010 : fMinUncertainty = std::min(fMinUncertainty, fVal);
1736 1440010 : fMaxUncertainty = std::max(fMaxUncertainty, fVal);
1737 : }
1738 : }
1739 : }
1740 : }
1741 :
1742 : H5OFFSET_TYPE offset[] = {
1743 159 : static_cast<H5OFFSET_TYPE>(iY) *
1744 159 : static_cast<H5OFFSET_TYPE>(nBlockYSize),
1745 159 : static_cast<H5OFFSET_TYPE>(iX) *
1746 159 : static_cast<H5OFFSET_TYPE>(nBlockXSize)};
1747 159 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
1748 159 : static_cast<hsize_t>(nReqCountX)};
1749 : GH5_HIDSpaceHolder hMemSpace(
1750 159 : H5_CHECK(H5Screate_simple(2, count, nullptr)));
1751 159 : bRet =
1752 159 : bRet &&
1753 159 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1754 159 : nullptr, count, nullptr)) >= 0 &&
1755 159 : hMemSpace &&
1756 159 : H5_CHECK(H5Dwrite(hDatasetID, hDataType, hMemSpace, hFileSpace,
1757 318 : H5P_DEFAULT, afValues.data())) >= 0 &&
1758 159 : pfnProgress((static_cast<double>(iY) * nXBlocks + iX + 1) /
1759 159 : (static_cast<double>(nXBlocks) * nYBlocks),
1760 : "", pProgressData) != 0;
1761 : }
1762 : }
1763 :
1764 16 : if (fMinDepth > fMaxDepth)
1765 : {
1766 0 : fMinDepth = fMaxDepth = fNoDataValue;
1767 : }
1768 16 : else if (!(fMinDepth >= MIN_DEPTH_VALUE && fMaxDepth <= MAX_DEPTH_VALUE))
1769 : {
1770 2 : CPLError(CE_Warning, CPLE_AppDefined,
1771 : "Range of depth in the dataset is [%f, %f] whereas the "
1772 : "allowed range is [%d, %d]",
1773 : fMinDepth, fMaxDepth, MIN_DEPTH_VALUE, MAX_DEPTH_VALUE);
1774 : }
1775 :
1776 16 : if (fMinUncertainty > fMaxUncertainty)
1777 : {
1778 13 : fMinUncertainty = fMaxUncertainty = fNoDataValue;
1779 : }
1780 3 : else if (fMinUncertainty < 0)
1781 : {
1782 1 : CPLError(CE_Warning, CPLE_AppDefined,
1783 : "Negative uncertainty value found, which is not allowed");
1784 : }
1785 :
1786 16 : return bRet &&
1787 16 : WriteFloat32Value(m_valuesGroup, "minimumDepth", fMinDepth) &&
1788 16 : WriteFloat32Value(m_valuesGroup, "maximumDepth", fMaxDepth) &&
1789 16 : WriteFloat32Value(m_valuesGroup, "minimumUncertainty",
1790 32 : fMinUncertainty) &&
1791 16 : WriteFloat32Value(m_valuesGroup, "maximumUncertainty",
1792 32 : fMaxUncertainty);
1793 : }
1794 :
1795 : /************************************************************************/
1796 : /* S102Creator::CopyQualityValues() */
1797 : /************************************************************************/
1798 :
1799 3 : bool S102Creator::CopyQualityValues(GDALDataset *poQualityDS,
1800 : const std::set<int> &oSetRATId,
1801 : GDALProgressFunc pfnProgress,
1802 : void *pProgressData)
1803 : {
1804 3 : CPLAssert(m_valuesGroup.get() >= 0);
1805 :
1806 3 : const int nYSize = poQualityDS->GetRasterYSize();
1807 3 : const int nXSize = poQualityDS->GetRasterXSize();
1808 :
1809 3 : hsize_t dims[] = {static_cast<hsize_t>(nYSize),
1810 3 : static_cast<hsize_t>(nXSize)};
1811 :
1812 6 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(2, dims, nullptr)));
1813 3 : bool bRet = hDataSpace;
1814 :
1815 : const bool bDeflate =
1816 3 : EQUAL(m_aosOptions.FetchNameValueDef("COMPRESS", "DEFLATE"), "DEFLATE");
1817 : const int nCompressionLevel =
1818 3 : atoi(m_aosOptions.FetchNameValueDef("ZLEVEL", "6"));
1819 : const int nBlockSize =
1820 3 : std::min(4096, std::max(100, atoi(m_aosOptions.FetchNameValueDef(
1821 3 : "BLOCK_SIZE", "100"))));
1822 3 : const int nBlockXSize = std::min(nXSize, nBlockSize);
1823 3 : const int nBlockYSize = std::min(nYSize, nBlockSize);
1824 3 : constexpr uint32_t nNoDataValue = 0;
1825 :
1826 3 : hsize_t chunk_size[] = {static_cast<hsize_t>(nBlockYSize),
1827 3 : static_cast<hsize_t>(nBlockXSize)};
1828 :
1829 6 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1830 3 : bRet = bRet && hParams &&
1831 3 : H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) >= 0 &&
1832 3 : H5_CHECK(H5Pset_fill_value(hParams, H5T_STD_U32LE, &nNoDataValue)) >=
1833 3 : 0 &&
1834 9 : H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0 &&
1835 3 : H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) >= 0;
1836 :
1837 3 : if (bRet && bDeflate)
1838 : {
1839 3 : bRet = H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) >= 0;
1840 : }
1841 :
1842 6 : GH5_HIDDatasetHolder hDatasetID;
1843 3 : if (bRet)
1844 : {
1845 3 : hDatasetID.reset(H5_CHECK(H5Dcreate(
1846 : m_valuesGroup, "values", H5T_STD_U32LE, hDataSpace, hParams)));
1847 3 : bRet = hDatasetID;
1848 : }
1849 :
1850 6 : GH5_HIDSpaceHolder hFileSpace(H5_CHECK(H5Dget_space(hDatasetID)));
1851 3 : bRet = bRet && hFileSpace;
1852 :
1853 3 : const int nYBlocks = static_cast<int>(DIV_ROUND_UP(nYSize, nBlockYSize));
1854 3 : const int nXBlocks = static_cast<int>(DIV_ROUND_UP(nXSize, nBlockXSize));
1855 3 : std::vector<uint32_t> anValues(static_cast<size_t>(nBlockYSize) *
1856 6 : nBlockXSize);
1857 3 : const bool bReverseY = m_gt[5] < 0;
1858 :
1859 3 : int bHasSrcNoData = FALSE;
1860 : const double dfSrcNoData =
1861 3 : poQualityDS->GetRasterBand(1)->GetNoDataValue(&bHasSrcNoData);
1862 3 : const uint32_t nSrcNoData = static_cast<uint32_t>(dfSrcNoData);
1863 :
1864 3 : std::set<int> oSetRATIdCopy(oSetRATId);
1865 6 : for (int iY = 0; iY < nYBlocks && bRet; iY++)
1866 : {
1867 : const int nSrcYOff = bReverseY
1868 3 : ? std::max(0, nYSize - (iY + 1) * nBlockYSize)
1869 3 : : iY * nBlockYSize;
1870 3 : const int nReqCountY = std::min(nBlockYSize, nYSize - iY * nBlockYSize);
1871 6 : for (int iX = 0; iX < nXBlocks && bRet; iX++)
1872 : {
1873 : const int nReqCountX =
1874 3 : std::min(nBlockXSize, nXSize - iX * nBlockXSize);
1875 :
1876 3 : bRet =
1877 3 : poQualityDS->GetRasterBand(1)->RasterIO(
1878 3 : GF_Read, iX * nBlockXSize, nSrcYOff, nReqCountX, nReqCountY,
1879 0 : bReverseY ? anValues.data() + (nReqCountY - 1) * nReqCountX
1880 3 : : anValues.data(),
1881 : nReqCountX, nReqCountY, GDT_UInt32, 0,
1882 : bReverseY ? -static_cast<GPtrDiff_t>(sizeof(uint32_t)) *
1883 0 : nReqCountX
1884 : : 0,
1885 : nullptr) == CE_None;
1886 :
1887 3 : if (bRet)
1888 : {
1889 15 : for (int i = 0; i < nReqCountY * nReqCountX; i++)
1890 : {
1891 12 : if (bHasSrcNoData && anValues[i] == nSrcNoData)
1892 : {
1893 0 : anValues[i] = nNoDataValue;
1894 : }
1895 24 : else if (anValues[i] != 0 &&
1896 12 : !cpl::contains(oSetRATIdCopy, anValues[i]))
1897 : {
1898 1 : CPLError(
1899 : CE_Warning, CPLE_AppDefined,
1900 : "Quality grid contains nodes with id %u, but there "
1901 : "is no such entry in the feature attribute table",
1902 1 : anValues[i]);
1903 1 : oSetRATIdCopy.insert(anValues[i]);
1904 : }
1905 : }
1906 : }
1907 :
1908 : H5OFFSET_TYPE offset[] = {
1909 3 : static_cast<H5OFFSET_TYPE>(iY) *
1910 3 : static_cast<H5OFFSET_TYPE>(nBlockYSize),
1911 3 : static_cast<H5OFFSET_TYPE>(iX) *
1912 3 : static_cast<H5OFFSET_TYPE>(nBlockXSize)};
1913 3 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
1914 3 : static_cast<hsize_t>(nReqCountX)};
1915 3 : GH5_HIDSpaceHolder hMemSpace(H5Screate_simple(2, count, nullptr));
1916 3 : bRet =
1917 3 : bRet && hMemSpace &&
1918 3 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1919 3 : nullptr, count, nullptr)) >= 0 &&
1920 3 : H5_CHECK(H5Dwrite(hDatasetID, H5T_NATIVE_UINT, hMemSpace,
1921 : hFileSpace, H5P_DEFAULT, anValues.data())) >=
1922 6 : 0 &&
1923 3 : pfnProgress((static_cast<double>(iY) * nXBlocks + iX + 1) /
1924 3 : (static_cast<double>(nXBlocks) * nYBlocks),
1925 : "", pProgressData) != 0;
1926 : }
1927 : }
1928 :
1929 6 : return bRet;
1930 : }
1931 :
1932 : /************************************************************************/
1933 : /* S102Dataset::CreateCopy() */
1934 : /************************************************************************/
1935 :
1936 : /* static */
1937 56 : GDALDataset *S102Dataset::CreateCopy(const char *pszFilename,
1938 : GDALDataset *poSrcDS, int /* bStrict*/,
1939 : char **papszOptions,
1940 : GDALProgressFunc pfnProgress,
1941 : void *pProgressData)
1942 : {
1943 112 : S102Creator creator(pszFilename, poSrcDS, papszOptions);
1944 56 : if (!creator.Create(pfnProgress, pProgressData))
1945 40 : return nullptr;
1946 :
1947 : VSIStatBufL sStatBuf;
1948 32 : if (VSIStatL(pszFilename, &sStatBuf) == 0 &&
1949 16 : sStatBuf.st_size > 10 * 1024 * 1024)
1950 : {
1951 1 : CPLError(CE_Warning, CPLE_AppDefined,
1952 : "%s file size exceeds 10 MB, which is the upper limit "
1953 : "suggested for wireless transmission to marine vessels",
1954 : pszFilename);
1955 : }
1956 :
1957 32 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
1958 16 : return Open(&oOpenInfo);
1959 : }
1960 :
1961 : /************************************************************************/
1962 : /* S102DatasetDriverUnload() */
1963 : /************************************************************************/
1964 :
1965 9 : static void S102DatasetDriverUnload(GDALDriver *)
1966 : {
1967 9 : HDF5UnloadFileDriver();
1968 9 : }
1969 :
1970 : /************************************************************************/
1971 : /* GDALRegister_S102() */
1972 : /************************************************************************/
1973 14 : void GDALRegister_S102()
1974 :
1975 : {
1976 14 : if (!GDAL_CHECK_VERSION("S102"))
1977 0 : return;
1978 :
1979 14 : if (GDALGetDriverByName(S102_DRIVER_NAME) != nullptr)
1980 0 : return;
1981 :
1982 14 : GDALDriver *poDriver = new GDALDriver();
1983 :
1984 14 : S102DriverSetCommonMetadata(poDriver);
1985 14 : poDriver->pfnOpen = S102Dataset::Open;
1986 14 : poDriver->pfnCreateCopy = S102Dataset::CreateCopy;
1987 14 : poDriver->pfnUnloadDriver = S102DatasetDriverUnload;
1988 :
1989 14 : GetGDALDriverManager()->RegisterDriver(poDriver);
1990 : }
|
