Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Hierarchical Data Format Release 5 (HDF5)
       4             :  * Purpose:  Read S100 bathymetric 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 "s100.h"
      14             : #include "hdf5dataset.h"
      15             : 
      16             : #include <algorithm>
      17             : #include <cmath>
      18             : 
      19             : /************************************************************************/
      20             : /*                       S100BaseDataset()                              */
      21             : /************************************************************************/
      22             : 
      23          29 : S100BaseDataset::S100BaseDataset(const std::string &osFilename)
      24          29 :     : m_osFilename(osFilename)
      25             : {
      26          29 : }
      27             : 
      28             : /************************************************************************/
      29             : /*                              Init()                                  */
      30             : /************************************************************************/
      31             : 
      32          29 : bool S100BaseDataset::Init()
      33             : {
      34             :     // Open the file as an HDF5 file.
      35          29 :     hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
      36          29 :     H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
      37          29 :     hid_t hHDF5 = H5Fopen(m_osFilename.c_str(), H5F_ACC_RDONLY, fapl);
      38          29 :     H5Pclose(fapl);
      39          29 :     if (hHDF5 < 0)
      40           0 :         return false;
      41             : 
      42          58 :     auto poSharedResources = GDAL::HDF5SharedResources::Create(m_osFilename);
      43          29 :     poSharedResources->m_hHDF5 = hHDF5;
      44             : 
      45          29 :     m_poRootGroup = HDF5Dataset::OpenGroup(poSharedResources);
      46          29 :     if (m_poRootGroup == nullptr)
      47           0 :         return false;
      48             : 
      49          29 :     S100ReadSRS(m_poRootGroup.get(), m_oSRS);
      50             : 
      51          29 :     S100ReadVerticalDatum(this, m_poRootGroup.get());
      52             : 
      53             :     m_osMetadataFile =
      54          29 :         S100ReadMetadata(this, m_osFilename, m_poRootGroup.get());
      55             : 
      56          29 :     return true;
      57             : }
      58             : 
      59             : /************************************************************************/
      60             : /*                          GetGeoTransform()                           */
      61             : /************************************************************************/
      62             : 
      63          12 : CPLErr S100BaseDataset::GetGeoTransform(double *padfGeoTransform)
      64             : 
      65             : {
      66          12 :     if (m_bHasGT)
      67             :     {
      68          12 :         memcpy(padfGeoTransform, m_adfGeoTransform, sizeof(double) * 6);
      69          12 :         return CE_None;
      70             :     }
      71             : 
      72           0 :     return GDALPamDataset::GetGeoTransform(padfGeoTransform);
      73             : }
      74             : 
      75             : /************************************************************************/
      76             : /*                         GetSpatialRef()                              */
      77             : /************************************************************************/
      78             : 
      79          10 : const OGRSpatialReference *S100BaseDataset::GetSpatialRef() const
      80             : {
      81          10 :     if (!m_oSRS.IsEmpty())
      82          10 :         return &m_oSRS;
      83           0 :     return GDALPamDataset::GetSpatialRef();
      84             : }
      85             : 
      86             : /************************************************************************/
      87             : /*                         GetFileList()                                */
      88             : /************************************************************************/
      89             : 
      90           4 : char **S100BaseDataset::GetFileList()
      91             : {
      92           4 :     char **papszFileList = GDALPamDataset::GetFileList();
      93           4 :     if (!m_osMetadataFile.empty())
      94           4 :         papszFileList = CSLAddString(papszFileList, m_osMetadataFile.c_str());
      95           4 :     return papszFileList;
      96             : }
      97             : 
      98             : /************************************************************************/
      99             : /*                            S100ReadSRS()                             */
     100             : /************************************************************************/
     101             : 
     102          55 : bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS)
     103             : {
     104             :     // Get SRS
     105         110 :     auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS");
     106          87 :     if (poHorizontalCRS &&
     107          87 :         poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC)
     108             :     {
     109             :         // horizontalCRS is v2.2
     110          32 :         oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     111          32 :         if (oSRS.importFromEPSG(poHorizontalCRS->ReadAsInt()) != OGRERR_NONE)
     112             :         {
     113           0 :             oSRS.Clear();
     114             :         }
     115             :     }
     116             :     else
     117             :     {
     118             :         auto poHorizontalDatumReference =
     119          69 :             poRootGroup->GetAttribute("horizontalDatumReference");
     120             :         auto poHorizontalDatumValue =
     121          69 :             poRootGroup->GetAttribute("horizontalDatumValue");
     122          23 :         if (poHorizontalDatumReference && poHorizontalDatumValue)
     123             :         {
     124             :             const char *pszAuthName =
     125          23 :                 poHorizontalDatumReference->ReadAsString();
     126          23 :             const char *pszAuthCode = poHorizontalDatumValue->ReadAsString();
     127          23 :             if (pszAuthName && pszAuthCode)
     128             :             {
     129          23 :                 oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     130          46 :                 if (oSRS.SetFromUserInput(
     131          46 :                         (std::string(pszAuthName) + ':' + pszAuthCode).c_str(),
     132             :                         OGRSpatialReference::
     133          23 :                             SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
     134             :                     OGRERR_NONE)
     135             :                 {
     136           0 :                     oSRS.Clear();
     137             :                 }
     138             :             }
     139             :         }
     140             :     }
     141         110 :     return !oSRS.IsEmpty();
     142             : }
     143             : 
     144             : /************************************************************************/
     145             : /*               S100GetNumPointsLongitudinalLatitudinal()              */
     146             : /************************************************************************/
     147             : 
     148          61 : bool S100GetNumPointsLongitudinalLatitudinal(const GDALGroup *poGroup,
     149             :                                              int &nNumPointsLongitudinal,
     150             :                                              int &nNumPointsLatitudinal)
     151             : {
     152         183 :     auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
     153         183 :     auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
     154             :     auto poNumPointsLongitudinal =
     155         183 :         poGroup->GetAttribute("numPointsLongitudinal");
     156         183 :     auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal");
     157          61 :     if (poSpacingX &&
     158         183 :         poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
     159         122 :         poSpacingY &&
     160         122 :         poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 &&
     161          61 :         poNumPointsLongitudinal &&
     162          61 :         GDALDataTypeIsInteger(
     163         122 :             poNumPointsLongitudinal->GetDataType().GetNumericDataType()) &&
     164         183 :         poNumPointsLatitudinal &&
     165          61 :         GDALDataTypeIsInteger(
     166          61 :             poNumPointsLatitudinal->GetDataType().GetNumericDataType()))
     167             :     {
     168          61 :         nNumPointsLongitudinal = poNumPointsLongitudinal->ReadAsInt();
     169          61 :         nNumPointsLatitudinal = poNumPointsLatitudinal->ReadAsInt();
     170             : 
     171             :         // Those are optional, but use them when available, to detect
     172             :         // potential inconsistency
     173         183 :         auto poEastBoundLongitude = poGroup->GetAttribute("eastBoundLongitude");
     174         183 :         auto poWestBoundLongitude = poGroup->GetAttribute("westBoundLongitude");
     175             :         auto poSouthBoundLongitude =
     176         183 :             poGroup->GetAttribute("southBoundLatitude");
     177         122 :         auto poNorthBoundLatitude = poGroup->GetAttribute("northBoundLatitude");
     178          61 :         if (poEastBoundLongitude &&
     179           0 :             GDALDataTypeIsFloating(
     180          61 :                 poEastBoundLongitude->GetDataType().GetNumericDataType()) &&
     181           0 :             poWestBoundLongitude &&
     182           0 :             GDALDataTypeIsFloating(
     183           0 :                 poWestBoundLongitude->GetDataType().GetNumericDataType()) &&
     184           0 :             poSouthBoundLongitude &&
     185           0 :             GDALDataTypeIsFloating(
     186           0 :                 poSouthBoundLongitude->GetDataType().GetNumericDataType()) &&
     187          61 :             poNorthBoundLatitude &&
     188           0 :             GDALDataTypeIsFloating(
     189           0 :                 poNorthBoundLatitude->GetDataType().GetNumericDataType()))
     190             :         {
     191           0 :             const double dfSpacingX = poSpacingX->ReadAsDouble();
     192           0 :             const double dfSpacingY = poSpacingY->ReadAsDouble();
     193             : 
     194           0 :             const double dfEast = poEastBoundLongitude->ReadAsDouble();
     195           0 :             const double dfWest = poWestBoundLongitude->ReadAsDouble();
     196           0 :             const double dfSouth = poSouthBoundLongitude->ReadAsDouble();
     197           0 :             const double dfNorth = poNorthBoundLatitude->ReadAsDouble();
     198           0 :             if (std::fabs((dfWest + nNumPointsLongitudinal * dfSpacingX) -
     199           0 :                           dfEast) < 5 * dfSpacingX &&
     200           0 :                 std::fabs((dfSouth + nNumPointsLatitudinal * dfSpacingY) -
     201           0 :                           dfNorth) < 5 * dfSpacingY)
     202             :             {
     203             :                 // We need up to 5 spacings for product
     204             :                 // S-111 Trial Data Set Release 1.1/111UK_20210401T000000Z_SolentAndAppr_dcf2.h5
     205             :             }
     206             :             else
     207             :             {
     208           0 :                 CPLError(
     209             :                     CE_Warning, CPLE_AppDefined,
     210             :                     "Caution: "
     211             :                     "eastBoundLongitude/westBoundLongitude/southBoundLatitude/"
     212             :                     "northBoundLatitude/gridSpacingLongitudinal/"
     213             :                     "gridSpacingLatitudinal/numPointsLongitudinal/"
     214             :                     "numPointsLatitudinal do not seem to be consistent");
     215           0 :                 CPLDebug("S100", "Computed east = %f. Actual = %f",
     216           0 :                          dfWest + nNumPointsLongitudinal * dfSpacingX, dfEast);
     217           0 :                 CPLDebug("S100", "Computed north = %f. Actual = %f",
     218           0 :                          dfSouth + nNumPointsLatitudinal * dfSpacingY, dfNorth);
     219             :             }
     220             :         }
     221             : 
     222          61 :         return true;
     223             :     }
     224           0 :     return false;
     225             : }
     226             : 
     227             : /************************************************************************/
     228             : /*                        S100GetGeoTransform()                         */
     229             : /************************************************************************/
     230             : 
     231          27 : bool S100GetGeoTransform(const GDALGroup *poGroup, double adfGeoTransform[6],
     232             :                          bool bNorthUp)
     233             : {
     234          81 :     auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
     235          81 :     auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
     236          81 :     auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
     237          81 :     auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
     238          27 :     if (poOriginX &&
     239          81 :         poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
     240          54 :         poOriginY &&
     241          54 :         poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
     242          54 :         poSpacingX &&
     243          54 :         poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
     244          81 :         poSpacingY &&
     245          27 :         poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64)
     246             :     {
     247          27 :         int nNumPointsLongitudinal = 0;
     248          27 :         int nNumPointsLatitudinal = 0;
     249          27 :         if (!S100GetNumPointsLongitudinalLatitudinal(
     250             :                 poGroup, nNumPointsLongitudinal, nNumPointsLatitudinal))
     251           0 :             return false;
     252             : 
     253          27 :         const double dfSpacingX = poSpacingX->ReadAsDouble();
     254          27 :         const double dfSpacingY = poSpacingY->ReadAsDouble();
     255             : 
     256          27 :         adfGeoTransform[0] = poOriginX->ReadAsDouble();
     257          27 :         adfGeoTransform[3] =
     258          27 :             poOriginY->ReadAsDouble() +
     259          27 :             (bNorthUp ? dfSpacingY * (nNumPointsLatitudinal - 1) : 0);
     260          27 :         adfGeoTransform[1] = dfSpacingX;
     261          27 :         adfGeoTransform[5] = bNorthUp ? -dfSpacingY : dfSpacingY;
     262             : 
     263             :         // From pixel-center convention to pixel-corner convention
     264          27 :         adfGeoTransform[0] -= adfGeoTransform[1] / 2;
     265          27 :         adfGeoTransform[3] -= adfGeoTransform[5] / 2;
     266             : 
     267          27 :         return true;
     268             :     }
     269           0 :     return false;
     270             : }
     271             : 
     272             : /************************************************************************/
     273             : /*                        S100GetDimensions()                           */
     274             : /************************************************************************/
     275             : 
     276          26 : bool S100GetDimensions(
     277             :     const GDALGroup *poGroup,
     278             :     std::vector<std::shared_ptr<GDALDimension>> &apoDims,
     279             :     std::vector<std::shared_ptr<GDALMDArray>> &apoIndexingVars)
     280             : {
     281          78 :     auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
     282          78 :     auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
     283          78 :     auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
     284          78 :     auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
     285          26 :     if (poOriginX &&
     286          78 :         poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
     287          52 :         poOriginY &&
     288          52 :         poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
     289          52 :         poSpacingX &&
     290          52 :         poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
     291          78 :         poSpacingY &&
     292          26 :         poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64)
     293             :     {
     294          26 :         int nNumPointsLongitudinal = 0;
     295          26 :         int nNumPointsLatitudinal = 0;
     296          26 :         if (!S100GetNumPointsLongitudinalLatitudinal(
     297             :                 poGroup, nNumPointsLongitudinal, nNumPointsLatitudinal))
     298           0 :             return false;
     299             : 
     300             :         {
     301             :             auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
     302          52 :                 std::string(), "Y", GDAL_DIM_TYPE_HORIZONTAL_Y, std::string(),
     303          52 :                 nNumPointsLatitudinal);
     304             :             auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
     305          78 :                 std::string(), poDim->GetName(), poDim,
     306         104 :                 poOriginY->ReadAsDouble(), poSpacingY->ReadAsDouble(), 0);
     307          26 :             poDim->SetIndexingVariable(poIndexingVar);
     308          26 :             apoDims.emplace_back(poDim);
     309          26 :             apoIndexingVars.emplace_back(poIndexingVar);
     310             :         }
     311             : 
     312             :         {
     313             :             auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
     314          52 :                 std::string(), "X", GDAL_DIM_TYPE_HORIZONTAL_X, std::string(),
     315          52 :                 nNumPointsLongitudinal);
     316             :             auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
     317          78 :                 std::string(), poDim->GetName(), poDim,
     318         104 :                 poOriginX->ReadAsDouble(), poSpacingX->ReadAsDouble(), 0);
     319          26 :             poDim->SetIndexingVariable(poIndexingVar);
     320          26 :             apoDims.emplace_back(poDim);
     321          26 :             apoIndexingVars.emplace_back(poIndexingVar);
     322             :         }
     323             : 
     324          26 :         return true;
     325             :     }
     326           0 :     return false;
     327             : }
     328             : 
     329             : /************************************************************************/
     330             : /*                       S100ReadVerticalDatum()                        */
     331             : /************************************************************************/
     332             : 
     333          29 : void S100ReadVerticalDatum(GDALDataset *poDS, const GDALGroup *poRootGroup)
     334             : {
     335             :     // https://iho.int/uploads/user/pubs/standards/s-100/S-100_5.0.0_Final_Clean_Web.pdf
     336             :     // Table S100_VerticalAndSoundingDatum page 20
     337             :     static const struct
     338             :     {
     339             :         int nCode;
     340             :         const char *pszMeaning;
     341             :         const char *pszAbbrev;
     342             :     } asVerticalDatums[] = {
     343             :         {1, "meanLowWaterSprings", "MLWS"},
     344             :         {2, "meanLowerLowWaterSprings", nullptr},
     345             :         {3, "meanSeaLevel", "MSL"},
     346             :         {4, "lowestLowWater", nullptr},
     347             :         {5, "meanLowWater", "MLW"},
     348             :         {6, "lowestLowWaterSprings", nullptr},
     349             :         {7, "approximateMeanLowWaterSprings", nullptr},
     350             :         {8, "indianSpringLowWater", nullptr},
     351             :         {9, "lowWaterSprings", nullptr},
     352             :         {10, "approximateLowestAstronomicalTide", nullptr},
     353             :         {11, "nearlyLowestLowWater", nullptr},
     354             :         {12, "meanLowerLowWater", "MLLW"},
     355             :         {13, "lowWater", "LW"},
     356             :         {14, "approximateMeanLowWater", nullptr},
     357             :         {15, "approximateMeanLowerLowWater", nullptr},
     358             :         {16, "meanHighWater", "MHW"},
     359             :         {17, "meanHighWaterSprings", "MHWS"},
     360             :         {18, "highWater", "HW"},
     361             :         {19, "approximateMeanSeaLevel", nullptr},
     362             :         {20, "highWaterSprings", nullptr},
     363             :         {21, "meanHigherHighWater", "MHHW"},
     364             :         {22, "equinoctialSpringLowWater", nullptr},
     365             :         {23, "lowestAstronomicalTide", "LAT"},
     366             :         {24, "localDatum", nullptr},
     367             :         {25, "internationalGreatLakesDatum1985", nullptr},
     368             :         {26, "meanWaterLevel", nullptr},
     369             :         {27, "lowerLowWaterLargeTide", nullptr},
     370             :         {28, "higherHighWaterLargeTide", nullptr},
     371             :         {29, "nearlyHighestHighWater", nullptr},
     372             :         {30, "highestAstronomicalTide", "HAT"},
     373             :         {44, "balticSeaChartDatum2000", nullptr},
     374             :         {46, "internationalGreatLakesDatum2020", nullptr},
     375             :     };
     376             : 
     377          87 :     auto poVerticalDatum = poRootGroup->GetAttribute("verticalDatum");
     378          58 :     if (poVerticalDatum &&
     379          58 :         poVerticalDatum->GetDataType().GetClass() == GEDTC_NUMERIC)
     380             :     {
     381          29 :         bool bFound = false;
     382          29 :         const auto nVal = poVerticalDatum->ReadAsInt();
     383         348 :         for (const auto &sVerticalDatum : asVerticalDatums)
     384             :         {
     385         348 :             if (sVerticalDatum.nCode == nVal)
     386             :             {
     387          29 :                 bFound = true;
     388          29 :                 poDS->GDALDataset::SetMetadataItem("VERTICAL_DATUM_MEANING",
     389          29 :                                                    sVerticalDatum.pszMeaning);
     390          29 :                 if (sVerticalDatum.pszAbbrev)
     391          29 :                     poDS->GDALDataset::SetMetadataItem(
     392          29 :                         "VERTICAL_DATUM_ABBREV", sVerticalDatum.pszAbbrev);
     393          29 :                 break;
     394             :             }
     395             :         }
     396          29 :         if (!bFound)
     397             :         {
     398           0 :             poDS->GDALDataset::SetMetadataItem("verticalDatum",
     399             :                                                CPLSPrintf("%d", nVal));
     400             :         }
     401             :     }
     402          29 : }
     403             : 
     404             : /************************************************************************/
     405             : /*                         S100ReadMetadata()                           */
     406             : /************************************************************************/
     407             : 
     408          29 : std::string S100ReadMetadata(GDALDataset *poDS, const std::string &osFilename,
     409             :                              const GDALGroup *poRootGroup)
     410             : {
     411          29 :     std::string osMetadataFile;
     412         297 :     for (const auto &poAttr : poRootGroup->GetAttributes())
     413             :     {
     414         268 :         const auto &osName = poAttr->GetName();
     415         268 :         if (osName == "metadata")
     416             :         {
     417          29 :             const char *pszVal = poAttr->ReadAsString();
     418          29 :             if (pszVal && pszVal[0])
     419             :             {
     420             :                 osMetadataFile = CPLFormFilename(CPLGetPath(osFilename.c_str()),
     421          29 :                                                  pszVal, nullptr);
     422             :                 VSIStatBufL sStat;
     423          29 :                 if (VSIStatL(osMetadataFile.c_str(), &sStat) != 0)
     424             :                 {
     425             :                     // Test products from https://data.admiralty.co.uk/portal/apps/sites/#/marine-data-portal/pages/s-100
     426             :                     // advertise a metadata filename starting with "MD_", per the spec,
     427             :                     // but the actual filename does not start with "MD_"...
     428           6 :                     if (STARTS_WITH(pszVal, "MD_"))
     429             :                     {
     430             :                         osMetadataFile =
     431             :                             CPLFormFilename(CPLGetPath(osFilename.c_str()),
     432           6 :                                             pszVal + strlen("MD_"), nullptr);
     433           6 :                         if (VSIStatL(osMetadataFile.c_str(), &sStat) != 0)
     434             :                         {
     435           6 :                             osMetadataFile.clear();
     436             :                         }
     437             :                     }
     438             :                     else
     439             :                     {
     440           0 :                         osMetadataFile.clear();
     441             :                     }
     442             :                 }
     443             :             }
     444             :         }
     445         460 :         else if (osName != "horizontalCRS" &&
     446         431 :                  osName != "horizontalDatumReference" &&
     447         409 :                  osName != "horizontalDatumValue" &&
     448         369 :                  osName != "productSpecification" &&
     449         340 :                  osName != "eastBoundLongitude" &&
     450         340 :                  osName != "northBoundLatitude" &&
     451         340 :                  osName != "southBoundLatitude" &&
     452         510 :                  osName != "westBoundLongitude" && osName != "extentTypeCode" &&
     453         456 :                  osName != "verticalCS" && osName != "verticalCoordinateBase" &&
     454         594 :                  osName != "verticalDatumReference" &&
     455         116 :                  osName != "verticalDatum")
     456             :         {
     457          87 :             const char *pszVal = poAttr->ReadAsString();
     458          87 :             if (pszVal && pszVal[0])
     459          87 :                 poDS->GDALDataset::SetMetadataItem(osName.c_str(), pszVal);
     460             :         }
     461             :     }
     462          29 :     return osMetadataFile;
     463             : }