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