update MDAL 0.3.1

external/mdal/frmts/mdal_2dm.cpp

@@ -13,6 +13,7 @@
 #include <map>
 #include <cassert>
 #include <limits>
+#include <algorithm>
 
 #include "mdal_2dm.hpp"
 #include "mdal.h"
@@ -64,6 +65,19 @@ 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,
           "2DM Mesh File",
@@ -146,6 +160,7 @@ 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 ) )
   {
@@ -205,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];

external/mdal/frmts/mdal_2dm.hpp

@@ -29,21 +29,26 @@ 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 and BASEMENET solvers
+   * 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)
@@ -60,6 +65,13 @@ namespace MDAL
    *
    * 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
   {

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,7 +336,10 @@ 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 );
@@ -324,9 +347,11 @@ void MDAL::DriverAsciiDat::readVertexTimestep(
 
     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 )
     {

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,

external/mdal/frmts/mdal_gdal.cpp

@@ -317,7 +317,28 @@ void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
 {
   assert( raster_band );
 
-  double nodata =  GDALGetRasterNoDataValue( raster_band, nullptr );
+  // nodata
+  int pbSuccess;
+  double nodata =  GDALGetRasterNoDataValue( raster_band, &pbSuccess );
+  if ( pbSuccess == 0 ) nodata = std::numeric_limits<double>::quiet_NaN();
+  bool hasNoData = !std::isnan( nodata );
+
+  // offset and scale
+  double offset = 0.0;
+  double scale = GDALGetRasterScale( raster_band, &pbSuccess );
+  if ( ( pbSuccess == 0 ) || MDAL::equals( scale, 0.0 ) || std::isnan( scale ) )
+  {
+    scale = 1.0;
+  }
+  else
+  {
+    offset = GDALGetRasterOffset( raster_band, &pbSuccess );
+    if ( ( pbSuccess == 0 ) || std::isnan( offset ) )
+    {
+      offset = 0.0;
+    }
+  }
+
   unsigned int mXSize = meshGDALDataset()->mXSize;
   unsigned int mYSize = meshGDALDataset()->mYSize;
 
@@ -349,8 +370,11 @@ void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
     {
       unsigned int idx = x + mXSize * y;
       double val = mPafScanline[x];
-      if ( !MDAL::equals( val, nodata ) )
+      if ( !hasNoData || !MDAL::equals( val, nodata ) )
       {
+        // Apply scale and offset
+        val = val * scale + offset;
+
         // values is prepolulated with NODATA values, so store only legal values
         if ( is_vector )
         {

external/mdal/frmts/mdal_xdmf.cpp

@@ -153,7 +153,7 @@ size_t MDAL::XdmfFunctionDataset::scalarData( size_t indexStart, size_t count, d
   assert( mType != FunctionType::Join );
 
   if ( mType == FunctionType::Subtract )
-    return substractFunction( indexStart, count, buffer );
+    return subtractFunction( indexStart, count, buffer );
 
   if ( mType == FunctionType::Flow )
     return flowFunction( indexStart, count, buffer );
@@ -176,7 +176,7 @@ size_t MDAL::XdmfFunctionDataset::activeData( size_t indexStart, size_t count, i
   return count;
 }
 
-size_t MDAL::XdmfFunctionDataset::substractFunction( size_t indexStart, size_t count, double *buffer )
+size_t MDAL::XdmfFunctionDataset::subtractFunction( size_t indexStart, size_t count, double *buffer )
 {
   std::vector<double> buf( 2 * count, std::numeric_limits<double>::quiet_NaN() );
   size_t copyVals = extractRawData( indexStart, count, 2, buf );

external/mdal/frmts/mdal_xdmf.hpp

@@ -86,7 +86,7 @@ namespace MDAL
    * Currently we do not use any fancy bison/flex based
    * expression parsing, just supporting few types of
    * most common function types:
-   *   - substraction (A-B)
+   *   - subtraction (A-B)
    *   - join ( [A, B] vector)
    *   - magnitude
    *
@@ -126,7 +126,7 @@ namespace MDAL
       size_t activeData( size_t indexStart, size_t count, int *buffer ) override;
 
     private:
-      size_t substractFunction( size_t indexStart, size_t count, double *buffer );
+      size_t subtractFunction( size_t indexStart, size_t count, double *buffer );
       size_t flowFunction( size_t indexStart, size_t count, double *buffer );
       size_t joinFunction( size_t indexStart, size_t count, double *buffer );
       size_t extractRawData( size_t indexStart, size_t count, size_t nDatasets, std::vector<double> &buf );

external/mdal/mdal.cpp

@@ -22,7 +22,7 @@ static MDAL_Status sLastStatus;
 
 const char *MDAL_Version()
 {
-  return "0.3.0";
+  return "0.3.1";
 }
 
 MDAL_Status MDAL_LastStatus()