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