From a4efedb838fd484d046a9bc3dd9b1fdcdec9c58f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 17 Oct 2023 15:58:14 +0200 Subject: [PATCH] S102: add multidimensional support --- autotest/gdrivers/data/s102/generate_test.py | 2 + autotest/gdrivers/data/s102/test_s102_v2.1.h5 | Bin 12792 -> 12984 bytes autotest/gdrivers/data/s102/test_s102_v2.2.h5 | Bin 12896 -> 13056 bytes autotest/gdrivers/s102.py | 27 +++ frmts/hdf5/CMakeLists.txt | 2 + frmts/hdf5/hdf5dataset.h | 5 +- frmts/hdf5/hdf5multidim.cpp | 47 +++++ frmts/hdf5/s100.cpp | 183 ++++++++++++++++++ frmts/hdf5/s100.h | 47 +++++ frmts/hdf5/s102dataset.cpp | 78 +------- 10 files changed, 320 insertions(+), 71 deletions(-) create mode 100644 frmts/hdf5/s100.cpp create mode 100644 frmts/hdf5/s100.h diff --git a/autotest/gdrivers/data/s102/generate_test.py b/autotest/gdrivers/data/s102/generate_test.py index a3a7ee321de5..861614e88290 100755 --- a/autotest/gdrivers/data/s102/generate_test.py +++ b/autotest/gdrivers/data/s102/generate_test.py @@ -63,6 +63,8 @@ def generate(filename, version): BathymetryCoverage_01.attrs["gridOriginLatitude"] = np.float64(48) BathymetryCoverage_01.attrs["gridSpacingLongitudinal"] = np.float64(0.4) BathymetryCoverage_01.attrs["gridSpacingLatitudinal"] = np.float64(0.5) + BathymetryCoverage_01.attrs["numPointsLongitudinal"] = np.uint32(values.shape[1]) + BathymetryCoverage_01.attrs["numPointsLatitudinal"] = np.uint32(values.shape[0]) f.create_group("Group_F") diff --git a/autotest/gdrivers/data/s102/test_s102_v2.1.h5 b/autotest/gdrivers/data/s102/test_s102_v2.1.h5 index 021894a01d788c8df0a62d30068e588065cfac59..ea3bfdc44d1ca2cadf57edf30d89fe81e7ba33a2 100644 GIT binary patch delta 334 zcmey7yd!mj2Gb6siCPx&ybO#C3=9%L3<3f`2m&kQ!4$)e$${MJ8&?#tO1+SQh=5ds zfPgH7HjtmV@f$NE1LtIIc3H+BAoZII*spLfS5DO1>@L8_Q_sWTz`z1j%*Y_dz{9}7 zkXM=;ke``XQtXqTm!4Ttnv$89m;*8qh(YQY7!<%H%qC_qj{&F`W{xPDIf*4m#-N)c z2=*Yz3?`^SlRacKxDf8&TqR@2#1tSud6!(A5G(*D7(x_lFmBe8pT-Du1INUTwu>hy X3QV3L$uU_$UTE?mg=Lc`7-j$f_4YXU delta 179 zcmdmy`XhOQ2GbA2iCPv4JPeEs3=9%L3<3f`2m%uFV2UAwX|fl)+QtOPG7u3)2*n^E3tTPxxVC1RiVQ^qz0V-x>5M$tB;9$rr%?-%U%quDO${E&baMoO7DE6Nl%DJ%qrru6 u|K=(gJ0_+8`N_NF@|X%3H*3jnV+7jCF>#~p<_U^Y0+Syp%$vNza0dX#{x!A$ delta 140 zcmZojdyq0ggDJsiqLv?*gggTnR4`5UVprR^qJUK@LIxtl2%#7np!@}L6E}XFY``rs z*@2ydF=Jw*`sM=mD;&%W^K>`63o!Cbz9FN*g;2biOV*Bw>44m1FZn#C3k;ig$!}wv Vyuxq+GXn$1=4M58fyrJ*+W@)+B7Xn? diff --git a/autotest/gdrivers/s102.py b/autotest/gdrivers/s102.py index fe8ca7fe7a05..0d7f46abe2e1 100755 --- a/autotest/gdrivers/s102.py +++ b/autotest/gdrivers/s102.py @@ -175,3 +175,30 @@ def test_s102_identify_fallback_through_HDF5_driver(): assert ds.GetDriver().GetDescription() == "S102" del ds assert not os.path.exists("data/s102/test_s102_v2.1.h5.aux.xml") + + +############################################################################### + + +def test_s102_multidim(): + ds = gdal.OpenEx("data/s102/test_s102_v2.1.h5", gdal.OF_MULTIDIM_RASTER) + rg = ds.GetRootGroup() + ar = rg.OpenMDArrayFromFullname( + "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values" + ) + assert ar.GetSpatialRef().GetAuthorityCode(None) == "4326" + + assert ar.GetDimensions()[0].GetName() == "Y" + y = ar.GetDimensions()[0].GetIndexingVariable() + y_data = struct.unpack("d" * y.GetDimensions()[0].GetSize(), y.Read()) + assert y_data[0] == 48.0 + assert y_data[-1] == 48.5 + + assert ar.GetDimensions()[1].GetName() == "X" + x = ar.GetDimensions()[1].GetIndexingVariable() + x_data = struct.unpack("d" * x.GetDimensions()[0].GetSize(), x.Read()) + assert x_data[0] == 2.0 + assert x_data[-1] == 2.8 + + # Check that it doesn't go into infinite recursion + gdal.MultiDimInfo(ds) diff --git a/frmts/hdf5/CMakeLists.txt b/frmts/hdf5/CMakeLists.txt index 100d13567202..81864f7e1116 100644 --- a/frmts/hdf5/CMakeLists.txt +++ b/frmts/hdf5/CMakeLists.txt @@ -10,6 +10,8 @@ set(SOURCE bagdataset.cpp hdf5multidim.cpp hdf5eosparser.cpp + s100.cpp + s100.h s102dataset.cpp ) diff --git a/frmts/hdf5/hdf5dataset.h b/frmts/hdf5/hdf5dataset.h index 6c255ed6577c..3e39520bb5ca 100644 --- a/frmts/hdf5/hdf5dataset.h +++ b/frmts/hdf5/hdf5dataset.h @@ -131,8 +131,7 @@ class HDF5SharedResources m_oMapEOSGridNameToDimensions{}; std::map>> m_oMapEOSSwathNameToDimensions{}; - std::map> - m_oRefKeeper; + std::map> m_oRefKeeper; explicit HDF5SharedResources(const std::string &osFilename); @@ -199,7 +198,7 @@ class HDF5SharedResources return oIter->second; } - void KeepRef(const std::shared_ptr &poArray) + void KeepRef(const std::shared_ptr &poArray) { m_oRefKeeper[poArray->GetFullName()] = poArray; } diff --git a/frmts/hdf5/hdf5multidim.cpp b/frmts/hdf5/hdf5multidim.cpp index c27aa74aa503..a66796a6b0e0 100644 --- a/frmts/hdf5/hdf5multidim.cpp +++ b/frmts/hdf5/hdf5multidim.cpp @@ -27,6 +27,7 @@ #include "hdf5dataset.h" #include "hdf5eosparser.h" +#include "s100.h" #include #include @@ -1233,6 +1234,52 @@ void HDF5Array::InstantiateDimensions(const std::string &osParentName, } return; } + + // Special case for S102 + if (nDims == 2 && + GetFullName() == + "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values") + { + auto poRootGroup = m_poShared->GetRootGroup(); + if (poRootGroup) + { + m_poSRS = std::make_shared(); + if (S100ReadSRS(poRootGroup.get(), *(m_poSRS.get()))) + { + if (m_poSRS->GetDataAxisToSRSAxisMapping() == + std::vector{2, 1}) + m_poSRS->SetDataAxisToSRSAxisMapping({1, 2}); + else + m_poSRS->SetDataAxisToSRSAxisMapping({2, 1}); + } + else + { + m_poSRS.reset(); + } + + auto poBathymetryCoverage01 = + poRootGroup->OpenGroupFromFullname( + "/BathymetryCoverage/BathymetryCoverage.01"); + if (poBathymetryCoverage01) + { + std::vector> apoIndexingVars; + if (S100GetDimensions(poBathymetryCoverage01.get(), m_dims, + apoIndexingVars) && + m_dims.size() == 2 && + m_dims[0]->GetSize() == anDimSizes[0] && + m_dims[1]->GetSize() == anDimSizes[1]) + { + for (const auto &poIndexingVar : apoIndexingVars) + m_poShared->KeepRef(poIndexingVar); + return; + } + else + { + m_dims.clear(); + } + } + } + } } std::map> oMapFullNameToDim; diff --git a/frmts/hdf5/s100.cpp b/frmts/hdf5/s100.cpp new file mode 100644 index 000000000000..31882160e8c4 --- /dev/null +++ b/frmts/hdf5/s100.cpp @@ -0,0 +1,183 @@ +/****************************************************************************** + * + * Project: Hierarchical Data Format Release 5 (HDF5) + * Purpose: Read S100 bathymetric datasets. + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2023, Even Rouault + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "s100.h" + +/************************************************************************/ +/* S100ReadSRS() */ +/************************************************************************/ + +bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS) +{ + // Get SRS + auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS"); + if (poHorizontalCRS && + poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC) + { + // horizontalCRS is v2.2 + oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + if (oSRS.importFromEPSG(poHorizontalCRS->ReadAsInt()) != OGRERR_NONE) + { + oSRS.Clear(); + } + } + else + { + auto poHorizontalDatumReference = + poRootGroup->GetAttribute("horizontalDatumReference"); + auto poHorizontalDatumValue = + poRootGroup->GetAttribute("horizontalDatumValue"); + if (poHorizontalDatumReference && poHorizontalDatumValue) + { + const char *pszAuthName = + poHorizontalDatumReference->ReadAsString(); + const char *pszAuthCode = poHorizontalDatumValue->ReadAsString(); + if (pszAuthName && pszAuthCode) + { + oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + if (oSRS.SetFromUserInput( + (std::string(pszAuthName) + ':' + pszAuthCode).c_str(), + OGRSpatialReference:: + SET_FROM_USER_INPUT_LIMITATIONS_get()) != + OGRERR_NONE) + { + oSRS.Clear(); + } + } + } + } + return !oSRS.IsEmpty(); +} + +/************************************************************************/ +/* S100GetGeoTransform() */ +/************************************************************************/ + +bool S100GetGeoTransform(const GDALGroup *poGroup, double adfGeoTransform[6], + bool bNorthUp) +{ + auto poOriginX = poGroup->GetAttribute("gridOriginLongitude"); + auto poOriginY = poGroup->GetAttribute("gridOriginLatitude"); + auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal"); + auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal"); + auto poNumPointsLongitudinal = + poGroup->GetAttribute("numPointsLongitudinal"); + auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal"); + if (poOriginX && + poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 && + poOriginY && + poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 && + poSpacingX && + poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 && + poSpacingY && + poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 && + poNumPointsLongitudinal && + GDALDataTypeIsInteger( + poNumPointsLongitudinal->GetDataType().GetNumericDataType()) && + poNumPointsLatitudinal && + GDALDataTypeIsInteger( + poNumPointsLatitudinal->GetDataType().GetNumericDataType())) + { + adfGeoTransform[0] = poOriginX->ReadAsDouble(); + adfGeoTransform[3] = + poOriginY->ReadAsDouble() + + (bNorthUp ? poSpacingY->ReadAsDouble() * + (poNumPointsLatitudinal->ReadAsInt() - 1) + : 0); + adfGeoTransform[1] = poSpacingX->ReadAsDouble(); + adfGeoTransform[5] = + bNorthUp ? -poSpacingY->ReadAsDouble() : poSpacingY->ReadAsDouble(); + + // From pixel-center convention to pixel-corner convention + adfGeoTransform[0] -= adfGeoTransform[1] / 2; + adfGeoTransform[3] -= adfGeoTransform[5] / 2; + + return true; + } + return false; +} + +/************************************************************************/ +/* S100GetDimensions() */ +/************************************************************************/ + +bool S100GetDimensions( + const GDALGroup *poGroup, + std::vector> &apoDims, + std::vector> &apoIndexingVars) +{ + auto poOriginX = poGroup->GetAttribute("gridOriginLongitude"); + auto poOriginY = poGroup->GetAttribute("gridOriginLatitude"); + auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal"); + auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal"); + auto poNumPointsLongitudinal = + poGroup->GetAttribute("numPointsLongitudinal"); + auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal"); + if (poOriginX && + poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 && + poOriginY && + poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 && + poSpacingX && + poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 && + poSpacingY && + poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 && + poNumPointsLongitudinal && + GDALDataTypeIsInteger( + poNumPointsLongitudinal->GetDataType().GetNumericDataType()) && + poNumPointsLatitudinal && + GDALDataTypeIsInteger( + poNumPointsLatitudinal->GetDataType().GetNumericDataType())) + { + { + auto poDim = std::make_shared( + std::string(), "Y", GDAL_DIM_TYPE_HORIZONTAL_Y, std::string(), + poNumPointsLatitudinal->ReadAsInt()); + auto poIndexingVar = GDALMDArrayRegularlySpaced::Create( + std::string(), poDim->GetName(), poDim, + poOriginY->ReadAsDouble(), poSpacingY->ReadAsDouble(), 0); + poDim->SetIndexingVariable(poIndexingVar); + apoDims.emplace_back(poDim); + apoIndexingVars.emplace_back(poIndexingVar); + } + + { + auto poDim = std::make_shared( + std::string(), "X", GDAL_DIM_TYPE_HORIZONTAL_X, std::string(), + poNumPointsLongitudinal->ReadAsInt()); + auto poIndexingVar = GDALMDArrayRegularlySpaced::Create( + std::string(), poDim->GetName(), poDim, + poOriginX->ReadAsDouble(), poSpacingX->ReadAsDouble(), 0); + poDim->SetIndexingVariable(poIndexingVar); + apoDims.emplace_back(poDim); + apoIndexingVars.emplace_back(poIndexingVar); + } + + return true; + } + return false; +} diff --git a/frmts/hdf5/s100.h b/frmts/hdf5/s100.h new file mode 100644 index 000000000000..722065ee31ed --- /dev/null +++ b/frmts/hdf5/s100.h @@ -0,0 +1,47 @@ +/****************************************************************************** + * + * Project: Hierarchical Data Format Release 5 (HDF5) + * Purpose: Read S100 bathymetric datasets. + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2023, Even Rouault + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#ifndef S100_H +#define S100_H + +#include "cpl_port.h" + +#include "gdal_priv.h" +#include "ogr_spatialref.h" + +bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS); + +bool S100GetDimensions( + const GDALGroup *poGroup, + std::vector> &apoDims, + std::vector> &apoIndexingVars); + +bool S100GetGeoTransform(const GDALGroup *poGroup, double adfGeoTransform[6], + bool bNorthUp); + +#endif // S100_H diff --git a/frmts/hdf5/s102dataset.cpp b/frmts/hdf5/s102dataset.cpp index aa7db5b5b669..2e1df112084c 100644 --- a/frmts/hdf5/s102dataset.cpp +++ b/frmts/hdf5/s102dataset.cpp @@ -29,6 +29,7 @@ #include "cpl_port.h" #include "hdf5dataset.h" #include "gh5_convenience.h" +#include "s100.h" #include "gdal_priv.h" #include "gdal_proxy.h" @@ -210,6 +211,11 @@ GDALDataset *S102Dataset::Open(GDALOpenInfo *poOpenInfo) HDF5_GLOBAL_LOCK(); + if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER) + { + return HDF5Dataset::OpenMultiDim(poOpenInfo); + } + // Confirm the requested access is supported. if (poOpenInfo->eAccess == GA_Update) { @@ -388,76 +394,11 @@ GDALDataset *S102Dataset::Open(GDALOpenInfo *poOpenInfo) poDS->SetBand(2, poUncertaintyBand); // Compute geotransform - auto poOriginX = - poBathymetryCoverage01->GetAttribute("gridOriginLongitude"); - auto poOriginY = poBathymetryCoverage01->GetAttribute("gridOriginLatitude"); - auto poSpacingX = - poBathymetryCoverage01->GetAttribute("gridSpacingLongitudinal"); - auto poSpacingY = - poBathymetryCoverage01->GetAttribute("gridSpacingLatitudinal"); - if (poOriginX && - poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 && - poOriginY && - poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 && - poSpacingX && - poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 && - poSpacingY && - poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64) - { - poDS->m_adfGeoTransform[0] = poOriginX->ReadAsDouble(); - poDS->m_adfGeoTransform[3] = - poOriginY->ReadAsDouble() + - (bNorthUp ? poSpacingY->ReadAsDouble() * (poDS->nRasterYSize - 1) - : 0); - poDS->m_adfGeoTransform[1] = poSpacingX->ReadAsDouble(); - poDS->m_adfGeoTransform[5] = - bNorthUp ? -poSpacingY->ReadAsDouble() : poSpacingY->ReadAsDouble(); - - // From pixel-center convention to pixel-corner convention - poDS->m_adfGeoTransform[0] -= poDS->m_adfGeoTransform[1] / 2; - poDS->m_adfGeoTransform[3] -= poDS->m_adfGeoTransform[5] / 2; - - poDS->m_bHasGT = true; - } + poDS->m_bHasGT = S100GetGeoTransform(poBathymetryCoverage01.get(), + poDS->m_adfGeoTransform, bNorthUp); // Get SRS - auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS"); - if (poHorizontalCRS && - poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC) - { - // horizontalCRS is v2.2 - poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); - if (poDS->m_oSRS.importFromEPSG(poHorizontalCRS->ReadAsInt()) != - OGRERR_NONE) - { - poDS->m_oSRS.Clear(); - } - } - else - { - auto poHorizontalDatumReference = - poRootGroup->GetAttribute("horizontalDatumReference"); - auto poHorizontalDatumValue = - poRootGroup->GetAttribute("horizontalDatumValue"); - if (poHorizontalDatumReference && poHorizontalDatumValue) - { - const char *pszAuthName = - poHorizontalDatumReference->ReadAsString(); - const char *pszAuthCode = poHorizontalDatumValue->ReadAsString(); - if (pszAuthName && pszAuthCode) - { - poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); - if (poDS->m_oSRS.SetFromUserInput( - (std::string(pszAuthName) + ':' + pszAuthCode).c_str(), - OGRSpatialReference:: - SET_FROM_USER_INPUT_LIMITATIONS_get()) != - OGRERR_NONE) - { - poDS->m_oSRS.Clear(); - } - } - } - } + S100ReadSRS(poRootGroup.get(), poDS->m_oSRS); // https://iho.int/uploads/user/pubs/standards/s-100/S-100_5.0.0_Final_Clean_Web.pdf // Table S100_VerticalAndSoundingDatum page 20 @@ -616,6 +557,7 @@ void GDALRegister_S102() poDriver->SetDescription("S102"); poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES"); + poDriver->SetMetadataItem(GDAL_DCAP_MULTIDIM_RASTER, "YES"); poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "S-102 Bathymetric Surface Product"); poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/s102.html");