update MDAL 0.3.1

qgis · Apr 18, 2019 · da7f655 · da7f655
1 parent 1fac5a6
commit da7f655
Show file tree

Hide file tree

Showing 24 changed files with 1,467 additions and 164 deletions.
diff --git a/external/mdal/cmake/FindNetCDF.cmake b/external/mdal/cmake/FindNetCDF.cmake
@@ -28,4 +28,3 @@ FIND_LIBRARY (NETCDF_LIBRARY
 INCLUDE (FindPackageHandleStandardArgs)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS (NetCDF
   DEFAULT_MSG NETCDF_LIBRARY NETCDF_INCLUDE_DIR)
-
diff --git a/external/mdal/cmake_templates/mdal_config.hpp.in b/external/mdal/cmake_templates/mdal_config.hpp.in
@@ -9,6 +9,8 @@
 
 #cmakedefine HAVE_NETCDF
 
+#cmakedefine HAVE_XML
+
 #endif // MDAL_CONFIG_HPP
 
 
diff --git a/external/mdal/frmts/mdal_2dm.cpp b/external/mdal/frmts/mdal_2dm.cpp
@@ -12,6 +12,8 @@
 #include <vector>
 #include <map>
 #include <cassert>
+#include <limits>
+#include <algorithm>
 
 #include "mdal_2dm.hpp"
 #include "mdal.h"
@@ -63,6 +65,18 @@ size_t MDAL::Mesh2dm::vertexIndex( size_t vertexID ) const
   return vertexID;
 }
 
