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