Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Hierarchical Data Format Release 5 (HDF5)
4 : * Purpose: Read/create S111 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 "hdf5dataset.h"
15 : #include "hdf5drivercore.h"
16 : #include "gh5_convenience.h"
17 : #include "rat.h"
18 : #include "s100.h"
19 :
20 : #include "gdal_frmts.h"
21 : #include "gdal_priv.h"
22 : #include "gdal_proxy.h"
23 : #include "gdal_rat.h"
24 :
25 : #include "cpl_time.h"
26 :
27 : #include <algorithm>
28 : #include <array>
29 : #include <charconv>
30 : #include <cmath>
31 : #include <limits>
32 : #include <map>
33 : #include <variant>
34 :
35 : /************************************************************************/
36 : /* S111Dataset */
37 : /************************************************************************/
38 :
39 162 : class S111Dataset final : public S100BaseDataset
40 : {
41 : public:
42 81 : explicit S111Dataset(const std::string &osFilename)
43 81 : : S100BaseDataset(osFilename)
44 : {
45 81 : }
46 :
47 : ~S111Dataset() override;
48 :
49 : static GDALDataset *Open(GDALOpenInfo *);
50 : static GDALDataset *CreateCopy(const char *pszFilename,
51 : GDALDataset *poSrcDS, int bStrict,
52 : char **papszOptions,
53 : GDALProgressFunc pfnProgress,
54 : void *pProgressData);
55 : };
56 :
57 : S111Dataset::~S111Dataset() = default;
58 :
59 : /************************************************************************/
60 : /* S111RasterBand */
61 : /************************************************************************/
62 :
63 : class S111RasterBand final : public GDALProxyRasterBand
64 : {
65 : friend class S111Dataset;
66 : std::unique_ptr<GDALDataset> m_poDS{};
67 : GDALRasterBand *m_poUnderlyingBand = nullptr;
68 : std::string m_osUnitType{};
69 : std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
70 :
71 : CPL_DISALLOW_COPY_ASSIGN(S111RasterBand)
72 :
73 : public:
74 72 : explicit S111RasterBand(std::unique_ptr<GDALDataset> &&poDSIn)
75 144 : : m_poDS(std::move(poDSIn)),
76 72 : m_poUnderlyingBand(m_poDS->GetRasterBand(1))
77 : {
78 72 : eDataType = m_poUnderlyingBand->GetRasterDataType();
79 72 : m_poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
80 72 : }
81 :
82 : GDALRasterBand *
83 : RefUnderlyingRasterBand(bool /*bForceOpen*/ = true) const override;
84 :
85 4 : const char *GetUnitType() override
86 : {
87 4 : return m_osUnitType.c_str();
88 : }
89 :
90 1 : GDALRasterAttributeTable *GetDefaultRAT() override
91 : {
92 1 : return m_poRAT.get();
93 : }
94 :
95 0 : char **GetMetadata(const char *pszDomain) override
96 : {
97 : // Short-circuit GDALProxyRasterBand...
98 0 : return GDALRasterBand::GetMetadata(pszDomain);
99 : }
100 : };
101 :
102 : GDALRasterBand *
103 50 : S111RasterBand::RefUnderlyingRasterBand(bool /*bForceOpen*/) const
104 : {
105 50 : return m_poUnderlyingBand;
106 : }
107 :
108 : /************************************************************************/
109 : /* Open() */
110 : /************************************************************************/
111 :
112 82 : GDALDataset *S111Dataset::Open(GDALOpenInfo *poOpenInfo)
113 :
114 : {
115 : // Confirm that this appears to be a S111 file.
116 82 : if (!S111DatasetIdentify(poOpenInfo))
117 0 : return nullptr;
118 :
119 : HDF5_GLOBAL_LOCK();
120 :
121 82 : if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
122 : {
123 1 : return HDF5Dataset::OpenMultiDim(poOpenInfo);
124 : }
125 :
126 : // Confirm the requested access is supported.
127 81 : if (poOpenInfo->eAccess == GA_Update)
128 : {
129 0 : ReportUpdateNotSupportedByDriver("S111");
130 0 : return nullptr;
131 : }
132 :
133 162 : std::string osFilename(poOpenInfo->pszFilename);
134 162 : std::string osFeatureInstance = "SurfaceCurrent.01";
135 162 : std::string osGroup;
136 81 : if (STARTS_WITH(poOpenInfo->pszFilename, "S111:"))
137 : {
138 : const CPLStringList aosTokens(
139 37 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
140 37 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
141 :
142 37 : if (aosTokens.size() == 2)
143 : {
144 1 : osFilename = aosTokens[1];
145 : }
146 36 : else if (aosTokens.size() == 3)
147 : {
148 28 : osFilename = aosTokens[1];
149 28 : osGroup = aosTokens[2];
150 : }
151 8 : else if (aosTokens.size() == 4)
152 : {
153 8 : osFilename = aosTokens[1];
154 8 : osFeatureInstance = aosTokens[2];
155 8 : osGroup = aosTokens[3];
156 : }
157 : else
158 : {
159 0 : return nullptr;
160 : }
161 : }
162 :
163 162 : auto poDS = std::make_unique<S111Dataset>(osFilename);
164 81 : if (!poDS->Init())
165 0 : return nullptr;
166 :
167 81 : const auto &poRootGroup = poDS->m_poRootGroup;
168 :
169 : // Read additional metadata
170 162 : for (const char *pszAttrName : {"depthTypeIndex"})
171 : {
172 243 : auto poAttr = poRootGroup->GetAttribute(pszAttrName);
173 81 : if (poAttr)
174 : {
175 75 : const char *pszVal = poAttr->ReadAsString();
176 75 : if (pszVal)
177 : {
178 75 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
179 : }
180 : }
181 : }
182 :
183 81 : if (const auto poDepthTypeIndex =
184 243 : poRootGroup->GetAttribute("depthTypeIndex"))
185 : {
186 75 : if (poDepthTypeIndex->GetDataType().GetClass() == GEDTC_NUMERIC)
187 : {
188 75 : const int nDepthType = poDepthTypeIndex->ReadAsInt();
189 75 : if (nDepthType == 1)
190 : {
191 16 : poDS->GDALDataset::SetMetadataItem("DEPTH_TYPE_INDEX_NAME",
192 : "heightOrDepth");
193 16 : poDS->GDALDataset::SetMetadataItem(
194 : "DEPTH_TYPE_INDEX_DEFINITION", "Height or depth");
195 : }
196 59 : else if (nDepthType == 2)
197 : {
198 59 : poDS->GDALDataset::SetMetadataItem("DEPTH_TYPE_INDEX_NAME",
199 : "layerAverage");
200 59 : poDS->GDALDataset::SetMetadataItem(
201 : "DEPTH_TYPE_INDEX_DEFINITION", "Layer average");
202 : }
203 150 : poDS->GDALDataset::SetMetadataItem(
204 150 : "depthTypeIndex", std::to_string(nDepthType).c_str());
205 : }
206 : }
207 :
208 243 : auto poVerticalCS = poRootGroup->GetAttribute("verticalCS");
209 81 : if (poVerticalCS && poVerticalCS->GetDataType().GetClass() == GEDTC_NUMERIC)
210 : {
211 18 : const int nVerticalCS = poVerticalCS->ReadAsInt();
212 18 : if (nVerticalCS == 6498)
213 18 : poDS->GDALDataset::SetMetadataItem(
214 : "VERTICAL_CS_DEFINITION", "depth, meters, orientation down");
215 0 : else if (nVerticalCS == 6499)
216 0 : poDS->GDALDataset::SetMetadataItem(
217 : "VERTICAL_CS_DEFINITION", "height, meters, orientation up");
218 :
219 36 : poDS->GDALDataset::SetMetadataItem("verticalCS",
220 36 : std::to_string(nVerticalCS).c_str());
221 : }
222 :
223 243 : auto poSurfaceCurrent = poRootGroup->OpenGroup("SurfaceCurrent");
224 81 : if (!poSurfaceCurrent)
225 : {
226 7 : CPLError(CE_Failure, CPLE_AppDefined,
227 : "Cannot find /SurfaceCurrent group");
228 7 : return nullptr;
229 : }
230 :
231 : auto poDataCodingFormat =
232 222 : poSurfaceCurrent->GetAttribute("dataCodingFormat");
233 148 : if (!poDataCodingFormat ||
234 74 : poDataCodingFormat->GetDataType().GetClass() != GEDTC_NUMERIC)
235 : {
236 0 : CPLError(CE_Failure, CPLE_AppDefined,
237 : "Cannot find /SurfaceCurrent/dataCodingFormat attribute");
238 0 : return nullptr;
239 : }
240 74 : const int nDataCodingFormat = poDataCodingFormat->ReadAsInt();
241 74 : if (nDataCodingFormat != 2)
242 : {
243 0 : CPLError(CE_Failure, CPLE_NotSupported,
244 : "dataCodingFormat=%d is not supported by the S111 driver",
245 : nDataCodingFormat);
246 0 : return nullptr;
247 : }
248 :
249 : // Read additional metadata
250 370 : for (const char *pszAttrName :
251 : {"methodCurrentsProduct", "minDatasetCurrentSpeed",
252 444 : "maxDatasetCurrentSpeed", "commonPointRule", "interpolationType"})
253 : {
254 1110 : auto poAttr = poSurfaceCurrent->GetAttribute(pszAttrName);
255 370 : if (poAttr)
256 : {
257 284 : const char *pszVal = poAttr->ReadAsString();
258 284 : if (pszVal)
259 : {
260 284 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
261 : }
262 : }
263 : }
264 :
265 74 : if (auto poCommonPointRule =
266 222 : poSurfaceCurrent->GetAttribute("commonPointRule"))
267 : {
268 70 : poDS->SetMetadataForCommonPointRule(poCommonPointRule.get());
269 : }
270 :
271 74 : if (auto poInterpolationType =
272 222 : poSurfaceCurrent->GetAttribute("interpolationType"))
273 : {
274 70 : poDS->SetMetadataForInterpolationType(poInterpolationType.get());
275 : }
276 :
277 74 : int nNumInstances = 1;
278 74 : if (osGroup.empty())
279 : {
280 114 : auto poNumInstances = poSurfaceCurrent->GetAttribute("numInstances");
281 75 : if (poNumInstances &&
282 75 : poNumInstances->GetDataType().GetClass() == GEDTC_NUMERIC)
283 : {
284 37 : nNumInstances = poNumInstances->ReadAsInt();
285 : }
286 : }
287 74 : if (nNumInstances != 1)
288 : {
289 8 : CPLStringList aosSubDSList;
290 4 : int iSubDS = 0;
291 8 : for (const std::string &featureInstanceName :
292 20 : poSurfaceCurrent->GetGroupNames())
293 : {
294 : auto poFeatureInstance =
295 16 : poSurfaceCurrent->OpenGroup(featureInstanceName);
296 8 : if (poFeatureInstance)
297 : {
298 16 : GDALMajorObject mo;
299 : // Read first vertical datum from root group and let the
300 : // coverage override it.
301 8 : S100ReadVerticalDatum(&mo, poRootGroup.get());
302 8 : S100ReadVerticalDatum(&mo, poFeatureInstance.get());
303 :
304 16 : const auto aosGroupNames = poFeatureInstance->GetGroupNames();
305 16 : for (const auto &osSubGroup : aosGroupNames)
306 : {
307 8 : if (auto poSubGroup =
308 16 : poFeatureInstance->OpenGroup(osSubGroup))
309 : {
310 8 : ++iSubDS;
311 : aosSubDSList.SetNameValue(
312 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
313 : CPLSPrintf("S111:\"%s\":%s:%s", osFilename.c_str(),
314 : featureInstanceName.c_str(),
315 8 : osSubGroup.c_str()));
316 :
317 16 : std::string verticalDatum;
318 : const char *pszValue =
319 8 : mo.GetMetadataItem(S100_VERTICAL_DATUM_NAME);
320 8 : if (pszValue)
321 : {
322 2 : verticalDatum = ", vertical datum ";
323 2 : verticalDatum += pszValue;
324 : pszValue =
325 2 : mo.GetMetadataItem(S100_VERTICAL_DATUM_ABBREV);
326 2 : if (pszValue)
327 : {
328 2 : verticalDatum += " (";
329 2 : verticalDatum += pszValue;
330 2 : verticalDatum += ')';
331 : }
332 : }
333 :
334 16 : std::string osSubDSDesc;
335 : const auto poTimePoint =
336 24 : poSubGroup->GetAttribute("timePoint");
337 8 : if (poTimePoint)
338 : {
339 8 : const char *pszVal = poTimePoint->ReadAsString();
340 8 : if (pszVal)
341 : {
342 8 : osSubDSDesc = "Values for feature instance ";
343 8 : osSubDSDesc += featureInstanceName;
344 8 : osSubDSDesc += verticalDatum;
345 8 : osSubDSDesc += " at timestamp ";
346 8 : osSubDSDesc += pszVal;
347 : }
348 : }
349 8 : if (osSubDSDesc.empty())
350 : {
351 0 : osSubDSDesc = "Values for feature instance ";
352 0 : osSubDSDesc += featureInstanceName;
353 0 : osSubDSDesc += verticalDatum;
354 0 : osSubDSDesc += " and group ";
355 0 : osSubDSDesc += osSubGroup;
356 : }
357 :
358 : aosSubDSList.SetNameValue(
359 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
360 8 : osSubDSDesc.c_str());
361 : }
362 : }
363 : }
364 : }
365 :
366 4 : poDS->GDALDataset::SetMetadata(aosSubDSList.List(), "SUBDATASETS");
367 :
368 : // Setup/check for pam .aux.xml.
369 4 : poDS->SetDescription(osFilename.c_str());
370 4 : poDS->TryLoadXML();
371 :
372 : // Setup overviews.
373 4 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
374 :
375 4 : return poDS.release();
376 : }
377 :
378 140 : auto poFeatureInstance = poSurfaceCurrent->OpenGroup(osFeatureInstance);
379 70 : if (!poFeatureInstance)
380 : {
381 0 : CPLError(CE_Failure, CPLE_AppDefined,
382 : "Cannot find /SurfaceCurrent/%s group",
383 : osFeatureInstance.c_str());
384 0 : return nullptr;
385 : }
386 :
387 : // Read additional metadata
388 350 : for (const char *pszAttrName :
389 : {"timeRecordInterval", "dateTimeOfFirstRecord", "dateTimeOfLastRecord",
390 420 : "numberOfTimes", "dataDynamicity"})
391 : {
392 1050 : auto poAttr = poFeatureInstance->GetAttribute(pszAttrName);
393 350 : if (poAttr)
394 : {
395 292 : const char *pszVal = poAttr->ReadAsString();
396 292 : if (pszVal)
397 : {
398 292 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
399 : }
400 : }
401 : }
402 :
403 70 : if (auto poDataDynamicity =
404 210 : poFeatureInstance->GetAttribute("dataDynamicity"))
405 : {
406 64 : poDS->SetMetadataForDataDynamicity(poDataDynamicity.get());
407 : }
408 :
409 : // Read optional uncertainty array
410 198 : if (auto poUncertainty = poFeatureInstance->OpenMDArray("uncertainty"))
411 : {
412 58 : auto &apoDims = poUncertainty->GetDimensions();
413 116 : if (poUncertainty->GetDataType().GetClass() == GEDTC_COMPOUND &&
414 116 : apoDims.size() == 1 && apoDims[0]->GetSize() == 2)
415 : {
416 : const auto &oComponents =
417 58 : poUncertainty->GetDataType().GetComponents();
418 116 : if (oComponents.size() == 2 &&
419 116 : oComponents[0]->GetName() == "name" &&
420 116 : oComponents[0]->GetType().GetClass() == GEDTC_STRING &&
421 175 : oComponents[1]->GetName() == "value" &&
422 58 : (oComponents[1]->GetType().GetNumericDataType() ==
423 1 : GDT_Float32 ||
424 1 : oComponents[1]->GetType().GetNumericDataType() == GDT_Float64))
425 : {
426 174 : auto poName = poUncertainty->GetView("[\"name\"]");
427 174 : auto poValue = poUncertainty->GetView("[\"value\"]");
428 58 : if (poName && poValue)
429 : {
430 58 : char *apszStr[2] = {nullptr, nullptr};
431 58 : double adfVals[2] = {0, 0};
432 58 : GUInt64 arrayStartIdx[] = {0};
433 58 : size_t count[] = {2};
434 58 : GInt64 arrayStep[] = {1};
435 58 : GPtrDiff_t bufferStride[] = {1};
436 58 : if (poName->Read(arrayStartIdx, count, arrayStep,
437 58 : bufferStride, oComponents[0]->GetType(),
438 174 : apszStr) &&
439 116 : poValue->Read(
440 : arrayStartIdx, count, arrayStep, bufferStride,
441 116 : GDALExtendedDataType::Create(GDT_Float64), adfVals))
442 : {
443 174 : for (int i = 0; i < 2; ++i)
444 : {
445 232 : std::string osName = apszStr[i];
446 116 : if (osName[0] >= 'a' && osName[0] <= 'z')
447 116 : osName[0] = osName[0] - 'a' + 'A';
448 116 : osName = "uncertainty" + osName;
449 116 : poDS->GDALDataset::SetMetadataItem(
450 : osName.c_str(), CPLSPrintf("%f", adfVals[i]));
451 : }
452 : }
453 58 : VSIFree(apszStr[0]);
454 58 : VSIFree(apszStr[1]);
455 : }
456 : }
457 : }
458 : }
459 :
460 140 : if (auto poStartSequence = poFeatureInstance->GetAttribute("startSequence"))
461 : {
462 70 : const char *pszStartSequence = poStartSequence->ReadAsString();
463 70 : if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
464 : {
465 0 : CPLError(CE_Failure, CPLE_AppDefined,
466 : "startSequence (=%s) != 0,0 is not supported",
467 : pszStartSequence);
468 0 : return nullptr;
469 : }
470 : }
471 :
472 70 : if (!S100GetNumPointsLongitudinalLatitudinal(
473 70 : poFeatureInstance.get(), poDS->nRasterXSize, poDS->nRasterYSize))
474 : {
475 0 : return nullptr;
476 : }
477 :
478 : // Potentially override vertical datum
479 70 : S100ReadVerticalDatum(poDS.get(), poFeatureInstance.get());
480 :
481 70 : const bool bNorthUp = CPLTestBool(
482 70 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
483 :
484 : // Compute geotransform
485 140 : poDS->m_bHasGT =
486 70 : S100GetGeoTransform(poFeatureInstance.get(), poDS->m_gt, bNorthUp);
487 :
488 70 : if (osGroup.empty())
489 : {
490 68 : const auto aosGroupNames = poFeatureInstance->GetGroupNames();
491 34 : int iSubDS = 1;
492 72 : for (const auto &osSubGroup : aosGroupNames)
493 : {
494 76 : if (auto poSubGroup = poFeatureInstance->OpenGroup(osSubGroup))
495 : {
496 38 : poDS->GDALDataset::SetMetadataItem(
497 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
498 : CPLSPrintf("S111:\"%s\":%s", osFilename.c_str(),
499 : osSubGroup.c_str()),
500 : "SUBDATASETS");
501 76 : std::string osSubDSDesc = "Values for group ";
502 38 : osSubDSDesc += osSubGroup;
503 76 : const auto poTimePoint = poSubGroup->GetAttribute("timePoint");
504 38 : if (poTimePoint)
505 : {
506 38 : const char *pszVal = poTimePoint->ReadAsString();
507 38 : if (pszVal)
508 : {
509 38 : osSubDSDesc = "Values at timestamp ";
510 38 : osSubDSDesc += pszVal;
511 : }
512 : }
513 38 : poDS->GDALDataset::SetMetadataItem(
514 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
515 : osSubDSDesc.c_str(), "SUBDATASETS");
516 38 : ++iSubDS;
517 : }
518 : }
519 : }
520 : else
521 : {
522 36 : auto poGroup = poFeatureInstance->OpenGroup(osGroup);
523 36 : if (!poGroup)
524 : {
525 1 : CPLError(CE_Failure, CPLE_AppDefined,
526 : "Cannot find /SurfaceCurrent/%s/%s group",
527 : osFeatureInstance.c_str(), osGroup.c_str());
528 1 : return nullptr;
529 : }
530 :
531 : // Read additional metadata
532 70 : for (const char *pszAttrName : {"timePoint"})
533 : {
534 105 : auto poAttr = poGroup->GetAttribute(pszAttrName);
535 35 : if (poAttr)
536 : {
537 35 : const char *pszVal = poAttr->ReadAsString();
538 35 : if (pszVal)
539 : {
540 35 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
541 : }
542 : }
543 : }
544 :
545 70 : auto poValuesArray = poGroup->OpenMDArray("values");
546 35 : if (!poValuesArray)
547 : {
548 0 : CPLError(CE_Failure, CPLE_AppDefined,
549 : "Cannot find /SurfaceCurrent/%s/%s/values array",
550 : osFeatureInstance.c_str(), osGroup.c_str());
551 0 : return nullptr;
552 : }
553 :
554 35 : if (poValuesArray->GetDimensionCount() != 2)
555 : {
556 0 : CPLError(CE_Failure, CPLE_AppDefined,
557 : "Wrong dimension count for %s",
558 0 : poValuesArray->GetFullName().c_str());
559 0 : return nullptr;
560 : }
561 :
562 35 : const auto &oType = poValuesArray->GetDataType();
563 35 : if (oType.GetClass() != GEDTC_COMPOUND)
564 : {
565 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
566 0 : poValuesArray->GetFullName().c_str());
567 0 : return nullptr;
568 : }
569 :
570 35 : const auto &oComponents = oType.GetComponents();
571 70 : if (!(oComponents.size() >= 2 &&
572 70 : ((oComponents[0]->GetName() == "surfaceCurrentSpeed" &&
573 35 : oComponents[0]->GetType().GetNumericDataType() == GDT_Float32 &&
574 70 : oComponents[1]->GetName() == "surfaceCurrentDirection" &&
575 35 : oComponents[1]->GetType().GetNumericDataType() ==
576 0 : GDT_Float32) ||
577 : // S111US_20170829.0100_W078.N44_F2_loofs_type2.h5 has direction first...
578 0 : (oComponents[0]->GetName() == "surfaceCurrentDirection" &&
579 0 : oComponents[0]->GetType().GetNumericDataType() == GDT_Float32 &&
580 0 : oComponents[1]->GetName() == "surfaceCurrentSpeed" &&
581 0 : oComponents[1]->GetType().GetNumericDataType() ==
582 : GDT_Float32))))
583 : {
584 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
585 0 : poValuesArray->GetFullName().c_str());
586 0 : return nullptr;
587 : }
588 :
589 35 : const auto &apoDims = poValuesArray->GetDimensions();
590 35 : if (apoDims[0]->GetSize() != static_cast<unsigned>(poDS->nRasterYSize))
591 : {
592 0 : CPLError(CE_Failure, CPLE_AppDefined,
593 : "numPointsLatitudinal(=%d) doesn't match first dimension "
594 : "size of %s (=%d)",
595 0 : poDS->nRasterYSize, poValuesArray->GetFullName().c_str(),
596 0 : static_cast<int>(apoDims[0]->GetSize()));
597 0 : return nullptr;
598 : }
599 35 : if (apoDims[1]->GetSize() != static_cast<unsigned>(poDS->nRasterXSize))
600 : {
601 0 : CPLError(CE_Failure, CPLE_AppDefined,
602 : "numPointsLongitudinal(=%d) doesn't match second "
603 : "dimension size of %s (=%d)",
604 0 : poDS->nRasterXSize, poValuesArray->GetFullName().c_str(),
605 0 : static_cast<int>(apoDims[1]->GetSize()));
606 0 : return nullptr;
607 : }
608 :
609 35 : if (bNorthUp)
610 34 : poValuesArray = poValuesArray->GetView("[::-1,...]");
611 :
612 : // Create surfaceCurrentSpeed band
613 : auto poSurfaceCurrentSpeed =
614 105 : poValuesArray->GetView("[\"surfaceCurrentSpeed\"]");
615 : auto poSurfaceCurrentSpeedDS = std::unique_ptr<GDALDataset>(
616 70 : poSurfaceCurrentSpeed->AsClassicDataset(1, 0));
617 : auto poSurfaceCurrentSpeedBand = std::make_unique<S111RasterBand>(
618 70 : std::move(poSurfaceCurrentSpeedDS));
619 35 : poSurfaceCurrentSpeedBand->SetDescription("surfaceCurrentSpeed");
620 35 : poSurfaceCurrentSpeedBand->m_osUnitType = "knots";
621 :
622 : // From S-111 v1.2 table 9.1 (Speed ranges) and 9.2 (Colour schemas)
623 70 : auto poRAT = std::make_unique<GDALDefaultRasterAttributeTable>();
624 35 : poRAT->CreateColumn("speed_band", GFT_Integer, GFU_Generic);
625 35 : poRAT->CreateColumn("min_speed", GFT_Real, GFU_Min);
626 35 : poRAT->CreateColumn("width_band", GFT_Real, GFU_Generic);
627 35 : poRAT->CreateColumn("color", GFT_String, GFU_Generic);
628 35 : poRAT->CreateColumn("red", GFT_Integer, GFU_RedMin);
629 35 : poRAT->CreateColumn("green", GFT_Integer, GFU_GreenMin);
630 35 : poRAT->CreateColumn("blue", GFT_Integer, GFU_BlueMin);
631 :
632 : const struct
633 : {
634 : int nSpeedBand;
635 : double dfMinSpeed;
636 : double dfWidthBand;
637 : const char *pszColor;
638 : int nRed;
639 : int nGreen;
640 : int nBlue;
641 35 : } aoRatValues[] = {
642 : {1, 0.0, 0.5, "purple", 118, 82, 226},
643 : {2, 0.5, 0.5, "dark blue", 72, 152, 211},
644 : {3, 1.0, 1.0, "light blue", 97, 203, 229},
645 : {4, 2.0, 1.0, "dark green", 109, 188, 69},
646 : {5, 3.0, 2.0, "light green", 180, 220, 0},
647 : {6, 5.0, 2.0, "yellow green", 205, 193, 0},
648 : {7, 7.0, 3.0, "orange", 248, 167, 24},
649 : {8, 10.0, 3.0, "pink", 247, 162, 157},
650 : {9, 13.0, 86.0, "red", 255, 30, 30},
651 : };
652 :
653 35 : int iRow = 0;
654 350 : for (const auto &oRecord : aoRatValues)
655 : {
656 315 : int iCol = 0;
657 315 : poRAT->SetValue(iRow, iCol++, oRecord.nSpeedBand);
658 315 : poRAT->SetValue(iRow, iCol++, oRecord.dfMinSpeed);
659 315 : poRAT->SetValue(iRow, iCol++, oRecord.dfWidthBand);
660 315 : poRAT->SetValue(iRow, iCol++, oRecord.pszColor);
661 315 : poRAT->SetValue(iRow, iCol++, oRecord.nRed);
662 315 : poRAT->SetValue(iRow, iCol++, oRecord.nGreen);
663 315 : poRAT->SetValue(iRow, iCol++, oRecord.nBlue);
664 315 : ++iRow;
665 : }
666 :
667 35 : poSurfaceCurrentSpeedBand->m_poRAT = std::move(poRAT);
668 :
669 35 : poDS->SetBand(1, poSurfaceCurrentSpeedBand.release());
670 :
671 : // Create surfaceCurrentDirection band
672 : auto poSurfaceCurrentDirection =
673 105 : poValuesArray->GetView("[\"surfaceCurrentDirection\"]");
674 : auto poSurfaceCurrentDirectionDS = std::unique_ptr<GDALDataset>(
675 70 : poSurfaceCurrentDirection->AsClassicDataset(1, 0));
676 : auto poSurfaceCurrentDirectionBand = std::make_unique<S111RasterBand>(
677 70 : std::move(poSurfaceCurrentDirectionDS));
678 35 : poSurfaceCurrentDirectionBand->SetDescription(
679 35 : "surfaceCurrentDirection");
680 35 : poSurfaceCurrentDirectionBand->m_osUnitType = "degree";
681 35 : poSurfaceCurrentDirectionBand->GDALRasterBand::SetMetadataItem(
682 : "ANGLE_CONVENTION", "From true north, clockwise");
683 35 : poDS->SetBand(2, poSurfaceCurrentDirectionBand.release());
684 :
685 37 : for (size_t i = 2; i < oComponents.size(); ++i)
686 : {
687 3 : if (oComponents[i]->GetName() == "speedUncertainty" &&
688 1 : oComponents[i]->GetType().GetNumericDataType() == GDT_Float32)
689 : {
690 : auto poSubArray =
691 3 : poValuesArray->GetView("[\"speedUncertainty\"]");
692 : auto poSubArrayDS = std::unique_ptr<GDALDataset>(
693 2 : poSubArray->AsClassicDataset(1, 0));
694 : auto poSubArrayBand =
695 2 : std::make_unique<S111RasterBand>(std::move(poSubArrayDS));
696 1 : poSubArrayBand->SetDescription("speedUncertainty");
697 1 : poSubArrayBand->m_osUnitType = "knot";
698 1 : poDS->SetBand(poDS->nBands + 1, poSubArrayBand.release());
699 : }
700 2 : else if (oComponents[i]->GetName() == "directionUncertainty" &&
701 1 : oComponents[i]->GetType().GetNumericDataType() ==
702 : GDT_Float32)
703 : {
704 : auto poSubArray =
705 3 : poValuesArray->GetView("[\"directionUncertainty\"]");
706 : auto poSubArrayDS = std::unique_ptr<GDALDataset>(
707 2 : poSubArray->AsClassicDataset(1, 0));
708 : auto poSubArrayBand =
709 2 : std::make_unique<S111RasterBand>(std::move(poSubArrayDS));
710 1 : poSubArrayBand->SetDescription("directionUncertainty");
711 1 : poSubArrayBand->m_osUnitType = "degree";
712 1 : poDS->SetBand(poDS->nBands + 1, poSubArrayBand.release());
713 : }
714 : }
715 : }
716 :
717 69 : poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
718 :
719 : // Setup/check for pam .aux.xml.
720 69 : if (osFilename != poOpenInfo->pszFilename)
721 : {
722 36 : poDS->SetSubdatasetName((osFeatureInstance + "/" + osGroup).c_str());
723 36 : poDS->SetPhysicalFilename(osFilename.c_str());
724 : }
725 69 : poDS->SetDescription(poOpenInfo->pszFilename);
726 69 : poDS->TryLoadXML();
727 :
728 : // Setup overviews.
729 69 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
730 :
731 69 : return poDS.release();
732 : }
733 :
734 : /************************************************************************/
735 : /* S111Creator */
736 : /************************************************************************/
737 :
738 : class S111Creator final : public S100BaseWriter
739 : {
740 : public:
741 69 : S111Creator(const char *pszDestFilename, GDALDataset *poSrcDS,
742 : CSLConstList papszOptions)
743 69 : : S100BaseWriter(pszDestFilename, poSrcDS, papszOptions)
744 : {
745 69 : }
746 :
747 : ~S111Creator() override;
748 :
749 : bool Create(GDALProgressFunc pfnProgress, void *pProgressData);
750 :
751 : static constexpr const char *FEATURE_TYPE = "SurfaceCurrent";
752 :
753 : protected:
754 102 : bool Close() override
755 : {
756 102 : return BaseClose();
757 : }
758 :
759 : private:
760 : bool WriteFeatureGroupAttributes();
761 : bool WriteUncertaintyDataset();
762 : bool FillFeatureInstanceGroup(
763 : const std::map<std::string, std::variant<GDALDataset *, std::string>>
764 : &oMapTimestampToDS,
765 : GDALProgressFunc pfnProgress, void *pProgressData);
766 : bool CopyValues(GDALDataset *poSrcDS, GDALProgressFunc pfnProgress,
767 : void *pProgressData);
768 : bool CreateGroupF();
769 : };
770 :
771 : /************************************************************************/
772 : /* S111Creator::~S111Creator() */
773 : /************************************************************************/
774 :
775 69 : S111Creator::~S111Creator()
776 : {
777 69 : S111Creator::Close();
778 69 : }
779 :
780 : /************************************************************************/
781 : /* S111Creator::Create() */
782 : /************************************************************************/
783 :
784 69 : bool S111Creator::Create(GDALProgressFunc pfnProgress, void *pProgressData)
785 : {
786 : CPLStringList aosDatasets(
787 138 : CSLTokenizeString2(m_aosOptions.FetchNameValue("DATASETS"), ",", 0));
788 69 : if (m_poSrcDS->GetRasterCount() == 0 && aosDatasets.empty())
789 : {
790 : // Deal with S111 -> S111 translation;
791 3 : CSLConstList papszSubdatasets = m_poSrcDS->GetMetadata("SUBDATASETS");
792 3 : if (papszSubdatasets)
793 : {
794 2 : int iSubDS = 0;
795 2 : std::string osFirstDataset;
796 2 : std::string osDatasets;
797 20 : for (const auto &[pszItem, pszValue] :
798 22 : cpl::IterateNameValue(papszSubdatasets))
799 : {
800 30 : if (STARTS_WITH(pszItem, "SUBDATASET_") &&
801 15 : cpl::ends_with(std::string_view(pszItem), "_NAME") &&
802 5 : STARTS_WITH(pszValue, "S111:"))
803 : {
804 5 : if (strstr(pszValue, ":SurfaceCurrent."))
805 : {
806 : auto poTmpDS =
807 : std::unique_ptr<GDALDataset>(GDALDataset::Open(
808 : pszValue,
809 2 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
810 2 : if (!poTmpDS)
811 0 : return false;
812 2 : CPLStringList aosOptions(m_aosOptions);
813 2 : if (iSubDS > 0)
814 1 : aosOptions.SetNameValue("APPEND_SUBDATASET", "YES");
815 : S111Creator oAuxCreator(m_osDestFilename.c_str(),
816 : poTmpDS.get(),
817 2 : aosOptions.List());
818 : const int nSubDSCount =
819 2 : ((CSLCount(papszSubdatasets) + 1) / 2);
820 : std::unique_ptr<void,
821 : decltype(&GDALDestroyScaledProgress)>
822 : pScaledProgressData(
823 : GDALCreateScaledProgress(
824 2 : static_cast<double>(iSubDS) / nSubDSCount,
825 2 : static_cast<double>(iSubDS + 1) /
826 : nSubDSCount,
827 : pfnProgress, pProgressData),
828 2 : GDALDestroyScaledProgress);
829 2 : ++iSubDS;
830 2 : if (!oAuxCreator.Create(GDALScaledProgress,
831 : pScaledProgressData.get()))
832 0 : return false;
833 : }
834 : else
835 : {
836 3 : if (osFirstDataset.empty())
837 1 : osFirstDataset = pszValue;
838 3 : if (!osDatasets.empty())
839 2 : osDatasets += ',';
840 3 : osDatasets += pszValue;
841 : }
842 : }
843 : }
844 2 : if (iSubDS > 0)
845 : {
846 1 : return true;
847 : }
848 1 : else if (!osDatasets.empty())
849 : {
850 : auto poTmpDS = std::unique_ptr<GDALDataset>(
851 : GDALDataset::Open(osFirstDataset.c_str(),
852 2 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
853 1 : if (!poTmpDS)
854 0 : return false;
855 2 : CPLStringList aosOptions(m_aosOptions);
856 1 : aosOptions.SetNameValue("DATASETS", osDatasets.c_str());
857 : S111Creator oAuxCreator(m_osDestFilename.c_str(), poTmpDS.get(),
858 2 : aosOptions.List());
859 1 : return oAuxCreator.Create(pfnProgress, pProgressData);
860 : }
861 : }
862 : }
863 :
864 67 : if (m_poSrcDS->GetRasterCount() != 2 && m_poSrcDS->GetRasterCount() != 4)
865 : {
866 17 : CPLError(CE_Failure, CPLE_NotSupported,
867 : "Source dataset %s must have 2 or 4 bands",
868 17 : m_poSrcDS->GetDescription());
869 17 : return false;
870 : }
871 :
872 50 : if (!BaseChecks("S111", /* crsMustBeEPSG = */ false,
873 : /* verticalDatumRequired = */ false))
874 6 : return false;
875 :
876 : std::map<std::string, std::variant<GDALDataset *, std::string>>
877 88 : oMapTimestampToDS;
878 : CPLStringList aosDatasetsTimePoint(CSLTokenizeString2(
879 88 : m_aosOptions.FetchNameValue("DATASETS_TIME_POINT"), ",", 0));
880 44 : if (!aosDatasets.empty())
881 : {
882 2 : if (!aosDatasetsTimePoint.empty() &&
883 0 : aosDatasetsTimePoint.size() != aosDatasets.size())
884 : {
885 0 : CPLError(CE_Failure, CPLE_AppDefined,
886 : "DATASETS_TIME_POINT does not have the same number of "
887 : "values as DATASETS");
888 0 : return false;
889 : }
890 2 : int i = 0;
891 8 : for (const char *pszDataset : aosDatasets)
892 : {
893 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
894 6 : pszDataset, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
895 6 : if (!poDS)
896 0 : return false;
897 12 : if (poDS->GetRasterXSize() != m_poSrcDS->GetRasterXSize() ||
898 6 : poDS->GetRasterYSize() != m_poSrcDS->GetRasterYSize())
899 : {
900 0 : CPLError(CE_Failure, CPLE_NotSupported,
901 : "Dataset %s does not have the same dimensions as %s",
902 0 : poDS->GetDescription(), m_poSrcDS->GetDescription());
903 0 : return false;
904 : }
905 6 : if (poDS->GetRasterCount() != m_poSrcDS->GetRasterCount())
906 : {
907 0 : CPLError(CE_Failure, CPLE_NotSupported,
908 : "Dataset %s must have %d bands",
909 0 : poDS->GetDescription(), m_poSrcDS->GetRasterCount());
910 0 : return false;
911 : }
912 6 : auto poSRS = poDS->GetSpatialRef();
913 6 : if (!poSRS || !poSRS->IsSame(m_poSRS))
914 : {
915 0 : CPLError(CE_Failure, CPLE_NotSupported,
916 : "Dataset %s does not have the same CRS as %s",
917 0 : poDS->GetDescription(), m_poSrcDS->GetDescription());
918 0 : return false;
919 : }
920 6 : GDALGeoTransform gt;
921 6 : if (poDS->GetGeoTransform(gt) != CE_None || gt != m_gt)
922 : {
923 0 : CPLError(CE_Failure, CPLE_NotSupported,
924 : "Dataset %s does not have the same geotransform as %s",
925 0 : poDS->GetDescription(), m_poSrcDS->GetDescription());
926 0 : return false;
927 : }
928 : const char *pszVerticalDatum =
929 6 : poDS->GetMetadataItem("VERTICAL_DATUM");
930 6 : if (pszVerticalDatum)
931 : {
932 : const int nVerticalDatum =
933 0 : S100GetVerticalDatumCodeFromNameOrAbbrev(pszVerticalDatum);
934 0 : if (nVerticalDatum != m_nVerticalDatum)
935 : {
936 0 : CPLError(CE_Failure, CPLE_NotSupported,
937 : "Dataset %s does not have the same vertical datum "
938 : "as %s",
939 0 : poDS->GetDescription(),
940 0 : m_poSrcDS->GetDescription());
941 0 : return false;
942 : }
943 : }
944 6 : const char *pszTimePoint = poDS->GetMetadataItem("timePoint");
945 6 : if (!pszTimePoint && !aosDatasetsTimePoint.empty())
946 0 : pszTimePoint = aosDatasetsTimePoint[i];
947 6 : if (!pszTimePoint)
948 : {
949 0 : CPLError(
950 : CE_Failure, CPLE_NotSupported,
951 : "Dataset %s does not have a timePoint metadata item, and "
952 : "the DATASETS_TIME_POINT creation option is not set",
953 0 : poDS->GetDescription());
954 0 : return false;
955 : }
956 6 : if (strlen(pszTimePoint) != strlen("YYYYMMDDTHHMMSSZ") ||
957 6 : pszTimePoint[8] != 'T' || pszTimePoint[15] != 'Z')
958 : {
959 0 : CPLError(CE_Failure, CPLE_AppDefined,
960 : "timePoint value for dataset %s is %s, but does not "
961 : "conform to a YYYYMMDDTHHMMSSZ datetime value.",
962 0 : poDS->GetDescription(), pszTimePoint);
963 0 : return false;
964 : }
965 6 : if (cpl::contains(oMapTimestampToDS, pszTimePoint))
966 : {
967 0 : CPLError(CE_Failure, CPLE_AppDefined,
968 : "Several datasets are at timePoint %s.", pszTimePoint);
969 0 : return false;
970 : }
971 6 : oMapTimestampToDS[pszTimePoint] = pszDataset;
972 6 : ++i;
973 : }
974 : }
975 :
976 : {
977 44 : const char *pszTimePoint = m_aosOptions.FetchNameValueDef(
978 44 : "TIME_POINT", m_poSrcDS->GetMetadataItem("timePoint"));
979 44 : if (!pszTimePoint)
980 : {
981 3 : CPLError(CE_Failure, CPLE_AppDefined,
982 : "TIME_POINT creation option value must "
983 : "be set, or source dataset must have a timePoint metadata "
984 : "item.");
985 3 : return false;
986 : }
987 41 : if (strlen(pszTimePoint) != strlen("YYYYMMDDTHHMMSSZ") ||
988 40 : pszTimePoint[8] != 'T' || pszTimePoint[15] != 'Z')
989 : {
990 1 : CPLError(CE_Failure, CPLE_AppDefined,
991 : "TIME_POINT creation option value must "
992 : "be set to a YYYYMMDDTHHMMSSZ datetime value.");
993 1 : return false;
994 : }
995 :
996 40 : if (oMapTimestampToDS.empty())
997 : {
998 38 : oMapTimestampToDS[pszTimePoint] = m_poSrcDS;
999 : }
1000 : else
1001 : {
1002 2 : const auto oIter = oMapTimestampToDS.find(pszTimePoint);
1003 3 : if (oIter != oMapTimestampToDS.end() &&
1004 3 : CPLString(std::get<std::string>(oIter->second))
1005 1 : .replaceAll('\\', '/') !=
1006 3 : CPLString(m_poSrcDS->GetDescription())
1007 1 : .replaceAll('\\', '/'))
1008 : {
1009 0 : CPLError(CE_Failure, CPLE_AppDefined,
1010 : "Several datasets are at timePoint %s (%s vs %s).",
1011 : pszTimePoint,
1012 0 : std::get<std::string>(oIter->second).c_str(),
1013 0 : m_poSrcDS->GetDescription());
1014 0 : return false;
1015 : }
1016 : }
1017 : }
1018 40 : if (oMapTimestampToDS.size() > 999)
1019 : {
1020 0 : CPLError(CE_Failure, CPLE_AppDefined,
1021 : "Only up to 999 datasets are supported for a same feature "
1022 : "instance group");
1023 0 : return false;
1024 : }
1025 :
1026 80 : if (m_poSRS->IsVertical() || m_poSRS->IsCompound() || m_poSRS->IsLocal() ||
1027 40 : m_poSRS->GetAxesCount() != 2)
1028 : {
1029 0 : CPLError(CE_Failure, CPLE_NotSupported,
1030 : "The CRS must be a geographic 2D or projected 2D CRS");
1031 0 : return false;
1032 : }
1033 :
1034 : const bool bAppendSubdataset =
1035 40 : CPLTestBool(m_aosOptions.FetchNameValueDef("APPEND_SUBDATASET", "NO"));
1036 40 : if (bAppendSubdataset)
1037 : {
1038 4 : GDALOpenInfo oOpenInfo(m_osDestFilename.c_str(), GA_ReadOnly);
1039 : auto poOriDS =
1040 4 : std::unique_ptr<GDALDataset>(S111Dataset::Open(&oOpenInfo));
1041 2 : if (!poOriDS)
1042 : {
1043 0 : CPLError(CE_Failure, CPLE_AppDefined,
1044 : "%s is not a valid existing S111 dataset",
1045 : m_osDestFilename.c_str());
1046 0 : return false;
1047 : }
1048 2 : const auto poOriSRS = poOriDS->GetSpatialRef();
1049 2 : if (!poOriSRS)
1050 : {
1051 : // shouldn't happen
1052 0 : return false;
1053 : }
1054 2 : if (!poOriSRS->IsSame(m_poSRS))
1055 : {
1056 0 : CPLError(CE_Failure, CPLE_AppDefined,
1057 : "CRS of %s is not the same as the one of %s",
1058 0 : m_osDestFilename.c_str(), m_poSrcDS->GetDescription());
1059 0 : return false;
1060 : }
1061 2 : poOriDS.reset();
1062 :
1063 2 : OGREnvelope sExtent;
1064 2 : if (m_poSrcDS->GetExtentWGS84LongLat(&sExtent) != OGRERR_NONE)
1065 : {
1066 0 : CPLError(CE_Failure, CPLE_AppDefined,
1067 : "Cannot get dataset extent in WGS84 longitude/latitude");
1068 0 : return false;
1069 : }
1070 :
1071 2 : bool ret = OpenFileUpdateMode();
1072 2 : if (ret)
1073 : {
1074 2 : m_featureGroup.reset(H5_CHECK(H5Gopen(m_hdf5, "SurfaceCurrent")));
1075 : }
1076 :
1077 2 : ret = ret && m_featureGroup;
1078 2 : double dfNumInstances = 0;
1079 2 : ret = ret && GH5_FetchAttribute(m_featureGroup, "numInstances",
1080 : dfNumInstances, true);
1081 2 : if (ret && !(dfNumInstances >= 1 && dfNumInstances <= 99 &&
1082 2 : std::round(dfNumInstances) == dfNumInstances))
1083 : {
1084 0 : CPLError(CE_Failure, CPLE_AppDefined,
1085 : "Invalid value for numInstances");
1086 0 : ret = false;
1087 : }
1088 2 : else if (ret && dfNumInstances == 99)
1089 : {
1090 0 : CPLError(CE_Failure, CPLE_AppDefined,
1091 : "Too many existing feature instances");
1092 0 : ret = false;
1093 : }
1094 : else
1095 : {
1096 2 : const int newNumInstances = static_cast<int>(dfNumInstances) + 1;
1097 2 : ret = ret && GH5_WriteAttribute(m_featureGroup, "numInstances",
1098 : newNumInstances);
1099 2 : ret = ret && CreateFeatureInstanceGroup(CPLSPrintf(
1100 : "SurfaceCurrent.%02d", newNumInstances));
1101 2 : ret = ret && FillFeatureInstanceGroup(oMapTimestampToDS,
1102 : pfnProgress, pProgressData);
1103 : }
1104 :
1105 2 : return Close() && ret;
1106 : }
1107 : else
1108 : {
1109 38 : bool ret = CreateFile();
1110 38 : ret = ret && WriteProductSpecification("INT.IHO.S-111.2.0");
1111 38 : ret = ret && WriteIssueDate();
1112 38 : ret = ret && WriteIssueTime(/* bAutogenerateFromCurrent = */ true);
1113 38 : ret = ret && WriteHorizontalCRS();
1114 38 : ret = ret && WriteTopLevelBoundingBox();
1115 :
1116 38 : const char *pszGeographicIdentifier = m_aosOptions.FetchNameValueDef(
1117 : "GEOGRAPHIC_IDENTIFIER",
1118 38 : m_poSrcDS->GetMetadataItem("geographicIdentifier"));
1119 38 : if (pszGeographicIdentifier)
1120 : {
1121 0 : ret =
1122 0 : ret && WriteVarLengthStringValue(m_hdf5, "geographicIdentifier",
1123 : pszGeographicIdentifier);
1124 : }
1125 :
1126 38 : const char *pszDepthTypeIndex = m_aosOptions.FetchNameValueDef(
1127 38 : "DEPTH_TYPE", m_poSrcDS->GetMetadataItem("depthTypeIndex"));
1128 38 : if (!pszDepthTypeIndex)
1129 : {
1130 1 : CPLError(CE_Failure, CPLE_AppDefined,
1131 : "DEPTH_TYPE creation option must be specified.");
1132 1 : return false;
1133 : }
1134 :
1135 : GH5_HIDTypeHolder hDepthTypeIndexEnumDataType(
1136 74 : H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1137 37 : ret = ret && hDepthTypeIndexEnumDataType;
1138 :
1139 : uint8_t val;
1140 37 : val = 1;
1141 37 : ret = ret && H5_CHECK(H5Tenum_insert(hDepthTypeIndexEnumDataType,
1142 : "heightOrDepth", &val)) >= 0;
1143 37 : val = 2;
1144 37 : ret = ret && H5_CHECK(H5Tenum_insert(hDepthTypeIndexEnumDataType,
1145 : "layerAverage", &val)) >= 0;
1146 :
1147 74 : const int nDepthTypeIndex = EQUAL(pszDepthTypeIndex, "heightOrDepth")
1148 74 : ? 1
1149 37 : : EQUAL(pszDepthTypeIndex, "layerAverage")
1150 37 : ? 2
1151 37 : : atoi(pszDepthTypeIndex);
1152 37 : if (nDepthTypeIndex != 1 && nDepthTypeIndex != 2)
1153 : {
1154 1 : CPLError(CE_Failure, CPLE_AppDefined,
1155 : "DEPTH_TYPE creation option must "
1156 : "be set to heightOrDepth/1 or layerAverage/2.");
1157 1 : return false;
1158 : }
1159 72 : ret = ret &&
1160 36 : GH5_CreateAttribute(m_hdf5, "depthTypeIndex",
1161 72 : hDepthTypeIndexEnumDataType) &&
1162 36 : GH5_WriteAttribute(m_hdf5, "depthTypeIndex", nDepthTypeIndex);
1163 :
1164 36 : const char *pszVerticalCS = m_aosOptions.FetchNameValueDef(
1165 36 : "VERTICAL_CS", m_poSrcDS->GetMetadataItem("verticalCS"));
1166 36 : if (!pszVerticalCS)
1167 : {
1168 26 : if (nDepthTypeIndex == 1)
1169 : {
1170 1 : CPLError(CE_Failure, CPLE_AppDefined,
1171 : "VERTICAL_CS creation option must be specified when "
1172 : "DEPTH_TYPE = heightOrDepth");
1173 1 : return false;
1174 : }
1175 : }
1176 : else
1177 : {
1178 20 : const int nVerticalCS = EQUAL(pszVerticalCS, "DEPTH") ? 6498
1179 10 : : EQUAL(pszVerticalCS, "HEIGHT")
1180 10 : ? 6499
1181 10 : : atoi(pszVerticalCS);
1182 10 : if (nVerticalCS != 6498 && nVerticalCS != 6499)
1183 : {
1184 1 : CPLError(
1185 : CE_Failure, CPLE_NotSupported,
1186 : "VERTICAL_CS creation option must be set either to 6498 "
1187 : "(depth/down, metre), or 6499 (height/up, metre)");
1188 1 : return false;
1189 : }
1190 :
1191 9 : ret = ret && WriteVerticalCS(nVerticalCS);
1192 : }
1193 :
1194 34 : ret = ret && WriteVerticalCoordinateBase(2); // verticalDatum
1195 :
1196 34 : if (m_nVerticalDatum > 0)
1197 : {
1198 : // 1=s100VerticalDatum, 2=EPSG
1199 16 : ret = ret && WriteVerticalDatumReference(
1200 8 : m_hdf5, m_nVerticalDatum <= 1024 ? 1 : 2);
1201 16 : ret = ret &&
1202 8 : WriteVerticalDatum(m_hdf5, H5T_STD_I32LE, m_nVerticalDatum);
1203 : }
1204 26 : else if (nDepthTypeIndex == 1)
1205 : {
1206 1 : CPLError(CE_Failure, CPLE_AppDefined,
1207 : "VERTICAL_DATUM creation option must be specified when "
1208 : "DEPTH_TYPE = heightOrDepth");
1209 1 : return false;
1210 : }
1211 :
1212 : const char *pszSurfaceCurrentTrendThreshold =
1213 33 : m_aosOptions.FetchNameValueDef(
1214 : "SURFACE_CURRENT_DEPTH",
1215 33 : m_poSrcDS->GetMetadataItem("surfaceCurrentDepth"));
1216 33 : if (!pszSurfaceCurrentTrendThreshold)
1217 : {
1218 1 : CPLError(CE_Failure, CPLE_AppDefined,
1219 : "SURFACE_CURRENT_DEPTH creation option must be "
1220 : "specified.");
1221 1 : return false;
1222 : }
1223 32 : if (CPLGetValueType(pszSurfaceCurrentTrendThreshold) ==
1224 : CPL_VALUE_STRING)
1225 : {
1226 1 : CPLError(CE_Failure, CPLE_AppDefined,
1227 : "SURFACE_CURRENT_DEPTH creation option value must "
1228 : "be a numeric value.");
1229 1 : return false;
1230 : }
1231 31 : ret =
1232 31 : ret && WriteFloat32Value(m_hdf5, "surfaceCurrentDepth",
1233 : CPLAtof(pszSurfaceCurrentTrendThreshold));
1234 :
1235 31 : const char *pszDatasetDeliveryInterval = m_aosOptions.FetchNameValueDef(
1236 : "DATASET_DELIVERY_INTERVAL",
1237 31 : m_poSrcDS->GetMetadataItem("datasetDeliveryInterval"));
1238 31 : if (pszDatasetDeliveryInterval)
1239 : {
1240 0 : ret = ret &&
1241 0 : WriteVarLengthStringValue(m_hdf5, "datasetDeliveryInterval",
1242 : pszDatasetDeliveryInterval);
1243 : }
1244 :
1245 : // SurfaceCurrent
1246 31 : ret = ret && CreateFeatureGroup(FEATURE_TYPE);
1247 31 : ret = ret && WriteFeatureGroupAttributes();
1248 31 : ret = ret && WriteAxisNames(m_featureGroup);
1249 :
1250 31 : ret = ret && CreateFeatureInstanceGroup("SurfaceCurrent.01");
1251 31 : ret = ret && FillFeatureInstanceGroup(oMapTimestampToDS, pfnProgress,
1252 : pProgressData);
1253 :
1254 31 : ret = ret && CreateGroupF();
1255 :
1256 31 : return Close() && ret;
1257 : }
1258 : }
1259 :
1260 : /************************************************************************/
1261 : /* S111Creator::WriteFeatureGroupAttributes() */
1262 : /************************************************************************/
1263 :
1264 31 : bool S111Creator::WriteFeatureGroupAttributes()
1265 : {
1266 31 : CPLAssert(m_featureGroup);
1267 :
1268 : // 3 = high (recommended)
1269 31 : const char *pszCommonPointRule = m_aosOptions.FetchNameValueDef(
1270 31 : "COMMON_POINT_RULE", m_poSrcDS->GetMetadataItem("commonPointRule"));
1271 31 : if (!pszCommonPointRule)
1272 29 : pszCommonPointRule = "3"; // all (recommended)
1273 31 : const int nCommonPointRule = atoi(pszCommonPointRule);
1274 31 : bool ret = WriteCommonPointRule(m_featureGroup, nCommonPointRule);
1275 31 : ret = ret && WriteDataCodingFormat(m_featureGroup, 2); // Regular grid
1276 31 : ret = ret && WriteDataOffsetCode(m_featureGroup, 5); // Center of cell
1277 31 : ret = ret && WriteDimension(m_featureGroup, 2);
1278 : const char *pszHorizontalPositionUncertainty =
1279 31 : m_aosOptions.FetchNameValueDef(
1280 : "HORIZONTAL_POSITION_UNCERTAINTY",
1281 31 : m_poSrcDS->GetMetadataItem("horizontalPositionUncertainty"));
1282 31 : ret =
1283 62 : ret &&
1284 31 : WriteHorizontalPositionUncertainty(
1285 : m_featureGroup,
1286 0 : pszHorizontalPositionUncertainty &&
1287 0 : pszHorizontalPositionUncertainty[0]
1288 0 : ? static_cast<float>(CPLAtof(pszHorizontalPositionUncertainty))
1289 : : -1.0f);
1290 31 : const char *pszVerticalUncertainty = m_aosOptions.FetchNameValueDef(
1291 : "VERTICAL_UNCERTAINTY",
1292 31 : m_poSrcDS->GetMetadataItem("verticalUncertainty"));
1293 31 : ret = ret && WriteVerticalUncertainty(
1294 : m_featureGroup,
1295 0 : pszVerticalUncertainty && pszVerticalUncertainty[0]
1296 0 : ? static_cast<float>(CPLAtof(pszVerticalUncertainty))
1297 : : -1.0f);
1298 31 : const char *pszTimeUncertainty = m_aosOptions.FetchNameValueDef(
1299 31 : "TIME_UNCERTAINTY", m_poSrcDS->GetMetadataItem("timeUncertainty"));
1300 31 : if (pszTimeUncertainty)
1301 0 : WriteFloat32Value(m_featureGroup, "timeUncertainty",
1302 : CPLAtof(pszTimeUncertainty));
1303 31 : const char *pszMethodCurrentsProduct = m_aosOptions.FetchNameValueDef(
1304 : "METHOD_CURRENTS_PRODUCT",
1305 31 : m_poSrcDS->GetMetadataItem("methodCurrentsProduct"));
1306 31 : if (pszMethodCurrentsProduct)
1307 0 : WriteVarLengthStringValue(m_featureGroup, "methodCurrentsProduct",
1308 : pszMethodCurrentsProduct);
1309 31 : ret = ret && WriteInterpolationType(m_featureGroup, 10); // discrete
1310 31 : ret = ret && WriteNumInstances(m_featureGroup, H5T_STD_U32LE, 1);
1311 62 : ret = ret && WriteSequencingRuleScanDirection(m_featureGroup,
1312 31 : m_poSRS->IsProjected()
1313 : ? "Easting, Northing"
1314 : : "Longitude, Latitude");
1315 31 : ret = ret && WriteSequencingRuleType(m_featureGroup, 1); // Linear
1316 31 : return ret;
1317 : }
1318 :
1319 : /************************************************************************/
1320 : /* S111Creator::WriteUncertaintyDataset() */
1321 : /************************************************************************/
1322 :
1323 27 : bool S111Creator::WriteUncertaintyDataset()
1324 : {
1325 27 : CPLAssert(m_featureInstanceGroup);
1326 :
1327 : GH5_HIDTypeHolder hDataType(
1328 54 : H5_CHECK(H5Tcreate(H5T_COMPOUND, sizeof(char *) + sizeof(float))));
1329 54 : GH5_HIDTypeHolder hVarLengthStringDataType(H5_CHECK(H5Tcopy(H5T_C_S1)));
1330 : bool bRet =
1331 54 : hVarLengthStringDataType &&
1332 54 : H5_CHECK(H5Tset_size(hVarLengthStringDataType, H5T_VARIABLE)) >= 0;
1333 54 : bRet = bRet && hVarLengthStringDataType &&
1334 27 : H5_CHECK(
1335 : H5Tset_strpad(hVarLengthStringDataType, H5T_STR_NULLTERM)) >= 0;
1336 27 : bRet = bRet && hDataType &&
1337 27 : H5_CHECK(H5Tinsert(hDataType, "name", 0,
1338 54 : hVarLengthStringDataType)) >= 0 &&
1339 27 : H5_CHECK(H5Tinsert(hDataType, "value", sizeof(char *),
1340 : H5T_IEEE_F32LE)) >= 0;
1341 27 : constexpr hsize_t NUM_ROWS = 2;
1342 27 : hsize_t dims[] = {NUM_ROWS};
1343 54 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(1, dims, nullptr)));
1344 54 : GH5_HIDDatasetHolder hDatasetID;
1345 54 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1346 27 : bRet = bRet && hParams;
1347 27 : if (bRet)
1348 : {
1349 27 : hDatasetID.reset(
1350 : H5_CHECK(H5Dcreate(m_featureInstanceGroup, "uncertainty", hDataType,
1351 : hDataSpace, hParams)));
1352 27 : bRet = hDatasetID;
1353 : }
1354 :
1355 27 : GH5_HIDSpaceHolder hFileSpace;
1356 27 : if (bRet)
1357 : {
1358 27 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
1359 27 : bRet = hFileSpace;
1360 : }
1361 : GByte abyValues[NUM_ROWS * (sizeof(char *) + sizeof(float))];
1362 : {
1363 27 : const char *pszName = "surfaceCurrentSpeed";
1364 27 : const char *pszVal = m_aosOptions.FetchNameValueDef(
1365 : "UNCERTAINTY_SPEED",
1366 27 : m_poSrcDS->GetMetadataItem("uncertaintySurfaceCurrentSpeed"));
1367 27 : float fVal = pszVal ? static_cast<float>(CPLAtof(pszVal)) : -1.0f;
1368 27 : CPL_LSBPTR32(&fVal);
1369 27 : memcpy(abyValues, &pszName, sizeof(char **));
1370 27 : memcpy(abyValues + sizeof(char *), &fVal, sizeof(fVal));
1371 : }
1372 : {
1373 27 : const char *pszName = "surfaceCurrentDirection";
1374 27 : const char *pszVal = m_aosOptions.FetchNameValueDef(
1375 : "UNCERTAINTY_DIRECTION",
1376 27 : m_poSrcDS->GetMetadataItem("uncertaintySurfaceCurrentDirection"));
1377 27 : float fVal = pszVal ? static_cast<float>(CPLAtof(pszVal)) : -1.0f;
1378 27 : CPL_LSBPTR32(&fVal);
1379 27 : memcpy(abyValues + sizeof(char *) + sizeof(float), &pszName,
1380 : sizeof(char *));
1381 27 : memcpy(abyValues + sizeof(char *) + sizeof(float) + sizeof(char *),
1382 : &fVal, sizeof(fVal));
1383 : }
1384 :
1385 27 : H5OFFSET_TYPE offset[] = {0};
1386 54 : bRet = bRet &&
1387 27 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1388 54 : nullptr, dims, nullptr)) >= 0 &&
1389 27 : H5_CHECK(H5Dwrite(hDatasetID, hDataType, hDataSpace, hFileSpace,
1390 : H5P_DEFAULT, abyValues)) >= 0;
1391 54 : return bRet;
1392 : }
1393 :
1394 : /************************************************************************/
1395 : /* S111Creator::FillFeatureInstanceGroup() */
1396 : /************************************************************************/
1397 :
1398 33 : bool S111Creator::FillFeatureInstanceGroup(
1399 : const std::map<std::string, std::variant<GDALDataset *, std::string>>
1400 : &oMapTimestampToDS,
1401 : GDALProgressFunc pfnProgress, void *pProgressData)
1402 : {
1403 33 : bool ret = WriteFIGGridRelatedParameters(m_featureInstanceGroup);
1404 :
1405 33 : const int numInstances = static_cast<int>(oMapTimestampToDS.size());
1406 :
1407 33 : ret =
1408 33 : ret && WriteNumGRP(m_featureInstanceGroup, H5T_STD_U32LE, numInstances);
1409 33 : ret = ret && WriteUInt32Value(m_featureInstanceGroup, "numberOfTimes",
1410 : numInstances);
1411 :
1412 : // Check if value groups are spaced at a regular time interval
1413 33 : GIntBig nLastInterval = 0;
1414 33 : GIntBig nLastTS = 0;
1415 70 : for (const auto &[key, value] : oMapTimestampToDS)
1416 : {
1417 37 : CPL_IGNORE_RET_VAL(value);
1418 : int nYear, nMonth, nDay, nHour, nMinute, nSecond;
1419 37 : if (sscanf(key.c_str(), "%04d%02d%02dT%02d%02d%02dZ", &nYear, &nMonth,
1420 37 : &nDay, &nHour, &nMinute, &nSecond) == 6)
1421 : {
1422 : struct tm brokenDown;
1423 37 : memset(&brokenDown, 0, sizeof(brokenDown));
1424 37 : brokenDown.tm_year = nYear - 1900;
1425 37 : brokenDown.tm_mon = nMonth - 1;
1426 37 : brokenDown.tm_mday = nDay;
1427 37 : brokenDown.tm_hour = nHour;
1428 37 : brokenDown.tm_min = nMinute;
1429 37 : brokenDown.tm_sec = nMinute;
1430 37 : const GIntBig nTS = CPLYMDHMSToUnixTime(&brokenDown);
1431 37 : if (nLastTS != 0)
1432 : {
1433 4 : if (nLastInterval == 0)
1434 : {
1435 2 : nLastInterval = nTS - nLastTS;
1436 : }
1437 2 : else if (nLastInterval != nTS - nLastTS)
1438 : {
1439 0 : nLastInterval = 0;
1440 0 : break;
1441 : }
1442 : }
1443 37 : nLastTS = nTS;
1444 : }
1445 : }
1446 :
1447 33 : const char *pszTimeRecordInterval = m_aosOptions.FetchNameValueDef(
1448 : "TIME_RECORD_INTERVAL",
1449 33 : m_poSrcDS->GetMetadataItem("timeRecordInterval"));
1450 33 : if (pszTimeRecordInterval)
1451 : {
1452 2 : ret = ret &&
1453 1 : WriteUInt16Value(m_featureInstanceGroup, "timeRecordInterval",
1454 : atoi(pszTimeRecordInterval));
1455 : }
1456 32 : else if (nLastInterval > 0 && nLastInterval < 65536)
1457 : {
1458 2 : ret = ret &&
1459 1 : WriteUInt16Value(m_featureInstanceGroup, "timeRecordInterval",
1460 : static_cast<int>(nLastInterval));
1461 : }
1462 :
1463 66 : ret = ret && WriteVarLengthStringValue(
1464 : m_featureInstanceGroup, "dateTimeOfFirstRecord",
1465 33 : oMapTimestampToDS.begin()->first.c_str());
1466 66 : ret = ret && WriteVarLengthStringValue(
1467 : m_featureInstanceGroup, "dateTimeOfLastRecord",
1468 33 : oMapTimestampToDS.rbegin()->first.c_str());
1469 :
1470 33 : const char *pszDataDynamicity = m_aosOptions.FetchNameValueDef(
1471 33 : "DATA_DYNAMICITY", m_poSrcDS->GetMetadataItem("dataDynamicity"));
1472 33 : if (!pszDataDynamicity)
1473 : {
1474 1 : CPLError(CE_Failure, CPLE_AppDefined,
1475 : "DATA_DYNAMICITY creation option must "
1476 : "be specified.");
1477 1 : return false;
1478 : }
1479 : {
1480 : GH5_HIDTypeHolder hDepthTypeIndexEnumDataType(
1481 32 : H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1482 32 : ret = ret && hDepthTypeIndexEnumDataType;
1483 :
1484 : uint8_t val;
1485 32 : val = 1;
1486 32 : ret = ret && H5_CHECK(H5Tenum_insert(hDepthTypeIndexEnumDataType,
1487 : "observation", &val)) >= 0;
1488 32 : val = 2;
1489 64 : ret = ret &&
1490 32 : H5_CHECK(H5Tenum_insert(hDepthTypeIndexEnumDataType,
1491 : "astronomicalPrediction", &val)) >= 0;
1492 32 : val = 3;
1493 32 : ret = ret && H5_CHECK(H5Tenum_insert(hDepthTypeIndexEnumDataType,
1494 : "analysisOrHybrid", &val)) >= 0;
1495 32 : val = 5;
1496 32 : ret =
1497 32 : ret && H5_CHECK(H5Tenum_insert(hDepthTypeIndexEnumDataType,
1498 : "hydrodynamicForecast", &val)) >= 0;
1499 :
1500 32 : const int nDataDynamicity =
1501 64 : EQUAL(pszDataDynamicity, "observation") ? 1
1502 64 : : EQUAL(pszDataDynamicity, "astronomicalPrediction") ? 2
1503 64 : : EQUAL(pszDataDynamicity, "analysisOrHybrid") ? 3
1504 32 : : EQUAL(pszDataDynamicity, "hydrodynamicForecast")
1505 32 : ? 5
1506 32 : : atoi(pszDataDynamicity);
1507 32 : if (nDataDynamicity != 1 && nDataDynamicity != 2 &&
1508 32 : nDataDynamicity != 3 && nDataDynamicity != 5)
1509 : {
1510 1 : CPLError(CE_Failure, CPLE_AppDefined,
1511 : "DATA_DYNAMICITY creation option must "
1512 : "be set to observation/1, astronomicalPrediction/2, "
1513 : "analysisOrHybrid/3 or hydrodynamicForecast/5.");
1514 1 : return false;
1515 : }
1516 31 : ret = ret &&
1517 31 : GH5_CreateAttribute(m_featureInstanceGroup, "dataDynamicity",
1518 62 : hDepthTypeIndexEnumDataType) &&
1519 31 : GH5_WriteAttribute(m_featureInstanceGroup, "dataDynamicity",
1520 : nDataDynamicity);
1521 : }
1522 :
1523 31 : if (m_poSrcDS->GetRasterCount() == 2 ||
1524 35 : m_aosOptions.FetchNameValue("UNCERTAINTY_SPEED") ||
1525 4 : m_aosOptions.FetchNameValue("UNCERTAINTY_DIRECTION"))
1526 : {
1527 27 : ret = ret && WriteUncertaintyDataset();
1528 : }
1529 :
1530 31 : int iInstance = 0;
1531 31 : double dfLastRatio = 0;
1532 66 : for (const auto &iter : oMapTimestampToDS)
1533 : {
1534 35 : ++iInstance;
1535 35 : ret = ret && CreateValuesGroup(CPLSPrintf("Group_%03d", iInstance));
1536 :
1537 35 : ret = ret && WriteVarLengthStringValue(m_valuesGroup, "timePoint",
1538 : iter.first.c_str());
1539 :
1540 0 : std::unique_ptr<GDALDataset> poTmpDSHolder;
1541 : GDALDataset *poSrcDS;
1542 35 : if (std::holds_alternative<std::string>(iter.second))
1543 : {
1544 6 : poTmpDSHolder.reset(
1545 6 : GDALDataset::Open(std::get<std::string>(iter.second).c_str(),
1546 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
1547 6 : if (!poTmpDSHolder)
1548 : {
1549 0 : return false;
1550 : }
1551 6 : poSrcDS = poTmpDSHolder.get();
1552 : }
1553 : else
1554 : {
1555 29 : CPLAssert(std::holds_alternative<GDALDataset *>(iter.second));
1556 29 : poSrcDS = std::get<GDALDataset *>(iter.second);
1557 : }
1558 :
1559 35 : const double dfNewRatio = static_cast<double>(iInstance) / numInstances;
1560 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>
1561 : pScaledProgressData(
1562 : GDALCreateScaledProgress(dfLastRatio, dfNewRatio, pfnProgress,
1563 : pProgressData),
1564 35 : GDALDestroyScaledProgress);
1565 35 : ret = ret && CopyValues(poSrcDS, GDALScaledProgress,
1566 : pScaledProgressData.get());
1567 35 : dfLastRatio = dfNewRatio;
1568 : }
1569 :
1570 31 : return ret;
1571 : }
1572 :
1573 : /************************************************************************/
1574 : /* S111Creator::CreateGroupF() */
1575 : /************************************************************************/
1576 :
1577 : // Per S-111 v2.0 spec
1578 :
1579 : constexpr float string_to_float(std::string_view str)
1580 : {
1581 : if (str.empty())
1582 : return 0.0f;
1583 :
1584 : size_t i = 0;
1585 : bool neg = false;
1586 : if (str[0] == '-')
1587 : {
1588 : neg = true;
1589 : ++i;
1590 : }
1591 : else if (str[0] == '+')
1592 : {
1593 : ++i;
1594 : }
1595 :
1596 : float int_part = 0.0f;
1597 : for (; i < str.size() && str[i] >= '0' && str[i] <= '9'; ++i)
1598 : {
1599 : int_part = int_part * 10.0f + (str[i] - '0');
1600 : }
1601 :
1602 : float frac_part = 0.0f;
1603 : float divisor = 1.0f;
1604 : if (i < str.size() && str[i] == '.')
1605 : {
1606 : for (++i; i < str.size() && str[i] >= '0' && str[i] <= '9'; ++i)
1607 : {
1608 : frac_part = frac_part * 10.0f + (str[i] - '0');
1609 : divisor *= 10.0f;
1610 : }
1611 : }
1612 :
1613 : float result = int_part + frac_part / divisor;
1614 : return neg ? -result : result;
1615 : }
1616 :
1617 : constexpr const char *MIN_SPEED_STR = "0.00";
1618 : constexpr float MIN_SPEED = string_to_float(MIN_SPEED_STR);
1619 : constexpr const char *MAX_SPEED_STR = "99.00";
1620 : constexpr float MAX_SPEED = string_to_float(MAX_SPEED_STR);
1621 : constexpr const char *NODATA_SPEED_STR = "-9999.00";
1622 : constexpr float NODATA_SPEED = string_to_float(NODATA_SPEED_STR);
1623 :
1624 : constexpr const char *MIN_DIR_STR = "0.0";
1625 : constexpr float MIN_DIR = string_to_float(MIN_DIR_STR);
1626 : constexpr const char *MAX_DIR_STR = "359.9";
1627 : constexpr float MAX_DIR = string_to_float(MAX_DIR_STR);
1628 : constexpr const char *NODATA_DIR_STR = "-9999.0";
1629 : constexpr float NODATA_DIR = string_to_float(NODATA_DIR_STR);
1630 :
1631 : constexpr const char *NODATA_UNCT_STR = "-1.0";
1632 : constexpr float NODATA_UNCT = string_to_float(NODATA_UNCT_STR);
1633 :
1634 29 : bool S111Creator::CreateGroupF()
1635 : {
1636 29 : bool ret = S100BaseWriter::CreateGroupF();
1637 :
1638 29 : CPLStringList aosFeatureCodes;
1639 29 : aosFeatureCodes.push_back(FEATURE_TYPE);
1640 58 : ret = ret && WriteOneDimensionalVarLengthStringArray(
1641 29 : m_GroupF, "featureCode", aosFeatureCodes.List());
1642 :
1643 : {
1644 : std::vector<std::array<const char *, GROUP_F_DATASET_FIELD_COUNT>> rows{
1645 : {"surfaceCurrentSpeed", "Surface Current Speed", "knot",
1646 : NODATA_SPEED_STR, "H5T_FLOAT", MIN_SPEED_STR, MAX_SPEED_STR,
1647 : "geSemiInterval"},
1648 : {"surfaceCurrentDirection", "Surface Current Direction", "degree",
1649 : "-9999.0", "H5T_FLOAT", MIN_DIR_STR, MAX_DIR_STR,
1650 : "closedInterval"},
1651 : {"speedUncertainty", "Speed Uncertainty", "knot", NODATA_UNCT_STR,
1652 : "H5T_FLOAT", MIN_SPEED_STR, MAX_SPEED_STR, "geSemiInterval"},
1653 : {"directionUncertainty", "Direction Uncertainty", "degree",
1654 : NODATA_UNCT_STR, "H5T_FLOAT", MIN_DIR_STR, MAX_DIR_STR,
1655 : "closedInterval"},
1656 29 : };
1657 29 : rows.resize(m_poSrcDS->GetRasterCount());
1658 29 : ret = ret && WriteGroupFDataset(FEATURE_TYPE, rows);
1659 : }
1660 :
1661 58 : return ret;
1662 : }
1663 :
1664 : /************************************************************************/
1665 : /* S111Creator::CopyValues() */
1666 : /************************************************************************/
1667 :
1668 35 : bool S111Creator::CopyValues(GDALDataset *poSrcDS, GDALProgressFunc pfnProgress,
1669 : void *pProgressData)
1670 : {
1671 35 : CPLAssert(m_valuesGroup.get() >= 0);
1672 :
1673 35 : const int nYSize = poSrcDS->GetRasterYSize();
1674 35 : const int nXSize = poSrcDS->GetRasterXSize();
1675 :
1676 35 : hsize_t dims[] = {static_cast<hsize_t>(nYSize),
1677 35 : static_cast<hsize_t>(nXSize)};
1678 :
1679 70 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(2, dims, nullptr)));
1680 35 : bool bRet = hDataSpace;
1681 :
1682 : const bool bDeflate =
1683 35 : EQUAL(m_aosOptions.FetchNameValueDef("COMPRESS", "DEFLATE"), "DEFLATE");
1684 : const int nCompressionLevel =
1685 35 : atoi(m_aosOptions.FetchNameValueDef("ZLEVEL", "6"));
1686 : const int nBlockSize =
1687 35 : std::min(4096, std::max(100, atoi(m_aosOptions.FetchNameValueDef(
1688 35 : "BLOCK_SIZE", "100"))));
1689 35 : const int nBlockXSize = std::min(nXSize, nBlockSize);
1690 35 : const int nBlockYSize = std::min(nYSize, nBlockSize);
1691 35 : const int nComponents = poSrcDS->GetRasterCount();
1692 35 : CPLAssert(nComponents == 2 || nComponents == 4);
1693 :
1694 : GH5_HIDTypeHolder hDataType(
1695 70 : H5_CHECK(H5Tcreate(H5T_COMPOUND, sizeof(float) * nComponents)));
1696 35 : bRet = bRet && hDataType &&
1697 35 : H5_CHECK(H5Tinsert(hDataType, "surfaceCurrentSpeed", 0,
1698 70 : H5T_IEEE_F32LE)) >= 0 &&
1699 35 : H5_CHECK(H5Tinsert(hDataType, "surfaceCurrentDirection",
1700 : sizeof(float), H5T_IEEE_F32LE)) >= 0;
1701 35 : if (nComponents == 4 && bRet)
1702 : {
1703 8 : bRet = H5_CHECK(H5Tinsert(hDataType, "speedUncertainty",
1704 8 : 2 * sizeof(float), H5T_IEEE_F32LE)) >= 0 &&
1705 4 : H5_CHECK(H5Tinsert(hDataType, "directionUncertainty",
1706 : 3 * sizeof(float), H5T_IEEE_F32LE)) >= 0;
1707 : }
1708 :
1709 35 : hsize_t chunk_size[] = {static_cast<hsize_t>(nBlockYSize),
1710 35 : static_cast<hsize_t>(nBlockXSize)};
1711 :
1712 70 : GH5_HIDParametersHolder hParams(H5_CHECK(H5Pcreate(H5P_DATASET_CREATE)));
1713 35 : bRet = bRet && hParams &&
1714 35 : H5_CHECK(H5Pset_fill_time(hParams, H5D_FILL_TIME_ALLOC)) >= 0 &&
1715 105 : H5_CHECK(H5Pset_layout(hParams, H5D_CHUNKED)) >= 0 &&
1716 35 : H5_CHECK(H5Pset_chunk(hParams, 2, chunk_size)) >= 0;
1717 :
1718 35 : if (bRet && bDeflate)
1719 : {
1720 34 : bRet = H5_CHECK(H5Pset_deflate(hParams, nCompressionLevel)) >= 0;
1721 : }
1722 :
1723 70 : GH5_HIDDatasetHolder hDatasetID;
1724 35 : if (bRet)
1725 : {
1726 35 : hDatasetID.reset(H5_CHECK(H5Dcreate(m_valuesGroup, "values", hDataType,
1727 : hDataSpace, hParams)));
1728 35 : bRet = hDatasetID;
1729 : }
1730 :
1731 70 : GH5_HIDSpaceHolder hFileSpace;
1732 35 : if (bRet)
1733 : {
1734 35 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
1735 35 : bRet = hFileSpace;
1736 : }
1737 :
1738 35 : const int nYBlocks = static_cast<int>(DIV_ROUND_UP(nYSize, nBlockYSize));
1739 35 : const int nXBlocks = static_cast<int>(DIV_ROUND_UP(nXSize, nBlockXSize));
1740 35 : std::vector<float> afValues(static_cast<size_t>(nBlockYSize) * nBlockXSize *
1741 35 : nComponents);
1742 35 : const bool bReverseY = m_gt[5] < 0;
1743 :
1744 35 : constexpr std::array<float, 4> afNoDataValue{NODATA_SPEED, NODATA_DIR,
1745 : NODATA_UNCT, NODATA_UNCT};
1746 :
1747 35 : constexpr float INF = std::numeric_limits<float>::infinity();
1748 35 : std::array<float, 4> afMin{INF, INF, INF, INF};
1749 35 : std::array<float, 4> afMax{-INF, -INF, -INF, -INF};
1750 :
1751 35 : std::array<int, 4> abHasNoDataBand{FALSE, FALSE, FALSE, FALSE};
1752 35 : std::array<float, 4> afSrcNoData{0, 0, 0, 0};
1753 113 : for (int i = 0; i < nComponents; ++i)
1754 : {
1755 78 : afSrcNoData[i] = static_cast<float>(
1756 78 : poSrcDS->GetRasterBand(i + 1)->GetNoDataValue(&abHasNoDataBand[i]));
1757 : }
1758 :
1759 81 : for (int iY = 0; iY < nYBlocks && bRet; iY++)
1760 : {
1761 : const int nSrcYOff = bReverseY
1762 46 : ? std::max(0, nYSize - (iY + 1) * nBlockYSize)
1763 41 : : iY * nBlockYSize;
1764 46 : const int nReqCountY = std::min(nBlockYSize, nYSize - iY * nBlockYSize);
1765 224 : for (int iX = 0; iX < nXBlocks && bRet; iX++)
1766 : {
1767 : const int nReqCountX =
1768 178 : std::min(nBlockXSize, nXSize - iX * nBlockXSize);
1769 :
1770 178 : bRet =
1771 178 : poSrcDS->RasterIO(
1772 178 : GF_Read, iX * nBlockXSize, nSrcYOff, nReqCountX, nReqCountY,
1773 5 : bReverseY ? afValues.data() +
1774 5 : (nReqCountY - 1) * nReqCountX * nComponents
1775 173 : : afValues.data(),
1776 : nReqCountX, nReqCountY, GDT_Float32, nComponents, nullptr,
1777 178 : static_cast<int>(sizeof(float)) * nComponents,
1778 : bReverseY ? -static_cast<GPtrDiff_t>(sizeof(float)) *
1779 5 : nComponents * nReqCountX
1780 : : 0,
1781 : sizeof(float), nullptr) == CE_None;
1782 :
1783 178 : if (bRet)
1784 : {
1785 1440470 : for (int i = 0; i < nReqCountY * nReqCountX; i++)
1786 : {
1787 7200920 : for (int iC = 0; iC < nComponents; ++iC)
1788 : {
1789 5760630 : float fVal = afValues[i * nComponents + iC];
1790 11521300 : if ((abHasNoDataBand[iC] && fVal == afSrcNoData[iC]) ||
1791 5760630 : std::isnan(fVal))
1792 : {
1793 0 : afValues[i * nComponents + iC] = afNoDataValue[iC];
1794 : }
1795 : else
1796 : {
1797 5760630 : afMin[iC] = std::min(afMin[iC], fVal);
1798 5760630 : afMax[iC] = std::max(afMax[iC], fVal);
1799 : }
1800 5760630 : CPL_LSBPTR32(&afValues[i * nComponents + iC]);
1801 : }
1802 : }
1803 : }
1804 :
1805 : H5OFFSET_TYPE offset[] = {
1806 178 : static_cast<H5OFFSET_TYPE>(iY) *
1807 178 : static_cast<H5OFFSET_TYPE>(nBlockYSize),
1808 178 : static_cast<H5OFFSET_TYPE>(iX) *
1809 178 : static_cast<H5OFFSET_TYPE>(nBlockXSize)};
1810 178 : hsize_t count[2] = {static_cast<hsize_t>(nReqCountY),
1811 178 : static_cast<hsize_t>(nReqCountX)};
1812 : GH5_HIDSpaceHolder hMemSpace(
1813 178 : H5_CHECK(H5Screate_simple(2, count, nullptr)));
1814 178 : bRet =
1815 178 : bRet &&
1816 178 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
1817 178 : nullptr, count, nullptr)) >= 0 &&
1818 178 : hMemSpace &&
1819 178 : H5_CHECK(H5Dwrite(hDatasetID, hDataType, hMemSpace, hFileSpace,
1820 356 : H5P_DEFAULT, afValues.data())) >= 0 &&
1821 178 : pfnProgress((static_cast<double>(iY) * nXBlocks + iX + 1) /
1822 178 : (static_cast<double>(nXBlocks) * nYBlocks),
1823 : "", pProgressData) != 0;
1824 : }
1825 : }
1826 :
1827 35 : constexpr int IDX_SPEED = 0;
1828 35 : constexpr int IDX_DIR = 1;
1829 35 : constexpr int IDX_UNC_SPEED = 2;
1830 35 : constexpr int IDX_UNC_DIR = 3;
1831 :
1832 35 : if (afMin[IDX_SPEED] > afMax[IDX_SPEED])
1833 : {
1834 0 : afMin[IDX_SPEED] = afMax[IDX_SPEED] = NODATA_SPEED;
1835 : }
1836 35 : else if (!(afMin[IDX_SPEED] >= MIN_SPEED && afMax[IDX_SPEED] <= MAX_SPEED))
1837 : {
1838 2 : CPLError(CE_Warning, CPLE_AppDefined,
1839 : "Range of surface current speed in the dataset is [%f, %f] "
1840 : "whereas the "
1841 : "allowed range is [%.2f, %.2f]",
1842 2 : afMin[IDX_SPEED], afMax[IDX_SPEED], MIN_SPEED, MAX_SPEED);
1843 : }
1844 :
1845 35 : if (afMin[IDX_DIR] > afMax[IDX_DIR])
1846 : {
1847 0 : afMin[IDX_DIR] = afMax[IDX_DIR] = NODATA_DIR;
1848 : }
1849 35 : else if (!(afMin[IDX_DIR] >= MIN_DIR && afMax[IDX_DIR] <= MAX_DIR))
1850 : {
1851 2 : CPLError(
1852 : CE_Warning, CPLE_AppDefined,
1853 : "Range of surface current direction in the dataset is [%f, %f] "
1854 : "whereas the "
1855 : "allowed range is [%.2f, %.2f]",
1856 2 : afMin[IDX_DIR], afMax[IDX_DIR], MIN_DIR, MAX_DIR);
1857 : }
1858 :
1859 39 : if (afMax[IDX_UNC_SPEED] >= afMin[IDX_UNC_SPEED] &&
1860 4 : afMin[IDX_UNC_SPEED] < 0)
1861 : {
1862 1 : CPLError(
1863 : CE_Warning, CPLE_AppDefined,
1864 : "Negative speed uncertainty value found (%f), which is not allowed "
1865 : "(except nodata value %s)",
1866 1 : afMin[IDX_UNC_SPEED], NODATA_UNCT_STR);
1867 : }
1868 :
1869 35 : if (afMax[IDX_UNC_DIR] >= afMin[IDX_UNC_DIR] && afMin[IDX_UNC_DIR] < 0)
1870 : {
1871 1 : CPLError(CE_Warning, CPLE_AppDefined,
1872 : "Negative direction uncertainty value found (%f), which is "
1873 : "not allowed (except nodata value %s)",
1874 1 : afMin[IDX_UNC_DIR], NODATA_UNCT_STR);
1875 : }
1876 :
1877 35 : if (bRet)
1878 : {
1879 35 : double prevMinSpeed = 0;
1880 35 : double prevMaxSpeed = 0;
1881 35 : if (GH5_FetchAttribute(m_featureGroup, "minDatasetCurrentSpeed",
1882 41 : prevMinSpeed) &&
1883 6 : GH5_FetchAttribute(m_featureGroup, "maxDatasetCurrentSpeed",
1884 : prevMaxSpeed))
1885 : {
1886 6 : if (afMin[IDX_SPEED] != NODATA_SPEED)
1887 : {
1888 6 : prevMinSpeed = std::min<double>(prevMinSpeed, afMin[IDX_SPEED]);
1889 6 : prevMaxSpeed = std::max<double>(prevMaxSpeed, afMax[IDX_SPEED]);
1890 6 : bRet =
1891 6 : GH5_WriteAttribute(m_featureGroup, "minDatasetCurrentSpeed",
1892 12 : prevMinSpeed) &&
1893 6 : GH5_WriteAttribute(m_featureGroup, "maxDatasetCurrentSpeed",
1894 : prevMaxSpeed);
1895 : }
1896 : }
1897 : else
1898 : {
1899 58 : bRet = WriteFloat64Value(m_featureGroup, "minDatasetCurrentSpeed",
1900 87 : afMin[IDX_SPEED]) &&
1901 29 : WriteFloat64Value(m_featureGroup, "maxDatasetCurrentSpeed",
1902 29 : afMax[IDX_SPEED]);
1903 : }
1904 : }
1905 :
1906 70 : return bRet;
1907 : }
1908 :
1909 : /************************************************************************/
1910 : /* S111Dataset::CreateCopy() */
1911 : /************************************************************************/
1912 :
1913 : /* static */
1914 66 : GDALDataset *S111Dataset::CreateCopy(const char *pszFilename,
1915 : GDALDataset *poSrcDS, int /* bStrict*/,
1916 : char **papszOptions,
1917 : GDALProgressFunc pfnProgress,
1918 : void *pProgressData)
1919 : {
1920 132 : S111Creator creator(pszFilename, poSrcDS, papszOptions);
1921 66 : if (!creator.Create(pfnProgress, pProgressData))
1922 36 : return nullptr;
1923 :
1924 : VSIStatBufL sStatBuf;
1925 60 : if (VSIStatL(pszFilename, &sStatBuf) == 0 &&
1926 30 : sStatBuf.st_size > 10 * 1024 * 1024)
1927 : {
1928 1 : CPLError(CE_Warning, CPLE_AppDefined,
1929 : "%s file size exceeds 10 MB, which is the upper limit "
1930 : "suggested for wireless transmission to marine vessels",
1931 : pszFilename);
1932 : }
1933 :
1934 60 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
1935 30 : return Open(&oOpenInfo);
1936 : }
1937 :
1938 : /************************************************************************/
1939 : /* S111DatasetDriverUnload() */
1940 : /************************************************************************/
1941 :
1942 9 : static void S111DatasetDriverUnload(GDALDriver *)
1943 : {
1944 9 : HDF5UnloadFileDriver();
1945 9 : }
1946 :
1947 : /************************************************************************/
1948 : /* GDALRegister_S111() */
1949 : /************************************************************************/
1950 14 : void GDALRegister_S111()
1951 :
1952 : {
1953 14 : if (!GDAL_CHECK_VERSION("S111"))
1954 0 : return;
1955 :
1956 14 : if (GDALGetDriverByName(S111_DRIVER_NAME) != nullptr)
1957 0 : return;
1958 :
1959 14 : GDALDriver *poDriver = new GDALDriver();
1960 :
1961 14 : S111DriverSetCommonMetadata(poDriver);
1962 14 : poDriver->pfnOpen = S111Dataset::Open;
1963 14 : poDriver->pfnCreateCopy = S111Dataset::CreateCopy;
1964 14 : poDriver->pfnUnloadDriver = S111DatasetDriverUnload;
1965 :
1966 14 : GetGDALDriverManager()->RegisterDriver(poDriver);
1967 : }
|
