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