[MDAL] update MDAL to 0.1.0 (new API)

This commit is contained in:
Peter Petrik 2018-12-04 17:28:05 +01:00
parent d6701f24eb
commit d43f6376eb
23 changed files with 1958 additions and 729 deletions

View File

@ -63,8 +63,9 @@ enum MDAL_Status
Warn_NodeNotUnique
};
//! Mesh
typedef void *MeshH;
typedef void *MeshVertexIteratorH;
typedef void *MeshFaceIteratorH;
typedef void *DatasetGroupH;
typedef void *DatasetH;
@ -78,7 +79,10 @@ MDAL_EXPORT MDAL_Status MDAL_LastStatus();
/// MESH
///////////////////////////////////////////////////////////////////////////////////////
//! Loads mesh file. On error see MDAL_LastStatus for error type This effectively loads whole mesh in-memory
/**
* Loads mesh file. On error see MDAL_LastStatus for error type
* This may effectively load whole mesh in-memory for some providers
*/
MDAL_EXPORT MeshH MDAL_LoadMesh( const char *meshFile );
//! Closes mesh, frees the memory
MDAL_EXPORT void MDAL_CloseMesh( MeshH mesh );
@ -88,20 +92,18 @@ MDAL_EXPORT void MDAL_CloseMesh( MeshH mesh );
* not thread-safe and valid only till next call
*/
MDAL_EXPORT const char *MDAL_M_projection( MeshH mesh );
/**
* Returns mesh extent in native projection
* Returns NaN on error
*/
MDAL_EXPORT void MDAL_M_extent( MeshH mesh, double *minX, double *maxX, double *minY, double *maxY );
//! Returns vertex count for the mesh
MDAL_EXPORT int MDAL_M_vertexCount( MeshH mesh );
//! Returns vertex X coord for the mesh
MDAL_EXPORT double MDAL_M_vertexXCoordinatesAt( MeshH mesh, int index );
//! Returns vertex Y coord for the mesh
MDAL_EXPORT double MDAL_M_vertexYCoordinatesAt( MeshH mesh, int index );
//! Returns vertex Z coord for the mesh
MDAL_EXPORT double MDAL_M_vertexZCoordinatesAt( MeshH mesh, int index );
//! Returns face count for the mesh
MDAL_EXPORT int MDAL_M_faceCount( MeshH mesh );
//! Returns number of vertices face consist of, e.g. 3 for triangle
MDAL_EXPORT int MDAL_M_faceVerticesCountAt( MeshH mesh, int index );
//! Returns vertex index for face
MDAL_EXPORT int MDAL_M_faceVerticesIndexAt( MeshH mesh, int face_index, int vertex_index );
//! Returns maximum number of vertices face can consist of, e.g. 4 for regular quad mesh
MDAL_EXPORT int MDAL_M_faceVerticesMaximumCount( MeshH mesh );
/**
* Loads dataset file. On error see MDAL_LastStatus for error type.
@ -110,17 +112,74 @@ MDAL_EXPORT int MDAL_M_faceVerticesIndexAt( MeshH mesh, int face_index, int vert
* can be freed manually with MDAL_CloseDataset if needed
*/
MDAL_EXPORT void MDAL_M_LoadDatasets( MeshH mesh, const char *datasetFile );
//! Returns dataset groups count
MDAL_EXPORT int MDAL_M_datasetGroupCount( MeshH mesh );
//! Returns dataset group handle
MDAL_EXPORT DatasetGroupH MDAL_M_datasetGroup( MeshH mesh, int index );
///////////////////////////////////////////////////////////////////////////////////////
/// MESH VERTICES
///////////////////////////////////////////////////////////////////////////////////////
/**
* Returns iterator to the mesh vertices
* For some formats this may effectively load all vertices in-memory until iterator is closed
*/
MDAL_EXPORT MeshVertexIteratorH MDAL_M_vertexIterator( MeshH mesh );
/**
* Returns vertices from iterator for the mesh
* \param iterator mesh data iterator
* \param verticesCount maximum number or vertices to be written to buffer
* \param coordinates must be allocated to 3* verticesCount items to store x1, y1, z1, ..., xN, yN, zN coordinates
* \returns number of vertices written in the buffer
*/
MDAL_EXPORT int MDAL_VI_next( MeshVertexIteratorH iterator, int verticesCount, double *coordinates );
//! Closes mesh data iterator, frees the memory
MDAL_EXPORT void MDAL_VI_close( MeshVertexIteratorH iterator );
///////////////////////////////////////////////////////////////////////////////////////
/// MESH FACES
///////////////////////////////////////////////////////////////////////////////////////
/**
* Returns iterator to the mesh faces
* For some formats this may effectively load all faces in-memory until iterator is closed
*/
MDAL_EXPORT MeshFaceIteratorH MDAL_M_faceIterator( MeshH mesh );
/**
* Returns next faces from iterator for the mesh
*
* Reading stops when vertexIndicesBuffer capacity is full / faceOffsetsBuffer
* capacity is full / end of faces is reached, whatever comes first
*
* \param iterator mesh data iterator
* \param faceOffsetsBufferLen size of faceOffsetsBuffer, minimum 1
* \param faceOffsetsBuffer allocated array to store face offset in vertexIndicesBuffer for given face.
* To find number of vertices of face i, calculate faceOffsetsBuffer[i] - faceOffsetsBuffer[i-1]
* \param vertexIndicesBufferLen size of vertexIndicesBuffer, minimum is MDAL_M_faceVerticesMaximumCount()
* \param vertexIndicesBuffer writes vertex indexes for faces
* faceOffsetsBuffer[i-1] is index where the vertices for face i begins,
* \returns number of faces written in the buffer
*/
MDAL_EXPORT int MDAL_FI_next( MeshFaceIteratorH iterator,
int faceOffsetsBufferLen,
int *faceOffsetsBuffer,
int vertexIndicesBufferLen,
int *vertexIndicesBuffer );
//! Closes mesh data iterator, frees the memory
MDAL_EXPORT void MDAL_FI_close( MeshFaceIteratorH iterator );
///////////////////////////////////////////////////////////////////////////////////////
/// DATASET GROUPS
///////////////////////////////////////////////////////////////////////////////////////
//! Returns dataset parent mesh
MDAL_EXPORT MeshH MDAL_G_mesh( DatasetGroupH group );
//! Returns dataset count in group
MDAL_EXPORT int MDAL_G_datasetCount( DatasetGroupH group );
@ -154,6 +213,12 @@ MDAL_EXPORT bool MDAL_G_hasScalarData( DatasetGroupH group );
//! Whether dataset is on vertices
MDAL_EXPORT bool MDAL_G_isOnVertices( DatasetGroupH group );
/**
* Returns the min and max values of the group
* Returns NaN on error
*/
MDAL_EXPORT void MDAL_G_minimumMaximum( DatasetGroupH group, double *min, double *max );
///////////////////////////////////////////////////////////////////////////////////////
/// DATASETS
///////////////////////////////////////////////////////////////////////////////////////
@ -167,33 +232,40 @@ MDAL_EXPORT double MDAL_D_time( DatasetH dataset );
//! Returns number of values
MDAL_EXPORT int MDAL_D_valueCount( DatasetH dataset );
/**
* Returns scalar value associated with the index from the dataset
* for nodata return numeric_limits<double>:quiet_NaN
*/
MDAL_EXPORT double MDAL_D_value( DatasetH dataset, int valueIndex );
/**
* Returns X value associated with the index from the vector dataset
* for nodata return numeric_limits<double>:quiet_NaN
*/
MDAL_EXPORT double MDAL_D_valueX( DatasetH dataset, int valueIndex );
/**
* Returns Y value associated with the index from the vector dataset
* for nodata return numeric_limits<double>:quiet_NaN
*/
MDAL_EXPORT double MDAL_D_valueY( DatasetH dataset, int valueIndex );
/**
* Whether element is active - should be taken into account
* Some formats support switching off the element for particular timestep
*/
MDAL_EXPORT bool MDAL_D_active( DatasetH dataset, int faceIndex );
//! Returns whether dataset is valid
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.
};
/**
* Populates buffer with values from the dataset
* for nodata, returned is numeric_limits<double>::quiet_NaN
*
* \param dataset handle to dataset
* \param indexStart index of face/vertex to start reading of values to the buffer
* \param count number of values to be written to the buffer
* \param dataType type of values to be written to the buffer
* \param buffer output array to be populated with the values. must be already allocated
* For SCALAR_DOUBLE, the minimum size must be valuesCount * size_of(double)
* 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)
* \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 );
/**
* Returns the min and max values of the dataset
* Returns NaN on error
*/
MDAL_EXPORT void MDAL_D_minimumMaximum( DatasetH dataset, double *min, double *max );
#ifdef __cplusplus
}
#endif

View File

@ -17,6 +17,46 @@
#include "mdal.h"
#include "mdal_utils.hpp"
MDAL::Mesh2dm::Mesh2dm( size_t verticesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
MDAL::BBox extent,
const std::string &uri,
const std::map<size_t, size_t> vertexIDtoIndex )
: MemoryMesh( verticesCount, facesCount, faceVerticesMaximumCount, extent, uri )
, mVertexIDtoIndex( vertexIDtoIndex )
{
}
MDAL::Mesh2dm::~Mesh2dm() = default;
bool _parse_vertex_id_gaps( std::map<size_t, size_t> &vertexIDtoIndex, size_t vertexIndex, size_t vertexID, MDAL_Status *status )
{
if ( vertexIndex == vertexID )
return false;
std::map<size_t, size_t>::iterator search = vertexIDtoIndex.find( vertexID );
if ( search != vertexIDtoIndex.end() )
{
if ( status ) *status = MDAL_Status::Warn_ElementNotUnique;
return true;
}
vertexIDtoIndex[vertexID] = vertexIndex;
return false;
}
size_t MDAL::Mesh2dm::vertexIndex( size_t vertexID ) const
{
auto ni2i = mVertexIDtoIndex.find( vertexID );
if ( ni2i != mVertexIDtoIndex.end() )
{
return ni2i->second; // convert from ID to index
}
return vertexID;
}
MDAL::Loader2dm::Loader2dm( const std::string &meshFile ):
mMeshFile( meshFile )
{
@ -71,7 +111,6 @@ std::unique_ptr<MDAL::Mesh> MDAL::Loader2dm::load( MDAL_Status *status )
size_t faceIndex = 0;
size_t vertexIndex = 0;
std::map<size_t, size_t> faceIDtoIndex;
std::map<size_t, size_t> vertexIDtoIndex;
while ( std::getline( in, line ) )
@ -81,20 +120,11 @@ std::unique_ptr<MDAL::Mesh> MDAL::Loader2dm::load( MDAL_Status *status )
chunks = split( line, " ", SplitBehaviour::SkipEmptyParts );
assert( faceIndex < faceCount );
size_t elemID = toSizeT( chunks[1] );
std::map<size_t, size_t>::iterator search = faceIDtoIndex.find( elemID );
if ( search != faceIDtoIndex.end() )
{
if ( status ) *status = MDAL_Status::Warn_ElementNotUnique;
continue;
}
faceIDtoIndex[elemID] = faceIndex;
Face &face = faces[faceIndex];
face.resize( 4 );
// Right now we just store node IDs here - we will convert them to node indices afterwards
for ( size_t i = 0; i < 4; ++i )
face[i] = toSizeT( chunks[i + 2] );
face[i] = toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1
faceIndex++;
}
@ -103,21 +133,12 @@ std::unique_ptr<MDAL::Mesh> MDAL::Loader2dm::load( MDAL_Status *status )
chunks = split( line, " ", SplitBehaviour::SkipEmptyParts );
assert( faceIndex < faceCount );
size_t elemID = toSizeT( chunks[1] );
std::map<size_t, size_t>::iterator search = faceIDtoIndex.find( elemID );
if ( search != faceIDtoIndex.end() )
{
if ( status ) *status = MDAL_Status::Warn_ElementNotUnique;
continue;
}
faceIDtoIndex[elemID] = faceIndex;
Face &face = faces[faceIndex];
face.resize( 3 );
// Right now we just store node IDs here - we will convert them to node indices afterwards
for ( size_t i = 0; i < 3; ++i )
{
face[i] = toSizeT( chunks[i + 2] );
face[i] = toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1
}
faceIndex++;
@ -132,15 +153,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::Loader2dm::load( MDAL_Status *status )
chunks = split( line, " ", SplitBehaviour::SkipEmptyParts );
assert( faceIndex < faceCount );
size_t elemID = toSizeT( chunks[1] );
std::map<size_t, size_t>::iterator search = faceIDtoIndex.find( elemID );
if ( search != faceIDtoIndex.end() )
{
if ( status ) *status = MDAL_Status::Warn_ElementNotUnique;
continue;
}
faceIDtoIndex[elemID] = faceIndex;
//size_t elemID = toSizeT( chunks[1] );
assert( false ); //TODO mark element as unusable
faceIndex++;
@ -148,15 +161,8 @@ std::unique_ptr<MDAL::Mesh> MDAL::Loader2dm::load( MDAL_Status *status )
else if ( startsWith( line, "ND" ) )
{
chunks = split( line, " ", SplitBehaviour::SkipEmptyParts );
size_t nodeID = toSizeT( chunks[1] );
std::map<size_t, size_t>::iterator search = vertexIDtoIndex.find( nodeID );
if ( search != vertexIDtoIndex.end() )
{
if ( status ) *status = MDAL_Status::Warn_NodeNotUnique;
continue;
}
vertexIDtoIndex[nodeID] = vertexIndex;
size_t nodeID = toSizeT( chunks[1] ) - 1; // 2dm is numbered from 1
_parse_vertex_id_gaps( vertexIDtoIndex, vertexIndex, nodeID, status );
assert( vertexIndex < vertexCount );
Vertex &vertex = vertices[vertexIndex];
vertex.x = toDouble( chunks[2] );
@ -178,25 +184,28 @@ std::unique_ptr<MDAL::Mesh> MDAL::Loader2dm::load( MDAL_Status *status )
{
face[nd] = ni2i->second; // convert from ID to index
}
else
else if ( vertices.size() < nodeID )
{
assert( false ); //TODO mark element as unusable
if ( status ) *status = MDAL_Status::Warn_ElementWithInvalidNode;
}
}
//TODO check validity of the face
//check that we have distinct nodes
}
std::unique_ptr< Mesh > mesh( new Mesh );
mesh->uri = mMeshFile;
std::unique_ptr< MemoryMesh > mesh(
new Mesh2dm(
vertices.size(),
faces.size(),
4, //maximum quads
computeExtent( vertices ),
mMeshFile,
vertexIDtoIndex
)
);
mesh->faces = faces;
mesh->vertices = vertices;
mesh->faceIDtoIndex = faceIDtoIndex;
mesh->vertexIDtoIndex = vertexIDtoIndex;
mesh->addBedElevationDataset();
mesh->addBedElevationDataset( vertices, faces );
return mesh;
}

