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_priv.h"
21 : #include "gdal_proxy.h"
22 : #include "gdal_rat.h"
23 :
24 : #include <limits>
25 :
26 : /************************************************************************/
27 : /* S104Dataset */
28 : /************************************************************************/
29 :
30 8 : class S104Dataset final : public S100BaseDataset
31 : {
32 : public:
33 4 : explicit S104Dataset(const std::string &osFilename)
34 4 : : S100BaseDataset(osFilename)
35 : {
36 4 : }
37 :
38 : ~S104Dataset() override;
39 :
40 : static GDALDataset *Open(GDALOpenInfo *);
41 : };
42 :
43 : S104Dataset::~S104Dataset() = default;
44 :
45 : /************************************************************************/
46 : /* S104RasterBand */
47 : /************************************************************************/
48 :
49 : class S104RasterBand final : public GDALProxyRasterBand
50 : {
51 : friend class S104Dataset;
52 : std::unique_ptr<GDALDataset> m_poDS{};
53 : GDALRasterBand *m_poUnderlyingBand = nullptr;
54 : std::string m_osUnitType{};
55 : std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
56 :
57 : CPL_DISALLOW_COPY_ASSIGN(S104RasterBand)
58 :
59 : public:
60 4 : explicit S104RasterBand(std::unique_ptr<GDALDataset> &&poDSIn)
61 8 : : m_poDS(std::move(poDSIn)),
62 4 : m_poUnderlyingBand(m_poDS->GetRasterBand(1))
63 : {
64 4 : eDataType = m_poUnderlyingBand->GetRasterDataType();
65 4 : m_poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
66 4 : }
67 :
68 : GDALRasterBand *
69 : RefUnderlyingRasterBand(bool /*bForceOpen*/ = true) const override;
70 :
71 2 : const char *GetUnitType() override
72 : {
73 2 : return m_osUnitType.c_str();
74 : }
75 :
76 1 : GDALRasterAttributeTable *GetDefaultRAT() override
77 : {
78 1 : return m_poRAT.get();
79 : }
80 : };
81 :
82 : GDALRasterBand *
83 8 : S104RasterBand::RefUnderlyingRasterBand(bool /*bForceOpen*/) const
84 : {
85 8 : return m_poUnderlyingBand;
86 : }
87 :
88 : /************************************************************************/
89 : /* Open() */
90 : /************************************************************************/
91 :
92 5 : GDALDataset *S104Dataset::Open(GDALOpenInfo *poOpenInfo)
93 :
94 : {
95 : // Confirm that this appears to be a S104 file.
96 5 : if (!S104DatasetIdentify(poOpenInfo))
97 0 : return nullptr;
98 :
99 : HDF5_GLOBAL_LOCK();
100 :
101 5 : if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
102 : {
103 1 : return HDF5Dataset::OpenMultiDim(poOpenInfo);
104 : }
105 :
106 : // Confirm the requested access is supported.
107 4 : if (poOpenInfo->eAccess == GA_Update)
108 : {
109 0 : ReportUpdateNotSupportedByDriver("S104");
110 0 : return nullptr;
111 : }
112 :
113 8 : std::string osFilename(poOpenInfo->pszFilename);
114 8 : std::string osGroup;
115 4 : if (STARTS_WITH(poOpenInfo->pszFilename, "S104:"))
116 : {
117 : const CPLStringList aosTokens(
118 3 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
119 3 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
120 :
121 3 : if (aosTokens.size() == 2)
122 : {
123 0 : osFilename = aosTokens[1];
124 : }
125 3 : else if (aosTokens.size() == 3)
126 : {
127 3 : osFilename = aosTokens[1];
128 3 : osGroup = aosTokens[2];
129 : }
130 : else
131 : {
132 0 : return nullptr;
133 : }
134 : }
135 :
136 8 : auto poDS = std::make_unique<S104Dataset>(osFilename);
137 4 : if (!poDS->Init())
138 0 : return nullptr;
139 :
140 4 : const auto &poRootGroup = poDS->m_poRootGroup;
141 :
142 12 : auto poWaterLevel = poRootGroup->OpenGroup("WaterLevel");
143 4 : if (!poWaterLevel)
144 : {
145 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find /WaterLevel group");
146 0 : return nullptr;
147 : }
148 :
149 12 : auto poDataCodingFormat = poWaterLevel->GetAttribute("dataCodingFormat");
150 8 : if (!poDataCodingFormat ||
151 4 : poDataCodingFormat->GetDataType().GetClass() != GEDTC_NUMERIC)
152 : {
153 0 : CPLError(CE_Failure, CPLE_AppDefined,
154 : "Cannot find /WaterLevel/dataCodingFormat attribute");
155 0 : return nullptr;
156 : }
157 4 : const int nDataCodingFormat = poDataCodingFormat->ReadAsInt();
158 4 : if (nDataCodingFormat != 2)
159 : {
160 0 : CPLError(CE_Failure, CPLE_NotSupported,
161 : "dataCodingFormat=%d is not supported by the S104 driver",
162 : nDataCodingFormat);
163 0 : return nullptr;
164 : }
165 :
166 : // Read additional metadata
167 12 : for (const char *pszAttrName :
168 16 : {"methodWaterLevelProduct", "minDatasetHeight", "maxDatasetHeight"})
169 : {
170 36 : auto poAttr = poWaterLevel->GetAttribute(pszAttrName);
171 12 : if (poAttr)
172 : {
173 8 : const char *pszVal = poAttr->ReadAsString();
174 8 : if (pszVal)
175 : {
176 8 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
177 : }
178 : }
179 : }
180 :
181 12 : auto poWaterLevel01 = poWaterLevel->OpenGroup("WaterLevel.01");
182 4 : if (!poWaterLevel01)
183 : {
184 0 : CPLError(CE_Failure, CPLE_AppDefined,
185 : "Cannot find /WaterLevel/WaterLevel.01 group");
186 0 : return nullptr;
187 : }
188 :
189 : // Read additional metadata
190 16 : for (const char *pszAttrName :
191 : {"timeRecordInterval", "dateTimeOfFirstRecord", "dateTimeOfLastRecord",
192 20 : "numberOfTimes"})
193 : {
194 48 : auto poAttr = poWaterLevel01->GetAttribute(pszAttrName);
195 16 : if (poAttr)
196 : {
197 16 : const char *pszVal = poAttr->ReadAsString();
198 16 : if (pszVal)
199 : {
200 16 : poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
201 : }
202 : }
203 : }
204 :
205 8 : if (auto poStartSequence = poWaterLevel01->GetAttribute("startSequence"))
206 : {
207 4 : const char *pszStartSequence = poStartSequence->ReadAsString();
208 4 : if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
209 : {
210 0 : CPLError(CE_Failure, CPLE_AppDefined,
211 : "startSequence (=%s) != 0,0 is not supported",
212 : pszStartSequence);
213 0 : return nullptr;
214 : }
215 : }
216 :
217 4 : if (!S100GetNumPointsLongitudinalLatitudinal(
218 4 : poWaterLevel01.get(), poDS->nRasterXSize, poDS->nRasterYSize))
219 : {
220 0 : return nullptr;
221 : }
222 :
223 4 : const bool bNorthUp = CPLTestBool(
224 4 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
225 :
226 : // Compute geotransform
227 8 : poDS->m_bHasGT =
228 4 : S100GetGeoTransform(poWaterLevel01.get(), poDS->m_gt, bNorthUp);
229 :
230 4 : if (osGroup.empty())
231 : {
232 2 : const auto aosGroupNames = poWaterLevel01->GetGroupNames();
233 1 : int iSubDS = 1;
234 2 : for (const auto &osSubGroup : aosGroupNames)
235 : {
236 2 : if (auto poSubGroup = poWaterLevel01->OpenGroup(osSubGroup))
237 : {
238 1 : poDS->GDALDataset::SetMetadataItem(
239 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
240 : CPLSPrintf("S104:\"%s\":%s", osFilename.c_str(),
241 : osSubGroup.c_str()),
242 : "SUBDATASETS");
243 2 : std::string osSubDSDesc = "Values for group ";
244 1 : osSubDSDesc += osSubGroup;
245 2 : const auto poTimePoint = poSubGroup->GetAttribute("timePoint");
246 1 : if (poTimePoint)
247 : {
248 1 : const char *pszVal = poTimePoint->ReadAsString();
249 1 : if (pszVal)
250 : {
251 1 : osSubDSDesc = "Values at timestamp ";
252 1 : osSubDSDesc += pszVal;
253 : }
254 : }
255 1 : poDS->GDALDataset::SetMetadataItem(
256 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
257 : osSubDSDesc.c_str(), "SUBDATASETS");
258 1 : ++iSubDS;
259 : }
260 : }
261 : }
262 : else
263 : {
264 3 : auto poGroup = poWaterLevel01->OpenGroup(osGroup);
265 3 : if (!poGroup)
266 : {
267 1 : CPLError(CE_Failure, CPLE_AppDefined,
268 : "Cannot find /WaterLevel/WaterLevel.01/%s group",
269 : osGroup.c_str());
270 1 : return nullptr;
271 : }
272 :
273 4 : auto poValuesArray = poGroup->OpenMDArray("values");
274 2 : if (!poValuesArray)
275 : {
276 0 : CPLError(CE_Failure, CPLE_AppDefined,
277 : "Cannot find /WaterLevel/WaterLevel.01/%s/values array",
278 : osGroup.c_str());
279 0 : return nullptr;
280 : }
281 :
282 2 : if (poValuesArray->GetDimensionCount() != 2)
283 : {
284 0 : CPLError(CE_Failure, CPLE_AppDefined,
285 : "Wrong dimension count for %s",
286 0 : poValuesArray->GetFullName().c_str());
287 0 : return nullptr;
288 : }
289 :
290 2 : const auto &oType = poValuesArray->GetDataType();
291 2 : if (oType.GetClass() != GEDTC_COMPOUND)
292 : {
293 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
294 0 : poValuesArray->GetFullName().c_str());
295 0 : return nullptr;
296 : }
297 :
298 2 : const auto &oComponents = oType.GetComponents();
299 2 : if (oComponents.size() != 2 ||
300 4 : oComponents[0]->GetName() != "waterLevelHeight" ||
301 2 : oComponents[0]->GetType().GetNumericDataType() != GDT_Float32 ||
302 6 : oComponents[1]->GetName() != "waterLevelTrend" ||
303 2 : (oComponents[1]->GetType().GetNumericDataType() != GDT_Byte &&
304 : // In theory should be Byte, but 104US00_ches_dcf2_20190606T12Z.h5 uses Int32
305 0 : oComponents[1]->GetType().GetNumericDataType() != GDT_Int32))
306 : {
307 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
308 0 : poValuesArray->GetFullName().c_str());
309 0 : return nullptr;
310 : }
311 :
312 2 : const auto &apoDims = poValuesArray->GetDimensions();
313 2 : if (apoDims[0]->GetSize() != static_cast<unsigned>(poDS->nRasterYSize))
314 : {
315 0 : CPLError(CE_Failure, CPLE_AppDefined,
316 : "numPointsLatitudinal(=%d) doesn't match first dimension "
317 : "size of %s (=%d)",
318 0 : poDS->nRasterYSize, poValuesArray->GetFullName().c_str(),
319 0 : static_cast<int>(apoDims[0]->GetSize()));
320 0 : return nullptr;
321 : }
322 2 : if (apoDims[1]->GetSize() != static_cast<unsigned>(poDS->nRasterXSize))
323 : {
324 0 : CPLError(CE_Failure, CPLE_AppDefined,
325 : "numPointsLongitudinal(=%d) doesn't match second "
326 : "dimension size of %s (=%d)",
327 0 : poDS->nRasterXSize, poValuesArray->GetFullName().c_str(),
328 0 : static_cast<int>(apoDims[1]->GetSize()));
329 0 : return nullptr;
330 : }
331 :
332 2 : if (bNorthUp)
333 1 : poValuesArray = poValuesArray->GetView("[::-1,...]");
334 :
335 : // Create waterLevelHeight band
336 : auto poWaterLevelHeight =
337 6 : poValuesArray->GetView("[\"waterLevelHeight\"]");
338 : auto poWaterLevelHeightDS = std::unique_ptr<GDALDataset>(
339 4 : poWaterLevelHeight->AsClassicDataset(1, 0));
340 : auto poWaterLevelHeightBand =
341 4 : std::make_unique<S104RasterBand>(std::move(poWaterLevelHeightDS));
342 2 : poWaterLevelHeightBand->SetDescription("waterLevelHeight");
343 2 : poWaterLevelHeightBand->m_osUnitType = "metre";
344 2 : poDS->SetBand(1, poWaterLevelHeightBand.release());
345 :
346 : // Create waterLevelTrend band
347 : auto poWaterLevelTrend =
348 6 : poValuesArray->GetView("[\"waterLevelTrend\"]");
349 : auto poWaterLevelTrendDS = std::unique_ptr<GDALDataset>(
350 4 : poWaterLevelTrend->AsClassicDataset(1, 0));
351 : auto poWaterLevelTrendBand =
352 4 : std::make_unique<S104RasterBand>(std::move(poWaterLevelTrendDS));
353 2 : poWaterLevelTrendBand->SetDescription("waterLevelTrend");
354 :
355 : // From D-5.3 Water Level Trend of S-101 v1.1 spec
356 4 : auto poRAT = std::make_unique<GDALDefaultRasterAttributeTable>();
357 2 : poRAT->CreateColumn("code", GFT_Integer, GFU_MinMax);
358 2 : poRAT->CreateColumn("label", GFT_String, GFU_Generic);
359 2 : poRAT->CreateColumn("definition", GFT_String, GFU_Generic);
360 :
361 : const struct
362 : {
363 : int nCode;
364 : const char *pszLabel;
365 : const char *pszDefinition;
366 2 : } aoRatValues[] = {
367 : {0, "Nodata", "No data"},
368 : {1, "Decreasing", "Becoming smaller in magnitude"},
369 : {2, "Increasing", "Becoming larger in magnitude"},
370 : {3, "Steady", "Constant"},
371 : };
372 :
373 2 : int iRow = 0;
374 10 : for (const auto &oRecord : aoRatValues)
375 : {
376 8 : int iCol = 0;
377 8 : poRAT->SetValue(iRow, iCol++, oRecord.nCode);
378 8 : poRAT->SetValue(iRow, iCol++, oRecord.pszLabel);
379 8 : poRAT->SetValue(iRow, iCol++, oRecord.pszDefinition);
380 8 : ++iRow;
381 : }
382 :
383 2 : poWaterLevelTrendBand->m_poRAT = std::move(poRAT);
384 :
385 2 : poDS->SetBand(2, poWaterLevelTrendBand.release());
386 : }
387 :
388 3 : poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
389 :
390 : // Setup/check for pam .aux.xml.
391 3 : poDS->SetDescription(osFilename.c_str());
392 3 : poDS->TryLoadXML();
393 :
394 : // Setup overviews.
395 3 : poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
396 :
397 3 : return poDS.release();
398 : }
399 :
400 : /************************************************************************/
401 : /* S104DatasetDriverUnload() */
402 : /************************************************************************/
403 :
404 6 : static void S104DatasetDriverUnload(GDALDriver *)
405 : {
406 6 : HDF5UnloadFileDriver();
407 6 : }
408 :
409 : /************************************************************************/
410 : /* GDALRegister_S104() */
411 : /************************************************************************/
412 11 : void GDALRegister_S104()
413 :
414 : {
415 11 : if (!GDAL_CHECK_VERSION("S104"))
416 0 : return;
417 :
418 11 : if (GDALGetDriverByName(S104_DRIVER_NAME) != nullptr)
419 0 : return;
420 :
421 11 : GDALDriver *poDriver = new GDALDriver();
422 :
423 11 : S104DriverSetCommonMetadata(poDriver);
424 11 : poDriver->pfnOpen = S104Dataset::Open;
425 11 : poDriver->pfnUnloadDriver = S104DatasetDriverUnload;
426 :
427 11 : GetGDALDriverManager()->RegisterDriver(poDriver);
428 : }
|