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