View File

@ -10,10 +10,36 @@
#include <memory>
#include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp"
#include "mdal.h"
namespace MDAL
{
class Mesh2dm: public MemoryMesh
{
public:
Mesh2dm( size_t verticesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
BBox extent,
const std::string &uri,
const std::map<size_t, size_t> vertexIDtoIndex
);
~Mesh2dm() override;
//! Some formats supports gaps in the vertex indexing, but we return continuos array from MDAL
//! for most of the formats this returns
//! \param vertexID internal index/ID of the vertex that native format uses
//! \returns index of the vertex in the continuous array of vertices we returned by readVertices()
virtual size_t vertexIndex( size_t vertexID ) const;
private:
// 2dm supports "gaps" in the mesh indexing
// Store only the indices that have different index and ID
// https://github.com/lutraconsulting/MDAL/issues/51
std::map<size_t, size_t> mVertexIDtoIndex;
};
class Loader2dm
{

View File

@ -33,11 +33,11 @@ MDAL::CFDimensions MDAL::Loader3Di::populateDimensions()
return dims;
}
void MDAL::Loader3Di::populateFacesAndVertices( MDAL::Mesh *mesh )
void MDAL::Loader3Di::populateFacesAndVertices( Vertices &vertices, Faces &faces )
{
assert( mesh );
assert( vertices.empty() );
size_t faceCount = mDimensions.size( CFDimensions::Face2D );
mesh->faces.resize( faceCount );
faces.resize( faceCount );
size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
size_t arrsize = faceCount * verticesInFace;
std::map<std::string, size_t> xyToVertex2DId;
@ -79,9 +79,9 @@ void MDAL::Loader3Di::populateFacesAndVertices( MDAL::Mesh *mesh )
if ( it == xyToVertex2DId.end() )
{
// new vertex
vertexId = mesh->vertices.size();
vertexId = vertices.size();
xyToVertex2DId[key] = vertexId;
mesh->vertices.push_back( vertex );
vertices.push_back( vertex );
}
else
{
@ -93,21 +93,21 @@ void MDAL::Loader3Di::populateFacesAndVertices( MDAL::Mesh *mesh )
}
mesh->faces[faceId] = face;
faces[faceId] = face;
}
// Only now we have number of vertices, since we identified vertices that
// are used in multiple faces
mDimensions.setDimension( CFDimensions::Vertex2D, mesh->vertices.size() );
mDimensions.setDimension( CFDimensions::Vertex2D, vertices.size() );
}
void MDAL::Loader3Di::addBedElevation( MDAL::Mesh *mesh )
{
assert( mesh );
if ( mesh->faces.empty() )
if ( 0 == mesh->facesCount() )
return;
size_t faceCount = mesh->faces.size();
size_t faceCount = mesh->facesCount();
// read Z coordinate of 3di computation nodes centers
int ncidZ = mNcFile.getVarId( "Mesh2DFace_zcc" );
@ -117,20 +117,24 @@ void MDAL::Loader3Di::addBedElevation( MDAL::Mesh *mesh )
return; //error reading the array
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >();
group->isOnVertices = false;
group->isScalar = true;
group->setName( "Bed Elevation" );
group->uri = mesh->uri;
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< Dataset >();
dataset->time = 0.0;
dataset->values.resize( faceCount );
dataset->active.resize( faceCount );
dataset->parent = group.get();
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
mesh,
mesh->uri(),
"Bed Elevation"
);
group->setIsOnVertices( false );
group->setIsScalar( true );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
dataset->setTime( 0.0 );
double *values = dataset->values();
for ( size_t i = 0; i < faceCount; ++i )
{
dataset->values[i].x = MDAL::safeValue( coordZ[i], fillZ );
values[i] = MDAL::safeValue( coordZ[i], fillZ );
}
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->setStatistics( MDAL::calculateStatistics( group ) );
group->datasets.push_back( dataset );
mesh->datasetGroups.push_back( group );
}

View File

@ -43,8 +43,8 @@ namespace MDAL
private:
CFDimensions populateDimensions() override;
void populateFacesAndVertices( MDAL::Mesh *mesh ) override;
void addBedElevation( MDAL::Mesh *mesh ) override;
void populateFacesAndVertices( Vertices &vertices, Faces &faces ) override;
void addBedElevation( Mesh *mesh ) override;
std::string getCoordinateSystemVariableName() override;
std::set<std::string> ignoreNetCDFVariables() override;
std::string nameSuffix( CFDimensions::Type type ) override;

View File

