Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Hierarchical Data Format Release 5 (HDF5)
       4             :  * Purpose:  Read S111 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             : /*                             S111Dataset                              */
      28             : /************************************************************************/
      29             : 
      30           8 : class S111Dataset final : public S100BaseDataset
      31             : {
      32             :   public:
      33           4 :     explicit S111Dataset(const std::string &osFilename)
      34           4 :         : S100BaseDataset(osFilename)
      35             :     {
      36           4 :     }
      37             : 
      38             :     ~S111Dataset() override;
      39             : 
      40             :     static GDALDataset *Open(GDALOpenInfo *);
      41             : };
      42             : 
      43             : S111Dataset::~S111Dataset() = default;
      44             : 
      45             : /************************************************************************/
      46             : /*                            S111RasterBand                            */
      47             : /************************************************************************/
      48             : 
      49             : class S111RasterBand final : public GDALProxyRasterBand
      50             : {
      51             :     friend class S111Dataset;
      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(S111RasterBand)
      58             : 
      59             :   public:
      60           4 :     explicit S111RasterBand(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           4 :     const char *GetUnitType() override
      72             :     {
      73           4 :         return m_osUnitType.c_str();
      74             :     }
      75             : 
      76           1 :     GDALRasterAttributeTable *GetDefaultRAT() override
      77             :     {
      78           1 :         return m_poRAT.get();
      79             :     }
      80             : 
      81           0 :     char **GetMetadata(const char *pszDomain) override
      82             :     {
      83             :         // Short-circuit GDALProxyRasterBand...
      84           0 :         return GDALRasterBand::GetMetadata(pszDomain);
      85             :     }
      86             : };
      87             : 
      88             : GDALRasterBand *
      89           8 : S111RasterBand::RefUnderlyingRasterBand(bool /*bForceOpen*/) const
      90             : {
      91           8 :     return m_poUnderlyingBand;
      92             : }
      93             : 
      94             : /************************************************************************/
      95             : /*                                Open()                                */
      96             : /************************************************************************/
      97             : 
      98           5 : GDALDataset *S111Dataset::Open(GDALOpenInfo *poOpenInfo)
      99             : 
     100             : {
     101             :     // Confirm that this appears to be a S111 file.
     102           5 :     if (!S111DatasetIdentify(poOpenInfo))
     103           0 :         return nullptr;
     104             : 
     105             :     HDF5_GLOBAL_LOCK();
     106             : 
     107           5 :     if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
     108             :     {
     109           1 :         return HDF5Dataset::OpenMultiDim(poOpenInfo);
     110             :     }
     111             : 
     112             :     // Confirm the requested access is supported.
     113           4 :     if (poOpenInfo->eAccess == GA_Update)
     114             :     {
     115           0 :         ReportUpdateNotSupportedByDriver("S111");
     116           0 :         return nullptr;
     117             :     }
     118             : 
     119           8 :     std::string osFilename(poOpenInfo->pszFilename);
     120           8 :     std::string osGroup;
     121           4 :     if (STARTS_WITH(poOpenInfo->pszFilename, "S111:"))
     122             :     {
     123             :         const CPLStringList aosTokens(
     124           3 :             CSLTokenizeString2(poOpenInfo->pszFilename, ":",
     125           3 :                                CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
     126             : 
     127           3 :         if (aosTokens.size() == 2)
     128             :         {
     129           0 :             osFilename = aosTokens[1];
     130             :         }
     131           3 :         else if (aosTokens.size() == 3)
     132             :         {
     133           3 :             osFilename = aosTokens[1];
     134           3 :             osGroup = aosTokens[2];
     135             :         }
     136             :         else
     137             :         {
     138           0 :             return nullptr;
     139             :         }
     140             :     }
     141             : 
     142           8 :     auto poDS = std::make_unique<S111Dataset>(osFilename);
     143           4 :     if (!poDS->Init())
     144           0 :         return nullptr;
     145             : 
     146           4 :     const auto &poRootGroup = poDS->m_poRootGroup;
     147             : 
     148          12 :     auto poSurfaceCurrent = poRootGroup->OpenGroup("SurfaceCurrent");
     149           4 :     if (!poSurfaceCurrent)
     150             :     {
     151           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     152             :                  "Cannot find /SurfaceCurrent group");
     153           0 :         return nullptr;
     154             :     }
     155             : 
     156             :     auto poDataCodingFormat =
     157          12 :         poSurfaceCurrent->GetAttribute("dataCodingFormat");
     158           8 :     if (!poDataCodingFormat ||
     159           4 :         poDataCodingFormat->GetDataType().GetClass() != GEDTC_NUMERIC)
     160             :     {
     161           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     162             :                  "Cannot find /SurfaceCurrent/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 S111 driver",
     170             :                  nDataCodingFormat);
     171           0 :         return nullptr;
     172             :     }
     173             : 
     174             :     // Read additional metadata
     175          12 :     for (const char *pszAttrName :
     176             :          {"methodCurrentsProduct", "minDatasetCurrentSpeed",
     177          16 :           "maxDatasetCurrentSpeed"})
     178             :     {
     179          36 :         auto poAttr = poSurfaceCurrent->GetAttribute(pszAttrName);
     180          12 :         if (poAttr)
     181             :         {
     182           8 :             const char *pszVal = poAttr->ReadAsString();
     183           8 :             if (pszVal)
     184             :             {
     185           8 :                 poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
     186             :             }
     187             :         }
     188             :     }
     189             : 
     190          12 :     auto poSurfaceCurrent01 = poSurfaceCurrent->OpenGroup("SurfaceCurrent.01");
     191           4 :     if (!poSurfaceCurrent01)
     192             :     {
     193           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     194             :                  "Cannot find /SurfaceCurrent/SurfaceCurrent.01 group");
     195           0 :         return nullptr;
     196             :     }
     197             : 
     198             :     // Read additional metadata
     199          16 :     for (const char *pszAttrName :
     200             :          {"timeRecordInterval", "dateTimeOfFirstRecord", "dateTimeOfLastRecord",
     201          20 :           "numberOfTimes"})
     202             :     {
     203          48 :         auto poAttr = poSurfaceCurrent01->GetAttribute(pszAttrName);
     204          16 :         if (poAttr)
     205             :         {
     206          16 :             const char *pszVal = poAttr->ReadAsString();
     207          16 :             if (pszVal)
     208             :             {
     209          16 :                 poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
     210             :             }
     211             :         }
     212             :     }
     213             : 
     214           4 :     if (auto poStartSequence =
     215           8 :             poSurfaceCurrent01->GetAttribute("startSequence"))
     216             :     {
     217           4 :         const char *pszStartSequence = poStartSequence->ReadAsString();
     218           4 :         if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
     219             :         {
     220           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     221             :                      "startSequence (=%s) != 0,0 is not supported",
     222             :                      pszStartSequence);
     223           0 :             return nullptr;
     224             :         }
     225             :     }
     226             : 
     227           4 :     if (!S100GetNumPointsLongitudinalLatitudinal(
     228           4 :             poSurfaceCurrent01.get(), poDS->nRasterXSize, poDS->nRasterYSize))
     229             :     {
     230           0 :         return nullptr;
     231             :     }
     232             : 
     233           4 :     const bool bNorthUp = CPLTestBool(
     234           4 :         CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
     235             : 
     236             :     // Compute geotransform
     237           8 :     poDS->m_bHasGT =
     238           4 :         S100GetGeoTransform(poSurfaceCurrent01.get(), poDS->m_gt, bNorthUp);
     239             : 
     240           4 :     if (osGroup.empty())
     241             :     {
     242           2 :         const auto aosGroupNames = poSurfaceCurrent01->GetGroupNames();
     243           1 :         int iSubDS = 1;
     244           2 :         for (const auto &osSubGroup : aosGroupNames)
     245             :         {
     246           2 :             if (auto poSubGroup = poSurfaceCurrent01->OpenGroup(osSubGroup))
     247             :             {
     248           1 :                 poDS->GDALDataset::SetMetadataItem(
     249             :                     CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
     250             :                     CPLSPrintf("S111:\"%s\":%s", osFilename.c_str(),
     251             :                                osSubGroup.c_str()),
     252             :                     "SUBDATASETS");
     253           2 :                 std::string osSubDSDesc = "Values for group ";
     254           1 :                 osSubDSDesc += osSubGroup;
     255           2 :                 const auto poTimePoint = poSubGroup->GetAttribute("timePoint");
     256           1 :                 if (poTimePoint)
     257             :                 {
     258           1 :                     const char *pszVal = poTimePoint->ReadAsString();
     259           1 :                     if (pszVal)
     260             :                     {
     261           1 :                         osSubDSDesc = "Values at timestamp ";
     262           1 :                         osSubDSDesc += pszVal;
     263             :                     }
     264             :                 }
     265           1 :                 poDS->GDALDataset::SetMetadataItem(
     266             :                     CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
     267             :                     osSubDSDesc.c_str(), "SUBDATASETS");
     268           1 :                 ++iSubDS;
     269             :             }
     270             :         }
     271             :     }
     272             :     else
     273             :     {
     274           3 :         auto poGroup = poSurfaceCurrent01->OpenGroup(osGroup);
     275           3 :         if (!poGroup)
     276             :         {
     277           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     278             :                      "Cannot find /SurfaceCurrent/SurfaceCurrent.01/%s group",
     279             :                      osGroup.c_str());
     280           1 :             return nullptr;
     281             :         }
     282             : 
     283           4 :         auto poValuesArray = poGroup->OpenMDArray("values");
     284           2 :         if (!poValuesArray)
     285             :         {
     286           0 :             CPLError(
     287             :                 CE_Failure, CPLE_AppDefined,
     288             :                 "Cannot find /SurfaceCurrent/SurfaceCurrent.01/%s/values array",
     289             :                 osGroup.c_str());
     290           0 :             return nullptr;
     291             :         }
     292             : 
     293           2 :         if (poValuesArray->GetDimensionCount() != 2)
     294             :         {
     295           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     296             :                      "Wrong dimension count for %s",
     297           0 :                      poValuesArray->GetFullName().c_str());
     298           0 :             return nullptr;
     299             :         }
     300             : 
     301           2 :         const auto &oType = poValuesArray->GetDataType();
     302           2 :         if (oType.GetClass() != GEDTC_COMPOUND)
     303             :         {
     304           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
     305           0 :                      poValuesArray->GetFullName().c_str());
     306           0 :             return nullptr;
     307             :         }
     308             : 
     309           2 :         const auto &oComponents = oType.GetComponents();
     310           4 :         if (!(oComponents.size() == 2 &&
     311           4 :               ((oComponents[0]->GetName() == "surfaceCurrentSpeed" &&
     312           2 :                 oComponents[0]->GetType().GetNumericDataType() == GDT_Float32 &&
     313           4 :                 oComponents[1]->GetName() == "surfaceCurrentDirection" &&
     314           2 :                 oComponents[1]->GetType().GetNumericDataType() ==
     315           0 :                     GDT_Float32) ||
     316             :                // S111US_20170829.0100_W078.N44_F2_loofs_type2.h5 has direction first...
     317           0 :                (oComponents[0]->GetName() == "surfaceCurrentDirection" &&
     318           0 :                 oComponents[0]->GetType().GetNumericDataType() == GDT_Float32 &&
     319           0 :                 oComponents[1]->GetName() == "surfaceCurrentSpeed" &&
     320           0 :                 oComponents[1]->GetType().GetNumericDataType() ==
     321             :                     GDT_Float32))))
     322             :         {
     323           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
     324           0 :                      poValuesArray->GetFullName().c_str());
     325           0 :             return nullptr;
     326             :         }
     327             : 
     328           2 :         const auto &apoDims = poValuesArray->GetDimensions();
     329           2 :         if (apoDims[0]->GetSize() != static_cast<unsigned>(poDS->nRasterYSize))
     330             :         {
     331           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     332             :                      "numPointsLatitudinal(=%d) doesn't match first dimension "
     333             :                      "size of %s (=%d)",
     334           0 :                      poDS->nRasterYSize, poValuesArray->GetFullName().c_str(),
     335           0 :                      static_cast<int>(apoDims[0]->GetSize()));
     336           0 :             return nullptr;
     337             :         }
     338           2 :         if (apoDims[1]->GetSize() != static_cast<unsigned>(poDS->nRasterXSize))
     339             :         {
     340           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     341             :                      "numPointsLongitudinal(=%d) doesn't match second "
     342             :                      "dimension size of %s (=%d)",
     343           0 :                      poDS->nRasterXSize, poValuesArray->GetFullName().c_str(),
     344           0 :                      static_cast<int>(apoDims[1]->GetSize()));
     345           0 :             return nullptr;
     346             :         }
     347             : 
     348           2 :         if (bNorthUp)
     349           1 :             poValuesArray = poValuesArray->GetView("[::-1,...]");
     350             : 
     351             :         // Create surfaceCurrentSpeed band
     352             :         auto poSurfaceCurrentSpeed =
     353           6 :             poValuesArray->GetView("[\"surfaceCurrentSpeed\"]");
     354             :         auto poSurfaceCurrentSpeedDS = std::unique_ptr<GDALDataset>(
     355           4 :             poSurfaceCurrentSpeed->AsClassicDataset(1, 0));
     356             :         auto poSurfaceCurrentSpeedBand = std::make_unique<S111RasterBand>(
     357           4 :             std::move(poSurfaceCurrentSpeedDS));
     358           2 :         poSurfaceCurrentSpeedBand->SetDescription("surfaceCurrentSpeed");
     359           2 :         poSurfaceCurrentSpeedBand->m_osUnitType = "knots";
     360             : 
     361             :         // From S-111 v1.2 table 9.1 (Speed ranges) and 9.2 (Colour schemas)
     362           4 :         auto poRAT = std::make_unique<GDALDefaultRasterAttributeTable>();
     363           2 :         poRAT->CreateColumn("speed_band", GFT_Integer, GFU_Generic);
     364           2 :         poRAT->CreateColumn("min_speed", GFT_Real, GFU_Min);
     365           2 :         poRAT->CreateColumn("width_band", GFT_Real, GFU_Generic);
     366           2 :         poRAT->CreateColumn("color", GFT_String, GFU_Generic);
     367           2 :         poRAT->CreateColumn("red", GFT_Integer, GFU_RedMin);
     368           2 :         poRAT->CreateColumn("green", GFT_Integer, GFU_GreenMin);
     369           2 :         poRAT->CreateColumn("blue", GFT_Integer, GFU_BlueMin);
     370             : 
     371             :         const struct
     372             :         {
     373             :             int nSpeedBand;
     374             :             double dfMinSpeed;
     375             :             double dfWidthBand;
     376             :             const char *pszColor;
     377             :             int nRed;
     378             :             int nGreen;
     379             :             int nBlue;
     380           2 :         } aoRatValues[] = {
     381             :             {1, 0.0, 0.5, "purple", 118, 82, 226},
     382             :             {2, 0.5, 0.5, "dark blue", 72, 152, 211},
     383             :             {3, 1.0, 1.0, "light blue", 97, 203, 229},
     384             :             {4, 2.0, 1.0, "dark green", 109, 188, 69},
     385             :             {5, 3.0, 2.0, "light green", 180, 220, 0},
     386             :             {6, 5.0, 2.0, "yellow green", 205, 193, 0},
     387             :             {7, 7.0, 3.0, "orange", 248, 167, 24},
     388             :             {8, 10.0, 3.0, "pink", 247, 162, 157},
     389             :             {9, 13.0, 86.0, "red", 255, 30, 30},
     390             :         };
     391             : 
     392           2 :         int iRow = 0;
     393          20 :         for (const auto &oRecord : aoRatValues)
     394             :         {
     395          18 :             int iCol = 0;
     396          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nSpeedBand);
     397          18 :             poRAT->SetValue(iRow, iCol++, oRecord.dfMinSpeed);
     398          18 :             poRAT->SetValue(iRow, iCol++, oRecord.dfWidthBand);
     399          18 :             poRAT->SetValue(iRow, iCol++, oRecord.pszColor);
     400          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nRed);
     401          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nGreen);
     402          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nBlue);
     403          18 :             ++iRow;
     404             :         }
     405             : 
     406           2 :         poSurfaceCurrentSpeedBand->m_poRAT = std::move(poRAT);
     407             : 
     408           2 :         poDS->SetBand(1, poSurfaceCurrentSpeedBand.release());
     409             : 
     410             :         // Create surfaceCurrentDirection band
     411             :         auto poSurfaceCurrentDirection =
     412           6 :             poValuesArray->GetView("[\"surfaceCurrentDirection\"]");
     413             :         auto poSurfaceCurrentDirectionDS = std::unique_ptr<GDALDataset>(
     414           4 :             poSurfaceCurrentDirection->AsClassicDataset(1, 0));
     415             :         auto poSurfaceCurrentDirectionBand = std::make_unique<S111RasterBand>(
     416           4 :             std::move(poSurfaceCurrentDirectionDS));
     417           2 :         poSurfaceCurrentDirectionBand->SetDescription(
     418           2 :             "surfaceCurrentDirection");
     419           2 :         poSurfaceCurrentDirectionBand->m_osUnitType = "degree";
     420           2 :         poSurfaceCurrentDirectionBand->GDALRasterBand::SetMetadataItem(
     421             :             "ANGLE_CONVENTION", "From true north, clockwise");
     422           2 :         poDS->SetBand(2, poSurfaceCurrentDirectionBand.release());
     423             :     }
     424             : 
     425           3 :     poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
     426             : 
     427             :     // Setup/check for pam .aux.xml.
     428           3 :     poDS->SetDescription(osFilename.c_str());
     429           3 :     poDS->TryLoadXML();
     430             : 
     431             :     // Setup overviews.
     432           3 :     poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
     433             : 
     434           3 :     return poDS.release();
     435             : }
     436             : 
     437             : /************************************************************************/
     438             : /*                      S111DatasetDriverUnload()                       */
     439             : /************************************************************************/
     440             : 
     441           6 : static void S111DatasetDriverUnload(GDALDriver *)
     442             : {
     443           6 :     HDF5UnloadFileDriver();
     444           6 : }
     445             : 
     446             : /************************************************************************/
     447             : /*                         GDALRegister_S111()                          */
     448             : /************************************************************************/
     449          11 : void GDALRegister_S111()
     450             : 
     451             : {
     452          11 :     if (!GDAL_CHECK_VERSION("S111"))
     453           0 :         return;
     454             : 
     455          11 :     if (GDALGetDriverByName(S111_DRIVER_NAME) != nullptr)
     456           0 :         return;
     457             : 
     458          11 :     GDALDriver *poDriver = new GDALDriver();
     459             : 
     460          11 :     S111DriverSetCommonMetadata(poDriver);
     461          11 :     poDriver->pfnOpen = S111Dataset::Open;
     462          11 :     poDriver->pfnUnloadDriver = S111DatasetDriverUnload;
     463             : 
     464          11 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     465             : }