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