@ -17,6 +17,7 @@
#include "mdal_ascii_dat.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
#include "mdal_2dm.hpp"
#include <math.h>
@ -69,10 +70,12 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
oldFormat = true;
isVector = ( line == "VECTOR" );
group = std::make_shared< DatasetGroup >();
group->uri = mDatFile;
group->setName( name );
group->isScalar = !isVector;
group = std::make_shared< DatasetGroup >(
mesh,
mDatFile,
name
);
group->setIsScalar( !isVector );
}
else
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
@ -83,10 +86,14 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
faceCentered = true;
if ( group )
group->isOnVertices = !faceCentered;
group->setIsOnVertices( !faceCentered );
while ( std::getline( in, line ) )
{
// Replace tabs by spaces,
// since basement v.2.8 uses tabs instead of spaces (e.g. 'TS 0\t0.0')
line = replace( line, "\t", " " );
std::vector<std::string> items = split( line, " ", SplitBehaviour::SkipEmptyParts );
if ( items.size() < 1 )
continue; // empty line?? let's skip it
@ -95,13 +102,13 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
if ( cardType == "ND" && items.size() >= 2 )
{
size_t fileNodeCount = toSizeT( items[1] );
if ( mesh->vertexIDtoIndex.size() != fileNodeCount )
if ( mesh->verticesCount() != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( !oldFormat && cardType == "NC" && items.size() >= 2 )
{
size_t fileElemCount = toSizeT( items[1] );
if ( mesh->faceIDtoIndex.size() != fileElemCount )
if ( mesh->facesCount() != fileElemCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( !oldFormat && cardType == "OBJTYPE" )
@ -118,11 +125,13 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
}
isVector = cardType == "BEGVEC";
group = std::make_shared< DatasetGroup >();
group->uri = mDatFile;
group->setName( name );
group->isScalar = !isVector;
group->isOnVertices = !faceCentered;
group = std::make_shared< DatasetGroup >(
mesh,
mDatFile,
name
);
group->setIsScalar( !isVector );
group->setIsOnVertices( !faceCentered );
}
else if ( !oldFormat && cardType == "ENDDS" )
{
@ -131,6 +140,7 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
debug( "ENDDS card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
group.reset();
}
@ -179,6 +189,7 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
if ( !group || group->datasets.size() == 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
group.reset();
}
@ -194,15 +205,13 @@ void MDAL::LoaderAsciiDat::readVertexTimestep(
{
assert( group );
size_t vertexCount = mesh->vertices.size();
size_t faceCount = mesh->faces.size();
size_t vertexCount = mesh->verticesCount();
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< MDAL::Dataset >();
dataset->time = t / 3600.; // TODO read TIMEUNITS
dataset->values.resize( vertexCount );
dataset->active.resize( faceCount );
dataset->parent = group.get();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
dataset->setTime( t / 3600. ); // TODO read TIMEUNITS
int *active = dataset->active();
// only for new format
for ( size_t i = 0; i < faceCount; ++i )
{
@ -210,49 +219,48 @@ void MDAL::LoaderAsciiDat::readVertexTimestep(
{
std::string line;
std::getline( stream, line );
dataset->active[i] = toBool( line );
active[i] = toBool( line );
}
else
dataset->active[i] = true;
}
for ( size_t i = 0; i < mesh->vertexIDtoIndex.size(); ++i )
const Mesh2dm *m2dm = dynamic_cast<const Mesh2dm *>( mesh );
double *values = dataset->values();
for ( size_t i = 0; i < mesh->verticesCount(); ++i )
{
std::string line;
std::getline( stream, line );
std::vector<std::string> tsItems = split( line, " ", SplitBehaviour::SkipEmptyParts );
auto idx = mesh->vertexIDtoIndex.find( i + 1 ); // ID in 2dm are numbered from 1
if ( idx == mesh->vertexIDtoIndex.end() )
continue; // node ID that does not exist in the mesh
size_t index = idx->second;
size_t index;
if ( m2dm )
index = m2dm->vertexIndex( i );
else
index = i;
if ( isVector )
{
if ( tsItems.size() >= 2 ) // BASEMENT files with vectors have 3 columns
{
dataset->values[index].x = toDouble( tsItems[0] );
dataset->values[index].y = toDouble( tsItems[1] );
values[2 * index] = toDouble( tsItems[0] );
values[2 * index + 1] = toDouble( tsItems[1] );
}
else
{
debug( "invalid timestep line" );
dataset->values[index].noData = true;
}
}
else
{
if ( tsItems.size() >= 1 )
dataset->values[index].x = toDouble( tsItems[0] );
values[index] = toDouble( tsItems[0] );
else
{
debug( "invalid timestep line" );
dataset->values[index].noData = true;
}
}
}
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
@ -265,13 +273,11 @@ void MDAL::LoaderAsciiDat::readFaceTimestep(
{
assert( group );
size_t faceCount = mesh->faces.size();
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< MDAL::Dataset >();
dataset->time = t / 3600.;
dataset->values.resize( faceCount );
dataset->parent = group.get();
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
dataset->setTime( t / 3600. );
double *values = dataset->values();
// TODO: hasStatus
for ( size_t index = 0; index < faceCount; ++index )
{
@ -283,26 +289,25 @@ void MDAL::LoaderAsciiDat::readFaceTimestep(
{
if ( tsItems.size() >= 2 ) // BASEMENT files with vectors have 3 columns
{
dataset->values[index].x = toDouble( tsItems[0] );
dataset->values[index].y = toDouble( tsItems[1] );
values[2 * index] = toDouble( tsItems[0] );
values[2 * index + 1] = toDouble( tsItems[1] );
}
else
{
debug( "invalid timestep line" );
dataset->values[index].noData = true;
}
}
else
{
if ( tsItems.size() >= 1 )
dataset->values[index].x = toDouble( tsItems[0] );
values[index] = toDouble( tsItems[0] );
else
{
debug( "invalid timestep line" );
dataset->values[index].noData = true;
}
}
}
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}

View File

@ -103,8 +103,8 @@ void MDAL::LoaderBinaryDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
if ( !in )
EXIT_WITH_ERROR( MDAL_Status::Err_FileNotFound ); // Couldn't open the file
size_t vertexCount = mesh->vertices.size();
size_t elemCount = mesh->faces.size();
size_t vertexCount = mesh->verticesCount();
size_t elemCount = mesh->facesCount();
int card = 0;
int version;
@ -125,15 +125,19 @@ void MDAL::LoaderBinaryDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
if ( version != CT_VERSION ) // Version should be 3000
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(); // DAT datasets
group->uri = mDatFile;
group->isOnVertices = true;
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
mesh,
mDatFile
); // DAT datasets
group->setIsOnVertices( true );
// in TUFLOW results there could be also a special timestep (99999) with maximums
// we will put it into a separate dataset
std::shared_ptr<DatasetGroup> groupMax = std::make_shared< DatasetGroup >();
groupMax->uri = mDatFile;
groupMax->isOnVertices = true;
std::shared_ptr<DatasetGroup> groupMax = std::make_shared< DatasetGroup >(
mesh,
mDatFile
);
groupMax->setIsOnVertices( true );
while ( card != CT_ENDDS )
{
@ -166,13 +170,13 @@ void MDAL::LoaderBinaryDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
break;
case CT_BEGSCL:
group->isScalar = true;
groupMax->isScalar = true;
group->setIsScalar( true );
groupMax->setIsScalar( true );
break;
case CT_BEGVEC:
group->isScalar = false;
groupMax->isScalar = false;
group->setIsScalar( false );
groupMax->setIsScalar( false );
break;
case CT_VECTYPE:
@ -231,8 +235,14 @@ void MDAL::LoaderBinaryDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
if ( !group || group->datasets.size() == 0 )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
if ( groupMax && groupMax->datasets.size() > 0 )
{
groupMax->setStatistics( MDAL::calculateStatistics( groupMax ) );
mesh->datasetGroups.push_back( groupMax );
}
}
bool MDAL::LoaderBinaryDat::readVertexTimestep( const MDAL::Mesh *mesh,
@ -243,17 +253,15 @@ bool MDAL::LoaderBinaryDat::readVertexTimestep( const MDAL::Mesh *mesh,
int sflg,
std::ifstream &in )
{
assert( group && groupMax && ( group->isScalar == groupMax->isScalar ) );
bool isScalar = group->isScalar;
assert( group && groupMax && ( group->isScalar() == groupMax->isScalar() ) );
bool isScalar = group->isScalar();
size_t vertexCount = mesh->vertices.size();
size_t faceCount = mesh->faces.size();
size_t vertexCount = mesh->verticesCount();
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< MDAL::Dataset >();
dataset->values.resize( vertexCount );
dataset->active.resize( faceCount );
dataset->parent = group.get();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
int *activeFlags = dataset->active();
bool active = true;
for ( size_t i = 0; i < faceCount; ++i )
{
@ -263,9 +271,10 @@ bool MDAL::LoaderBinaryDat::readVertexTimestep( const MDAL::Mesh *mesh,
return true; //error
}
dataset->active[i] = active;
activeFlags[i] = active;
}
double *values = dataset->values();
for ( size_t i = 0; i < vertexCount; ++i )
{
if ( !isScalar )
@ -277,8 +286,8 @@ bool MDAL::LoaderBinaryDat::readVertexTimestep( const MDAL::Mesh *mesh,
if ( read( in, reinterpret_cast< char * >( &y ), 4 ) )
return true; //error
dataset->values[i].x = static_cast< double >( x );
dataset->values[i].y = static_cast< double >( y );
values[2 * i] = static_cast< double >( x );
values[2 * i + 1] = static_cast< double >( y );
}
else
{
@ -287,18 +296,20 @@ bool MDAL::LoaderBinaryDat::readVertexTimestep( const MDAL::Mesh *mesh,
if ( read( in, reinterpret_cast< char * >( &scalar ), 4 ) )
return true; //error
dataset->values[i].x = static_cast< double >( scalar );
values[i] = static_cast< double >( scalar );
}
}
if ( MDAL::equals( time, 99999.0 ) ) // Special TUFLOW dataset with maximus
{
dataset->time = time;
dataset->setTime( time );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
groupMax->datasets.push_back( dataset );
}
else
{
dataset->time = time; // TODO read TIMEUNITS
dataset->setTime( time ); // TODO read TIMEUNITS
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
return false; //OK

View File

@ -132,30 +132,22 @@ MDAL::cfdataset_info_map MDAL::LoaderCF::parseDatasetGroupInfo()
return dsinfo_map;
}
static void populate_vals( bool is_vector, std::vector<MDAL::Value> &vals, size_t i,
static void populate_vals( bool is_vector, double *vals, size_t i,
const std::vector<double> &vals_x, const std::vector<double> &vals_y,
size_t idx, double fill_val_x, double fill_val_y )
{
vals[i].x = MDAL::safeValue( vals_x[idx], fill_val_x );
if ( is_vector )
{
vals[i].y = MDAL::safeValue( vals_y[idx], fill_val_y );
vals[2 * i] = MDAL::safeValue( vals_x[idx], fill_val_x );
vals[2 * i + 1] = MDAL::safeValue( vals_y[idx], fill_val_y );
}
}
static void populate_nodata( std::vector<MDAL::Value> &vals, size_t from_i, size_t to_i )
{
for ( size_t i = from_i; i < to_i; ++i )
else
{
vals[i].noData = true;
vals[i].x = std::numeric_limits<double>::quiet_NaN();
vals[i].y = std::numeric_limits<double>::quiet_NaN();
vals[i] = MDAL::safeValue( vals_x[idx], fill_val_x );
}
}
std::shared_ptr<MDAL::Dataset> MDAL::LoaderCF::createFace2DDataset( size_t ts, const MDAL::CFDatasetGroupInfo &dsi,
std::shared_ptr<MDAL::Dataset> MDAL::LoaderCF::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 )
{
@ -163,18 +155,13 @@ std::shared_ptr<MDAL::Dataset> MDAL::LoaderCF::createFace2DDataset( size_t ts, c
size_t nFaces2D = mDimensions.size( CFDimensions::Face2D );
size_t nLine1D = mDimensions.size( CFDimensions::Line1D );
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared<MDAL::Dataset>();
dataset->values.resize( mDimensions.faceCount() );
populate_nodata( dataset->values,
0,
nLine1D );
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,
dataset->values(),
nLine1D + i,
vals_x,
vals_y,
@ -194,10 +181,13 @@ void MDAL::LoaderCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
{
const CFDatasetGroupInfo dsi = it.second;
// Create a dataset group
std::shared_ptr<MDAL::DatasetGroup> group = std::make_shared<MDAL::DatasetGroup>();
group->uri = mFileName;
group->setName( dsi.name );
group->isScalar = !dsi.is_vector;
std::shared_ptr<MDAL::DatasetGroup> group = std::make_shared<MDAL::DatasetGroup>(
mesh,
mFileName,
dsi.name
);
group->setIsScalar( !dsi.is_vector );
group->setIsOnVertices( false );
// read X data
double fill_val_x = mNcFile.getFillValue( dsi.ncid_x );
@ -221,18 +211,20 @@ void MDAL::LoaderCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
if ( dsi.outputType == CFDimensions::Face2D )
{
group->isOnVertices = false;
dataset = createFace2DDataset( ts, dsi, vals_x, vals_y, fill_val_x, fill_val_y );
dataset = createFace2DDataset( group, ts, dsi, vals_x, vals_y, fill_val_x, fill_val_y );
}
dataset->parent = group.get();
dataset->time = time;
dataset->setTime( time );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
// Add to mesh
if ( !group->datasets.empty() )
{
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
}
}
}
@ -299,9 +291,6 @@ std::unique_ptr< MDAL::Mesh > MDAL::LoaderCF::load( MDAL_Status *status )
{
if ( status ) *status = MDAL_Status::None;
std::unique_ptr< MDAL::Mesh > mesh( new MDAL::Mesh );
mesh->uri = mFileName;
//Dimensions dims;
std::vector<double> times;
@ -314,7 +303,21 @@ std::unique_ptr< MDAL::Mesh > MDAL::LoaderCF::load( MDAL_Status *status )
mDimensions = populateDimensions();
// Create mMesh
populateFacesAndVertices( mesh.get() );
Faces faces;
Vertices vertices;
populateFacesAndVertices( vertices, faces );
std::unique_ptr< MemoryMesh > mesh(
new MemoryMesh(
vertices.size(),
faces.size(),
mDimensions.MaxVerticesInFace,
computeExtent( vertices ),
mFileName
)
);
mesh->faces = faces;
mesh->vertices = vertices;
addBedElevation( mesh.get() );
setProjection( mesh.get() );
@ -326,14 +329,14 @@ std::unique_ptr< MDAL::Mesh > MDAL::LoaderCF::load( MDAL_Status *status )
// Create datasets
addDatasetGroups( mesh.get(), times, dsinfo_map );
return mesh;
}
catch ( MDAL_Status error )
{
if ( status ) *status = ( error );
if ( mesh ) mesh.reset();
return std::unique_ptr<Mesh>();
}
return mesh;
}
//////////////////////////////////////////////////////////////////////////////////////

View File

@ -76,7 +76,7 @@ namespace MDAL
protected:
virtual CFDimensions populateDimensions() = 0;
virtual void populateFacesAndVertices( MDAL::Mesh *mesh ) = 0;
virtual void populateFacesAndVertices( Vertices &vertices, Faces &faces ) = 0;
virtual void addBedElevation( MDAL::Mesh *mesh ) = 0;
virtual std::string getCoordinateSystemVariableName() = 0;
virtual std::set<std::string> ignoreNetCDFVariables() = 0;
@ -87,10 +87,13 @@ namespace MDAL
void setProjection( MDAL::Mesh *m );
cfdataset_info_map parseDatasetGroupInfo();
void parseTime( std::vector<double> &times );
std::shared_ptr<MDAL::Dataset> createFace2DDataset( 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> 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 );
void addDatasetGroups( Mesh *mesh,
const std::vector<double> &times,

View File

@ -258,7 +258,7 @@ void MDAL::LoaderGdal::parseRasterBands( const MDAL::GdalDataset *cfGDALDataset
}
}
void MDAL::LoaderGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<Dataset> tos, bool is_vector, bool is_x )
void MDAL::LoaderGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x )
{
assert( raster_band );
@ -266,6 +266,8 @@ void MDAL::LoaderGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
unsigned int mXSize = meshGDALDataset()->mXSize;
unsigned int mYSize = meshGDALDataset()->mYSize;
double *values = tos->values();
for ( unsigned int y = 0; y < mYSize; ++y )
{
// buffering per-line
@ -292,53 +294,65 @@ void MDAL::LoaderGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
{
unsigned int idx = x + mXSize * y;
double val = mPafScanline[x];
bool noData = false;
if ( MDAL::equals( val, nodata ) )
if ( !MDAL::equals( val, nodata ) )
{
// store all nodata value as this hardcoded number
val = MDAL_NODATA;
noData = true;
}
if ( is_vector )
{
if ( is_x )
// values is prepolulated with NODATA values, so store only legal values
if ( is_vector )
{
tos->values[idx].x = val;
tos->values[idx].noData = noData;
if ( is_x )
{
values[2 * idx] = val;
}
else
{
values[2 * idx + 1] = val;
}
}
else
{
tos->values[idx].y = val;
tos->values[idx].noData = noData;
values[idx] = val;
}
}
else
{
tos->values[idx].x = val;
tos->values[idx].noData = noData;
}
}
}
}
void MDAL::LoaderGdal::activateFaces( std::shared_ptr<Dataset> tos )
void MDAL::LoaderGdal::activateFaces( std::shared_ptr<MemoryDataset> tos )
{
// only for data on vertices
if ( !tos->group()->isOnVertices() )
return;
bool isScalar = tos->group()->isScalar();
// Activate only Faces that do all Vertex's outputs with some data
int *active = tos->active();
const double *values = tos->constValues();
for ( unsigned int idx = 0; idx < meshGDALDataset()->mNVolumes; ++idx )
{
Face elem = mMesh->faces.at( idx );
if ( tos->values[elem[0]].noData ||
tos->values[elem[1]].noData ||
tos->values[elem[2]].noData ||
tos->values[elem[3]].noData )
for ( size_t i = 0; i < 4; ++i )
{
tos->active[idx] = 0; //NOT ACTIVE
}
else
{
tos->active[idx] = 1; //ACTIVE
if ( isScalar )
{
double val = values[elem[i]];
if ( std::isnan( val ) )
{
active[idx] = 0; //NOT ACTIVE
break;
}
}
else
{
double x = values[elem[2 * i]];
double y = values[elem[2 * i + 1]];
if ( std::isnan( x ) || std::isnan( y ) )
{
active[idx] = 0; //NOT ACTIVE
break;
}
}
}
}
}
@ -348,32 +362,35 @@ void MDAL::LoaderGdal::addDatasetGroups()
// Add dataset to mMesh
for ( data_hash::const_iterator band = mBands.begin(); band != mBands.end(); band++ )
{
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >();
group->uri = mFileName;
group->setName( band->first );
group->isOnVertices = true;
if ( band->second.empty() )
continue;
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
mMesh.get(),
mFileName,
band->first
);
group->setIsOnVertices( true );
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;
bool is_vector = ( raster_bands.size() > 1 );
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< MDAL::Dataset >();
group->isScalar = !is_vector;
dataset->time = time_step->first;
dataset->values.resize( meshGDALDataset()->mNPoints );
dataset->active.resize( meshGDALDataset()->mNVolumes );
dataset->parent = group.get();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
dataset->setTime( time_step->first );
for ( std::vector<GDALRasterBandH>::size_type i = 0; i < raster_bands.size(); ++i )
{
addDataToOutput( raster_bands[i], dataset, is_vector, i == 0 );
}
activateFaces( dataset );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
// TODO use GDALComputeRasterMinMax
group->setStatistics( MDAL::calculateStatistics( group ) );
mMesh->datasetGroups.push_back( group );
}
}
@ -386,7 +403,14 @@ void MDAL::LoaderGdal::createMesh()
Faces faces( meshGDALDataset()->mNVolumes );
initFaces( vertices, faces, is_longitude_shifted );
mMesh.reset( new Mesh() );
mMesh.reset( new MemoryMesh(
vertices.size(),
faces.size(),
4, //maximum quads
computeExtent( vertices ),
mFileName
)
);
mMesh->vertices = vertices;
mMesh->faces = faces;
bool proj_added = addSrcProj();

View File

@ -82,9 +82,9 @@ 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<Dataset> tos, bool is_vector, bool is_x );
void addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x );
bool addSrcProj();
void activateFaces( std::shared_ptr<Dataset> tos );
void activateFaces( std::shared_ptr<MemoryDataset> tos );
void addDatasetGroups();
void createMesh();
void parseRasterBands( const GdalDataset *cfGDALDataset );
@ -92,7 +92,7 @@ namespace MDAL
const std::string mFileName;
const std::string mDriverName; /* GDAL driver name */
double *mPafScanline; /* temporary buffer for reading one raster line */
std::unique_ptr< Mesh > mMesh;
std::unique_ptr< MemoryMesh > mMesh;
gdal_datasets_vector gdal_datasets;
data_hash mBands; /* raster bands GDAL handle */
};

241
external/mdal/frmts/mdal_hdf5.cpp vendored Normal file
View File

@ -0,0 +1,241 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2018 Lutra Consulting Limited
*/
#include "mdal_hdf5.hpp"
HdfFile::HdfFile( const std::string &path )
: d( std::make_shared< Handle >( H5Fopen( path.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) ) )
{
}
bool HdfFile::isValid() const { return d->id >= 0; }
hid_t HdfFile::id() const { return d->id; }
HdfGroup::HdfGroup( hid_t file, const std::string &path )
: d( std::make_shared< Handle >( H5Gopen( file, path.c_str() ) ) )
{}
bool HdfGroup::isValid() const { return d->id >= 0; }
hid_t HdfGroup::id() const { return d->id; }
hid_t HdfGroup::file_id() const { return H5Iget_file_id( d->id ); }
std::string HdfGroup::name() const
{
char name[HDF_MAX_NAME];
H5Iget_name( d->id, name, HDF_MAX_NAME );
return std::string( name );
}
std::vector<std::string> HdfGroup::groups() const { return objects( H5G_GROUP ); }
std::vector<std::string> HdfGroup::datasets() const { return objects( H5G_DATASET ); }
std::vector<std::string> HdfGroup::objects() const { return objects( H5G_UNKNOWN ); }
std::string HdfGroup::childPath( const std::string &childName ) const { return name() + "/" + childName; }
std::vector<std::string> HdfGroup::objects( H5G_obj_t type ) const
{
std::vector<std::string> lst;
hsize_t nobj;
H5Gget_num_objs( d->id, &nobj );
for ( hsize_t i = 0; i < nobj; ++i )
{
if ( type == H5G_UNKNOWN || H5Gget_objtype_by_idx( d->id, i ) == type )
{
char name[HDF_MAX_NAME];
H5Gget_objname_by_idx( d->id, i, name, ( size_t )HDF_MAX_NAME );
lst.push_back( std::string( name ) );
}
}
return lst;
}
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 ) ) )
{}
bool HdfAttribute::isValid() const { return d->id >= 0; }
hid_t HdfAttribute::id() const { return d->id; }
std::string HdfAttribute::readString() const
{
char name[HDF_MAX_NAME];
hid_t datatype = H5Tcopy( H5T_C_S1 );
H5Tset_size( datatype, HDF_MAX_NAME );
herr_t status = H5Aread( d->id, datatype, name );
if ( status < 0 )
{
//MDAL::debug("Failed to read data!");
return std::string();
}
H5Tclose( datatype );
return std::string( name );
}
HdfDataset::HdfDataset( hid_t file, const std::string &path )
: d( std::make_shared< Handle >( H5Dopen2( file, path.c_str(), H5P_DEFAULT ) ) )
{}
bool HdfDataset::isValid() const { return d->id >= 0; }
hid_t HdfDataset::id() const { return d->id; }
std::vector<hsize_t> HdfDataset::dims() const
{
hid_t sid = H5Dget_space( d->id );
std::vector<hsize_t> d( H5Sget_simple_extent_ndims( sid ) );
H5Sget_simple_extent_dims( sid, d.data(), NULL );
H5Sclose( sid );
return d;
}
hsize_t HdfDataset::elementCount() const
{
hsize_t count = 1;
for ( hsize_t dsize : dims() )
count *= dsize;
return count;
}
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;
}
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 ); }
std::vector<float> HdfDataset::readArray( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const { return readArray<float>( H5T_NATIVE_FLOAT, offsets, counts ); }
std::vector<double> HdfDataset::readArrayDouble( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const { return readArray<double>( H5T_NATIVE_DOUBLE, offsets, counts ); }
std::vector<int> HdfDataset::readArrayInt( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const { return readArray<int>( H5T_NATIVE_INT, offsets, counts ); }
std::vector<uchar> HdfDataset::readArrayUint8() const { return readArray<uchar>( H5T_NATIVE_UINT8 ); }
std::vector<float> HdfDataset::readArray() const { return readArray<float>( H5T_NATIVE_FLOAT ); }
std::vector<double> HdfDataset::readArrayDouble() const { return readArray<double>( H5T_NATIVE_DOUBLE ); }
std::vector<int> HdfDataset::readArrayInt() const { return readArray<int>( H5T_NATIVE_INT ); }
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 );
for ( const HdfString &str : arr )
{
std::string dat = std::string( str.data );
ret.push_back( MDAL::trim( dat ) );
}
return ret;
}
float HdfDataset::readFloat() const
{
if ( elementCount() != 1 )
{
MDAL::debug( "Not scalar!" );
return 0;
}
float value;
herr_t status = H5Dread( d->id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value );
if ( status < 0 )
{
MDAL::debug( "Failed to read data!" );
return 0;
}
return value;
}
std::string HdfDataset::readString() const
{
if ( elementCount() != 1 )
{
MDAL::debug( "Not scalar!" );
return std::string();
}
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 );
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(
dims.size(),
dims.data(),
dims.data()
)
)
)
{
}
HdfDataspace::HdfDataspace( hid_t dataset )
: d( std::make_shared< Handle >( H5Dget_space( dataset ) ) )
{
}
void HdfDataspace::selectHyperslab( hsize_t start, hsize_t count )
{
// this function works only for 1D arrays
assert( H5Sget_simple_extent_ndims( d->id ) == 1 );
herr_t status = H5Sselect_hyperslab( d->id, H5S_SELECT_SET, &start, NULL, &count, NULL );
if ( status < 0 )
{
MDAL::debug( "Failed to select 1D hyperslab!" );
}
}
void HdfDataspace::selectHyperslab( const std::vector<hsize_t> offsets,
const std::vector<hsize_t> counts )
{
assert( H5Sget_simple_extent_ndims( d->id ) == static_cast<int>( offsets.size() ) );
assert( offsets.size() == counts.size() );
herr_t status = H5Sselect_hyperslab( d->id,
H5S_SELECT_SET,
offsets.data(),
NULL,
counts.data(),
NULL );
if ( status < 0 )
{
MDAL::debug( "Failed to select 1D hyperslab!" );
}
}
bool HdfDataspace::isValid() const { return d->id >= 0; }
hid_t HdfDataspace::id() const { return d->id; }

