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 :         CPLError(CE_Failure, CPLE_NotSupported,
     107             :                  "The S111 driver does not support update access.");
     108           0 :         return nullptr;
     109             :     }
     110             : 
     111           8 :     std::string osFilename(poOpenInfo->pszFilename);
     112           8 :     std::string osGroup;
     113           4 :     if (STARTS_WITH(poOpenInfo->pszFilename, "S111:"))
     114             :     {
     115             :         const CPLStringList aosTokens(
     116           3 :             CSLTokenizeString2(poOpenInfo->pszFilename, ":",
     117           3 :                                CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
     118             : 
     119           3 :         if (aosTokens.size() == 2)
     120             :         {
     121           0 :             osFilename = aosTokens[1];
     122             :         }
     123           3 :         else if (aosTokens.size() == 3)
     124             :         {
     125           3 :             osFilename = aosTokens[1];
     126           3 :             osGroup = aosTokens[2];
     127             :         }
     128             :         else
     129             :         {
     130           0 :             return nullptr;
     131             :         }
     132             :     }
     133             : 
     134           8 :     auto poDS = std::make_unique<S111Dataset>(osFilename);
     135           4 :     if (!poDS->Init())
     136           0 :         return nullptr;
     137             : 
     138           4 :     const auto &poRootGroup = poDS->m_poRootGroup;
     139             : 
     140          12 :     auto poSurfaceCurrent = poRootGroup->OpenGroup("SurfaceCurrent");
     141           4 :     if (!poSurfaceCurrent)
     142             :     {
     143           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     144             :                  "Cannot find /SurfaceCurrent group");
     145           0 :         return nullptr;
     146             :     }
     147             : 
     148             :     auto poDataCodingFormat =
     149          12 :         poSurfaceCurrent->GetAttribute("dataCodingFormat");
     150           8 :     if (!poDataCodingFormat ||
     151           4 :         poDataCodingFormat->GetDataType().GetClass() != GEDTC_NUMERIC)
     152             :     {
     153           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     154             :                  "Cannot find /SurfaceCurrent/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 S111 driver",
     162             :                  nDataCodingFormat);
     163           0 :         return nullptr;
     164             :     }
     165             : 
     166             :     // Read additional metadata
     167          12 :     for (const char *pszAttrName :
     168             :          {"methodCurrentsProduct", "minDatasetCurrentSpeed",
     169          16 :           "maxDatasetCurrentSpeed"})
     170             :     {
     171          36 :         auto poAttr = poSurfaceCurrent->GetAttribute(pszAttrName);
     172          12 :         if (poAttr)
     173             :         {
     174           8 :             const char *pszVal = poAttr->ReadAsString();
     175           8 :             if (pszVal)
     176             :             {
     177           8 :                 poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
     178             :             }
     179             :         }
     180             :     }
     181             : 
     182          12 :     auto poSurfaceCurrent01 = poSurfaceCurrent->OpenGroup("SurfaceCurrent.01");
     183           4 :     if (!poSurfaceCurrent01)
     184             :     {
     185           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     186             :                  "Cannot find /SurfaceCurrent/SurfaceCurrent.01 group");
     187           0 :         return nullptr;
     188             :     }
     189             : 
     190             :     // Read additional metadata
     191          16 :     for (const char *pszAttrName :
     192             :          {"timeRecordInterval", "dateTimeOfFirstRecord", "dateTimeOfLastRecord",
     193          20 :           "numberOfTimes"})
     194             :     {
     195          48 :         auto poAttr = poSurfaceCurrent01->GetAttribute(pszAttrName);
     196          16 :         if (poAttr)
     197             :         {
     198          16 :             const char *pszVal = poAttr->ReadAsString();
     199          16 :             if (pszVal)
     200             :             {
     201          16 :                 poDS->GDALDataset::SetMetadataItem(pszAttrName, pszVal);
     202             :             }
     203             :         }
     204             :     }
     205             : 
     206           4 :     if (auto poStartSequence =
     207           8 :             poSurfaceCurrent01->GetAttribute("startSequence"))
     208             :     {
     209           4 :         const char *pszStartSequence = poStartSequence->ReadAsString();
     210           4 :         if (pszStartSequence && !EQUAL(pszStartSequence, "0,0"))
     211             :         {
     212           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     213             :                      "startSequence (=%s) != 0,0 is not supported",
     214             :                      pszStartSequence);
     215           0 :             return nullptr;
     216             :         }
     217             :     }
     218             : 
     219           4 :     if (!S100GetNumPointsLongitudinalLatitudinal(
     220           4 :             poSurfaceCurrent01.get(), poDS->nRasterXSize, poDS->nRasterYSize))
     221             :     {
     222           0 :         return nullptr;
     223             :     }
     224             : 
     225           4 :     const bool bNorthUp = CPLTestBool(
     226           4 :         CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES"));
     227             : 
     228             :     // Compute geotransform
     229           4 :     poDS->m_bHasGT = S100GetGeoTransform(poSurfaceCurrent01.get(),
     230           4 :                                          poDS->m_adfGeoTransform, bNorthUp);
     231             : 
     232           4 :     if (osGroup.empty())
     233             :     {
     234           2 :         const auto aosGroupNames = poSurfaceCurrent01->GetGroupNames();
     235           1 :         int iSubDS = 1;
     236           2 :         for (const auto &osSubGroup : aosGroupNames)
     237             :         {
     238           2 :             if (auto poSubGroup = poSurfaceCurrent01->OpenGroup(osSubGroup))
     239             :             {
     240           1 :                 poDS->GDALDataset::SetMetadataItem(
     241             :                     CPLSPrintf("SUBDATASET_%d_NAME", iSubDS),
     242             :                     CPLSPrintf("S111:\"%s\":%s", osFilename.c_str(),
     243             :                                osSubGroup.c_str()),
     244             :                     "SUBDATASETS");
     245           2 :                 std::string osSubDSDesc = "Values for group ";
     246           1 :                 osSubDSDesc += osSubGroup;
     247           2 :                 const auto poTimePoint = poSubGroup->GetAttribute("timePoint");
     248           1 :                 if (poTimePoint)
     249             :                 {
     250           1 :                     const char *pszVal = poTimePoint->ReadAsString();
     251           1 :                     if (pszVal)
     252             :                     {
     253           1 :                         osSubDSDesc = "Values at timestamp ";
     254           1 :                         osSubDSDesc += pszVal;
     255             :                     }
     256             :                 }
     257           1 :                 poDS->GDALDataset::SetMetadataItem(
     258             :                     CPLSPrintf("SUBDATASET_%d_DESC", iSubDS),
     259             :                     osSubDSDesc.c_str(), "SUBDATASETS");
     260           1 :                 ++iSubDS;
     261             :             }
     262             :         }
     263             :     }
     264             :     else
     265             :     {
     266           3 :         auto poGroup = poSurfaceCurrent01->OpenGroup(osGroup);
     267           3 :         if (!poGroup)
     268             :         {
     269           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     270             :                      "Cannot find /SurfaceCurrent/SurfaceCurrent.01/%s group",
     271             :                      osGroup.c_str());
     272           1 :             return nullptr;
     273             :         }
     274             : 
     275           4 :         auto poValuesArray = poGroup->OpenMDArray("values");
     276           2 :         if (!poValuesArray)
     277             :         {
     278           0 :             CPLError(
     279             :                 CE_Failure, CPLE_AppDefined,
     280             :                 "Cannot find /SurfaceCurrent/SurfaceCurrent.01/%s/values array",
     281             :                 osGroup.c_str());
     282           0 :             return nullptr;
     283             :         }
     284             : 
     285           2 :         if (poValuesArray->GetDimensionCount() != 2)
     286             :         {
     287           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     288             :                      "Wrong dimension count for %s",
     289           0 :                      poValuesArray->GetFullName().c_str());
     290           0 :             return nullptr;
     291             :         }
     292             : 
     293           2 :         const auto &oType = poValuesArray->GetDataType();
     294           2 :         if (oType.GetClass() != GEDTC_COMPOUND)
     295             :         {
     296           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Wrong data type for %s",
     297           0 :                      poValuesArray->GetFullName().c_str());
     298           0 :             return nullptr;
     299             :         }
     300             : 
     301           2 :         const auto &oComponents = oType.GetComponents();
     302           4 :         if (!(oComponents.size() == 2 &&
     303           4 :               ((oComponents[0]->GetName() == "surfaceCurrentSpeed" &&
     304           2 :                 oComponents[0]->GetType().GetNumericDataType() == GDT_Float32 &&
     305           4 :                 oComponents[1]->GetName() == "surfaceCurrentDirection" &&
     306           2 :                 oComponents[1]->GetType().GetNumericDataType() ==
     307           0 :                     GDT_Float32) ||
     308             :                // S111US_20170829.0100_W078.N44_F2_loofs_type2.h5 has direction first...
     309           0 :                (oComponents[0]->GetName() == "surfaceCurrentDirection" &&
     310           0 :                 oComponents[0]->GetType().GetNumericDataType() == GDT_Float32 &&
     311           0 :                 oComponents[1]->GetName() == "surfaceCurrentSpeed" &&
     312           0 :                 oComponents[1]->GetType().GetNumericDataType() ==
     313             :                     GDT_Float32))))
     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 surfaceCurrentSpeed band
     344             :         auto poSurfaceCurrentSpeed =
     345           6 :             poValuesArray->GetView("[\"surfaceCurrentSpeed\"]");
     346             :         auto poSurfaceCurrentSpeedDS = std::unique_ptr<GDALDataset>(
     347           4 :             poSurfaceCurrentSpeed->AsClassicDataset(1, 0));
     348             :         auto poSurfaceCurrentSpeedBand = std::make_unique<S111RasterBand>(
     349           4 :             std::move(poSurfaceCurrentSpeedDS));
     350           2 :         poSurfaceCurrentSpeedBand->SetDescription("surfaceCurrentSpeed");
     351           2 :         poSurfaceCurrentSpeedBand->m_osUnitType = "knots";
     352             : 
     353             :         // From S-111 v1.2 table 9.1 (Speed ranges) and 9.2 (Colour schemas)
     354           4 :         auto poRAT = std::make_unique<GDALDefaultRasterAttributeTable>();
     355           2 :         poRAT->CreateColumn("speed_band", GFT_Integer, GFU_Generic);
     356           2 :         poRAT->CreateColumn("min_speed", GFT_Real, GFU_Min);
     357           2 :         poRAT->CreateColumn("width_band", GFT_Real, GFU_Generic);
     358           2 :         poRAT->CreateColumn("color", GFT_String, GFU_Generic);
     359           2 :         poRAT->CreateColumn("red", GFT_Integer, GFU_RedMin);
     360           2 :         poRAT->CreateColumn("green", GFT_Integer, GFU_GreenMin);
     361           2 :         poRAT->CreateColumn("blue", GFT_Integer, GFU_BlueMin);
     362             : 
     363             :         const struct
     364             :         {
     365             :             int nSpeedBand;
     366             :             double dfMinSpeed;
     367             :             double dfWidthBand;
     368             :             const char *pszColor;
     369             :             int nRed;
     370             :             int nGreen;
     371             :             int nBlue;
     372           2 :         } aoRatValues[] = {
     373             :             {1, 0.0, 0.5, "purple", 118, 82, 226},
     374             :             {2, 0.5, 0.5, "dark blue", 72, 152, 211},
     375             :             {3, 1.0, 1.0, "light blue", 97, 203, 229},
     376             :             {4, 2.0, 1.0, "dark green", 109, 188, 69},
     377             :             {5, 3.0, 2.0, "light green", 180, 220, 0},
     378             :             {6, 5.0, 2.0, "yellow green", 205, 193, 0},
     379             :             {7, 7.0, 3.0, "orange", 248, 167, 24},
     380             :             {8, 10.0, 3.0, "pink", 247, 162, 157},
     381             :             {9, 13.0, 86.0, "red", 255, 30, 30},
     382             :         };
     383             : 
     384           2 :         int iRow = 0;
     385          20 :         for (const auto &oRecord : aoRatValues)
     386             :         {
     387          18 :             int iCol = 0;
     388          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nSpeedBand);
     389          18 :             poRAT->SetValue(iRow, iCol++, oRecord.dfMinSpeed);
     390          18 :             poRAT->SetValue(iRow, iCol++, oRecord.dfWidthBand);
     391          18 :             poRAT->SetValue(iRow, iCol++, oRecord.pszColor);
     392          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nRed);
     393          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nGreen);
     394          18 :             poRAT->SetValue(iRow, iCol++, oRecord.nBlue);
     395          18 :             ++iRow;
     396             :         }
     397             : 
     398           2 :         poSurfaceCurrentSpeedBand->m_poRAT = std::move(poRAT);
     399             : 
     400           2 :         poDS->SetBand(1, poSurfaceCurrentSpeedBand.release());
     401             : 
     402             :         // Create surfaceCurrentDirection band
     403             :         auto poSurfaceCurrentDirection =
     404           6 :             poValuesArray->GetView("[\"surfaceCurrentDirection\"]");
     405             :         auto poSurfaceCurrentDirectionDS = std::unique_ptr<GDALDataset>(
     406           4 :             poSurfaceCurrentDirection->AsClassicDataset(1, 0));
     407             :         auto poSurfaceCurrentDirectionBand = std::make_unique<S111RasterBand>(
     408           4 :             std::move(poSurfaceCurrentDirectionDS));
     409           2 :         poSurfaceCurrentDirectionBand->SetDescription(
     410           2 :             "surfaceCurrentDirection");
     411           2 :         poSurfaceCurrentDirectionBand->m_osUnitType = "degree";
     412           2 :         poSurfaceCurrentDirectionBand->GDALRasterBand::SetMetadataItem(
     413             :             "ANGLE_CONVENTION", "From true north, clockwise");
     414           2 :         poDS->SetBand(2, poSurfaceCurrentDirectionBand.release());
     415             :     }
     416             : 
     417           3 :     poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
     418             : 
     419             :     // Setup/check for pam .aux.xml.
     420           3 :     poDS->SetDescription(osFilename.c_str());
     421           3 :     poDS->TryLoadXML();
     422             : 
     423             :     // Setup overviews.
     424           3 :     poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str());
     425             : 
     426           3 :     return poDS.release();
     427             : }
     428             : 
     429             : /************************************************************************/
     430             : /*                      S111DatasetDriverUnload()                       */
     431             : /************************************************************************/
     432             : 
     433           6 : static void S111DatasetDriverUnload(GDALDriver *)
     434             : {
     435           6 :     HDF5UnloadFileDriver();
     436           6 : }
     437             : 
     438             : /************************************************************************/
     439             : /*                         GDALRegister_S111()                          */
     440             : /************************************************************************/
     441          10 : void GDALRegister_S111()
     442             : 
     443             : {
     444          10 :     if (!GDAL_CHECK_VERSION("S111"))
     445           0 :         return;
     446             : 
     447          10 :     if (GDALGetDriverByName(S111_DRIVER_NAME) != nullptr)
     448           0 :         return;
     449             : 
     450          10 :     GDALDriver *poDriver = new GDALDriver();
     451             : 
     452          10 :     S111DriverSetCommonMetadata(poDriver);
     453          10 :     poDriver->pfnOpen = S111Dataset::Open;
     454          10 :     poDriver->pfnUnloadDriver = S111DatasetDriverUnload;
     455             : 
     456          10 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     457             : }