sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-10-08 00:05:09 -04:00
Code Issues Packages Projects Releases Wiki Activity

update MDAL to 0.1.4 (RC1 for QGIS 3.6)

Browse Source
This commit is contained in:
Peter Petrik 2019-01-22 10:29:53 +01:00 committed by Nyall Dawson
parent 4048d37a46
commit 928a559aa9
18 changed files with 1213 additions and 711 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/mdal/frmts/mdal_3di.cpp vendored
View File

@ -4,6 +4,8 @@
*/ */
#include "mdal_3di.hpp" #include "mdal_3di.hpp"
#include <netcdf.h>
#include <assert.h>
MDAL::Driver3Di::Driver3Di() MDAL::Driver3Di::Driver3Di()
: DriverCF( : DriverCF(

351
external/mdal/frmts/mdal_ascii_dat.cpp vendored
View File

@ -49,15 +49,203 @@ bool MDAL::DriverAsciiDat::canRead( const std::string &uri )
} }
line = trim( line ); line = trim( line );
if ( line != "DATASET" && return canReadNewFormat( line ) || canReadOldFormat( line );
line != "SCALAR" &&
line != "VECTOR" )
{
return false;
}
return true;
} }
bool MDAL::DriverAsciiDat::canReadOldFormat( const std::string &line ) const
{
return MDAL::contains( line, "SCALAR" ) ||
MDAL::contains( line, "VECTOR" ) ||
MDAL::contains( line, "TS" );
}
bool MDAL::DriverAsciiDat::canReadNewFormat( const std::string &line ) const
{
return line == "DATASET";
}
void MDAL::DriverAsciiDat::loadOldFormat( std::ifstream &in,
Mesh *mesh,
MDAL_Status *status ) const
{
std::shared_ptr<DatasetGroup> group; // DAT outputs data
std::string groupName( MDAL::baseName( mDatFile ) );
std::string line;
std::getline( in, line );
// Read the first line
bool isVector = MDAL::contains( line, "VECTOR" );
group = std::make_shared< DatasetGroup >(
name(),
mesh,
mDatFile,
groupName
);
group->setIsScalar( !isVector );
group->setIsOnVertices( true );
do
{
// Replace tabs by spaces,
// since basement v.2.8 uses tabs instead of spaces (e.g. 'TS 0\t0.0')
line = replace( line, "\t", " " );
// Trim string for cases when file has inconsistent new line symbols
line = MDAL::trim( line );
// Split to tokens
std::vector<std::string> items = split( line, " ", SplitBehaviour::SkipEmptyParts );
if ( items.size() < 1 )
continue; // empty line?? let's skip it
std::string cardType = items[0];
if ( cardType == "ND" && items.size() >= 2 )
{
size_t fileNodeCount = toSizeT( items[1] );
if ( mesh->verticesCount() != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( cardType == "SCALAR" || cardType == "VECTOR" )
{
// just ignore - we know the type from earlier...
}
else if ( cardType == "TS" && items.size() >= 2 )
{
double t = toDouble( items[ 1 ] );
readVertexTimestep( mesh, group, t, isVector, false, in );
}
else
{
std::stringstream str;
str << " Unknown card:" << line;
debug( str.str() );
}
}
while ( std::getline( in, line ) );
if ( !group || group->datasets.size() == 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
group.reset();
}
void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
Mesh *mesh,
MDAL_Status *status ) const
{
bool isVector = false;
std::shared_ptr<DatasetGroup> group; // DAT outputs data
std::string groupName( MDAL::baseName( mDatFile ) );
std::string line;
// see if it contains face-centered results - supported by BASEMENT
bool faceCentered = false;
if ( contains( groupName, "_els_" ) )
faceCentered = true;
if ( group )
group->setIsOnVertices( !faceCentered );
while ( std::getline( in, line ) )
{
// Replace tabs by spaces,
// since basement v.2.8 uses tabs instead of spaces (e.g. 'TS 0\t0.0')
line = replace( line, "\t", " " );
// Trim string for cases when file has inconsistent new line symbols
line = MDAL::trim( line );
// Split to tokens
std::vector<std::string> items = split( line, " ", SplitBehaviour::SkipEmptyParts );
if ( items.size() < 1 )
continue; // empty line?? let's skip it
std::string cardType = items[0];
if ( cardType == "ND" && items.size() >= 2 )
{
size_t fileNodeCount = toSizeT( items[1] );
if ( mesh->verticesCount() != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( cardType == "NC" && items.size() >= 2 )
{
size_t fileElemCount = toSizeT( items[1] );
if ( mesh->facesCount() != fileElemCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( cardType == "OBJTYPE" )
{
if ( items[1] != "mesh2d" && items[1] != "\"mesh2d\"" )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
else if ( cardType == "BEGSCL" || cardType == "BEGVEC" )
{
if ( group )
{
debug( "New dataset while previous one is still active!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
isVector = cardType == "BEGVEC";
group = std::make_shared< DatasetGroup >(
name(),
mesh,
mDatFile,
groupName
);
group->setIsScalar( !isVector );
group->setIsOnVertices( !faceCentered );
}
else if ( cardType == "ENDDS" )
{
if ( !group )
{
debug( "ENDDS card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
group.reset();
}
else if ( cardType == "NAME" && items.size() >= 2 )
{
if ( !group )
{
debug( "NAME card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
size_t quoteIdx1 = line.find( '\"' );
size_t quoteIdx2 = line.find( '\"', quoteIdx1 + 1 );
if ( quoteIdx1 != std::string::npos && quoteIdx2 != std::string::npos )
group->setName( line.substr( quoteIdx1 + 1, quoteIdx2 - quoteIdx1 - 1 ) );
}
else if ( cardType == "TS" && items.size() >= 3 )
{
double t = toDouble( items[2] );
if ( faceCentered )
{
readFaceTimestep( mesh, group, t, isVector, in );
}
else
{
bool hasStatus = ( toBool( items[1] ) );
readVertexTimestep( mesh, group, t, isVector, hasStatus, in );
}
}
else
{
std::stringstream str;
str << " Unknown card:" << line;
debug( str.str() );
}
}
}
/** /**
* The DAT format contains "datasets" and each dataset has N-outputs. One output * The DAT format contains "datasets" and each dataset has N-outputs. One output
@ -85,147 +273,18 @@ void MDAL::DriverAsciiDat::load( const std::string &datFile, MDAL::Mesh *mesh, M
return; return;
} }
line = trim( line ); line = trim( line );
if ( canReadNewFormat( line ) )
// http://www.xmswiki.com/xms/SMS:ASCII_Dataset_Files_*.dat
// Apart from the format specified above, there is an older supported format used in BASEMENT (and SMS?)
// which is simpler (has only one dataset in one file, no status flags etc)
bool oldFormat;
bool isVector = false;
std::shared_ptr<DatasetGroup> group; // DAT outputs data
std::string groupName( MDAL::baseName( mDatFile ) );
if ( line == "DATASET" )
oldFormat = false;
else if ( line == "SCALAR" || line == "VECTOR" )
{ {
oldFormat = true; // we do not need to parse first line again
isVector = ( line == "VECTOR" ); loadNewFormat( in, mesh, status );
group = std::make_shared< DatasetGroup >(
name(),
mesh,
mDatFile,
groupName
);
group->setIsScalar( !isVector );
} }
else else
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
// see if it contains face-centered results - supported by BASEMENT
bool faceCentered = false;
if ( !oldFormat && contains( groupName, "_els_" ) )
faceCentered = true;
if ( group )
group->setIsOnVertices( !faceCentered );
while ( std::getline( in, line ) )
{ {
// Replace tabs by spaces, // we need to parse first line again to see
// since basement v.2.8 uses tabs instead of spaces (e.g. 'TS 0\t0.0') // scalar/vector flag or timestep flag
line = replace( line, "\t", " " ); in.clear();
in.seekg( 0 );
std::vector<std::string> items = split( line, " ", SplitBehaviour::SkipEmptyParts ); loadOldFormat( in, mesh, status );
if ( items.size() < 1 )
continue; // empty line?? let's skip it
std::string cardType = items[0];
if ( cardType == "ND" && items.size() >= 2 )
{
size_t fileNodeCount = toSizeT( items[1] );
if ( mesh->verticesCount() != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( !oldFormat && cardType == "NC" && items.size() >= 2 )
{
size_t fileElemCount = toSizeT( items[1] );
if ( mesh->facesCount() != fileElemCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( !oldFormat && cardType == "OBJTYPE" )
{
if ( items[1] != "mesh2d" && items[1] != "\"mesh2d\"" )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
else if ( !oldFormat && ( cardType == "BEGSCL" || cardType == "BEGVEC" ) )
{
if ( group )
{
debug( "New dataset while previous one is still active!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
isVector = cardType == "BEGVEC";
group = std::make_shared< DatasetGroup >(
name(),
mesh,
mDatFile,
groupName
);
group->setIsScalar( !isVector );
group->setIsOnVertices( !faceCentered );
}
else if ( !oldFormat && cardType == "ENDDS" )
{
if ( !group )
{
debug( "ENDDS card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
group.reset();
}
else if ( !oldFormat && cardType == "NAME" && items.size() >= 2 )
{
if ( !group )
{
debug( "NAME card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
size_t quoteIdx1 = line.find( '\"' );
size_t quoteIdx2 = line.find( '\"', quoteIdx1 + 1 );
if ( quoteIdx1 != std::string::npos && quoteIdx2 != std::string::npos )
group->setName( line.substr( quoteIdx1 + 1, quoteIdx2 - quoteIdx1 - 1 ) );
}
else if ( oldFormat && ( cardType == "SCALAR" || cardType == "VECTOR" ) )
{
// just ignore - we know the type from earlier...
}
else if ( cardType == "TS" && items.size() >= ( oldFormat ? 2 : 3 ) )
{
double t = toDouble( items[oldFormat ? 1 : 2] );
if ( faceCentered )
{
readFaceTimestep( mesh, group, t, isVector, in );
}
else
{
bool hasStatus = ( oldFormat ? false : toBool( items[1] ) );
readVertexTimestep( mesh, group, t, isVector, hasStatus, in );
}
}
else
{
std::stringstream str;
str << " Unknown card:" << line;
debug( str.str() );
}
}
if ( oldFormat )
{
if ( !group || group->datasets.size() == 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
group.reset();
} }
} }
@ -235,7 +294,7 @@ void MDAL::DriverAsciiDat::readVertexTimestep(
double t, double t,
bool isVector, bool isVector,
bool hasStatus, bool hasStatus,
std::ifstream &stream ) std::ifstream &stream ) const
{ {
assert( group ); assert( group );
size_t faceCount = mesh->facesCount(); size_t faceCount = mesh->facesCount();
@ -301,7 +360,7 @@ void MDAL::DriverAsciiDat::readFaceTimestep(
std::shared_ptr<DatasetGroup> group, std::shared_ptr<DatasetGroup> group,
double t, double t,
bool isVector, bool isVector,
std::ifstream &stream ) std::ifstream &stream ) const
{ {
assert( group ); assert( group );

34
external/mdal/frmts/mdal_ascii_dat.hpp vendored
View File

@ -20,6 +20,29 @@
namespace MDAL namespace MDAL
{ {
/**
* ASCII Dat format is used by various solvers and the output
* from various solvers can have slightly different header.
* The format is used by TUFLOW, BASEMENT and HYDRO_AS-2D solvers.
*
* The most frequent form is based on official SMS documentation
* https://www.xmswiki.com/wiki/SMS:ASCII_Dataset_Files_*.dat
* The official format only supports data defined on vertices,
* both scalar vector. The new format is recognized by keyword
* "DATASET" on the first line of the file.
*
* BASEMENT solver also stores data defined on faces. To recognize
* such dataset, the dataset name contains "_els_" substring
* (e.g. depth_els_1.dat)
*
* In one file, there is always one dataset group stored.
*
* Sometime the "older" datasets may have some part of the
* header missing, e.g. the file starts with SCALAR or VECTOR or TS
* keyword. The older format does not have "active" flags for faces
* and does not recognize most of the keywords. Old format data
* are always defined on vertices
*/
class DriverAsciiDat: public Driver class DriverAsciiDat: public Driver
{ {
public: public:
@ -30,22 +53,27 @@ namespace MDAL
bool canRead( const std::string &uri ) override; bool canRead( const std::string &uri ) override;
void load( const std::string &datFile, Mesh *mesh, MDAL_Status *status ) override; void load( const std::string &datFile, Mesh *mesh, MDAL_Status *status ) override;
private: private:
bool canReadOldFormat( const std::string &line ) const;
bool canReadNewFormat( const std::string &line ) const;
void loadOldFormat( std::ifstream &in, Mesh *mesh, MDAL_Status *status ) const;
void loadNewFormat( std::ifstream &in, Mesh *mesh, MDAL_Status *status ) const;
void readVertexTimestep( void readVertexTimestep(
const Mesh *mesh, const Mesh *mesh,
std::shared_ptr<DatasetGroup> group, std::shared_ptr<DatasetGroup> group,
double t, double t,
bool isVector, bool isVector,
bool hasStatus, bool hasStatus,
std::ifstream &stream ); std::ifstream &stream ) const;
void readFaceTimestep( void readFaceTimestep(
const Mesh *mesh, const Mesh *mesh,
std::shared_ptr<DatasetGroup> group, std::shared_ptr<DatasetGroup> group,
double t, double t,
bool isVector, bool isVector,
std::ifstream &stream ); std::ifstream &stream ) const;
std::string mDatFile; std::string mDatFile;
}; };

7
external/mdal/frmts/mdal_cf.cpp vendored
View File

@ -5,14 +5,15 @@
#include <vector> #include <vector>
#include <string> #include <string>
#include <netcdf.h>
#include "math.h"
#include <stdlib.h>
#include <assert.h>
#include "mdal_data_model.hpp" #include "mdal_data_model.hpp"
#include "mdal_cf.hpp" #include "mdal_cf.hpp"
#include "mdal_utils.hpp" #include "mdal_utils.hpp"
#include "math.h"
#include <stdlib.h>
#define CF_THROW_ERR throw MDAL_Status::Err_UnknownFormat #define CF_THROW_ERR throw MDAL_Status::Err_UnknownFormat
MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo() MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()

67
external/mdal/frmts/mdal_gdal.cpp vendored
View File

@ -259,6 +259,60 @@ void MDAL::DriverGdal::parseRasterBands( const MDAL::GdalDataset *cfGDALDataset
} }
} }
void MDAL::DriverGdal::fixRasterBands()
{
// go through all bands and find out if both x and y dataset are valid
// if such pair if found, make the group scalar by removal of null band pointer
// this can happen is parseBandInfo() returns false-positive in vector detection
for ( data_hash::iterator band = mBands.begin(); band != mBands.end(); band++ )
{
if ( band->second.empty() )
continue;
// scalars are always ok
bool is_scalar = ( band->second.begin()->second.size() == 1 );
if ( is_scalar )
continue;
// check if we have some null bands
int needs_fix = false;
for ( timestep_map::const_iterator time_step = band->second.begin(); time_step != band->second.end(); time_step++ )
{
std::vector<GDALRasterBandH> raster_bands = time_step->second;
if ( !raster_bands[0] )
{
needs_fix = true;
break;
}
if ( !raster_bands[1] )
{
needs_fix = true;
break;
}
}
// convert this vector to scalar
if ( needs_fix )
{
for ( timestep_map::iterator time_step = band->second.begin(); time_step != band->second.end(); time_step++ )
{
std::vector<GDALRasterBandH> &raster_bands = time_step->second;
if ( !raster_bands[0] )
{
raster_bands[0] = raster_bands[1];
}
// get rid of y-coord
raster_bands.resize( 1 );
// not we should have only 1 band and valid
assert( raster_bands[0] );
}
}
}
}
void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x ) void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x )
{ {
assert( raster_band ); assert( raster_band );
@ -335,9 +389,10 @@ void MDAL::DriverGdal::activateFaces( std::shared_ptr<MemoryDataset> tos )
Face elem = mMesh->faces.at( idx ); Face elem = mMesh->faces.at( idx );
for ( size_t i = 0; i < 4; ++i ) for ( size_t i = 0; i < 4; ++i )
{ {
const size_t vertexIndex = elem[i];
if ( isScalar ) if ( isScalar )
{ {
double val = values[elem[i]]; double val = values[vertexIndex];
if ( std::isnan( val ) ) if ( std::isnan( val ) )
{ {
active[idx] = 0; //NOT ACTIVE active[idx] = 0; //NOT ACTIVE
@ -346,8 +401,8 @@ void MDAL::DriverGdal::activateFaces( std::shared_ptr<MemoryDataset> tos )
} }
else else
{ {
double x = values[elem[2 * i]]; double x = values[2 * vertexIndex];
double y = values[elem[2 * i + 1]]; double y = values[2 * vertexIndex + 1];
if ( std::isnan( x ) || std::isnan( y ) ) if ( std::isnan( x ) || std::isnan( y ) )
{ {
active[idx] = 0; //NOT ACTIVE active[idx] = 0; //NOT ACTIVE
@ -551,6 +606,12 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverGdal::load( const std::string &fileName,
} }
} }
// Fix consistency of groups
// It can happen that we thought that the
// group is vector based on name, but it could be just coicidence
// or clash in naming
fixRasterBands();
// Create MDAL datasets // Create MDAL datasets
addDatasetGroups(); addDatasetGroups();
} }

1
external/mdal/frmts/mdal_gdal.hpp vendored
View File

@ -94,6 +94,7 @@ namespace MDAL
void addDatasetGroups(); void addDatasetGroups();
void createMesh(); void createMesh();
void parseRasterBands( const GdalDataset *cfGDALDataset ); void parseRasterBands( const GdalDataset *cfGDALDataset );
void fixRasterBands();
std::string mFileName; std::string mFileName;
const std::string mGdalDriverName; /* GDAL driver name */ const std::string mGdalDriverName; /* GDAL driver name */

112
external/mdal/frmts/mdal_hec2d.cpp vendored
View File

@ -164,7 +164,7 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
{ {
for ( size_t c = 0; c < 2; ++c ) for ( size_t c = 0; c < 2; ++c )
{ {
size_t cell_idx = face2Cells[2 * i + c] + areaElemStartIndex[nArea]; size_t cell_idx = static_cast<size_t>( face2Cells[2 * i + c] ) + areaElemStartIndex[nArea];
// Take just maximum // Take just maximum
if ( std::isnan( values[cell_idx] ) || values[cell_idx] < val ) if ( std::isnan( values[cell_idx] ) || values[cell_idx] < val )
{ {
@ -356,7 +356,7 @@ void MDAL::DriverHec2D::readElemResults(
); );
} }
std::vector<std::string> read2DFlowAreasNames( HdfGroup gGeom2DFlowAreas ) std::vector<std::string> MDAL::DriverHec2D::read2DFlowAreasNamesOld( HdfGroup gGeom2DFlowAreas ) const
{ {
HdfDataset dsNames = openHdfDataset( gGeom2DFlowAreas, "Names" ); HdfDataset dsNames = openHdfDataset( gGeom2DFlowAreas, "Names" );
std::vector<std::string> names = dsNames.readArrayString(); std::vector<std::string> names = dsNames.readArrayString();
@ -367,6 +367,85 @@ std::vector<std::string> read2DFlowAreasNames( HdfGroup gGeom2DFlowAreas )
return names; return names;
} }
/**
For 5.0.5+ format
DATATYPE H5T_COMPOUND {
H5T_STRING {
STRSIZE 16;
STRPAD H5T_STR_NULLTERM;
CSET H5T_CSET_ASCII;
CTYPE H5T_C_S1;
} "Name";
H5T_IEEE_F32LE "Mann";
H5T_IEEE_F32LE "Cell Vol Tol";
H5T_IEEE_F32LE "Cell Min Area Fraction";
H5T_IEEE_F32LE "Face Profile Tol";
H5T_IEEE_F32LE "Face Area Tol";
H5T_IEEE_F32LE "Face Conv Ratio";
H5T_IEEE_F32LE "Laminar Depth";
H5T_IEEE_F32LE "Spacing dx";
H5T_IEEE_F32LE "Spacing dy";
H5T_IEEE_F32LE "Shift dx";
H5T_IEEE_F32LE "Shift dy";
H5T_STD_I32LE "Cell Count";
}
*/
typedef struct FlowAreasAttribute505
{
char name[HDF_MAX_NAME];
float mann;
float cellVolTol;
float cellMinAreaFraction;
float faceProfileTol;
float faceAreaTol;
float faceConvRatio;
float laminarDepth;
float spacingDx;
float spacingDy;
float shifyDx;
float shifyDy;
int cellCount;
} FlowAreasAttribute505;
std::vector<std::string> MDAL::DriverHec2D::read2DFlowAreasNames505( HdfGroup gGeom2DFlowAreas ) const
{
HdfDataset dsAttributes = openHdfDataset( gGeom2DFlowAreas, "Attributes" );
hid_t attributeHID = H5Tcreate( H5T_COMPOUND, sizeof( FlowAreasAttribute505 ) );
hid_t stringHID = H5Tcopy( H5T_C_S1 );
H5Tset_size( stringHID, HDF_MAX_NAME );
H5Tinsert( attributeHID, "Name", HOFFSET( FlowAreasAttribute505, name ), stringHID );
H5Tinsert( attributeHID, "Mann", HOFFSET( FlowAreasAttribute505, mann ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Cell Vol Tol", HOFFSET( FlowAreasAttribute505, cellVolTol ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Cell Min Area Fraction", HOFFSET( FlowAreasAttribute505, cellMinAreaFraction ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Face Profile Tol", HOFFSET( FlowAreasAttribute505, faceProfileTol ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Face Area Tol", HOFFSET( FlowAreasAttribute505, faceAreaTol ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Face Conv Ratio", HOFFSET( FlowAreasAttribute505, faceConvRatio ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Laminar Depth", HOFFSET( FlowAreasAttribute505, laminarDepth ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Spacing dx", HOFFSET( FlowAreasAttribute505, spacingDx ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Spacing dy", HOFFSET( FlowAreasAttribute505, spacingDy ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Shift dx", HOFFSET( FlowAreasAttribute505, shifyDx ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Shift dy", HOFFSET( FlowAreasAttribute505, shifyDy ), H5T_NATIVE_FLOAT );
H5Tinsert( attributeHID, "Cell Count", HOFFSET( FlowAreasAttribute505, cellCount ), H5T_NATIVE_INT );
std::vector<FlowAreasAttribute505> attributes = dsAttributes.readArray<FlowAreasAttribute505>( attributeHID );
H5Tclose( attributeHID );
H5Tclose( stringHID );
std::vector<std::string> names;
if ( attributes.empty() )
{
throw MDAL_Status::Err_InvalidData;
}
for ( const auto &attr : attributes )
{
std::string dat = std::string( attr.name );
names.push_back( MDAL::trim( dat ) );
}
return names;
}
void MDAL::DriverHec2D::setProjection( HdfFile hdfFile ) void MDAL::DriverHec2D::setProjection( HdfFile hdfFile )
{ {
try try
@ -474,16 +553,22 @@ bool MDAL::DriverHec2D::canRead( const std::string &uri )
{ {
HdfFile hdfFile = openHdfFile( uri ); HdfFile hdfFile = openHdfFile( uri );
std::string fileType = openHdfAttribute( hdfFile, "File Type" ); std::string fileType = openHdfAttribute( hdfFile, "File Type" );
if ( fileType != "HEC-RAS Results" ) return canReadOldFormat( fileType ) || canReadFormat505( fileType );
{
return false;
}
} }
catch ( MDAL_Status ) catch ( MDAL_Status )
{ {
return false; return false;
} }
return true; }
bool MDAL::DriverHec2D::canReadOldFormat( const std::string &fileType ) const
{
return fileType == "HEC-RAS Results";
}
bool MDAL::DriverHec2D::canReadFormat505( const std::string &fileType ) const
{
return fileType == "HEC-RAS Geometry";
} }
std::unique_ptr<MDAL::Mesh> MDAL::DriverHec2D::load( const std::string &resultsFile, MDAL_Status *status ) std::unique_ptr<MDAL::Mesh> MDAL::DriverHec2D::load( const std::string &resultsFile, MDAL_Status *status )
@ -498,15 +583,17 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverHec2D::load( const std::string &resultsF
// Verify it is correct file // Verify it is correct file
std::string fileType = openHdfAttribute( hdfFile, "File Type" ); std::string fileType = openHdfAttribute( hdfFile, "File Type" );
if ( fileType != "HEC-RAS Results" ) bool oldFormat = canReadOldFormat( fileType );
{
throw MDAL_Status::Err_UnknownFormat;
}
HdfGroup gGeom = openHdfGroup( hdfFile, "Geometry" ); HdfGroup gGeom = openHdfGroup( hdfFile, "Geometry" );
HdfGroup gGeom2DFlowAreas = openHdfGroup( gGeom, "2D Flow Areas" ); HdfGroup gGeom2DFlowAreas = openHdfGroup( gGeom, "2D Flow Areas" );
std::vector<std::string> flowAreaNames = read2DFlowAreasNames( gGeom2DFlowAreas ); std::vector<std::string> flowAreaNames;
if ( oldFormat )
flowAreaNames = read2DFlowAreasNamesOld( gGeom2DFlowAreas );
else
flowAreaNames = read2DFlowAreasNames505( gGeom2DFlowAreas );
std::vector<size_t> areaElemStartIndex( flowAreaNames.size() + 1 ); std::vector<size_t> areaElemStartIndex( flowAreaNames.size() + 1 );
parseMesh( gGeom2DFlowAreas, areaElemStartIndex, flowAreaNames ); parseMesh( gGeom2DFlowAreas, areaElemStartIndex, flowAreaNames );
@ -520,6 +607,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverHec2D::load( const std::string &resultsF
// Face centered Values // Face centered Values
readFaceResults( hdfFile, areaElemStartIndex, flowAreaNames ); readFaceResults( hdfFile, areaElemStartIndex, flowAreaNames );
} }
catch ( MDAL_Status error ) catch ( MDAL_Status error )
{ {

26
external/mdal/frmts/mdal_hec2d.hpp vendored
View File

@ -16,6 +16,17 @@
namespace MDAL namespace MDAL
{ {
/**
* HEC-RAS 2D format.
*
* This is a HDF5-based format to store mesh and datasets
* in a single file. The format supports meshes is multiple
* (disconnected) areas.
*
* There is a small change in the format in HEC-RAS 5.0.5+, where
* - Header File Type is different (HEC-RAS Results vs HEC-RAS Geometry)
* - Names or areas are stored in different place (names array vs attributes array)
*/
class DriverHec2D: public Driver class DriverHec2D: public Driver
{ {
public: public:
@ -30,6 +41,15 @@ namespace MDAL
std::unique_ptr< MDAL::MemoryMesh > mMesh; std::unique_ptr< MDAL::MemoryMesh > mMesh;
std::string mFileName; std::string mFileName;
// Pre 5.0.5 format
bool canReadOldFormat( const std::string &fileType ) const;
std::vector<std::string> read2DFlowAreasNamesOld( HdfGroup gGeom2DFlowAreas ) const;
// 5.0.5 + format
bool canReadFormat505( const std::string &fileType ) const;
std::vector<std::string> read2DFlowAreasNames505( HdfGroup gGeom2DFlowAreas ) const;
// Common functions
void readFaceOutput( const HdfFile &hdfFile, void readFaceOutput( const HdfFile &hdfFile,
const HdfGroup &rootGroup, const HdfGroup &rootGroup,
const std::vector<size_t> &areaElemStartIndex, const std::vector<size_t> &areaElemStartIndex,
@ -37,9 +57,11 @@ namespace MDAL
const std::string rawDatasetName, const std::string rawDatasetName,
const std::string datasetName, const std::string datasetName,
const std::vector<float> &times ); const std::vector<float> &times );
void readFaceResults( const HdfFile &hdfFile, void readFaceResults( const HdfFile &hdfFile,
const std::vector<size_t> &areaElemStartIndex, const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames ); const std::vector<std::string> &flowAreaNames );
std::shared_ptr<MDAL::MemoryDataset> readElemOutput( std::shared_ptr<MDAL::MemoryDataset> readElemOutput(
const HdfGroup &rootGroup, const HdfGroup &rootGroup,
const std::vector<size_t> &areaElemStartIndex, const std::vector<size_t> &areaElemStartIndex,
@ -48,14 +70,18 @@ namespace MDAL
const std::string datasetName, const std::string datasetName,
const std::vector<float> &times, const std::vector<float> &times,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation ); std::shared_ptr<MDAL::MemoryDataset> bed_elevation );
std::shared_ptr<MDAL::MemoryDataset> readBedElevation( std::shared_ptr<MDAL::MemoryDataset> readBedElevation(
const HdfGroup &gGeom2DFlowAreas, const HdfGroup &gGeom2DFlowAreas,
const std::vector<size_t> &areaElemStartIndex, const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames ); const std::vector<std::string> &flowAreaNames );
void setProjection( HdfFile hdfFile ); void setProjection( HdfFile hdfFile );
void parseMesh( HdfGroup gGeom2DFlowAreas, void parseMesh( HdfGroup gGeom2DFlowAreas,
std::vector<size_t> &areaElemStartIndex, std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames ); const std::vector<std::string> &flowAreaNames );
void readElemResults( void readElemResults(
const HdfFile &hdfFile, const HdfFile &hdfFile,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation, std::shared_ptr<MDAL::MemoryDataset> bed_elevation,

206
external/mdal/frmts/mdal_netcdf.cpp vendored Normal file
View File

@ -0,0 +1,206 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#include <string>
#include <vector>
#include <assert.h>
#include <netcdf.h>
#include "mdal_netcdf.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
NetCDFFile::NetCDFFile(): mNcid( 0 ) {}
NetCDFFile::~NetCDFFile()
{
if ( mNcid != 0 ) nc_close( mNcid );
}
int NetCDFFile::handle() const
{
return mNcid;
}
void NetCDFFile::openFile( const std::string &fileName )
{
int res = nc_open( fileName.c_str(), NC_NOWRITE, &mNcid );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_UnknownFormat;
}
}
bool NetCDFFile::hasVariable( const std::string &name ) const
{
assert( mNcid != 0 );
int varid;
return ( nc_inq_varid( mNcid, name.c_str(), &varid ) == NC_NOERR );
}
std::vector<int> NetCDFFile::readIntArr( const std::string &name, size_t dim ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
std::vector<int> arr_val( dim );
if ( nc_get_var_int( mNcid, arr_id, arr_val.data() ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}
std::vector<double> NetCDFFile::readDoubleArr( const std::string &name, size_t dim ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
std::vector<double> arr_val( dim );
if ( nc_get_var_double( mNcid, arr_id, arr_val.data() ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}
bool NetCDFFile::hasArr( const std::string &name ) const
{
assert( mNcid != 0 );
int arr_id;
return nc_inq_varid( mNcid, name.c_str(), &arr_id ) == NC_NOERR;
}
std::vector<std::string> NetCDFFile::readArrNames() const
{
assert( mNcid != 0 );
std::vector<std::string> res;
int nvars;
if ( nc_inq_varids( mNcid, &nvars, nullptr ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
std::vector<int> varids( static_cast<size_t>( nvars ) );
if ( nc_inq_varids( mNcid, &nvars, varids.data() ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
for ( size_t i = 0; i < static_cast<size_t>( nvars ); ++i )
{
std::vector<char> cname( NC_MAX_NAME + 1 );
if ( nc_inq_varname( mNcid, varids[i], cname.data() ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
res.push_back( cname.data() );
}
return res;
}
int NetCDFFile::getAttrInt( const std::string &name, const std::string &attr_name ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
int val;
if ( nc_get_att_int( mNcid, arr_id, attr_name.c_str(), &val ) )
{
throw MDAL_Status::Err_UnknownFormat;
}
return val;
}
std::string NetCDFFile::getAttrStr( const std::string &name, const std::string &attr_name ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) ) throw MDAL_Status::Err_UnknownFormat;
return getAttrStr( attr_name, arr_id );
}
std::string NetCDFFile::getAttrStr( const std::string &name, int varid ) const
{
assert( mNcid != 0 );
size_t attlen = 0;
if ( nc_inq_attlen( mNcid, varid, name.c_str(), &attlen ) )
{
// attribute is missing
return std::string();
}
char *string_attr;
string_attr = ( char * ) malloc( attlen + 1 );
if ( nc_get_att_text( mNcid, varid, name.c_str(), string_attr ) ) throw MDAL_Status::Err_UnknownFormat;
string_attr[attlen] = '\0';
std::string res( string_attr );
free( string_attr );
return res;
}
double NetCDFFile::getFillValue( int varid ) const
{
return getAttrDouble( varid, "_FillValue" );
}
double NetCDFFile::getAttrDouble( int varid, const std::string &attr_name ) const
{
double res;
if ( nc_get_att_double( mNcid, varid, attr_name.c_str(), &res ) )
res = std::numeric_limits<double>::quiet_NaN(); // not present/set
return res;
}
int NetCDFFile::getVarId( const std::string &name )
{
int ncid_val;
if ( nc_inq_varid( mNcid, name.c_str(), &ncid_val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
return ncid_val;
}
void NetCDFFile::getDimension( const std::string &name, size_t *val, int *ncid_val ) const
{
assert( mNcid != 0 );
if ( nc_inq_dimid( mNcid, name.c_str(), ncid_val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_dimlen( mNcid, *ncid_val, val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
void NetCDFFile::getDimensionOptional( const std::string &name, size_t *val, int *ncid_val ) const
{
assert( mNcid != 0 );
try
{
getDimension( name, val, ncid_val );
}
catch ( MDAL_Status )
{
*ncid_val = -1;
*val = 0;
}
}

175
external/mdal/frmts/mdal_netcdf.hpp vendored
View File

@ -8,168 +8,35 @@
#include <string> #include <string>
#include <vector> #include <vector>
#include <assert.h>
#include <netcdf.h>
#include "mdal_utils.hpp"
#include "mdal.h"
//! C++ Wrapper around netcdf C library
class NetCDFFile class NetCDFFile
{ {
public: public:
NetCDFFile(): mNcid( 0 ) {} //! Create file with invalid handle
~NetCDFFile() NetCDFFile();
{ //! Closes the file
if ( mNcid != 0 ) nc_close( mNcid ); ~NetCDFFile();
}
int handle() const int handle() const;
{ void openFile( const std::string &fileName );
return mNcid; bool hasVariable( const std::string &name ) const;
}
void openFile( const std::string &fileName ) std::vector<int> readIntArr( const std::string &name, size_t dim ) const;
{ std::vector<double> readDoubleArr( const std::string &name, size_t dim ) const;
int res = nc_open( fileName.c_str(), NC_NOWRITE, &mNcid ); bool hasArr( const std::string &name ) const;
if ( res != NC_NOERR ) std::vector<std::string> readArrNames() const;
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_UnknownFormat;
}
}
bool hasVariable( const std::string &name ) const int getAttrInt( const std::string &name, const std::string &attr_name ) const;
{ double getAttrDouble( int varid, const std::string &attr_name ) const;
assert( mNcid != 0 ); std::string getAttrStr( const std::string &name, const std::string &attr_name ) const;
int varid; std::string getAttrStr( const std::string &name, int varid ) const;
return ( nc_inq_varid( mNcid, name.c_str(), &varid ) == NC_NOERR ); double getFillValue( int varid ) const;
} int getVarId( const std::string &name );
void getDimension( const std::string &name, size_t *val, int *ncid_val ) const;
std::vector<int> readIntArr( const std::string &name, size_t dim ) const void getDimensionOptional( const std::string &name, size_t *val, int *ncid_val ) const;
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
std::vector<int> arr_val( dim );
if ( nc_get_var_int( mNcid, arr_id, arr_val.data() ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}
std::vector<double> readDoubleArr( const std::string &name, size_t dim ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
std::vector<double> arr_val( dim );
if ( nc_get_var_double( mNcid, arr_id, arr_val.data() ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}
int getAttrInt( const std::string &name, const std::string &attr_name ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
throw MDAL_Status::Err_UnknownFormat;
}
int val;
if ( nc_get_att_int( mNcid, arr_id, attr_name.c_str(), &val ) )
{
throw MDAL_Status::Err_UnknownFormat;
}
return val;
}
std::string getAttrStr( const std::string &name, const std::string &attr_name ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) ) throw MDAL_Status::Err_UnknownFormat;
return getAttrStr( attr_name, arr_id );
}
std::string getAttrStr( const std::string &name, int varid ) const
{
assert( mNcid != 0 );
size_t attlen = 0;
if ( nc_inq_attlen( mNcid, varid, name.c_str(), &attlen ) )
{
// attribute is missing
return std::string();
}
char *string_attr;
string_attr = ( char * ) malloc( attlen + 1 );
if ( nc_get_att_text( mNcid, varid, name.c_str(), string_attr ) ) throw MDAL_Status::Err_UnknownFormat;
string_attr[attlen] = '\0';
std::string res( string_attr );
free( string_attr );
return res;
}
double getFillValue( int varid ) const
{
double fill;
if ( nc_get_att_double( mNcid, varid, "_FillValue", &fill ) )
fill = std::numeric_limits<double>::quiet_NaN(); // not present/set
return fill;
}
int getVarId( const std::string &name )
{
int ncid_val;
if ( nc_inq_varid( mNcid, name.c_str(), &ncid_val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
return ncid_val;
}
void getDimension( const std::string &name, size_t *val, int *ncid_val ) const
{
assert( mNcid != 0 );
if ( nc_inq_dimid( mNcid, name.c_str(), ncid_val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_dimlen( mNcid, *ncid_val, val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
void getDimensionOptional( const std::string &name, size_t *val, int *ncid_val ) const
{
assert( mNcid != 0 );
try
{
getDimension( name, val, ncid_val );
}
catch ( MDAL_Status )
{
*ncid_val = -1;
*val = 0;
}
}
private: private:
int mNcid; int mNcid; // C handle to the file
}; };
#endif // MDAL_NETCDF_HPP #endif // MDAL_NETCDF_HPP

761
external/mdal/frmts/mdal_sww.cpp vendored
View File

@ -7,14 +7,11 @@
#include <string> #include <string>
#include <string.h> #include <string.h>
#include <vector> #include <vector>
#include <netcdf.h>
#include <set>
#include "mdal_netcdf.hpp"
#include "mdal_sww.hpp" #include "mdal_sww.hpp"
#include "mdal_utils.hpp"
// threshold for determining whether an element is active (wet)
// the format does not explicitly store that information so we
// determine that when loading data
#define DEPTH_THRESHOLD 0.0001 // in meters
MDAL::DriverSWW::DriverSWW() MDAL::DriverSWW::DriverSWW()
: Driver( "SWW", : Driver( "SWW",
@ -29,350 +26,442 @@ MDAL::DriverSWW *MDAL::DriverSWW::create()
return new DriverSWW(); return new DriverSWW();
} }
bool MDAL::DriverSWW::canRead( const std::string &uri ) size_t MDAL::DriverSWW::getVertexCount( const NetCDFFile &ncFile ) const
{ {
int ncid; int nPointsId;
int res; size_t res;
ncFile.getDimension( "number_of_points", &res, &nPointsId );
// open return res;
res = nc_open( uri.c_str(), NC_NOWRITE, &ncid );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
nc_close( ncid );
return false;
}
// get dimensions
int nVolumesId, nVerticesId, nPointsId, nTimestepsId;
if ( nc_inq_dimid( ncid, "number_of_volumes", &nVolumesId ) != NC_NOERR ||
nc_inq_dimid( ncid, "number_of_vertices", &nVerticesId ) != NC_NOERR ||
nc_inq_dimid( ncid, "number_of_points", &nPointsId ) != NC_NOERR ||
nc_inq_dimid( ncid, "number_of_timesteps", &nTimestepsId ) != NC_NOERR )
{
nc_close( ncid );
return false;
}
return true;
} }
std::unique_ptr<MDAL::Mesh> MDAL::DriverSWW::load( const std::string &resultsFile, bool MDAL::DriverSWW::canRead( const std::string &uri )
MDAL_Status *status ) {
NetCDFFile ncFile;
try
{
ncFile.openFile( uri );
getVertexCount( ncFile );
}
catch ( MDAL_Status )
{
return false;
}
return true;
}
std::vector<double> MDAL::DriverSWW::readZCoords( const NetCDFFile &ncFile ) const
{
size_t nPoints = getVertexCount( ncFile );
std::vector<double> pz;
// newer sww files does have elevation array that is time-dependent
if ( ncFile.hasArr( "z" ) )
{
pz = ncFile.readDoubleArr( "z", nPoints );
}
else if ( ncFile.hasArr( "elevation" ) )
{
int zDims = 0;
int zid;
if ( nc_inq_varid( ncFile.handle(), "elevation", &zid ) == NC_NOERR &&
nc_inq_varndims( ncFile.handle(), zid, &zDims ) == NC_NOERR )
{
if ( zDims == 1 )
{
// just one elevation for all times, treat as z coord
pz = ncFile.readDoubleArr( "elevation", nPoints );
}
else
{
pz.resize( nPoints );
// fetching "elevation" data for the first timestep,
// and threat it as z coord
size_t start[2], count[2];
const ptrdiff_t stride[2] = {1, 1};
start[0] = 0; // t = 0
start[1] = 0;
count[0] = 1;
count[1] = nPoints;
nc_get_vars_double( ncFile.handle(), zid, start, count, stride, pz.data() );
}
}
}
return pz;
}
void MDAL::DriverSWW::addBedElevation( const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times
) const
{
if ( ncFile.hasArr( "elevation" ) )
{
std::shared_ptr<MDAL::DatasetGroup> grp = readScalarGroup( ncFile,
mesh,
times,
"Bed Elevation",
"elevation" );
mesh->datasetGroups.push_back( grp );
}
else
{
MDAL::addBedElevationDatasetGroup( mesh, mesh->vertices );
}
}
MDAL::Vertices MDAL::DriverSWW::readVertices( const NetCDFFile &ncFile ) const
{
size_t nPoints = getVertexCount( ncFile );
// load mesh data
std::vector<double> px = ncFile.readDoubleArr( "x", nPoints );
std::vector<double> py = ncFile.readDoubleArr( "y", nPoints );
std::vector<double> pz = readZCoords( ncFile );
// we may need to apply a shift to the X,Y coordinates
double xLLcorner = ncFile.getAttrDouble( NC_GLOBAL, "xllcorner" );
double yLLcorner = ncFile.getAttrDouble( NC_GLOBAL, "yllcorner" );
MDAL::Vertices vertices( nPoints );
Vertex *vertexPtr = vertices.data();
for ( size_t i = 0; i < nPoints; ++i, ++vertexPtr )
{
vertexPtr->x = px[i] + xLLcorner ;
vertexPtr->y = py[i] + yLLcorner ;
if ( !pz.empty() ) // missing both "z" and "elevation"
vertexPtr->z = pz[i];
}
return vertices;
}
MDAL::Faces MDAL::DriverSWW::readFaces( const NetCDFFile &ncFile ) const
{
// get dimensions
int nVolumesId, nVerticesId;
size_t nVolumes, nVertices;
ncFile.getDimension( "number_of_volumes", &nVolumes, &nVolumesId );
ncFile.getDimension( "number_of_vertices", &nVertices, &nVerticesId );
if ( nVertices != 3 )
throw MDAL_Status::Err_UnknownFormat;
std::vector<int> pvolumes = ncFile.readIntArr( "volumes", nVertices * nVolumes );
MDAL::Faces faces( nVolumes );
for ( size_t i = 0; i < nVolumes; ++i )
{
faces[i].resize( 3 );
faces[i][0] = static_cast<size_t>( pvolumes[3 * i + 0] );
faces[i][1] = static_cast<size_t>( pvolumes[3 * i + 1] );
faces[i][2] = static_cast<size_t>( pvolumes[3 * i + 2] );
}
return faces;
}
std::vector<double> MDAL::DriverSWW::readTimes( const NetCDFFile &ncFile ) const
{
size_t nTimesteps;
int nTimestepsId;
ncFile.getDimension( "number_of_timesteps", &nTimesteps, &nTimestepsId );
std::vector<double> times = ncFile.readDoubleArr( "time", nTimesteps );
return times;
}
void MDAL::DriverSWW::readDatasetGroups(
const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times
) const
{
std::set<std::string> parsedVariableNames; // already parsed arrays somewhere else
parsedVariableNames.insert( "x" );
parsedVariableNames.insert( "y" );
parsedVariableNames.insert( "z" );
parsedVariableNames.insert( "volumes" );
parsedVariableNames.insert( "time" );
std::vector<std::string> names = ncFile.readArrNames();
std::set<std::string> namesSet( names.begin(), names.end() );
// Add bed elevation group
parsedVariableNames.insert( "elevations" );
addBedElevation( ncFile, mesh, times );
for ( const std::string &name : names )
{
// skip already parsed variables
if ( parsedVariableNames.find( name ) == parsedVariableNames.cend() )
{
std::string xName, yName, groupName( name );
bool isVector = parseGroupName( groupName, xName, yName );
std::shared_ptr<MDAL::DatasetGroup> grp;
if ( isVector && ncFile.hasArr( xName ) && ncFile.hasArr( yName ) )
{
// vector dataset group
grp = readVectorGroup(
ncFile,
mesh,
times,
groupName,
xName,
yName );
parsedVariableNames.insert( xName );
parsedVariableNames.insert( yName );
}
else
{
// scalar dataset group
grp = readScalarGroup(
ncFile,
mesh,
times,
groupName,
name );
parsedVariableNames.insert( name );
}
if ( grp )
mesh->datasetGroups.push_back( grp );
}
}
}
bool MDAL::DriverSWW::parseGroupName( std::string &groupName,
std::string &xName,
std::string &yName ) const
{
bool isVector = false;
std::string baseName( groupName );
// X and Y variables
if ( groupName.size() > 1 )
{
if ( MDAL::startsWith( groupName, "x" ) )
{
baseName = groupName.substr( 1, groupName.size() - 1 );
xName = groupName;
yName = "y" + baseName;
isVector = true;
}
else if ( MDAL::startsWith( groupName, "y" ) )
{
baseName = groupName.substr( 1, groupName.size() - 1 );
xName = "x" + baseName;
yName = groupName;
isVector = true;
}
}
// Maximums
groupName = baseName;
if ( MDAL::endsWith( baseName, "_range" ) )
{
groupName = groupName.substr( 0, baseName.size() - 6 ) + "/Maximums";
}
return isVector;
}
std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readScalarGroup(
const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times,
const std::string groupName,
const std::string arrName
) const
{
size_t nPoints = getVertexCount( ncFile );
std::shared_ptr<MDAL::DatasetGroup> mds;
int varxid;
if ( nc_inq_varid( ncFile.handle(), arrName.c_str(), &varxid ) == NC_NOERR )
{
mds = std::make_shared<MDAL::DatasetGroup> (
name(),
mesh,
mFileName,
groupName );
mds->setIsOnVertices( true );
mds->setIsScalar( true );
int zDimsX = 0;
if ( nc_inq_varndims( ncFile.handle(), varxid, &zDimsX ) != NC_NOERR )
throw MDAL_Status::Err_UnknownFormat;
if ( zDimsX == 1 )
{
// TIME INDEPENDENT
std::shared_ptr<MDAL::MemoryDataset> o = std::make_shared<MDAL::MemoryDataset>( mds.get() );
o->setTime( 0.0 );
double *values = o->values();
std::vector<double> valuesX = ncFile.readDoubleArr( arrName, nPoints );
for ( size_t i = 0; i < nPoints; ++i )
{
values[i] = valuesX[i];
}
o->setStatistics( MDAL::calculateStatistics( o ) );
mds->datasets.push_back( o );
}
else
{
// TIME DEPENDENT
for ( size_t t = 0; t < times.size(); ++t )
{
std::shared_ptr<MDAL::MemoryDataset> mto = std::make_shared<MDAL::MemoryDataset>( mds.get() );
mto->setTime( static_cast<double>( times[t] ) / 3600. );
double *values = mto->values();
// fetching data for one timestep
size_t start[2], count[2];
const ptrdiff_t stride[2] = {1, 1};
start[0] = t;
start[1] = 0;
count[0] = 1;
count[1] = nPoints;
nc_get_vars_double( ncFile.handle(), varxid, start, count, stride, values );
mto->setStatistics( MDAL::calculateStatistics( mto ) );
mds->datasets.push_back( mto );
}
}
mds->setStatistics( MDAL::calculateStatistics( mds ) );
}
return mds;
}
std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readVectorGroup(
const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times,
const std::string groupName,
const std::string arrXName,
const std::string arrYName
) const
{
size_t nPoints = getVertexCount( ncFile );
std::shared_ptr<MDAL::DatasetGroup> mds;
int varxid, varyid;
if ( nc_inq_varid( ncFile.handle(), arrXName.c_str(), &varxid ) == NC_NOERR &&
nc_inq_varid( ncFile.handle(), arrYName.c_str(), &varyid ) == NC_NOERR )
{
mds = std::make_shared<MDAL::DatasetGroup> (
name(),
mesh,
mFileName,
groupName );
mds->setIsOnVertices( true );
mds->setIsScalar( false );
int zDimsX = 0;
int zDimsY = 0;
if ( nc_inq_varndims( ncFile.handle(), varxid, &zDimsX ) != NC_NOERR )
throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_varndims( ncFile.handle(), varyid, &zDimsY ) != NC_NOERR )
throw MDAL_Status::Err_UnknownFormat;
if ( zDimsX != zDimsY )
throw MDAL_Status::Err_UnknownFormat;
std::vector<double> valuesX( nPoints ), valuesY( nPoints );
if ( zDimsX == 1 )
{
// TIME INDEPENDENT
std::shared_ptr<MDAL::MemoryDataset> o = std::make_shared<MDAL::MemoryDataset>( mds.get() );
o->setTime( 0.0 );
double *values = o->values();
std::vector<double> valuesX = ncFile.readDoubleArr( arrXName, nPoints );
std::vector<double> valuesY = ncFile.readDoubleArr( arrYName, nPoints );
for ( size_t i = 0; i < nPoints; ++i )
{
values[2 * i] = valuesX[i];
values[2 * i + 1] = valuesY[i];
}
o->setStatistics( MDAL::calculateStatistics( o ) );
mds->datasets.push_back( o );
}
else
{
// TIME DEPENDENT
for ( size_t t = 0; t < times.size(); ++t )
{
std::shared_ptr<MDAL::MemoryDataset> mto = std::make_shared<MDAL::MemoryDataset>( mds.get() );
mto->setTime( static_cast<double>( times[t] ) / 3600. );
double *values = mto->values();
// fetching data for one timestep
size_t start[2], count[2];
const ptrdiff_t stride[2] = {1, 1};
start[0] = t;
start[1] = 0;
count[0] = 1;
count[1] = nPoints;
nc_get_vars_double( ncFile.handle(), varxid, start, count, stride, valuesX.data() );
nc_get_vars_double( ncFile.handle(), varyid, start, count, stride, valuesY.data() );
for ( size_t i = 0; i < nPoints; ++i )
{
values[2 * i] = static_cast<double>( valuesX[i] );
values[2 * i + 1] = static_cast<double>( valuesY[i] );
}
mto->setStatistics( MDAL::calculateStatistics( mto ) );
mds->datasets.push_back( mto );
}
}
mds->setStatistics( MDAL::calculateStatistics( mds ) );
}
return mds;
}
std::unique_ptr<MDAL::Mesh> MDAL::DriverSWW::load(
const std::string &resultsFile,
MDAL_Status *status )
{ {
mFileName = resultsFile; mFileName = resultsFile;
if ( status ) *status = MDAL_Status::None; if ( status ) *status = MDAL_Status::None;
int ncid; NetCDFFile ncFile;
int res;
res = nc_open( mFileName.c_str(), NC_NOWRITE, &ncid ); try
if ( res != NC_NOERR )
{ {
MDAL::debug( nc_strerror( res ) ); // Open file for reading
nc_close( ncid ); ncFile.openFile( mFileName );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
// get dimensions // Read mesh
int nVolumesId, nVerticesId, nPointsId, nTimestepsId; MDAL::Vertices vertices = readVertices( ncFile );
size_t nVolumes, nVertices, nPoints, nTimesteps; MDAL::Faces faces = readFaces( ncFile );
if ( nc_inq_dimid( ncid, "number_of_volumes", &nVolumesId ) != NC_NOERR || std::unique_ptr< MDAL::MemoryMesh > mesh(
nc_inq_dimid( ncid, "number_of_vertices", &nVerticesId ) != NC_NOERR || new MemoryMesh(
nc_inq_dimid( ncid, "number_of_points", &nPointsId ) != NC_NOERR ||
nc_inq_dimid( ncid, "number_of_timesteps", &nTimestepsId ) != NC_NOERR )
{
nc_close( ncid );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
if ( nc_inq_dimlen( ncid, nVolumesId, &nVolumes ) != NC_NOERR ||
nc_inq_dimlen( ncid, nVerticesId, &nVertices ) != NC_NOERR ||
nc_inq_dimlen( ncid, nPointsId, &nPoints ) != NC_NOERR ||
nc_inq_dimlen( ncid, nTimestepsId, &nTimesteps ) != NC_NOERR )
{
nc_close( ncid );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
if ( nVertices != 3 )
{
MDAL::debug( "Expecting triangular elements!" );
nc_close( ncid );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
int xid, yid, zid, volumesid, timeid, stageid;
if ( nc_inq_varid( ncid, "x", &xid ) != NC_NOERR ||
nc_inq_varid( ncid, "y", &yid ) != NC_NOERR ||
nc_inq_varid( ncid, "volumes", &volumesid ) != NC_NOERR ||
nc_inq_varid( ncid, "time", &timeid ) != NC_NOERR ||
nc_inq_varid( ncid, "stage", &stageid ) != NC_NOERR )
{
nc_close( ncid );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
// load mesh data
std::vector<float> px( nPoints ), py( nPoints ), pz( nPoints );
std::vector<int> pvolumes( nVertices * nVolumes );
if ( nc_get_var_float( ncid, xid, px.data() ) != NC_NOERR ||
nc_get_var_float( ncid, yid, py.data() ) != NC_NOERR ||
nc_get_var_int( ncid, volumesid, pvolumes.data() ) != NC_NOERR )
{
nc_close( ncid );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
// we may need to apply a shift to the X,Y coordinates
float xLLcorner = 0, yLLcorner = 0;
nc_get_att_float( ncid, NC_GLOBAL, "xllcorner", &xLLcorner );
nc_get_att_float( ncid, NC_GLOBAL, "yllcorner", &yLLcorner );
MDAL::Vertices nodes( nPoints );
Vertex *nodesPtr = nodes.data();
for ( size_t i = 0; i < nPoints; ++i, ++nodesPtr )
{
nodesPtr->x = static_cast<double>( px[i] + xLLcorner );
nodesPtr->y = static_cast<double>( py[i] + yLLcorner );
}
std::vector<float> times( nTimesteps );
nc_get_var_float( ncid, timeid, times.data() );
int zDims = 0;
if ( nc_inq_varid( ncid, "z", &zid ) == NC_NOERR &&
nc_get_var_float( ncid, zid, pz.data() ) == NC_NOERR )
{
// older SWW format: elevation is constant over time
zDims = 1;
}
else if ( nc_inq_varid( ncid, "elevation", &zid ) == NC_NOERR &&
nc_inq_varndims( ncid, zid, &zDims ) == NC_NOERR &&
( ( zDims == 1 && nc_get_var_float( ncid, zid, pz.data() ) == NC_NOERR ) || zDims == 2 ) )
{
// we're good
}
else
{
// neither "z" nor "elevation" are present -> something is going wrong
nc_close( ncid );
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
return std::unique_ptr< MDAL::Mesh >();
}
MDAL::Faces elements( nVolumes );
for ( size_t i = 0; i < nVolumes; ++i )
{
elements[i].resize( 3 );
elements[i][0] = static_cast<size_t>( pvolumes[3 * i + 0] );
elements[i][1] = static_cast<size_t>( pvolumes[3 * i + 1] );
elements[i][2] = static_cast<size_t>( pvolumes[3 * i + 2] );
}
std::unique_ptr< MDAL::MemoryMesh > mesh(
new MemoryMesh(
name(),
nodes.size(),
elements.size(),
3, // triangles
computeExtent( nodes ),
mFileName
)
);
mesh->faces = elements;
mesh->vertices = nodes;
// Create a dataset for the bed elevation
std::shared_ptr<MDAL::DatasetGroup> bedDs = std::make_shared<MDAL::DatasetGroup> (
name(), name(),
mesh.get(), vertices.size(),
mFileName, faces.size(),
"Bed Elevation" ); 3, // triangles
bedDs->setIsOnVertices( true ); computeExtent( vertices ),
bedDs->setIsScalar( true ); mFileName
)
);
mesh->faces = faces;
mesh->vertices = vertices;
// read bed elevations // Read times
std::vector<std::shared_ptr<MDAL::MemoryDataset>> elevationOutputs; std::vector<double> times = readTimes( ncFile );
if ( zDims == 1 )
{ // Create a dataset(s)
// either "z" or "elevation" with 1 dimension readDatasetGroups( ncFile, mesh.get(), times );
std::shared_ptr<MDAL::MemoryDataset> o = std::make_shared<MDAL::MemoryDataset>( bedDs.get() );
o->setTime( 0.0 ); // Success!
double *values = o->values(); return std::unique_ptr<Mesh>( mesh.release() );
for ( size_t i = 0; i < nPoints; ++i )
{
double z = static_cast<double>( pz[i] );
values[i] = z;
mesh->vertices[i].z = z;
}
o->setStatistics( MDAL::calculateStatistics( o ) );
bedDs->datasets.push_back( o );
elevationOutputs.push_back( o );
} }
else if ( zDims == 2 ) catch ( MDAL_Status err )
{ {
// newer SWW format: elevation may change over time if ( status ) *status = err;
for ( size_t t = 0; t < nTimesteps; ++t ) return std::unique_ptr< MDAL::Mesh >();
{
std::shared_ptr<MDAL::MemoryDataset> toe = std::make_shared<MDAL::MemoryDataset>( bedDs.get() );
toe->setTime( static_cast<double>( times[t] ) / 3600. );
double *elev = toe->values();
// fetching "elevation" data for one timestep
size_t start[2], count[2];
const ptrdiff_t stride[2] = {1, 1};
start[0] = t;
start[1] = 0;
count[0] = 1;
count[1] = nPoints;
std::vector<float> buffer( nPoints );
nc_get_vars_float( ncid, zid, start, count, stride, buffer.data() );
for ( size_t i = 0; i < nPoints; ++i )
{
double val = static_cast<double>( buffer[i] );
elev[i] = val;
}
toe->setStatistics( MDAL::calculateStatistics( toe ) );
bedDs->datasets.push_back( toe );
elevationOutputs.push_back( toe );
}
} }
bedDs->setStatistics( MDAL::calculateStatistics( bedDs ) );
mesh->datasetGroups.push_back( bedDs );
// load results
std::shared_ptr<MDAL::DatasetGroup> dss = std::make_shared<MDAL::DatasetGroup> (
name(),
mesh.get(),
mFileName,
"Stage" );
dss->setIsOnVertices( true );
dss->setIsScalar( true );
std::shared_ptr<MDAL::DatasetGroup> dsd = std::make_shared<MDAL::DatasetGroup> (
name(),
mesh.get(),
mFileName,
"Depth" );
dsd->setIsOnVertices( true );
dsd->setIsScalar( true );
for ( size_t t = 0; t < nTimesteps; ++t )
{
const std::shared_ptr<MDAL::MemoryDataset> elevO = elevationOutputs.size() > 1 ? elevationOutputs[t] : elevationOutputs[0];
const double *elev = elevO->constValues();
std::shared_ptr<MDAL::MemoryDataset> tos = std::make_shared<MDAL::MemoryDataset>( dss.get() );
tos->setTime( static_cast<double>( times[t] ) / 3600. );
double *values = tos->values();
// fetching "stage" data for one timestep
size_t start[2], count[2];
const ptrdiff_t stride[2] = {1, 1};
start[0] = t;
start[1] = 0;
count[0] = 1;
count[1] = nPoints;
std::vector<float> buffer( nPoints );
nc_get_vars_float( ncid, stageid, start, count, stride, buffer.data() );
for ( size_t i = 0; i < nPoints; ++i )
{
double val = static_cast<double>( buffer[i] );
values[i] = val;
}
// derived data: depth = stage - elevation
std::shared_ptr<MDAL::MemoryDataset> tod = std::make_shared<MDAL::MemoryDataset>( dsd.get() );
tod->setTime( tos->time() );
double *depths = tod->values();
int *activeTos = tos->active();
int *activeTod = tod->active();
for ( size_t j = 0; j < nPoints; ++j )
depths[j] = values[j] - elev[j];
// determine which elements are active (wet)
for ( size_t elemidx = 0; elemidx < nVolumes; ++elemidx )
{
const Face &elem = mesh->faces[elemidx];
double v0 = depths[elem[0]];
double v1 = depths[elem[1]];
double v2 = depths[elem[2]];
activeTos[elemidx] = v0 > DEPTH_THRESHOLD && v1 > DEPTH_THRESHOLD && v2 > DEPTH_THRESHOLD;
activeTod[elemidx] = activeTos[elemidx];
}
tos->setStatistics( MDAL::calculateStatistics( tos ) );
dss->datasets.push_back( tos );
tod->setStatistics( MDAL::calculateStatistics( tod ) );
dsd->datasets.push_back( tod );
}
dss->setStatistics( MDAL::calculateStatistics( dss ) );
mesh->datasetGroups.push_back( dss );
dsd->setStatistics( MDAL::calculateStatistics( dsd ) );
mesh->datasetGroups.push_back( dsd );
int momentumxid, momentumyid;
if ( nc_inq_varid( ncid, "xmomentum", &momentumxid ) == NC_NOERR &&
nc_inq_varid( ncid, "ymomentum", &momentumyid ) == NC_NOERR )
{
std::shared_ptr<MDAL::DatasetGroup> mds = std::make_shared<MDAL::DatasetGroup> (
name(),
mesh.get(),
mFileName,
"Momentum" );
mds->setIsOnVertices( true );
mds->setIsScalar( false );
std::vector<float> valuesX( nPoints ), valuesY( nPoints );
for ( size_t t = 0; t < nTimesteps; ++t )
{
std::shared_ptr<MDAL::MemoryDataset> mto = std::make_shared<MDAL::MemoryDataset>( mds.get() );
mto->setTime( static_cast<double>( times[t] ) / 3600. );
double *values = mto->values();
std::shared_ptr<MDAL::MemoryDataset> mto0 = std::static_pointer_cast<MDAL::MemoryDataset>( dsd->datasets[t] );
memcpy( mto->active(), mto0->active(), mesh->facesCount() * sizeof( int ) );
// fetching "stage" data for one timestep
size_t start[2], count[2];
const ptrdiff_t stride[2] = {1, 1};
start[0] = t;
start[1] = 0;
count[0] = 1;
count[1] = nPoints;
nc_get_vars_float( ncid, momentumxid, start, count, stride, valuesX.data() );
nc_get_vars_float( ncid, momentumyid, start, count, stride, valuesY.data() );
for ( size_t i = 0; i < nPoints; ++i )
{
values[2 * i] = static_cast<double>( valuesX[i] );
values[2 * i + 1] = static_cast<double>( valuesY[i] );
}
mto->setStatistics( MDAL::calculateStatistics( mto ) );
mds->datasets.push_back( mto );
}
mds->setStatistics( MDAL::calculateStatistics( mds ) );
mesh->datasetGroups.push_back( mds );
}
nc_close( ncid );
return std::unique_ptr<Mesh>( mesh.release() );
} }

55
external/mdal/frmts/mdal_sww.hpp vendored
View File

@ -7,15 +7,32 @@
#define MDAL_SWW_HPP #define MDAL_SWW_HPP
#include <string> #include <string>
#include <memory>
#include <vector>
#include <map>
#include <utility>
#include "mdal_data_model.hpp" #include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp" #include "mdal_memory_data_model.hpp"
#include "mdal.h" #include "mdal.h"
#include "mdal_driver.hpp" #include "mdal_driver.hpp"
#include "mdal_netcdf.hpp"
namespace MDAL namespace MDAL
{ {
// AnuGA format with extension .SWW /**
* AnuGA format with extension .SWW
*
* The format is based on NetCDF storage
*
* Bed Elevation can be static (== one bed elevation for all timesteps, stored
* in "z" variable or "elevation" variable)
* or dynamic (each timestep has its own elevation data, stored in "elevation"
* variable with multiple dimensions)
*
* Vector data are recognized by prefix "x" and "y" in the name
* Maximums data are recognized by suffix "_range" in the name
*/
class DriverSWW: public Driver class DriverSWW: public Driver
{ {
public: public:
@ -25,7 +42,43 @@ namespace MDAL
std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override; std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override;
bool canRead( const std::string &uri ) override; bool canRead( const std::string &uri ) override;
private: private:
size_t getVertexCount( const NetCDFFile &ncFile ) const;
std::vector<double> readZCoords( const NetCDFFile &ncFile ) const;
MDAL::Vertices readVertices( const NetCDFFile &ncFile ) const;
MDAL::Faces readFaces( const NetCDFFile &ncFile ) const;
std::vector<double> readTimes( const NetCDFFile &ncFile ) const;
/**
* Finds all variables (arrays) in netcdf file and base on the name add it as
* vector or scalar dataset group
*/
void readDatasetGroups( const NetCDFFile &ncFile, MDAL::MemoryMesh *mesh, const std::vector<double> &times ) const;
bool parseGroupName( std::string &groupName, std::string &xName, std::string &yName ) const;
std::shared_ptr<MDAL::DatasetGroup> readScalarGroup(
const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times,
const std::string variableBaseName,
const std::string arrName
) const;
std::shared_ptr<MDAL::DatasetGroup> readVectorGroup(
const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times,
const std::string variableBaseName,
const std::string arrXName,
const std::string arrYName
) const;
void addBedElevation(
const NetCDFFile &ncFile,
MDAL::MemoryMesh *mesh,
const std::vector<double> &times
) const;
std::string mFileName; std::string mFileName;
}; };
} // namespace MDAL } // namespace MDAL

52
external/mdal/frmts/mdal_xmdf.cpp vendored
View File

@ -139,6 +139,7 @@ void MDAL::DriverXmdf::load( const std::string &datFile, MDAL::Mesh *mesh, MDAL
size_t vertexCount = mesh->verticesCount(); size_t vertexCount = mesh->verticesCount();
size_t faceCount = mesh->facesCount(); size_t faceCount = mesh->facesCount();
std::vector<std::string> rootGroups = file.groups(); std::vector<std::string> rootGroups = file.groups();
if ( rootGroups.size() != 1 ) if ( rootGroups.size() != 1 )
{ {
@ -147,44 +148,21 @@ void MDAL::DriverXmdf::load( const std::string &datFile, MDAL::Mesh *mesh, MDAL
return; return;
} }
HdfGroup gMesh = file.group( rootGroups[0] ); HdfGroup gMesh = file.group( rootGroups[0] );
// TODO: read Times group (e.g. time of peak velocity)
DatasetGroups groups; // DAT outputs data DatasetGroups groups; // DAT outputs data
if ( gMesh.pathExists( "Temporal" ) ) for ( const std::string &groupName : gMesh.groups() )
{ {
HdfGroup gTemporal = gMesh.group( "Temporal" ); HdfGroup gGroup = gMesh.group( groupName );
if ( gTemporal.isValid() ) if ( gGroup.isValid() )
{ {
addDatasetGroupsFromXmdfGroup( groups, gTemporal, vertexCount, faceCount ); if ( groupName == "Maximums" )
}
}
if ( gMesh.pathExists( "Temporal" ) )
{
HdfGroup gMaximums = gMesh.group( "Maximums" );
if ( gMaximums.isValid() )
{
for ( const std::string &groupName : gMaximums.groups() )
{ {
HdfGroup g = gMaximums.group( groupName ); addDatasetGroupsFromXmdfGroup( groups, gGroup, "/Maximums", vertexCount, faceCount );
std::shared_ptr<MDAL::DatasetGroup> maxGroup = readXmdfGroupAsDatasetGroup( g, groupName + "/Maximums", vertexCount, faceCount ); }
if ( !maxGroup || maxGroup->datasets.size() != 1 ) else
MDAL::debug( "Maximum dataset should have just one timestep!" ); {
else addDatasetGroupsFromXmdfGroup( groups, gGroup, "", vertexCount, faceCount );
groups.push_back( maxGroup );
} }
}
}
// res_to_res.exe (TUFLOW utiity tool)
if ( gMesh.pathExists( "Difference" ) )
{
HdfGroup gDifference = gMesh.group( "Difference" );
if ( gDifference.isValid() )
{
addDatasetGroupsFromXmdfGroup( groups, gDifference, vertexCount, faceCount );
} }
} }
@ -195,14 +173,20 @@ void MDAL::DriverXmdf::load( const std::string &datFile, MDAL::Mesh *mesh, MDAL
); );
} }
void MDAL::DriverXmdf::addDatasetGroupsFromXmdfGroup( DatasetGroups &groups, const HdfGroup &rootGroup, size_t vertexCount, size_t faceCount ) void MDAL::DriverXmdf::addDatasetGroupsFromXmdfGroup( DatasetGroups &groups,
const HdfGroup &rootGroup,
const std::string &nameSuffix,
size_t vertexCount,
size_t faceCount )
{ {
for ( const std::string &groupName : rootGroup.groups() ) for ( const std::string &groupName : rootGroup.groups() )
{ {
HdfGroup g = rootGroup.group( groupName ); HdfGroup g = rootGroup.group( groupName );
std::shared_ptr<DatasetGroup> ds = readXmdfGroupAsDatasetGroup( g, groupName, vertexCount, faceCount ); std::shared_ptr<DatasetGroup> ds = readXmdfGroupAsDatasetGroup( g, groupName + nameSuffix, vertexCount, faceCount );
if ( ds && ds->datasets.size() > 0 ) if ( ds && ds->datasets.size() > 0 )
{
groups.push_back( ds ); groups.push_back( ds );
}
} }
} }

21
external/mdal/frmts/mdal_xmdf.hpp vendored
View File

@ -58,6 +58,26 @@ namespace MDAL
class DriverXmdf: public Driver class DriverXmdf: public Driver
{ {
public: public:
/**
* Driver for XMDF Files
*
* Structure of the TUFLOW file. Groups are optional since it depends
* on tools which groups are created.
* - root
* - Temporal
* - Depth
* - Velocity
* - ..
* - Maximums
* - Depth
* - Velocity
* - ..
* - Difference (res_to_res.exe TUFLOW utility tool)
* - ..
* - Times (e.g. time of peak velocity)
* - ..
* - ...
*/
DriverXmdf(); DriverXmdf();
~DriverXmdf( ) override = default; ~DriverXmdf( ) override = default;
DriverXmdf *create() override; DriverXmdf *create() override;
@ -77,6 +97,7 @@ namespace MDAL
void addDatasetGroupsFromXmdfGroup( void addDatasetGroupsFromXmdfGroup(
DatasetGroups &groups, DatasetGroups &groups,
const HdfGroup &rootGroup, const HdfGroup &rootGroup,
const std::string &nameSuffix,
size_t vertexCount, size_t vertexCount,
size_t faceCount ); size_t faceCount );
}; };

2
external/mdal/mdal.cpp vendored
View File

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

24
external/mdal/mdal_utils.cpp vendored
View File

@ -24,6 +24,9 @@ bool MDAL::fileExists( const std::string &filename )
bool MDAL::startsWith( const std::string &str, const std::string &substr, ContainsBehaviour behaviour ) bool MDAL::startsWith( const std::string &str, const std::string &substr, ContainsBehaviour behaviour )
{ {
if ( str.size() < substr.size() )
return false;
if ( behaviour == ContainsBehaviour::CaseSensitive ) if ( behaviour == ContainsBehaviour::CaseSensitive )
return str.rfind( substr, 0 ) == 0; return str.rfind( substr, 0 ) == 0;
else else
@ -32,6 +35,9 @@ bool MDAL::startsWith( const std::string &str, const std::string &substr, Contai
bool MDAL::endsWith( const std::string &str, const std::string &substr, ContainsBehaviour behaviour ) bool MDAL::endsWith( const std::string &str, const std::string &substr, ContainsBehaviour behaviour )
{ {
if ( str.size() < substr.size() )
return false;
if ( behaviour == ContainsBehaviour::CaseSensitive ) if ( behaviour == ContainsBehaviour::CaseSensitive )
return str.rfind( substr ) == ( str.size() - substr.size() ); return str.rfind( substr ) == ( str.size() - substr.size() );
else else
@ -237,6 +243,24 @@ std::string MDAL::replace( const std::string &str, const std::string &substr, co
return res; return res;
} }
// http://www.cplusplus.com/faq/sequences/strings/trim/
std::string MDAL::trim( const std::string &s, const std::string &delimiters )
{
return ltrim( rtrim( s, delimiters ), delimiters );
}
// http://www.cplusplus.com/faq/sequences/strings/trim/
std::string MDAL::ltrim( const std::string &s, const std::string &delimiters )
{
return s.substr( s.find_first_not_of( delimiters ) );
}
// http://www.cplusplus.com/faq/sequences/strings/trim/
std::string MDAL::rtrim( const std::string &s, const std::string &delimiters )
{
return s.substr( 0, s.find_last_not_of( delimiters ) + 1 );
}
MDAL::BBox MDAL::computeExtent( const MDAL::Vertices &vertices ) MDAL::BBox MDAL::computeExtent( const MDAL::Vertices &vertices )
{ {
BBox b; BBox b;

27
external/mdal/mdal_utils.hpp vendored
View File

@ -67,29 +67,20 @@ namespace MDAL
std::vector<std::string> split( const std::string &str, const std::string &delimiter, SplitBehaviour behaviour ); std::vector<std::string> split( const std::string &str, const std::string &delimiter, SplitBehaviour behaviour );
std::string join( const std::vector<std::string> parts, const std::string &delimiter ); std::string join( const std::vector<std::string> parts, const std::string &delimiter );
// http://www.cplusplus.com/faq/sequences/strings/trim/ //! Right trim
inline std::string rtrim( std::string rtrim(
const std::string &s, const std::string &s,
const std::string &delimiters = " \f\n\r\t\v" ) const std::string &delimiters = " \f\n\r\t\v" );
{
return s.substr( 0, s.find_last_not_of( delimiters ) + 1 );
}
// http://www.cplusplus.com/faq/sequences/strings/trim/ //! Left trim
inline std::string ltrim( std::string ltrim(
const std::string &s, const std::string &s,
const std::string &delimiters = " \f\n\r\t\v" ) const std::string &delimiters = " \f\n\r\t\v" );
{
return s.substr( s.find_first_not_of( delimiters ) );
}
// http://www.cplusplus.com/faq/sequences/strings/trim/ //! Right and left trim
inline std::string trim( std::string trim(
const std::string &s, const std::string &s,
const std::string &delimiters = " \f\n\r\t\v" ) const std::string &delimiters = " \f\n\r\t\v" );
{
return ltrim( rtrim( s, delimiters ), delimiters );
}
// extent // extent
BBox computeExtent( const Vertices &vertices ); BBox computeExtent( const Vertices &vertices );

1
src/providers/mdal/CMakeLists.txt
View File

@ -88,6 +88,7 @@ IF (WITH_INTERNAL_MDAL)
SET(MDAL_LIB_SRCS ${MDAL_LIB_SRCS} SET(MDAL_LIB_SRCS ${MDAL_LIB_SRCS}
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_cf.cpp ${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_cf.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_3di.cpp ${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_3di.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_netcdf.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_sww.cpp ${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_sww.cpp
) )
SET(MDAL_LIB_HDRS ${MDAL_LIB_HDRS} SET(MDAL_LIB_HDRS ${MDAL_LIB_HDRS}