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 : CSLConstList 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 : // 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
663 0 : oType.GetNumericDataType() == GDT_Int32))
664 : {
665 : // ok
666 : }
667 4 : else if (oType.GetClass() == GEDTC_COMPOUND &&
668 4 : oType.GetComponents().size() == 1 &&
669 2 : oType.GetComponents()[0]->GetType().GetClass() ==
670 4 : GEDTC_NUMERIC &&
671 2 : oType.GetComponents()[0]->GetType().GetNumericDataType() ==
672 : GDT_UInt32)
673 : {
674 : // seen in a S102 v3 product (102DE00CA22_UNC_MD.H5), although
675 : // I believe this is non-conformant.
676 :
677 : // Escape potentials single-quote and double-quote with back-slash
678 2 : CPLString osEscapedCompName(oType.GetComponents()[0]->GetName());
679 4 : osEscapedCompName.replaceAll("\\", "\\\\")
680 4 : .replaceAll("'", "\\'")
681 2 : .replaceAll("\"", "\\\"");
682 :
683 : // Gets a view with that single component extracted.
684 4 : poValuesArray = poValuesArray->GetView(
685 4 : std::string("['").append(osEscapedCompName).append("']"));
686 2 : if (!poValuesArray)
687 0 : return false;
688 : }
689 : else
690 : {
691 0 : CPLError(CE_Failure, CPLE_NotSupported,
692 : "Unsupported data type for %s",
693 0 : poValuesArray->GetFullName().c_str());
694 0 : return false;
695 : }
696 : }
697 :
698 7 : if (poValuesArray->GetDimensionCount() != 2)
699 : {
700 0 : CPLError(CE_Failure, CPLE_NotSupported,
701 : "Unsupported number of dimensions for %s",
702 0 : poValuesArray->GetFullName().c_str());
703 0 : return false;
704 : }
705 :
706 : auto poFeatureAttributeTable =
707 21 : poGroupQuality->OpenMDArray("featureAttributeTable");
708 7 : if (!poFeatureAttributeTable)
709 : {
710 0 : CPLError(CE_Failure, CPLE_AppDefined,
711 : "Cannot find array /%s/featureAttributeTable",
712 : pszNameOfQualityGroup);
713 0 : return false;
714 : }
715 :
716 : {
717 7 : const auto &oType = poFeatureAttributeTable->GetDataType();
718 7 : if (oType.GetClass() != GEDTC_COMPOUND)
719 : {
720 0 : CPLError(CE_Failure, CPLE_NotSupported,
721 : "Unsupported data type for %s",
722 0 : poFeatureAttributeTable->GetFullName().c_str());
723 0 : return false;
724 : }
725 :
726 7 : const auto &poComponents = oType.GetComponents();
727 7 : if (poComponents.size() >= 1 && poComponents[0]->GetName() != "id")
728 : {
729 0 : CPLError(CE_Failure, CPLE_AppDefined,
730 : "Missing 'id' component in %s",
731 0 : poFeatureAttributeTable->GetFullName().c_str());
732 0 : return false;
733 : }
734 : }
735 :
736 7 : if (bNorthUp)
737 5 : poValuesArray = poValuesArray->GetView("[::-1,...]");
738 :
739 : auto poDS =
740 14 : std::unique_ptr<GDALDataset>(poValuesArray->AsClassicDataset(1, 0));
741 7 : if (!poDS)
742 0 : return false;
743 :
744 7 : nRasterXSize = poDS->GetRasterXSize();
745 7 : nRasterYSize = poDS->GetRasterYSize();
746 :
747 : auto poRAT =
748 14 : HDF5CreateRAT(poFeatureAttributeTable, /* bFirstColIsMinMax = */ true);
749 : auto poBand = std::make_unique<S102GeoreferencedMetadataRasterBand>(
750 7 : std::move(poDS), std::move(poRAT));
751 7 : SetBand(1, poBand.release());
752 :
753 7 : return true;
754 : }
755 :
756 : /************************************************************************/
757 : /* S102Creator */
758 : /************************************************************************/
759 :
760 : class S102Creator final : public S100BaseWriter
761 : {
762 : public:
763 56 : S102Creator(const char *pszDestFilename, GDALDataset *poSrcDS,
764 : CSLConstList papszOptions)
765 56 : : S100BaseWriter(pszDestFilename, poSrcDS, papszOptions)
766 : {
767 56 : }
768 :
769 : ~S102Creator() override;
770 :
771 : bool Create(GDALProgressFunc pfnProgress, void *pProgressData);
772 :
773 : // From the S102 spec
774 : static constexpr float NODATA = 1000000.0f;
775 : static constexpr const char *FEATURE_TYPE = "BathymetryCoverage";
776 : static constexpr const char *QUALITY_FEATURE_TYPE =
777 : "QualityOfBathymetryCoverage";
778 :
779 : protected:
780 73 : bool Close() override
781 : {
782 73 : return BaseClose();
783 : }
784 :
785 : private:
786 : bool WriteFeatureGroupAttributes(bool isQuality);
787 : bool CopyValues(GDALProgressFunc pfnProgress, void *pProgressData);
788 : bool CopyQualityValues(GDALDataset *poQualityDS,
789 : const std::set<int> &oSetRATId,
790 : GDALProgressFunc pfnProgress, void *pProgressData);
791 : bool WriteFeatureAttributeTable(const GDALRasterAttributeTable *poRAT);
792 : bool CreateGroupF(bool hasQualityOfBathymetryCoverage);
793 : };
794 :
795 : /************************************************************************/
796 : /* S102Creator::~S102Creator() */
797 : /************************************************************************/
798 :
799 56 : S102Creator::~S102Creator()
800 : {
801 56 : S102Creator::Close();
802 56 : }
803 :
804 : /************************************************************************/
805 : /* S102Creator::Create() */
806 : /************************************************************************/
807 :
808 : // S102 v3.0 Table 10-8 - Elements of featureAttributeTable compound datatype
809 : static const struct
810 : {
811 : const char *pszName;
812 : const char *pszType;
813 : } gasFeatureAttributeTableMembers[] = {
814 : {"id", "uint32"},
815 : {"dataAssessment", "uint8"},
816 : {"featuresDetected.leastDepthOfDetectedFeaturesMeasured", "boolean"},
817 : {"featuresDetected.significantFeaturesDetected", "boolean"},
818 : {"featuresDetected.sizeOfFeaturesDetected", "float32"},
819 : {"featureSizeVar", "float32"},
820 : {"fullSeafloorCoverageAchieved", "boolean"},
821 : {"bathyCoverage", "boolean"},
822 : {"zoneOfConfidence.horizontalPositionUncertainty.uncertaintyFixed",
823 : "float32"},
824 : {"zoneOfConfidence.horizontalPositionUncertainty.uncertaintyVariableFactor",
825 : "float32"},
826 : {"surveyDateRange.dateStart", "date"},
827 : {"surveyDateRange.dateEnd", "date"},
828 : {"sourceSurveyID", "string"},
829 : {"surveyAuthority", "string"},
830 : {"typeOfBathymetricEstimationUncertainty", "enumeration"},
831 : };
832 :
833 56 : bool S102Creator::Create(GDALProgressFunc pfnProgress, void *pProgressData)
834 : {
835 56 : if (m_poSrcDS->GetRasterCount() != 1 && m_poSrcDS->GetRasterCount() != 2)
836 : {
837 4 : CPLError(CE_Failure, CPLE_NotSupported,
838 : "Source dataset must have one or two bands");
839 4 : return false;
840 : }
841 :
842 52 : if (!BaseChecks("S102", /* crsMustBeEPSG = */ true,
843 : /* verticalDatumRequired = */ true))
844 21 : return false;
845 :
846 : const bool bAppendSubdataset =
847 31 : CPLTestBool(m_aosOptions.FetchNameValueDef("APPEND_SUBDATASET", "NO"));
848 :
849 31 : std::unique_ptr<GDALDataset> poQualityDS;
850 : const char *pszQualityDataset =
851 31 : m_aosOptions.FetchNameValue("QUALITY_DATASET");
852 31 : const GDALRasterAttributeTable *poRAT = nullptr;
853 31 : if (!pszQualityDataset && !bAppendSubdataset)
854 : {
855 : const char *pszSubDSName =
856 13 : m_poSrcDS->GetMetadataItem("SUBDATASET_2_NAME", "SUBDATASETS");
857 13 : if (pszSubDSName &&
858 13 : cpl::starts_with(std::string_view(pszSubDSName), "S102:") &&
859 13 : cpl::ends_with(std::string_view(pszSubDSName),
860 : ":QualityOfBathymetryCoverage"))
861 : {
862 0 : pszQualityDataset = pszSubDSName;
863 : }
864 : }
865 :
866 62 : std::set<int> oSetRATId;
867 31 : if (pszQualityDataset)
868 : {
869 15 : if (bAppendSubdataset)
870 : {
871 0 : CPLError(CE_Failure, CPLE_NotSupported,
872 : "Quality dataset can only be set on initial creation");
873 12 : return false;
874 : }
875 15 : poQualityDS.reset(GDALDataset::Open(
876 : pszQualityDataset, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
877 : nullptr, nullptr));
878 15 : if (!poQualityDS)
879 1 : return false;
880 :
881 14 : if (poQualityDS->GetRasterCount() != 1)
882 : {
883 1 : CPLError(CE_Failure, CPLE_AppDefined,
884 : "%s does not have a single band.", pszQualityDataset);
885 1 : return false;
886 : }
887 13 : if (!GDALDataTypeIsInteger(
888 : poQualityDS->GetRasterBand(1)->GetRasterDataType()))
889 : {
890 1 : CPLError(CE_Failure, CPLE_AppDefined,
891 : "%s band is not of an integer data type.",
892 : pszQualityDataset);
893 1 : return false;
894 : }
895 23 : if (poQualityDS->GetRasterXSize() != m_poSrcDS->GetRasterXSize() ||
896 11 : poQualityDS->GetRasterYSize() != m_poSrcDS->GetRasterYSize())
897 : {
898 2 : CPLError(CE_Failure, CPLE_AppDefined,
899 : "%s does not have the same dimensions as %s.",
900 2 : pszQualityDataset, m_poSrcDS->GetDescription());
901 2 : return false;
902 : }
903 :
904 10 : const auto poQualityDS_SRS = poQualityDS->GetSpatialRef();
905 10 : if (!poQualityDS_SRS || !poQualityDS_SRS->IsSame(m_poSRS))
906 : {
907 1 : CPLError(CE_Failure, CPLE_AppDefined,
908 : "%s does not have the same CRS as %s.", pszQualityDataset,
909 1 : m_poSrcDS->GetDescription());
910 1 : return false;
911 : }
912 :
913 9 : GDALGeoTransform gt;
914 9 : if (poQualityDS->GetGeoTransform(gt) != CE_None || gt != m_gt)
915 : {
916 1 : CPLError(CE_Failure, CPLE_AppDefined,
917 : "%s does not have the same geotransform as %s.",
918 1 : pszQualityDataset, m_poSrcDS->GetDescription());
919 1 : return false;
920 : }
921 :
922 8 : poRAT = poQualityDS->GetRasterBand(1)->GetDefaultRAT();
923 8 : if (!poRAT)
924 : {
925 1 : CPLError(CE_Failure, CPLE_AppDefined,
926 : "%s does not have a raster attribute table.",
927 1 : poQualityDS->GetDescription());
928 1 : return false;
929 : }
930 :
931 7 : const int nRATColumnCount = poRAT->GetColumnCount();
932 7 : std::set<std::string_view> setKnownColumnNames;
933 112 : for (const auto &entry : gasFeatureAttributeTableMembers)
934 105 : setKnownColumnNames.insert(entry.pszName);
935 7 : int iRATIdField = -1;
936 24 : for (int i = 0; i < nRATColumnCount; ++i)
937 : {
938 17 : const char *pszColName = poRAT->GetNameOfCol(i);
939 17 : if (strcmp(pszColName, "id") == 0)
940 : {
941 6 : iRATIdField = i;
942 : }
943 11 : else if (!cpl::contains(setKnownColumnNames, pszColName))
944 : {
945 5 : CPLError(CE_Warning, CPLE_AppDefined,
946 : "'%s' is not a valid S102 feature attribute table "
947 : "column name.",
948 : pszColName);
949 : }
950 : }
951 7 : if (iRATIdField < 0)
952 : {
953 1 : CPLError(
954 : CE_Failure, CPLE_AppDefined,
955 : "Input raster attribute table lacks an integer 'id' field");
956 1 : return false;
957 : }
958 6 : const int nRATRowCount = poRAT->GetRowCount();
959 11 : for (int i = 0; i < nRATRowCount; ++i)
960 : {
961 8 : const int nID = poRAT->GetValueAsInt(i, iRATIdField);
962 8 : if (nID == 0)
963 : {
964 1 : CPLError(CE_Failure, CPLE_AppDefined,
965 : "id=0 is not allowed in input raster attribute table");
966 3 : return false;
967 : }
968 7 : else if (nID < 0)
969 : {
970 1 : CPLError(CE_Failure, CPLE_AppDefined,
971 : "Negative id is not allowed in input raster attribute "
972 : "table");
973 1 : return false;
974 : }
975 6 : else if (cpl::contains(oSetRATId, nID))
976 : {
977 1 : CPLError(
978 : CE_Failure, CPLE_AppDefined,
979 : "Several rows of input raster attribute table have id=%d",
980 : nID);
981 1 : return false;
982 : }
983 5 : oSetRATId.insert(nID);
984 : }
985 : }
986 :
987 19 : if (!((m_nVerticalDatum >= 1 && m_nVerticalDatum <= 30) ||
988 0 : m_nVerticalDatum == 44))
989 : {
990 0 : CPLError(CE_Warning, CPLE_AppDefined,
991 : "VERTICAL_DATUM=%d value is a valid S100 value but not "
992 : "allowed in S102. Valid values are [1, 30] or 44",
993 : m_nVerticalDatum);
994 : }
995 :
996 19 : if (!(m_nEPSGCode == 4326 || m_nEPSGCode == 5041 || m_nEPSGCode == 5042 ||
997 19 : (m_nEPSGCode >= 32601 && m_nEPSGCode <= 32660) ||
998 1 : (m_nEPSGCode >= 32701 && m_nEPSGCode <= 32760)))
999 : {
1000 1 : CPLError(CE_Warning, CPLE_NotSupported,
1001 : "The EPSG code of the CRS is %d. "
1002 : "Only EPSG codes 4326, 5041, 5042, [32601, 32660], "
1003 : "[32701, 32760] are officially supported. "
1004 : "The dataset may not be recognized by other software",
1005 : m_nEPSGCode);
1006 : }
1007 :
1008 19 : if (bAppendSubdataset)
1009 : {
1010 6 : GDALOpenInfo oOpenInfo(m_osDestFilename.c_str(), GA_ReadOnly);
1011 : auto poOriDS =
1012 6 : std::unique_ptr<GDALDataset>(S102Dataset::Open(&oOpenInfo));
1013 3 : if (!poOriDS)
1014 : {
1015 1 : CPLError(CE_Failure, CPLE_AppDefined,
1016 : "%s is not a valid existing S102 dataset",
1017 : m_osDestFilename.c_str());
1018 1 : return false;
1019 : }
1020 2 : const auto poOriSRS = poOriDS->GetSpatialRef();
1021 2 : if (!poOriSRS)
1022 : {
1023 : // shouldn't happen
1024 0 : return false;
1025 : }
1026 2 : if (!poOriSRS->IsSame(m_poSRS))
1027 : {
1028 1 : CPLError(CE_Failure, CPLE_AppDefined,
1029 : "CRS of %s is not the same as the one of %s",
1030 1 : m_osDestFilename.c_str(), m_poSrcDS->GetDescription());
1031 1 : return false;
1032 : }
1033 1 : poOriDS.reset();
1034 :
1035 1 : OGREnvelope sExtent;
1036 1 : if (m_poSrcDS->GetExtentWGS84LongLat(&sExtent) != OGRERR_NONE)
1037 : {
1038 0 : CPLError(CE_Failure, CPLE_AppDefined,
1039 : "Cannot get dataset extent in WGS84 longitude/latitude");
1040 0 : return false;
1041 : }
1042 :
1043 1 : bool ret = OpenFileUpdateMode();
1044 1 : if (ret)
1045 : {
1046 1 : m_featureGroup.reset(
1047 : H5_CHECK(H5Gopen(m_hdf5, "BathymetryCoverage")));
1048 : }
1049 :
1050 1 : ret = ret && m_featureGroup;
1051 1 : double dfNumInstances = 0;
1052 1 : ret = ret && GH5_FetchAttribute(m_featureGroup, "numInstances",
1053 : dfNumInstances, true);
1054 1 : if (ret && !(dfNumInstances >= 1 && dfNumInstances <= 99 &&
1055 1 : std::round(dfNumInstances) == dfNumInstances))
1056 : {
1057 0 : CPLError(CE_Failure, CPLE_AppDefined,
1058 : "Invalid value for numInstances");
1059 0 : ret = false;
1060 : }
1061 1 : else if (ret && dfNumInstances == 99)
1062 : {
1063 0 : CPLError(CE_Failure, CPLE_AppDefined,
1064 : "Too many existing feature instances");
1065 0 : ret = false;
1066 : }
1067 : else
1068 : {
1069 1 : double dfMainVerticalDatum = 0;
1070 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "verticalDatum",
1071 : dfMainVerticalDatum, true);
1072 :
1073 1 : const int newNumInstances = static_cast<int>(dfNumInstances) + 1;
1074 1 : ret = ret && GH5_WriteAttribute(m_featureGroup, "numInstances",
1075 : newNumInstances);
1076 1 : ret = ret && CreateFeatureInstanceGroup(CPLSPrintf(
1077 : "BathymetryCoverage.%02d", newNumInstances));
1078 1 : ret = ret && WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1079 1 : if (dfMainVerticalDatum != m_nVerticalDatum)
1080 : {
1081 2 : ret = ret &&
1082 1 : GH5_CreateAttribute(m_featureInstanceGroup,
1083 : "verticalDatumReference",
1084 2 : H5T_STD_U8LE) &&
1085 : // s100VerticalDatum
1086 1 : GH5_WriteAttribute(m_featureInstanceGroup,
1087 : "verticalDatumReference", 1);
1088 1 : ret =
1089 2 : ret && WriteVerticalDatum(m_featureInstanceGroup,
1090 1 : H5T_STD_U16LE, m_nVerticalDatum);
1091 : }
1092 :
1093 1 : ret = ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U8LE, 1);
1094 :
1095 1 : ret = ret && CreateValuesGroup("Group_001");
1096 1 : ret = ret && WriteVarLengthStringValue(m_valuesGroup, "timePoint",
1097 : "00010101T000000Z");
1098 1 : ret = ret && CopyValues(pfnProgress, pProgressData);
1099 : }
1100 :
1101 : // Update global bounding box
1102 1 : OGREnvelope sExistingExtent;
1103 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "westBoundLongitude",
1104 : sExistingExtent.MinX, true);
1105 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "southBoundLatitude",
1106 : sExistingExtent.MinY, true);
1107 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "eastBoundLongitude",
1108 : sExistingExtent.MaxX, true);
1109 1 : ret = ret && GH5_FetchAttribute(m_hdf5, "northBoundLatitude",
1110 : sExistingExtent.MaxY, true);
1111 :
1112 1 : sExtent.Merge(sExistingExtent);
1113 2 : ret = ret &&
1114 1 : GH5_WriteAttribute(m_hdf5, "westBoundLongitude", sExtent.MinX);
1115 2 : ret = ret &&
1116 1 : GH5_WriteAttribute(m_hdf5, "southBoundLatitude", sExtent.MinY);
1117 2 : ret = ret &&
1118 1 : GH5_WriteAttribute(m_hdf5, "eastBoundLongitude", sExtent.MaxX);
1119 2 : ret = ret &&
1120 1 : GH5_WriteAttribute(m_hdf5, "northBoundLatitude", sExtent.MaxY);
1121 :
1122 1 : return Close() && ret;
1123 : }
1124 : else
1125 : {
1126 16 : bool ret = CreateFile();
1127 16 : ret = ret && WriteProductSpecification("INT.IHO.S-102.3.0.0");
1128 16 : ret = ret && WriteIssueDate();
1129 16 : ret = ret && WriteIssueTime(/* bAutogenerateFromCurrent = */ false);
1130 16 : ret = ret && WriteHorizontalCRS();
1131 16 : ret = ret && WriteTopLevelBoundingBox();
1132 16 : ret = ret && WriteVerticalCS(6498); // Depth, metre, down
1133 16 : ret = ret && WriteVerticalCoordinateBase(2); // verticalDatum
1134 : // s100VerticalDatum
1135 16 : ret = ret && WriteVerticalDatumReference(m_hdf5, 1);
1136 16 : ret =
1137 16 : ret && WriteVerticalDatum(m_hdf5, H5T_STD_U16LE, m_nVerticalDatum);
1138 :
1139 : // BathymetryCoverage
1140 16 : ret = ret && CreateFeatureGroup(FEATURE_TYPE);
1141 16 : ret = ret && WriteFeatureGroupAttributes(/* isQuality = */ false);
1142 16 : ret = ret && WriteAxisNames(m_featureGroup);
1143 :
1144 16 : ret = ret && CreateFeatureInstanceGroup("BathymetryCoverage.01");
1145 16 : ret = ret && WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1146 16 : ret = ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U8LE, 1);
1147 :
1148 16 : ret = ret && CreateValuesGroup("Group_001");
1149 :
1150 16 : ret = ret && WriteVarLengthStringValue(m_valuesGroup, "timePoint",
1151 : "00010101T000000Z");
1152 :
1153 : const double dfIntermediatePct =
1154 16 : m_poSrcDS->GetRasterCount() /
1155 16 : (m_poSrcDS->GetRasterCount() + (poQualityDS ? 1.0 : 0.0));
1156 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>
1157 : pScaledProgressData(GDALCreateScaledProgress(0.0, dfIntermediatePct,
1158 : pfnProgress,
1159 : pProgressData),
1160 16 : GDALDestroyScaledProgress);
1161 16 : ret = ret && CopyValues(GDALScaledProgress, pScaledProgressData.get());
1162 :
1163 16 : if (poQualityDS)
1164 : {
1165 : // QualityOfBathymetryCoverage group
1166 3 : ret = ret && CreateFeatureGroup(QUALITY_FEATURE_TYPE);
1167 3 : ret = ret && WriteFeatureGroupAttributes(/* isQuality = */ true);
1168 3 : ret = ret && WriteAxisNames(m_featureGroup);
1169 3 : ret = ret && WriteFeatureAttributeTable(poRAT);
1170 :
1171 6 : ret = ret &&
1172 3 : CreateFeatureInstanceGroup("QualityOfBathymetryCoverage.01");
1173 3 : ret = ret && WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1174 3 : ret = ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U8LE, 1);
1175 :
1176 3 : ret = ret && CreateValuesGroup("Group_001");
1177 3 : pScaledProgressData.reset(GDALCreateScaledProgress(
1178 : dfIntermediatePct, 1.0, pfnProgress, pProgressData));
1179 3 : ret = ret && CopyQualityValues(poQualityDS.get(), oSetRATId,
1180 : GDALScaledProgress,
1181 : pScaledProgressData.get());
1182 : }
1183 :
1184 16 : ret = ret && CreateGroupF(poQualityDS != nullptr);
1185 :
1186 16 : return Close() && ret;
1187 : }
1188 : }
1189 :
1190 : /************************************************************************/
1191 : /* S102Creator::WriteFeatureGroupAttributes() */
1192 : /************************************************************************/
1193 :
1194 18 : bool S102Creator::WriteFeatureGroupAttributes(bool isQuality)
1195 : {
1196 18 : CPLAssert(m_featureGroup);
1197 :
1198 18 : bool ret = WriteCommonPointRule(m_featureGroup, 2); // low
1199 18 : if (isQuality)
1200 : {
1201 : // Feature oriented Regular Grid
1202 3 : ret = ret && WriteDataCodingFormat(m_featureGroup, 9);
1203 : }
1204 : else
1205 : {
1206 15 : ret = ret && WriteDataCodingFormat(m_featureGroup, 2); // Regular grid
1207 : }
1208 18 : ret = ret && WriteDataOffsetCode(m_featureGroup, 5); // Center of cell
1209 18 : ret = ret && WriteDimension(m_featureGroup, 2);
1210 : const char *pszHorizontalPositionUncertainty =
1211 18 : m_aosOptions.FetchNameValue("HORIZONTAL_POSITION_UNCERTAINTY");
1212 18 : ret =
1213 36 : ret &&
1214 18 : WriteHorizontalPositionUncertainty(
1215 : m_featureGroup,
1216 0 : pszHorizontalPositionUncertainty &&
1217 0 : pszHorizontalPositionUncertainty[0]
1218 0 : ? static_cast<float>(CPLAtof(pszHorizontalPositionUncertainty))
1219 : : -1.0f);
1220 : const char *pszVerticalUncertainty =
1221 18 : m_aosOptions.FetchNameValue("VERTICAL_UNCERTAINTY");
1222 18 : ret = ret && WriteVerticalUncertainty(
1223 : m_featureGroup,
1224 0 : pszVerticalUncertainty && pszVerticalUncertainty[0]
1225 0 : ? static_cast<float>(CPLAtof(pszVerticalUncertainty))
1226 : : -1.0f);
1227 18 : ret = ret && WriteInterpolationType(m_featureGroup, 1); // Nearest neighbor
1228 18 : ret = ret && WriteNumInstances(m_featureGroup, H5T_STD_U8LE, 1);
1229 36 : ret = ret && WriteSequencingRuleScanDirection(m_featureGroup,
1230 18 : m_poSRS->IsProjected()
1231 : ? "Easting, Northing"
1232 : : "Longitude, Latitude");
1233 18 : ret = ret && WriteSequencingRuleType(m_featureGroup, 1); // Linear
1234 18 : return ret;
1235 : }
1236 :
1237 : /************************************************************************/
1238 : /* S102Creator::WriteFeatureAttributeTable() */
1239 : /************************************************************************/
1240 :
1241 3 : bool S102Creator::WriteFeatureAttributeTable(
1242 : const GDALRasterAttributeTable *poRAT)
1243 : {
1244 3 : CPLAssert(m_featureGroup);
1245 :
1246 6 : std::map<std::string_view, const char *> mapKnownColumns;
1247 48 : for (const auto &entry : gasFeatureAttributeTableMembers)
1248 45 : mapKnownColumns[entry.pszName] = entry.pszType;
1249 :
1250 3 : const int nColCount = poRAT->GetColumnCount();
1251 :
1252 3 : size_t nCompoundSize = 0;
1253 3 : size_t nMEMCompoundSize = 0;
1254 16 : for (int i = 0; i < nColCount; ++i)
1255 : {
1256 13 : const char *pszColName = poRAT->GetNameOfCol(i);
1257 13 : const auto iter = mapKnownColumns.find(pszColName);
1258 13 : size_t nMemberSize = sizeof(char *);
1259 13 : if (iter != mapKnownColumns.end())
1260 : {
1261 9 : const char *pszType = iter->second;
1262 9 : if (strcmp(pszType, "uint8") == 0 ||
1263 8 : strcmp(pszType, "boolean") == 0 ||
1264 7 : strcmp(pszType, "enumeration") == 0)
1265 : {
1266 3 : nMemberSize = sizeof(uint8_t);
1267 : }
1268 6 : else if (strcmp(pszType, "uint32") == 0)
1269 : {
1270 3 : nMemberSize = sizeof(uint32_t);
1271 : }
1272 3 : else if (strcmp(pszType, "float32") == 0)
1273 : {
1274 1 : nMemberSize = sizeof(float);
1275 : }
1276 2 : else if (strcmp(pszType, "string") == 0 ||
1277 1 : strcmp(pszType, "date") == 0)
1278 : {
1279 2 : nMemberSize = sizeof(char *);
1280 : }
1281 : else
1282 : {
1283 0 : CPLAssert(false);
1284 : }
1285 : }
1286 : else
1287 : {
1288 4 : GDALRATFieldType eType = poRAT->GetTypeOfCol(i);
1289 4 : switch (eType)
1290 : {
1291 1 : case GFT_Integer:
1292 1 : nMemberSize = sizeof(int32_t);
1293 1 : break;
1294 1 : case GFT_Real:
1295 1 : nMemberSize = sizeof(double);
1296 1 : break;
1297 1 : case GFT_Boolean:
1298 1 : nMemberSize = sizeof(uint8_t);
1299 1 : break;
1300 1 : case GFT_String:
1301 : case GFT_DateTime:
1302 : case GFT_WKBGeometry:
1303 1 : nMemberSize = sizeof(char *);
1304 1 : break;
1305 : }
1306 : }
1307 13 : nCompoundSize += nMemberSize;
1308 13 : if ((nMEMCompoundSize % nMemberSize) != 0)
1309 2 : nMEMCompoundSize += nMemberSize - (nMEMCompoundSize % nMemberSize);
1310 13 : nMEMCompoundSize += nMemberSize;
1311 : }
1312 :
1313 : GH5_HIDTypeHolder hDataType(
1314 6 : H5_CHECK(H5Tcreate(H5T_COMPOUND, nCompoundSize)));
1315 : GH5_HIDTypeHolder hDataTypeMEM(
1316 6 : H5_CHECK(H5Tcreate(H5T_COMPOUND, nMEMCompoundSize)));
1317 6 : GH5_HIDTypeHolder hVarLengthType(H5_CHECK(H5Tcopy(H5T_C_S1)));
1318 6 : bool bRet = hDataType && hDataTypeMEM && hVarLengthType &&
1319 9 : H5_CHECK(H5Tset_size(hVarLengthType, H5T_VARIABLE)) >= 0 &&
1320 3 : H5_CHECK(H5Tset_strpad(hVarLengthType, H5T_STR_NULLTERM)) >= 0;
1321 :
1322 6 : GH5_HIDTypeHolder hEnumType;
1323 6 : std::vector<const char *> apszTypes;
1324 :
1325 3 : size_t nOffset = 0;
1326 3 : size_t nMEMOffset = 0;
1327 6 : std::vector<size_t> anMEMOffsets;
1328 16 : for (int i = 0; i < nColCount && bRet; ++i)
1329 : {
1330 13 : const char *pszColName = poRAT->GetNameOfCol(i);
1331 13 : const auto iter = mapKnownColumns.find(pszColName);
1332 13 : hid_t hMemberType = hVarLengthType.get();
1333 13 : hid_t hMemberNativeType = hVarLengthType.get();
1334 13 : if (iter != mapKnownColumns.end())
1335 : {
1336 9 : const char *pszType = iter->second;
1337 9 : if (strcmp(pszType, "uint8") == 0 ||
1338 8 : strcmp(pszType, "boolean") == 0)
1339 : {
1340 2 : hMemberType = H5T_STD_U8LE;
1341 2 : hMemberNativeType = H5T_NATIVE_UCHAR;
1342 : }
1343 7 : else if (strcmp(pszType, "uint32") == 0)
1344 : {
1345 3 : hMemberType = H5T_STD_U32LE;
1346 3 : hMemberNativeType = H5T_NATIVE_UINT;
1347 : }
1348 4 : else if (strcmp(pszType, "float32") == 0)
1349 : {
1350 1 : hMemberType = H5T_IEEE_F32LE;
1351 1 : hMemberNativeType = H5T_NATIVE_FLOAT;
1352 : }
1353 3 : else if (strcmp(pszType, "string") == 0 ||
1354 2 : strcmp(pszType, "date") == 0)
1355 : {
1356 2 : hMemberType = hVarLengthType.get();
1357 2 : hMemberNativeType = hVarLengthType.get();
1358 : }
1359 1 : else if (strcmp(pszType, "enumeration") == 0 &&
1360 1 : strcmp(pszColName,
1361 : "typeOfBathymetricEstimationUncertainty") == 0)
1362 : {
1363 1 : hEnumType.reset(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1364 1 : bRet = hEnumType;
1365 1 : if (bRet)
1366 : {
1367 : uint8_t val;
1368 1 : val = 1;
1369 2 : bRet = bRet &&
1370 1 : H5_CHECK(H5Tenum_insert(
1371 : hEnumType, "rawStandardDeviation", &val)) >= 0;
1372 1 : val = 2;
1373 2 : bRet = bRet &&
1374 1 : H5_CHECK(H5Tenum_insert(
1375 : hEnumType, "cUBEStandardDeviation", &val)) >= 0;
1376 1 : val = 3;
1377 2 : bRet = bRet &&
1378 1 : H5_CHECK(H5Tenum_insert(
1379 : hEnumType, "productUncertainty", &val)) >= 0;
1380 1 : val = 4;
1381 1 : bRet = bRet && H5_CHECK(H5Tenum_insert(
1382 : hEnumType, "historicalStandardDeviation",
1383 : &val)) >= 0;
1384 :
1385 1 : hMemberType = hEnumType.get();
1386 1 : hMemberNativeType = hEnumType.get();
1387 1 : }
1388 : }
1389 : else
1390 : {
1391 0 : CPLAssert(false);
1392 : }
1393 9 : apszTypes.push_back(pszType);
1394 : }
1395 : else
1396 : {
1397 4 : GDALRATFieldType eType = poRAT->GetTypeOfCol(i);
1398 4 : switch (eType)
1399 : {
1400 1 : case GFT_Integer:
1401 1 : hMemberType = H5T_STD_I32LE;
1402 1 : hMemberNativeType = H5T_NATIVE_INT;
1403 1 : apszTypes.push_back("int32");
1404 1 : break;
1405 1 : case GFT_Real:
1406 1 : hMemberType = H5T_IEEE_F64LE;
1407 1 : hMemberNativeType = H5T_NATIVE_DOUBLE;
1408 1 : apszTypes.push_back("float64");
1409 1 : break;
1410 1 : case GFT_Boolean:
1411 1 : hMemberType = H5T_STD_U8LE;
1412 1 : hMemberNativeType = H5T_NATIVE_UCHAR;
1413 1 : apszTypes.push_back("boolean");
1414 1 : break;
1415 1 : case GFT_String:
1416 : case GFT_DateTime:
1417 : case GFT_WKBGeometry:
1418 1 : apszTypes.push_back("string");
1419 1 : break;
1420 : }
1421 : }
1422 :
1423 13 : CPLAssert(H5Tget_size(hMemberType) == H5Tget_size(hMemberNativeType));
1424 :
1425 13 : bRet = bRet && H5_CHECK(H5Tinsert(hDataType, pszColName, nOffset,
1426 : hMemberType)) >= 0;
1427 :
1428 13 : const size_t nMemberSize = H5Tget_size(hMemberType);
1429 13 : if ((nMEMOffset % nMemberSize) != 0)
1430 2 : nMEMOffset += nMemberSize - (nMEMOffset % nMemberSize);
1431 13 : anMEMOffsets.push_back(nMEMOffset);
1432 13 : bRet = bRet && H5_CHECK(H5Tinsert(hDataTypeMEM, pszColName, nMEMOffset,
1433 : hMemberNativeType)) >= 0;
1434 13 : nOffset += nMemberSize;
1435 13 : nMEMOffset += nMemberSize;
1436 : }
1437 3 : CPLAssert(nOffset == nCompoundSize);
1438 3 : CPLAssert(nMEMOffset == nMEMCompoundSize);
1439 :
1440 3 : CPLAssert(apszTypes.size() == static_cast<size_t>(nColCount));
1441 :
1442 3 : const int nRowCount = poRAT->GetRowCount();
1443 3 : hsize_t dims[] = {static_cast<hsize_t>(nRowCount)};
1444 6 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(1, dims, nullptr)));
1445 3 : bRet = bRet && hDataSpace;
1446 6 : GH5_HIDDatasetHolder hDatasetID;
1447 6 : GH5_HIDSpaceHolder hFileSpace;
1448 6 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1449 3 : bRet = bRet && hParams;
1450 3 : if (bRet)
1451 : {
1452 3 : bRet = H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0;
1453 3 : hsize_t chunk_size[] = {static_cast<hsize_t>(1)};
1454 3 : bRet = bRet && H5_CHECK(H5Pset_chunk(hParams, 1, chunk_size)) >= 0;
1455 3 : hDatasetID.reset(
1456 : H5_CHECK(H5Dcreate(m_featureGroup, "featureAttributeTable",
1457 : hDataType, hDataSpace, hParams)));
1458 3 : bRet = bRet && hDatasetID;
1459 : }
1460 3 : if (bRet)
1461 : {
1462 3 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
1463 3 : bRet = hFileSpace;
1464 : }
1465 :
1466 3 : hsize_t count[] = {1};
1467 6 : GH5_HIDSpaceHolder hMemSpace(H5_CHECK(H5Screate_simple(1, count, nullptr)));
1468 3 : bRet = bRet && hMemSpace;
1469 :
1470 6 : std::vector<GByte> abyBuffer(nMEMCompoundSize);
1471 3 : std::vector<std::string> asBuffers(nColCount);
1472 7 : for (int iRow = 0; iRow < nRowCount && bRet; ++iRow)
1473 : {
1474 24 : for (int iCol = 0; iCol < nColCount && bRet; ++iCol)
1475 : {
1476 20 : const char *const pszType = apszTypes[iCol];
1477 20 : GByte *const pabyDst = abyBuffer.data() + anMEMOffsets[iCol];
1478 20 : if (strcmp(pszType, "uint8") == 0 ||
1479 18 : strcmp(pszType, "boolean") == 0 ||
1480 15 : strcmp(pszType, "enumeration") == 0)
1481 : {
1482 : const uint8_t nVal =
1483 7 : static_cast<uint8_t>(poRAT->GetValueAsInt(iRow, iCol));
1484 7 : *pabyDst = nVal;
1485 : }
1486 13 : else if (strcmp(pszType, "int32") == 0 ||
1487 12 : strcmp(pszType, "uint32") == 0)
1488 : {
1489 5 : const int nVal = poRAT->GetValueAsInt(iRow, iCol);
1490 5 : memcpy(pabyDst, &nVal, sizeof(nVal));
1491 : }
1492 8 : else if (strcmp(pszType, "float32") == 0)
1493 : {
1494 : const float fVal =
1495 2 : static_cast<float>(poRAT->GetValueAsDouble(iRow, iCol));
1496 2 : memcpy(pabyDst, &fVal, sizeof(fVal));
1497 : }
1498 6 : else if (strcmp(pszType, "float64") == 0)
1499 : {
1500 1 : const double dfVal = poRAT->GetValueAsDouble(iRow, iCol);
1501 1 : memcpy(pabyDst, &dfVal, sizeof(dfVal));
1502 : }
1503 5 : else if (strcmp(pszType, "string") == 0)
1504 : {
1505 3 : asBuffers[iCol] = poRAT->GetValueAsString(iRow, iCol);
1506 3 : const char *pszStr = asBuffers[iCol].c_str();
1507 3 : memcpy(pabyDst, &pszStr, sizeof(pszStr));
1508 : }
1509 2 : else if (strcmp(pszType, "date") == 0)
1510 : {
1511 2 : asBuffers[iCol] = poRAT->GetValueAsString(iRow, iCol);
1512 2 : if (asBuffers[iCol].size() != 8)
1513 : {
1514 : OGRField sField;
1515 1 : if (OGRParseDate(asBuffers[iCol].c_str(), &sField, 0))
1516 : {
1517 2 : asBuffers[iCol] = CPLString().Printf(
1518 1 : "%04d%02d%02d", sField.Date.Year, sField.Date.Month,
1519 1 : sField.Date.Day);
1520 : }
1521 : }
1522 2 : const char *pszStr = asBuffers[iCol].c_str();
1523 2 : memcpy(pabyDst, &pszStr, sizeof(pszStr));
1524 : }
1525 : else
1526 : {
1527 0 : CPLAssert(false);
1528 : }
1529 : }
1530 :
1531 4 : H5OFFSET_TYPE offset[] = {static_cast<H5OFFSET_TYPE>(iRow)};
1532 4 : bRet =
1533 4 : bRet &&
1534 4 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1535 8 : nullptr, count, nullptr)) >= 0 &&
1536 4 : H5_CHECK(H5Dwrite(hDatasetID, hDataTypeMEM, hMemSpace, hFileSpace,
1537 : H5P_DEFAULT, abyBuffer.data())) >= 0;
1538 : }
1539 :
1540 6 : return bRet;
1541 : }
1542 :
1543 : /************************************************************************/
1544 : /* S102Creator::CreateGroupF() */
1545 : /************************************************************************/
1546 :
1547 : // Per S-102 v3.0 spec
1548 : #define MIN_DEPTH_VALUE -14
1549 : #define MAX_DEPTH_VALUE 11050
1550 :
1551 : #define STRINGIFY(x) #x
1552 : #define XSTRINGIFY(x) STRINGIFY(x)
1553 :
1554 15 : bool S102Creator::CreateGroupF(bool hasQualityOfBathymetryCoverage)
1555 : {
1556 15 : bool ret = S100BaseWriter::CreateGroupF();
1557 :
1558 15 : CPLStringList aosFeatureCodes;
1559 15 : aosFeatureCodes.push_back(FEATURE_TYPE);
1560 15 : if (hasQualityOfBathymetryCoverage)
1561 3 : aosFeatureCodes.push_back(QUALITY_FEATURE_TYPE);
1562 30 : ret = ret && WriteOneDimensionalVarLengthStringArray(
1563 15 : m_GroupF, "featureCode", aosFeatureCodes.List());
1564 :
1565 : {
1566 : std::vector<std::array<const char *, GROUP_F_DATASET_FIELD_COUNT>> rows{
1567 : {"depth", "depth", "metres", "1000000", "H5T_FLOAT",
1568 : XSTRINGIFY(MIN_DEPTH_VALUE), XSTRINGIFY(MAX_DEPTH_VALUE),
1569 : "closedInterval"},
1570 : {"uncertainty", "uncertainty", "metres", "1000000", "H5T_FLOAT",
1571 15 : "0", "", "geSemiInterval"}};
1572 15 : rows.resize(m_poSrcDS->GetRasterCount());
1573 15 : ret = ret && WriteGroupFDataset(FEATURE_TYPE, rows);
1574 : }
1575 : {
1576 : std::vector<std::array<const char *, GROUP_F_DATASET_FIELD_COUNT>> rows{
1577 15 : {"iD", "ID", "", "0", "H5T_INTEGER", "1", "", "geSemiInterval"}};
1578 15 : ret = ret && WriteGroupFDataset(QUALITY_FEATURE_TYPE, rows);
1579 : }
1580 :
1581 30 : return ret;
1582 : }
1583 :
1584 : /************************************************************************/
1585 : /* S102Creator::CopyValues() */
1586 : /************************************************************************/
1587 :
1588 16 : bool S102Creator::CopyValues(GDALProgressFunc pfnProgress, void *pProgressData)
1589 : {
1590 16 : CPLAssert(m_valuesGroup.get() >= 0);
1591 :
1592 16 : const int nYSize = m_poSrcDS->GetRasterYSize();
1593 16 : const int nXSize = m_poSrcDS->GetRasterXSize();
1594 :
1595 16 : hsize_t dims[] = {static_cast<hsize_t>(nYSize),
1596 16 : static_cast<hsize_t>(nXSize)};
1597 :
1598 32 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(2, dims, nullptr)));
1599 16 : bool bRet = hDataSpace;
1600 :
1601 : const bool bDeflate =
1602 16 : EQUAL(m_aosOptions.FetchNameValueDef("COMPRESS", "DEFLATE"), "DEFLATE");
1603 : const int nCompressionLevel =
1604 16 : atoi(m_aosOptions.FetchNameValueDef("ZLEVEL", "6"));
1605 : const int nBlockSize =
1606 16 : std::min(4096, std::max(100, atoi(m_aosOptions.FetchNameValueDef(
1607 16 : "BLOCK_SIZE", "100"))));
1608 16 : const int nBlockXSize = std::min(nXSize, nBlockSize);
1609 16 : const int nBlockYSize = std::min(nYSize, nBlockSize);
1610 16 : const float fNoDataValue = NODATA;
1611 16 : const int nComponents = m_poSrcDS->GetRasterCount();
1612 :
1613 : GH5_HIDTypeHolder hDataType(
1614 32 : H5_CHECK(H5Tcreate(H5T_COMPOUND, nComponents * sizeof(float))));
1615 16 : bRet = bRet && hDataType &&
1616 35 : H5_CHECK(H5Tinsert(hDataType, "depth", 0, H5T_IEEE_F32LE)) >= 0 &&
1617 3 : (nComponents == 1 ||
1618 3 : H5_CHECK(H5Tinsert(hDataType, "uncertainty", sizeof(float),
1619 : H5T_IEEE_F32LE)) >= 0);
1620 :
1621 16 : hsize_t chunk_size[] = {static_cast<hsize_t>(nBlockYSize),
1622 16 : static_cast<hsize_t>(nBlockXSize)};
1623 :
1624 16 : const float afFillValue[] = {fNoDataValue, fNoDataValue};
1625 32 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1626 16 : bRet = bRet && hParams &&
1627 16 : H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) >= 0 &&
1628 16 : H5_CHECK(H5Pset_fill_value(hParams, hDataType, afFillValue)) >= 0 &&
1629 48 : H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0 &&
1630 16 : H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) >= 0;
1631 :
1632 16 : if (bRet && bDeflate)
1633 : {
1634 15 : bRet = H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) >= 0;
1635 : }
1636 :
1637 32 : GH5_HIDDatasetHolder hDatasetID;
1638 16 : if (bRet)
1639 : {
1640 16 : hDatasetID.reset(H5_CHECK(H5Dcreate(m_valuesGroup, "values", hDataType,
1641 : hDataSpace, hParams)));
1642 16 : bRet = hDatasetID;
1643 : }
1644 :
1645 32 : GH5_HIDSpaceHolder hFileSpace;
1646 16 : if (bRet)
1647 : {
1648 16 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
1649 16 : bRet = hFileSpace;
1650 : }
1651 :
1652 16 : const int nYBlocks = static_cast<int>(DIV_ROUND_UP(nYSize, nBlockYSize));
1653 16 : const int nXBlocks = static_cast<int>(DIV_ROUND_UP(nXSize, nBlockXSize));
1654 16 : std::vector<float> afValues(static_cast<size_t>(nBlockYSize) * nBlockXSize *
1655 16 : nComponents);
1656 16 : const bool bReverseY = m_gt.yscale < 0;
1657 :
1658 16 : float fMinDepth = std::numeric_limits<float>::infinity();
1659 16 : float fMaxDepth = -std::numeric_limits<float>::infinity();
1660 16 : float fMinUncertainty = std::numeric_limits<float>::infinity();
1661 16 : float fMaxUncertainty = -std::numeric_limits<float>::infinity();
1662 :
1663 16 : int bHasNoDataBand1 = FALSE;
1664 : const char *pszFirstBandDesc =
1665 16 : m_poSrcDS->GetRasterBand(1)->GetDescription();
1666 16 : const float fMulFactor =
1667 16 : EQUAL(pszFirstBandDesc, "elevation") ? -1.0f : 1.0f;
1668 16 : if (fMulFactor < 0.0f)
1669 : {
1670 1 : CPLError(CE_Warning, CPLE_AppDefined,
1671 : "Automatically convert from elevation to depth by negating "
1672 : "elevation values");
1673 : }
1674 : const double dfSrcNoDataBand1 =
1675 16 : m_poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoDataBand1);
1676 16 : const float fSrcNoDataBand1 = static_cast<float>(dfSrcNoDataBand1);
1677 16 : int bHasNoDataBand2 = FALSE;
1678 : const double dfSrcNoDataBand2 =
1679 : nComponents == 2
1680 16 : ? m_poSrcDS->GetRasterBand(2)->GetNoDataValue(&bHasNoDataBand2)
1681 16 : : 0.0;
1682 16 : const float fSrcNoDataBand2 = static_cast<float>(dfSrcNoDataBand2);
1683 :
1684 43 : for (int iY = 0; iY < nYBlocks && bRet; iY++)
1685 : {
1686 : const int nSrcYOff = bReverseY
1687 27 : ? std::max(0, nYSize - (iY + 1) * nBlockYSize)
1688 19 : : iY * nBlockYSize;
1689 27 : const int nReqCountY = std::min(nBlockYSize, nYSize - iY * nBlockYSize);
1690 186 : for (int iX = 0; iX < nXBlocks && bRet; iX++)
1691 : {
1692 : const int nReqCountX =
1693 159 : std::min(nBlockXSize, nXSize - iX * nBlockXSize);
1694 :
1695 159 : bRet =
1696 159 : m_poSrcDS->RasterIO(
1697 159 : GF_Read, iX * nBlockXSize, nSrcYOff, nReqCountX, nReqCountY,
1698 8 : bReverseY ? afValues.data() +
1699 8 : (nReqCountY - 1) * nReqCountX * nComponents
1700 151 : : afValues.data(),
1701 : nReqCountX, nReqCountY, GDT_Float32, nComponents, nullptr,
1702 159 : static_cast<int>(sizeof(float)) * nComponents,
1703 : bReverseY ? -static_cast<GPtrDiff_t>(sizeof(float)) *
1704 8 : nComponents * nReqCountX
1705 : : 0,
1706 : sizeof(float), nullptr) == CE_None;
1707 :
1708 159 : if (bRet)
1709 : {
1710 1441810 : for (int i = 0; i < nReqCountY * nReqCountX; i++)
1711 : {
1712 : {
1713 1441650 : float fVal = afValues[i * nComponents];
1714 2883290 : if ((bHasNoDataBand1 && fVal == fSrcNoDataBand1) ||
1715 1441650 : std::isnan(fVal))
1716 : {
1717 8 : afValues[i * nComponents] = fNoDataValue;
1718 : }
1719 : else
1720 : {
1721 1441640 : fVal *= fMulFactor;
1722 1441640 : afValues[i * nComponents] = fVal;
1723 1441640 : fMinDepth = std::min(fMinDepth, fVal);
1724 1441640 : fMaxDepth = std::max(fMaxDepth, fVal);
1725 : }
1726 : }
1727 1441650 : if (nComponents == 2)
1728 : {
1729 1440010 : const float fVal = afValues[i * nComponents + 1];
1730 2880020 : if ((bHasNoDataBand2 && fVal == fSrcNoDataBand2) ||
1731 1440010 : std::isnan(fVal))
1732 : {
1733 1 : afValues[i * nComponents + 1] = fNoDataValue;
1734 : }
1735 : else
1736 : {
1737 1440010 : fMinUncertainty = std::min(fMinUncertainty, fVal);
1738 1440010 : fMaxUncertainty = std::max(fMaxUncertainty, fVal);
1739 : }
1740 : }
1741 : }
1742 : }
1743 :
1744 : H5OFFSET_TYPE offset[] = {
1745 159 : static_cast<H5OFFSET_TYPE>(iY) *
1746 159 : static_cast<H5OFFSET_TYPE>(nBlockYSize),
1747 159 : static_cast<H5OFFSET_TYPE>(iX) *
1748 159 : static_cast<H5OFFSET_TYPE>(nBlockXSize)};
1749 159 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
1750 159 : static_cast<hsize_t>(nReqCountX)};
1751 : GH5_HIDSpaceHolder hMemSpace(
1752 159 : H5_CHECK(H5Screate_simple(2, count, nullptr)));
1753 159 : bRet =
1754 159 : bRet &&
1755 159 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1756 159 : nullptr, count, nullptr)) >= 0 &&
1757 159 : hMemSpace &&
1758 159 : H5_CHECK(H5Dwrite(hDatasetID, hDataType, hMemSpace, hFileSpace,
1759 318 : H5P_DEFAULT, afValues.data())) >= 0 &&
1760 159 : pfnProgress((static_cast<double>(iY) * nXBlocks + iX + 1) /
1761 159 : (static_cast<double>(nXBlocks) * nYBlocks),
1762 : "", pProgressData) != 0;
1763 : }
1764 : }
1765 :
1766 16 : if (fMinDepth > fMaxDepth)
1767 : {
1768 0 : fMinDepth = fMaxDepth = fNoDataValue;
1769 : }
1770 16 : else if (!(fMinDepth >= MIN_DEPTH_VALUE && fMaxDepth <= MAX_DEPTH_VALUE))
1771 : {
1772 2 : CPLError(CE_Warning, CPLE_AppDefined,
1773 : "Range of depth in the dataset is [%f, %f] whereas the "
1774 : "allowed range is [%d, %d]",
1775 : fMinDepth, fMaxDepth, MIN_DEPTH_VALUE, MAX_DEPTH_VALUE);
1776 : }
1777 :
1778 16 : if (fMinUncertainty > fMaxUncertainty)
1779 : {
1780 13 : fMinUncertainty = fMaxUncertainty = fNoDataValue;
1781 : }
1782 3 : else if (fMinUncertainty < 0)
1783 : {
1784 1 : CPLError(CE_Warning, CPLE_AppDefined,
1785 : "Negative uncertainty value found, which is not allowed");
1786 : }
1787 :
1788 16 : return bRet &&
1789 16 : WriteFloat32Value(m_valuesGroup, "minimumDepth", fMinDepth) &&
1790 16 : WriteFloat32Value(m_valuesGroup, "maximumDepth", fMaxDepth) &&
1791 16 : WriteFloat32Value(m_valuesGroup, "minimumUncertainty",
1792 32 : fMinUncertainty) &&
1793 16 : WriteFloat32Value(m_valuesGroup, "maximumUncertainty",
1794 32 : fMaxUncertainty);
1795 : }
1796 :
1797 : /************************************************************************/
1798 : /* S102Creator::CopyQualityValues() */
1799 : /************************************************************************/
1800 :
1801 3 : bool S102Creator::CopyQualityValues(GDALDataset *poQualityDS,
1802 : const std::set<int> &oSetRATId,
1803 : GDALProgressFunc pfnProgress,
1804 : void *pProgressData)
1805 : {
1806 3 : CPLAssert(m_valuesGroup.get() >= 0);
1807 :
1808 3 : const int nYSize = poQualityDS->GetRasterYSize();
1809 3 : const int nXSize = poQualityDS->GetRasterXSize();
1810 :
1811 3 : hsize_t dims[] = {static_cast<hsize_t>(nYSize),
1812 3 : static_cast<hsize_t>(nXSize)};
1813 :
1814 6 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(2, dims, nullptr)));
1815 3 : bool bRet = hDataSpace;
1816 :
1817 : const bool bDeflate =
1818 3 : EQUAL(m_aosOptions.FetchNameValueDef("COMPRESS", "DEFLATE"), "DEFLATE");
1819 : const int nCompressionLevel =
1820 3 : atoi(m_aosOptions.FetchNameValueDef("ZLEVEL", "6"));
1821 : const int nBlockSize =
1822 3 : std::min(4096, std::max(100, atoi(m_aosOptions.FetchNameValueDef(
1823 3 : "BLOCK_SIZE", "100"))));
1824 3 : const int nBlockXSize = std::min(nXSize, nBlockSize);
1825 3 : const int nBlockYSize = std::min(nYSize, nBlockSize);
1826 3 : constexpr uint32_t nNoDataValue = 0;
1827 :
1828 3 : hsize_t chunk_size[] = {static_cast<hsize_t>(nBlockYSize),
1829 3 : static_cast<hsize_t>(nBlockXSize)};
1830 :
1831 6 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1832 3 : bRet = bRet && hParams &&
1833 3 : H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) >= 0 &&
1834 3 : H5_CHECK(H5Pset_fill_value(hParams, H5T_STD_U32LE, &nNoDataValue)) >=
1835 3 : 0 &&
1836 9 : H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0 &&
1837 3 : H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) >= 0;
1838 :
1839 3 : if (bRet && bDeflate)
1840 : {
1841 3 : bRet = H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) >= 0;
1842 : }
1843 :
1844 6 : GH5_HIDDatasetHolder hDatasetID;
1845 3 : if (bRet)
1846 : {
1847 3 : hDatasetID.reset(H5_CHECK(H5Dcreate(
1848 : m_valuesGroup, "values", H5T_STD_U32LE, hDataSpace, hParams)));
1849 3 : bRet = hDatasetID;
1850 : }
1851 :
1852 6 : GH5_HIDSpaceHolder hFileSpace(H5_CHECK(H5Dget_space(hDatasetID)));
1853 3 : bRet = bRet && hFileSpace;
1854 :
1855 3 : const int nYBlocks = static_cast<int>(DIV_ROUND_UP(nYSize, nBlockYSize));
1856 3 : const int nXBlocks = static_cast<int>(DIV_ROUND_UP(nXSize, nBlockXSize));
1857 3 : std::vector<uint32_t> anValues(static_cast<size_t>(nBlockYSize) *
1858 6 : nBlockXSize);
1859 3 : const bool bReverseY = m_gt.yscale < 0;
1860 :
1861 3 : int bHasSrcNoData = FALSE;
1862 : const double dfSrcNoData =
1863 3 : poQualityDS->GetRasterBand(1)->GetNoDataValue(&bHasSrcNoData);
1864 3 : const uint32_t nSrcNoData = static_cast<uint32_t>(dfSrcNoData);
1865 :
1866 3 : std::set<int> oSetRATIdCopy(oSetRATId);
1867 6 : for (int iY = 0; iY < nYBlocks && bRet; iY++)
1868 : {
1869 : const int nSrcYOff = bReverseY
1870 3 : ? std::max(0, nYSize - (iY + 1) * nBlockYSize)
1871 3 : : iY * nBlockYSize;
1872 3 : const int nReqCountY = std::min(nBlockYSize, nYSize - iY * nBlockYSize);
1873 6 : for (int iX = 0; iX < nXBlocks && bRet; iX++)
1874 : {
1875 : const int nReqCountX =
1876 3 : std::min(nBlockXSize, nXSize - iX * nBlockXSize);
1877 :
1878 3 : bRet =
1879 3 : poQualityDS->GetRasterBand(1)->RasterIO(
1880 3 : GF_Read, iX * nBlockXSize, nSrcYOff, nReqCountX, nReqCountY,
1881 0 : bReverseY ? anValues.data() + (nReqCountY - 1) * nReqCountX
1882 3 : : anValues.data(),
1883 : nReqCountX, nReqCountY, GDT_UInt32, 0,
1884 : bReverseY ? -static_cast<GPtrDiff_t>(sizeof(uint32_t)) *
1885 0 : nReqCountX
1886 : : 0,
1887 : nullptr) == CE_None;
1888 :
1889 3 : if (bRet)
1890 : {
1891 15 : for (int i = 0; i < nReqCountY * nReqCountX; i++)
1892 : {
1893 12 : if (bHasSrcNoData && anValues[i] == nSrcNoData)
1894 : {
1895 0 : anValues[i] = nNoDataValue;
1896 : }
1897 24 : else if (anValues[i] != 0 &&
1898 12 : !cpl::contains(oSetRATIdCopy, anValues[i]))
1899 : {
1900 1 : CPLError(
1901 : CE_Warning, CPLE_AppDefined,
1902 : "Quality grid contains nodes with id %u, but there "
1903 : "is no such entry in the feature attribute table",
1904 1 : anValues[i]);
1905 1 : oSetRATIdCopy.insert(anValues[i]);
1906 : }
1907 : }
1908 : }
1909 :
1910 : H5OFFSET_TYPE offset[] = {
1911 3 : static_cast<H5OFFSET_TYPE>(iY) *
1912 3 : static_cast<H5OFFSET_TYPE>(nBlockYSize),
1913 3 : static_cast<H5OFFSET_TYPE>(iX) *
1914 3 : static_cast<H5OFFSET_TYPE>(nBlockXSize)};
1915 3 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
1916 3 : static_cast<hsize_t>(nReqCountX)};
1917 3 : GH5_HIDSpaceHolder hMemSpace(H5Screate_simple(2, count, nullptr));
1918 3 : bRet =
1919 3 : bRet && hMemSpace &&
1920 3 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1921 3 : nullptr, count, nullptr)) >= 0 &&
1922 3 : H5_CHECK(H5Dwrite(hDatasetID, H5T_NATIVE_UINT, hMemSpace,
1923 : hFileSpace, H5P_DEFAULT, anValues.data())) >=
1924 6 : 0 &&
1925 3 : pfnProgress((static_cast<double>(iY) * nXBlocks + iX + 1) /
1926 3 : (static_cast<double>(nXBlocks) * nYBlocks),
1927 : "", pProgressData) != 0;
1928 : }
1929 : }
1930 :
1931 6 : return bRet;
1932 : }
1933 :
1934 : /************************************************************************/
1935 : /* S102Dataset::CreateCopy() */
1936 : /************************************************************************/
1937 :
1938 : /* static */
1939 56 : GDALDataset *S102Dataset::CreateCopy(const char *pszFilename,
1940 : GDALDataset *poSrcDS, int /* bStrict*/,
1941 : CSLConstList papszOptions,
1942 : GDALProgressFunc pfnProgress,
1943 : void *pProgressData)
1944 : {
1945 112 : S102Creator creator(pszFilename, poSrcDS, papszOptions);
1946 56 : if (!creator.Create(pfnProgress, pProgressData))
1947 40 : return nullptr;
1948 :
1949 : VSIStatBufL sStatBuf;
1950 32 : if (VSIStatL(pszFilename, &sStatBuf) == 0 &&
1951 16 : sStatBuf.st_size > 10 * 1024 * 1024)
1952 : {
1953 1 : CPLError(CE_Warning, CPLE_AppDefined,
1954 : "%s file size exceeds 10 MB, which is the upper limit "
1955 : "suggested for wireless transmission to marine vessels",
1956 : pszFilename);
1957 : }
1958 :
1959 32 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
1960 16 : return Open(&oOpenInfo);
1961 : }
1962 :
1963 : /************************************************************************/
1964 : /* S102DatasetDriverUnload() */
1965 : /************************************************************************/
1966 :
1967 9 : static void S102DatasetDriverUnload(GDALDriver *)
1968 : {
1969 9 : HDF5UnloadFileDriver();
1970 9 : }
1971 :
1972 : /************************************************************************/
1973 : /* GDALRegister_S102() */
1974 : /************************************************************************/
1975 14 : void GDALRegister_S102()
1976 :
1977 : {
1978 14 : if (!GDAL_CHECK_VERSION("S102"))
1979 0 : return;
1980 :
1981 14 : if (GDALGetDriverByName(S102_DRIVER_NAME) != nullptr)
1982 0 : return;
1983 :
1984 14 : GDALDriver *poDriver = new GDALDriver();
1985 :
1986 14 : S102DriverSetCommonMetadata(poDriver);
1987 14 : poDriver->pfnOpen = S102Dataset::Open;
1988 14 : poDriver->pfnCreateCopy = S102Dataset::CreateCopy;
1989 14 : poDriver->pfnUnloadDriver = S102DatasetDriverUnload;
1990 :
1991 14 : GetGDALDriverManager()->RegisterDriver(poDriver);
1992 : }
|
