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, 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 <limits>
26 : #include <map>
27 :
28 : /************************************************************************/
29 : /* S104Dataset */
30 : /************************************************************************/
31 :
32 42 : class S104Dataset final : public S100BaseDataset
33 : {
34 : public:
35 21 : explicit S104Dataset(const std::string &osFilename)
36 21 : : S100BaseDataset(osFilename)
37 : {
38 21 : }
39 :
40 : ~S104Dataset() override;
41 :
42 : static GDALDataset *Open(GDALOpenInfo *);
43 : };
44 :
45 : S104Dataset::~S104Dataset() = default;
46 :
47 : /************************************************************************/
48 : /* S104RasterBand */
49 : /************************************************************************/
50 :
51 : class S104RasterBand final : public GDALProxyRasterBand
52 : {
53 : friend class S104Dataset;
54 : std::unique_ptr<GDALDataset> m_poDS{};
55 : GDALRasterBand *m_poUnderlyingBand = nullptr;
56 : std::string m_osUnitType{};
57 : std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
58 :
59 : CPL_DISALLOW_COPY_ASSIGN(S104RasterBand)
60 :
61 : public:
62 8 : explicit S104RasterBand(std::unique_ptr<GDALDataset> &&poDSIn)
63 16 : : m_poDS(std::move(poDSIn)),
64 8 : m_poUnderlyingBand(m_poDS->GetRasterBand(1))
65 : {
66 8 : eDataType = m_poUnderlyingBand->GetRasterDataType();
67 8 : m_poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
68 8 : }
69 :
70 : GDALRasterBand *
71 : RefUnderlyingRasterBand(bool /*bForceOpen*/ = true) const override;
72 :
73 2 : const char *GetUnitType() override
74 : {
75 2 : return m_osUnitType.c_str();
76 : }
77 :
78 1 : GDALRasterAttributeTable *GetDefaultRAT() override
79 : {
80 1 : return m_poRAT.get();
81 : }
82 : };
83 :
84 : GDALRasterBand *
85 12 : S104RasterBand::RefUnderlyingRasterBand(bool /*bForceOpen*/) const
86 : {
87 12 : return m_poUnderlyingBand;
88 : }
89 :
90 : /************************************************************************/
91 : /* Open() */
92 : /************************************************************************/
93 :
94 22 : GDALDataset *S104Dataset::Open(GDALOpenInfo *poOpenInfo)
95 :
96 : {
97 : // Confirm that this appears to be a S104 file.
98 22 : if (!S104DatasetIdentify(poOpenInfo))
99 0 : return nullptr;
100 :
101 : HDF5_GLOBAL_LOCK();
102 :
103 22 : if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
104 : {
105 1 : return HDF5Dataset::OpenMultiDim(poOpenInfo);
106 : }
107 :
108 : // Confirm the requested access is supported.
109 21 : if (poOpenInfo->eAccess == GA_Update)
110 : {
111 0 : ReportUpdateNotSupportedByDriver("S104");
112 0 : return nullptr;
113 : }
114 :
115 42 : std::string osFilename(poOpenInfo->pszFilename);
116 42 : std::string osFeatureInstance = "WaterLevel.01";
117 42 : std::string osGroup;
118 21 : if (STARTS_WITH(poOpenInfo->pszFilename, "S104:"))
119 : {
120 : const CPLStringList aosTokens(
121 5 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
122 5 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
123 :
124 5 : if (aosTokens.size() == 2)
125 : {
126 0 : osFilename = aosTokens[1];
127 : }
128 5 : else if (aosTokens.size() == 3)
129 : {
130 3 : osFilename = aosTokens[1];
131 3 : osGroup = aosTokens[2];
132 : }
133 2 : else if (aosTokens.size() == 4)
134 : {
135 2 : osFilename = aosTokens[1];
136 2 : osFeatureInstance = aosTokens[2];
137 2 : osGroup = aosTokens[3];
138 : }
139 : else
140 : {
141 0 : return nullptr;
142 : }
143 : }
144 :
145 42 : auto poDS = std::make_unique<S104Dataset>(osFilename);
146 21 : if (!poDS->Init())
147 0 : return nullptr;
148 :
149 21 : const auto &poRootGroup = poDS->m_poRootGroup;
150 :
151 63 : auto poWaterLevel = poRootGroup->OpenGroup("WaterLevel");
152 21 : if (!poWaterLevel)
153 : {
154 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find /WaterLevel group");
155 0 : return nullptr;
156 : }
157 :
158 63 : auto poDataCodingFormat = poWaterLevel->GetAttribute("dataCodingFormat");
159 42 : if (!poDataCodingFormat ||
160 21 : poDataCodingFormat->GetDataType().GetClass() != GEDTC_NUMERIC)
161 : {
162 0 : CPLError(CE_Failure, CPLE_AppDefined,
163 : "Cannot find /WaterLevel/dataCodingFormat attribute");
164 0 : return nullptr;
165 : }
166 21 : const int nDataCodingFormat = poDataCodingFormat->ReadAsInt();
167 21 : if (nDataCodingFormat != 2)
168 : {
169 0 : CPLError(CE_Failure, CPLE_NotSupported,
170 : "dataCodingFormat=%d is not supported by the S104 driver",
171 : nDataCodingFormat);
172 0 : return nullptr;
173 : }
174 :
175 : // Read additional metadata
176 63 : for (const char *pszAttrName :
177 84 : {"methodWaterLevelProduct", "minDatasetHeight", "maxDatasetHeight"})
178 : {
179 189 : auto poAttr = poWaterLevel->GetAttribute(pszAttrName);
180 63 : if (poAttr)
181 : {
182 42 : const char *pszVal = poAttr->ReadAsString();
183 42 : if (pszVal)
184 : {
185 42 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
186 : }
187 : }
188 : }
189 :
190 21 : int nNumInstances = 1;
191 21 : if (osGroup.empty())
192 : {
193 48 : auto poNumInstances = poWaterLevel->GetAttribute("numInstances");
194 17 : if (poNumInstances &&
195 17 : poNumInstances->GetDataType().GetClass() == GEDTC_NUMERIC)
196 : {
197 1 : nNumInstances = poNumInstances->ReadAsInt();
198 : }
199 : }
200 21 : if (nNumInstances != 1)
201 : {
202 2 : CPLStringList aosSubDSList;
203 1 : int iSubDS = 0;
204 2 : for (const std::string &featureInstanceName :
205 5 : poWaterLevel->GetGroupNames())
206 : {
207 : auto poFeatureInstance =
208 4 : poWaterLevel->OpenGroup(featureInstanceName);
209 2 : if (poFeatureInstance)
210 : {
211 4 : GDALMajorObject mo;
212 : // Read first vertical datum from root group and let the
213 : // coverage override it.
214 2 : S100ReadVerticalDatum(&mo, poRootGroup.get());
215 2 : S100ReadVerticalDatum(&mo, poFeatureInstance.get());
216 :
217 4 : const auto aosGroupNames = poFeatureInstance->GetGroupNames();
218 4 : for (const auto &osSubGroup : aosGroupNames)
219 : {
220 2 : if (auto poSubGroup =
221 4 : poFeatureInstance->OpenGroup(osSubGroup))
222 : {
223 2 : ++iSubDS;
224 : aosSubDSList.SetNameValue(
225 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
226 : CPLSPrintf("S104:\"%s\":%s:%s", osFilename.c_str(),
227 : featureInstanceName.c_str(),
228 2 : osSubGroup.c_str()));
229 :
230 4 : std::string verticalDatum;
231 : const char *pszValue =
232 2 : mo.GetMetadataItem(S100_VERTICAL_DATUM_MEANING);
233 2 : if (pszValue)
234 : {
235 2 : verticalDatum = ", vertical datum ";
236 2 : verticalDatum += pszValue;
237 : pszValue =
238 2 : mo.GetMetadataItem(S100_VERTICAL_DATUM_ABBREV);
239 2 : if (pszValue)
240 : {
241 2 : verticalDatum += " (";
242 2 : verticalDatum += pszValue;
243 2 : verticalDatum += ')';
244 : }
245 : }
246 : else
247 : {
248 : pszValue =
249 0 : mo.GetMetadataItem(S100_VERTICAL_DATUM_NAME);
250 0 : if (pszValue)
251 : {
252 0 : verticalDatum = ", vertical datum ";
253 0 : verticalDatum += pszValue;
254 : }
255 : }
256 :
257 4 : std::string osSubDSDesc;
258 : const auto poTimePoint =
259 6 : poSubGroup->GetAttribute("timePoint");
260 2 : if (poTimePoint)
261 : {
262 2 : const char *pszVal = poTimePoint->ReadAsString();
263 2 : if (pszVal)
264 : {
265 2 : osSubDSDesc = "Values for feature instance ";
266 2 : osSubDSDesc += featureInstanceName;
267 2 : osSubDSDesc += verticalDatum;
268 2 : osSubDSDesc += " at timestamp ";
269 2 : osSubDSDesc += pszVal;
270 : }
271 : }
272 2 : if (osSubDSDesc.empty())
273 : {
274 0 : osSubDSDesc = "Values for feature instance ";
275 0 : osSubDSDesc += featureInstanceName;
276 0 : osSubDSDesc += verticalDatum;
277 0 : osSubDSDesc += " and group ";
278 0 : osSubDSDesc += osSubGroup;
279 : }
280 :
281 : aosSubDSList.SetNameValue(
282 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
283 2 : osSubDSDesc.c_str());
284 : }
285 : }
286 : }
287 : }
288 :
289 1 : poDS->GDALDataset::SetMetadata(aosSubDSList.List(), "SUBDATASETS");
290 :
291 : // Setup/check for pam .aux.xml.
292 1 : poDS->SetDescription(osFilename.c_str());
293 1 : poDS->TryLoadXML();
294 :
295 : // Setup overviews.
296 1 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
297 :
298 1 : return poDS.release();
299 : }
300 :
301 40 : auto poFeatureInstance = poWaterLevel->OpenGroup(osFeatureInstance);
302 20 : if (!poFeatureInstance)
303 : {
304 0 : CPLError(CE_Failure, CPLE_AppDefined,
305 : "Cannot find /WaterLevel/%s group", osFeatureInstance.c_str());
306 0 : return nullptr;
307 : }
308 :
309 : // Read additional metadata
310 100 : for (const char *pszAttrName :
311 : {"timeRecordInterval", "dateTimeOfFirstRecord", "dateTimeOfLastRecord",
312 120 : "numberOfTimes", "dataDynamicity"})
313 : {
314 300 : auto poAttr = poFeatureInstance->GetAttribute(pszAttrName);
315 100 : if (poAttr)
316 : {
317 80 : const char *pszVal = poAttr->ReadAsString();
318 80 : if (pszVal)
319 : {
320 80 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
321 : }
322 : }
323 : }
324 :
325 20 : if (auto poDataDynamicity =
326 60 : poFeatureInstance->GetAttribute("dataDynamicity"))
327 : {
328 0 : if (poDataDynamicity->GetDataType().GetClass() == GEDTC_NUMERIC)
329 : {
330 0 : const int nVal = poDataDynamicity->ReadAsInt();
331 : const std::map<int, const char *> oDataDynamicityMap = {
332 : {1, "Observation"},
333 : {2, "Astronomical prediction"},
334 : {3, "Analysis or hybrid method"},
335 : {5, "Hydrodynamic model forecast"},
336 0 : };
337 0 : const auto oIter = oDataDynamicityMap.find(nVal);
338 0 : if (oIter != oDataDynamicityMap.end())
339 0 : poDS->GDALDataset::SetMetadataItem("DATA_DYNAMICITY_MEANING",
340 0 : oIter->second);
341 : }
342 : }
343 :
344 40 : if (auto poStartSequence = poFeatureInstance->GetAttribute("startSequence"))
345 : {
346 20 : const char *pszStartSequence = poStartSequence->ReadAsString();
347 20 : if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
348 : {
349 0 : CPLError(CE_Failure, CPLE_AppDefined,
350 : "startSequence (=%s) != 0,0 is not supported",
351 : pszStartSequence);
352 0 : return nullptr;
353 : }
354 : }
355 :
356 20 : if (!S100GetNumPointsLongitudinalLatitudinal(
357 20 : poFeatureInstance.get(), poDS->nRasterXSize, poDS->nRasterYSize))
358 : {
359 0 : return nullptr;
360 : }
361 :
362 : // Potentially override vertical datum
363 20 : S100ReadVerticalDatum(poDS.get(), poFeatureInstance.get());
364 :
365 20 : const bool bNorthUp = CPLTestBool(
366 20 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
367 :
368 : // Compute geotransform
369 40 : poDS->m_bHasGT =
370 20 : S100GetGeoTransform(poFeatureInstance.get(), poDS->m_gt, bNorthUp);
371 :
372 20 : if (osGroup.empty())
373 : {
374 30 : const auto aosGroupNames = poFeatureInstance->GetGroupNames();
375 15 : int iSubDS = 1;
376 30 : for (const auto &osSubGroup : aosGroupNames)
377 : {
378 30 : if (auto poSubGroup = poFeatureInstance->OpenGroup(osSubGroup))
379 : {
380 15 : poDS->GDALDataset::SetMetadataItem(
381 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
382 : CPLSPrintf("S104:\"%s\":%s", osFilename.c_str(),
383 : osSubGroup.c_str()),
384 : "SUBDATASETS");
385 30 : std::string osSubDSDesc = "Values for group ";
386 15 : osSubDSDesc += osSubGroup;
387 30 : const auto poTimePoint = poSubGroup->GetAttribute("timePoint");
388 15 : if (poTimePoint)
389 : {
390 15 : const char *pszVal = poTimePoint->ReadAsString();
391 15 : if (pszVal)
392 : {
393 15 : osSubDSDesc = "Values at timestamp ";
394 15 : osSubDSDesc += pszVal;
395 : }
396 : }
397 15 : poDS->GDALDataset::SetMetadataItem(
398 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
399 : osSubDSDesc.c_str(), "SUBDATASETS");
400 15 : ++iSubDS;
401 : }
402 : }
403 : }
404 : else
405 : {
406 5 : auto poGroup = poFeatureInstance->OpenGroup(osGroup);
407 5 : if (!poGroup)
408 : {
409 1 : CPLError(CE_Failure, CPLE_AppDefined,
410 : "Cannot find /WaterLevel/%s/%s group",
411 : osFeatureInstance.c_str(), osGroup.c_str());
412 1 : return nullptr;
413 : }
414 :
415 8 : auto poValuesArray = poGroup->OpenMDArray("values");
416 4 : if (!poValuesArray)
417 : {
418 0 : CPLError(CE_Failure, CPLE_AppDefined,
419 : "Cannot find /WaterLevel/%s/%s/values array",
420 : osFeatureInstance.c_str(), osGroup.c_str());
421 0 : return nullptr;
422 : }
423 :
424 4 : if (poValuesArray->GetDimensionCount() != 2)
425 : {
426 0 : CPLError(CE_Failure, CPLE_AppDefined,
427 : "Wrong dimension count for %s",
428 0 : poValuesArray->GetFullName().c_str());
429 0 : return nullptr;
430 : }
431 :
432 4 : const auto &oType = poValuesArray->GetDataType();
433 4 : if (oType.GetClass() != GEDTC_COMPOUND)
434 : {
435 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
436 0 : poValuesArray->GetFullName().c_str());
437 0 : return nullptr;
438 : }
439 :
440 4 : const auto &oComponents = oType.GetComponents();
441 4 : if (oComponents.size() != 2 ||
442 8 : oComponents[0]->GetName() != "waterLevelHeight" ||
443 4 : oComponents[0]->GetType().GetNumericDataType() != GDT_Float32 ||
444 12 : oComponents[1]->GetName() != "waterLevelTrend" ||
445 4 : (oComponents[1]->GetType().GetNumericDataType() != GDT_Byte &&
446 : // In theory should be Byte, but 104US00_ches_dcf2_20190606T12Z.h5 uses Int32
447 0 : oComponents[1]->GetType().GetNumericDataType() != GDT_Int32))
448 : {
449 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
450 0 : poValuesArray->GetFullName().c_str());
451 0 : return nullptr;
452 : }
453 :
454 4 : const auto &apoDims = poValuesArray->GetDimensions();
455 4 : if (apoDims[0]->GetSize() != static_cast<unsigned>(poDS->nRasterYSize))
456 : {
457 0 : CPLError(CE_Failure, CPLE_AppDefined,
458 : "numPointsLatitudinal(=%d) doesn't match first dimension "
459 : "size of %s (=%d)",
460 0 : poDS->nRasterYSize, poValuesArray->GetFullName().c_str(),
461 0 : static_cast<int>(apoDims[0]->GetSize()));
462 0 : return nullptr;
463 : }
464 4 : if (apoDims[1]->GetSize() != static_cast<unsigned>(poDS->nRasterXSize))
465 : {
466 0 : CPLError(CE_Failure, CPLE_AppDefined,
467 : "numPointsLongitudinal(=%d) doesn't match second "
468 : "dimension size of %s (=%d)",
469 0 : poDS->nRasterXSize, poValuesArray->GetFullName().c_str(),
470 0 : static_cast<int>(apoDims[1]->GetSize()));
471 0 : return nullptr;
472 : }
473 :
474 4 : if (bNorthUp)
475 3 : poValuesArray = poValuesArray->GetView("[::-1,...]");
476 :
477 : // Create waterLevelHeight band
478 : auto poWaterLevelHeight =
479 12 : poValuesArray->GetView("[\"waterLevelHeight\"]");
480 : auto poWaterLevelHeightDS = std::unique_ptr<GDALDataset>(
481 8 : poWaterLevelHeight->AsClassicDataset(1, 0));
482 : auto poWaterLevelHeightBand =
483 8 : std::make_unique<S104RasterBand>(std::move(poWaterLevelHeightDS));
484 4 : poWaterLevelHeightBand->SetDescription("waterLevelHeight");
485 4 : poWaterLevelHeightBand->m_osUnitType = "metre";
486 4 : poDS->SetBand(1, poWaterLevelHeightBand.release());
487 :
488 : // Create waterLevelTrend band
489 : auto poWaterLevelTrend =
490 12 : poValuesArray->GetView("[\"waterLevelTrend\"]");
491 : auto poWaterLevelTrendDS = std::unique_ptr<GDALDataset>(
492 8 : poWaterLevelTrend->AsClassicDataset(1, 0));
493 : auto poWaterLevelTrendBand =
494 8 : std::make_unique<S104RasterBand>(std::move(poWaterLevelTrendDS));
495 4 : poWaterLevelTrendBand->SetDescription("waterLevelTrend");
496 :
497 : // From D-5.3 Water Level Trend of S-101 v1.1 spec
498 8 : auto poRAT = std::make_unique<GDALDefaultRasterAttributeTable>();
499 4 : poRAT->CreateColumn("code", GFT_Integer, GFU_MinMax);
500 4 : poRAT->CreateColumn("label", GFT_String, GFU_Generic);
501 4 : poRAT->CreateColumn("definition", GFT_String, GFU_Generic);
502 :
503 : const struct
504 : {
505 : int nCode;
506 : const char *pszLabel;
507 : const char *pszDefinition;
508 4 : } aoRatValues[] = {
509 : {0, "Nodata", "No data"},
510 : {1, "Decreasing", "Becoming smaller in magnitude"},
511 : {2, "Increasing", "Becoming larger in magnitude"},
512 : {3, "Steady", "Constant"},
513 : };
514 :
515 4 : int iRow = 0;
516 20 : for (const auto &oRecord : aoRatValues)
517 : {
518 16 : int iCol = 0;
519 16 : poRAT->SetValue(iRow, iCol++, oRecord.nCode);
520 16 : poRAT->SetValue(iRow, iCol++, oRecord.pszLabel);
521 16 : poRAT->SetValue(iRow, iCol++, oRecord.pszDefinition);
522 16 : ++iRow;
523 : }
524 :
525 4 : poWaterLevelTrendBand->m_poRAT = std::move(poRAT);
526 :
527 4 : poDS->SetBand(2, poWaterLevelTrendBand.release());
528 : }
529 :
530 19 : poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
531 :
532 : // Setup/check for pam .aux.xml.
533 19 : poDS->SetDescription(osFilename.c_str());
534 19 : poDS->TryLoadXML();
535 :
536 : // Setup overviews.
537 19 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
538 :
539 19 : return poDS.release();
540 : }
541 :
542 : /************************************************************************/
543 : /* S104DatasetDriverUnload() */
544 : /************************************************************************/
545 :
546 323 : static void S104DatasetDriverUnload(GDALDriver *)
547 : {
548 323 : HDF5UnloadFileDriver();
549 323 : }
550 :
551 : /************************************************************************/
552 : /* GDALRegister_S104() */
553 : /************************************************************************/
554 355 : void GDALRegister_S104()
555 :
556 : {
557 355 : if (!GDAL_CHECK_VERSION("S104"))
558 0 : return;
559 :
560 355 : if (GDALGetDriverByName(S104_DRIVER_NAME) != nullptr)
561 0 : return;
562 :
563 355 : GDALDriver *poDriver = new GDALDriver();
564 :
565 355 : S104DriverSetCommonMetadata(poDriver);
566 355 : poDriver->pfnOpen = S104Dataset::Open;
567 355 : poDriver->pfnUnloadDriver = S104DatasetDriverUnload;
568 :
569 355 : GetGDALDriverManager()->RegisterDriver(poDriver);
570 : }
|
