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