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

update MDAL to 0.4.91 (alpha release of 0.5.0 for QGIS 3.12) - stacked meshes 3d

Browse Source
This commit is contained in:
Peter Petrik 2019-11-29 15:15:01 +01:00
parent a30bda6c3d
commit 8aa0c28070
46 changed files with 2622 additions and 665 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
external/mdal/api/mdal.h vendored
View File

@ -6,12 +6,6 @@
#ifndef MDAL_H
#define MDAL_H
/**********************************************************************/
/**********************************************************************/
/* API is considered EXPERIMENTAL and can be changed without a notice */
/**********************************************************************/
/**********************************************************************/
#ifdef MDAL_STATIC
# define MDAL_EXPORT
#else
@ -65,12 +59,28 @@ enum MDAL_Status
Warn_NodeNotUnique
};
/**
* Specifies where the data is defined
*/
enum MDAL_DataLocation
{
//! Unknown/Invalid location
DataInvalidLocation = 0,
//! Data is defined on vertices of 2D mesh
DataOnVertices2D,
//! Data is defined on face centres of 2D mesh
DataOnFaces2D,
//! Data is defined on volume centres of 3D mesh
DataOnVolumes3D
};
typedef void *MeshH;
typedef void *MeshVertexIteratorH;
typedef void *MeshFaceIteratorH;
typedef void *DatasetGroupH;
typedef void *DatasetH;
typedef void *DriverH;
typedef void *AveragingMethodH;
//! Returns MDAL version
MDAL_EXPORT const char *MDAL_Version();
@ -104,10 +114,10 @@ MDAL_EXPORT DriverH MDAL_driverFromName( const char *name );
MDAL_EXPORT bool MDAL_DR_meshLoadCapability( DriverH driver );
//! Returns whether driver has capability to write/edit dataset (groups)
MDAL_EXPORT bool MDAL_DR_writeDatasetsCapability( DriverH driver );
MDAL_EXPORT bool MDAL_DR_writeDatasetsCapability( DriverH driver, MDAL_DataLocation location );
//! Returns whether driver has capability to save mesh
MDAL_EXPORT bool MDAL_DR_SaveMeshCapability( DriverH driver );
MDAL_EXPORT bool MDAL_DR_saveMeshCapability( DriverH driver );
/**
* Returns name of MDAL driver
@ -188,14 +198,14 @@ MDAL_EXPORT DatasetGroupH MDAL_M_datasetGroup( MeshH mesh, int index );
* \param mesh mesh handle
* \param driver the driver to use for storing the data
* \param name dataset group name
* \param isOnVertices whether data is defined on vertices
* \param dataLocation location of data (face, vertex, volume)
* \param hasScalarData whether data is scalar (false = vector data)
* \param datasetGroupFile file to store the new dataset group
* \returns empty pointer if not possible to create group, otherwise handle to new group
*/
MDAL_EXPORT DatasetGroupH MDAL_M_addDatasetGroup( MeshH mesh,
const char *name,
bool isOnVertices,
MDAL_DataLocation dataLocation,
bool hasScalarData,
DriverH driver,
const char *datasetGroupFile );
@ -312,7 +322,7 @@ MDAL_EXPORT const char *MDAL_G_driverName( DatasetGroupH group );
MDAL_EXPORT bool MDAL_G_hasScalarData( DatasetGroupH group );
//! Whether dataset is on vertices
MDAL_EXPORT bool MDAL_G_isOnVertices( DatasetGroupH group );
MDAL_EXPORT MDAL_DataLocation MDAL_G_dataLocation( DatasetGroupH group );
/**
* Returns the minimum and maximum values of the group
@ -372,7 +382,18 @@ MDAL_EXPORT DatasetGroupH MDAL_D_group( DatasetH dataset );
//! Returns dataset time
MDAL_EXPORT double MDAL_D_time( DatasetH dataset );
//! Returns number of values
//! Returns volumes count for the mesh (for 3D meshes)
MDAL_EXPORT int MDAL_D_volumesCount( DatasetH dataset );
//! Returns maximum number of vertical levels (for 3D meshes)
MDAL_EXPORT int MDAL_D_maximumVerticalLevelCount( DatasetH dataset );
/**
* Returns number of values
* For dataset with data location DataOnVertices2D returns vertex count
* For dataset with data location DataOnFaces2D returns face count
* For dataset with data location DataOnVolumes3D returns volumes count
*/
MDAL_EXPORT int MDAL_D_valueCount( DatasetH dataset );
//! Returns whether dataset is valid
@ -381,9 +402,15 @@ MDAL_EXPORT bool MDAL_D_isValid( DatasetH dataset );
//! Data type to be returned by MDAL_D_data
enum MDAL_DataType
{
SCALAR_DOUBLE, //!< Double value for scalar datasets
VECTOR_2D_DOUBLE, //!< Double, double value for vector datasets
ACTIVE_INTEGER //!< Integer, active flag for dataset faces. Some formats support switching off the element for particular timestep.
SCALAR_DOUBLE = 0, //!< Double value for scalar datasets (DataOnVertices2D or DataOnFaces2D)
VECTOR_2D_DOUBLE, //!< Double, double value for vector datasets (DataOnVertices2D or DataOnFaces2D)
ACTIVE_INTEGER, //!< Integer, active flag for dataset faces. Some formats support switching off the element for particular timestep (DataOnVertices2D or DataOnFaces2D)
VERTICAL_LEVEL_COUNT_INTEGER, //!< Number of vertical level for particular mesh's face in 3D Stacked Meshes (DataOnVolumes3D)
VERTICAL_LEVEL_DOUBLE, //!< Vertical level extrusion for particular mesh's face in 3D Stacked Meshes (DataOnVolumes3D)
FACE_INDEX_TO_VOLUME_INDEX_INTEGER, //!< The first index of 3D volume for particular mesh's face in 3D Stacked Meshes (DataOnVolumes3D)
SCALAR_VOLUMES_DOUBLE, //!< Double scalar values for volumes in 3D Stacked Meshes (DataOnVolumes3D)
VECTOR_2D_VOLUMES_DOUBLE, //!< Double, double value for volumes in 3D Stacked Meshes (DataOnVolumes3D)
ACTIVE_VOLUMES_INTEGER //!< Integer, active flag for dataset volumes. Some formats support switching off the element for particular timestep (DataOnVolumes3D)
};
/**
@ -399,6 +426,12 @@ enum MDAL_DataType
* For VECTOR_2D_DOUBLE, the minimum size must be valuesCount * 2 * size_of(double).
* Values are returned as x1, y1, x2, y2, ..., xN, yN
* For ACTIVE_INTEGER, the minimum size must be valuesCount * size_of(int)
* For VERTICAL_LEVEL_COUNT_INTEGER, the minimum size must be faceCount * size_of(int)
* For VERTICAL_LEVEL_DOUBLE, the minimum size must be (faceCount + volumesCount) * size_of(double)
* For FACE_INDEX_TO_VOLUME_INDEX_INTEGER, the minimum size must be faceCount * size_of(int)
* For SCALAR_VOLUMES_DOUBLE, the minimum size must be volumesCount * size_of(double)
* For VECTOR_2D_VOLUMES_DOUBLE, the minimum size must be 2 * volumesCount * size_of(double)
* For ACTIVE_VOLUMES_INTEGER, , the minimum size must be volumesCount * size_of(int)
* \returns number of values written to buffer. If return value != count requested, see MDAL_LastStatus() for error type
*/
MDAL_EXPORT int MDAL_D_data( DatasetH dataset, int indexStart, int count, MDAL_DataType dataType, void *buffer );

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

@ -94,7 +94,7 @@ MDAL::Driver2dm *MDAL::Driver2dm::create()
MDAL::Driver2dm::~Driver2dm() = default;
bool MDAL::Driver2dm::canRead( const std::string &uri )
bool MDAL::Driver2dm::canReadMesh( const std::string &uri )
{
std::ifstream in( uri, std::ifstream::in );
std::string line;

2
external/mdal/frmts/mdal_2dm.hpp vendored
View File

@ -86,7 +86,7 @@ namespace MDAL
int faceVerticesMaximumCount() const override
{return MAX_VERTICES_PER_FACE_2DM;}
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile, MDAL_Status *status ) override;
void save( const std::string &uri, Mesh *mesh, MDAL_Status *status ) override;

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