View File

@ -6,6 +6,7 @@
#ifndef MDAL_HDF5_HPP
#define MDAL_HDF5_HPP
/** A simple C++ wrapper around HDF5 library API */
// for compatibility (older hdf5 version in Travis)
@ -18,6 +19,8 @@ typedef unsigned char uchar;
#include <memory>
#include <vector>
#include <string>
#include <numeric>
#include <assert.h>
#include "stdlib.h"
@ -31,7 +34,8 @@ template <int TYPE> inline void hdfClose( hid_t id ) { MDAL_UNUSED( id ); assert
template <> inline void hdfClose<H5I_FILE>( hid_t id ) { H5Fclose( id ); }
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 ) { 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 <int TYPE>
class HdfH
@ -47,25 +51,24 @@ class HdfH
class HdfGroup;
class HdfDataset;
class HdfAttribute;
class HdfDataspace;
class HdfFile
{
public:
typedef HdfH<H5I_FILE> Handle;
HdfFile( const std::string &path )
: d( std::make_shared< Handle >( H5Fopen( path.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) ) )
{
}
HdfFile( const std::string &path );
bool isValid() const { return d->id >= 0; }
hid_t id() const { return d->id; }
bool isValid() const;
hid_t id() const;
inline std::vector<std::string> groups() const;
inline HdfGroup group( const std::string &path ) const;
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;
protected:
std::shared_ptr<Handle> d;
@ -76,49 +79,27 @@ class HdfGroup
public:
typedef HdfH<H5I_GROUP> Handle;
HdfGroup( hid_t file, const std::string &path )
: d( std::make_shared< Handle >( H5Gopen( file, path.c_str() ) ) )
{}
HdfGroup( hid_t file, const std::string &path );
bool isValid() const { return d->id >= 0; }
hid_t id() const { return d->id; }
hid_t file_id() const { return H5Iget_file_id( d->id ); }
bool isValid() const;
hid_t id() const;
hid_t file_id() const;
std::string name() const
{
char name[HDF_MAX_NAME];
H5Iget_name( d->id, name, HDF_MAX_NAME );
return std::string( name );
}
std::string name() const;
std::vector<std::string> groups() const { return objects( H5G_GROUP ); }
std::vector<std::string> datasets() const { return objects( H5G_DATASET ); }
std::vector<std::string> objects() const { return objects( H5G_UNKNOWN ); }
std::vector<std::string> groups() const;
std::vector<std::string> datasets() const;
std::vector<std::string> objects() const;
std::string childPath( const std::string &childName ) const { return name() + "/" + childName; }
std::string childPath( const std::string &childName ) const;
inline HdfGroup group( const std::string &groupName ) const;
inline HdfDataset dataset( const std::string &dsName ) const;
inline HdfAttribute attribute( const std::string &attr_name ) const;
inline bool pathExists( const std::string &path ) const;
protected:
std::vector<std::string> objects( H5G_obj_t type ) const
{
std::vector<std::string> lst;
hsize_t nobj;
H5Gget_num_objs( d->id, &nobj );
for ( hsize_t i = 0; i < nobj; ++i )
{
if ( type == H5G_UNKNOWN || H5Gget_objtype_by_idx( d->id, i ) == type )
{
char name[HDF_MAX_NAME];
H5Gget_objname_by_idx( d->id, i, name, ( size_t )HDF_MAX_NAME );
lst.push_back( std::string( name ) );
}
}
return lst;
}
std::vector<std::string> objects( H5G_obj_t type ) const;
protected:
std::shared_ptr<Handle> d;
@ -130,27 +111,33 @@ class HdfAttribute
public:
typedef HdfH<H5I_ATTR> Handle;
HdfAttribute( hid_t obj_id, const std::string &attr_name )
: d( std::make_shared< Handle >( H5Aopen( obj_id, attr_name.c_str(), H5P_DEFAULT ) ) )
{}
HdfAttribute( hid_t obj_id, const std::string &attr_name );
bool isValid() const { return d->id >= 0; }
hid_t id() const { return d->id; }
bool isValid() const;
hid_t id() const;
std::string readString() const;
protected:
std::shared_ptr<Handle> d;
};
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 );
//! select from 1D array
void selectHyperslab( hsize_t start, hsize_t count );
//! select from N-D array
void selectHyperslab( const std::vector<hsize_t> offsets,
const std::vector<hsize_t> counts );
bool isValid() const;
hid_t id() const;
std::string readString() const
{
char name[HDF_MAX_NAME];
hid_t datatype = H5Tcopy( H5T_C_S1 );
H5Tset_size( datatype, HDF_MAX_NAME );
herr_t status = H5Aread( d->id, datatype, name );
if ( status < 0 )
{
//MDAL::debug("Failed to read data!");
return std::string();
}
H5Tclose( datatype );
return std::string( name );
}
protected:
std::shared_ptr<Handle> d;
};
@ -160,65 +147,34 @@ class HdfDataset
public:
typedef HdfH<H5I_DATASET> Handle;
HdfDataset( hid_t file, const std::string &path )
: d( std::make_shared< Handle >( H5Dopen2( file, path.c_str(), H5P_DEFAULT ) ) )
{}
HdfDataset( hid_t file, const std::string &path );
bool isValid() const { return d->id >= 0; }
hid_t id() const { return d->id; }
bool isValid() const;
hid_t id() const;
std::vector<hsize_t> dims() const
{
hid_t sid = H5Dget_space( d->id );
std::vector<hsize_t> d( H5Sget_simple_extent_ndims( sid ) );
H5Sget_simple_extent_dims( sid, d.data(), NULL );
H5Sclose( sid );
return d;
}
std::vector<hsize_t> dims() const;
hsize_t elementCount() const
{
hsize_t count = 1;
for ( hsize_t dsize : dims() )
count *= dsize;
return count;
}
hsize_t elementCount() const;
H5T_class_t type() const
{
hid_t tid = H5Dget_type( d->id );
H5T_class_t t_class = H5Tget_class( tid );
H5Tclose( tid );
return t_class;
}
H5T_class_t type() const;
std::vector<uchar> readArrayUint8() const { return readArray<uchar>( H5T_NATIVE_UINT8 ); }
//! Reads full array into vector
//! Array can have any number of dimenstions
//! and it is fully read into 1D vector
std::vector<uchar> readArrayUint8() const;
std::vector<float> readArray() const;
std::vector<double> readArrayDouble() const;
std::vector<int> readArrayInt() const;
std::vector<std::string> readArrayString() const;
std::vector<float> readArray() const { return readArray<float>( H5T_NATIVE_FLOAT ); }
std::vector<double> readArrayDouble() const { return readArray<double>( H5T_NATIVE_DOUBLE ); }
std::vector<int> readArrayInt() const { return readArray<int>( H5T_NATIVE_INT ); }
std::vector<std::string> 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 );
for ( const HdfString &str : arr )
{
std::string dat = std::string( str.data );
ret.push_back( MDAL::trim( dat ) );
}
return ret;
}
//! Reads part of the N-D array into vector,
//! for each dimension specified by offset and count
//! size of offsets and counts must be same as rank (number of dims) of dataset
//! the results array is 1D
std::vector<uchar> readArrayUint8( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const;
std::vector<float> readArray( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const;
std::vector<double> readArrayDouble( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const;
std::vector<int> readArrayInt( const std::vector<hsize_t> offsets, const std::vector<hsize_t> counts ) const;
template <typename T> std::vector<T> readArray( hid_t mem_type_id ) const
@ -234,44 +190,34 @@ class HdfDataset
return data;
}
float readFloat() const
template <typename T> std::vector<T> readArray( hid_t mem_type_id,
const std::vector<hsize_t> offsets,
const std::vector<hsize_t> counts ) const
{
if ( elementCount() != 1 )
{
MDAL::debug( "Not scalar!" );
return 0;
}
HdfDataspace dataspace( d->id );
dataspace.selectHyperslab( offsets, counts );
float value;
herr_t status = H5Dread( d->id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value );
hsize_t totalItems = 1;
for ( auto it = counts.begin(); it != counts.end(); ++it )
totalItems *= *it;
std::vector<hsize_t> dims = {totalItems};
HdfDataspace memspace( dims );
memspace.selectHyperslab( 0, totalItems );
std::vector<T> data( totalItems );
herr_t status = H5Dread( d->id, mem_type_id, memspace.id(), dataspace.id(), H5P_DEFAULT, data.data() );
if ( status < 0 )
{
MDAL::debug( "Failed to read data!" );
return 0;
return std::vector<T>();
}
return value;
return data;
}
std::string readString() const
{
if ( elementCount() != 1 )
{
MDAL::debug( "Not scalar!" );
return std::string();
}
float readFloat() 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 );
if ( status < 0 )
{
MDAL::debug( "Failed to read data!" );
return std::string();
}
H5Tclose( datatype );
return std::string( name );
}
std::string readString() const;
protected:
std::shared_ptr<Handle> d;
@ -291,4 +237,8 @@ inline HdfAttribute HdfFile::attribute( const std::string &attr_name ) const { r
inline HdfAttribute HdfGroup::attribute( const std::string &attr_name ) const { return HdfAttribute( d->id, attr_name ); }
inline bool HdfFile::pathExists( const std::string &path ) const { return H5Lexists( d->id, path.c_str(), H5P_DEFAULT ) > 0; }
inline bool HdfGroup::pathExists( const std::string &path ) const { return H5Lexists( d->id, path.c_str(), H5P_DEFAULT ) > 0; }
#endif // MDAL_HDF5_HPP

View File

@ -12,7 +12,78 @@
#include <string>
#include <vector>
#include <memory>
#include <algorithm>
MDAL::XmdfDataset::~XmdfDataset() = default;
MDAL::XmdfDataset::XmdfDataset( DatasetGroup *grp, const HdfDataset &valuesDs, const HdfDataset &activeDs, hsize_t timeIndex )
: Dataset( grp )
, mHdf5DatasetValues( valuesDs )
, mHdf5DatasetActive( activeDs )
, mTimeIndex( timeIndex )
{
}
const HdfDataset &MDAL::XmdfDataset::dsValues() const
{
return mHdf5DatasetValues;
}
const HdfDataset &MDAL::XmdfDataset::dsActive() const
{
return mHdf5DatasetActive;
}
hsize_t MDAL::XmdfDataset::timeIndex() const
{
return mTimeIndex;
}
size_t MDAL::XmdfDataset::scalarData( size_t indexStart, size_t count, double *buffer )
{
assert( group()->isScalar() ); //checked in C API interface
std::vector<hsize_t> offsets = {timeIndex(), indexStart};
std::vector<hsize_t> counts = {1, count};
std::vector<float> values = dsValues().readArray( offsets, counts );
const float *input = values.data();
for ( size_t j = 0; j < count; ++j )
{
buffer[j] = double( input[j] );
}
return count;
}
size_t MDAL::XmdfDataset::vectorData( size_t indexStart, size_t count, double *buffer )
{
assert( !group()->isScalar() ); //checked in C API interface
std::vector<hsize_t> offsets = {timeIndex(), indexStart, 0};
std::vector<hsize_t> counts = {1, count, 2};
std::vector<float> values = dsValues().readArray( offsets, counts );
const float *input = values.data();
for ( size_t j = 0; j < count; ++j )
{
buffer[2 * j] = double( input[2 * j] );
buffer[2 * j + 1] = double( input[2 * j + 1] );
}
return count;
}
size_t MDAL::XmdfDataset::activeData( size_t indexStart, size_t count, int *buffer )
{
std::vector<hsize_t> offsets = {timeIndex(), indexStart};
std::vector<hsize_t> counts = {1, count};
std::vector<uchar> active = dsActive().readArrayUint8( offsets, counts );
const uchar *input = active.data();
for ( size_t j = 0; j < count; ++j )
{
buffer[j] = bool( input[ j ] );
}
return count;
}
///////////////////////////////////////////////////////////////////////////////////////
MDAL::LoaderXmdf::LoaderXmdf( const std::string &datFile )
: mDatFile( datFile )
@ -20,6 +91,7 @@ MDAL::LoaderXmdf::LoaderXmdf( const std::string &datFile )
void MDAL::LoaderXmdf::load( MDAL::Mesh *mesh, MDAL_Status *status )
{
mMesh = mesh;
if ( status ) *status = MDAL_Status::None;
HdfFile file( mDatFile );
@ -38,8 +110,8 @@ void MDAL::LoaderXmdf::load( MDAL::Mesh *mesh, MDAL_Status *status )
// TODO: check version?
size_t vertexCount = mesh->vertices.size();
size_t faceCount = mesh->faces.size();
size_t vertexCount = mesh->verticesCount();
size_t faceCount = mesh->facesCount();
std::vector<std::string> rootGroups = file.groups();
if ( rootGroups.size() != 1 )
{
@ -53,31 +125,40 @@ void MDAL::LoaderXmdf::load( MDAL::Mesh *mesh, MDAL_Status *status )
DatasetGroups groups; // DAT outputs data
HdfGroup gTemporal = gMesh.group( "Temporal" );
if ( gTemporal.isValid() )
if ( gMesh.pathExists( "Temporal" ) )
{
addDatasetGroupsFromXmdfGroup( groups, gTemporal, vertexCount, faceCount );
HdfGroup gTemporal = gMesh.group( "Temporal" );
if ( gTemporal.isValid() )
{
addDatasetGroupsFromXmdfGroup( groups, gTemporal, vertexCount, faceCount );
}
}
HdfGroup gMaximums = gMesh.group( "Maximums" );
if ( gMaximums.isValid() )
if ( gMesh.pathExists( "Temporal" ) )
{
for ( const std::string &name : gMaximums.groups() )
HdfGroup gMaximums = gMesh.group( "Maximums" );
if ( gMaximums.isValid() )
{
HdfGroup g = gMaximums.group( name );
std::shared_ptr<MDAL::DatasetGroup> maxGroup = readXmdfGroupAsDatasetGroup( g, name + "/Maximums", vertexCount, faceCount );
if ( maxGroup->datasets.size() != 1 )
MDAL::debug( "Maximum dataset should have just one timestep!" );
else
groups.push_back( maxGroup );
for ( const std::string &name : gMaximums.groups() )
{
HdfGroup g = gMaximums.group( name );
std::shared_ptr<MDAL::DatasetGroup> maxGroup = readXmdfGroupAsDatasetGroup( g, name + "/Maximums", vertexCount, faceCount );
if ( !maxGroup || maxGroup->datasets.size() != 1 )
MDAL::debug( "Maximum dataset should have just one timestep!" );
else
groups.push_back( maxGroup );
}
}
}
// res_to_res.exe (TUFLOW utiity tool)
HdfGroup gDifference = gMesh.group( "Difference" );
if ( gDifference.isValid() )
if ( gMesh.pathExists( "Difference" ) )
{
addDatasetGroupsFromXmdfGroup( groups, gDifference, vertexCount, faceCount );
HdfGroup gDifference = gMesh.group( "Difference" );
if ( gDifference.isValid() )
{
addDatasetGroupsFromXmdfGroup( groups, gDifference, vertexCount, faceCount );
}
}
mesh->datasetGroups.insert(
@ -93,18 +174,21 @@ void MDAL::LoaderXmdf::addDatasetGroupsFromXmdfGroup( DatasetGroups &groups, con
{
HdfGroup g = rootGroup.group( name );
std::shared_ptr<DatasetGroup> ds = readXmdfGroupAsDatasetGroup( g, name, vertexCount, faceCount );
groups.push_back( ds );
if ( ds && ds->datasets.size() > 0 )
groups.push_back( ds );
}
}
std::shared_ptr<MDAL::DatasetGroup> MDAL::LoaderXmdf::readXmdfGroupAsDatasetGroup(
const HdfGroup &rootGroup, const std::string &name, size_t vertexCount, size_t faceCount )
{
std::shared_ptr<DatasetGroup> group( new DatasetGroup() );
std::shared_ptr<DatasetGroup> group;
std::vector<std::string> gDataNames = rootGroup.datasets();
if ( !MDAL::contains( gDataNames, "Times" ) ||
!MDAL::contains( gDataNames, "Values" ) ||
!MDAL::contains( gDataNames, "Active" ) )
!MDAL::contains( gDataNames, "Active" ) ||
!MDAL::contains( gDataNames, "Mins" ) ||
!MDAL::contains( gDataNames, "Maxs" ) )
{
MDAL::debug( "ignoring dataset " + name + " - not having required arrays" );
return group;
@ -113,19 +197,31 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::LoaderXmdf::readXmdfGroupAsDatasetGrou
HdfDataset dsTimes = rootGroup.dataset( "Times" );
HdfDataset dsValues = rootGroup.dataset( "Values" );
HdfDataset dsActive = rootGroup.dataset( "Active" );
HdfDataset dsMins = rootGroup.dataset( "Mins" );
HdfDataset dsMaxs = rootGroup.dataset( "Maxs" );
std::vector<hsize_t> dimTimes = dsTimes.dims();
std::vector<hsize_t> dimValues = dsValues.dims();
std::vector<hsize_t> dimActive = dsActive.dims();
std::vector<hsize_t> dimMins = dsMins.dims();
std::vector<hsize_t> dimMaxs = dsMaxs.dims();
if ( dimTimes.size() != 1 || ( dimValues.size() != 2 && dimValues.size() != 3 ) || dimActive.size() != 2 )
if ( dimTimes.size() != 1 ||
( dimValues.size() != 2 && dimValues.size() != 3 ) ||
dimActive.size() != 2 ||
dimMins.size() != 1 ||
dimMaxs.size() != 1
)
{
MDAL::debug( "ignoring dataset " + name + " - arrays not having correct dimension counts" );
return group;
}
hsize_t nTimeSteps = dimTimes[0];
if ( dimValues[0] != nTimeSteps || dimActive[0] != nTimeSteps )
if ( dimValues[0] != nTimeSteps ||
dimActive[0] != nTimeSteps ||
dimMins[0] != nTimeSteps ||
dimMaxs[0] != nTimeSteps )
{
MDAL::debug( "ignoring dataset " + name + " - arrays not having correct dimension sizes" );
return group;
@ -139,44 +235,32 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::LoaderXmdf::readXmdfGroupAsDatasetGrou
bool isVector = dimValues.size() == 3;
std::vector<float> times = dsTimes.readArray();
std::vector<float> values = dsValues.readArray();
std::vector<uchar> active = dsActive.readArrayUint8();
group->setName( name );
group->isScalar = !isVector;
group->isOnVertices = true;
group->uri = mDatFile;
// all fine, set group and return
group = std::make_shared<MDAL::DatasetGroup>(
mMesh,
mDatFile,
name
);
group->setIsScalar( !isVector );
group->setIsOnVertices( true );
// lazy loading of min and max of the dataset group
std::vector<float> mins = dsMins.readArray();
std::vector<float> maxs = dsMaxs.readArray();
Statistics stats;
stats.minimum = static_cast<double>( *std::min_element( mins.begin(), mins.end() ) );
stats.maximum = static_cast<double>( *std::max_element( maxs.begin(), maxs.end() ) );
group->setStatistics( stats );
for ( hsize_t i = 0; i < nTimeSteps; ++i )
{
std::shared_ptr<Dataset> dataset( new Dataset() );
dataset->values.resize( vertexCount );
dataset->active.resize( faceCount );
dataset->parent = group.get();
dataset->time = double( times[i] );
if ( isVector )
{
const float *input = values.data() + 2 * i * vertexCount;
for ( size_t j = 0; j < vertexCount; ++j )
{
dataset->values[j].x = double( input[2 * j] );
dataset->values[j].y = double( input[2 * j + 1] );
}
}
else
{
const float *input = values.data() + i * vertexCount;
for ( size_t j = 0; j < vertexCount; ++j )
{
dataset->values[j].x = double( input[j] );
}
}
const uchar *input = active.data() + i * faceCount;
for ( size_t j = 0; j < faceCount; ++j )
{
dataset->active[j] = input[j];
}
std::shared_ptr<XmdfDataset> dataset( new XmdfDataset( group.get(), dsValues, dsActive, i ) );
dataset->setTime( double( times[i] ) );
Statistics stats;
stats.minimum = static_cast<double>( mins[i] );
stats.maximum = static_cast<double>( maxs[i] );
dataset->setStatistics( stats );
group->datasets.push_back( dataset );
}

View File

@ -20,6 +20,40 @@
namespace MDAL
{
/**
* The XmdfDataset reads the data directly from HDF5 file
* by usage of hyperslabs retrieval
*
* basically all (timesteps) data for one particular dataset groups
* are stored in single
* 3D arrays (time, x, y) for vector datasets
* 2D arrays (time, x) for scalar datasets
* 2D arrays (time, active) for active flags
*/
class XmdfDataset: public Dataset
{
public:
XmdfDataset( DatasetGroup *grp,
const HdfDataset &valuesDs,
const HdfDataset &activeDs,
hsize_t timeIndex );
~XmdfDataset() 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;
size_t activeData( size_t indexStart, size_t count, int *buffer ) override;
const HdfDataset &dsValues() const;
const HdfDataset &dsActive() const;
hsize_t timeIndex() const;
private:
HdfDataset mHdf5DatasetValues;
HdfDataset mHdf5DatasetActive;
// index or row where the data for this timestep begins
hsize_t mTimeIndex;
};
class LoaderXmdf
{
public:
@ -27,6 +61,7 @@ namespace MDAL
void load( Mesh *mesh, MDAL_Status *status );
private:
MDAL::Mesh *mMesh = nullptr;
std::string mDatFile;
std::shared_ptr<MDAL::DatasetGroup> readXmdfGroupAsDatasetGroup(
const HdfGroup &rootGroup,

421
external/mdal/mdal.cpp vendored
View File

@ -6,19 +6,21 @@
#include <string>
#include <stddef.h>
#include <limits>
#include <assert.h>
#include "mdal.h"
#include "mdal_loader.hpp"
#include "mdal_data_model.hpp"
#define NODATA std::numeric_limits<double>::quiet_NaN()
static const char *EMPTY_STR = "";
static MDAL_Status sLastStatus;
const char *MDAL_Version()
{
return "0.0.10";
return "0.1.0";
}
MDAL_Status MDAL_LastStatus()
@ -70,9 +72,29 @@ const char *MDAL_M_projection( MeshH mesh )
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
return _return_str( m->crs );
return _return_str( m->crs() );
}
void MDAL_M_extent( MeshH mesh, double *minX, double *maxX, double *minY, double *maxY )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
*minX = std::numeric_limits<double>::quiet_NaN();
*maxX = std::numeric_limits<double>::quiet_NaN();
*minY = std::numeric_limits<double>::quiet_NaN();
*maxY = std::numeric_limits<double>::quiet_NaN();
}
else
{
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
const MDAL::BBox extent = m->extent();
*minX = extent.minX;
*maxX = extent.maxX;
*minY = extent.minY;
*maxY = extent.maxY;
}
}
int MDAL_M_vertexCount( MeshH mesh )
{
@ -83,76 +105,10 @@ int MDAL_M_vertexCount( MeshH mesh )
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
int len = static_cast<int>( m->vertices.size() );
int len = static_cast<int>( m->verticesCount() );
return len;
}
double MDAL_M_vertexXCoordinatesAt( MeshH mesh, int index )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
if ( index < 0 )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
size_t i = static_cast<size_t>( index );
if ( m->vertices.size() <= i )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
return m->vertices[i].x;
}
double MDAL_M_vertexYCoordinatesAt( MeshH mesh, int index )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
if ( index < 0 )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
size_t i = static_cast<size_t>( index );
if ( m->vertices.size() <= i )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
return m->vertices[i].y;
}
double MDAL_M_vertexZCoordinatesAt( MeshH mesh, int index )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
if ( index < 0 )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
size_t i = static_cast<size_t>( index );
if ( m->vertices.size() <= i )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return NODATA;
}
return m->vertices[i].z;
}
int MDAL_M_faceCount( MeshH mesh )
{
if ( !mesh )
@ -161,11 +117,11 @@ int MDAL_M_faceCount( MeshH mesh )
return 0;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
int len = static_cast<int>( m->faces.size() );
int len = static_cast<int>( m->facesCount() );
return len;
}
int MDAL_M_faceVerticesCountAt( MeshH mesh, int index )
int MDAL_M_faceVerticesMaximumCount( MeshH mesh )
{
if ( !mesh )
{
@ -173,52 +129,7 @@ int MDAL_M_faceVerticesCountAt( MeshH mesh, int index )
return 0;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
if ( index < 0 )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
size_t i = static_cast<size_t>( index );
if ( m->faces.size() <= i )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
int len = static_cast<int>( m->faces[i].size() );
return len;
}
int MDAL_M_faceVerticesIndexAt( MeshH mesh, int face_index, int vertex_index )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
if ( face_index < 0 )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
size_t fi = static_cast<size_t>( face_index );
if ( m->faces.size() <= fi )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
if ( vertex_index < 0 )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
size_t vi = static_cast<size_t>( vertex_index );
if ( m->faces[fi].size() <= vi )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
int len = static_cast<int>( m->faces[fi][vi] );
int len = static_cast<int>( m->faceVerticesMaximumCount() );
return len;
}
@ -279,10 +190,109 @@ DatasetGroupH MDAL_M_datasetGroup( MeshH mesh, int index )
return static_cast< DatasetH >( m->datasetGroups[i].get() );
}
///////////////////////////////////////////////////////////////////////////////////////
/// MESH VERTICES
///////////////////////////////////////////////////////////////////////////////////////
MeshVertexIteratorH MDAL_M_vertexIterator( MeshH mesh )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return nullptr;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
std::unique_ptr<MDAL::MeshVertexIterator> it = m->readVertices();
return static_cast< MeshVertexIteratorH >( it.release() );
}
int MDAL_VI_next( MeshVertexIteratorH iterator, int verticesCount, double *coordinates )
{
if ( !iterator )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
MDAL::MeshVertexIterator *it = static_cast< MDAL::MeshVertexIterator * >( iterator );
size_t size = static_cast<size_t>( verticesCount );
if ( size == 0 )
{
return 0;
}
size_t ret = it->next( size, coordinates );
return static_cast<int>( ret );
}
void MDAL_VI_close( MeshVertexIteratorH iterator )
{
if ( iterator )
{
MDAL::MeshVertexIterator *it = static_cast< MDAL::MeshVertexIterator * >( iterator );
delete it;
}
}
///////////////////////////////////////////////////////////////////////////////////////
/// MESH FACES
///////////////////////////////////////////////////////////////////////////////////////
MeshFaceIteratorH MDAL_M_faceIterator( MeshH mesh )
{
if ( !mesh )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return nullptr;
}
MDAL::Mesh *m = static_cast< MDAL::Mesh * >( mesh );
std::unique_ptr<MDAL::MeshFaceIterator > it = m->readFaces();
return static_cast< MeshFaceIteratorH >( it.release() );
}
int MDAL_FI_next( MeshFaceIteratorH iterator,
int faceOffsetsBufferLen,
int *faceOffsetsBuffer,
int vertexIndicesBufferLen,
int *vertexIndicesBuffer )
{
if ( !iterator )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return 0;
}
MDAL::MeshFaceIterator *it = static_cast< MDAL::MeshFaceIterator * >( iterator );
size_t ret = it->next( static_cast<size_t>( faceOffsetsBufferLen ),
faceOffsetsBuffer,
static_cast<size_t>( vertexIndicesBufferLen ),
vertexIndicesBuffer );
return static_cast<int>( ret );
}
void MDAL_FI_close( MeshFaceIteratorH iterator )
{
if ( iterator )
{
MDAL::MeshFaceIterator *it = static_cast< MDAL::MeshFaceIterator * >( iterator );
delete it;
}
}
///////////////////////////////////////////////////////////////////////////////////////
/// DATASET GROUPS
///////////////////////////////////////////////////////////////////////////////////////
MeshH MDAL_G_mesh( DatasetGroupH group )
{
if ( !group )
{
sLastStatus = MDAL_Status::Err_IncompatibleDatasetGroup;
return nullptr;
}
MDAL::DatasetGroup *g = static_cast< MDAL::DatasetGroup * >( group );
MDAL::Mesh *m = g->mesh();
return static_cast< MeshH >( m );
}
int MDAL_G_datasetCount( DatasetGroupH group )
{
@ -388,7 +398,7 @@ bool MDAL_G_hasScalarData( DatasetGroupH group )
return true;
}
MDAL::DatasetGroup *g = static_cast< MDAL::DatasetGroup * >( group );
return g->isScalar;
return g->isScalar();
}
bool MDAL_G_isOnVertices( DatasetGroupH group )
@ -399,9 +409,32 @@ bool MDAL_G_isOnVertices( DatasetGroupH group )
return true;
}
MDAL::DatasetGroup *g = static_cast< MDAL::DatasetGroup * >( group );
return g->isOnVertices;
return g->isOnVertices();
}
void MDAL_G_minimumMaximum( DatasetGroupH group, double *min, double *max )
{
if ( !min || !max )
{
sLastStatus = MDAL_Status::Err_InvalidData;
return;
}
if ( !group )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
*min = NODATA;
*max = NODATA;
return;
}
MDAL::DatasetGroup *g = static_cast< MDAL::DatasetGroup * >( group );
MDAL::Statistics stats = g->statistics();
*min = stats.minimum;
*max = stats.maximum;
}
///////////////////////////////////////////////////////////////////////////////////////
/// DATASETS
///////////////////////////////////////////////////////////////////////////////////////
@ -414,7 +447,7 @@ DatasetGroupH MDAL_D_group( DatasetH dataset )
return nullptr;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
return static_cast< MDAL::DatasetGroup * >( d->parent );
return static_cast< MDAL::DatasetGroup * >( d->group() );
}
double MDAL_D_time( DatasetH dataset )
@ -425,7 +458,7 @@ double MDAL_D_time( DatasetH dataset )
return NODATA;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
return d->time;
return d->time();
}
@ -437,61 +470,10 @@ int MDAL_D_valueCount( DatasetH dataset )
return 0;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
int len = static_cast<int>( d->values.size() );
int len = static_cast<int>( d->valuesCount() );
return len;
}
double MDAL_D_value( DatasetH dataset, int valueIndex )
{
return MDAL_D_valueX( dataset, valueIndex );
}
double MDAL_D_valueX( DatasetH dataset, int valueIndex )
{
if ( !dataset )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return NODATA;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
int len = static_cast<int>( d->values.size() );
if ( len <= valueIndex )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return NODATA;
}
size_t i = static_cast<size_t>( valueIndex );
if ( d->values[i].noData )
{
return NODATA;
}
else
return d->values[i].x;
}
double MDAL_D_valueY( DatasetH dataset, int valueIndex )
{
if ( !dataset )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return NODATA;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
int len = static_cast<int>( d->values.size() );
if ( len <= valueIndex )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return NODATA;
}
size_t i = static_cast<size_t>( valueIndex );
if ( d->values[i].noData )
{
return NODATA;
}
else
return d->values[i].y;
}
bool MDAL_D_isValid( DatasetH dataset )
{
if ( !dataset )
@ -500,17 +482,100 @@ bool MDAL_D_isValid( DatasetH dataset )
return false;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
return d->isValid;
return d->isValid();
}
bool MDAL_D_active( DatasetH dataset, int faceIndex )
int MDAL_D_data( DatasetH dataset, int indexStart, int count, MDAL_DataType dataType, void *buffer )
{
if ( !dataset )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return false;
return 0;
}
MDAL::Dataset *d = static_cast< MDAL::Dataset * >( dataset );
size_t i = static_cast<size_t>( faceIndex );
return d->isActive( i );
size_t indexStartSizeT = static_cast<size_t>( indexStart );
size_t countSizeT = static_cast<size_t>( count );
MDAL::DatasetGroup *g = d->group();
assert( g );
MDAL::Mesh *m = d->mesh();
assert( m );
size_t valuesCount;
// Check that we are requesting correct 1D/2D for given dataset
switch ( dataType )
{
case MDAL_DataType::SCALAR_DOUBLE:
if ( !g->isScalar() )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = d->valuesCount();
break;
case MDAL_DataType::VECTOR_2D_DOUBLE:
if ( g->isScalar() )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
valuesCount = d->valuesCount();
break;
case MDAL_DataType::ACTIVE_INTEGER:
valuesCount = m->facesCount();
break;
}
// Check that we are not reaching out of values limit
if ( valuesCount <= indexStartSizeT )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
if ( valuesCount < indexStartSizeT + countSizeT )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return 0;
}
// Request data
size_t writtenValuesCount;
switch ( dataType )
{
case MDAL_DataType::SCALAR_DOUBLE:
writtenValuesCount = d->scalarData( indexStartSizeT, countSizeT, static_cast<double *>( buffer ) );
break;
case MDAL_DataType::VECTOR_2D_DOUBLE:
writtenValuesCount = d->vectorData( indexStartSizeT, countSizeT, static_cast<double *>( buffer ) );
break;
case MDAL_DataType::ACTIVE_INTEGER:
writtenValuesCount = d->activeData( indexStartSizeT, countSizeT, static_cast<int *>( buffer ) );
break;
}
return static_cast<int>( writtenValuesCount );
}
void MDAL_D_minimumMaximum( DatasetH dataset, double *min, double *max )
{
if ( !min || !max )
{
sLastStatus = MDAL_Status::Err_InvalidData;
return;
}
if ( !dataset )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
*min = NODATA;
*max = NODATA;
return;
}
MDAL::Dataset *ds = static_cast< MDAL::Dataset * >( dataset );
MDAL::Statistics stats = ds->statistics();
*min = stats.minimum;
*max = stats.maximum;
}

