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