XMDF fixes (#492)

* rootGroup might contain dataset, load it if it is there * read mesh from xmdf * add scalar dataset group * new tests for h5 * add info * use buildUri and findMeshesNames to identify existing groups with Mesh Data ... in the same way UGRID driver does that * fix test data * file with multiple mesh datasets * add h5 as extension for XMDF * only store faces that have more than 2 vertices * fix messages Co-authored-by: Stefanos Natsis <[email protected]> * vertexDims size check Co-authored-by: Stefanos Natsis <[email protected]>
lutraconsulting · Nov 6, 2024 · e781a53 · e781a53
1 parent 9cf7a70
commit e781a53
Show file tree

Hide file tree

Showing 8 changed files with 409 additions and 3 deletions.
diff --git a/README.md b/README.md
@@ -125,7 +125,9 @@ QGIS contains internal copy of MDAL library in following versions:
 | 3.30.0  | 1.0.2   | |
 | 3.36.0  | 1.1.0   | Mike21 format support read/write |
 | 3.38.0  | 1.2.0   | Groundwater / surface water meshes for 3Di format |
-| 3.42.0  | 1.3.0   | Fix 2dm format coordinates saving, Support Dataset Group Removal From Mesh |
+| 3.42.0  | 1.3.0   | Fix 2dm format coordinates saving |
+|         |         | XMDF loading Dataset Group fix, support for Mesh in XMDF files (as 2DMeshModule) |
+|         |         | Support Dataset Group Removal From Mesh |
 
 versions `X.Y.9Z` are development versions or alpha/beta releases (e.g. `0.4.90`, `0.4.91`, ...)
 

diff --git a/docs/source/drivers/xmdf.rst b/docs/source/drivers/xmdf.rst
@@ -10,5 +10,8 @@ XMDF -- eXtensible Model Data Format
 
 MDAL supports reading of the XMDF format generated by TUFLOW_, HYDRO_AS-2D_ and other hydraulic modelling software applications.
 
+Since 1.3 MDAL supports reading mesh stored in XMDF format according to specification - http://xmdf.aquaveo.com/doc2.2/html/modules.html. 
+
 .. _TUFLOW: https://www.tuflow.com/
 .. _HYDRO_AS-2D : https://www.hydroas-2d.com/
+.. _Mesh Format: http://xmdf.aquaveo.com/doc2.2/html/group__d2d4d1d.html
diff --git a/mdal/frmts/mdal_xmdf.cpp b/mdal/frmts/mdal_xmdf.cpp
@@ -92,8 +92,8 @@ size_t MDAL::XmdfDataset::activeData( size_t indexStart, size_t count, int *buff
 MDAL::DriverXmdf::DriverXmdf()
   : Driver( "XMDF",
             "TUFLOW XMDF",
-            "*.xmdf",
-            Capability::ReadDatasets )
+            "*.xmdf;;*.h5",
+            Capability::ReadDatasets | Capability::ReadMesh )
 {
 }
 
@@ -206,6 +206,17 @@ void MDAL::DriverXmdf::addDatasetGroupsFromXmdfGroup( DatasetGroups &groups,
     size_t vertexCount,
     size_t faceCount ) const
 {
+  // check if this root group can be loaded as a dataset group and if so, then load it
+  std::vector<std::string> gDataNames = rootGroup.datasets();
+  if ( MDAL::contains( gDataNames, "Times" ) &&
+       MDAL::contains( gDataNames, "Values" ) &&
+       MDAL::contains( gDataNames, "Mins" ) &&
+       MDAL::contains( gDataNames, "Maxs" ) )
+  {
+    std::shared_ptr<DatasetGroup> ds = readXmdfGroupAsDatasetGroup( rootGroup, rootGroup.name() + nameSuffix, vertexCount, faceCount );
+    groups.push_back( ds );
+  }
+
   for ( const std::string &groupName : rootGroup.groups() )
   {
     HdfGroup g = rootGroup.group( groupName );
@@ -331,3 +342,241 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXmdf::readXmdfGroupAsDatasetGrou
 
   return group;
 }
+
+bool MDAL::DriverXmdf::canReadMesh( const std::string &uri )
+{
+  HdfFile file( uri, HdfFile::ReadOnly );
+  if ( !file.isValid() )
+  {
+    return false;
+  }
+
+  HdfDataset dsFileType = file.dataset( "/File Type" );
+  if ( dsFileType.readString() != "Xmdf" )
+  {
+    return false;
+  }
+
+  std::vector<std::string> meshPaths = meshGroupPaths( file );
+
+  return !meshPaths.empty();
+}
+
+std::string MDAL::DriverXmdf::buildUri( const std::string &meshFile )
+{
+  mDatFile = meshFile;
+
+  std::vector<std::string> meshNames = findMeshesNames();
+
+  return MDAL::buildAndMergeMeshUris( meshFile, meshNames, name() );
+}
+
+std::vector<std::string> MDAL::DriverXmdf::findMeshesNames() const
+{
+  std::vector<std::string> meshesInFile;
+
+  HdfFile file( mDatFile, HdfFile::ReadOnly );
+  if ( !file.isValid() )
+  {
+    return meshesInFile;
+  }
+
+  meshesInFile = meshGroupPaths( file );
+
+  return meshesInFile;
+}
+
+std::vector<std::string> MDAL::DriverXmdf::meshGroupPaths( const HdfFile &file ) const
+{
+  std::vector<std::string> meshPaths;
+
+  std::vector<std::string> rootGroups = file.groups();
+
+  for ( const std::string &groupName : rootGroups )
+  {
+    HdfGroup g = file.group( groupName );
+    std::vector<std::string> paths = meshGroupPaths( g );
+    meshPaths.insert( meshPaths.end(), paths.begin(), paths.end() );
+  }
+
+  return meshPaths;
+}
+
+std::vector<std::string> MDAL::DriverXmdf::meshGroupPaths( const HdfGroup &group ) const
+{
+  std::vector<std::string> meshPaths;
+
+  std::vector<std::string> gDataNames = group.groups();
+
+  if ( MDAL::contains( gDataNames, "Nodes" ) ||
+       MDAL::contains( gDataNames, "Elements" ) )
+  {
+    meshPaths.push_back( group.name() );
+  }
+
+  for ( const std::string &groupName : gDataNames )
+  {
+    HdfGroup g = group.group( groupName );
+    std::vector<std::string> paths = meshGroupPaths( g );
+    meshPaths.insert( meshPaths.end(), paths.begin(), paths.end() );
+  }
+
+  return meshPaths;
+}
+
+std::unique_ptr< MDAL::Mesh > MDAL::DriverXmdf::load( const std::string &meshFile, const std::string &meshName )
+{
+  mDatFile = meshFile;
+
+  MDAL::Log::resetLastStatus();
+
+  HdfFile file( mDatFile, HdfFile::ReadOnly );
+  if ( !file.isValid() )
+  {
+    MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), "File " + mDatFile + " is not valid" );
+    return nullptr;
+  }
+
+  HdfDataset dsFileType = file.dataset( "/File Type" );
+  if ( dsFileType.readString() != "Xmdf" )
+  {
+    MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), "Unknown dataset file type" );
+    return nullptr;
+  }
+
+  std::vector<std::string> meshNames = findMeshesNames();
+
+  if ( meshNames.empty() )
+  {
+    MDAL::Log::error( MDAL_Status::Err_IncompatibleMesh, name(), "No meshes found in file " + mDatFile );
+    return nullptr;
+  }
+
+  std::string meshNameToLoad = meshName;
+
+  if ( meshNameToLoad.empty() )
+  {
+    meshNameToLoad = meshNames[0];
+  }
+
+  if ( !MDAL::contains( meshNames, meshNameToLoad ) )
+  {
+    MDAL::Log::error( MDAL_Status::Err_IncompatibleMesh, name(), "No meshes with name " + meshNameToLoad + " found in file " + mDatFile );
+    return nullptr;
+  }
+
+  HdfGroup groupMeshModule = file.group( meshNameToLoad );
+
+  std::vector<std::string> gDataNames = groupMeshModule.groups();
+
+  HdfGroup gNodes = groupMeshModule.group( "Nodes" );
+
+  std::vector<std::string> namesNodes = gNodes.datasets();
+  HdfDataset nodes = gNodes.dataset( namesNodes[0] );
+
+  std::vector<hsize_t> nodesDims = nodes.dims();
+  hsize_t nodesRows = nodesDims[0];
+  size_t vertexDims = nodesDims[1];
+
+  if ( vertexDims < 2 || vertexDims > 3 )
+  {
+    MDAL::Log::error( MDAL_Status::Err_IncompatibleMesh, name(), "Vertices have unsupported number of dimensions " + std::to_string( vertexDims ) + " only 2 (X,Y) or 3 (X, Y, Z) dimensions are allowed." );
+    return nullptr;
+  }
+
+  std::vector<double> nodesData = nodes.readArrayDouble();
+
+  Vertices vertices( nodesRows );
+
+  size_t currentVertexIndex = 0;
+  size_t i = 0;
+  while ( i < nodesData.size() )
+  {
+    Vertex &vertex = vertices[currentVertexIndex];
+
+    vertex.x = nodesData[i];
+    i++;
+    vertex.y = nodesData[i];
+    i++;
+    if ( vertexDims == 3 )
+    {
+      vertex.z = nodesData[i];
+      i++;
+    }
+    currentVertexIndex++;
+  }
+
+  nodesData.clear();
+
+  HdfGroup gElements = groupMeshModule.group( "Elements" );
+
+  std::vector<std::string> namesElements = gElements.datasets();
+  HdfDataset elements = gElements.dataset( namesElements[0] );
+
+  std::vector<hsize_t> elementsDims = elements.dims();
+  hsize_t elementsRows = elementsDims[0];
+  int elementsRowsDims = elementsDims[1];
+
+  std::vector<int> facesData = elements.readArrayInt();
+
+  Faces faces( elementsRows );
+  int maxVerticesPerFace = 0;
+
+  size_t currentFaceIndex = 0;
+  i = 0;
+  while ( i < facesData.size() )
+  {
+    std::vector<size_t> tempFace;
+    for ( int j = 0; j < elementsRowsDims; j++ )
+    {
+      int vertexIndex = facesData[i];
+      if ( vertexIndex > 0 )
+      {
+        // XMDF is 1-based, MDAL is 0-based
+        tempFace.push_back( facesData[i] - 1 );
+      }
+      i++;
+    }
+
+    // only store faces with more than 2 vertices
+    if ( tempFace.size() > 2 )
+    {
+      Face &face = faces[currentFaceIndex];
+      std::copy( tempFace.begin(), tempFace.end(), std::back_inserter( face ) );
+
+      if ( tempFace.size() > maxVerticesPerFace )
+      {
+        maxVerticesPerFace = tempFace.size();
+      }
+
+      currentFaceIndex++;
+    }
+  }
+
+  facesData.clear();
+
+  // copy only the faces that have been properly filled
+  faces = Faces( faces.begin(), faces.begin() + currentFaceIndex );
+
+  // create the mesh and set the required data
+  std::unique_ptr< MemoryMesh > mesh(
+    new MemoryMesh(
+      name(),
+      maxVerticesPerFace,
+      mDatFile
+    )
+  );
+
+  std::vector<double> values( vertices.size() );
+  for ( size_t i = 0; i < vertices.size(); ++i )
+  {
+    values[i] = vertices[i].z;
+  }
+
+  mesh->setFaces( std::move( faces ) );
+  mesh->setVertices( std::move( vertices ) );
+
+  addVertexScalarDatasetGroup( mesh.get(), values, "Z-Values" );
+
+  return mesh;
+}
diff --git a/mdal/frmts/mdal_xmdf.hpp b/mdal/frmts/mdal_xmdf.hpp
@@ -89,6 +89,9 @@ namespace MDAL
       bool canReadDatasets( const std::string &uri ) override;
       void load( const std::string &datFile, Mesh *mesh ) override;
 
+      bool canReadMesh( const std::string &uri ) override;
+      std::unique_ptr< Mesh > load( const std::string &meshFile, const std::string &meshName = "" ) override;
+
     private:
       MDAL::Mesh *mMesh = nullptr;
       std::string mDatFile;
@@ -111,6 +114,11 @@ namespace MDAL
                            size_t vertexCount,
                            size_t faceCount ) const;
 
+      std::string buildUri( const std::string &meshFile ) override;
+      std::vector<std::string> findMeshesNames() const;
+
+      std::vector<std::string> meshGroupPaths( const HdfGroup &group ) const;
+      std::vector<std::string> meshGroupPaths( const HdfFile &file ) const;
   };
 
 } // namespace MDAL

diff --git a/tests/data/xmdf/withMesh/data.h5 b/tests/data/xmdf/withMesh/data.h5
diff --git a/tests/data/xmdf/withMesh/mesh.h5 b/tests/data/xmdf/withMesh/mesh.h5
diff --git a/tests/data/xmdf/withMesh/multiple_meshes.h5 b/tests/data/xmdf/withMesh/multiple_meshes.h5