View File

@ -5,25 +5,87 @@
#include "mdal_data_model.hpp"
#include <assert.h>
#include <math.h>
#include <algorithm>
#include "mdal_utils.hpp"
bool MDAL::Dataset::isActive( size_t faceIndex )
MDAL::Dataset::~Dataset() = default;
MDAL::Dataset::Dataset( MDAL::DatasetGroup *parent )
: mParent( parent )
{
assert( parent );
if ( parent->isOnVertices )
assert( mParent );
}
size_t MDAL::Dataset::valuesCount() const
{
if ( group()->isOnVertices() )
{
if ( active.size() > faceIndex )
return active[faceIndex];
else
return false;
return mesh()->verticesCount();
}
else
{
return true;
return mesh()->facesCount();
}
}
MDAL::Statistics MDAL::Dataset::statistics() const
{
return mStatistics;
}
void MDAL::Dataset::setStatistics( const MDAL::Statistics &statistics )
{
mStatistics = statistics;
}
MDAL::DatasetGroup *MDAL::Dataset::group() const
{
return mParent;
}
MDAL::Mesh *MDAL::Dataset::mesh() const
{
return mParent->mesh();
}
double MDAL::Dataset::time() const
{
return mTime;
}
void MDAL::Dataset::setTime( double time )
{
mTime = time;
}
bool MDAL::Dataset::isValid() const
{
return mIsValid;
}
void MDAL::Dataset::setIsValid( bool isValid )
{
mIsValid = isValid;
}
MDAL::DatasetGroup::DatasetGroup( MDAL::Mesh *parent,
const std::string &uri,
const std::string &name )
: mParent( parent )
, mUri( uri )
{
assert( mParent );
setName( name );
}
MDAL::DatasetGroup::DatasetGroup( MDAL::Mesh *parent, const std::string &uri )
: mParent( parent )
, mUri( uri )
{
assert( mParent );
}
std::string MDAL::DatasetGroup::getMetadata( const std::string &key )
{
for ( auto &pair : metadata )
@ -61,9 +123,66 @@ void MDAL::DatasetGroup::setName( const std::string &name )
setMetadata( "name", name );
}
std::string MDAL::DatasetGroup::uri() const
{
return mUri;
}
MDAL::Statistics MDAL::DatasetGroup::statistics() const
{
return mStatistics;
}
void MDAL::DatasetGroup::setStatistics( const Statistics &statistics )
{
mStatistics = statistics;
}
MDAL::Mesh *MDAL::DatasetGroup::mesh() const
{
return mParent;
}
bool MDAL::DatasetGroup::isOnVertices() const
{
return mIsOnVertices;
}
void MDAL::DatasetGroup::setIsOnVertices( bool isOnVertices )
{
// 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;
}
bool MDAL::DatasetGroup::isScalar() const
{
return mIsScalar;
}
void MDAL::DatasetGroup::setIsScalar( bool isScalar )
{
// datasets are initialized (e.g. values array, active array) based
// on this property. Do not allow to modify later on.
assert( datasets.empty() );
mIsScalar = isScalar;
}
MDAL::Mesh::Mesh( size_t verticesCount, size_t facesCount, size_t faceVerticesMaximumCount, MDAL::BBox extent, const std::string &uri )
: mVerticesCount( verticesCount )
, mFacesCount( facesCount )
, mFaceVerticesMaximumCount( faceVerticesMaximumCount )
, mExtent( extent )
, mUri( uri )
{
}
MDAL::Mesh::~Mesh() = default;
void MDAL::Mesh::setSourceCrs( const std::string &str )
{
crs = MDAL::trim( str );
mCrs = MDAL::trim( str );
}
void MDAL::Mesh::setSourceCrsFromWKT( const std::string &wkt )
@ -76,26 +195,56 @@ void MDAL::Mesh::setSourceCrsFromEPSG( int code )
setSourceCrs( std::string( "EPSG:" ) + std::to_string( code ) );
}
void MDAL::Mesh::addBedElevationDataset()
void MDAL::Mesh::setExtent( const BBox &extent )
{
if ( faces.empty() )
return;
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >();
group->isOnVertices = true;
group->isScalar = true;
group->setName( "Bed Elevation" );
group->uri = uri;
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< Dataset >();
dataset->time = 0.0;
dataset->values.resize( vertices.size() );
dataset->active.resize( faces.size() );
dataset->parent = group.get();
std::fill( dataset->active.begin(), dataset->active.end(), 1 );
for ( size_t i = 0; i < vertices.size(); ++i )
{
dataset->values[i].x = vertices[i].z;
}
group->datasets.push_back( dataset );
datasetGroups.push_back( group );
mExtent = extent;
}
void MDAL::Mesh::setFaceVerticesMaximumCount( size_t faceVerticesMaximumCount )
{
mFaceVerticesMaximumCount = faceVerticesMaximumCount;
}
void MDAL::Mesh::setFacesCount( size_t facesCount )
{
mFacesCount = facesCount;
}
void MDAL::Mesh::setVerticesCount( size_t verticesCount )
{
mVerticesCount = verticesCount;
}
size_t MDAL::Mesh::verticesCount() const
{
return mVerticesCount;
}
size_t MDAL::Mesh::facesCount() const
{
return mFacesCount;
}
std::string MDAL::Mesh::uri() const
{
return mUri;
}
MDAL::BBox MDAL::Mesh::extent() const
{
return mExtent;
}
std::string MDAL::Mesh::crs() const
{
return mCrs;
}
size_t MDAL::Mesh::faceVerticesMaximumCount() const
{
return mFaceVerticesMaximumCount;
}
MDAL::MeshVertexIterator::~MeshVertexIterator() = default;
MDAL::MeshFaceIterator::~MeshFaceIterator() = default;