@ -11,7 +11,8 @@ MDAL::Driver3Di::Driver3Di()
: DriverCF(
"3Di",
"3Di Results",
"results_3di.nc" )
"results_3di.nc",
Capability::ReadMesh )
{
}
@ -27,17 +28,17 @@ MDAL::CFDimensions MDAL::Driver3Di::populateDimensions( )
int ncid;
// 2D Mesh
mNcFile.getDimension( "nMesh2D_nodes", &count, &ncid );
mNcFile->getDimension( "nMesh2D_nodes", &count, &ncid );
dims.setDimension( CFDimensions::Face2D, count, ncid );
mNcFile.getDimension( "nCorner_Nodes", &count, &ncid );
mNcFile->getDimension( "nCorner_Nodes", &count, &ncid );
dims.setDimension( CFDimensions::MaxVerticesInFace, count, ncid );
// Vertices count is populated later in populateFacesAndVertices
// it is not known from the array dimensions
// Time
mNcFile.getDimension( "time", &count, &ncid );
mNcFile->getDimension( "time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );
return dims;
@ -53,16 +54,16 @@ void MDAL::Driver3Di::populateFacesAndVertices( Vertices &vertices, Faces &faces
std::map<std::string, size_t> xyToVertex2DId;
// X coordinate
int ncidX = mNcFile.getVarId( "Mesh2DContour_x" );
double fillX = mNcFile.getFillValue( ncidX );
int ncidX = mNcFile->getVarId( "Mesh2DContour_x" );
double fillX = mNcFile->getFillValue( ncidX );
std::vector<double> faceVerticesX( arrsize );
if ( nc_get_var_double( mNcFile.handle(), ncidX, faceVerticesX.data() ) ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_get_var_double( mNcFile->handle(), ncidX, faceVerticesX.data() ) ) throw MDAL_Status::Err_UnknownFormat;
// Y coordinate
int ncidY = mNcFile.getVarId( "Mesh2DContour_y" );
double fillY = mNcFile.getFillValue( ncidY );
int ncidY = mNcFile->getVarId( "Mesh2DContour_y" );
double fillY = mNcFile->getFillValue( ncidY );
std::vector<double> faceVerticesY( arrsize );
if ( nc_get_var_double( mNcFile.handle(), ncidY, faceVerticesY.data() ) ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_get_var_double( mNcFile->handle(), ncidY, faceVerticesY.data() ) ) throw MDAL_Status::Err_UnknownFormat;
// now populate create faces and backtrack which vertices
// are used in multiple faces
@ -120,10 +121,10 @@ void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )
size_t faceCount = mesh->facesCount();
// read Z coordinate of 3di computation nodes centers
int ncidZ = mNcFile.getVarId( "Mesh2DFace_zcc" );
double fillZ = mNcFile.getFillValue( ncidZ );
int ncidZ = mNcFile->getVarId( "Mesh2DFace_zcc" );
double fillZ = mNcFile->getFillValue( ncidZ );
std::vector<double> coordZ( faceCount );
if ( nc_get_var_double( mNcFile.handle(), ncidZ, coordZ.data() ) )
if ( nc_get_var_double( mNcFile->handle(), ncidZ, coordZ.data() ) )
return; //error reading the array
@ -134,10 +135,10 @@ void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )
"Bed Elevation"
);
group->setIsOnVertices( false );
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
group->setIsScalar( true );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
dataset->setTime( 0.0 );
double *values = dataset->values();
for ( size_t i = 0; i < faceCount; ++i )
@ -155,6 +156,11 @@ std::string MDAL::Driver3Di::getCoordinateSystemVariableName()
return "projected_coordinate_system";
}
std::string MDAL::Driver3Di::getTimeVariableName() const
{
return "time";
}
std::set<std::string> MDAL::Driver3Di::ignoreNetCDFVariables()
{
std::set<std::string> ignore_variables;
@ -194,10 +200,10 @@ void MDAL::Driver3Di::parseNetCDFVariableMetadata( int varid, const std::string
*is_vector = false;
*is_x = true;
std::string long_name = mNcFile.getAttrStr( "long_name", varid );
std::string long_name = mNcFile->getAttrStr( "long_name", varid );
if ( long_name.empty() )
{
std::string standard_name = mNcFile.getAttrStr( "standard_name", varid );
std::string standard_name = mNcFile->getAttrStr( "standard_name", varid );
if ( standard_name.empty() )
{
name = variableName;

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

@ -48,6 +48,7 @@ namespace MDAL
void populateFacesAndVertices( Vertices &vertices, Faces &faces ) override;
void addBedElevation( MemoryMesh *mesh ) override;
std::string getCoordinateSystemVariableName() override;
std::string getTimeVariableName() const override;
std::set<std::string> ignoreNetCDFVariables() override;
void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) override;

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

@ -28,7 +28,7 @@ MDAL::DriverAsciiDat::DriverAsciiDat( ):
Driver( "ASCII_DAT",
"DAT",
"*.dat",
Capability::ReadDatasets
Capability::ReadDatasets | Capability::WriteDatasetsOnFaces2D | Capability::WriteDatasetsOnVertices2D
)
{
}
@ -40,7 +40,7 @@ MDAL::DriverAsciiDat *MDAL::DriverAsciiDat::create()
MDAL::DriverAsciiDat::~DriverAsciiDat( ) = default;
bool MDAL::DriverAsciiDat::canRead( const std::string &uri )
bool MDAL::DriverAsciiDat::canReadDatasets( const std::string &uri )
{
std::ifstream in( uri, std::ifstream::in );
std::string line;
@ -84,7 +84,7 @@ void MDAL::DriverAsciiDat::loadOldFormat( std::ifstream &in,
groupName
);
group->setIsScalar( !isVector );
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
do
{
@ -106,7 +106,9 @@ void MDAL::DriverAsciiDat::loadOldFormat( std::ifstream &in,
size_t fileNodeCount = toSizeT( items[1] );
size_t meshIdCount = maximumId( mesh ) + 1;
if ( meshIdCount != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
{
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh )
}
}
else if ( cardType == "SCALAR" || cardType == "VECTOR" )
{
@ -127,16 +129,19 @@ void MDAL::DriverAsciiDat::loadOldFormat( std::ifstream &in,
while ( std::getline( in, line ) );
if ( !group || group->datasets.size() == 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
{
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
void MDAL::DriverAsciiDat::loadNewFormat(
std::ifstream &in,
Mesh *mesh,
MDAL_Status *status ) const
{
bool isVector = false;
std::shared_ptr<DatasetGroup> group; // DAT outputs data
@ -145,12 +150,9 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
std::string referenceTime;
// see if it contains face-centered results - supported by BASEMENT
bool faceCentered = false;
if ( contains( groupName, "_els_" ) )
if ( contains( groupName, "_els" ) )
faceCentered = true;
if ( group )
group->setIsOnVertices( !faceCentered );
while ( std::getline( in, line ) )
{
// Replace tabs by spaces,
@ -171,25 +173,25 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
size_t fileNodeCount = toSizeT( items[1] );
size_t meshIdCount = maximumId( mesh ) + 1;
if ( meshIdCount != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
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 );
}
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 );
}
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 );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat )
}
isVector = cardType == "BEGVEC";
@ -200,7 +202,7 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
groupName
);
group->setIsScalar( !isVector );
group->setIsOnVertices( !faceCentered );
group->setDataLocation( faceCentered ? MDAL_DataLocation::DataOnFaces2D : MDAL_DataLocation::DataOnVertices2D );
group->setReferenceTime( referenceTime );
}
else if ( cardType == "ENDDS" )
@ -208,7 +210,7 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
if ( !group )
{
debug( "ENDDS card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat )
}
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
@ -219,7 +221,7 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
if ( !group )
{
debug( "NAME card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat )
}
size_t quoteIdx1 = line.find( '\"' );
@ -236,7 +238,7 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
if ( !group )
{
debug( "TIMEUNITS card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat )
}
group->setMetadata( "TIMEUNITS", items[1] );
@ -334,8 +336,9 @@ void MDAL::DriverAsciiDat::readVertexTimestep(
{
assert( group );
size_t faceCount = mesh->facesCount();
size_t vertexCount = mesh->verticesCount();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MDAL::MemoryDataset2D >( group.get() );
dataset->setTime( t );
int *active = dataset->active();
@ -353,7 +356,7 @@ void MDAL::DriverAsciiDat::readVertexTimestep(
const Mesh2dm *m2dm = dynamic_cast<const Mesh2dm *>( mesh );
double *values = dataset->values();
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 )
{
@ -407,7 +410,7 @@ void MDAL::DriverAsciiDat::readFaceTimestep(
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MDAL::MemoryDataset2D >( group.get() );
dataset->setTime( t );
double *values = dataset->values();
// TODO: hasStatus
@ -444,6 +447,94 @@ void MDAL::DriverAsciiDat::readFaceTimestep(
group->datasets.push_back( dataset );
}
bool MDAL::DriverAsciiDat::persist( MDAL::DatasetGroup *group )
{
assert( ( group->dataLocation() == MDAL_DataLocation::DataOnFaces2D ) ||
( group->dataLocation() == MDAL_DataLocation::DataOnVertices2D ) );
const bool isScalar = group->isScalar();
const bool isOnVertices = group->dataLocation() == MDAL_DataLocation::DataOnVertices2D;
std::string uri = group->uri();
if ( !MDAL::contains( uri, "_els" ) && isOnVertices == false )
{
// Should contain _els in name but it does not
uri.insert( uri.size() - 4, "_els" );
}
std::ofstream out( uri, std::ofstream::out );
// implementation based on information from:
// https://www.xmswiki.com/wiki/SMS:ASCII_Dataset_Files_*.dat
if ( !out )
return true; // Couldn't open the file
const Mesh *mesh = group->mesh();
size_t nodeCount = mesh->verticesCount();
size_t elemCount = mesh->facesCount();
out << "DATASET\n";
out << "OBJTYPE \"mesh2d\"\n";
if ( isScalar )
out << "BEGSCL\n";
else
out << "BEGVEC\n";
out << "ND " << nodeCount << "\n";
out << "NC " << elemCount << "\n";
out << "NAME " "\"" << group->name() << "\"" "\n";
std::string referenceTimeStr = group->referenceTime();
if ( !referenceTimeStr.empty() )
{
// Cutting of the JULIAN prefix
std::vector<std::string> referenceTimeStrWords = split( referenceTimeStr, ' ' );
if ( referenceTimeStrWords.size() > 1 )
out << "RT_JULIAN " << referenceTimeStrWords[1] << "\n";
else
out << "RT_JULIAN " << referenceTimeStr << "\n";
}
out << "TIMEUNITS " << 0 << "\n";
for ( size_t time_index = 0; time_index < group->datasets.size(); ++ time_index )
{
const std::shared_ptr<MDAL::MemoryDataset2D> dataset
= std::dynamic_pointer_cast<MDAL::MemoryDataset2D>( group->datasets[time_index] );
bool hasActiveStatus = isOnVertices && dataset->active();
out << "TS " << hasActiveStatus << " " << std::to_string( dataset->time() ) << "\n";
if ( hasActiveStatus )
{
// Fill the active data
for ( size_t i = 0; i < elemCount; ++i )
{
int active = dataset->active()[i];
out << ( active == 1 ? true : false ) << "\n";
}
}
size_t valuesToWrite = isOnVertices ? nodeCount : elemCount;
for ( size_t i = 0; i < valuesToWrite; ++i )
{
// Read values flags
if ( isScalar )
out << dataset->values()[i] << "\n";
else
{
out << dataset->values()[2 * i] << " " << dataset->values()[2 * i + 1 ] << "\n";
}
}
}
out << "ENDDS";
return false;
}
double MDAL::DriverAsciiDat::convertTimeDataToHours( double time, const std::string &originalTimeDataUnit ) const
{

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

@ -53,8 +53,9 @@ namespace MDAL
~DriverAsciiDat( ) override;
DriverAsciiDat *create() override;
bool canRead( const std::string &uri ) override;
bool canReadDatasets( const std::string &uri ) override;
void load( const std::string &datFile, Mesh *mesh, MDAL_Status *status ) override;
bool persist( DatasetGroup *group ) override;
private:
bool canReadOldFormat( const std::string &line ) const;

81
external/mdal/frmts/mdal_binary_dat.cpp vendored
View File

@ -40,6 +40,13 @@ static const int CT_FLOAT_SIZE = 4;
static const int CF_FLAG_SIZE = 1;
static const int CF_FLAG_INT_SIZE = 4;
static void exit_with_error( MDAL_Status *status, MDAL_Status error, const std::string msg )
{
MDAL::debug( "BINARY DAT: " + msg );
if ( status ) *status = ( error );
}
#define EXIT_WITH_ERROR(error) { if (status) *status = (error); return; }
static bool read( std::ifstream &in, char *s, int n )
@ -75,7 +82,7 @@ MDAL::DriverBinaryDat::DriverBinaryDat():
Driver( "BINARY_DAT",
"Binary DAT",
"*.dat",
Capability::ReadDatasets | Capability::WriteDatasets
Capability::ReadDatasets | Capability::WriteDatasetsOnVertices2D
)
{
}
@ -87,7 +94,7 @@ MDAL::DriverBinaryDat *MDAL::DriverBinaryDat::create()
MDAL::DriverBinaryDat::~DriverBinaryDat() = default;
bool MDAL::DriverBinaryDat::canRead( const std::string &uri )
bool MDAL::DriverBinaryDat::canReadDatasets( const std::string &uri )
{
std::ifstream in( uri, std::ifstream::in | std::ifstream::binary );
int version;
@ -126,8 +133,7 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
// implementation based on information from:
// http://www.xmswiki.com/wiki/SMS:Binary_Dataset_Files_*.dat
if ( !in )
EXIT_WITH_ERROR( MDAL_Status::Err_FileNotFound ); // Couldn't open the file
if ( !in ) return exit_with_error( status, MDAL_Status::Err_FileNotFound, "Couldn't open the file" );
size_t vertexCount = mesh->verticesCount();
size_t elemCount = mesh->facesCount();
@ -148,18 +154,16 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
char istat;
float time;
if ( read( in, reinterpret_cast< char * >( &version ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &version ), 4 ) ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Unable to read version" );
if ( version != CT_VERSION ) // Version should be 3000
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( version != CT_VERSION ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Invalid version " );
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
name(),
mesh,
mDatFile
); // DAT datasets
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
// in TUFLOW results there could be also a special timestep (99999) with maximums
// we will put it into a separate dataset
@ -168,7 +172,7 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
mesh,
mDatFile
);
groupMax->setIsOnVertices( true );
groupMax->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
while ( card != CT_ENDDS )
{
@ -180,24 +184,21 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
switch ( card )
{
case CT_OBJTYPE:
// Object type
if ( read( in, reinterpret_cast< char * >( &objecttype ), 4 ) || objecttype != CT_2D_MESHES )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &objecttype ), 4 ) || objecttype != CT_2D_MESHES ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Invalid object type" );
break;
case CT_SFLT:
// Float size
if ( read( in, reinterpret_cast< char * >( &sflt ), 4 ) || sflt != CT_FLOAT_SIZE )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &sflt ), 4 ) || sflt != CT_FLOAT_SIZE ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read float size" );
break;
case CT_SFLG:
// Flag size
if ( read( in, reinterpret_cast< char * >( &sflg ), 4 ) )
if ( sflg != CF_FLAG_SIZE && sflg != CF_FLAG_INT_SIZE )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read flag size" );
break;
case CT_BEGSCL:
@ -212,36 +213,29 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
case CT_VECTYPE:
// Vector type
if ( read( in, reinterpret_cast< char * >( &vectype ), 4 ) || vectype != 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &vectype ), 4 ) || vectype != 0 ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read vector type" );
break;
case CT_OBJID:
// Object id
if ( read( in, reinterpret_cast< char * >( &objid ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &objid ), 4 ) ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read object id" );
break;
case CT_NUMDATA:
// Num data
if ( read( in, reinterpret_cast< char * >( &numdata ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( numdata != static_cast< int >( vertexCount ) )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
if ( read( in, reinterpret_cast< char * >( &numdata ), 4 ) ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read num data" );
if ( numdata != static_cast< int >( vertexCount ) ) return exit_with_error( status, MDAL_Status::Err_IncompatibleMesh, "invalid num data" );
break;
case CT_NUMCELLS:
// Num data
if ( read( in, reinterpret_cast< char * >( &numcells ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( numcells != static_cast< int >( elemCount ) )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
if ( read( in, reinterpret_cast< char * >( &numcells ), 4 ) ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read num cells" );
if ( numcells != static_cast< int >( elemCount ) ) return exit_with_error( status, MDAL_Status::Err_IncompatibleMesh, "invalid num cells" );
break;
case CT_NAME:
// Name
if ( read( in, reinterpret_cast< char * >( &groupName ), 40 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &groupName ), 40 ) ) return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read name" );
if ( groupName[39] != 0 )
groupName[39] = 0;
group->setName( trim( std::string( groupName ) ) );
@ -251,10 +245,10 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
case CT_RT_JULIAN:
// Reference time
if ( readIStat( in, sflg, &istat ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read reference time" );
if ( read( in, reinterpret_cast< char * >( &time ), 8 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "unable to read reference time" );
referenceTime = static_cast<double>( time );
group->setReferenceTime( "JULIAN " + std::to_string( referenceTime ) );
@ -263,7 +257,7 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
case CT_TIMEUNITS:
// Time unit
if ( read( in, reinterpret_cast< char * >( &timeUnit ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Unable to read time units" );
switch ( timeUnit )
{
@ -289,23 +283,24 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
case CT_TS:
// Time step!
if ( readIStat( in, sflg, &istat ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Invalid time step" );
if ( read( in, reinterpret_cast< char * >( &time ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Invalid time step" );
double t = static_cast<double>( time );
t = convertTimeDataToHours( t, timeUnit );
if ( readVertexTimestep( mesh, group, groupMax, t, istat, sflg, in ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "Unable to read vertex timestep" );
break;
}
}
if ( !group || group->datasets.size() == 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
return exit_with_error( status, MDAL_Status::Err_UnknownFormat, "No datasets" );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
@ -330,7 +325,7 @@ bool MDAL::DriverBinaryDat::readVertexTimestep( const MDAL::Mesh *mesh,
size_t vertexCount = mesh->verticesCount();
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MDAL::MemoryDataset2D >( group.get() );
int *activeFlags = dataset->active();
bool active = true;
@ -401,6 +396,8 @@ static bool writeRawData( std::ofstream &out, const char *s, int n )
bool MDAL::DriverBinaryDat::persist( MDAL::DatasetGroup *group )
{
assert( group->dataLocation() == MDAL_DataLocation::DataOnVertices2D );
std::ofstream out( group->uri(), std::ofstream::out | std::ofstream::binary );
// implementation based on information from:
@ -412,12 +409,6 @@ bool MDAL::DriverBinaryDat::persist( MDAL::DatasetGroup *group )
size_t nodeCount = mesh->verticesCount();
size_t elemCount = mesh->facesCount();
if ( !group->isOnVertices() )
{
// Element outputs not supported in the format
return true;
}
// version card
writeRawData( out, reinterpret_cast< const char * >( &CT_VERSION ), 4 );
@ -465,7 +456,7 @@ bool MDAL::DriverBinaryDat::persist( MDAL::DatasetGroup *group )
for ( size_t time_index = 0; time_index < group->datasets.size(); ++ time_index )
{
const std::shared_ptr<MDAL::MemoryDataset> dataset = std::dynamic_pointer_cast<MDAL::MemoryDataset>( group->datasets[time_index] );
const std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::dynamic_pointer_cast<MDAL::MemoryDataset2D>( group->datasets[time_index] );
writeRawData( out, reinterpret_cast< const char * >( &CT_TS ), 4 );
writeRawData( out, reinterpret_cast< const char * >( &istat ), 1 );

2
external/mdal/frmts/mdal_binary_dat.hpp vendored
View File

@ -27,7 +27,7 @@ namespace MDAL
~DriverBinaryDat( ) override;
DriverBinaryDat *create() override;
bool canRead( const std::string &uri ) override;
bool canReadDatasets( const std::string &uri ) override;
void load( const std::string &datFile, Mesh *mesh, MDAL_Status *status ) override;
bool persist( DatasetGroup *group ) override;

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

@ -9,6 +9,7 @@
#include "math.h"
#include <stdlib.h>
#include <assert.h>
#include <cstring>
#include "mdal_data_model.hpp"
#include "mdal_cf.hpp"
@ -47,22 +48,23 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
// get variable name
char variable_name_c[NC_MAX_NAME];
if ( nc_inq_varname( mNcFile.handle(), varid, variable_name_c ) ) break; // probably we are at the end of available arrays, quit endless loop
if ( nc_inq_varname( mNcFile->handle(), varid, variable_name_c ) ) break; // probably we are at the end of available arrays, quit endless loop
std::string variable_name( variable_name_c );
if ( ignoreVariables.find( variable_name ) == ignoreVariables.end() )
{
// get number of dimensions
int ndims;
if ( nc_inq_varndims( mNcFile.handle(), varid, &ndims ) ) continue;
if ( nc_inq_varndims( mNcFile->handle(), varid, &ndims ) ) continue;
// we parse either time-dependent or time-independent (e.g. Bed/Maximums)
if ( ( ndims < 1 ) || ( ndims > 2 ) ) continue;
int dimids[2];
if ( nc_inq_vardimid( mNcFile.handle(), varid, dimids ) ) continue;
if ( nc_inq_vardimid( mNcFile->handle(), varid, dimids ) ) continue;
int dimid;
size_t nTimesteps;
bool time_is_first_dimension = true;
if ( ndims == 1 )
{
@ -77,10 +79,12 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
*/
if ( mDimensions.type( dimids[0] ) == CFDimensions::Time )
{
time_is_first_dimension = true;
dimid = dimids[1];
}
else if ( mDimensions.type( dimids[1] ) == CFDimensions::Time )
{
time_is_first_dimension = false;
dimid = dimids[0];
}
else
@ -94,8 +98,6 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
if ( !mDimensions.isDatasetType( mDimensions.type( dimid ) ) )
continue;
size_t arr_size = mDimensions.size( mDimensions.type( dimid ) ) * nTimesteps;
// Get name, if it is vector and if it is x or y
std::string name;
bool is_vector = true;
@ -131,7 +133,8 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
}
dsInfo.outputType = mDimensions.type( dimid );
dsInfo.name = name;
dsInfo.arr_size = arr_size;
dsInfo.nValues = mDimensions.size( mDimensions.type( dimid ) );
dsInfo.time_first_dim = time_is_first_dimension;
dsinfo_map[name] = dsInfo;
}
}
@ -158,61 +161,6 @@ static void populate_vals( bool is_vector, double *vals, size_t i,
}
}
std::shared_ptr<MDAL::Dataset> MDAL::DriverCF::createFace2DDataset( std::shared_ptr<DatasetGroup> group, size_t ts, const MDAL::CFDatasetGroupInfo &dsi,
const std::vector<double> &vals_x, const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y )
{
assert( dsi.outputType == CFDimensions::Face2D );
size_t nFaces2D = mDimensions.size( CFDimensions::Face2D );
size_t nLine1D = mDimensions.size( CFDimensions::Line1D );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared<MDAL::MemoryDataset>( group.get() );
for ( size_t i = 0; i < nFaces2D; ++i )
{
size_t idx = ts * nFaces2D + i;
populate_vals( dsi.is_vector,
dataset->values(),
nLine1D + i,
vals_x,
vals_y,
idx,
fill_val_x,
fill_val_y );
}
return dataset;
}
std::shared_ptr<MDAL::Dataset> MDAL::DriverCF::createVertex2DDataset( std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts, const MDAL::CFDatasetGroupInfo &dsi,
const std::vector<double> &vals_x,
const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y )
{
assert( dsi.outputType == CFDimensions::Vertex2D );
size_t nVertices2D = mDimensions.size( CFDimensions::Vertex2D );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared<MDAL::MemoryDataset>( group.get() );
for ( size_t i = 0; i < nVertices2D; ++i )
{
size_t idx = ts * nVertices2D + i;
populate_vals( dsi.is_vector,
dataset->values(),
i,
vals_x,
vals_y,
idx,
fill_val_x,
fill_val_y );
}
return dataset;
}
void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<double> &times, const MDAL::cfdataset_info_map &dsinfo_map )
{
/* PHASE 2 - add dataset groups */
@ -229,9 +177,11 @@ void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
group->setIsScalar( !dsi.is_vector );
if ( dsi.outputType == CFDimensions::Vertex2D )
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
else if ( dsi.outputType == CFDimensions::Face2D )
group->setIsOnVertices( false );
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
else if ( dsi.outputType == CFDimensions::Volume3D )
group->setDataLocation( MDAL_DataLocation::DataOnVolumes3D );
else
{
// unsupported
@ -239,38 +189,38 @@ void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
}
// read X data
double fill_val_x = mNcFile.getFillValue( dsi.ncid_x );
std::vector<double> vals_x( dsi.arr_size, std::numeric_limits<double>::quiet_NaN() );
if ( nc_get_var_double( mNcFile.handle(), dsi.ncid_x, vals_x.data() ) ) CF_THROW_ERR;
double fill_val_x = mNcFile->getFillValue( dsi.ncid_x );
// read Y data if vector
double fill_val_y = std::numeric_limits<double>::quiet_NaN();
std::vector<double> vals_y;
if ( dsi.is_vector )
{
fill_val_y = mNcFile.getFillValue( dsi.ncid_y );
vals_y.resize( dsi.arr_size, std::numeric_limits<double>::quiet_NaN() );
if ( nc_get_var_double( mNcFile.handle(), dsi.ncid_y, vals_y.data() ) ) CF_THROW_ERR;
fill_val_y = mNcFile->getFillValue( dsi.ncid_y );
}
// Create dataset
for ( size_t ts = 0; ts < dsi.nTimesteps; ++ts )
{
double time = times[ts];
std::shared_ptr<MDAL::Dataset> dataset;
double time = times[ts];
if ( dsi.outputType == CFDimensions::Volume3D )
{
dataset = create3DDataset(
group, ts, dsi, fill_val_x, fill_val_y );
}
else
{
dataset = create2DDataset(
group,
ts,
dsi, fill_val_x, fill_val_y
);
}
if ( dsi.outputType == CFDimensions::Face2D )
{
dataset = createFace2DDataset( group, ts, dsi, vals_x, vals_y, fill_val_x, fill_val_y );
}
else // Vertex2D
{
dataset = createVertex2DDataset( group, ts, dsi, vals_x, vals_y, fill_val_x, fill_val_y );
}
if ( dataset )
{
dataset->setTime( time );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
}
@ -295,8 +245,9 @@ void MDAL::DriverCF::parseTime( std::vector<double> &times )
times = std::vector<double>( 1, 0 );
return;
}
times = mNcFile.readDoubleArr( "time", nTimesteps );
std::string units = mNcFile.getAttrStr( "time", "units" );
const std::string timeArrName = getTimeVariableName();
times = mNcFile->readDoubleArr( timeArrName, nTimesteps );
std::string units = mNcFile->getAttrStr( timeArrName, "units" );
double div_by = MDAL::parseTimeUnits( units );
for ( size_t i = 0; i < nTimesteps; ++i )
{
@ -304,21 +255,49 @@ void MDAL::DriverCF::parseTime( std::vector<double> &times )
}
}
std::shared_ptr<MDAL::Dataset> MDAL::DriverCF::create2DDataset( std::shared_ptr<MDAL::DatasetGroup> group, size_t ts, const MDAL::CFDatasetGroupInfo &dsi, double fill_val_x, double fill_val_y )
{
std::shared_ptr<MDAL::CFDataset2D> dataset = std::make_shared<MDAL::CFDataset2D>(
group.get(),
fill_val_x,
fill_val_y,
dsi.ncid_x,
dsi.ncid_y,
dsi.time_first_dim,
dsi.nTimesteps,
dsi.nValues,
ts,
mNcFile
);
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
return std::move( dataset );
}
std::shared_ptr<MDAL::Dataset> MDAL::DriverCF::create3DDataset( std::shared_ptr<MDAL::DatasetGroup>,
size_t, const MDAL::CFDatasetGroupInfo &,
double, double )
{
std::shared_ptr<MDAL::Dataset> dataset;
return dataset;
}
MDAL::DriverCF::DriverCF( const std::string &name,
const std::string &longName,
const std::string &filters ):
Driver( name, longName, filters, Capability::ReadMesh )
const std::string &filters,
const int capabilities ):
Driver( name, longName, filters, capabilities )
{
}
bool MDAL::DriverCF::canRead( const std::string &uri )
MDAL::DriverCF::~DriverCF() = default;
bool MDAL::DriverCF::canReadMesh( const std::string &uri )
{
try
{
NetCDFFile ncFile;
ncFile.openFile( uri );
mNcFile = ncFile;
mNcFile.reset( new NetCDFFile );
mNcFile->openFile( uri );
populateDimensions( );
}
catch ( MDAL_Status )
@ -337,13 +316,13 @@ void MDAL::DriverCF::setProjection( MDAL::Mesh *mesh )
if ( !coordinate_system_variable.empty() )
{
std::string wkt = mNcFile.getAttrStr( coordinate_system_variable, "wkt" );
std::string wkt = mNcFile->getAttrStr( coordinate_system_variable, "wkt" );
if ( wkt.empty() )
{
std::string epsg_code = mNcFile.getAttrStr( coordinate_system_variable, "EPSG_code" );
std::string epsg_code = mNcFile->getAttrStr( coordinate_system_variable, "EPSG_code" );
if ( epsg_code.empty() )
{
int epsg = mNcFile.getAttrInt( coordinate_system_variable, "epsg" );
int epsg = mNcFile->getAttrInt( coordinate_system_variable, "epsg" );
if ( epsg != 0 )
{
mesh->setSourceCrsFromEPSG( epsg );
@ -370,6 +349,8 @@ void MDAL::DriverCF::setProjection( MDAL::Mesh *mesh )
std::unique_ptr< MDAL::Mesh > MDAL::DriverCF::load( const std::string &fileName, MDAL_Status *status )
{
mNcFile.reset( new NetCDFFile );
mFileName = fileName;
if ( status ) *status = MDAL_Status::None;
@ -380,7 +361,7 @@ std::unique_ptr< MDAL::Mesh > MDAL::DriverCF::load( const std::string &fileName,
try
{
// Open file
mNcFile.openFile( mFileName );
mNcFile->openFile( mFileName );
// Parse dimensions
mDimensions = populateDimensions( );
@ -394,14 +375,13 @@ std::unique_ptr< MDAL::Mesh > MDAL::DriverCF::load( const std::string &fileName,
name(),
vertices.size(),
faces.size(),
mDimensions.MaxVerticesInFace,
mDimensions.size( mDimensions.MaxVerticesInFace ),
computeExtent( vertices ),
mFileName
)
);
mesh->faces = faces;
mesh->vertices = vertices;
addBedElevation( mesh.get() );
setProjection( mesh.get() );
@ -457,6 +437,117 @@ bool MDAL::CFDimensions::isDatasetType( MDAL::CFDimensions::Type type ) const
( type == Vertex2D ) ||
( type == Line1D ) ||
( type == Face2DEdge ) ||
( type == Face2D )
( type == Face2D ) ||
( type == Volume3D )
);
}
//////////////////////////////////////////////////////////////////////////////////////
MDAL::CFDataset2D::CFDataset2D( MDAL::DatasetGroup *parent,
double fill_val_x, double fill_val_y,
int ncid_x, int ncid_y, bool time_first_dim,
size_t timesteps, size_t values, size_t ts, std::shared_ptr<NetCDFFile> ncFile )
: Dataset2D( parent )
, mFillValX( fill_val_x )
, mFillValY( fill_val_y )
, mNcidX( ncid_x )
, mNcidY( ncid_y )
, mTimeFirstDim( time_first_dim )
, mTimesteps( timesteps )
, mValues( values )
, mTs( ts )
, mNcFile( ncFile )
{
}
MDAL::CFDataset2D::~CFDataset2D() = default;
size_t MDAL::CFDataset2D::scalarData( size_t indexStart, size_t count, double *buffer )
{
assert( group()->isScalar() ); //checked in C API interface
if ( ( count < 1 ) || ( indexStart >= mValues ) )
return 0;
if ( mTs >= mTimesteps )
return 0;
size_t copyValues = std::min( mValues - indexStart, count );
size_t start_dim1 = mTimeFirstDim ? mTs : indexStart;
size_t start_dim2 = mTimeFirstDim ? indexStart : mTs;
size_t count_dim1 = mTimeFirstDim ? 1 : copyValues;
size_t count_dim2 = mTimeFirstDim ? copyValues : 1;
std::vector<double> values_x = mNcFile->readDoubleArr(
mNcidX,
start_dim1,
start_dim2,
count_dim1,
count_dim2
);
for ( size_t i = 0; i < copyValues; ++i )
{
populate_vals( false,
buffer,
i,
values_x,
std::vector<double>(),
i,
mFillValX,
mFillValY );
}
return copyValues;
}
size_t MDAL::CFDataset2D::vectorData( size_t indexStart, size_t count, double *buffer )
{
assert( !group()->isScalar() ); //checked in C API interface
if ( ( count < 1 ) || ( indexStart >= mValues ) )
return 0;
if ( mTs >= mTimesteps )
return 0;
size_t copyValues = std::min( mValues - indexStart, count );
size_t start_dim1 = mTimeFirstDim ? mTs : indexStart;
size_t start_dim2 = mTimeFirstDim ? indexStart : mTs;
size_t count_dim1 = mTimeFirstDim ? 1 : copyValues;
size_t count_dim2 = mTimeFirstDim ? copyValues : 1;
std::vector<double> values_x = mNcFile->readDoubleArr(
mNcidX,
start_dim1,
start_dim2,
count_dim1,
count_dim2
);
std::vector<double> values_y = mNcFile->readDoubleArr(
mNcidY,
start_dim1,
start_dim2,
count_dim1,
count_dim2
);
for ( size_t i = 0; i < copyValues; ++i )
{
populate_vals( true,
buffer,
i,
values_x,
values_y,
i,
mFillValX,
mFillValY );
}
return copyValues;
}
size_t MDAL::CFDataset2D::activeData( size_t indexStart, size_t count, int *buffer )
{
MDAL_UNUSED( indexStart )
memset( buffer, 1, count * sizeof( int ) );
return count;
}

74
external/mdal/frmts/mdal_cf.hpp vendored
View File

@ -31,6 +31,8 @@ namespace MDAL
Line1D, //!< Line joining 1D vertices
Face2DEdge, //!< Edge of 2D Face
Face2D, //!< 2D (Polygon) Face
Volume3D, //!< 3D (stacked) volumes
StackedFace3D, //! 3D (stacked) Faces
Time, //!< Time steps
MaxVerticesInFace //!< Maximum number of vertices in a face
};
@ -54,13 +56,45 @@ namespace MDAL
std::string name; //!< Dataset group name
CFDimensions::Type outputType;
bool is_vector;
bool time_first_dim;
size_t nTimesteps;
size_t nValues;
int ncid_x; //!< NetCDF variable id
int ncid_y; //!< NetCDF variable id
size_t arr_size;
};
typedef std::map<std::string, CFDatasetGroupInfo> cfdataset_info_map; // name -> DatasetInfo
class CFDataset2D: public Dataset2D
{
public:
CFDataset2D( DatasetGroup *parent,
double fill_val_x,
double fill_val_y,
int ncid_x,
int ncid_y,
bool time_first_dim,
size_t timesteps,
size_t values,
size_t ts,
std::shared_ptr<NetCDFFile> ncFile
);
virtual ~CFDataset2D() override;
virtual size_t scalarData( size_t indexStart, size_t count, double *buffer ) override;
virtual size_t vectorData( size_t indexStart, size_t count, double *buffer ) override;
virtual size_t activeData( size_t indexStart, size_t count, int *buffer ) override;
private:
double mFillValX;
double mFillValY;
int mNcidX; //!< NetCDF variable id
int mNcidY; //!< NetCDF variable id
bool mTimeFirstDim;
size_t mTimesteps;
size_t mValues;
size_t mTs;
std::shared_ptr<NetCDFFile> mNcFile;
};
//! NetCDF Climate and Forecast (CF) Metadata Conventions
//! http://cfconventions.org
//! and http://ugrid-conventions.github.io/ugrid-conventions/
@ -69,9 +103,10 @@ namespace MDAL
public:
DriverCF( const std::string &name,
const std::string &longName,
const std::string &filters );
virtual ~DriverCF() override = default;
bool canRead( const std::string &uri ) override;
const std::string &filters,
const int capabilities );
virtual ~DriverCF() override;
bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &fileName, MDAL_Status *status ) override;
protected:
@ -82,33 +117,28 @@ namespace MDAL
virtual std::set<std::string> ignoreNetCDFVariables() = 0;
virtual void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) = 0;
virtual std::string getTimeVariableName() const = 0;
virtual std::shared_ptr<MDAL::Dataset> create2DDataset(
std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
double fill_val_x, double fill_val_y );
virtual std::shared_ptr<MDAL::Dataset> create3DDataset(
std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
double fill_val_x, double fill_val_y );
void setProjection( MDAL::Mesh *m );
cfdataset_info_map parseDatasetGroupInfo();
void parseTime( std::vector<double> &times );
std::shared_ptr<MDAL::Dataset> createFace2DDataset(
std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
const std::vector<double> &vals_x,
const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y );
std::shared_ptr<MDAL::Dataset> createVertex2DDataset(
std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
const std::vector<double> &vals_x,
const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y );
void addDatasetGroups( Mesh *mesh,
const std::vector<double> &times,
const cfdataset_info_map &dsinfo_map );
std::string mFileName;
NetCDFFile mNcFile;
std::shared_ptr<NetCDFFile> mNcFile;
CFDimensions mDimensions;
};

62
external/mdal/frmts/mdal_driver.cpp vendored
View File

@ -17,7 +17,6 @@ MDAL::Driver::Driver( const std::string &name,
, mFilters( filters )
, mCapabilityFlags( capabilityFlags )
{
}
MDAL::Driver::~Driver() = default;
@ -42,34 +41,34 @@ bool MDAL::Driver::hasCapability( MDAL::Capability capability ) const
return capability == ( mCapabilityFlags & capability );
}
int MDAL::Driver::faceVerticesMaximumCount() const
bool MDAL::Driver::canReadMesh( const std::string & ) { return false; }
bool MDAL::Driver::canReadDatasets( const std::string & ) { return false; }
bool MDAL::Driver::hasWriteDatasetCapability( MDAL_DataLocation location ) const
{
return -1;
switch ( location )
{
case MDAL_DataLocation::DataOnVertices2D:
return hasCapability( MDAL::Capability::WriteDatasetsOnVertices2D );
case MDAL_DataLocation::DataOnFaces2D:
return hasCapability( MDAL::Capability::WriteDatasetsOnFaces2D );
case MDAL_DataLocation::DataOnVolumes3D:
return hasCapability( MDAL::Capability::WriteDatasetsOnVolumes3D );
default:
return false;
}
}
std::unique_ptr< MDAL::Mesh > MDAL::Driver::load( const std::string &uri, MDAL_Status *status )
{
MDAL_UNUSED( uri );
MDAL_UNUSED( status );
return std::unique_ptr< MDAL::Mesh >();
}
int MDAL::Driver::faceVerticesMaximumCount() const { return -1; }
void MDAL::Driver::load( const std::string &uri, Mesh *mesh, MDAL_Status *status )
{
MDAL_UNUSED( uri );
MDAL_UNUSED( mesh );
MDAL_UNUSED( status );
return;
}
std::unique_ptr< MDAL::Mesh > MDAL::Driver::load( const std::string &, MDAL_Status * ) { return std::unique_ptr< MDAL::Mesh >(); }
void MDAL::Driver::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Status *status )
{
MDAL_UNUSED( uri );
MDAL_UNUSED( mesh );
MDAL_UNUSED( status );
}
void MDAL::Driver::load( const std::string &, Mesh *, MDAL_Status * ) {}
void MDAL::Driver::createDatasetGroup( MDAL::Mesh *mesh, const std::string &groupName, bool isOnVertices, bool hasScalarData, const std::string &datasetGroupFile )
void MDAL::Driver::save( const std::string &, MDAL::Mesh *, MDAL_Status * ) {}
void MDAL::Driver::createDatasetGroup( MDAL::Mesh *mesh, const std::string &groupName, MDAL_DataLocation dataLocation, bool hasScalarData, const std::string &datasetGroupFile )
{
std::shared_ptr<MDAL::DatasetGroup> grp(
new MDAL::DatasetGroup( name(),
@ -77,7 +76,7 @@ void MDAL::Driver::createDatasetGroup( MDAL::Mesh *mesh, const std::string &grou
datasetGroupFile )
);
grp->setName( groupName );
grp->setIsOnVertices( isOnVertices );
grp->setDataLocation( dataLocation );
grp->setIsScalar( hasScalarData );
grp->startEditing();
mesh->datasetGroups.push_back( grp );
@ -85,17 +84,18 @@ void MDAL::Driver::createDatasetGroup( MDAL::Mesh *mesh, const std::string &grou
void MDAL::Driver::createDataset( MDAL::DatasetGroup *group, double time, const double *values, const int *active )
{
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group );
dataset->setTime( time );
memcpy( dataset->values(), values, sizeof( double ) * dataset->valuesCount() );
size_t count = dataset->valuesCount();
if ( !group->isScalar() )
count *= 2;
memcpy( dataset->values(), values, sizeof( double ) * count );
if ( active && dataset->active() )
memcpy( dataset->active(), active, sizeof( int ) * dataset->mesh()->facesCount() );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
bool MDAL::Driver::persist( MDAL::DatasetGroup *group )
{
MDAL_UNUSED( group );
return true; // failure
}
bool MDAL::Driver::persist( MDAL::DatasetGroup * ) { return true; } // failure

19
external/mdal/frmts/mdal_driver.hpp vendored
View File

@ -14,11 +14,13 @@ namespace MDAL
{
enum Capability
{
None = 0,
ReadMesh = 1 << 0, //! Can read mesh and all datasets stored in the mesh file
ReadDatasets = 1 << 1, //! Can read only datasets (groups) from existing mesh
WriteDatasets = 1 << 2, //! Can write datasets (groups)
SaveMesh = 1 << 3, //! Can save the mesh
None = 0,
ReadMesh = 1 << 0, //! Can read mesh and all datasets stored in the mesh file
SaveMesh = 1 << 1, //! Can save the mesh
ReadDatasets = 1 << 2, //! Can read only datasets (groups) from existing mesh
WriteDatasetsOnVertices2D = 1 << 3, //! Can write datasets (groups) on MDAL_DataLocation::DataOnVertices2D
WriteDatasetsOnFaces2D = 1 << 4, //! Can write datasets (groups) on MDAL_DataLocation::DataOnFaces2D
WriteDatasetsOnVolumes3D = 1 << 5, //! Can write datasets (groups) on MDAL_DataLocation::DataOnVolumes3D
};
class Driver
@ -37,8 +39,10 @@ namespace MDAL
std::string longName() const;
std::string filters() const;
bool hasCapability( Capability capability ) const;
bool hasWriteDatasetCapability( MDAL_DataLocation location ) const;
virtual bool canRead( const std::string &uri ) = 0;
virtual bool canReadMesh( const std::string &uri );
virtual bool canReadDatasets( const std::string &uri );
//! returns the maximum vertices per face
virtual int faceVerticesMaximumCount() const;
@ -47,7 +51,6 @@ namespace MDAL
virtual std::unique_ptr< Mesh > load( const std::string &uri, MDAL_Status *status );
// loads datasets
virtual void load( const std::string &uri, Mesh *mesh, MDAL_Status *status );
// save mesh
virtual void save( const std::string &uri, Mesh *mesh, MDAL_Status *status );
@ -55,7 +58,7 @@ namespace MDAL
virtual void createDatasetGroup(
Mesh *mesh,
const std::string &groupName,
bool isOnVertices,
MDAL_DataLocation dataLocation,
bool hasScalarData,
const std::string &datasetGroupFile );

294
external/mdal/frmts/mdal_esri_tin.cpp vendored
View File

@ -23,13 +23,155 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverEsriTin::load( const std::string &uri, M
try
{
std::list<int> superpointIndexes;
Vertices vertices;
Faces faces;
bool isNativeLittleEndian = MDAL::isNativeLittleEndian();
readSuperpoints( uri, superpointIndexes );
populateVertices( uri, vertices, superpointIndexes );
populateFaces( uri, faces, superpointIndexes );
//read the total number of vertices (including superpoints and isolated vertices)
int32_t totalIndexesCount32;
std::ifstream inDenv( denvFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( !inDenv.is_open() )
{
inDenv.open( denv9File( uri ), std::ifstream::in | std::ifstream::binary );
if ( !inDenv.is_open() )
throw MDAL_Status::Err_UnknownFormat;
}
readValue( totalIndexesCount32, inDenv, isNativeLittleEndian );
size_t totalIndexesCount = static_cast<size_t>( totalIndexesCount32 );
/* Round 1 :populates faces with raw indexes from the file
* rawAndCorrectedIndexesMap is used to map raw indexes from the files and corrected indexes
* Corrected indexes take into account the unwanted vertexes (superpoints, isolated vertices)
* Wanted vertices are associated with the corrected index
* Unwanted vertices are associated with the totalIndexesCount value
*/
std::vector<size_t> rawAndCorrectedIndexesMap( totalIndexesCount, totalIndexesCount );
std::ifstream inFaces( faceFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inMsk( mskFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inMsx( msxFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! inFaces.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inMsk.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inMsx.is_open() )
throw MDAL_Status::Err_FileNotFound;
//Find the beginning of data in mskFile
inMsx.seekg( -4, std::ios::end );
int32_t mskBegin;
if ( ! readValue( mskBegin, inMsx, true ) )
throw MDAL_Status::Err_UnknownFormat;
//read information in mskFile
inMsk.seekg( -mskBegin * 2, std::ios::end );
int32_t maskIntergerCount;
if ( ! readValue( maskIntergerCount, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
inMsk.ignore( 4 ); //unused 4 bytes
int32_t maskBitsCount;
if ( ! readValue( maskBitsCount, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
int c = 0;
int32_t maskInt = 0;
while ( true )
{
//read mask file
if ( c % 32 == 0 && c < maskBitsCount ) //first bit in the mask array have to be used-->read next maskInt
if ( ! readValue( maskInt, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
Face f;
for ( int i = 0; i < 3; ++i )
{
int32_t index;
if ( ! readValue( index, inFaces, isNativeLittleEndian ) )
break;
f.push_back( static_cast<size_t>( index - 1 ) );
}
if ( f.size() == 0 ) //that's mean this is the end of the file
break;
if ( f.size() < 3 ) //that's mean the face is not complete
throw MDAL_Status::Err_UnknownFormat;
//exclude masked face
if ( !( maskInt & 0x01 ) )
{
faces.push_back( f );
//fill raw indexes
for ( auto ri : f )
{
if ( ri >= totalIndexesCount )
throw MDAL_Status::Err_UnknownFormat;
rawAndCorrectedIndexesMap[ri] = 1;
}
}
c++;
maskInt = maskInt >> 1;
}
inFaces.close();
inMsk.close();
inMsx.close();
//Round 2 :count the number of wanted indexes and fill the rawIndexes value with the correctedIndex
size_t correctedIndexCount = 0;
for ( size_t i = 0; i < rawAndCorrectedIndexesMap.size(); ++i )
{
if ( rawAndCorrectedIndexesMap.at( i ) < totalIndexesCount )
{
rawAndCorrectedIndexesMap[i] = correctedIndexCount;
correctedIndexCount++;
}
}
//Round 3: populate vertices
Vertices vertices( correctedIndexCount );
std::ifstream inXY( xyFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inZ( zFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! inXY.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inZ.is_open() )
throw MDAL_Status::Err_FileNotFound;
size_t rawIndex = 0;
while ( rawIndex < rawAndCorrectedIndexesMap.size() )
{
Vertex vert;
if ( !readValue( vert.x, inXY, isNativeLittleEndian ) )
break;
if ( !readValue( vert.y, inXY, isNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
float zValue;
if ( !readValue( zValue, inZ, isNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
vert.z = double( zValue );
// store the vertex only if it is a wanted index
if ( rawAndCorrectedIndexesMap[rawIndex] < totalIndexesCount )
vertices[rawAndCorrectedIndexesMap[rawIndex]] = vert ;
rawIndex++;
}
inXY.close();
inZ.close();
//Round 4 :apply correction to the face's indexes
for ( auto &face : faces )
for ( auto &fi : face )
fi = rawAndCorrectedIndexesMap[fi];
//create the memory mesh
std::unique_ptr< MemoryMesh > mesh(
new MemoryMesh(
name(),
@ -41,9 +183,11 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverEsriTin::load( const std::string &uri, M
)
);
//move the faces and the vertices in the mesh
mesh->faces = std::move( faces );
mesh->vertices = std::move( vertices );
//create the "Altitude" dataset
addBedElevationDatasetGroup( mesh.get(), mesh->vertices );
mesh->datasetGroups.back()->setName( "Altitude" );
@ -60,7 +204,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverEsriTin::load( const std::string &uri, M
}
}
bool MDAL::DriverEsriTin::canRead( const std::string &uri )
bool MDAL::DriverEsriTin::canReadMesh( const std::string &uri )
{
std::string zFileName = zFile( uri );
@ -115,6 +259,16 @@ std::string MDAL::DriverEsriTin::hullFile( const std::string &uri ) const
return pathJoin( dirName( uri ), "thul.adf" );
}
std::string MDAL::DriverEsriTin::denvFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tdenv.adf" );
}
std::string MDAL::DriverEsriTin::denv9File( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tdenv9.adf" );
}
std::string MDAL::DriverEsriTin::crsFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "prj.adf" );
@ -152,114 +306,6 @@ std::string MDAL::DriverEsriTin::getTinName( const std::string &uri ) const
return tinName;
}
void MDAL::DriverEsriTin::populateVertices( const std::string &uri, MDAL::Vertices &vertices, const std::list<int> &superpointIndexes ) const
{
//first, read the (x,y) coordinates
bool isNativeLittleEndian = MDAL::isNativeLittleEndian();
std::ifstream inXY( xyFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inZ( zFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! inXY.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inZ.is_open() )
throw MDAL_Status::Err_FileNotFound;
int fileIndex = 1;
while ( true )
{
Vertex vert;
if ( !readValue( vert.x, inXY, isNativeLittleEndian ) )
break;
if ( !readValue( vert.y, inXY, isNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
float zValue;
if ( !readValue( zValue, inZ, isNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
vert.z = double( zValue );
if ( correctedIndex( fileIndex, superpointIndexes ) >= 0 ) //store the vertex only if it is not a superpoint
vertices.push_back( vert );
fileIndex++;
}
inXY.close();
inZ.close();
}
void MDAL::DriverEsriTin::populateFaces( const std::string &uri, MDAL::Faces &faces, const std::list<int> &superpointIndexes ) const
{
bool isNativeLittleEndian = MDAL::isNativeLittleEndian();
std::ifstream inFaces( faceFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inMsk( mskFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inMsx( msxFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! inFaces.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inMsk.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inMsx.is_open() )
throw MDAL_Status::Err_FileNotFound;
//Find the beginning of data in mskFile
inMsx.seekg( -4, std::ios::end );
int32_t mskBegin;
if ( ! readValue( mskBegin, inMsx, true ) )
throw MDAL_Status::Err_UnknownFormat;
//read information in mskFile
inMsk.seekg( -mskBegin * 2, std::ios::end );
int32_t maskIntergerCount;
if ( ! readValue( maskIntergerCount, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
inMsk.ignore( 4 ); //unused 4 bytes
int32_t maskBitsCount;
if ( ! readValue( maskBitsCount, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
int c = 0;
int32_t maskInt = 0;
while ( true )
{
//read mask file
if ( c % 32 == 0 && c < maskBitsCount ) //first bit in the mask array have to be used-->read next maskInt
if ( ! readValue( maskInt, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
Face f;
for ( int i = 0; i < 3; ++i )
{
int32_t index;
if ( ! readValue( index, inFaces, isNativeLittleEndian ) )
break;
index = correctedIndex( index, superpointIndexes );
f.push_back( static_cast<size_t>( index ) );
}
if ( f.size() == 0 ) //that's mean this is the end of the file
break;
if ( f.size() < 3 ) //that's mean the face is not complete
throw MDAL_Status::Err_UnknownFormat;
//exclude masked face
if ( !( maskInt & 0x01 ) )
faces.push_back( f );
c++;
maskInt = maskInt >> 1;
}
inFaces.close();
inMsk.close();
inMsx.close();
}
std::string MDAL::DriverEsriTin::getCrsWkt( const std::string &uri ) const
{
@ -276,25 +322,3 @@ std::string MDAL::DriverEsriTin::getCrsWkt( const std::string &uri ) const
return crsWkt;
}
int MDAL::DriverEsriTin::correctedIndex( int i, const std::list<int> &superPointsIndexes ) const
{
//usualy i is bigger than the biggest superpoint's index
if ( i > superPointsIndexes.back() )
return i - int( superPointsIndexes.size() ) - 1;
//correction of index depending of the position of the superpoints
int correctedIndex = i - 1;
for ( auto si : superPointsIndexes )
{
if ( i == si )
return -1; //return -1 if the index is associated with superpoint
else if ( i > si )
{
correctedIndex--;
}
else
break;
}
return correctedIndex;
}

11
external/mdal/frmts/mdal_esri_tin.hpp vendored
View File

@ -73,7 +73,7 @@ namespace MDAL
virtual std::unique_ptr< Mesh > load( const std::string &uri, MDAL_Status *status ) override;
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
private:
std::string xyFile( const std::string &uri ) const;
@ -82,18 +82,17 @@ namespace MDAL
std::string mskFile( const std::string &uri ) const;
std::string msxFile( const std::string &uri ) const;
std::string hullFile( const std::string &uri ) const;
std::string denvFile( const std::string &uri ) const;
std::string denv9File( const std::string &uri ) const;
std::string crsFile( const std::string &uri ) const;
//* can be used to read superpoints indexes, currently unused in MDAL
void readSuperpoints( const std::string &uri, std::list<int> &superpointsIndexes ) const;
void populateVertices( const std::string &uri, Vertices &vertices, const std::list<int> &superpointIndexes ) const;
void populateFaces( const std::string &uri, Faces &faces, const std::list<int> &superpointIndexes ) const;
std::string getCrsWkt( const std::string &uri ) const;
std::string getTinName( const std::string &uri ) const;
// correction of vertex index because indexing in ESRI TIN files begins with 1 and to take account superpoints position
// return -1 if i is the index of a superpoint
int correctedIndex( int i, const std::list<int> &superPointsIndexes ) const;
};
}

278
external/mdal/frmts/mdal_flo2d.cpp vendored
View File

@ -1,5 +1,5 @@
/*
MDAL - mMesh Data Abstraction Library (MIT License)
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2016 Lutra Consulting
Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
*/
@ -14,10 +14,13 @@
#include <string>
#include <cmath>
#include <cstring>
#include <assert.h>
#include "mdal_utils.hpp"
#include "mdal_hdf5.hpp"
#define FLO2D_NAN 0.0
struct VertexCompare
{
bool operator()( const MDAL::Vertex &lhs, const MDAL::Vertex &rhs ) const
@ -42,7 +45,7 @@ static std::string fileNameFromDir( const std::string &mainFileName, const std::
static double getDouble( double val )
{
if ( MDAL::equals( val, 0.0, 1e-8 ) )
if ( MDAL::equals( val, FLO2D_NAN, 1e-8 ) )
{
return MDAL_NAN;
}
@ -52,6 +55,18 @@ static double getDouble( double val )
}
}
static double toFlo2DDouble( double val )
{
if ( std::isnan( val ) )
{
return FLO2D_NAN;
}
else
{
return val;
}
}
static double getDouble( const std::string &val )
{
double valF = MDAL::toDouble( val );
@ -69,10 +84,10 @@ void MDAL::DriverFlo2D::addStaticDataset(
datFileName,
groupName
);
group->setIsOnVertices( false );
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
group->setIsScalar( true );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
assert( vals.size() == dataset->valuesCount() );
dataset->setTime( 0.0 );
double *values = dataset->values();
@ -143,7 +158,7 @@ void MDAL::DriverFlo2D::parseFPLAINFile( std::vector<double> &elevations,
}
}
static void addDatasetToGroup( std::shared_ptr<MDAL::DatasetGroup> group, std::shared_ptr<MDAL::MemoryDataset> dataset )
static void addDatasetToGroup( std::shared_ptr<MDAL::DatasetGroup> group, std::shared_ptr<MDAL::MemoryDataset2D> dataset )
{
if ( group && dataset && dataset->valuesCount() > 0 )
{
@ -181,7 +196,7 @@ void MDAL::DriverFlo2D::parseTIMDEPFile( const std::string &datFileName, const s
datFileName,
"Depth"
);
depthDsGroup->setIsOnVertices( false );
depthDsGroup->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
depthDsGroup->setIsScalar( true );
@ -191,7 +206,7 @@ void MDAL::DriverFlo2D::parseTIMDEPFile( const std::string &datFileName, const s
datFileName,
"Water Level"
);
waterLevelDsGroup->setIsOnVertices( false );
waterLevelDsGroup->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
waterLevelDsGroup->setIsScalar( true );
std::shared_ptr<DatasetGroup> flowDsGroup = std::make_shared< DatasetGroup >(
@ -200,12 +215,12 @@ void MDAL::DriverFlo2D::parseTIMDEPFile( const std::string &datFileName, const s
datFileName,
"Velocity"
);
flowDsGroup->setIsOnVertices( false );
flowDsGroup->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
flowDsGroup->setIsScalar( false );
std::shared_ptr<MDAL::MemoryDataset> flowDataset;
std::shared_ptr<MDAL::MemoryDataset> depthDataset;
std::shared_ptr<MDAL::MemoryDataset> waterLevelDataset;
std::shared_ptr<MDAL::MemoryDataset2D> flowDataset;
std::shared_ptr<MDAL::MemoryDataset2D> depthDataset;
std::shared_ptr<MDAL::MemoryDataset2D> waterLevelDataset;
while ( std::getline( inStream, line ) )
{
@ -220,9 +235,9 @@ void MDAL::DriverFlo2D::parseTIMDEPFile( const std::string &datFileName, const s
if ( flowDataset ) addDatasetToGroup( flowDsGroup, flowDataset );
if ( waterLevelDataset ) addDatasetToGroup( waterLevelDsGroup, waterLevelDataset );
depthDataset = std::make_shared< MemoryDataset >( depthDsGroup.get() );
flowDataset = std::make_shared< MemoryDataset >( flowDsGroup.get() );
waterLevelDataset = std::make_shared< MemoryDataset >( waterLevelDsGroup.get() );
depthDataset = std::make_shared< MemoryDataset2D >( depthDsGroup.get() );
flowDataset = std::make_shared< MemoryDataset2D >( flowDsGroup.get() );
waterLevelDataset = std::make_shared< MemoryDataset2D >( waterLevelDsGroup.get() );
depthDataset->setTime( time );
flowDataset->setTime( time );
@ -495,16 +510,14 @@ void MDAL::DriverFlo2D::createMesh( const std::vector<CellCenter> &cells, double
mMesh->vertices = vertices;
}
bool MDAL::DriverFlo2D::parseHDF5Datasets( const std::string &datFileName )
bool MDAL::DriverFlo2D::parseHDF5Datasets( MemoryMesh *mesh, const std::string &timedepFileName )
{
//return true on error
size_t nFaces = mMesh->facesCount();
std::string timedepFileName = fileNameFromDir( datFileName, "TIMDEP.HDF5" );
size_t nFaces = mesh->facesCount();
if ( !fileExists( timedepFileName ) ) return true;
HdfFile file( timedepFileName );
HdfFile file( timedepFileName, HdfFile::ReadOnly );
if ( !file.isValid() ) return true;
HdfGroup timedataGroup = file.group( "TIMDEP NETCDF OUTPUT RESULTS" );
@ -538,7 +551,7 @@ bool MDAL::DriverFlo2D::parseHDF5Datasets( const std::string &datFileName )
bool isVector = MDAL::contains( groupType.readString(), "vector", ContainsBehaviour::CaseInsensitive );
// Some sanity checks
size_t expectedSize = mMesh->facesCount() * timesteps;
size_t expectedSize = mesh->facesCount() * timesteps;
if ( isVector ) expectedSize *= 2;
if ( valuesDs.elementCount() != expectedSize ) return true;
@ -549,16 +562,16 @@ bool MDAL::DriverFlo2D::parseHDF5Datasets( const std::string &datFileName )
// Create dataset now
std::shared_ptr<DatasetGroup> ds = std::make_shared< DatasetGroup >(
name(),
mMesh.get(),
datFileName,
mesh,
timedepFileName,
grpName
);
ds->setIsOnVertices( false );
ds->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
ds->setIsScalar( !isVector );
for ( size_t ts = 0; ts < timesteps; ++ts )
{
std::shared_ptr< MemoryDataset > output = std::make_shared< MemoryDataset >( ds.get() );
std::shared_ptr< MemoryDataset2D > output = std::make_shared< MemoryDataset2D >( ds.get() );
output->setTime( times[ts] );
if ( isVector )
@ -588,7 +601,7 @@ bool MDAL::DriverFlo2D::parseHDF5Datasets( const std::string &datFileName )
// TODO use mins & maxs arrays
ds->setStatistics( MDAL::calculateStatistics( ds ) );
mMesh->datasetGroups.push_back( ds );
mesh->datasetGroups.push_back( ds );
}
@ -612,7 +625,7 @@ MDAL::DriverFlo2D::DriverFlo2D()
"FLO2D",
"Flo2D",
"*.nc",
Capability::ReadMesh )
Capability::ReadMesh | Capability::ReadDatasets | Capability::WriteDatasetsOnFaces2D )
{
}
@ -622,7 +635,7 @@ MDAL::DriverFlo2D *MDAL::DriverFlo2D::create()
return new DriverFlo2D();
}
bool MDAL::DriverFlo2D::canRead( const std::string &uri )
bool MDAL::DriverFlo2D::canReadMesh( const std::string &uri )
{
std::string cadptsFile( fileNameFromDir( uri, "CADPTS.DAT" ) );
if ( !MDAL::fileExists( cadptsFile ) )
@ -639,6 +652,44 @@ bool MDAL::DriverFlo2D::canRead( const std::string &uri )
return true;
}
bool MDAL::DriverFlo2D::canReadDatasets( const std::string &uri )
{
if ( !fileExists( uri ) ) return false;
HdfFile file( uri, HdfFile::ReadOnly );
if ( !file.isValid() ) return false;
HdfGroup timedataGroup = file.group( "TIMDEP NETCDF OUTPUT RESULTS" );
if ( !timedataGroup.isValid() ) return false;
return true;
}
void MDAL::DriverFlo2D::load( const std::string &uri, MDAL::Mesh *mesh, MDAL_Status *status )
{
if ( status ) *status = MDAL_Status::None;
MDAL::MemoryMesh *memoryMesh = dynamic_cast<MDAL::MemoryMesh *>( mesh );
if ( !memoryMesh )
{
if ( status ) *status = MDAL_Status::Err_IncompatibleMesh;
return;
}
if ( !MDAL::fileExists( uri ) )
{
if ( status ) *status = MDAL_Status::Err_FileNotFound;
return;
}
bool err = parseHDF5Datasets( memoryMesh, uri );
if ( err )
{
// TODO better error message?
if ( status ) *status = MDAL_Status::Err_InvalidData;
}
}
std::unique_ptr< MDAL::Mesh > MDAL::DriverFlo2D::load( const std::string &resultsFile, MDAL_Status *status )
{
mDatFileName = resultsFile;
@ -660,7 +711,9 @@ std::unique_ptr< MDAL::Mesh > MDAL::DriverFlo2D::load( const std::string &result
// create output for bed elevation
addStaticDataset( elevations, "Bed Elevation", mDatFileName );
if ( parseHDF5Datasets( mDatFileName ) )
// check if we have HDF5 file
std::string TIMDEPFileName = fileNameFromDir( mDatFileName, "TIMDEP.HDF5" );
if ( parseHDF5Datasets( mMesh.get(), TIMDEPFileName ) )
{
// some problem with HDF5 data, try text files
parseOUTDatasets( mDatFileName, elevations );
@ -675,3 +728,172 @@ std::unique_ptr< MDAL::Mesh > MDAL::DriverFlo2D::load( const std::string &result
return std::unique_ptr<Mesh>( mMesh.release() );
}
bool MDAL::DriverFlo2D::addToHDF5File( DatasetGroup *group )
{
assert( MDAL::fileExists( group->uri() ) );
HdfFile file( group->uri(), HdfFile::ReadWrite );
if ( !file.isValid() ) return true;
HdfGroup timedataGroup = file.group( "TIMDEP NETCDF OUTPUT RESULTS" );
if ( !timedataGroup.isValid() ) return true;
return appendGroup( file, group, timedataGroup );
}
bool MDAL::DriverFlo2D::saveNewHDF5File( DatasetGroup *dsGroup )
{
// Create file
HdfFile file( dsGroup->uri(), HdfFile::Create );
// Unable to create
if ( !file.isValid() ) return true;
// Create float dataset File Version
HdfDataset dsFileVersion( file.id(), "/File Version", H5T_NATIVE_FLOAT );
dsFileVersion.write( 1.0f );
// Create string dataset File Type
HdfDataset dsFileType( file.id(), "/File Type", HdfDataType::createString() );
dsFileType.write( "Xmdf" );
// Create group TIMDEP NETCDF OUTPUT RESULTS
HdfGroup groupTNOR = HdfGroup::create( file.id(), "/TIMDEP NETCDF OUTPUT RESULTS" );
// Create attribute
HdfAttribute attTNORGrouptype( groupTNOR.id(), "Grouptype", HdfDataType::createString() );
// Write string value to attribute
attTNORGrouptype.write( "Generic" );
return appendGroup( file, dsGroup, groupTNOR );
}
bool MDAL::DriverFlo2D::appendGroup( HdfFile &file, MDAL::DatasetGroup *dsGroup, HdfGroup &groupTNOR )
{
assert( dsGroup->dataLocation() == MDAL_DataLocation::DataOnFaces2D );
HdfDataType dtMaxString = HdfDataType::createString();
std::string dsGroupName = dsGroup->name();
const size_t timesCount = dsGroup->datasets.size();
const size_t facesCount = dsGroup->mesh()->facesCount();
size_t valCount = facesCount;
HdfDataspace dscValues;
if ( dsGroup->isScalar() )
{
std::vector<hsize_t> dimsForScalarValues = { timesCount, facesCount };
dscValues = HdfDataspace( dimsForScalarValues );
}
else
{
valCount *= 2;
std::vector<hsize_t> dimsForVectorValues = { timesCount, facesCount, 2 };
dscValues = HdfDataspace( dimsForVectorValues );
}
std::vector<hsize_t> timesCountVec = {timesCount};
HdfDataspace dscTimes( timesCountVec );
std::vector<float> maximums( timesCount );
std::vector<float> minimums( timesCount );
std::vector<double> times( timesCount );
std::vector<float> values( timesCount * valCount );
// prepare data
for ( size_t i = 0; i < dsGroup->datasets.size(); i++ )
{
const std::shared_ptr<Dataset> &dataset = dsGroup->datasets.at( i );
std::vector<double> singleRowValues( valCount );
if ( dsGroup->isScalar() )
{
dataset->scalarData( 0, facesCount, singleRowValues.data() );
}
else
{
dataset->vectorData( 0, facesCount, singleRowValues.data() );
}
for ( size_t j = 0; j < valCount; j++ )
{
double doubleValue = toFlo2DDouble( singleRowValues[j] );
values[i * valCount + j] = static_cast<float>( doubleValue );
}
const Statistics st = dataset->statistics();
maximums[i] = static_cast<float>( st.maximum );
minimums[i] = static_cast<float>( st.minimum );
times.push_back( dataset->time() );
}
// store data
int i = 0;
while ( file.pathExists( "/TIMDEP NETCDF OUTPUT RESULTS/" + dsGroupName ) )
{
dsGroupName = dsGroup->name() + "_" + std::to_string( i ); // make sure we have unique group name
}
HdfGroup group = HdfGroup::create( groupTNOR.id(), "/TIMDEP NETCDF OUTPUT RESULTS/" + dsGroupName );
HdfAttribute attDataType( group.id(), "Data Type", H5T_NATIVE_INT );
attDataType.write( 0 );
HdfAttribute attDatasetCompression( group.id(), "DatasetCompression", H5T_NATIVE_INT );
attDatasetCompression.write( -1 );
/*
HdfDataspace dscDatasetUnits( dimsSingle );
HdfAttribute attDatasetUnits( group.id(), "DatasetUnits", true );
attDatasetUnits.writeString( dscDatasetUnits.id(), "unknown" );
*/
HdfAttribute attGrouptype( group.id(), "Grouptype", dtMaxString );
if ( dsGroup->isScalar() )
attGrouptype.write( "DATASET SCALAR" );
else
attGrouptype.write( "DATASET VECTOR" );
HdfAttribute attTimeUnits( group.id(), "TimeUnits", dtMaxString );
attTimeUnits.write( "Hours" );
HdfDataset dsMaxs( file.id(), "/TIMDEP NETCDF OUTPUT RESULTS/" + dsGroupName + "/Maxs", H5T_NATIVE_FLOAT, timesCountVec );
dsMaxs.write( maximums );
HdfDataset dsMins( file.id(), "/TIMDEP NETCDF OUTPUT RESULTS/" + dsGroupName + "/Mins", H5T_NATIVE_FLOAT, timesCountVec );
dsMins.write( minimums );
HdfDataset dsTimes( file.id(), "/TIMDEP NETCDF OUTPUT RESULTS/" + dsGroupName + "/Times", H5T_NATIVE_DOUBLE, timesCountVec );
dsTimes.write( times );
HdfDataset dsValues( file.id(), "/TIMDEP NETCDF OUTPUT RESULTS/" + dsGroupName + "/Values", H5T_NATIVE_FLOAT, dscValues );
dsValues.write( values );
return false; //OK
}
bool MDAL::DriverFlo2D::persist( DatasetGroup *group )
{
if ( !group || ( group->dataLocation() != MDAL_DataLocation::DataOnFaces2D ) )
{
MDAL::debug( "flo-2d can store only 2D face datasets" );
return true;
}
try
{
// Return true on error
if ( MDAL::fileExists( group->uri() ) )
{
// Add dataset to a existing file
return addToHDF5File( group );
}
else
{
// Create new HDF5 file with Flow2D structure
return saveNewHDF5File( group );
}
}
catch ( MDAL_Status error )
{
MDAL::debug( "Error status: " + std::to_string( error ) );
return true;
}
}

17
external/mdal/frmts/mdal_flo2d.hpp vendored
View File

@ -13,6 +13,9 @@
#include "mdal.h"
#include "mdal_driver.hpp"
class HdfGroup;
class HdfFile;
namespace MDAL
{
class DriverFlo2D: public Driver
@ -22,8 +25,12 @@ namespace MDAL
~DriverFlo2D( ) override = default;
DriverFlo2D *create() override;
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
bool canReadDatasets( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override;
void load( const std::string &uri, Mesh *mesh, MDAL_Status *status ) override;
bool persist( DatasetGroup *group ) override;
private:
struct CellCenter
@ -39,7 +46,7 @@ namespace MDAL
void createMesh( const std::vector<CellCenter> &cells, double half_cell_size );
void parseOUTDatasets( const std::string &datFileName, const std::vector<double> &elevations );
bool parseHDF5Datasets( const std::string &datFileName );
bool parseHDF5Datasets( MDAL::MemoryMesh *mesh, const std::string &timedepFileName );
void parseVELFPVELOCFile( const std::string &datFileName );
void parseDEPTHFile( const std::string &datFileName, const std::vector<double> &elevations );
void parseTIMDEPFile( const std::string &datFileName, const std::vector<double> &elevations );
@ -48,6 +55,12 @@ namespace MDAL
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 );
// Write API
bool addToHDF5File( DatasetGroup *group );
bool saveNewHDF5File( DatasetGroup *group );
bool appendGroup( HdfFile &file, DatasetGroup *dsGroup, HdfGroup &groupTNOR );
};
} // namespace MDAL

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

@ -313,7 +313,7 @@ void MDAL::DriverGdal::fixRasterBands()
}
}
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<MemoryDataset2D> tos, bool is_vector, bool is_x )
{
assert( raster_band );
@ -410,14 +410,14 @@ void MDAL::DriverGdal::addDatasetGroups()
mFileName,
band->first
);
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
bool is_vector = ( band->second.begin()->second.size() > 1 );
group->setIsScalar( !is_vector );
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;
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MDAL::MemoryDataset2D >( group.get() );
dataset->setTime( time_step->first );
for ( std::vector<GDALRasterBandH>::size_type i = 0; i < raster_bands.size(); ++i )
@ -528,7 +528,7 @@ MDAL::DriverGdal::DriverGdal( const std::string &name,
mPafScanline( nullptr )
{}
bool MDAL::DriverGdal::canRead( const std::string &uri )
bool MDAL::DriverGdal::canReadMesh( const std::string &uri )
{
try
{

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

@ -57,7 +57,7 @@ namespace MDAL
const std::string &gdalDriverName );
virtual ~DriverGdal() override = default;
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &fileName, MDAL_Status *status ) override;
protected:
@ -88,7 +88,7 @@ namespace MDAL
bool meshes_equals( const GdalDataset *ds1, const GdalDataset *ds2 ) const;
metadata_hash parseMetadata( GDALMajorObjectH gdalBand, const char *pszDomain = nullptr );
void addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x );
void addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset2D> tos, bool is_vector, bool is_x );
bool addSrcProj();
void addDatasetGroups();
void createMesh();

243
external/mdal/frmts/mdal_hdf5.cpp vendored
View File

@ -5,20 +5,53 @@
#include "mdal_hdf5.hpp"
#include <cstring>
#include <algorithm>
HdfFile::HdfFile( const std::string &path )
: d( std::make_shared< Handle >( H5Fopen( path.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) ) )
HdfFile::HdfFile( const std::string &path, HdfFile::Mode mode )
: mPath( path )
{
switch ( mode )
{
case HdfFile::ReadOnly:
if ( H5Fis_hdf5( mPath.c_str() ) > 0 )
d = std::make_shared< Handle >( H5Fopen( path.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) );
break;
case HdfFile::ReadWrite:
if ( H5Fis_hdf5( mPath.c_str() ) > 0 )
d = std::make_shared< Handle >( H5Fopen( path.c_str(), H5F_ACC_RDWR, H5P_DEFAULT ) );
break;
case HdfFile::Create:
d = std::make_shared< Handle >( H5Fcreate( path.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT ) );
break;
}
}
bool HdfFile::isValid() const { return d->id >= 0; }
HdfFile::~HdfFile() = default;
bool HdfFile::isValid() const { return d && ( d->id >= 0 ); }
hid_t HdfFile::id() const { return d->id; }
std::string HdfFile::filePath() const
{
return mPath;
}
HdfGroup HdfGroup::create( hid_t file, const std::string &path )
{
auto d = std::make_shared< Handle >( H5Gcreate2( file, path.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT ) );
return HdfGroup( d );
}
HdfGroup::HdfGroup( hid_t file, const std::string &path )
: d( std::make_shared< Handle >( H5Gopen( file, path.c_str() ) ) )
{}
{
d = std::make_shared< Handle >( H5Gopen( file, path.c_str() ) );
}
HdfGroup::HdfGroup( std::shared_ptr<Handle> handle )
: d( handle )
{
}
bool HdfGroup::isValid() const { return d->id >= 0; }
@ -59,9 +92,21 @@ std::vector<std::string> HdfGroup::objects( H5G_obj_t type ) const
return lst;
}
HdfAttribute::HdfAttribute( hid_t obj_id, const std::string &attr_name, HdfDataType type )
: mType( type )
{
std::vector<hsize_t> dimsSingle = {1};
HdfDataspace dsc( dimsSingle );
d = std::make_shared< Handle >( H5Acreate2( obj_id, attr_name.c_str(), type.id(), dsc.id(), H5P_DEFAULT, H5P_DEFAULT ) );
}
HdfAttribute::HdfAttribute( hid_t obj_id, const std::string &attr_name )
: d( std::make_shared< Handle >( H5Aopen( obj_id, attr_name.c_str(), H5P_DEFAULT ) ) )
{}
: m_objId( obj_id ), m_name( attr_name )
{
d = std::make_shared< Handle >( H5Aopen( obj_id, attr_name.c_str(), H5P_DEFAULT ) );
}
HdfAttribute::~HdfAttribute() = default;
bool HdfAttribute::isValid() const { return d->id >= 0; }
@ -69,26 +114,65 @@ hid_t HdfAttribute::id() const { return d->id; }
std::string HdfAttribute::readString() const
{
hid_t datatype = H5Aget_type( id() );
HdfDataType datatype( H5Aget_type( id() ) );
char name[HDF_MAX_NAME + 1];
std::memset( name, '\0', HDF_MAX_NAME + 1 );
herr_t status = H5Aread( d->id, datatype, name );
herr_t status = H5Aread( d->id, datatype.id(), name );
if ( status < 0 )
{
//MDAL::debug("Failed to read data!");
H5Tclose( datatype );
return std::string();
}
H5Tclose( datatype );
std::string res( name );
res = MDAL::trim( res );
return res;
}
void HdfAttribute::write( const std::string &value )
{
if ( !isValid() || !mType.isValid() )
throw MDAL_Status::Err_FailToWriteToDisk;
// make sure you do not store more than it is possible
std::vector<char> buf( HDF_MAX_NAME + 1, '\0' );
size_t size = value.size() < HDF_MAX_NAME ? value.size() : HDF_MAX_NAME;
memcpy( buf.data(), value.c_str(), size );
if ( H5Awrite( d->id, mType.id(), buf.data() ) < 0 )
throw MDAL_Status::Err_FailToWriteToDisk;
}
void HdfAttribute::write( int value )
{
if ( !isValid() || !mType.isValid() )
throw MDAL_Status::Err_FailToWriteToDisk;
if ( H5Awrite( d->id, mType.id(), &value ) < 0 )
throw MDAL_Status::Err_FailToWriteToDisk;
}
HdfDataset::HdfDataset( hid_t file, const std::string &path, HdfDataType dtype, size_t nItems )
: mType( dtype )
{
// Crete dataspace for attribute
std::vector<hsize_t> dimsSingle = {nItems};
HdfDataspace dsc( dimsSingle );
d = std::make_shared< Handle >( H5Dcreate2( file, path.c_str(), dtype.id(), dsc.id(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT ) );
}
HdfDataset::HdfDataset( hid_t file, const std::string &path, HdfDataType dtype, HdfDataspace dataspace )
: mType( dtype )
{
d = std::make_shared< Handle >( H5Dcreate2( file, path.c_str(), dtype.id(), dataspace.id(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT ) );
}
HdfDataset::HdfDataset( hid_t file, const std::string &path )
: d( std::make_shared< Handle >( H5Dopen2( file, path.c_str(), H5P_DEFAULT ) ) )
{}
{
}
HdfDataset::~HdfDataset() = default;
bool HdfDataset::isValid() const { return d->id >= 0; }
@ -113,10 +197,13 @@ hsize_t HdfDataset::elementCount() const
H5T_class_t HdfDataset::type() const
{
hid_t tid = H5Dget_type( d->id );
H5T_class_t t_class = H5Tget_class( tid );
H5Tclose( tid );
return t_class;
if ( mType.isValid() )
return H5Tget_class( mType.id() );
else
{
HdfDataType dt( H5Dget_type( d->id ) );
return H5Tget_class( dt.id() );
}
}
std::vector<uchar> HdfDataset::readArrayUint8( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const { return readArray<uchar>( H5T_NATIVE_UINT8, offsets, counts ); }
@ -139,12 +226,8 @@ std::vector<std::string> HdfDataset::readArrayString() const
{
std::vector<std::string> ret;
hid_t datatype = H5Tcopy( H5T_C_S1 );
H5Tset_size( datatype, HDF_MAX_NAME );
std::vector<HdfString> arr = readArray<HdfString>( datatype );
H5Tclose( datatype );
HdfDataType datatype = HdfDataType::createString();
std::vector<HdfString> arr = readArray<HdfString>( datatype.id() );
for ( const HdfString &str : arr )
{
@ -173,6 +256,51 @@ float HdfDataset::readFloat() const
return value;
}
void HdfDataset::write( std::vector<float> &value )
{
if ( !isValid() || !mType.isValid() )
throw MDAL_Status::Err_FailToWriteToDisk;
// Write float array to dataset
if ( H5Dwrite( d->id, mType.id(), H5S_ALL, H5S_ALL, H5P_DEFAULT, value.data() ) < 0 )
throw MDAL_Status::Err_FailToWriteToDisk;
}
void HdfDataset::write( float value )
{
if ( !isValid() || !mType.isValid() )
throw MDAL_Status::Err_FailToWriteToDisk;
// Write float array to dataset
if ( H5Dwrite( d->id, mType.id(), H5S_ALL, H5S_ALL, H5P_DEFAULT, &value ) < 0 )
throw MDAL_Status::Err_FailToWriteToDisk;
}
void HdfDataset::write( std::vector<double> &value )
{
if ( !isValid() || !mType.isValid() )
throw MDAL_Status::Err_FailToWriteToDisk;
// Write double array to dataset.
if ( H5Dwrite( d->id, mType.id(), H5S_ALL, H5S_ALL, H5P_DEFAULT, value.data() ) < 0 )
throw MDAL_Status::Err_FailToWriteToDisk;
}
void HdfDataset::write( const std::string &value )
{
if ( !isValid() || !mType.isValid() )
throw MDAL_Status::Err_FailToWriteToDisk;
// make sure you do not store more than it is possible
std::vector<char> buf( HDF_MAX_NAME + 1, '\0' );
size_t size = value.size() < HDF_MAX_NAME ? value.size() : HDF_MAX_NAME;
memcpy( buf.data(), value.c_str(), size );
// Write string to dataset.
if ( H5Dwrite( d->id, mType.id(), H5S_ALL, H5S_ALL, H5P_DEFAULT, buf.data() ) < 0 )
throw MDAL_Status::Err_FailToWriteToDisk;
}
std::string HdfDataset::readString() const
{
if ( elementCount() != 1 )
@ -182,34 +310,34 @@ std::string HdfDataset::readString() const
}
char name[HDF_MAX_NAME];
hid_t datatype = H5Tcopy( H5T_C_S1 );
H5Tset_size( datatype, HDF_MAX_NAME );
herr_t status = H5Dread( d->id, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, name );
HdfDataType datatype = HdfDataType::createString();
herr_t status = H5Dread( d->id, datatype.id(), H5S_ALL, H5S_ALL, H5P_DEFAULT, name );
if ( status < 0 )
{
MDAL::debug( "Failed to read data!" );
return std::string();
}
H5Tclose( datatype );
return std::string( name );
}
HdfDataspace::HdfDataspace( const std::vector<hsize_t> &dims )
: d( std::make_shared< Handle >( H5Screate_simple(
static_cast<int>( dims.size() ),
dims.data(),
dims.data()
)
)
)
{
d = std::make_shared< Handle >( H5Screate_simple(
static_cast<int>( dims.size() ),
dims.data(),
dims.data()
) );
}
HdfDataspace::HdfDataspace( hid_t dataset )
: d( std::make_shared< Handle >( H5Dget_space( dataset ) ) )
{
if ( dataset >= 0 )
d = std::make_shared< Handle >( H5Dget_space( dataset ) );
}
HdfDataspace::~HdfDataspace() = default;
void HdfDataspace::selectHyperslab( hsize_t start, hsize_t count )
{
// this function works only for 1D arrays
@ -243,3 +371,44 @@ void HdfDataspace::selectHyperslab( const std::vector<hsize_t> offsets,
bool HdfDataspace::isValid() const { return d->id >= 0; }
hid_t HdfDataspace::id() const { return d->id; }
HdfDataType::HdfDataType() = default;
HdfDataType::HdfDataType( hid_t type, bool isNativeType )
{
if ( isNativeType )
mNativeId = type;
else
d = std::make_shared< Handle >( type );
}
HdfDataType HdfDataType::createString( int size )
{
assert( size > 0 );
if ( size > HDF_MAX_NAME )
size = HDF_MAX_NAME;
hid_t atype = H5Tcopy( H5T_C_S1 );
H5Tset_size( atype, static_cast<size_t>( size ) );
H5Tset_strpad( atype, H5T_STR_NULLTERM );
return HdfDataType( atype, false );
}
HdfDataType::~HdfDataType() = default;
bool HdfDataType::isValid() const
{
if ( d )
return d->id >= 0;
else
return mNativeId >= 0;
}
hid_t HdfDataType::id() const
{
if ( d )
return d->id;
else
return mNativeId;
}

78
external/mdal/frmts/mdal_hdf5.hpp vendored
View File

@ -36,6 +36,7 @@ template <> inline void hdfClose<H5I_GROUP>( hid_t id ) { H5Gclose( id ); }
template <> inline void hdfClose<H5I_DATASET>( hid_t id ) { H5Dclose( id ); }
template <> inline void hdfClose<H5I_ATTR>( hid_t id ) { H5Aclose( id ); }
template <> inline void hdfClose<H5I_DATASPACE>( hid_t id ) { H5Sclose( id ); }
template <> inline void hdfClose<H5I_DATATYPE>( hid_t id ) { H5Tclose( id ); }
template <int TYPE>
class HdfH
@ -52,14 +53,21 @@ class HdfGroup;
class HdfDataset;
class HdfAttribute;
class HdfDataspace;
class HdfDataType;
class HdfFile
{
public:
enum Mode
{
ReadOnly,
ReadWrite,
Create
};
typedef HdfH<H5I_FILE> Handle;
HdfFile( const std::string &path );
HdfFile( const std::string &path, HdfFile::Mode mode );
~HdfFile();
bool isValid() const;
hid_t id() const;
@ -69,9 +77,30 @@ class HdfFile
inline HdfDataset dataset( const std::string &path ) const;
inline HdfAttribute attribute( const std::string &attr_name ) const;
inline bool pathExists( const std::string &path ) const;
std::string filePath() const;
protected:
std::shared_ptr<Handle> d;
std::string mPath;
};
class HdfDataType
{
public:
typedef HdfH<H5I_DATATYPE> Handle;
HdfDataType();
HdfDataType( hid_t type, bool isNativeType = true );
~HdfDataType();
// Creates new string type with size, use HDF_MAX_NAME for maximum length
static HdfDataType createString( int size = HDF_MAX_NAME );
bool isValid() const;
hid_t id() const;
protected:
std::shared_ptr<Handle> d;
hid_t mNativeId = -1;
};
class HdfGroup
@ -79,7 +108,9 @@ class HdfGroup
public:
typedef HdfH<H5I_GROUP> Handle;
static HdfGroup create( hid_t file, const std::string &path );
HdfGroup( hid_t file, const std::string &path );
HdfGroup( std::shared_ptr<Handle> handle );
bool isValid() const;
hid_t id() const;
@ -104,30 +135,44 @@ class HdfGroup
std::shared_ptr<Handle> d;
};
class HdfAttribute
{
public:
typedef HdfH<H5I_ATTR> Handle;
HdfAttribute( hid_t obj_id, const std::string &attr_name );
//! Create new attribute for writing for 1 item
HdfAttribute( hid_t obj_id, const std::string &attr_name, HdfDataType type );
//! Open existing attribute for reading
HdfAttribute( hid_t obj_id, const std::string &attr_name );
~HdfAttribute();
bool isValid() const;
hid_t id() const;
std::string readString() const;
void write( const std::string &value );
void write( int value );
protected:
std::shared_ptr<Handle> d;
hid_t m_objId;
std::string m_name;
HdfDataType mType; // when in write mode
};
class HdfDataspace
{
public:
typedef HdfH<H5I_DATASPACE> Handle;
//! memory dataspace for simple N-D array
HdfDataspace( const std::vector<hsize_t> &dims );
//! dataspace of the dataset
HdfDataspace( hid_t dataset );
HdfDataspace( hid_t dataset = -1 );
~HdfDataspace( );
//! select from 1D array
void selectHyperslab( hsize_t start, hsize_t count );
//! select from N-D array
@ -146,8 +191,13 @@ class HdfDataset
public:
typedef HdfH<H5I_DATASET> Handle;
//! Create new, simple 1 dimensional dataset
HdfDataset( hid_t file, const std::string &path, HdfDataType dtype, size_t nItems = 1 );
//! Create new dataset with custom dimensions
HdfDataset( hid_t file, const std::string &path, HdfDataType dtype, HdfDataspace dataspace );
//! Opens dataset for reading
HdfDataset( hid_t file, const std::string &path );
~HdfDataset();
bool isValid() const;
hid_t id() const;
@ -216,12 +266,28 @@ class HdfDataset
return data;
}
//! Reads float value
float readFloat() const;
//! Reads string value
std::string readString() const;
//! Writes string dataset with single entry
void write( const std::string &value );
//! Writes array of float data
void write( std::vector<float> &value );
void write( float value );
//! Writes array of double data
void write( std::vector<double> &value );
protected:
std::shared_ptr<Handle> d;
hid_t m_fileId;
std::string m_path;
HdfDataType mType; // when in write mode
};
inline std::vector<std::string> HdfFile::groups() const { return group( "/" ).groups(); }

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

@ -17,7 +17,7 @@
static HdfFile openHdfFile( const std::string &fileName )
{
HdfFile file( fileName );
HdfFile file( fileName, HdfFile::ReadOnly );
if ( !file.isValid() ) throw MDAL_Status::Err_UnknownFormat;
return file;
}
@ -168,21 +168,21 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
mFileName,
datasetName
);
group->setIsOnVertices( false );
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
group->setIsScalar( true );
group->setReferenceTime( referenceTime );
std::vector<std::shared_ptr<MDAL::MemoryDataset>> datasets;
std::vector<std::shared_ptr<MDAL::MemoryDataset2D>> datasets;
for ( size_t tidx = 0; tidx < times.size(); ++tidx )
{
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
double time = static_cast<double>( times[tidx] );
dataset->setTime( time );
datasets.push_back( dataset );
}
std::shared_ptr<MDAL::MemoryDataset> firstDataset;
std::shared_ptr<MDAL::MemoryDataset2D> firstDataset;
for ( size_t nArea = 0; nArea < flowAreaNames.size(); ++nArea )
{
@ -197,7 +197,7 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
for ( size_t tidx = 0; tidx < times.size(); ++tidx )
{
std::shared_ptr<MDAL::MemoryDataset> dataset = datasets[tidx];
std::shared_ptr<MDAL::MemoryDataset2D> dataset = datasets[tidx];
double *values = dataset->values();
for ( size_t i = 0; i < nFaces; ++i )
@ -251,13 +251,13 @@ void MDAL::DriverHec2D::readFaceResults( const HdfFile &hdfFile,
}
std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readElemOutput( const HdfGroup &rootGroup,
std::shared_ptr<MDAL::MemoryDataset2D> MDAL::DriverHec2D::readElemOutput( const HdfGroup &rootGroup,
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames,
const std::string rawDatasetName,
const std::string datasetName,
const std::vector<float> &times,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation,
std::shared_ptr<MDAL::MemoryDataset2D> bed_elevation,
const std::string &referenceTime )
{
double eps = std::numeric_limits<double>::min();
@ -268,15 +268,15 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readElemOutput( const Hd
mFileName,
datasetName
);
group->setIsOnVertices( false );
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
group->setIsScalar( true );
group->setReferenceTime( referenceTime );
std::vector<std::shared_ptr<MDAL::MemoryDataset>> datasets;
std::vector<std::shared_ptr<MDAL::MemoryDataset2D>> datasets;
for ( size_t tidx = 0; tidx < times.size(); ++tidx )
{
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
double time = static_cast<double>( times[tidx] );
dataset->setTime( time );
datasets.push_back( dataset );
@ -293,7 +293,7 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readElemOutput( const Hd
for ( size_t tidx = 0; tidx < times.size(); ++tidx )
{
std::shared_ptr<MDAL::MemoryDataset> dataset = datasets[tidx];
std::shared_ptr<MDAL::MemoryDataset2D> dataset = datasets[tidx];
double *values = dataset->values();
for ( size_t i = 0; i < nAreaElements; ++i )
@ -343,7 +343,7 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readElemOutput( const Hd
return datasets[0];
}
std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readBedElevation(
std::shared_ptr<MDAL::MemoryDataset2D> MDAL::DriverHec2D::readBedElevation(
const HdfGroup &gGeom2DFlowAreas,
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames )
@ -358,14 +358,14 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readBedElevation(
"Cells Minimum Elevation",
"Bed Elevation",
times,
std::shared_ptr<MDAL::MemoryDataset>(),
std::shared_ptr<MDAL::MemoryDataset2D>(),
referenceTime
);
}
void MDAL::DriverHec2D::readElemResults(
const HdfFile &hdfFile,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation,
std::shared_ptr<MDAL::MemoryDataset2D> bed_elevation,
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames )
{
@ -596,7 +596,7 @@ MDAL::DriverHec2D *MDAL::DriverHec2D::create()
return new DriverHec2D();
}
bool MDAL::DriverHec2D::canRead( const std::string &uri )
bool MDAL::DriverHec2D::canReadMesh( const std::string &uri )
{
try
{
@ -649,7 +649,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverHec2D::load( const std::string &resultsF
setProjection( hdfFile );
//Elevation
std::shared_ptr<MDAL::MemoryDataset> bed_elevation = readBedElevation( gGeom2DFlowAreas, areaElemStartIndex, flowAreaNames );
std::shared_ptr<MDAL::MemoryDataset2D> bed_elevation = readBedElevation( gGeom2DFlowAreas, areaElemStartIndex, flowAreaNames );
// Element centered Values
readElemResults( hdfFile, bed_elevation, areaElemStartIndex, flowAreaNames );

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

@ -42,7 +42,7 @@ namespace MDAL
~DriverHec2D( ) override = default;
DriverHec2D *create() override;
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override;
private:
@ -71,17 +71,17 @@ namespace MDAL
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames );
std::shared_ptr<MDAL::MemoryDataset> readElemOutput(
std::shared_ptr<MDAL::MemoryDataset2D> readElemOutput(
const HdfGroup &rootGroup,
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames,
const std::string rawDatasetName,
const std::string datasetName,
const std::vector<float> &times,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation,
std::shared_ptr<MDAL::MemoryDataset2D> bed_elevation,
const std::string &referenceTime );
std::shared_ptr<MDAL::MemoryDataset> readBedElevation(
std::shared_ptr<MDAL::MemoryDataset2D> readBedElevation(
const HdfGroup &gGeom2DFlowAreas,
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames );
@ -94,7 +94,7 @@ namespace MDAL
void readElemResults(
const HdfFile &hdfFile,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation,
std::shared_ptr<MDAL::MemoryDataset2D> bed_elevation,
const std::vector<size_t> &areaElemStartIndex,
const std::vector<std::string> &flowAreaNames );
};

218
external/mdal/frmts/mdal_netcdf.cpp vendored
View File

@ -7,6 +7,8 @@
#include <vector>
#include <assert.h>
#include <netcdf.h>
#include <cmath>
#include "mdal_netcdf.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
@ -15,7 +17,11 @@ NetCDFFile::NetCDFFile(): mNcid( 0 ) {}
NetCDFFile::~NetCDFFile()
{
if ( mNcid != 0 ) nc_close( mNcid );
if ( mNcid != 0 )
{
nc_close( mNcid );
mNcid = 0;
}
}
int NetCDFFile::handle() const
@ -43,6 +49,34 @@ std::vector<int> NetCDFFile::readIntArr( const std::string &name, size_t dim ) c
return arr_val;
}
std::vector<int> NetCDFFile::readIntArr( int arr_id, size_t start_dim1, size_t start_dim2, size_t count_dim1, size_t count_dim2 ) const
{
assert( mNcid != 0 );
const std::vector<size_t> startp = {start_dim1, start_dim2};
const std::vector<size_t> countp = {count_dim1, count_dim2};
const std::vector<ptrdiff_t> stridep = {1, 1};
std::vector<int> arr_val( count_dim1 * count_dim2 );
int res = nc_get_vars_int( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val.data() );
if ( res != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
return arr_val;
}
std::vector<int> NetCDFFile::readIntArr( int arr_id, size_t start_dim, size_t count_dim ) const
{
assert( mNcid != 0 );
const std::vector<size_t> startp = {start_dim};
const std::vector<size_t> countp = {count_dim};
const std::vector<ptrdiff_t> stridep = {1};
std::vector<int> arr_val( count_dim );
int res = nc_get_vars_int( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val.data() );
if ( res != 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 );
@ -50,7 +84,84 @@ std::vector<double> NetCDFFile::readDoubleArr( const std::string &name, size_t d
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;
nc_type typep;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_vartype( mNcid, arr_id, &typep ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( typep == NC_FLOAT )
{
std::vector<float> arr_val_f( dim );
if ( nc_get_var_float( mNcid, arr_id, arr_val_f.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
for ( size_t i = 0; i < dim; ++i )
{
const float val = arr_val_f[i];
if ( std::isnan( val ) )
arr_val[i] = std::numeric_limits<double>::quiet_NaN();
else
arr_val[i] = static_cast<double>( val );
}
}
else if ( typep == NC_DOUBLE )
{
if ( nc_get_var_double( mNcid, arr_id, arr_val.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
else
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}
std::vector<double> NetCDFFile::readDoubleArr( int arr_id,
size_t start_dim1, size_t start_dim2,
size_t count_dim1, size_t count_dim2 ) const
{
assert( mNcid != 0 );
const std::vector<size_t> startp = {start_dim1, start_dim2};
const std::vector<size_t> countp = {count_dim1, count_dim2};
const std::vector<ptrdiff_t> stridep = {1, 1};
std::vector<double> arr_val( count_dim1 * count_dim2 );
nc_type typep;
if ( nc_inq_vartype( mNcid, arr_id, &typep ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( typep == NC_FLOAT )
{
std::vector<float> arr_val_f( count_dim1 * count_dim2 );
if ( nc_get_vars_float( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val_f.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
for ( size_t i = 0; i < count_dim1 * count_dim2; ++i )
{
const float val = arr_val_f[i];
if ( std::isnan( val ) )
arr_val[i] = std::numeric_limits<double>::quiet_NaN();
else
arr_val[i] = static_cast<double>( val );
}
}
else if ( typep == NC_BYTE )
{
std::vector<unsigned char> arr_val_b( count_dim1 * count_dim2 );
if ( nc_get_vars_ubyte( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val_b.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
for ( size_t i = 0; i < count_dim1 * count_dim2; ++i )
{
const unsigned char val = arr_val_b[i];
if ( val == 129 )
arr_val[i] = std::numeric_limits<double>::quiet_NaN();
else
arr_val[i] = double( int( val ) );
}
}
else if ( typep == NC_DOUBLE )
{
if ( nc_get_vars_double( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
else
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}
@ -61,6 +172,16 @@ bool NetCDFFile::hasArr( const std::string &name ) const
return nc_inq_varid( mNcid, name.c_str(), &arr_id ) == NC_NOERR;
}
int NetCDFFile::arrId( const std::string &name ) const
{
int arr_id = -1;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR )
{
arr_id = -1;
}
return arr_id;
}
std::vector<std::string> NetCDFFile::readArrNames() const
{
assert( mNcid != 0 );
@ -186,7 +307,6 @@ void NetCDFFile::getDimensions( const std::string &variableName, std::vector<siz
{
nc_inq_dimlen( mNcid, dimensionIds[size_t( i )], &dimensions[size_t( i )] );
}
}
bool NetCDFFile::hasDimension( const std::string &name ) const
@ -194,3 +314,95 @@ bool NetCDFFile::hasDimension( const std::string &name ) const
int ncid_val;
return nc_inq_dimid( mNcid, name.c_str(), &ncid_val ) == NC_NOERR;
}
void NetCDFFile::createFile( const std::string &fileName )
{
int res = nc_create( fileName.c_str(), NC_CLOBBER, &mNcid );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
}
int NetCDFFile::defineDimension( const std::string &name, size_t size )
{
int dimId = 0;
int res = nc_def_dim( mNcid, name.c_str(), size, &dimId );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
return dimId;
}
int NetCDFFile::defineVar( const std::string &varName,
int ncType, int dimensionCount, const int *dimensions )
{
int varIdp;
int res = nc_def_var( mNcid, varName.c_str(), ncType, dimensionCount, dimensions, &varIdp );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
return varIdp;
}
void NetCDFFile::putAttrStr( int varId, const std::string &attrName, const std::string &value )
{
int res = nc_put_att_text( mNcid, varId, attrName.c_str(), value.size(), value.c_str() );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
}
void NetCDFFile::putAttrInt( int varId, const std::string &attrName, int value )
{
int res = nc_put_att_int( mNcid, varId, attrName.c_str(), NC_INT, 1, &value );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
}
void NetCDFFile::putAttrDouble( int varId, const std::string &attrName, double value )
{
int res = nc_put_att_double( mNcid, varId, attrName.c_str(), NC_DOUBLE, 1, &value );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
}
void NetCDFFile::putDataDouble( int varId, const size_t index, const double value )
{
int res = nc_put_var1_double( mNcid, varId, &index, &value );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
}
void NetCDFFile::putDataArrayInt( int varId, size_t line, size_t faceVerticesMax, int *values )
{
// Configuration of these two vectors determines how is value array read and stored in the file
// https://www.unidata.ucar.edu/software/netcdf/docs/programming_notes.html#specify_hyperslabfileNameToSave
const size_t start[] = { line, 0 };
const size_t count[] = { 1, faceVerticesMax };
int res = nc_put_vara_int( mNcid, varId, start, count, values );
if ( res != NC_NOERR )
{
MDAL::debug( nc_strerror( res ) );
throw MDAL_Status::Err_FailToWriteToDisk;
}
}

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

@ -22,8 +22,33 @@ class NetCDFFile
void openFile( const std::string &fileName );
std::vector<int> readIntArr( const std::string &name, size_t dim ) const;
/** Reads hyperslap from double variable with 2 dimensions*/
std::vector<int> readIntArr( int arr_id,
size_t start_dim1,
size_t start_dim2,
size_t count_dim1,
size_t count_dim2
) const;
/** Reads hyperslap from double variable with 1 dimension */
std::vector<int> readIntArr( int arr_id,
size_t start_dim,
size_t count_dim
) const;
std::vector<double> readDoubleArr( const std::string &name, size_t dim ) const;
/** Reads hyperslap from double variable */
std::vector<double> readDoubleArr( int arr_id,
size_t start_dim1,
size_t start_dim2,
size_t count_dim1,
size_t count_dim2
) const;
bool hasArr( const std::string &name ) const;
int arrId( const std::string &name ) const;
std::vector<std::string> readArrNames() const;
bool hasAttrInt( const std::string &name, const std::string &attr_name ) const;
@ -37,12 +62,22 @@ class NetCDFFile
*/
std::string getAttrStr( const std::string &name, const std::string &attr_name ) const;
std::string getAttrStr( const std::string &attr_name, int varid ) const;
double getFillValue( int varid ) const;
int getVarId( const std::string &name );
void getDimension( const std::string &name, size_t *val, int *ncid_val ) const;
void getDimensions( const std::string &variableName, std::vector<size_t> &dimensionsId, std::vector<int> &dimensionIds );
bool hasDimension( const std::string &name ) const;
void createFile( const std::string &fileName );
int defineDimension( const std::string &name, size_t size );
int defineVar( const std::string &varName, int ncType, int dimensionCount, const int *dimensions );
void putAttrStr( int varId, const std::string &attrName, const std::string &value );
void putAttrInt( int varId, const std::string &attrName, int value );
void putAttrDouble( int varId, const std::string &attrName, double value );
void putDataDouble( int varId, const size_t index, const double value );
void putDataArrayInt( int varId, size_t line, size_t faceVerticesMax, int *values );
private:
int mNcid; // C handle to the file
};

14
external/mdal/frmts/mdal_selafin.cpp vendored
View File

@ -321,7 +321,7 @@ void MDAL::DriverSelafin::parseFile( std::vector<std::string> &var_names,
gives the numbering of boundary points for the others
*/
std::vector<int> iPointBoundary = mReader.read_int_arr( *nPoint );
MDAL_UNUSED( iPointBoundary );
MDAL_UNUSED( iPointBoundary )
/* 1 record containing table X (real array of dimension NPOIN containing the
abscissae of the points)
@ -459,7 +459,7 @@ void MDAL::DriverSelafin::addData( const std::vector<std::string> &var_names, co
mMesh->uri(),
var_name
);
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
group->setIsScalar( !is_vector );
mMesh->datasetGroups.push_back( group );
@ -468,14 +468,14 @@ void MDAL::DriverSelafin::addData( const std::vector<std::string> &var_names, co
size_t i = 0;
for ( timestep_map::const_iterator it = data[nName].begin(); it != data[nName].end(); ++it, ++i )
{
std::shared_ptr<MDAL::MemoryDataset> dataset;
std::shared_ptr<MDAL::MemoryDataset2D> dataset;
if ( group->datasets.size() > i )
{
dataset = std::dynamic_pointer_cast<MDAL::MemoryDataset>( group->datasets[i] );
dataset = std::dynamic_pointer_cast<MDAL::MemoryDataset2D>( group->datasets[i] );
}
else
{
dataset = std::make_shared< MemoryDataset >( group.get() );
dataset = std::make_shared< MemoryDataset2D >( group.get() );
dataset->setTime( it->first );
group->datasets.push_back( dataset );
}
@ -512,7 +512,7 @@ void MDAL::DriverSelafin::addData( const std::vector<std::string> &var_names, co
{
for ( auto dataset : group->datasets )
{
std::shared_ptr<MDAL::MemoryDataset> dts = std::dynamic_pointer_cast<MDAL::MemoryDataset>( dataset );
std::shared_ptr<MDAL::MemoryDataset2D> dts = std::dynamic_pointer_cast<MDAL::MemoryDataset2D>( dataset );
MDAL::activateFaces( mMesh.get(), dts );
MDAL::Statistics stats = MDAL::calculateStatistics( dataset );
@ -524,7 +524,7 @@ void MDAL::DriverSelafin::addData( const std::vector<std::string> &var_names, co
}
}
bool MDAL::DriverSelafin::canRead( const std::string &uri )
bool MDAL::DriverSelafin::canReadMesh( const std::string &uri )
{
if ( !MDAL::fileExists( uri ) ) return false;

2
external/mdal/frmts/mdal_selafin.hpp vendored
View File

@ -62,7 +62,7 @@ namespace MDAL
~DriverSelafin() override;
DriverSelafin *create() override;
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile, MDAL_Status *status ) override;
private:

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

@ -35,7 +35,7 @@ size_t MDAL::DriverSWW::getVertexCount( const NetCDFFile &ncFile ) const
}
bool MDAL::DriverSWW::canRead( const std::string &uri )
bool MDAL::DriverSWW::canReadMesh( const std::string &uri )
{
NetCDFFile ncFile;
@ -288,7 +288,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readScalarGroup(
mesh,
mFileName,
groupName );
mds->setIsOnVertices( true );
mds->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
mds->setIsScalar( true );
int zDimsX = 0;
@ -298,7 +298,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readScalarGroup(
if ( zDimsX == 1 )
{
// TIME INDEPENDENT
std::shared_ptr<MDAL::MemoryDataset> o = std::make_shared<MDAL::MemoryDataset>( mds.get() );
std::shared_ptr<MDAL::MemoryDataset2D> o = std::make_shared<MDAL::MemoryDataset2D>( mds.get() );
o->setTime( 0.0 );
double *values = o->values();
std::vector<double> valuesX = ncFile.readDoubleArr( arrName, nPoints );
@ -314,7 +314,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readScalarGroup(
// TIME DEPENDENT
for ( size_t t = 0; t < times.size(); ++t )
{
std::shared_ptr<MDAL::MemoryDataset> mto = std::make_shared<MDAL::MemoryDataset>( mds.get() );
std::shared_ptr<MDAL::MemoryDataset2D> mto = std::make_shared<MDAL::MemoryDataset2D>( mds.get() );
mto->setTime( static_cast<double>( times[t] ) / 3600. );
double *values = mto->values();
@ -357,7 +357,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readVectorGroup(
mesh,
mFileName,
groupName );
mds->setIsOnVertices( true );
mds->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
mds->setIsScalar( false );
int zDimsX = 0;
@ -376,7 +376,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readVectorGroup(
if ( zDimsX == 1 )
{
// TIME INDEPENDENT
std::shared_ptr<MDAL::MemoryDataset> o = std::make_shared<MDAL::MemoryDataset>( mds.get() );
std::shared_ptr<MDAL::MemoryDataset2D> o = std::make_shared<MDAL::MemoryDataset2D>( mds.get() );
o->setTime( 0.0 );
double *values = o->values();
std::vector<double> valuesX = ncFile.readDoubleArr( arrXName, nPoints );
@ -394,7 +394,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverSWW::readVectorGroup(
// TIME DEPENDENT
for ( size_t t = 0; t < times.size(); ++t )
{
std::shared_ptr<MDAL::MemoryDataset> mto = std::make_shared<MDAL::MemoryDataset>( mds.get() );
std::shared_ptr<MDAL::MemoryDataset2D> mto = std::make_shared<MDAL::MemoryDataset2D>( mds.get() );
mto->setTime( static_cast<double>( times[t] ) / 3600. );
double *values = mto->values();

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

@ -41,7 +41,7 @@ namespace MDAL
DriverSWW *create() override;
std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override;
bool canRead( const std::string &uri ) override;
bool canReadMesh( const std::string &uri ) override;
private:
size_t getVertexCount( const NetCDFFile &ncFile ) const;

397
external/mdal/frmts/mdal_tuflowfv.cpp vendored Normal file
View File

@ -0,0 +1,397 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#include "mdal_tuflowfv.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
#include "mdal_netcdf.hpp"
#include <math.h>
#include <assert.h>
#include <cstring>
MDAL::TuflowFVDataset3D::TuflowFVDataset3D(
MDAL::DatasetGroup *parent,
int ncid_x,
int ncid_y,
size_t timesteps,
size_t volumesCount,
size_t facesCount,
size_t levelFacesCount,
size_t ts,
size_t maximumLevelsCount,
std::shared_ptr<NetCDFFile> ncFile )
: MDAL::Dataset3D( parent, volumesCount, maximumLevelsCount )
, mNcidX( ncid_x )
, mNcidY( ncid_y )
, mTimesteps( timesteps )
, mFacesCount( facesCount )
, mLevelFacesCount( levelFacesCount )
, mTs( ts )
, mNcFile( ncFile )
{
if ( ncFile )
{
mNcidVerticalLevels = ncFile->arrId( "NL" );
mNcidVerticalLevelsZ = ncFile->arrId( "layerface_Z" );
mNcidActive2D = ncFile->arrId( "stat" );
mNcid3DTo2D = ncFile->arrId( "idx2" );
mNcid2DTo3D = ncFile->arrId( "idx3" );
}
}
MDAL::TuflowFVDataset3D::~TuflowFVDataset3D() = default;
size_t MDAL::TuflowFVDataset3D::verticalLevelCountData( size_t indexStart, size_t count, int *buffer )
{
if ( ( count < 1 ) || ( indexStart >= mFacesCount ) )
return 0;
if ( mNcidVerticalLevels < 0 )
return 0;
size_t copyValues = std::min( mFacesCount - indexStart, count );
std::vector<int> vals = mNcFile->readIntArr(
mNcidVerticalLevels,
indexStart,
copyValues
);
memcpy( buffer, vals.data(), copyValues * sizeof( int ) );
return copyValues;
}
size_t MDAL::TuflowFVDataset3D::verticalLevelData( size_t indexStart, size_t count, double *buffer )
{
if ( ( count < 1 ) || ( indexStart >= mLevelFacesCount ) )
return 0;
if ( mTs >= mTimesteps )
return 0;
if ( mNcidVerticalLevelsZ < 0 )
return 0;
size_t copyValues = std::min( mLevelFacesCount - indexStart, count );
std::vector<double> vals = mNcFile->readDoubleArr(
mNcidVerticalLevelsZ,
mTs,
indexStart,
1,
copyValues
);
memcpy( buffer, vals.data(), copyValues * sizeof( double ) );
return copyValues;
}
size_t MDAL::TuflowFVDataset3D::faceToVolumeData( size_t indexStart, size_t count, int *buffer )
{
if ( ( count < 1 ) || ( indexStart >= mFacesCount ) )
return 0;
if ( mNcid2DTo3D < 0 )
return 0;
size_t copyValues = std::min( mFacesCount - indexStart, count );
std::vector<int> vals = mNcFile->readIntArr(
mNcid2DTo3D,
indexStart,
copyValues
);
// indexed from 1 in FV, from 0 in MDAL
for ( auto &element : vals )
element -= 1;
memcpy( buffer, vals.data(), copyValues * sizeof( int ) );
return copyValues;
}
size_t MDAL::TuflowFVDataset3D::scalarVolumesData( size_t indexStart, size_t count, double *buffer )
{
if ( ( count < 1 ) || ( indexStart >= volumesCount() ) )
return 0;
if ( mTs >= mTimesteps )
return 0;
size_t copyValues = std::min( volumesCount() - indexStart, count );
std::vector<double> vals = mNcFile->readDoubleArr(
mNcidX,
mTs,
indexStart,
1,
copyValues
);
memcpy( buffer, vals.data(), copyValues * sizeof( double ) );
return copyValues;
}
size_t MDAL::TuflowFVDataset3D::vectorVolumesData( size_t indexStart, size_t count, double *buffer )
{
if ( ( count < 1 ) || ( indexStart >= volumesCount() ) )
return 0;
if ( mTs >= mTimesteps )
return 0;
size_t copyValues = std::min( volumesCount() - indexStart, count );
std::vector<double> vals_x = mNcFile->readDoubleArr(
mNcidX,
mTs,
indexStart,
1,
copyValues
);
std::vector<double> vals_y = mNcFile->readDoubleArr(
mNcidY,
mTs,
indexStart,
1,
copyValues
);
for ( size_t i = 0; i < copyValues; ++i )
{
buffer[2 * i] = vals_x[i];
buffer[2 * i + 1] = vals_y[i];
}
return copyValues;
}
size_t MDAL::TuflowFVDataset3D::activeVolumesData( size_t indexStart, size_t count, int *buffer )
{
// TODO
MDAL_UNUSED( indexStart )
memset( buffer, 1, count * sizeof( int ) );
return count;
}
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
MDAL::DriverTuflowFV::DriverTuflowFV():
DriverCF( "TUFLOWFV",
"TUFLOW FV",
"*.nc",
Capability::ReadMesh
)
{
}
MDAL::DriverTuflowFV::~DriverTuflowFV() = default;
MDAL::DriverTuflowFV *MDAL::DriverTuflowFV::create()
{
return new DriverTuflowFV();
}
MDAL::CFDimensions MDAL::DriverTuflowFV::populateDimensions( )
{
CFDimensions dims;
size_t count;
int ncid;
// 2D Mesh
mNcFile->getDimension( "NumCells2D", &count, &ncid );
dims.setDimension( CFDimensions::Face2D, count, ncid );
mNcFile->getDimension( "MaxNumCellVert", &count, &ncid );
dims.setDimension( CFDimensions::MaxVerticesInFace, count, ncid );
mNcFile->getDimension( "NumVert2D", &count, &ncid );
dims.setDimension( CFDimensions::Vertex2D, count, ncid );
// 3D Mesh
mNcFile->getDimension( "NumCells3D", &count, &ncid );
dims.setDimension( CFDimensions::Volume3D, count, ncid );
mNcFile->getDimension( "NumLayerFaces3D", &count, &ncid );
dims.setDimension( CFDimensions::StackedFace3D, count, ncid );
// Time
mNcFile->getDimension( "Time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );
return dims;
}
void MDAL::DriverTuflowFV::populateFacesAndVertices( Vertices &vertices, Faces &faces )
{
populateVertices( vertices );
populateFaces( faces );
}
void MDAL::DriverTuflowFV::populateVertices( MDAL::Vertices &vertices )
{
assert( vertices.empty() );
size_t vertexCount = mDimensions.size( CFDimensions::Vertex2D );
vertices.resize( vertexCount );
Vertex *vertexPtr = vertices.data();
// Parse 2D Mesh
const std::vector<double> vertices2D_x = mNcFile->readDoubleArr( "node_X", vertexCount );
const std::vector<double> vertices2D_y = mNcFile->readDoubleArr( "node_Y", vertexCount );
const std::vector<double> vertices2D_z = mNcFile->readDoubleArr( "node_Zb", vertexCount );
for ( size_t i = 0; i < vertexCount; ++i, ++vertexPtr )
{
vertexPtr->x = vertices2D_x[i];
vertexPtr->y = vertices2D_y[i];
vertexPtr->z = vertices2D_z[i];
}
}
void MDAL::DriverTuflowFV::populateFaces( MDAL::Faces &faces )
{
assert( faces.empty() );
size_t faceCount = mDimensions.size( CFDimensions::Face2D );
size_t vertexCount = mDimensions.size( CFDimensions::Vertex2D );
faces.resize( faceCount );
// Parse 2D Mesh
size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
std::vector<int> face_nodes_conn = mNcFile->readIntArr( "cell_node", faceCount * verticesInFace );
std::vector<int> face_vertex_counts = mNcFile->readIntArr( "cell_Nvert", faceCount );
for ( size_t i = 0; i < faceCount; ++i )
{
size_t nVertices = static_cast<size_t>( face_vertex_counts[i] );
std::vector<size_t> idxs;
for ( size_t j = 0; j < nVertices; ++j )
{
size_t idx = verticesInFace * i + j;
size_t val = static_cast<size_t>( face_nodes_conn[idx] - 1 ); //indexed from 1
assert( val < vertexCount );
idxs.push_back( val );
}
faces[i] = idxs;
}
}
void MDAL::DriverTuflowFV::calculateMaximumLevelCount()
{
if ( mMaximumLevelsCount < 0 )
{
mMaximumLevelsCount = 0;
int ncidVerticalLevels = mNcFile->arrId( "NL" );
if ( ncidVerticalLevels < 0 )
return;
const size_t maxBufferLength = 1000;
size_t indexStart = 0;
size_t facesCount = mDimensions.size( CFDimensions::Face2D );
while ( true )
{
size_t copyValues = std::min( facesCount - indexStart, maxBufferLength );
if ( copyValues <= 0 ) break;
std::vector<int> vals = mNcFile->readIntArr(
ncidVerticalLevels,
indexStart,
copyValues
);
mMaximumLevelsCount = std::max( mMaximumLevelsCount, *std::max_element( vals.begin(), vals.end() ) );
indexStart += copyValues;
}
}
}
void MDAL::DriverTuflowFV::addBedElevation( MDAL::MemoryMesh *mesh )
{
MDAL::addBedElevationDatasetGroup( mesh, mesh->vertices );
}
std::string MDAL::DriverTuflowFV::getCoordinateSystemVariableName()
{
return "";
}
std::set<std::string> MDAL::DriverTuflowFV::ignoreNetCDFVariables()
{
std::set<std::string> ignore_variables;
ignore_variables.insert( getTimeVariableName() );
ignore_variables.insert( "NL" );
ignore_variables.insert( "cell_Nvert" );
ignore_variables.insert( "cell_node" );
ignore_variables.insert( "idx2" );
ignore_variables.insert( "idx3" );
ignore_variables.insert( "cell_X" );
ignore_variables.insert( "cell_Y" );
ignore_variables.insert( "cell_Zb" );
ignore_variables.insert( "cell_A" );
ignore_variables.insert( "node_X" );
ignore_variables.insert( "node_Y" );
ignore_variables.insert( "node_Zb" );
ignore_variables.insert( "layerface_Z" );
ignore_variables.insert( "stat" );
return ignore_variables;
}
void MDAL::DriverTuflowFV::parseNetCDFVariableMetadata( int varid, const std::string &variableName, std::string &name, bool *is_vector, bool *is_x )
{
*is_vector = false;
*is_x = true;
std::string long_name = mNcFile->getAttrStr( "long_name", varid );
if ( long_name.empty() )
{
name = variableName;
}
else
{
if ( MDAL::startsWith( long_name, "maximum value of " ) )
long_name = MDAL::replace( long_name, "maximum value of ", "" ) + "/Maximums";
if ( MDAL::startsWith( long_name, "minimum value of " ) )
long_name = MDAL::replace( long_name, "minimum value of ", "" ) + "/Minimums";
if ( MDAL::startsWith( long_name, "time at maximum value of " ) )
long_name = MDAL::replace( long_name, "time at maximum value of ", "" ) + "/Time at Maximums";
if ( MDAL::startsWith( long_name, "time at minimum value of " ) )
long_name = MDAL::replace( long_name, "time at minimum value of ", "" ) + "/Time at Minimums";
if ( MDAL::startsWith( long_name, "x_" ) )
{
*is_vector = true;
name = MDAL::replace( long_name, "x_", "" );
}
else if ( MDAL::startsWith( long_name, "y_" ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( long_name, "y_", "" );
}
else
{
name = long_name;
}
}
}
std::string MDAL::DriverTuflowFV::getTimeVariableName() const
{
return "ResTime";
}
std::shared_ptr<MDAL::Dataset> MDAL::DriverTuflowFV::create3DDataset( std::shared_ptr<MDAL::DatasetGroup> group, size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
double, double )
{
calculateMaximumLevelCount();
std::shared_ptr<MDAL::TuflowFVDataset3D> dataset = std::make_shared<MDAL::TuflowFVDataset3D>(
group.get(),
dsi.ncid_x,
dsi.ncid_y,
dsi.nTimesteps,
mDimensions.size( CFDimensions::Type::Volume3D ),
mDimensions.size( CFDimensions::Type::Face2D ),
mDimensions.size( CFDimensions::Type::StackedFace3D ),
ts,
mMaximumLevelsCount,
mNcFile
);
// TODO use "Maximums" from file
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
return std::move( dataset );
}

106
external/mdal/frmts/mdal_tuflowfv.hpp vendored Normal file
View File

@ -0,0 +1,106 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#ifndef MDAL_TUFLOWFV_HPP
#define MDAL_TUFLOWFV_HPP
#include <string>
#include <memory>
#include <map>
#include <iostream>
#include <fstream>
#include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp"
#include "mdal.h"
#include "mdal_driver.hpp"
#include "mdal_cf.hpp"
namespace MDAL
{
class TuflowFVDataset3D: public Dataset3D
{
public:
TuflowFVDataset3D(
DatasetGroup *parent,
int ncid_x,
int ncid_y,
size_t timesteps,
size_t volumesCount,
size_t facesCount,
size_t levelFacesCount,
size_t ts,
size_t maximumLevelsCount,
std::shared_ptr<NetCDFFile> ncFile
);
virtual ~TuflowFVDataset3D() override;
size_t verticalLevelCountData( size_t indexStart, size_t count, int *buffer ) override;
size_t verticalLevelData( size_t indexStart, size_t count, double *buffer ) override;
size_t faceToVolumeData( size_t indexStart, size_t count, int *buffer ) override;
size_t scalarVolumesData( size_t indexStart, size_t count, double *buffer ) override;
size_t vectorVolumesData( size_t indexStart, size_t count, double *buffer ) override;
size_t activeVolumesData( size_t indexStart, size_t count, int *buffer ) override;
private:
int mNcidX; //!< NetCDF variable id
int mNcidY; //!< NetCDF variable id
size_t mTimesteps;
size_t mFacesCount;
size_t mLevelFacesCount;
size_t mTs;
std::shared_ptr<NetCDFFile> mNcFile;
int mNcidVerticalLevels = -1; //! variable id of int NL(NumCells2D) ;
int mNcidVerticalLevelsZ = -1; //! variable id of float layerface_Z(Time, NumLayerFaces3D) ;
int mNcidActive2D = -1; //! variable id of int stat(Time, NumCells2D) ;
int mNcid3DTo2D = -1; //! variable id of int idx2(NumCells3D) ;
int mNcid2DTo3D = -1; //! variable id of int idx3(NumCells2D) ;
};
/**
* TUFLOW FV format
*
* Binary NetCDF format with structure similar to UGRID stored as
* 3D Layered Mesh (https://github.com/qgis/QGIS-Enhancement-Proposals/issues/158)
*
* Both mesh and dataset is stored in single file.
*/
class DriverTuflowFV: public DriverCF
{
public:
DriverTuflowFV();
~DriverTuflowFV() override;
DriverTuflowFV *create() override;
private:
CFDimensions populateDimensions( ) override;
void populateFacesAndVertices( Vertices &vertices, Faces &faces ) override;
void addBedElevation( MemoryMesh *mesh ) override;
std::string getCoordinateSystemVariableName() override;
std::set<std::string> ignoreNetCDFVariables() override;
void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) override;
std::string getTimeVariableName() const override;
// TODO 2d dataset with active flag
// TODO CRS from prj file
std::shared_ptr<MDAL::Dataset> create3DDataset(
std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
double fill_val_x, double fill_val_y ) override;
void addBedElevationDatasetOnFaces();
void populateVertices( MDAL::Vertices &vertices );
void populateFaces( MDAL::Faces &faces );
void calculateMaximumLevelCount();
int mMaximumLevelsCount = -1;
};
} // namespace MDAL
#endif //MDAL_TUFLOWFV_HPP

294
external/mdal/frmts/mdal_ugrid.cpp vendored
View File

@ -8,13 +8,17 @@
#include <netcdf.h>
#include <assert.h>
#include <algorithm>
#include <cmath>
MDAL::DriverUgrid::DriverUgrid()
: DriverCF(
"Ugrid",
"UGRID Results",
"*.nc" )
"*.nc",
Capability::ReadMesh | Capability::SaveMesh )
{
}
MDAL::DriverUgrid *MDAL::DriverUgrid::create()
@ -24,13 +28,13 @@ MDAL::DriverUgrid *MDAL::DriverUgrid::create()
std::string MDAL::DriverUgrid::findMeshName( int dimension, bool optional ) const
{
const std::vector<std::string> variables = mNcFile.readArrNames();
const std::vector<std::string> variables = mNcFile->readArrNames();
for ( const std::string &varName : variables )
{
const std::string cf_role = mNcFile.getAttrStr( varName, "cf_role" );
const std::string cf_role = mNcFile->getAttrStr( varName, "cf_role" );
if ( cf_role == "mesh_topology" )
{
int topology_dimension = mNcFile.getAttrInt( varName, "topology_dimension" );
int topology_dimension = mNcFile->getAttrInt( varName, "topology_dimension" );
if ( topology_dimension == dimension )
{
return varName;
@ -45,12 +49,12 @@ std::string MDAL::DriverUgrid::findMeshName( int dimension, bool optional ) cons
std::string MDAL::DriverUgrid::nodeZVariableName() const
{
const std::vector<std::string> variables = mNcFile.readArrNames();
const std::vector<std::string> variables = mNcFile->readArrNames();
for ( const std::string &varName : variables )
{
const std::string stdName = mNcFile.getAttrStr( varName, "standard_name" );
const std::string meshName = mNcFile.getAttrStr( varName, "mesh" );
const std::string location = mNcFile.getAttrStr( varName, "location" );
const std::string stdName = mNcFile->getAttrStr( varName, "standard_name" );
const std::string meshName = mNcFile->getAttrStr( varName, "mesh" );
const std::string location = mNcFile->getAttrStr( varName, "location" );
if ( stdName == "altitude" && meshName == mMesh2dName && location == "node" )
{
@ -76,13 +80,13 @@ MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
//node dimension location is retrieved from the node variable
std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile.getAttrStr( mMesh2dName, "node_coordinates" ), ' ' );
std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile->getAttrStr( mMesh2dName, "node_coordinates" ), ' ' );
if ( nodeVariablesName.size() < 2 )
throw MDAL_Status::Err_UnknownFormat;
std::vector<size_t> nodeDimension;
std::vector<int> nodeDimensionId;
mNcFile.getDimensions( nodeVariablesName.at( 0 ), nodeDimension, nodeDimensionId );
mNcFile->getDimensions( nodeVariablesName.at( 0 ), nodeDimension, nodeDimensionId );
if ( nodeDimension.size() != 1 )
throw MDAL_Status::Err_UnknownFormat;
@ -92,8 +96,8 @@ MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
//face dimension location is retrieved from the face_node_connectivity variable
//if face_dimension is defined as attribute, the dimension at this location help to desambiguate vertex per faces and number of faces
std::string faceConnectivityVariablesName = mNcFile.getAttrStr( mMesh2dName, "face_node_connectivity" );
std::string faceDimensionLocation = mNcFile.getAttrStr( mMesh2dName, "face_dimension" );
std::string faceConnectivityVariablesName = mNcFile->getAttrStr( mMesh2dName, "face_node_connectivity" );
std::string faceDimensionLocation = mNcFile->getAttrStr( mMesh2dName, "face_dimension" );
if ( faceConnectivityVariablesName == "" )
throw MDAL_Status::Err_UnknownFormat;
@ -104,13 +108,13 @@ MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
std::vector<int> faceDimensionId;
int facesIndexDimensionId;
int maxVerticesPerFaceDimensionId;
mNcFile.getDimensions( faceConnectivityVariablesName, faceDimension, faceDimensionId );
mNcFile->getDimensions( faceConnectivityVariablesName, faceDimension, faceDimensionId );
if ( faceDimension.size() != 2 )
throw MDAL_Status::Err_UnknownFormat;
if ( faceDimensionLocation != "" )
{
mNcFile.getDimension( faceDimensionLocation, &facesCount, &ncid );
mNcFile->getDimension( faceDimensionLocation, &facesCount, &ncid );
if ( facesCount == faceDimension.at( 0 ) )
{
facesIndexDimensionId = faceDimensionId.at( 0 );
@ -138,10 +142,10 @@ MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
// number of edges in the mesh, not required for UGRID format
const std::string mesh2dEdge = mNcFile.getAttrStr( mMesh2dName, "edge_dimension" );
if ( mNcFile.hasDimension( mesh2dEdge ) )
const std::string mesh2dEdge = mNcFile->getAttrStr( mMesh2dName, "edge_dimension" );
if ( mNcFile->hasDimension( mesh2dEdge ) )
{
mNcFile.getDimension( mesh2dEdge, &count, &ncid );
mNcFile->getDimension( mesh2dEdge, &count, &ncid );
dims.setDimension( CFDimensions::Face2DEdge, count, ncid );
}
else
@ -150,9 +154,9 @@ MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
}
// Time not required for UGRID format
if ( mNcFile.hasDimension( "time" ) )
if ( mNcFile->hasDimension( "time" ) )
{
mNcFile.getDimension( "time", &count, &ncid );
mNcFile->getDimension( "time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );
}
else
@ -180,13 +184,13 @@ void MDAL::DriverUgrid::populateVertices( MDAL::Vertices &vertices )
// node_coordinates should be something like Mesh2D_node_x Mesh2D_node_y
std::string verticesXName, verticesYName;
parse2VariablesFromAttribute( mMesh2dName, "node_coordinates", verticesXName, verticesYName, false );
const std::vector<double> vertices2D_x = mNcFile.readDoubleArr( verticesXName, vertexCount );
const std::vector<double> vertices2D_y = mNcFile.readDoubleArr( verticesYName, vertexCount );
const std::vector<double> vertices2D_x = mNcFile->readDoubleArr( verticesXName, vertexCount );
const std::vector<double> vertices2D_y = mNcFile->readDoubleArr( verticesYName, vertexCount );
std::vector<double> vertices2D_z;
if ( mNcFile.hasArr( nodeZVariableName() ) )
if ( mNcFile->hasArr( nodeZVariableName() ) )
{
vertices2D_z = mNcFile.readDoubleArr( nodeZVariableName(), vertexCount );
vertices2D_z = mNcFile->readDoubleArr( nodeZVariableName(), vertexCount );
}
for ( size_t i = 0; i < vertexCount; ++i, ++vertexPtr )
@ -206,14 +210,14 @@ void MDAL::DriverUgrid::populateFaces( MDAL::Faces &faces )
// Parse 2D Mesh
// face_node_connectivity is usually something like Mesh2D_face_nodes
const std::string mesh2dFaceNodeConnectivity = mNcFile.getAttrStr( mMesh2dName, "face_node_connectivity" );
const std::string mesh2dFaceNodeConnectivity = mNcFile->getAttrStr( mMesh2dName, "face_node_connectivity" );
size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
int fill_val = -1;
if ( mNcFile.hasAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" ) )
fill_val = mNcFile.getAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" );
int start_index = mNcFile.getAttrInt( mesh2dFaceNodeConnectivity, "start_index" );
std::vector<int> face_nodes_conn = mNcFile.readIntArr( mesh2dFaceNodeConnectivity, faceCount * verticesInFace );
if ( mNcFile->hasAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" ) )
fill_val = mNcFile->getAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" );
int start_index = mNcFile->getAttrInt( mesh2dFaceNodeConnectivity, "start_index" );
std::vector<int> face_nodes_conn = mNcFile->readIntArr( mesh2dFaceNodeConnectivity, faceCount * verticesInFace );
for ( size_t i = 0; i < faceCount; ++i )
{
@ -244,7 +248,7 @@ void MDAL::DriverUgrid::populateFaces( MDAL::Faces &faces )
void MDAL::DriverUgrid::addBedElevation( MDAL::MemoryMesh *mesh )
{
if ( mNcFile.hasArr( nodeZVariableName() ) ) MDAL::addBedElevationDatasetGroup( mesh, mesh->vertices );
if ( mNcFile->hasArr( nodeZVariableName() ) ) MDAL::addBedElevationDatasetGroup( mesh, mesh->vertices );
}
std::string MDAL::DriverUgrid::getCoordinateSystemVariableName()
@ -252,12 +256,12 @@ std::string MDAL::DriverUgrid::getCoordinateSystemVariableName()
std::string coordinate_system_variable;
// first try to get the coordinate system variable from grid definition
std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile.getAttrStr( mMesh2dName, "node_coordinates" ), ' ' );
std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile->getAttrStr( mMesh2dName, "node_coordinates" ), ' ' );
if ( nodeVariablesName.size() > 1 )
{
if ( mNcFile.hasArr( nodeVariablesName[0] ) )
if ( mNcFile->hasArr( nodeVariablesName[0] ) )
{
coordinate_system_variable = mNcFile.getAttrStr( nodeVariablesName[0], "grid_mapping" );
coordinate_system_variable = mNcFile->getAttrStr( nodeVariablesName[0], "grid_mapping" );
}
}
@ -265,9 +269,9 @@ std::string MDAL::DriverUgrid::getCoordinateSystemVariableName()
// if automatic discovery fails, try to check some hardcoded common variables that store projection
if ( coordinate_system_variable.empty() )
{
if ( mNcFile.hasArr( "projected_coordinate_system" ) )
if ( mNcFile->hasArr( "projected_coordinate_system" ) )
coordinate_system_variable = "projected_coordinate_system";
else if ( mNcFile.hasArr( "wgs84" ) )
else if ( mNcFile->hasArr( "wgs84" ) )
coordinate_system_variable = "wgs84";
}
@ -284,7 +288,7 @@ std::set<std::string> MDAL::DriverUgrid::ignoreNetCDFVariables()
ignore_variables.insert( "timestep" );
std::vector<std::string> meshes;
if ( mNcFile.hasArr( mMesh1dName ) )
if ( mNcFile->hasArr( mMesh1dName ) )
meshes.push_back( mMesh1dName );
meshes.push_back( mMesh2dName );
@ -295,31 +299,31 @@ std::set<std::string> MDAL::DriverUgrid::ignoreNetCDFVariables()
parse2VariablesFromAttribute( mesh, "node_coordinates", xName, yName, true );
ignore_variables.insert( xName );
ignore_variables.insert( yName );
ignore_variables.insert( mNcFile.getAttrStr( mesh, "edge_node_connectivity" ) );
ignore_variables.insert( mNcFile->getAttrStr( mesh, "edge_node_connectivity" ) );
parse2VariablesFromAttribute( mesh, "edge_coordinates", xName, yName, true );
if ( !xName.empty() )
{
ignore_variables.insert( xName );
ignore_variables.insert( mNcFile.getAttrStr( xName, "bounds" ) );
ignore_variables.insert( mNcFile->getAttrStr( xName, "bounds" ) );
}
if ( !yName.empty() )
{
ignore_variables.insert( yName );
ignore_variables.insert( mNcFile.getAttrStr( yName, "bounds" ) );
ignore_variables.insert( mNcFile->getAttrStr( yName, "bounds" ) );
}
ignore_variables.insert( mNcFile.getAttrStr( mesh, "face_node_connectivity" ) );
ignore_variables.insert( mNcFile->getAttrStr( mesh, "face_node_connectivity" ) );
parse2VariablesFromAttribute( mesh, "face_coordinates", xName, yName, true );
if ( !xName.empty() )
{
ignore_variables.insert( xName );
ignore_variables.insert( mNcFile.getAttrStr( xName, "bounds" ) );
ignore_variables.insert( mNcFile->getAttrStr( xName, "bounds" ) );
}
if ( !yName.empty() )
{
ignore_variables.insert( yName );
ignore_variables.insert( mNcFile.getAttrStr( yName, "bounds" ) );
ignore_variables.insert( mNcFile->getAttrStr( yName, "bounds" ) );
}
ignore_variables.insert( mNcFile.getAttrStr( mesh, "edge_face_connectivity" ) );
ignore_variables.insert( mNcFile->getAttrStr( mesh, "edge_face_connectivity" ) );
}
return ignore_variables;
@ -330,10 +334,10 @@ void MDAL::DriverUgrid::parseNetCDFVariableMetadata( int varid, const std::strin
*is_vector = false;
*is_x = true;
std::string long_name = mNcFile.getAttrStr( "long_name", varid );
std::string long_name = mNcFile->getAttrStr( "long_name", varid );
if ( long_name.empty() )
{
std::string standard_name = mNcFile.getAttrStr( "standard_name", varid );
std::string standard_name = mNcFile->getAttrStr( "standard_name", varid );
if ( standard_name.empty() )
{
name = variableName;
@ -379,10 +383,15 @@ void MDAL::DriverUgrid::parseNetCDFVariableMetadata( int varid, const std::strin
}
}
std::string MDAL::DriverUgrid::getTimeVariableName() const
{
return "time";
}
void MDAL::DriverUgrid::parse2VariablesFromAttribute( const std::string &name, const std::string &attr_name,
std::string &var1, std::string &var2, bool optional ) const
{
const std::string mesh2dNodeCoordinates = mNcFile.getAttrStr( name, attr_name );
const std::string mesh2dNodeCoordinates = mNcFile->getAttrStr( name, attr_name );
const std::vector<std::string> chunks = MDAL::split( mesh2dNodeCoordinates, ' ' );
if ( chunks.size() != 2 )
@ -401,3 +410,198 @@ void MDAL::DriverUgrid::parse2VariablesFromAttribute( const std::string &name, c
var2 = chunks[1];
}
}
void MDAL::DriverUgrid::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Status *status )
{
mFileName = uri;
try
{
// Create file
mNcFile.reset( new NetCDFFile );
mNcFile->createFile( mFileName );
// Write globals
writeGlobals( );
// Write variables
writeVariables( mesh );
}
catch ( MDAL_Status error )
{
if ( status ) *status = ( error );
}
}
void MDAL::DriverUgrid::writeVariables( MDAL::Mesh *mesh )
{
// Global dimensions
int dimNodeCountId = mNcFile->defineDimension( "nmesh2d_node", mesh->verticesCount() );
int dimFaceCountId = mNcFile->defineDimension( "nmesh2d_face", mesh->facesCount() );
mNcFile->defineDimension( "nmesh2d_edge", 1 ); // no data on edges, cannot be 0, since 0==NC_UNLIMITED
int dimTimeId = mNcFile->defineDimension( "time", NC_UNLIMITED );
int dimMaxNodesPerFaceId = mNcFile->defineDimension( "max_nmesh2d_face_nodes",
mesh->faceVerticesMaximumCount() );
// Mesh 2D Definition
int mesh2dId = mNcFile->defineVar( "mesh2d", NC_INT, 0, nullptr );
mNcFile->putAttrStr( mesh2dId, "cf_role", "mesh_topology" );
mNcFile->putAttrStr( mesh2dId, "long_name", "Topology data of 2D network" );
mNcFile->putAttrInt( mesh2dId, "topology_dimension", 2 );
mNcFile->putAttrStr( mesh2dId, "node_coordinates", "mesh2d_node_x mesh2d_node_y" );
mNcFile->putAttrStr( mesh2dId, "node_dimension", "nmesh2d_node" );
mNcFile->putAttrStr( mesh2dId, "edge_dimension", "nmesh2d_edge" );
mNcFile->putAttrStr( mesh2dId, "max_face_nodes_dimension", "max_nmesh2d_face_nodes" );
mNcFile->putAttrStr( mesh2dId, "face_node_connectivity", "mesh2d_face_nodes" );
mNcFile->putAttrStr( mesh2dId, "face_dimension", "nmesh2d_face" );
// Nodes X coordinate
int mesh2dNodeXId = mNcFile->defineVar( "mesh2d_node_x", NC_DOUBLE, 1, &dimNodeCountId );
mNcFile->putAttrStr( mesh2dNodeXId, "standard_name", "projection_x_coordinate" );
mNcFile->putAttrStr( mesh2dNodeXId, "long_name", "x-coordinate of mesh nodes" );
mNcFile->putAttrStr( mesh2dNodeXId, "mesh", "mesh2d" );
mNcFile->putAttrStr( mesh2dNodeXId, "location", "node" );
// Nodes Y coordinate
int mesh2dNodeYId = mNcFile->defineVar( "mesh2d_node_y", NC_DOUBLE, 1, &dimNodeCountId );
mNcFile->putAttrStr( mesh2dNodeYId, "standard_name", "projection_y_coordinate" );
mNcFile->putAttrStr( mesh2dNodeYId, "long_name", "y-coordinate of mesh nodes" );
mNcFile->putAttrStr( mesh2dNodeYId, "mesh", "mesh2d" );
mNcFile->putAttrStr( mesh2dNodeYId, "location", "node" );
// Nodes Z coordinate
int mesh2dNodeZId = mNcFile->defineVar( "mesh2d_node_z", NC_DOUBLE, 1, &dimNodeCountId );
mNcFile->putAttrStr( mesh2dNodeZId, "mesh", "mesh2d" );
mNcFile->putAttrStr( mesh2dNodeZId, "location", "node" );
mNcFile->putAttrStr( mesh2dNodeZId, "coordinates", "mesh2d_node_x mesh2d_node_y" );
mNcFile->putAttrStr( mesh2dNodeZId, "standard_name", "altitude" );
mNcFile->putAttrStr( mesh2dNodeZId, "long_name", "z-coordinate of mesh nodes" );
mNcFile->putAttrStr( mesh2dNodeZId, "grid_mapping", "projected_coordinate_system" );
double fillNodeZCoodVal = -999.0;
mNcFile->putAttrDouble( mesh2dNodeZId, "_FillValue", fillNodeZCoodVal );
// Faces 2D Variable
int mesh2FaceNodesId_dimIds [] { dimFaceCountId, dimMaxNodesPerFaceId };
int mesh2FaceNodesId = mNcFile->defineVar( "mesh2d_face_nodes", NC_INT, 2, mesh2FaceNodesId_dimIds );
mNcFile->putAttrStr( mesh2FaceNodesId, "cf_role", "face_node_connectivity" );
mNcFile->putAttrStr( mesh2FaceNodesId, "mesh", "mesh2d" );
mNcFile->putAttrStr( mesh2FaceNodesId, "location", "face" );
mNcFile->putAttrStr( mesh2FaceNodesId, "long_name", "Mapping from every face to its corner nodes (counterclockwise)" );
mNcFile->putAttrInt( mesh2FaceNodesId, "start_index", 0 );
int fillFace2DVertexValue = -999;
mNcFile->putAttrInt( mesh2FaceNodesId, "_FillValue", fillFace2DVertexValue );
// Projected Coordinate System
int pcsId = mNcFile->defineVar( "projected_coordinate_system", NC_INT, 0, nullptr );
if ( mesh->crs() == "" )
{
mNcFile->putAttrInt( pcsId, "epsg", 0 );
mNcFile->putAttrStr( pcsId, "EPSG_code", "epgs:0" );
}
else
{
std::vector<std::string> words = MDAL::split( mesh->crs(), ":" );
if ( words[0] == "EPSG" && words.size() > 1 )
{
mNcFile->putAttrInt( pcsId, "epsg", std::stoi( words[1] ) );
mNcFile->putAttrStr( pcsId, "EPSG_code", mesh->crs() );
}
else
{
mNcFile->putAttrStr( pcsId, "wkt", mesh->crs() );
}
}
// Time array
int timeId = mNcFile->defineVar( "time", NC_DOUBLE, 1, &dimTimeId );
mNcFile->putAttrStr( timeId, "units", "hours since 2000-01-01 00:00:00" );
// Turning off define mode - allows data write
nc_enddef( mNcFile->handle() );
// Write vertices
const size_t maxBufferSize = 1000;
const size_t bufferSize = std::min( mesh->verticesCount(), maxBufferSize );
const size_t verticesCoordCount = bufferSize * 3;
std::vector<double> verticesCoordinates( verticesCoordCount );
std::unique_ptr<MDAL::MeshVertexIterator> vertexIterator = mesh->readVertices();
size_t vertexIndex = 0;
size_t vertexFileIndex = 0;
while ( vertexIndex < mesh->verticesCount() )
{
size_t verticesRead = vertexIterator->next( bufferSize, verticesCoordinates.data() );
if ( verticesRead == 0 )
break;
for ( size_t i = 0; i < verticesRead; i++ )
{
mNcFile->putDataDouble( mesh2dNodeXId, vertexFileIndex, verticesCoordinates[3 * i] );
mNcFile->putDataDouble( mesh2dNodeYId, vertexFileIndex, verticesCoordinates[3 * i + 1] );
if ( std::isnan( verticesCoordinates[3 * i + 2] ) )
mNcFile->putDataDouble( mesh2dNodeZId, vertexFileIndex, fillNodeZCoodVal );
else
mNcFile->putDataDouble( mesh2dNodeZId, vertexFileIndex, verticesCoordinates[3 * i + 2] );
vertexFileIndex++;
}
vertexIndex += verticesRead;
}
// Write faces
std::unique_ptr<MDAL::MeshFaceIterator> faceIterator = mesh->readFaces();
const size_t faceVerticesMax = mesh->faceVerticesMaximumCount();
const size_t facesCount = mesh->facesCount();
const size_t faceOffsetsBufferLen = std::min( facesCount, maxBufferSize );
const size_t vertexIndicesBufferLen = faceOffsetsBufferLen * faceVerticesMax;
std::vector<int> faceOffsetsBuffer( faceOffsetsBufferLen );
std::vector<int> vertexIndicesBuffer( vertexIndicesBufferLen );
size_t faceIndex = 0;
while ( faceIndex < facesCount )
{
size_t facesRead = faceIterator->next(
faceOffsetsBufferLen,
faceOffsetsBuffer.data(),
vertexIndicesBufferLen,
vertexIndicesBuffer.data() );
if ( facesRead == 0 )
break;
for ( size_t i = 0; i < facesRead; i++ )
{
std::vector<int> verticesFaceData( faceVerticesMax, fillFace2DVertexValue );
int startIndex = 0;
if ( i > 0 )
startIndex = faceOffsetsBuffer[ i - 1 ];
int endIndex = faceOffsetsBuffer[ i ];
size_t k = 0;
for ( int j = startIndex; j < endIndex; ++j )
{
int vertexIndex = vertexIndicesBuffer[ static_cast<size_t>( j ) ];
verticesFaceData[k++] = vertexIndex;
}
mNcFile->putDataArrayInt( mesh2FaceNodesId, faceIndex + i, faceVerticesMax, verticesFaceData.data() );
}
faceIndex += facesRead;
}
// Time values (not implemented)
mNcFile->putDataDouble( timeId, 0, 0.0 );
// Turning on define mode
nc_redef( mNcFile->handle() );
}
void MDAL::DriverUgrid::writeGlobals()
{
mNcFile->putAttrStr( NC_GLOBAL, "source", "MDAL " + std::string( MDAL_Version() ) );
mNcFile->putAttrStr( NC_GLOBAL, "date_created", MDAL::getCurrentTimeStamp() );
mNcFile->putAttrStr( NC_GLOBAL, "Conventions", "CF-1.6 UGRID-1.0" );
}

8
external/mdal/frmts/mdal_ugrid.hpp vendored
View File

@ -26,6 +26,7 @@ namespace MDAL
DriverUgrid();
~DriverUgrid() override = default;
DriverUgrid *create() override;
void save( const std::string &uri, Mesh *mesh, MDAL_Status *status ) override;
private:
CFDimensions populateDimensions( ) override;
@ -37,6 +38,7 @@ namespace MDAL
std::set<std::string> ignoreNetCDFVariables() override;
void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) override;
std::string getTimeVariableName() const override;
void parse2VariablesFromAttribute( const std::string &name, const std::string &attr_name,
std::string &var1, std::string &var2,
@ -45,6 +47,12 @@ namespace MDAL
std::string mMesh2dName;
std::string mMesh1dName;
std::string nodeZVariableName() const;
void writeDimensions( MDAL::Mesh *mesh );
void writeVariables( MDAL::Mesh *mesh );
void writeGlobals();
int faceVerticesMaximumCount() const override
{ return std::numeric_limits<int>::max(); }
};
} // namespace MDAL

16
external/mdal/frmts/mdal_xdmf.cpp vendored
View File

@ -21,7 +21,7 @@
#include "mdal_xml.hpp"
MDAL::XdmfDataset::XdmfDataset( MDAL::DatasetGroup *grp, const MDAL::HyperSlab &slab, const HdfDataset &valuesDs, double time )
: MDAL::Dataset( grp )
: MDAL::Dataset2D( grp )
, mHdf5DatasetValues( valuesDs )
, mHyperSlab( slab )
{
@ -103,7 +103,7 @@ size_t MDAL::XdmfDataset::vectorData( size_t indexStart, size_t count, double *b
size_t MDAL::XdmfDataset::activeData( size_t indexStart, size_t count, int *buffer )
{
MDAL_UNUSED( indexStart );
MDAL_UNUSED( indexStart )
memset( buffer, 1, count * sizeof( int ) );
return count;
}
@ -117,13 +117,13 @@ MDAL::XdmfFunctionDataset::XdmfFunctionDataset(
MDAL::DatasetGroup *grp,
MDAL::XdmfFunctionDataset::FunctionType type,
double time )
: MDAL::Dataset( grp )
: MDAL::Dataset2D( grp )
, mType( type )
, mBaseReferenceGroup( "XDMF", grp->mesh(), grp->uri() )
{
setTime( time );
mBaseReferenceGroup.setIsScalar( true );
mBaseReferenceGroup.setIsOnVertices( grp->isOnVertices() );
mBaseReferenceGroup.setDataLocation( grp->dataLocation() );
mBaseReferenceGroup.setName( "Base group for reference datasets" );
}
@ -171,7 +171,7 @@ size_t MDAL::XdmfFunctionDataset::vectorData( size_t indexStart, size_t count, d
size_t MDAL::XdmfFunctionDataset::activeData( size_t indexStart, size_t count, int *buffer )
{
MDAL_UNUSED( indexStart );
MDAL_UNUSED( indexStart )
memset( buffer, 1, count * sizeof( int ) );
return count;
}
@ -341,7 +341,7 @@ HdfDataset MDAL::DriverXdmf::parseHdf5Node( const XMLFile &xmfFile, xmlNodePtr n
std::shared_ptr<HdfFile> hdfFile;
if ( mHdfFiles.count( hdf5Name ) == 0 )
{
hdfFile = std::make_shared<HdfFile>( hdf5Name );
hdfFile = std::make_shared<HdfFile>( hdf5Name, HdfFile::ReadOnly );
mHdfFiles[hdf5Name] = hdfFile;
}
else
@ -595,7 +595,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXdmf::findGroup( std::map<std::s
groupName
);
group->setIsScalar( isScalar );
group->setIsOnVertices( false ); //only center-based implemented
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D ); //only center-based implemented
groups[groupName] = group;
}
else
@ -626,7 +626,7 @@ MDAL::DriverXdmf *MDAL::DriverXdmf::create()
return new DriverXdmf();
}
bool MDAL::DriverXdmf::canRead( const std::string &uri )
bool MDAL::DriverXdmf::canReadDatasets( const std::string &uri )
{
XMLFile xmfFile;
try

6
external/mdal/frmts/mdal_xdmf.hpp vendored
View File

@ -55,7 +55,7 @@ namespace MDAL
* </DataItem>
* </Attribute>
*/
class XdmfDataset: public Dataset
class XdmfDataset: public Dataset2D
{
public:
XdmfDataset( DatasetGroup *grp,
@ -99,7 +99,7 @@ namespace MDAL
* </DataItem>
* </Attribute>
*/
class XdmfFunctionDataset: public Dataset
class XdmfFunctionDataset: public Dataset2D
{
public:
enum FunctionType
@ -159,7 +159,7 @@ namespace MDAL
~DriverXdmf( ) override;
DriverXdmf *create() override;
bool canRead( const std::string &uri ) override;
bool canReadDatasets( const std::string &uri ) override;
void load( const std::string &datFile, Mesh *mesh, MDAL_Status *status ) override;
private:

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

@ -17,7 +17,7 @@
MDAL::XmdfDataset::~XmdfDataset() = default;
MDAL::XmdfDataset::XmdfDataset( DatasetGroup *grp, const HdfDataset &valuesDs, const HdfDataset &activeDs, hsize_t timeIndex )
: Dataset( grp )
: Dataset2D( grp )
, mHdf5DatasetValues( valuesDs )
, mHdf5DatasetActive( activeDs )
, mTimeIndex( timeIndex )
@ -98,9 +98,9 @@ MDAL::DriverXmdf *MDAL::DriverXmdf::create()
return new DriverXmdf();
}
bool MDAL::DriverXmdf::canRead( const std::string &uri )
bool MDAL::DriverXmdf::canReadDatasets( const std::string &uri )
{
HdfFile file( uri );
HdfFile file( uri, HdfFile::ReadOnly );
if ( !file.isValid() )
{
return false;
@ -121,7 +121,7 @@ void MDAL::DriverXmdf::load( const std::string &datFile, MDAL::Mesh *mesh, MDAL
mMesh = mesh;
if ( status ) *status = MDAL_Status::None;
HdfFile file( mDatFile );
HdfFile file( mDatFile, HdfFile::ReadOnly );
if ( !file.isValid() )
{
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
@ -257,7 +257,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXmdf::readXmdfGroupAsDatasetGrou
groupName
);
group->setIsScalar( !isVector );
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
group->setMetadata( "TIMEUNITS", timeUnit );
// lazy loading of min and max of the dataset group

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

@ -31,7 +31,7 @@ namespace MDAL
* 2D arrays (time, x) for scalar datasets
* 2D arrays (time, active) for active flags
*/
class XmdfDataset: public Dataset
class XmdfDataset: public Dataset2D
{
public:
XmdfDataset( DatasetGroup *grp,
@ -82,7 +82,7 @@ namespace MDAL
~DriverXmdf( ) override = default;
DriverXmdf *create() override;
bool canRead( const std::string &uri ) override;
bool canReadDatasets( const std::string &uri ) override;
void load( const std::string &datFile, Mesh *mesh, MDAL_Status *status ) override;
private:

142
external/mdal/mdal.cpp vendored
View File

@ -22,7 +22,7 @@ static MDAL_Status sLastStatus;
const char *MDAL_Version()
{
return "0.4.0";
return "0.4.90";
}
MDAL_Status MDAL_LastStatus()
@ -81,18 +81,20 @@ bool MDAL_DR_meshLoadCapability( DriverH driver )
return d->hasCapability( MDAL::Capability::ReadMesh );
}
bool MDAL_DR_writeDatasetsCapability( DriverH driver )
bool MDAL_DR_writeDatasetsCapability( DriverH driver, MDAL_DataLocation location )
{
if ( !driver )
{
sLastStatus = MDAL_Status::Err_MissingDriver;
return false;
}
MDAL::Driver *d = static_cast< MDAL::Driver * >( driver );
return d->hasCapability( MDAL::Capability::WriteDatasets );
return d->hasWriteDatasetCapability( location );
}
bool MDAL_DR_SaveMeshCapability( DriverH driver )
bool MDAL_DR_saveMeshCapability( DriverH driver )
{
if ( !driver )
{
@ -328,7 +330,7 @@ DatasetGroupH MDAL_M_datasetGroup( MeshH mesh, int index )
DatasetGroupH MDAL_M_addDatasetGroup(
MeshH mesh,
const char *name,
bool isOnVertices,
MDAL_DataLocation dataLocation,
bool hasScalarData,
DriverH driver,
const char *datasetGroupFile )
@ -360,7 +362,7 @@ DatasetGroupH MDAL_M_addDatasetGroup(
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
MDAL::Driver *dr = static_cast< MDAL::Driver * >( driver );
if ( !dr->hasCapability( MDAL::Capability::WriteDatasets ) )
if ( !dr->hasWriteDatasetCapability( dataLocation ) )
{
sLastStatus = MDAL_Status::Err_MissingDriverCapability;
return nullptr;
@ -369,7 +371,7 @@ DatasetGroupH MDAL_M_addDatasetGroup(
const size_t index = m->datasetGroups.size();
dr->createDatasetGroup( m,
name,
isOnVertices,
dataLocation,
hasScalarData,
datasetGroupFile
);
@ -602,15 +604,15 @@ bool MDAL_G_hasScalarData( DatasetGroupH group )
return g->isScalar();
}
bool MDAL_G_isOnVertices( DatasetGroupH group )
MDAL_DataLocation MDAL_G_dataLocation( DatasetGroupH group )
{
if ( !group )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return true;
return DataInvalidLocation;
}
MDAL::DatasetGroup *g = static_cast< MDAL::DatasetGroup * >( group );
return g->isOnVertices();
return g->dataLocation();
}
void MDAL_G_minimumMaximum( DatasetGroupH group, double *min, double *max )
@ -664,7 +666,7 @@ DatasetH MDAL_G_addDataset( DatasetGroupH group, double time, const double *valu
return nullptr;
}
if ( !dr->hasCapability( MDAL::Capability::WriteDatasets ) )
if ( !dr->hasWriteDatasetCapability( g->dataLocation() ) )
{
sLastStatus = MDAL_Status::Err_MissingDriverCapability;
return nullptr;
@ -718,7 +720,7 @@ void MDAL_G_closeEditMode( DatasetGroupH group )
return;
}
if ( !dr->hasCapability( MDAL::Capability::WriteDatasets ) )
if ( !dr->hasWriteDatasetCapability( g->dataLocation() ) )
{
sLastStatus = MDAL_Status::Err_MissingDriverCapability;
return;
@ -805,6 +807,30 @@ double MDAL_D_time( DatasetH dataset )
}
int MDAL_D_volumesCount( DatasetH dataset )
{
if ( !dataset )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
int len = static_cast<int>( d->volumesCount() );
return len;
}
int MDAL_D_maximumVerticalLevelCount( DatasetH dataset )
{
if ( !dataset )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
int len = static_cast<int>( d->maximumVerticalLevelsCount() );
return len;
}
int MDAL_D_valueCount( DatasetH dataset )
{
if ( !dataset )
@ -855,6 +881,11 @@ int MDAL_D_data( DatasetH dataset, int indexStart, int count, MDAL_DataType data
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
if ( ( g->dataLocation() != MDAL_DataLocation::DataOnVertices2D ) && ( g->dataLocation() != MDAL_DataLocation::DataOnFaces2D ) )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = d->valuesCount();
break;
case MDAL_DataType::VECTOR_2D_DOUBLE:
@ -863,11 +894,79 @@ int MDAL_D_data( DatasetH dataset, int indexStart, int count, MDAL_DataType data
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
if ( ( g->dataLocation() != MDAL_DataLocation::DataOnVertices2D ) && ( g->dataLocation() != MDAL_DataLocation::DataOnFaces2D ) )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = d->valuesCount();
break;
case MDAL_DataType::ACTIVE_INTEGER:
if ( ( g->dataLocation() != MDAL_DataLocation::DataOnVertices2D ) && ( g->dataLocation() != MDAL_DataLocation::DataOnFaces2D ) )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = m->facesCount();
break;
case MDAL_DataType::VERTICAL_LEVEL_COUNT_INTEGER:
if ( g->dataLocation() != MDAL_DataLocation::DataOnVolumes3D )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = m->facesCount();
break;
case MDAL_DataType::VERTICAL_LEVEL_DOUBLE:
if ( g->dataLocation() != MDAL_DataLocation::DataOnVolumes3D )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = m->facesCount() + d->volumesCount();
break;
case MDAL_DataType::FACE_INDEX_TO_VOLUME_INDEX_INTEGER:
if ( g->dataLocation() != MDAL_DataLocation::DataOnVolumes3D )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = m->facesCount();
break;
case MDAL_DataType::SCALAR_VOLUMES_DOUBLE:
if ( g->dataLocation() != MDAL_DataLocation::DataOnVolumes3D )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
if ( !g->isScalar() )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = d->volumesCount();
break;
case MDAL_DataType::VECTOR_2D_VOLUMES_DOUBLE:
if ( g->dataLocation() != MDAL_DataLocation::DataOnVolumes3D )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
if ( g->isScalar() )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = 2 * d->volumesCount();
break;
case MDAL_DataType::ACTIVE_VOLUMES_INTEGER:
if ( g->dataLocation() != MDAL_DataLocation::DataOnVolumes3D )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = d->volumesCount();
break;
}
// Check that we are not reaching out of values limit
@ -896,6 +995,24 @@ int MDAL_D_data( DatasetH dataset, int indexStart, int count, MDAL_DataType data
case MDAL_DataType::ACTIVE_INTEGER:
writtenValuesCount = d->activeData( indexStartSizeT, countSizeT, static_cast<int *>( buffer ) );
break;
case MDAL_DataType::VERTICAL_LEVEL_COUNT_INTEGER:
writtenValuesCount = d->verticalLevelCountData( indexStartSizeT, countSizeT, static_cast<int *>( buffer ) );
break;
case MDAL_DataType::VERTICAL_LEVEL_DOUBLE:
writtenValuesCount = d->verticalLevelData( indexStartSizeT, countSizeT, static_cast<double *>( buffer ) );
break;
case MDAL_DataType::FACE_INDEX_TO_VOLUME_INDEX_INTEGER:
writtenValuesCount = d->faceToVolumeData( indexStartSizeT, countSizeT, static_cast<int *>( buffer ) );
break;
case MDAL_DataType::SCALAR_VOLUMES_DOUBLE:
writtenValuesCount = d->scalarVolumesData( indexStartSizeT, countSizeT, static_cast<double *>( buffer ) );
break;
case MDAL_DataType::VECTOR_2D_VOLUMES_DOUBLE:
writtenValuesCount = d->vectorVolumesData( indexStartSizeT, countSizeT, static_cast<double *>( buffer ) );
break;
case MDAL_DataType::ACTIVE_VOLUMES_INTEGER:
writtenValuesCount = d->activeVolumesData( indexStartSizeT, countSizeT, static_cast<int *>( buffer ) );
break;
}
return static_cast<int>( writtenValuesCount );
@ -922,4 +1039,3 @@ void MDAL_D_minimumMaximum( DatasetH dataset, double *min, double *max )
*min = stats.minimum;
*max = stats.maximum;
}

74
external/mdal/mdal_data_model.cpp vendored
View File

@ -17,20 +17,16 @@ MDAL::Dataset::Dataset( MDAL::DatasetGroup *parent )
assert( mParent );
}
std::string MDAL::Dataset::driverName() const
{
return group()->driverName();
}
size_t MDAL::Dataset::valuesCount() const
{
if ( group()->isOnVertices() )
const MDAL_DataLocation location = group()->dataLocation();
switch ( location )
{
return mesh()->verticesCount();
}
else
{
return mesh()->facesCount();
case MDAL_DataLocation::DataOnVertices2D: return mesh()->verticesCount();
case MDAL_DataLocation::DataOnFaces2D: return mesh()->facesCount();
case MDAL_DataLocation::DataOnVolumes3D: return volumesCount();
default: return 0;
}
}
@ -69,11 +65,55 @@ bool MDAL::Dataset::isValid() const
return mIsValid;
}
void MDAL::Dataset::setIsValid( bool isValid )
MDAL::Dataset2D::Dataset2D( MDAL::DatasetGroup *parent )
: Dataset( parent )
{
mIsValid = isValid;
}
MDAL::Dataset2D::~Dataset2D() = default;
size_t MDAL::Dataset2D::volumesCount() const { return 0; }
size_t MDAL::Dataset2D::maximumVerticalLevelsCount() const { return 0; }
size_t MDAL::Dataset2D::verticalLevelCountData( size_t, size_t, int * ) { return 0; }
size_t MDAL::Dataset2D::verticalLevelData( size_t, size_t, double * ) { return 0; }
size_t MDAL::Dataset2D::faceToVolumeData( size_t, size_t, int * ) { return 0; }
size_t MDAL::Dataset2D::scalarVolumesData( size_t, size_t, double * ) { return 0; }
size_t MDAL::Dataset2D::vectorVolumesData( size_t, size_t, double * ) { return 0; }
size_t MDAL::Dataset2D::activeVolumesData( size_t, size_t, int * ) { return 0; }
MDAL::Dataset3D::Dataset3D( MDAL::DatasetGroup *parent, size_t volumes, size_t maxVerticalLevelCount )
: Dataset( parent )
, mVolumesCount( volumes )
, mMaximumVerticalLevelsCount( maxVerticalLevelCount )
{
}
MDAL::Dataset3D::~Dataset3D() = default;
size_t MDAL::Dataset3D::volumesCount() const
{
return mVolumesCount;
}
size_t MDAL::Dataset3D::maximumVerticalLevelsCount() const
{
return mMaximumVerticalLevelsCount;
}
size_t MDAL::Dataset3D::scalarData( size_t, size_t, double * ) { return 0; }
size_t MDAL::Dataset3D::vectorData( size_t, size_t, double * ) { return 0; }
size_t MDAL::Dataset3D::activeData( size_t, size_t, int * ) { return 0; }
MDAL::DatasetGroup::DatasetGroup( const std::string &driverName,
MDAL::Mesh *parent,
const std::string &uri,
@ -185,17 +225,17 @@ void MDAL::DatasetGroup::stopEditing()
mInEditMode = false;
}
bool MDAL::DatasetGroup::isOnVertices() const
MDAL_DataLocation MDAL::DatasetGroup::dataLocation() const
{
return mIsOnVertices;
return mDataLocation;
}
void MDAL::DatasetGroup::setIsOnVertices( bool isOnVertices )
void MDAL::DatasetGroup::setDataLocation( MDAL_DataLocation dataLocation )
{
// datasets are initialized (e.g. values array, active array) based
// on this property. Do not allow to modify later on.
assert( datasets.empty() );
mIsOnVertices = isOnVertices;
mDataLocation = dataLocation;
}
bool MDAL::DatasetGroup::isScalar() const

70
external/mdal/mdal_data_model.hpp vendored
View File

@ -44,18 +44,35 @@ namespace MDAL
Dataset( DatasetGroup *parent );
virtual ~Dataset();
std::string driverName() const;
size_t valuesCount() const;
//! For DataOnVertices2D or DataOnFaces2D
virtual size_t scalarData( size_t indexStart, size_t count, double *buffer ) = 0;
//! For DataOnVertices2D or DataOnFaces2D
virtual size_t vectorData( size_t indexStart, size_t count, double *buffer ) = 0;
//! For DataOnVertices2D or DataOnFaces2D
virtual size_t activeData( size_t indexStart, size_t count, int *buffer ) = 0;
//! For DataOnVolumes3D
virtual size_t verticalLevelCountData( size_t indexStart, size_t count, int *buffer ) = 0;
//! For DataOnVolumes3D
virtual size_t verticalLevelData( size_t indexStart, size_t count, double *buffer ) = 0;
//! For DataOnVolumes3D
virtual size_t faceToVolumeData( size_t indexStart, size_t count, int *buffer ) = 0;
//! For DataOnVolumes3D
virtual size_t scalarVolumesData( size_t indexStart, size_t count, double *buffer ) = 0;
//! For DataOnVolumes3D
virtual size_t vectorVolumesData( size_t indexStart, size_t count, double *buffer ) = 0;
//! For DataOnVolumes3D
virtual size_t activeVolumesData( size_t indexStart, size_t count, int *buffer ) = 0;
virtual size_t volumesCount() const = 0;
virtual size_t maximumVerticalLevelsCount() const = 0;
Statistics statistics() const;
void setStatistics( const Statistics &statistics );
bool isValid() const;
void setIsValid( bool isValid );
DatasetGroup *group() const;
Mesh *mesh() const;
@ -70,6 +87,45 @@ namespace MDAL
Statistics mStatistics;
};
class Dataset2D: public Dataset
{
public:
Dataset2D( DatasetGroup *parent );
virtual ~Dataset2D() override;
size_t verticalLevelCountData( size_t indexStart, size_t count, int *buffer ) override;
size_t verticalLevelData( size_t indexStart, size_t count, double *buffer ) override;
size_t faceToVolumeData( size_t indexStart, size_t count, int *buffer ) override;
size_t scalarVolumesData( size_t indexStart, size_t count, double *buffer ) override;
size_t vectorVolumesData( size_t indexStart, size_t count, double *buffer ) override;
size_t activeVolumesData( size_t indexStart, size_t count, int *buffer ) override;
size_t volumesCount() const override;
size_t maximumVerticalLevelsCount() const override;
};
class Dataset3D: public Dataset
{
public:
Dataset3D(
DatasetGroup *parent,
size_t volumes,
size_t maxVerticalLevelCount
);
virtual ~Dataset3D() override;
virtual size_t scalarData( size_t indexStart, size_t count, double *buffer ) override;
virtual size_t vectorData( size_t indexStart, size_t count, double *buffer ) override;
virtual size_t activeData( size_t indexStart, size_t count, int *buffer ) override;
size_t volumesCount() const override;
size_t maximumVerticalLevelsCount() const override;
private:
size_t mVolumesCount = 0;
size_t mMaximumVerticalLevelsCount = 0;
};
typedef std::vector<std::shared_ptr<Dataset>> Datasets;
class DatasetGroup
@ -101,8 +157,8 @@ namespace MDAL
bool isScalar() const;
void setIsScalar( bool isScalar );
bool isOnVertices() const;
void setIsOnVertices( bool isOnVertices );
MDAL_DataLocation dataLocation() const;
void setDataLocation( MDAL_DataLocation dataLocation );
std::string uri() const;
@ -124,7 +180,7 @@ namespace MDAL
const std::string mDriverName;
Mesh *mParent = nullptr;
bool mIsScalar = true;
bool mIsOnVertices = true;
MDAL_DataLocation mDataLocation = MDAL_DataLocation::DataOnVertices2D;
std::string mUri; // file/uri from where it came
Statistics mStatistics;
std::string mReferenceTime;
@ -154,12 +210,14 @@ namespace MDAL
class Mesh
{
public:
//! Constructs 2D Mesh
Mesh( const std::string &driverName,
size_t verticesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
BBox extent,
const std::string &uri );
virtual ~Mesh();
std::string driverName() const;

6
external/mdal/mdal_driver_manager.cpp vendored
View File

@ -26,6 +26,7 @@
#include "frmts/mdal_ugrid.hpp"
#include "frmts/mdal_3di.hpp"
#include "frmts/mdal_sww.hpp"
#include "frmts/mdal_tuflowfv.hpp"
#endif
#if defined HAVE_GDAL && defined HAVE_NETCDF
@ -49,7 +50,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverManager::load( const std::string &meshFi
for ( const auto &driver : mDrivers )
{
if ( ( driver->hasCapability( Capability::ReadMesh ) ) &&
driver->canRead( meshFile ) )
driver->canReadMesh( meshFile ) )
{
std::unique_ptr<Driver> drv( driver->create() );
mesh = drv->load( meshFile, status );
@ -81,7 +82,7 @@ void MDAL::DriverManager::loadDatasets( Mesh *mesh, const std::string &datasetFi
for ( const auto &driver : mDrivers )
{
if ( driver->hasCapability( Capability::ReadDatasets ) &&
driver->canRead( datasetFile ) )
driver->canReadDatasets( datasetFile ) )
{
std::unique_ptr<Driver> drv( driver->create() );
drv->load( datasetFile, mesh, status );
@ -142,6 +143,7 @@ MDAL::DriverManager::DriverManager()
#endif
#ifdef HAVE_NETCDF
mDrivers.push_back( std::make_shared<MDAL::DriverTuflowFV>() );
mDrivers.push_back( std::make_shared<MDAL::Driver3Di>() );
mDrivers.push_back( std::make_shared<MDAL::DriverSWW>() );
mDrivers.push_back( std::make_shared<MDAL::DriverUgrid>() );

24
external/mdal/mdal_memory_data_model.cpp vendored
View File

@ -11,40 +11,40 @@
#include <iterator>
#include "mdal_utils.hpp"
MDAL::MemoryDataset::MemoryDataset( MDAL::DatasetGroup *grp )
: Dataset( grp )
MDAL::MemoryDataset2D::MemoryDataset2D( MDAL::DatasetGroup *grp )
: Dataset2D( grp )
, mValues( group()->isScalar() ? valuesCount() : 2 * valuesCount(),
std::numeric_limits<double>::quiet_NaN() )
{
if ( group()->isOnVertices() )
if ( group()->dataLocation() == MDAL_DataLocation::DataOnVertices2D )
mActive = std::vector<int>( mesh()->facesCount(), 1 );
}
MDAL::MemoryDataset::~MemoryDataset() = default;
MDAL::MemoryDataset2D::~MemoryDataset2D() = default;
int *MDAL::MemoryDataset::active()
int *MDAL::MemoryDataset2D::active()
{
return mActive.data();
}
double *MDAL::MemoryDataset::values()
double *MDAL::MemoryDataset2D::values()
{
return mValues.data();
}
const int *MDAL::MemoryDataset::constActive() const
const int *MDAL::MemoryDataset2D::constActive() const
{
return mActive.data();
}
const double *MDAL::MemoryDataset::constValues() const
const double *MDAL::MemoryDataset2D::constValues() const
{
return mValues.data();
}
size_t MDAL::MemoryDataset::activeData( size_t indexStart, size_t count, int *buffer )
size_t MDAL::MemoryDataset2D::activeData( size_t indexStart, size_t count, int *buffer )
{
if ( group()->isOnVertices() )
if ( group()->dataLocation() == MDAL_DataLocation::DataOnVertices2D )
{
size_t nValues = mActive.size();
@ -63,7 +63,7 @@ size_t MDAL::MemoryDataset::activeData( size_t indexStart, size_t count, int *bu
return count;
}
size_t MDAL::MemoryDataset::scalarData( size_t indexStart, size_t count, double *buffer )
size_t MDAL::MemoryDataset2D::scalarData( size_t indexStart, size_t count, double *buffer )
{
assert( group()->isScalar() ); //checked in C API interface
size_t nValues = valuesCount();
@ -77,7 +77,7 @@ size_t MDAL::MemoryDataset::scalarData( size_t indexStart, size_t count, double
return copyValues;
}
size_t MDAL::MemoryDataset::vectorData( size_t indexStart, size_t count, double *buffer )
size_t MDAL::MemoryDataset2D::vectorData( size_t indexStart, size_t count, double *buffer )
{
assert( !group()->isScalar() ); //checked in C API interface
size_t nValues = valuesCount();

7
external/mdal/mdal_memory_data_model.hpp vendored
View File

@ -31,11 +31,11 @@ namespace MDAL
/**
* The MemoryDataset stores all the data in the memory
*/
class MemoryDataset: public Dataset
class MemoryDataset2D: public Dataset2D
{
public:
MemoryDataset( DatasetGroup *grp );
~MemoryDataset() override;
MemoryDataset2D( DatasetGroup *grp );
~MemoryDataset2D() override;
size_t scalarData( size_t indexStart, size_t count, double *buffer ) override;
size_t vectorData( size_t indexStart, size_t count, double *buffer ) override;
@ -65,6 +65,7 @@ namespace MDAL
* - vertex count * 2 if isOnVertices & isVector
*/
std::vector<double> mValues;
/**
* Active flag, whether the face is active or not (disabled)
* Only make sense for dataset defined on vertices with size == face count

84
external/mdal/mdal_utils.cpp vendored
View File

@ -14,6 +14,7 @@
#include <cmath>
#include <string.h>
#include <stdio.h>
#include <ctime>
bool MDAL::fileExists( const std::string &filename )
{
@ -333,30 +334,50 @@ double MDAL::parseTimeUnits( const std::string &units )
{
double divBy = 1;
// We are trying to parse strings like
//
// "seconds since 2001-05-05 00:00:00"
// "hours since 1900-01-01 00:00:0.0"
// "days since 1961-01-01 00:00:00"
//
// or simply
// hours, days, seconds, ...
const std::vector<std::string> units_list = MDAL::split( units, " since " );
if ( units_list.size() == 2 )
std::string unit_definition = units;
if ( !units_list.empty() )
{
// Give me hours
if ( units_list[0] == "seconds" )
{
divBy = 3600.0;
}
else if ( units_list[0] == "minutes" )
{
divBy = 60.0;
}
else if ( units_list[0] == "days" )
{
divBy = 1.0 / 24.0;
}
unit_definition = units_list[0];
}
// Give me hours
if ( units_list[0] == "seconds" )
{
divBy = 3600.0;
}
else if ( units_list[0] == "minutes" )
{
divBy = 60.0;
}
else if ( units_list[0] == "days" )
{
divBy = 1.0 / 24.0;
}
return divBy;
}
std::string MDAL::getCurrentTimeStamp()
{
time_t t ;
struct tm *tmp ;
char MY_TIME[50];
time( &t );
tmp = localtime( &t );
strftime( MY_TIME, sizeof( MY_TIME ), "%Y-%m-%dT%H:%M:%S%z", tmp );
std::string s = MDAL::trim( MY_TIME );
return s;
}
MDAL::Statistics _calculateStatistics( const std::vector<double> &values, size_t count, bool isVector )
{
MDAL::Statistics ret;
@ -434,6 +455,7 @@ MDAL::Statistics MDAL::calculateStatistics( std::shared_ptr<Dataset> dataset )
return ret;
bool isVector = !dataset->group()->isScalar();
bool is3D = dataset->group()->dataLocation() == MDAL_DataLocation::DataOnVolumes3D;
size_t bufLen = 2000;
std::vector<double> buffer( isVector ? bufLen * 2 : bufLen );
@ -441,13 +463,27 @@ MDAL::Statistics MDAL::calculateStatistics( std::shared_ptr<Dataset> dataset )
while ( i < dataset->valuesCount() )
{
size_t valsRead;
if ( isVector )
if ( is3D )
{
valsRead = dataset->vectorData( i, bufLen, buffer.data() );
if ( isVector )
{
valsRead = dataset->vectorVolumesData( i, bufLen, buffer.data() );
}
else
{
valsRead = dataset->scalarVolumesData( i, bufLen, buffer.data() );
}
}
else
{
valsRead = dataset->scalarData( i, bufLen, buffer.data() );
if ( isVector )
{
valsRead = dataset->vectorData( i, bufLen, buffer.data() );
}
else
{
valsRead = dataset->scalarData( i, bufLen, buffer.data() );
}
}
if ( valsRead == 0 )
return ret;
@ -489,10 +525,10 @@ void MDAL::addBedElevationDatasetGroup( MDAL::Mesh *mesh, const Vertices &vertic
mesh->uri(),
"Bed Elevation"
);
group->setIsOnVertices( true );
group->setDataLocation( MDAL_DataLocation::DataOnVertices2D );
group->setIsScalar( true );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
dataset->setTime( 0.0 );
double *vals = dataset->values();
for ( size_t i = 0; i < vertices.size(); ++i )
@ -526,10 +562,10 @@ void MDAL::addFaceScalarDatasetGroup( MDAL::Mesh *mesh,
mesh->uri(),
name
);
group->setIsOnVertices( false );
group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
group->setIsScalar( true );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
dataset->setTime( 0.0 );
memcpy( dataset->values(), values.data(), sizeof( double )*values.size() );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
@ -538,10 +574,10 @@ void MDAL::addFaceScalarDatasetGroup( MDAL::Mesh *mesh,
mesh->datasetGroups.push_back( group );
}
void MDAL::activateFaces( MDAL::MemoryMesh *mesh, std::shared_ptr<MemoryDataset> dataset )
void MDAL::activateFaces( MDAL::MemoryMesh *mesh, std::shared_ptr<MemoryDataset2D> dataset )
{
// only for data on vertices
if ( !dataset->group()->isOnVertices() )
if ( !( dataset->group()->dataLocation() == MDAL_DataLocation::DataOnVertices2D ) )
return;
bool isScalar = dataset->group()->isScalar();
@ -624,5 +660,3 @@ std::string MDAL::doubleToString( double value, int precision )
oss << value;
return oss.str();
}

4
external/mdal/mdal_utils.hpp vendored
View File

@ -110,6 +110,8 @@ namespace MDAL
// time
//! Returns a delimiter to get time in hours
MDAL_TEST_EXPORT double parseTimeUnits( const std::string &units );
//! Returns current time stamp in YYYY-mm-ddTHH:MM:SSzoneOffset
std::string getCurrentTimeStamp();
// statistics
void combineStatistics( Statistics &main, const Statistics &other );
@ -127,7 +129,7 @@ namespace MDAL
//! Adds altitude dataset group to mesh
void addFaceScalarDatasetGroup( MDAL::Mesh *mesh, const std::vector<double> &values, const std::string &name );
//! Loop through all faces and activate those which has all 4 values on vertices valid
void activateFaces( MDAL::MemoryMesh *mesh, std::shared_ptr<MemoryDataset> dataset );
void activateFaces( MDAL::MemoryMesh *mesh, std::shared_ptr<MemoryDataset2D> dataset );
//! function used to read all of type of value. Option to change the endianness is provided
template<typename T>

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

@ -90,6 +90,7 @@ IF (WITH_INTERNAL_MDAL)
${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_ugrid.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_tuflowfv.cpp
)
SET(MDAL_LIB_HDRS ${MDAL_LIB_HDRS}
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_cf.hpp
@ -97,6 +98,7 @@ IF (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_netcdf.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_sww.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ugrid.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_tuflowfv.hpp
)
SET (HAVE_NETCDF TRUE)
ENDIF(NETCDF_FOUND)