+size_t MDAL::Mesh2dm::maximumVertexId() const
+{
+  size_t maxIndex = verticesCount() - 1;
+  if ( mVertexIDtoIndex.empty() )
+    return maxIndex;
+  else
+  {
+    // std::map is sorted!
+    size_t maxID = mVertexIDtoIndex.rbegin()->first;
+    return std::max( maxIndex, maxID );
+  }
+}
 
 MDAL::Driver2dm::Driver2dm():
   Driver( DRIVER_NAME,
@@ -135,6 +149,9 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
   std::vector<Vertex> vertices( vertexCount );
   std::vector<Face> faces( faceCount );
 
+  // Basement 3.x supports definition of elevation for cell centers
+  std::vector<double> elementCenteredElevation;
+
   in.clear();
   in.seekg( 0, std::ios::beg );
 
@@ -143,33 +160,44 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
   size_t faceIndex = 0;
   size_t vertexIndex = 0;
   std::map<size_t, size_t> vertexIDtoIndex;
+  size_t lastVertexID = 0;
 
   while ( std::getline( in, line ) )
   {
-    if ( startsWith( line, "E4Q" ) )
+    if ( startsWith( line, "E4Q" ) ||
+         startsWith( line, "E3T" )
+       )
     {
       chunks = split( line,  ' ' );
       assert( faceIndex < faceCount );
 
+      const size_t faceVertexCount = MDAL::toSizeT( line[1] );
+      assert( ( faceVertexCount == 3 ) || ( faceVertexCount == 4 ) );
+
       Face &face = faces[faceIndex];
-      face.resize( 4 );
+      face.resize( faceVertexCount );
+
+      // chunks format here
+      // E** id vertex_id1, vertex_id2, ... material_id (elevation - optional)
+      // vertex ids are numbered from 1
       // Right now we just store node IDs here - we will convert them to node indices afterwards
-      for ( size_t i = 0; i < 4; ++i )
-        face[i] = toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1
+      assert( chunks.size() > faceVertexCount + 1 );
 
-      faceIndex++;
-    }
-    else if ( startsWith( line, "E3T" ) )
-    {
-      chunks = split( line,  ' ' );
-      assert( faceIndex < faceCount );
+      for ( size_t i = 0; i < faceVertexCount; ++i )
+        face[i] = MDAL::toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1
 
-      Face &face = faces[faceIndex];
-      face.resize( 3 );
-      // Right now we just store node IDs here - we will convert them to node indices afterwards
-      for ( size_t i = 0; i < 3; ++i )
+      // OK, now find out if there is optional cell elevation (BASEMENT 3.x)
+      if ( chunks.size() == faceVertexCount + 4 )
       {
-        face[i] = toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1
+
+        // initialize dataset if it is still empty
+        if ( elementCenteredElevation.empty() )
+        {
+          elementCenteredElevation = std::vector<double>( faceCount, std::numeric_limits<double>::quiet_NaN() );
+        }
+
+        // add Bed Elevation (Face) value
+        elementCenteredElevation[faceIndex] = MDAL::toDouble( chunks[ faceVertexCount + 3 ] );
       }
 
       faceIndex++;
@@ -192,7 +220,22 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
     else if ( startsWith( line, "ND" ) )
     {
       chunks = split( line,  ' ' );
-      size_t nodeID = toSizeT( chunks[1] ) - 1; // 2dm is numbered from 1
+      size_t nodeID = toSizeT( chunks[1] );
+
+      if ( nodeID != 0 )
+      {
+        // specification of 2DM states that ID should be positive integer numbered from 1
+        // but it seems some formats do not respect that
+        if ( ( lastVertexID != 0 ) && ( nodeID <= lastVertexID ) )
+        {
+          // the algorithm requires that the file has NDs orderer by index
+          if ( status ) *status = MDAL_Status::Err_InvalidData;
+          return nullptr;
+        }
+        lastVertexID = nodeID;
+      }
+      nodeID -= 1; // 2dm is numbered from 1
+
       _parse_vertex_id_gaps( vertexIDtoIndex, vertexIndex, nodeID, status );
       assert( vertexIndex < vertexCount );
       Vertex &vertex = vertices[vertexIndex];
@@ -236,6 +279,10 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
   );
   mesh->faces = faces;
   mesh->vertices = vertices;
+
+  // Add Bed Elevations
+  MDAL::addFaceScalarDatasetGroup( mesh.get(), elementCenteredElevation, "Bed Elevation (Face)" );
   MDAL::addBedElevationDatasetGroup( mesh.get(), vertices );
+
   return std::unique_ptr<Mesh>( mesh.release() );
 }
diff --git a/external/mdal/frmts/mdal_2dm.hpp b/external/mdal/frmts/mdal_2dm.hpp
@@ -29,19 +29,50 @@ namespace MDAL
       ~Mesh2dm() override;
 
 
-      //! Some formats supports gaps in the vertex indexing, but we return continuos array from MDAL
-      //! for most of the formats this returns
+      //! HYDRO_AS-2D supports gaps in the vertex indexing,
+      //! but we use continuos array of vertices in MDAL
       //! \param vertexID internal index/ID of the vertex that native format uses
-      //! \returns index of the vertex in the continuous array of vertices we returned by readVertices()
+      //! \returns index of the vertex in the continuous array of vertices we returned by readVertices().
+      //!          For invalid vertexID it is returned index that is out of vertices array bounds.
       virtual size_t vertexIndex( size_t vertexID ) const;
 
+      //! Returns maximum vertex ID.
+      //! For meshes without gaps in vertex indexing, it is vertex count - 1
+      virtual size_t maximumVertexId() const;
+
     private:
-      // 2dm supports "gaps" in the mesh indexing
-      // Store only the indices that have different index and ID
-      // https://github.com/lutraconsulting/MDAL/issues/51
+      //! 2dm supports "gaps" in the mesh indexing
+      //! Store only the indices that have different index and ID
+      //! https://github.com/lutraconsulting/MDAL/issues/51
       std::map<size_t, size_t> mVertexIDtoIndex;
   };
 
+  /**
+   * 2DM format specification used in TUFLOW, HYDRO_AS-2D and BASEMENET solvers
+   * Text file format representing mesh vertices (ND) and faces (E**)
+   * ND id x y z
+   * Supports lines, triangles and polygons up to 9 vertices (implemented triangles and quads)
+   * E3T id 1 2 3 mat_id -> face type, id, vertex indices ..., material index
+   *
+   * full specification here: https://www.xmswiki.com/wiki/SMS:2D_Mesh_Files_*.2dm
+   *
+   * Exception for the official specification is for recognition of cell-centered
+   * elevation values supported by BASEMENT 3.x releases
+   * If face definition has extra column, it is parsed and recognized as
+   * elevation, e.g. format for triangle
+   * E3T id 1 2 3 mat_id elevation
+   * and added automatically as "Bed Elevation (Face)"
+   *
+   * Note that some 2dm formats do have some extra columns after mat_id column with
+   * data with unknown origin/name (e.g. tests/data/2dm/regular_grid.2dm)
+   *
+   * HYDRO_AS-2D also allows gaps in vertex indexing. In this case we support only files
+   * where the vertices are sorted by ID in the source file (limitation of the implementation)
+   *
+   * Vertex/Face IDs should be indexed from 1. We support indexing from 0 for datasets in xmdf format,
+   * but not for ascii dat format (since the loop is from 0 to maximumVertexId() which is in this case
+   * numberical_limits<size_t>::max(); (limitation of the implementation)
+   */
   class Driver2dm: public Driver
   {
     public:

diff --git a/external/mdal/frmts/mdal_ascii_dat.cpp b/external/mdal/frmts/mdal_ascii_dat.cpp
@@ -13,6 +13,7 @@
 #include <map>
 #include <cassert>
 #include <memory>
+#include <limits>
 
 #include "mdal_ascii_dat.hpp"
 #include "mdal.h"
@@ -103,7 +104,8 @@ void MDAL::DriverAsciiDat::loadOldFormat( std::ifstream &in,
     if ( cardType == "ND" && items.size() >= 2 )
     {
       size_t fileNodeCount = toSizeT( items[1] );
-      if ( mesh->verticesCount() != fileNodeCount )
+      size_t meshIdCount = maximumId( mesh ) + 1;
+      if ( meshIdCount != fileNodeCount )
         EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
     }
     else if ( cardType == "SCALAR" || cardType == "VECTOR" )
@@ -167,7 +169,8 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
     if ( cardType == "ND" && items.size() >= 2 )
     {
       size_t fileNodeCount = toSizeT( items[1] );
-      if ( mesh->verticesCount() != fileNodeCount )
+      size_t meshIdCount = maximumId( mesh ) + 1;
+      if ( meshIdCount != fileNodeCount )
         EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
     }
     else if ( cardType == "NC" && items.size() >= 2 )
@@ -247,6 +250,15 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
   }
 }
 
+size_t MDAL::DriverAsciiDat::maximumId( const MDAL::Mesh *mesh ) const
+{
+  const Mesh2dm *m2dm = dynamic_cast<const Mesh2dm *>( mesh );
+  if ( m2dm )
+    return m2dm->maximumVertexId();
+  else
+    return mesh->verticesCount() - 1;
+}
+
 /**
  * The DAT format contains "datasets" and each dataset has N-outputs. One output
  * represents data for all vertices/faces for one timestep
@@ -265,6 +277,14 @@ void MDAL::DriverAsciiDat::load( const std::string &datFile, MDAL::Mesh *mesh, M
     return;
   }
 
+  size_t mID = maximumId( mesh );
+  if ( mID == std::numeric_limits<size_t>::max() )
+  {
+    // This happens when mesh is 2DM and vertices are numbered from 0
+    if ( status ) *status = MDAL_Status::Err_IncompatibleMesh;
+    return;
+  }
+
   std::ifstream in( mDatFile, std::ifstream::in );
   std::string line;
   if ( !std::getline( in, line ) )
@@ -316,17 +336,22 @@ void MDAL::DriverAsciiDat::readVertexTimestep(
 
   const Mesh2dm *m2dm = dynamic_cast<const Mesh2dm *>( mesh );
   double *values = dataset->values();
-  for ( size_t i = 0; i < mesh->verticesCount(); ++i )
+  size_t meshIdCount = maximumId( mesh ) + 1; // these are native format indexes (IDs). For formats without gaps it equals vertex array index
+  size_t vertexCount = mesh->verticesCount();
+
+  for ( size_t id = 0; id < meshIdCount; ++id )
   {
     std::string line;
     std::getline( stream, line );
     std::vector<std::string> tsItems = split( line,  ' ' );
 
     size_t index;
     if ( m2dm )
-      index = m2dm->vertexIndex( i );
+      index = m2dm->vertexIndex( id ); //this index may be out of values array
     else
-      index = i;
+      index = id;
+
+    if ( index >= vertexCount ) continue;
 
     if ( isVector )
     {

diff --git a/external/mdal/frmts/mdal_ascii_dat.hpp b/external/mdal/frmts/mdal_ascii_dat.hpp
@@ -35,6 +35,9 @@ namespace MDAL
    * such dataset, the dataset name contains "_els_" substring
    * (e.g. depth_els_1.dat)
    *
+   * HYDRO_AS-2D solver can have mesh that has numbering gaps, but
+   * speficies values for even missing indexes in dataset file
+   *
    * In one file, there is always one dataset group stored.
    *
    * Sometime the "older" datasets may have some part of the
@@ -60,6 +63,12 @@ namespace MDAL
       void loadOldFormat( std::ifstream &in, Mesh *mesh, MDAL_Status *status ) const;
       void loadNewFormat( std::ifstream &in, Mesh *mesh, MDAL_Status *status ) const;
 
+      //! Gets maximum (native) index.
+      //! For meshes without indexing gap it is vertexCount - 1
+      //! For some HYDRO_AS-2D meshes with indexing gaps, it returns
+      //! maximum native index of the vertex in defined in the mesh
+      size_t maximumId( const Mesh *mesh ) const;
+
       void readVertexTimestep(
         const Mesh *mesh,
         std::shared_ptr<DatasetGroup> group,

diff --git a/external/mdal/frmts/mdal_cf.cpp b/external/mdal/frmts/mdal_cf.cpp
@@ -125,10 +125,7 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
   }
   while ( true );
 
-  if ( dsinfo_map.size() == 0 )
-  {
-    throw MDAL_Status::Err_InvalidData;
-  }
+  if ( dsinfo_map.size() == 0 ) throw MDAL_Status::Err_InvalidData;
 
   return dsinfo_map;
 }
@@ -388,16 +385,6 @@ void MDAL::CFDimensions::setDimension( MDAL::CFDimensions::Type type,
   mCount[type] = count;
 }
 
-size_t MDAL::CFDimensions::faceCount() const
-{
-  return size( Face2D ) + size( Line1D );
-}
-
-size_t MDAL::CFDimensions::vertexCount() const
-{
-  return size( Vertex1D ) + size( Vertex2D );
-}
-
 bool MDAL::CFDimensions::isDatasetType( MDAL::CFDimensions::Type type ) const
 {
   return (

diff --git a/external/mdal/frmts/mdal_cf.hpp b/external/mdal/frmts/mdal_cf.hpp
@@ -41,10 +41,6 @@ namespace MDAL
       size_t size( Type type ) const;
       //! Sets a dimension
       void setDimension( Type type, size_t count, int ncid = -1 );
-      //! Returns number of faces (1D + 2D)
-      size_t faceCount() const;
-      //! Returns number of vertices (1D + 2D)
-      size_t vertexCount() const;
       //! Returns whether the type is one that case be used for datasets definition
       bool isDatasetType( Type type ) const;
 

diff --git a/external/mdal/frmts/mdal_flo2d.cpp b/external/mdal/frmts/mdal_flo2d.cpp
@@ -502,23 +502,6 @@ void MDAL::DriverFlo2D::createMesh( const std::vector<CellCenter> &cells, double
   mMesh->vertices = vertices;
 }
 
-bool MDAL::DriverFlo2D::isFlo2DFile( const std::string &fileName )
-{
-  std::vector<std::string> required_files =
-  {
-    "CADPTS.DAT",
-    "FPLAIN.DAT"
-  };
-
-  for ( const std::string &str : required_files )
-  {
-    std::string fn( fileNameFromDir( fileName, str ) );
-    if ( !fileExists( fn ) )
-      return false;
-  }
-  return true;
-}
-
 bool MDAL::DriverFlo2D::parseHDF5Datasets( const std::string &datFileName )
 {
   //return true on error

diff --git a/external/mdal/frmts/mdal_flo2d.hpp b/external/mdal/frmts/mdal_flo2d.hpp
@@ -25,8 +25,6 @@ namespace MDAL
       bool canRead( const std::string &uri ) override;
       std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override;
 
-      static bool isFlo2DFile( const std::string &fileName );
-
     private:
       struct CellCenter
       {