View File

@ -3,18 +3,20 @@
Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
*/
#ifndef MDAL_DEFINES_HPP
#define MDAL_DEFINES_HPP
#ifndef MDAL_DATA_MODEL_HPP
#define MDAL_DATA_MODEL_HPP
#include <stddef.h>
#include <vector>
#include <memory>
#include <map>
#include <string>
#include "mdal.h"
namespace MDAL
{
class DatasetGroup;
class Mesh;
struct BBox
{
@ -27,45 +29,42 @@ namespace MDAL
double maxY;
};
typedef struct
{
double x;
double y;
double z; // Bed elevation
} Vertex;
typedef std::vector<size_t> Face;
typedef std::vector<Vertex> Vertices;
typedef std::vector<Face> Faces;
typedef struct
{
double x;
double y;
bool noData = false;
} Value; //Dataset Value
double minimum = std::numeric_limits<double>::quiet_NaN();
double maximum = std::numeric_limits<double>::quiet_NaN();
} Statistics;
typedef std::vector< std::pair< std::string, std::string > > Metadata;
class Dataset
{
public:
double time;
Dataset( DatasetGroup *parent );
virtual ~Dataset();
/**
* size - face count if !isOnVertices
* size - vertex count if isOnVertices
*/
std::vector<Value> values;
std::vector<bool> active; // size - face count. Whether the output for this is active...
size_t valuesCount() const;
virtual size_t scalarData( size_t indexStart, size_t count, double *buffer ) = 0;
virtual size_t vectorData( size_t indexStart, size_t count, double *buffer ) = 0;
virtual size_t activeData( size_t indexStart, size_t count, int *buffer ) = 0;
bool isValid = true;
DatasetGroup *parent = nullptr;
Statistics statistics() const;
void setStatistics( const Statistics &statistics );
bool isActive( size_t faceIndex );
bool isValid() const;
void setIsValid( bool isValid );
DatasetGroup *group() const;
Mesh *mesh() const;
double time() const;
void setTime( double time );
private:
double mTime = std::numeric_limits<double>::quiet_NaN();
bool mIsValid = true;
DatasetGroup *mParent = nullptr;
Statistics mStatistics;
};
typedef std::vector<std::shared_ptr<Dataset>> Datasets;
@ -73,43 +72,103 @@ namespace MDAL
class DatasetGroup
{
public:
std::string getMetadata( const std::string &key );
DatasetGroup( Mesh *parent,
const std::string &uri );
DatasetGroup( Mesh *parent,
const std::string &uri,
const std::string &name );
std::string getMetadata( const std::string &key );
void setMetadata( const std::string &key, const std::string &val );
std::string name();
void setName( const std::string &name );
Metadata metadata;
bool isScalar = true;
bool isOnVertices = true;
Datasets datasets;
std::string uri; // file/uri from where it came
bool isScalar() const;
void setIsScalar( bool isScalar );
bool isOnVertices() const;
void setIsOnVertices( bool isOnVertices );
std::string uri() const;
Statistics statistics() const;
void setStatistics( const Statistics &statistics );
Mesh *mesh() const;
private:
Mesh *mParent = nullptr;
bool mIsScalar = true;
bool mIsOnVertices = true;
std::string mUri; // file/uri from where it came
Statistics mStatistics;
};
typedef std::vector<std::shared_ptr<DatasetGroup>> DatasetGroups;
struct Mesh
class MeshVertexIterator
{
std::string uri; // file/uri from where it came
std::string crs;
public:
virtual ~MeshVertexIterator();
Vertices vertices;
std::map<size_t, size_t> vertexIDtoIndex; // only for 2DM and DAT files
Faces faces;
std::map<size_t, size_t> faceIDtoIndex; // only for 2DM and DAT files
DatasetGroups datasetGroups;
void setSourceCrs( const std::string &str );
void setSourceCrsFromWKT( const std::string &wkt );
void setSourceCrsFromEPSG( int code );
void addBedElevationDataset();
virtual size_t next( size_t vertexCount, double *coordinates ) = 0;
};
} // namespace MDAL
#endif //MDAL_DEFINES_HPP
class MeshFaceIterator
{
public:
virtual ~MeshFaceIterator();
virtual size_t next( size_t faceOffsetsBufferLen,
int *faceOffsetsBuffer,
size_t vertexIndicesBufferLen,
int *vertexIndicesBuffer ) = 0;
};
class Mesh
{
public:
Mesh( size_t verticesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
BBox extent,
const std::string &uri );
virtual ~Mesh();
void setSourceCrs( const std::string &str );
void setSourceCrsFromWKT( const std::string &wkt );
void setSourceCrsFromEPSG( int code );
void setVerticesCount( size_t verticesCount );
void setFacesCount( size_t facesCount );
void setFaceVerticesMaximumCount( size_t faceVerticesMaximumCount );
void setExtent( const BBox &extent );
virtual std::unique_ptr<MDAL::MeshVertexIterator> readVertices() = 0;
virtual std::unique_ptr<MDAL::MeshFaceIterator> readFaces() = 0;
DatasetGroups datasetGroups;
size_t verticesCount() const;
size_t facesCount() const;
std::string uri() const;
BBox extent() const;
std::string crs() const;
size_t faceVerticesMaximumCount() const;
private:
size_t mVerticesCount = 0;
size_t mFacesCount = 0;
size_t mFaceVerticesMaximumCount = 0; //typically 3 or 4, sometimes up to 9
BBox mExtent;
const std::string mUri; // file/uri from where it came
std::string mCrs;
};
} // namespace MDAL
#endif //MDAL_DATA_MODEL_HPP

