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