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             : }