207
external/mdal/mdal_memory_data_model.cpp vendored Normal file
View File

@ -0,0 +1,207 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
*/
#include "mdal_memory_data_model.hpp"
#include <assert.h>
#include <math.h>
#include <cstring>
#include <algorithm>
#include <iterator>
#include "mdal_utils.hpp"
MDAL::MemoryDataset::MemoryDataset( MDAL::DatasetGroup *grp )
: Dataset( grp )
, mValues( group()->isScalar() ? valuesCount() : 2 * valuesCount(),
std::numeric_limits<double>::quiet_NaN() )
, mActive( group()->isOnVertices() ? mesh()->facesCount() : 0,
1 )
{
}
MDAL::MemoryDataset::~MemoryDataset() = default;
int *MDAL::MemoryDataset::active()
{
return mActive.data();
}
double *MDAL::MemoryDataset::values()
{
return mValues.data();
}
const int *MDAL::MemoryDataset::constActive() const
{
return mActive.data();
}
const double *MDAL::MemoryDataset::constValues() const
{
return mValues.data();
}
size_t MDAL::MemoryDataset::activeData( size_t indexStart, size_t count, int *buffer )
{
if ( group()->isOnVertices() )
{
size_t nValues = mActive.size();
if ( ( count < 1 ) || ( indexStart >= nValues ) )
return 0;
size_t copyValues = std::min( nValues - indexStart, count );
memcpy( buffer, constActive() + indexStart, copyValues * sizeof( int ) );
return copyValues;
}
else
{
memset( buffer, true, count * sizeof( int ) );
}
return count;
}
size_t MDAL::MemoryDataset::scalarData( size_t indexStart, size_t count, double *buffer )
{
assert( group()->isScalar() ); //checked in C API interface
size_t nValues = valuesCount();
assert( mValues.size() == nValues );
if ( ( count < 1 ) || ( indexStart >= nValues ) )
return 0;
size_t copyValues = std::min( nValues - indexStart, count );
memcpy( buffer, constValues() + indexStart, copyValues * sizeof( double ) );
return copyValues;
}
size_t MDAL::MemoryDataset::vectorData( size_t indexStart, size_t count, double *buffer )
{
assert( !group()->isScalar() ); //checked in C API interface
size_t nValues = valuesCount();
assert( mValues.size() == nValues * 2 );
if ( ( count < 1 ) || ( indexStart >= nValues ) )
return 0;
size_t copyValues = std::min( nValues - indexStart, count );
memcpy( buffer, constValues() + 2 * indexStart, 2 * copyValues * sizeof( double ) );
return copyValues;
}
MDAL::MemoryMesh::MemoryMesh( size_t verticesCount, size_t facesCount, size_t faceVerticesMaximumCount, MDAL::BBox extent, const std::string &uri )
: MDAL::Mesh( verticesCount, facesCount, faceVerticesMaximumCount, extent, uri )
{
}
std::unique_ptr<MDAL::MeshVertexIterator> MDAL::MemoryMesh::readVertices()
{
std::unique_ptr<MDAL::MemoryMeshVertexIterator> it( new MemoryMeshVertexIterator( this ) );
return it;
}
std::unique_ptr<MDAL::MeshFaceIterator> MDAL::MemoryMesh::readFaces()
{
std::unique_ptr<MDAL::MemoryMeshFaceIterator> it( new MemoryMeshFaceIterator( this ) );
return it;
}
void MDAL::MemoryMesh::addBedElevationDataset( const MDAL::Vertices &vertices, const MDAL::Faces &faces )
{
MDAL::addBedElevationDatasetGroup( this, vertices, faces );
}
MDAL::MemoryMesh::~MemoryMesh() = default;
MDAL::MemoryMeshVertexIterator::MemoryMeshVertexIterator( const MDAL::MemoryMesh *mesh )
: mMemoryMesh( mesh )
{
}
MDAL::MemoryMeshVertexIterator::~MemoryMeshVertexIterator() = default;
size_t MDAL::MemoryMeshVertexIterator::next( size_t vertexCount, double *coordinates )
{
assert( mMemoryMesh );
assert( coordinates );
size_t maxVertices = mMemoryMesh->verticesCount();
if ( vertexCount > maxVertices )
return 0;
if ( mLastVertexIndex >= maxVertices )
return 0;
size_t i = 0;
while ( true )
{
if ( mLastVertexIndex + i >= maxVertices )
break;
if ( i >= vertexCount )
break;
const Vertex v = mMemoryMesh->vertices[mLastVertexIndex + i];
coordinates[3 * i] = v.x;
coordinates[3 * i + 1] = v.y;
coordinates[3 * i + 2] = v.z;
++i;
}
mLastVertexIndex += i;
return i;
}
MDAL::MemoryMeshFaceIterator::MemoryMeshFaceIterator( const MDAL::MemoryMesh *mesh )
: mMemoryMesh( mesh )
{
}
MDAL::MemoryMeshFaceIterator::~MemoryMeshFaceIterator() = default;
size_t MDAL::MemoryMeshFaceIterator::next(
size_t faceOffsetsBufferLen, int *faceOffsetsBuffer,
size_t vertexIndicesBufferLen, int *vertexIndicesBuffer )
{
assert( mMemoryMesh );
assert( faceOffsetsBuffer );
assert( vertexIndicesBuffer );
size_t maxFaces = mMemoryMesh->facesCount();
size_t faceVerticesMaximumCount = mMemoryMesh->faceVerticesMaximumCount();
size_t vertexIndex = 0;
size_t faceIndex = 0;
while ( true )
{
if ( vertexIndex + faceVerticesMaximumCount > vertexIndicesBufferLen )
break;
if ( faceIndex >= faceOffsetsBufferLen )
break;
if ( mLastFaceIndex + faceIndex >= maxFaces )
break;
const Face f = mMemoryMesh->faces[mLastFaceIndex + faceIndex];
for ( size_t faceVertexIndex = 0; faceVertexIndex < f.size(); ++faceVertexIndex )
{
assert( vertexIndex < vertexIndicesBufferLen );
vertexIndicesBuffer[vertexIndex] = static_cast<int>( f[faceVertexIndex] );
++vertexIndex;
}
assert( faceIndex < faceOffsetsBufferLen );
faceOffsetsBuffer[faceIndex] = static_cast<int>( vertexIndex );
++faceIndex;
}
mLastFaceIndex += faceIndex;
return faceIndex;
}

