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