Merge pull request #124 from lutraconsulting/fix_flo2d_max

fix #123 maximum datasets for flo-2d
lutraconsulting · May 27, 2019 · c37947d · c37947d
2 parents d837b4e + f4cfde5
commit c37947d
Showing 5 changed files with 68 additions and 46 deletions.
diff --git a/mdal/frmts/mdal_flo2d.cpp b/mdal/frmts/mdal_flo2d.cpp
@@ -34,11 +34,6 @@ struct VertexCompare
   }
 };
 
-static inline bool is_nodata( double val )
-{
-  return MDAL::equals( val, -9999.0, 1e-8 );
-}
-
 static std::string fileNameFromDir( const std::string &mainFileName, const std::string &name )
 {
   std::string dir = MDAL::dirName( mainFileName );
@@ -47,7 +42,7 @@ static std::string fileNameFromDir( const std::string &mainFileName, const std::
 
 static double getDouble( double val )
 {
-  if ( is_nodata( val ) )
+  if ( MDAL::equals( val, 0.0, 1e-8 ) )
   {
     return MDAL_NAN;
   }
@@ -59,19 +54,11 @@ static double getDouble( double val )
 
 static double getDouble( const std::string &val )
 {
-  if ( MDAL::isNumber( val ) )
-  {
-    double valF = MDAL::toDouble( val );
-    return getDouble( valF );
-  }
-  else
-  {
-    return MDAL_NAN;
-  }
+  double valF = MDAL::toDouble( val );
+  return getDouble( valF );
 }
 
 void MDAL::DriverFlo2D::addStaticDataset(
-  bool isOnVertices,
   std::vector<double> &vals,
   const std::string &groupName,
   const std::string &datFileName )
@@ -82,7 +69,7 @@ void MDAL::DriverFlo2D::addStaticDataset(
                                           datFileName,
                                           groupName
                                         );
-  group->setIsOnVertices( isOnVertices );
+  group->setIsOnVertices( false );
   group->setIsScalar( true );
 
   std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
@@ -257,7 +244,7 @@ void MDAL::DriverFlo2D::parseTIMDEPFile( const std::string &datFileName, const s
       double depth = getDouble( lineParts[1] );
       depthDataset->values()[face_idx] = depth;
 
-      if ( !is_nodata( depth ) ) depth += elevations[face_idx];
+      if ( !std::isnan( depth ) ) depth += elevations[face_idx];
       waterLevelDataset->values()[face_idx] = depth;
 
       face_idx ++;
@@ -295,17 +282,17 @@ void MDAL::DriverFlo2D::parseDEPTHFile( const std::string &datFileName, const st
   std::ifstream depthStream( depthFile, std::ifstream::in );
   std::string line;
 
-  size_t nVertices = mMesh->verticesCount();
-  std::vector<double> maxDepth( nVertices );
-  std::vector<double> maxWaterLevel( nVertices );
+  size_t nFaces = mMesh->facesCount();
+  std::vector<double> maxDepth( nFaces );
+  std::vector<double> maxWaterLevel( nFaces );
 
   size_t vertex_idx = 0;
 
   // DEPTH.OUT - COORDINATES (ELEM NUM, X, Y, MAX DEPTH)
   while ( std::getline( depthStream, line ) )
   {
     line = MDAL::rtrim( line );
-    if ( vertex_idx == nVertices ) throw MDAL_Status::Err_IncompatibleMesh;
+    if ( vertex_idx == nFaces ) throw MDAL_Status::Err_IncompatibleMesh;
 
     std::vector<std::string> lineParts = MDAL::split( line, ' ' );
     if ( lineParts.size() != 4 )
@@ -317,23 +304,23 @@ void MDAL::DriverFlo2D::parseDEPTHFile( const std::string &datFileName, const st
     maxDepth[vertex_idx] = val;
 
     //water level
-    if ( !is_nodata( val ) ) val += elevations[vertex_idx];
+    if ( !std::isnan( val ) ) val += elevations[vertex_idx];
     maxWaterLevel[vertex_idx] = val;
 
 
     vertex_idx++;
   }
 
-  addStaticDataset( true, maxDepth, "Depth/Maximums", datFileName );
-  addStaticDataset( true, maxWaterLevel, "Water Level/Maximums", datFileName );
+  addStaticDataset( maxDepth, "Depth/Maximums", datFileName );
+  addStaticDataset( maxWaterLevel, "Water Level/Maximums", datFileName );
 }
 
 
 void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
 {
   // these files are optional, so if not present, reading is skipped
-  size_t nVertices = mMesh->verticesCount();
-  std::vector<double> maxVel( nVertices );
+  size_t nFaces = mMesh->facesCount();
+  std::vector<double> maxVel( nFaces );
 
   {
     std::string velocityFile( fileNameFromDir( datFileName, "VELFP.OUT" ) );
@@ -350,7 +337,7 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
     // VELFP.OUT - COORDINATES (ELEM NUM, X, Y, MAX VEL) - Maximum floodplain flow velocity;
     while ( std::getline( velocityStream, line ) )
     {
-      if ( vertex_idx == nVertices ) throw MDAL_Status::Err_IncompatibleMesh;
+      if ( vertex_idx == nFaces ) throw MDAL_Status::Err_IncompatibleMesh;
 
       line = MDAL::rtrim( line );
       std::vector<std::string> lineParts = MDAL::split( line, ' ' );
@@ -381,7 +368,7 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
     // VELOC.OUT - COORDINATES (ELEM NUM, X, Y, MAX VEL)  - Maximum channel flow velocity
     while ( std::getline( velocityStream, line ) )
     {
-      if ( vertex_idx == nVertices ) throw MDAL_Status::Err_IncompatibleMesh;
+      if ( vertex_idx == nFaces ) throw MDAL_Status::Err_IncompatibleMesh;
 
       line = MDAL::rtrim( line );
       std::vector<std::string> lineParts = MDAL::split( line, ' ' );
@@ -391,7 +378,7 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
       }
 
       double val = getDouble( lineParts[3] );
-      if ( !is_nodata( val ) ) // overwrite value from VELFP if it is not 0
+      if ( !std::isnan( val ) )  // overwrite value from VELFP if it is not 0
       {
         maxVel[vertex_idx] = val;
       }
@@ -400,7 +387,7 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
     }
   }
 
-  addStaticDataset( true, maxVel, "Velocity/Maximums", datFileName );
+  addStaticDataset( maxVel, "Velocity/Maximums", datFileName );
 }
 
 double MDAL::DriverFlo2D::calcCellSize( const std::vector<CellCenter> &cells )
@@ -671,7 +658,7 @@ std::unique_ptr< MDAL::Mesh > MDAL::DriverFlo2D::load( const std::string &result
     createMesh( cells, cell_size / 2.0 );
 
     // create output for bed elevation
-    addStaticDataset( false, elevations, "Bed Elevation", mDatFileName );
+    addStaticDataset( elevations, "Bed Elevation", mDatFileName );
 
     if ( parseHDF5Datasets( mDatFileName ) )
     {

diff --git a/mdal/frmts/mdal_flo2d.hpp b/mdal/frmts/mdal_flo2d.hpp
@@ -45,7 +45,7 @@ namespace MDAL
       void parseTIMDEPFile( const std::string &datFileName, const std::vector<double> &elevations );
       void parseFPLAINFile( std::vector<double> &elevations, const std::string &datFileName, std::vector<CellCenter> &cells );
       void parseCADPTSFile( const std::string &datFileName, std::vector<CellCenter> &cells );
-      void addStaticDataset( bool isOnVertices, std::vector<double> &vals, const std::string &groupName, const std::string &datFileName );
+      void addStaticDataset( std::vector<double> &vals, const std::string &groupName, const std::string &datFileName );
       static MDAL::Vertex createVertex( size_t position, double half_cell_size, const CellCenter &cell );
       static double calcCellSize( const std::vector<CellCenter> &cells );
   };

diff --git a/mdal/mdal_utils.cpp b/mdal/mdal_utils.cpp
@@ -113,12 +113,6 @@ double MDAL::toDouble( const std::string &str )
   return atof( str.c_str() );
 }
 
-bool MDAL::isNumber( const std::string &str )
-{
-  // https://stackoverflow.com/a/16465826/2838364
-  return ( strspn( str.c_str(), "-.0123456789" ) == str.size() );
-}
-
 int MDAL::toInt( const std::string &str )
 {
   return atoi( str.c_str() );

diff --git a/mdal/mdal_utils.hpp b/mdal/mdal_utils.hpp
@@ -65,7 +65,6 @@ namespace MDAL
   int toInt( const std::string &str );
   double toDouble( const std::string &str );
   bool toBool( const std::string &str );
-  bool isNumber( const std::string &str );
 
   /**
    * Splits by deliminer and skips empty parts.

diff --git a/tests/test_flo2d.cpp b/tests/test_flo2d.cpp
@@ -162,24 +162,24 @@ TEST( MeshFlo2dTest, BarnHDF5 )
   EXPECT_EQ( false, onVertices );
 
   ASSERT_EQ( 20, MDAL_G_datasetCount( g ) );
-  ds = MDAL_G_dataset( g, 5 );
+  ds = MDAL_G_dataset( g, 17 );
   ASSERT_NE( ds, nullptr );
 
   valid = MDAL_D_isValid( ds );
   EXPECT_EQ( true, valid );
 
-  active = getActive( ds, 0 );
+  active = getActive( ds, 300 );
   EXPECT_EQ( true, active );
 
   count = MDAL_D_valueCount( ds );
   ASSERT_EQ( 521, count );
 
-  value = getValueX( ds, 0 );
-  EXPECT_DOUBLE_EQ( 0, value );
+  value = getValueX( ds, 234 );
+  EXPECT_DOUBLE_EQ( 0.27071857452392578, value );
 
   MDAL_D_minimumMaximum( ds, &min, &max );
-  EXPECT_DOUBLE_EQ( 0, min );
-  EXPECT_DOUBLE_EQ( 0, max );
+  EXPECT_DOUBLE_EQ( 0.1241119660936652, min );
+  EXPECT_DOUBLE_EQ( 2.847882132344469, max );
 
   MDAL_CloseMesh( m );
 }
@@ -321,6 +321,48 @@ TEST( MeshFlo2dTest, basic )
     double time = MDAL_D_time( ds );
     EXPECT_DOUBLE_EQ( 0.5, time );
 
+    // ///////////
+    // Max Depth
+    // ///////////
+    g = MDAL_M_datasetGroup( m, 4 );
+    ASSERT_NE( g, nullptr );
+
+    meta_count = MDAL_G_metadataCount( g );
+    ASSERT_EQ( 1, meta_count );
+
+    name = MDAL_G_name( g );
+    EXPECT_EQ( std::string( "Depth/Maximums" ), std::string( name ) );
+
+    scalar = MDAL_G_hasScalarData( g );
+    EXPECT_EQ( true, scalar );
+
+    onVertices = MDAL_G_isOnVertices( g );
+    EXPECT_EQ( false, onVertices );
+
+    ASSERT_EQ( 1, MDAL_G_datasetCount( g ) );
+    ds = MDAL_G_dataset( g, 0 );
+    ASSERT_NE( ds, nullptr );
+
+    valid = MDAL_D_isValid( ds );
+    EXPECT_EQ( true, valid );
+
+    active = getActive( ds, 0 );
+    EXPECT_EQ( true, active );
+
+    count = MDAL_D_valueCount( ds );
+    ASSERT_EQ( 9, count );
+
+    value = getValue( ds, 1 );
+    EXPECT_DOUBLE_EQ( 2, value );
+
+    MDAL_D_minimumMaximum( ds, &min, &max );
+    EXPECT_DOUBLE_EQ( 2, min );
+    EXPECT_DOUBLE_EQ( 4, max );
+
+    MDAL_G_minimumMaximum( g, &min, &max );
+    EXPECT_DOUBLE_EQ( 2, min );
+    EXPECT_DOUBLE_EQ( 4, max );
+
     MDAL_CloseMesh( m );
   }
 }