125
external/mdal/mdal_memory_data_model.hpp vendored Normal file
View File

@ -0,0 +1,125 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
*/
#ifndef MDAL_MEMORY_DATA_MODEL_HPP
#define MDAL_MEMORY_DATA_MODEL_HPP
#include <stddef.h>
#include <vector>
#include <memory>
#include <map>
#include <string>
#include "mdal.h"
#include "mdal_data_model.hpp"
namespace MDAL
{
typedef struct
{
double x;
double y;
double z; // Bed elevation
} Vertex;
typedef std::vector<size_t> Face;
typedef std::vector<Vertex> Vertices;
typedef std::vector<Face> Faces;
/**
* The MemoryDataset stores all the data in the memory
*/
class MemoryDataset: public Dataset
{
public:
MemoryDataset( DatasetGroup *grp );
~MemoryDataset() 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;
size_t activeData( size_t indexStart, size_t count, int *buffer ) override;
/**
* valid pointer only for for dataset defined on vertices
*/
int *active();
double *values();
const int *constActive() const;
const double *constValues() const;
private:
/**
* Stores vector2d/scalar data for dataset in form
* scalars: x1, x2, x3, ..., xN
* vector2D: x1, y1, x2, y2, x3, y3, .... , xN, yN
*
* all values are initialized to std::numerical_limits<double>::quiet_NaN (==NODATA)
*
* size:
* - face count if isOnFaces & isScalar
* - vertex count if isOnVertices & isScalar
* - face count * 2 if isOnFaces & isVector
* - 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
* For dataset defined on faces, this is empty vector
*
* Values are initialized by default to 1 (active)
*/
std::vector<int> mActive;
};
class MemoryMesh: public Mesh
{
public:
MemoryMesh( size_t verticesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
BBox extent,
const std::string &uri );
~MemoryMesh() override;
std::unique_ptr<MDAL::MeshVertexIterator> readVertices() override;
std::unique_ptr<MDAL::MeshFaceIterator> readFaces() override;
void addBedElevationDataset( const Vertices &vertices, const Faces &faces );
Vertices vertices;
Faces faces;
};
class MemoryMeshVertexIterator: public MeshVertexIterator
{
public:
MemoryMeshVertexIterator( const MemoryMesh *mesh );
~MemoryMeshVertexIterator() override;
size_t next( size_t vertexCount, double *coordinates ) override;
const MemoryMesh *mMemoryMesh;
size_t mLastVertexIndex = 0;
};
class MemoryMeshFaceIterator: public MeshFaceIterator
{
public:
MemoryMeshFaceIterator( const MemoryMesh *mesh );
~MemoryMeshFaceIterator() override;
size_t next( size_t faceOffsetsBufferLen,
int *faceOffsetsBuffer,
size_t vertexIndicesBufferLen,
int *vertexIndicesBuffer ) override;
const MemoryMesh *mMemoryMesh;
size_t mLastFaceIndex = 0;
};
} // namespace MDAL
#endif //MDAL_MEMORY_DATA_MODEL_HPP

View File

@ -10,6 +10,7 @@
#include <algorithm>
#include <sstream>
#include <math.h>
#include <assert.h>
bool MDAL::fileExists( const std::string &filename )
{
@ -254,3 +255,142 @@ double MDAL::parseTimeUnits( const std::string &units )
return divBy;
}
MDAL::Statistics _calculateStatistics( const std::vector<double> &values, size_t count, bool isVector )
{
MDAL::Statistics ret;
double min = std::numeric_limits<double>::quiet_NaN();
double max = std::numeric_limits<double>::quiet_NaN();
bool firstIteration = true;
for ( size_t i = 0; i < count; ++i )
{
double magnitude;
if ( isVector )
{
double x = values[2 * i];
double y = values[2 * i + 1];
if ( isnan( x ) || isnan( y ) )
continue;
magnitude = sqrt( x * x + y * y );
}
else
{
double x = values[i];
if ( isnan( x ) )
continue;
magnitude = x;
}
if ( firstIteration )
{
firstIteration = false;
min = magnitude;
max = magnitude;
}
else
{
if ( magnitude < min )
{
min = magnitude;
}
if ( magnitude > max )
{
max = magnitude;
}
}
}
ret.minimum = min;
ret.maximum = max;
return ret;
}
MDAL::Statistics MDAL::calculateStatistics( std::shared_ptr<DatasetGroup> grp )
{
Statistics ret;
if ( !grp )
return ret;
for ( std::shared_ptr<Dataset> ds : grp->datasets )
{
MDAL::Statistics dsStats = ds->statistics();
combineStatistics( ret, dsStats );
}
return ret;
}
MDAL::Statistics MDAL::calculateStatistics( std::shared_ptr<Dataset> dataset )
{
Statistics ret;
if ( !dataset )
return ret;
bool isVector = !dataset->group()->isScalar();
size_t bufLen = 2000;
std::vector<double> buffer( isVector ? bufLen * 2 : bufLen );
size_t i = 0;
while ( i < dataset->valuesCount() )
{
size_t valsRead;
if ( isVector )
{
valsRead = dataset->vectorData( i, bufLen, buffer.data() );
}
else
{
valsRead = dataset->scalarData( i, bufLen, buffer.data() );
}
MDAL::Statistics dsStats = _calculateStatistics( buffer, valsRead, isVector );
combineStatistics( ret, dsStats );
i += valsRead;
}
return ret;
}
void MDAL::combineStatistics( MDAL::Statistics &main, const MDAL::Statistics &other )
{
if ( isnan( main.minimum ) ||
( !isnan( other.minimum ) && ( main.minimum > other.minimum ) ) )
{
main.minimum = other.minimum;
}
if ( isnan( main.maximum ) ||
( !isnan( other.maximum ) && ( main.maximum < other.maximum ) ) )
{
main.maximum = other.maximum;
}
}
void MDAL::addBedElevationDatasetGroup( MDAL::Mesh *mesh, const Vertices &vertices, const Faces &faces )
{
if ( !mesh )
return;
if ( 0 == mesh->facesCount() )
return;
std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
mesh,
mesh->uri(),
"Bed Elevation"
);
group->setIsOnVertices( true );
group->setIsScalar( true );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MemoryDataset >( group.get() );
dataset->setTime( 0.0 );
double *vals = dataset->values();
for ( size_t i = 0; i < vertices.size(); ++i )
{
vals[i] = vertices[i].z;
}
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
}

View File

@ -12,9 +12,11 @@
#include <limits>
#include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp"
// avoid unused variable warnings
#define MDAL_UNUSED(x) (void)x;
#define MDAL_NAN std::numeric_limits<double>::quiet_NaN()
namespace MDAL
{
@ -88,9 +90,21 @@ namespace MDAL
BBox computeExtent( const Vertices &vertices );
// time
//! Returns a delimiter to get time in hours
double parseTimeUnits( const std::string &units );
// statistics
void combineStatistics( Statistics &main, const Statistics &other );
//! Calculates statistics for dataset group
Statistics calculateStatistics( std::shared_ptr<DatasetGroup> grp );
//! Calculates statistics for dataset
Statistics calculateStatistics( std::shared_ptr<Dataset> dataset );
// mesh & datasets
//! Add bed elevatiom dataset group to mesh
void addBedElevationDatasetGroup( MDAL::Mesh *mesh, const Vertices &vertices, const Faces &faces );
} // namespace MDAL
#endif //MDAL_UTILS_HPP

View File

@ -36,6 +36,7 @@ IF (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/mdal_utils.cpp
${CMAKE_SOURCE_DIR}/external/mdal/mdal_loader.cpp
${CMAKE_SOURCE_DIR}/external/mdal/mdal_data_model.cpp
${CMAKE_SOURCE_DIR}/external/mdal/mdal_memory_data_model.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_2dm.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ascii_dat.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_binary_dat.cpp
@ -46,6 +47,7 @@ IF (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/mdal_utils.hpp
${CMAKE_SOURCE_DIR}/external/mdal/mdal_loader.hpp
${CMAKE_SOURCE_DIR}/external/mdal/mdal_data_model.hpp
${CMAKE_SOURCE_DIR}/external/mdal/mdal_memory_data_model.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_2dm.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ascii_dat.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_binary_dat.hpp
@ -53,6 +55,7 @@ IF (WITH_INTERNAL_MDAL)
IF(HDF5_FOUND)
SET(MDAL_LIB_SRCS ${MDAL_LIB_SRCS}
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_hdf5.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_xmdf.cpp
)
SET(MDAL_LIB_HDRS ${MDAL_LIB_HDRS}