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

update MDAL 0.4.0

Browse Source
This commit is contained in:
Peter Petrik 2019-10-14 09:19:14 +02:00
parent 94fa3c44e0
commit 679e75e13c
33 changed files with 1148 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -56,6 +56,7 @@ enum MDAL_Status
Err_IncompatibleDatasetGroup,
Err_MissingDriver,
Err_MissingDriverCapability,
Err_FailToWriteToDisk,
// Warnings
Warn_UnsupportedElement,
Warn_InvalidElements,
@ -105,6 +106,9 @@ MDAL_EXPORT bool MDAL_DR_meshLoadCapability( DriverH driver );
//! Returns whether driver has capability to write/edit dataset (groups)
MDAL_EXPORT bool MDAL_DR_writeDatasetsCapability( DriverH driver );
//! Returns whether driver has capability to save mesh
MDAL_EXPORT bool MDAL_DR_SaveMeshCapability( DriverH driver );
/**
* Returns name of MDAL driver
* not thread-safe and valid only till next call
@ -134,9 +138,13 @@ MDAL_EXPORT const char *MDAL_DR_filters( DriverH driver );
* Caller must free memory with MDAL_CloseMesh() afterwards
*/
MDAL_EXPORT MeshH MDAL_LoadMesh( const char *meshFile );
//! Closes mesh, frees the memory
MDAL_EXPORT void MDAL_CloseMesh( MeshH mesh );
//! Saves mesh (only mesh structure) on a file with the specified driver. On error see MDAL_LastStatus for error type.
MDAL_EXPORT void MDAL_SaveMesh( MeshH mesh, const char *meshFile, const char *driver );
/**
* Returns mesh projection
* not thread-safe and valid only till next call
@ -348,6 +356,11 @@ MDAL_EXPORT bool MDAL_G_isInEditMode( DatasetGroupH group );
*/
MDAL_EXPORT void MDAL_G_closeEditMode( DatasetGroupH group );
/**
* Returns reference time for dataset group
* If returned value begins with word JULIAN, following number represents date in Julian format
*/
MDAL_EXPORT const char *MDAL_G_referenceTime( DatasetGroupH group );
///////////////////////////////////////////////////////////////////////////////////////
/// DATASETS

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

@ -82,7 +82,7 @@ MDAL::Driver2dm::Driver2dm():
Driver( DRIVER_NAME,
"2DM Mesh File",
"*.2dm",
Capability::ReadMesh
Capability::ReadMesh | Capability::SaveMesh
)
{
}
@ -271,7 +271,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
new Mesh2dm(
vertices.size(),
faces.size(),
4, //maximum quads
MAX_VERTICES_PER_FACE_2DM,
computeExtent( vertices ),
mMeshFile,
vertexIDtoIndex
@ -286,3 +286,66 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
return std::unique_ptr<Mesh>( mesh.release() );
}
void MDAL::Driver2dm::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Status *status )
{
if ( status ) *status = MDAL_Status::None;
std::ofstream file( uri, std::ofstream::out );
if ( !file.is_open() )
{
if ( status ) *status = MDAL_Status::Err_FailToWriteToDisk;
}
std::string line = "MESH2D";
file << line << std::endl;
//write vertices
std::unique_ptr<MDAL::MeshVertexIterator> vertexIterator = mesh->readVertices();
double vertex[3];
for ( size_t i = 0; i < mesh->verticesCount(); ++i )
{
vertexIterator->next( 1, vertex );
line = "ND ";
line.append( std::to_string( i + 1 ) );
for ( size_t j = 0; j < 2; ++j )
{
line.append( " " );
line.append( MDAL::coordinateToString( vertex[j] ) );
}
line.append( " " );
line.append( MDAL::doubleToString( vertex[2] ) );
file << line << std::endl;
}
//write faces
std::unique_ptr<MDAL::MeshFaceIterator> faceIterator = mesh->readFaces();
for ( size_t i = 0; i < mesh->facesCount(); ++i )
{
int faceOffsets[1];
int vertexIndices[4]; //max 4 vertices for a face
faceIterator->next( 1, faceOffsets, 4, vertexIndices );
if ( faceOffsets[0] > 2 && faceOffsets[0] < 5 )
{
if ( faceOffsets[0] == 3 )
line = "E3T ";
if ( faceOffsets[0] == 4 )
line = "E4Q ";
line.append( std::to_string( i + 1 ) );
for ( int j = 0; j < faceOffsets[0]; ++j )
{
line.append( " " );
line.append( std::to_string( vertexIndices[j] + 1 ) );
}
}
file << line << std::endl;
}
file.close();
}

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

@ -14,6 +14,9 @@
#include "mdal.h"
#include "mdal_driver.hpp"
#define MAX_VERTICES_PER_FACE_2DM 4
namespace MDAL
{
class Mesh2dm: public MemoryMesh
@ -80,8 +83,12 @@ namespace MDAL
~Driver2dm() override;
Driver2dm *create() override;
int faceVerticesMaximumCount() const override
{return MAX_VERTICES_PER_FACE_2DM;}
bool canRead( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile, MDAL_Status *status ) override;
void save( const std::string &uri, Mesh *mesh, MDAL_Status *status ) override;
private:
std::string mMeshFile;

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

@ -142,7 +142,7 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
std::shared_ptr<DatasetGroup> group; // DAT outputs data
std::string groupName( MDAL::baseName( mDatFile ) );
std::string line;
std::string referenceTime;
// see if it contains face-centered results - supported by BASEMENT
bool faceCentered = false;
if ( contains( groupName, "_els_" ) )
@ -201,6 +201,7 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
);
group->setIsScalar( !isVector );
group->setIsOnVertices( !faceCentered );
group->setReferenceTime( referenceTime );
}
else if ( cardType == "ENDDS" )
{
@ -226,9 +227,24 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
if ( quoteIdx1 != std::string::npos && quoteIdx2 != std::string::npos )
group->setName( line.substr( quoteIdx1 + 1, quoteIdx2 - quoteIdx1 - 1 ) );
}
else if ( cardType == "RT_JULIAN" && items.size() >= 2 )
{
referenceTime = "JULIAN " + items[1];
}
else if ( cardType == "TIMEUNITS" && items.size() >= 2 )
{
if ( !group )
{
debug( "TIMEUNITS card for no active dataset!" );
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
}
group->setMetadata( "TIMEUNITS", items[1] );
}
else if ( cardType == "TS" && items.size() >= 3 )
{
double t = toDouble( items[2] );
t = convertTimeDataToHours( t, group->getMetadata( "TIMEUNITS" ) );
if ( faceCentered )
{
@ -320,7 +336,7 @@ void MDAL::DriverAsciiDat::readVertexTimestep(
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
dataset->setTime( t / 3600. ); // TODO read TIMEUNITS
dataset->setTime( t );
int *active = dataset->active();
// only for new format
@ -392,7 +408,7 @@ void MDAL::DriverAsciiDat::readFaceTimestep(
size_t faceCount = mesh->facesCount();
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared< MDAL::MemoryDataset >( group.get() );
dataset->setTime( t / 3600. );
dataset->setTime( t );
double *values = dataset->values();
// TODO: hasStatus
for ( size_t index = 0; index < faceCount; ++index )
@ -427,3 +443,22 @@ void MDAL::DriverAsciiDat::readFaceTimestep(
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}
double MDAL::DriverAsciiDat::convertTimeDataToHours( double time, const std::string &originalTimeDataUnit ) const
{
if ( originalTimeDataUnit == "se" || originalTimeDataUnit == "2" || originalTimeDataUnit == "Seconds"
|| originalTimeDataUnit.empty() )
{
time /= 3600.0;
}
else if ( originalTimeDataUnit == "mi" || originalTimeDataUnit == "1" || originalTimeDataUnit == "Minutes" )
{
time /= 60.0;
}
else if ( originalTimeDataUnit == "days" )
{
time *= 24;
}
return time;
}

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

@ -84,6 +84,8 @@ namespace MDAL
bool isVector,
std::ifstream &stream ) const;
double convertTimeDataToHours( double time, const std::string &originalTimeDataUnit ) const;
std::string mDatFile;
};

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

@ -33,8 +33,8 @@ static const int CT_NUMCELLS = 180;
static const int CT_NAME = 190;
static const int CT_TS = 200;
static const int CT_ENDDS = 210;
//static const int CT_RT_JULIAN = 240;
//static const int CT_TIMEUNITS = 250;
static const int CT_RT_JULIAN = 240;
static const int CT_TIMEUNITS = 250;
static const int CT_2D_MESHES = 3;
static const int CT_FLOAT_SIZE = 4;
static const int CF_FLAG_SIZE = 1;
@ -142,6 +142,9 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
int numdata;
int numcells;
char groupName[40];
double referenceTime;
int timeUnit = 0;
std::string timeUnitStr;
char istat;
float time;
@ -245,6 +248,44 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
groupMax->setName( group->name() + "/Maximums" );
break;
case CT_RT_JULIAN:
// Reference time
if ( readIStat( in, sflg, &istat ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
if ( read( in, reinterpret_cast< char * >( &time ), 8 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
referenceTime = static_cast<double>( time );
group->setReferenceTime( "JULIAN " + std::to_string( referenceTime ) );
break;
case CT_TIMEUNITS:
// Time unit
if ( read( in, reinterpret_cast< char * >( &timeUnit ), 4 ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
switch ( timeUnit )
{
case 0:
timeUnitStr = "hours";
break;
case 1:
timeUnitStr = "minutes";
break;
case 2:
timeUnitStr = "seconds";
break;
case 4:
timeUnitStr = "days";
break;
default:
timeUnitStr = "unknown";
break;
}
group->setMetadata( "TIMEUNITS", timeUnitStr );
break;
case CT_TS:
// Time step!
if ( readIStat( in, sflg, &istat ) )
@ -254,6 +295,8 @@ void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh,
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
double t = static_cast<double>( time );
t = convertTimeDataToHours( t, timeUnit );
if ( readVertexTimestep( mesh, group, groupMax, t, istat, sflg, in ) )
EXIT_WITH_ERROR( MDAL_Status::Err_UnknownFormat );
@ -461,3 +504,23 @@ bool MDAL::DriverBinaryDat::persist( MDAL::DatasetGroup *group )
return false;
}
double MDAL::DriverBinaryDat::convertTimeDataToHours( double time, int originalTimeDataUnit )
{
switch ( originalTimeDataUnit )
{
case 1:
time /= 60.0;
break;
case 2:
time /= 3600.0;
break;
case 4:
time *= 24;
break;
case 0:
default:
break;
}
return time;
}

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

@ -40,6 +40,7 @@ namespace MDAL
int sflg,
std::ifstream &in );
double convertTimeDataToHours( double time, int originalTimeDataUnit );
std::string mDatFile;
};

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

@ -71,8 +71,24 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
}
else
{
/*UGRID convention says "The order of dimensions on a data variable is arbitrary",
*So time dimension is not necessary the first dimension, event if the convention recommands it.
*And prevent that any of the the two dimensions is not time dimension
*/
if ( mDimensions.type( dimids[0] ) == CFDimensions::Time )
{
dimid = dimids[1];
}
else if ( mDimensions.type( dimids[1] ) == CFDimensions::Time )
{
dimid = dimids[0];
}
else
{
continue;
}
nTimesteps = mDimensions.size( CFDimensions::Time );
dimid = dimids[1];
}
if ( !mDimensions.isDatasetType( mDimensions.type( dimid ) ) )
@ -169,6 +185,34 @@ std::shared_ptr<MDAL::Dataset> MDAL::DriverCF::createFace2DDataset( std::shared_
return dataset;
}
std::shared_ptr<MDAL::Dataset> MDAL::DriverCF::createVertex2DDataset( std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts, const MDAL::CFDatasetGroupInfo &dsi,
const std::vector<double> &vals_x,
const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y )
{
assert( dsi.outputType == CFDimensions::Vertex2D );
size_t nVertices2D = mDimensions.size( CFDimensions::Vertex2D );
std::shared_ptr<MDAL::MemoryDataset> dataset = std::make_shared<MDAL::MemoryDataset>( group.get() );
for ( size_t i = 0; i < nVertices2D; ++i )
{
size_t idx = ts * nVertices2D + i;
populate_vals( dsi.is_vector,
dataset->values(),
i,
vals_x,
vals_y,
idx,
fill_val_x,
fill_val_y );
}
return dataset;
}
void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<double> &times, const MDAL::cfdataset_info_map &dsinfo_map )
{
/* PHASE 2 - add dataset groups */
@ -183,7 +227,16 @@ void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
dsi.name
);
group->setIsScalar( !dsi.is_vector );
group->setIsOnVertices( false );
if ( dsi.outputType == CFDimensions::Vertex2D )
group->setIsOnVertices( true );
else if ( dsi.outputType == CFDimensions::Face2D )
group->setIsOnVertices( false );
else
{
// unsupported
continue;
}
// read X data
double fill_val_x = mNcFile.getFillValue( dsi.ncid_x );
@ -209,6 +262,13 @@ void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
if ( dsi.outputType == CFDimensions::Face2D )
{
dataset = createFace2DDataset( group, ts, dsi, vals_x, vals_y, fill_val_x, fill_val_y );
}
else // Vertex2D
{
dataset = createVertex2DDataset( group, ts, dsi, vals_x, vals_y, fill_val_x, fill_val_y );
}
if ( dataset )
{
dataset->setTime( time );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
@ -228,6 +288,13 @@ void MDAL::DriverCF::parseTime( std::vector<double> &times )
{
size_t nTimesteps = mDimensions.size( CFDimensions::Time );
if ( 0 == nTimesteps )
{
//if no time dimension is present creates only one time step to store the potential time-independent variable
nTimesteps = 1;
times = std::vector<double>( 1, 0 );
return;
}
times = mNcFile.readDoubleArr( "time", nTimesteps );
std::string units = mNcFile.getAttrStr( "time", "units" );
double div_by = MDAL::parseTimeUnits( units );

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

@ -94,6 +94,14 @@ namespace MDAL
const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y );
std::shared_ptr<MDAL::Dataset> createVertex2DDataset(
std::shared_ptr<MDAL::DatasetGroup> group,
size_t ts,
const MDAL::CFDatasetGroupInfo &dsi,
const std::vector<double> &vals_x,
const std::vector<double> &vals_y,
double fill_val_x, double fill_val_y );
void addDatasetGroups( Mesh *mesh,
const std::vector<double> &times,
const cfdataset_info_map &dsinfo_map );

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

@ -42,6 +42,11 @@ bool MDAL::Driver::hasCapability( MDAL::Capability capability ) const
return capability == ( mCapabilityFlags & capability );
}
int MDAL::Driver::faceVerticesMaximumCount() const
{
return -1;
}
std::unique_ptr< MDAL::Mesh > MDAL::Driver::load( const std::string &uri, MDAL_Status *status )
{
MDAL_UNUSED( uri );
@ -57,6 +62,13 @@ void MDAL::Driver::load( const std::string &uri, Mesh *mesh, MDAL_Status *status
return;
}
void MDAL::Driver::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Status *status )
{
MDAL_UNUSED( uri );
MDAL_UNUSED( mesh );
MDAL_UNUSED( status );
}
void MDAL::Driver::createDatasetGroup( MDAL::Mesh *mesh, const std::string &groupName, bool isOnVertices, bool hasScalarData, const std::string &datasetGroupFile )
{
std::shared_ptr<MDAL::DatasetGroup> grp(

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

@ -18,6 +18,7 @@ namespace MDAL
ReadMesh = 1 << 0, //! Can read mesh and all datasets stored in the mesh file
ReadDatasets = 1 << 1, //! Can read only datasets (groups) from existing mesh
WriteDatasets = 1 << 2, //! Can write datasets (groups)
SaveMesh = 1 << 3, //! Can save the mesh
};
class Driver
@ -39,11 +40,17 @@ namespace MDAL
virtual bool canRead( const std::string &uri ) = 0;
//! returns the maximum vertices per face
virtual int faceVerticesMaximumCount() const;
// loads mesh
virtual std::unique_ptr< Mesh > load( const std::string &uri, MDAL_Status *status );
// loads datasets
virtual void load( const std::string &uri, Mesh *mesh, MDAL_Status *status );
// save mesh
virtual void save( const std::string &uri, Mesh *mesh, MDAL_Status *status );
// create new dataset group
virtual void createDatasetGroup(
Mesh *mesh,

300
external/mdal/frmts/mdal_esri_tin.cpp vendored Normal file
View File

@ -0,0 +1,300 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Vincent Cloarec (vcloarec at gmail dot com)
*/
#include "mdal_esri_tin.hpp"
MDAL::DriverEsriTin::DriverEsriTin(): Driver( "ESRI_TIN",
"Esri TIN",
"*.adf",
Capability::ReadMesh )
{}
MDAL::Driver *MDAL::DriverEsriTin::create()
{
return new DriverEsriTin();
}
std::unique_ptr<MDAL::Mesh> MDAL::DriverEsriTin::load( const std::string &uri, MDAL_Status *status )
{
if ( status ) *status = MDAL_Status::None;
try
{
std::list<int> superpointIndexes;
Vertices vertices;
Faces faces;
readSuperpoints( uri, superpointIndexes );
populateVertices( uri, vertices, superpointIndexes );
populateFaces( uri, faces, superpointIndexes );
std::unique_ptr< MemoryMesh > mesh(
new MemoryMesh(
name(),
vertices.size(),
faces.size(),
3,
computeExtent( vertices ),
uri
)
);
mesh->faces = std::move( faces );
mesh->vertices = std::move( vertices );
addBedElevationDatasetGroup( mesh.get(), mesh->vertices );
mesh->datasetGroups.back()->setName( "Altitude" );
std::string crs = getCrsWkt( uri );
if ( ! crs.empty() )
mesh->setSourceCrsFromWKT( crs );
return std::unique_ptr<Mesh>( mesh.release() );
}
catch ( MDAL_Status error )
{
if ( status ) *status = ( error );
return std::unique_ptr<Mesh>();
}
}
bool MDAL::DriverEsriTin::canRead( const std::string &uri )
{
std::string zFileName = zFile( uri );
std::string faceFileName = faceFile( uri );
std::ifstream xyIn( xyFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! xyIn.is_open() )
return false;
std::ifstream zIn( zFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! zIn.is_open() )
return false;
std::ifstream faceIn( faceFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! faceIn.is_open() )
return false;
std::ifstream hullIn( hullFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! hullIn.is_open() )
return false;
return true;
}
std::string MDAL::DriverEsriTin::xyFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tnxy.adf" );
}
std::string MDAL::DriverEsriTin::zFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tnz.adf" );
}
std::string MDAL::DriverEsriTin::faceFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tnod.adf" );
}
std::string MDAL::DriverEsriTin::mskFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tmsk.adf" );
}
std::string MDAL::DriverEsriTin::msxFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "tmsx.adf" );
}
std::string MDAL::DriverEsriTin::hullFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "thul.adf" );
}
std::string MDAL::DriverEsriTin::crsFile( const std::string &uri ) const
{
return pathJoin( dirName( uri ), "prj.adf" );
}
void MDAL::DriverEsriTin::readSuperpoints( const std::string &uri, std::list<int> &superpointsIndexes ) const
{
superpointsIndexes.clear();
bool isNativeLittleEndian = MDAL::isNativeLittleEndian();
std::ifstream inHull( hullFile( uri ), std::ifstream::in | std::ifstream::binary );
int32_t index;
while ( readValue( index, inHull, isNativeLittleEndian ) && index != -1 )
superpointsIndexes.push_back( index );
superpointsIndexes.sort();
}
std::string MDAL::DriverEsriTin::getTinName( const std::string &uri ) const
{
std::string tinName = uri;
size_t last_slash_idx = tinName.find_last_of( "\\/" );
if ( last_slash_idx == std::string::npos )
return "";
tinName.erase( last_slash_idx, tinName.size() - last_slash_idx );
last_slash_idx = tinName.find_last_of( "\\/" );
if ( last_slash_idx == std::string::npos )
return "";
tinName.erase( 0, last_slash_idx + 1 );
return tinName;
}
void MDAL::DriverEsriTin::populateVertices( const std::string &uri, MDAL::Vertices &vertices, const std::list<int> &superpointIndexes ) const
{
//first, read the (x,y) coordinates
bool isNativeLittleEndian = MDAL::isNativeLittleEndian();
std::ifstream inXY( xyFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inZ( zFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! inXY.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inZ.is_open() )
throw MDAL_Status::Err_FileNotFound;
int fileIndex = 1;
while ( true )
{
Vertex vert;
if ( !readValue( vert.x, inXY, isNativeLittleEndian ) )
break;
if ( !readValue( vert.y, inXY, isNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
float zValue;
if ( !readValue( zValue, inZ, isNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
vert.z = double( zValue );
if ( correctedIndex( fileIndex, superpointIndexes ) >= 0 ) //store the vertex only if it is not a superpoint
vertices.push_back( vert );
fileIndex++;
}
inXY.close();
inZ.close();
}
void MDAL::DriverEsriTin::populateFaces( const std::string &uri, MDAL::Faces &faces, const std::list<int> &superpointIndexes ) const
{
bool isNativeLittleEndian = MDAL::isNativeLittleEndian();
std::ifstream inFaces( faceFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inMsk( mskFile( uri ), std::ifstream::in | std::ifstream::binary );
std::ifstream inMsx( msxFile( uri ), std::ifstream::in | std::ifstream::binary );
if ( ! inFaces.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inMsk.is_open() )
throw MDAL_Status::Err_FileNotFound;
if ( ! inMsx.is_open() )
throw MDAL_Status::Err_FileNotFound;
//Find the beginning of data in mskFile
inMsx.seekg( -4, std::ios::end );
int32_t mskBegin;
if ( ! readValue( mskBegin, inMsx, true ) )
throw MDAL_Status::Err_UnknownFormat;
//read information in mskFile
inMsk.seekg( -mskBegin * 2, std::ios::end );
int32_t maskIntergerCount;
if ( ! readValue( maskIntergerCount, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
inMsk.ignore( 4 ); //unused 4 bytes
int32_t maskBitsCount;
if ( ! readValue( maskBitsCount, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
int c = 0;
int32_t maskInt = 0;
while ( true )
{
//read mask file
if ( c % 32 == 0 && c < maskBitsCount ) //first bit in the mask array have to be used-->read next maskInt
if ( ! readValue( maskInt, inMsk, true ) )
throw MDAL_Status::Err_UnknownFormat;
Face f;
for ( int i = 0; i < 3; ++i )
{
int32_t index;
if ( ! readValue( index, inFaces, isNativeLittleEndian ) )
break;
index = correctedIndex( index, superpointIndexes );
f.push_back( static_cast<size_t>( index ) );
}
if ( f.size() == 0 ) //that's mean this is the end of the file
break;
if ( f.size() < 3 ) //that's mean the face is not complete
throw MDAL_Status::Err_UnknownFormat;
//exclude masked face
if ( !( maskInt & 0x01 ) )
faces.push_back( f );
c++;
maskInt = maskInt >> 1;
}
inFaces.close();
inMsk.close();
inMsx.close();
}
std::string MDAL::DriverEsriTin::getCrsWkt( const std::string &uri ) const
{
std::ifstream inCRS( crsFile( uri ), std::ifstream::in );
if ( ! inCRS.is_open() )
return std::string();
std::string crsWkt;
std::getline( inCRS, crsWkt );
if ( crsWkt == "{B286C06B-0879-11D2-AACA-00C04FA33C20}" ) //com class id of the esri UnknownCoordinateSystem class
crsWkt = "";
return crsWkt;
}
int MDAL::DriverEsriTin::correctedIndex( int i, const std::list<int> &superPointsIndexes ) const
{
//usualy i is bigger than the biggest superpoint's index
if ( i > superPointsIndexes.back() )
return i - int( superPointsIndexes.size() ) - 1;
//correction of index depending of the position of the superpoints
int correctedIndex = i - 1;
for ( auto si : superPointsIndexes )
{
if ( i == si )
return -1; //return -1 if the index is associated with superpoint
else if ( i > si )
{
correctedIndex--;
}
else
break;
}
return correctedIndex;
}

100
external/mdal/frmts/mdal_esri_tin.hpp vendored Normal file
View File

@ -0,0 +1,100 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Vincent Cloarec (vcloarec at gmail dot com)
*/
#ifndef MDAL_ESRI_TIN_H
#define MDAL_ESRI_TIN_H
#include <string>
#include <vector>
#include <list>
#include <memory>
#include <iosfwd>
#include <iostream>
#include <fstream>
#include "mdal_data_model.hpp"
#include "mdal.h"
#include "mdal_driver.hpp"
#include "mdal_utils.hpp"
namespace MDAL
{
/** *********************************************
* Structure of files
* https://en.wikipedia.org/wiki/Esri_TIN
************************************************
* thul.adf
*
* file containing only int (big-endianness)
* first, list of superpoint indexes (infinity point to North, South East and West)
* used by ArcGIS but unused here
* This list finished with value -1
* Then indexes of the outer boundary indexes or inner boundary indexes (holes)
* Those list are separated with value 0
*
* ***************************************
* tnxy.adf
*
* vertices position X,Y : for each value 8 bytes : double with big-endianness
*
* ***************************************
* tnz.adf : z value of each vertices
*
* vertices Z value: for each value 4 bytes : float with big endianness
*
* ***************************************
* tnod.adf : faces (triangles) of the TIN
*
* indexes of the 3 vertices for each faces (int32 big-endianness)
*
* ***************************************
* tmsk.adf and tmsx.adf
*
* files used to store information about masked faces. Those faces are not displayed
*
*****************************************
* prj.adf
*
* text file containing the WKT CRS
*
*/
class DriverEsriTin: public Driver
{
public:
DriverEsriTin();
~DriverEsriTin() override {}
Driver *create() override;
virtual std::unique_ptr< Mesh > load( const std::string &uri, MDAL_Status *status ) override;
bool canRead( const std::string &uri ) override;
private:
std::string xyFile( const std::string &uri ) const;
std::string zFile( const std::string &uri ) const;
std::string faceFile( const std::string &uri ) const;
std::string mskFile( const std::string &uri ) const;
std::string msxFile( const std::string &uri ) const;
std::string hullFile( const std::string &uri ) const;
std::string crsFile( const std::string &uri ) const;
void readSuperpoints( const std::string &uri, std::list<int> &superpointsIndexes ) const;
void populateVertices( const std::string &uri, Vertices &vertices, const std::list<int> &superpointIndexes ) const;
void populateFaces( const std::string &uri, Faces &faces, const std::list<int> &superpointIndexes ) const;
std::string getCrsWkt( const std::string &uri ) const;
std::string getTinName( const std::string &uri ) const;
// correction of vertex index because indexing in ESRI TIN files begins with 1 and to take account superpoints position
// return -1 if i is the index of a superpoint
int correctedIndex( int i, const std::list<int> &superPointsIndexes ) const;
};
}
#endif // MDAL_ESRI_TIN_H

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

@ -4,6 +4,7 @@
*/
#include "mdal_hdf5.hpp"
#include <cstring>
HdfFile::HdfFile( const std::string &path )
@ -68,17 +69,21 @@ 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 );
hid_t datatype = H5Aget_type( id() );
char name[HDF_MAX_NAME + 1];
std::memset( name, '\0', HDF_MAX_NAME + 1 );
herr_t status = H5Aread( d->id, datatype, name );
if ( status < 0 )
{
//MDAL::debug("Failed to read data!");
H5Tclose( datatype );
return std::string();
}
H5Tclose( datatype );
return std::string( name );
std::string res( name );
res = MDAL::trim( res );
return res;
}
HdfDataset::HdfDataset( hid_t file, const std::string &path )
@ -189,7 +194,6 @@ std::string HdfDataset::readString() const
return std::string( name );
}
HdfDataspace::HdfDataspace( const std::vector<hsize_t> &dims )
: d( std::make_shared< Handle >( H5Screate_simple(
static_cast<int>( dims.size() ),

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

@ -97,7 +97,6 @@ class HdfGroup
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;
@ -176,6 +175,8 @@ class HdfDataset
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;
inline bool hasAttribute( const std::string &attr_name ) const;
inline HdfAttribute attribute( const std::string &attr_name ) const;
template <typename T> std::vector<T> readArray( hid_t mem_type_id ) const
{
@ -233,10 +234,18 @@ inline HdfGroup HdfGroup::group( const std::string &groupName ) const { return H
inline HdfDataset HdfGroup::dataset( const std::string &dsName ) const { return HdfDataset( file_id(), childPath( dsName ) ); }
inline bool HdfDataset::hasAttribute( const std::string &attr_name ) const
{
htri_t res = H5Aexists( d->id, attr_name.c_str() );
return res > 0 ;
}
inline HdfAttribute HdfFile::attribute( const std::string &attr_name ) const { return HdfAttribute( d->id, attr_name ); }
inline HdfAttribute HdfGroup::attribute( const std::string &attr_name ) const { return HdfAttribute( d->id, attr_name ); }
inline HdfAttribute HdfDataset::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; }

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

@ -43,7 +43,6 @@ static HdfDataset openHdfDataset( const HdfGroup &hdfGroup, const std::string &n
return dsFileType;
}
static std::string openHdfAttribute( const HdfFile &hdfFile, const std::string &name )
{
HdfAttribute attr = hdfFile.attribute( name );
@ -51,6 +50,13 @@ static std::string openHdfAttribute( const HdfFile &hdfFile, const std::string &
return attr.readString();
}
static std::string openHdfAttribute( const HdfDataset &hdfDataset, const std::string &name )
{
HdfAttribute attr = hdfDataset.attribute( name );
if ( !attr.isValid() ) throw MDAL_Status::Err_UnknownFormat;
return attr.readString();
}
static HdfGroup getBaseOutputGroup( const HdfFile &hdfFile )
{
HdfGroup gResults = openHdfGroup( hdfFile, "Results" );
@ -69,12 +75,64 @@ static HdfGroup get2DFlowAreasGroup( const HdfFile &hdfFile, const std::string l
return g2DFlowRes;
}
static std::string getDataTimeUnit( HdfDataset &dsTime )
{
// Initially we expect data to be in hours
std::string dataTimeUnit = "Hours";
// First we look for the Time attribute
if ( dsTime.hasAttribute( "Time" ) )
{
dataTimeUnit = openHdfAttribute( dsTime, "Time" );
return dataTimeUnit;
}
// Some variants of HEC_RAS have time unit stored in Variables attribute
// With read values looking like "Time|Days "
if ( dsTime.hasAttribute( "Variables" ) )
{
dataTimeUnit = openHdfAttribute( dsTime, "Variables" );
// Removing the "Time|" prefix
dataTimeUnit = MDAL::replace( dataTimeUnit, "Time|", "" );
}
return dataTimeUnit;
}
static void convertTimeDataToHours( std::vector<float> &times, const std::string &originalTimeDataUnit )
{
if ( originalTimeDataUnit != "Hours" )
{
for ( size_t i = 0; i < times.size(); i++ )
{
if ( originalTimeDataUnit == "Seconds" ) { times[i] /= 3600.0f; }
else if ( originalTimeDataUnit == "Minutes" ) { times[i] /= 60.0f; }
else if ( originalTimeDataUnit == "Days" ) { times[i] *= 24; }
}
}
}
static std::string readReferenceTime( const HdfFile &hdfFile )
{
std::string refTime;
HdfGroup gBaseO = getBaseOutputGroup( hdfFile );
HdfGroup gUnsteadTS = openHdfGroup( gBaseO, "Unsteady Time Series" );
HdfDataset dsTimeDateStamp = openHdfDataset( gUnsteadTS, "Time Date Stamp" );
std::vector<std::string> timeStamps = dsTimeDateStamp.readArrayString();
if ( timeStamps.size() > 0 )
refTime = timeStamps[0];
return refTime;
}
static std::vector<float> readTimes( const HdfFile &hdfFile )
{
HdfGroup gBaseO = getBaseOutputGroup( hdfFile );
HdfGroup gUnsteadTS = openHdfGroup( gBaseO, "Unsteady Time Series" );
HdfDataset dsTimes = openHdfDataset( gUnsteadTS, "Time" );
std::vector<float> times = dsTimes.readArray();
HdfDataset dsTime = openHdfDataset( gUnsteadTS, "Time" );
std::string dataTimeUnits = getDataTimeUnit( dsTime );
std::vector<float> times = dsTime.readArray();
convertTimeDataToHours( times, dataTimeUnits );
return times;
}
@ -99,7 +157,8 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
const std::vector<std::string> &flowAreaNames,
const std::string rawDatasetName,
const std::string datasetName,
const std::vector<float> &times )
const std::vector<float> &times,
const std::string &referenceTime )
{
double eps = std::numeric_limits<double>::min();
@ -111,6 +170,7 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
);
group->setIsOnVertices( false );
group->setIsScalar( true );
group->setReferenceTime( referenceTime );
std::vector<std::shared_ptr<MDAL::MemoryDataset>> datasets;
@ -177,17 +237,17 @@ void MDAL::DriverHec2D::readFaceResults( const HdfFile &hdfFile,
// UNSTEADY
HdfGroup flowGroup = get2DFlowAreasGroup( hdfFile, "Unsteady Time Series" );
std::vector<float> times = readTimes( hdfFile );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Shear Stress", "Face Shear Stress", times );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Velocity", "Face Velocity", times );
const std::string referenceTime = readReferenceTime( hdfFile );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Shear Stress", "Face Shear Stress", times, referenceTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Velocity", "Face Velocity", times, referenceTime );
// SUMMARY
flowGroup = get2DFlowAreasGroup( hdfFile, "Summary Output" );
times.clear();
times.push_back( 0.0f );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Shear Stress", "Face Shear Stress/Maximums", times );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Velocity", "Face Velocity/Maximums", times );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Shear Stress", "Face Shear Stress/Maximums", times, referenceTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Velocity", "Face Velocity/Maximums", times, referenceTime );
}
@ -197,7 +257,8 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readElemOutput( const Hd
const std::string rawDatasetName,
const std::string datasetName,
const std::vector<float> &times,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation )
std::shared_ptr<MDAL::MemoryDataset> bed_elevation,
const std::string &referenceTime )
{
double eps = std::numeric_limits<double>::min();
@ -209,6 +270,7 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readElemOutput( const Hd
);
group->setIsOnVertices( false );
group->setIsScalar( true );
group->setReferenceTime( referenceTime );
std::vector<std::shared_ptr<MDAL::MemoryDataset>> datasets;
@ -287,6 +349,8 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readBedElevation(
const std::vector<std::string> &flowAreaNames )
{
std::vector<float> times( 1, 0.0f );
std::string referenceTime;
return readElemOutput(
gGeom2DFlowAreas,
areaElemStartIndex,
@ -294,7 +358,8 @@ std::shared_ptr<MDAL::MemoryDataset> MDAL::DriverHec2D::readBedElevation(
"Cells Minimum Elevation",
"Bed Elevation",
times,
std::shared_ptr<MDAL::MemoryDataset>()
std::shared_ptr<MDAL::MemoryDataset>(),
referenceTime
);
}
@ -307,6 +372,7 @@ void MDAL::DriverHec2D::readElemResults(
// UNSTEADY
HdfGroup flowGroup = get2DFlowAreasGroup( hdfFile, "Unsteady Time Series" );
std::vector<float> times = readTimes( hdfFile );
std::string referenceTime = readReferenceTime( hdfFile );
readElemOutput(
flowGroup,
@ -315,7 +381,8 @@ void MDAL::DriverHec2D::readElemResults(
"Water Surface",
"Water Surface",
times,
bed_elevation );
bed_elevation,
referenceTime );
readElemOutput(
flowGroup,
areaElemStartIndex,
@ -323,12 +390,14 @@ void MDAL::DriverHec2D::readElemResults(
"Depth",
"Depth",
times,
bed_elevation );
bed_elevation,
referenceTime );
// SUMMARY
flowGroup = get2DFlowAreasGroup( hdfFile, "Summary Output" );
times.clear();
times.push_back( 0.0f );
referenceTime.clear();
readElemOutput(
flowGroup,
@ -337,7 +406,8 @@ void MDAL::DriverHec2D::readElemResults(
"Maximum Water Surface",
"Water Surface/Maximums",
times,
bed_elevation
bed_elevation,
referenceTime
);
}

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

@ -26,6 +26,14 @@ namespace MDAL
* There is a small change in the format in HEC-RAS 5.0.5+, where
* - Header File Type is different (HEC-RAS Results vs HEC-RAS Geometry)
* - Names or areas are stored in different place (names array vs attributes array)
*
* Time data unit should be present in Time dataset and Time or Variable attribute for given dataset root,
* Since MDAL API is reporting times in float hours, the original values need to be corrected
* based on value found in the Time attribute.
*
* All reference times can be found in Time Data Stamp dataset.
* First value in the dataset is reported by MDAL as reference time
*
*/
class DriverHec2D: public Driver
{
@ -56,7 +64,8 @@ namespace MDAL
const std::vector<std::string> &flowAreaNames,
const std::string rawDatasetName,
const std::string datasetName,
const std::vector<float> &times );
const std::vector<float> &times,
const std::string &referenceTime );
void readFaceResults( const HdfFile &hdfFile,
const std::vector<size_t> &areaElemStartIndex,
@ -69,7 +78,8 @@ namespace MDAL
const std::string rawDatasetName,
const std::string datasetName,
const std::vector<float> &times,
std::shared_ptr<MDAL::MemoryDataset> bed_elevation );
std::shared_ptr<MDAL::MemoryDataset> bed_elevation,
const std::string &referenceTime );
std::shared_ptr<MDAL::MemoryDataset> readBedElevation(
const HdfGroup &gGeom2DFlowAreas,

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

@ -82,6 +82,19 @@ std::vector<std::string> NetCDFFile::readArrNames() const
return res;
}
bool NetCDFFile::hasAttrInt( const std::string &name, const std::string &attr_name ) const
{
assert( mNcid != 0 );
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR ) return false;
int val;
if ( nc_get_att_int( mNcid, arr_id, attr_name.c_str(), &val ) ) return false;
return true;
}
int NetCDFFile::getAttrInt( const std::string &name, const std::string &attr_name ) const
{
assert( mNcid != 0 );
@ -155,6 +168,27 @@ void NetCDFFile::getDimension( const std::string &name, size_t *val, int *ncid_v
if ( nc_inq_dimlen( mNcid, *ncid_val, val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
void NetCDFFile::getDimensions( const std::string &variableName, std::vector<size_t> &dimensions, std::vector<int> &dimensionIds )
{
assert( mNcid != 0 );
int n;
int varId;
if ( nc_inq_varid( mNcid, variableName.c_str(), &varId ) ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_varndims( mNcid, varId, &n ) ) throw MDAL_Status::Err_UnknownFormat;
dimensionIds.resize( size_t( n ) );
dimensions.resize( size_t( n ) );
if ( nc_inq_vardimid( mNcid, varId, dimensionIds.data() ) ) throw MDAL_Status::Err_UnknownFormat;
for ( int i = 0; i < n; ++i )
{
nc_inq_dimlen( mNcid, dimensionIds[size_t( i )], &dimensions[size_t( i )] );
}
}
bool NetCDFFile::hasDimension( const std::string &name ) const
{
int ncid_val;

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

@ -26,6 +26,7 @@ class NetCDFFile
bool hasArr( const std::string &name ) const;
std::vector<std::string> readArrNames() const;
bool hasAttrInt( const std::string &name, const std::string &attr_name ) const;
int getAttrInt( const std::string &name, const std::string &attr_name ) const;
double getAttrDouble( int varid, const std::string &attr_name ) const;
/**
@ -39,6 +40,7 @@ class NetCDFFile
double getFillValue( int varid ) const;
int getVarId( const std::string &name );
void getDimension( const std::string &name, size_t *val, int *ncid_val ) const;
void getDimensions( const std::string &variableName, std::vector<size_t> &dimensionsId, std::vector<int> &dimensionIds );
bool hasDimension( const std::string &name ) const;
private:

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

@ -145,35 +145,18 @@ std::string MDAL::SerafinStreamReader::read_string_without_length( size_t len )
double MDAL::SerafinStreamReader::read_double( )
{
double ret;
unsigned char buffer[8];
if ( mStreamInFloatPrecision )
{
float ret_f;
if ( mIn.read( reinterpret_cast< char * >( &buffer ), 4 ) )
if ( !mIn )
throw MDAL_Status::Err_UnknownFormat;
if ( mIsNativeLittleEndian )
{
std::reverse( reinterpret_cast< char * >( &buffer ), reinterpret_cast< char * >( &buffer ) + 4 );
}
memcpy( reinterpret_cast< char * >( &ret_f ),
reinterpret_cast< char * >( &buffer ),
4 );
if ( !readValue( ret_f, mIn, mIsNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
ret = static_cast<double>( ret_f );
}
else
{
if ( mIn.read( reinterpret_cast< char * >( &buffer ), 8 ) )
if ( !mIn )
throw MDAL_Status::Err_UnknownFormat;
if ( mIsNativeLittleEndian )
{
std::reverse( reinterpret_cast< char * >( &buffer ), reinterpret_cast< char * >( &buffer ) + 8 );
}
memcpy( reinterpret_cast< char * >( &ret ),
reinterpret_cast< char * >( &buffer ),
8 );
if ( !readValue( ret, mIn, mIsNativeLittleEndian ) )
throw MDAL_Status::Err_UnknownFormat;
}
return ret;
}
@ -453,6 +436,20 @@ void MDAL::DriverSelafin::addData( const std::vector<std::string> &var_names, co
var_name = MDAL::replace( var_name, "along y", "" );
}
if ( MDAL::contains( var_name, "vitesse u" ) || MDAL::contains( var_name, "suivant x" ) )
{
is_vector = true;
var_name = MDAL::replace( var_name, "vitesse u", "vitesse" );
var_name = MDAL::replace( var_name, "suivant x", "" );
}
else if ( MDAL::contains( var_name, "vitesse v" ) || MDAL::contains( var_name, "suivant y" ) )
{
is_vector = true;
is_x = false;
var_name = MDAL::replace( var_name, "vitesse v", "vitesse" );
var_name = MDAL::replace( var_name, "suivant y", "" );
}
std::shared_ptr<DatasetGroup> group = mMesh->group( var_name );
if ( !group )
{
@ -601,6 +598,5 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverSelafin::load( const std::string &meshFi
if ( status ) *status = ( error );
mMesh.reset();
}
return std::unique_ptr<Mesh>( mMesh.release() );
}

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

@ -190,6 +190,11 @@ void MDAL::DriverSWW::readDatasetGroups(
for ( const std::string &name : names )
{
// currently we do not support variables like elevation_c, friction_c, stage_c, xmomentum_c, ymomentum_c
// which contain values per volume instead of per vertex
if ( MDAL::endsWith( name, "_c" ) )
continue;
// skip already parsed variables
if ( parsedVariableNames.find( name ) == parsedVariableNames.cend() )
{

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

@ -45,8 +45,21 @@ std::string MDAL::DriverUgrid::findMeshName( int dimension, bool optional ) cons
std::string MDAL::DriverUgrid::nodeZVariableName() const
{
// looks like mesh attributes does not have node_z array name
// reference
const std::vector<std::string> variables = mNcFile.readArrNames();
for ( const std::string &varName : variables )
{
const std::string stdName = mNcFile.getAttrStr( varName, "standard_name" );
const std::string meshName = mNcFile.getAttrStr( varName, "mesh" );
const std::string location = mNcFile.getAttrStr( varName, "location" );
if ( stdName == "altitude" && meshName == mMesh2dName && location == "node" )
{
return varName;
}
}
// not found, the file in non UGRID standard conforming,
// but lets try the common name
return mMesh2dName + "_node_z";
}
@ -61,33 +74,91 @@ MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
// 2D Mesh
// node_dimension is usually something like nMesh2D_node
// number of nodes/vertices in the mesh
const std::string mesh2dNode = mNcFile.getAttrStr( mMesh2dName, "node_dimension" );
mNcFile.getDimension( mesh2dNode, &count, &ncid );
dims.setDimension( CFDimensions::Vertex2D, count, ncid );
//node dimension location is retrieved from the node variable
// face_dimension is usually something like nMesh2D_face
// number of faces in the mesh
const std::string mesh2dFace = mNcFile.getAttrStr( mMesh2dName, "face_dimension" );
mNcFile.getDimension( mesh2dFace, &count, &ncid );
dims.setDimension( CFDimensions::Face2D, count, ncid );
std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile.getAttrStr( mMesh2dName, "node_coordinates" ), ' ' );
if ( nodeVariablesName.size() < 2 )
throw MDAL_Status::Err_UnknownFormat;
// edge_dimension is usually something like nMesh2D_edge
// number of edges in the mesh
std::vector<size_t> nodeDimension;
std::vector<int> nodeDimensionId;
mNcFile.getDimensions( nodeVariablesName.at( 0 ), nodeDimension, nodeDimensionId );
if ( nodeDimension.size() != 1 )
throw MDAL_Status::Err_UnknownFormat;
dims.setDimension( CFDimensions::Vertex2D, nodeDimension.at( 0 ), nodeDimensionId.at( 0 ) );
//face dimension location is retrieved from the face_node_connectivity variable
//if face_dimension is defined as attribute, the dimension at this location help to desambiguate vertex per faces and number of faces
std::string faceConnectivityVariablesName = mNcFile.getAttrStr( mMesh2dName, "face_node_connectivity" );
std::string faceDimensionLocation = mNcFile.getAttrStr( mMesh2dName, "face_dimension" );
if ( faceConnectivityVariablesName == "" )
throw MDAL_Status::Err_UnknownFormat;
size_t facesCount;
size_t maxVerticesPerFace;
std::vector<size_t> faceDimension;
std::vector<int> faceDimensionId;
int facesIndexDimensionId;
int maxVerticesPerFaceDimensionId;
mNcFile.getDimensions( faceConnectivityVariablesName, faceDimension, faceDimensionId );
if ( faceDimension.size() != 2 )
throw MDAL_Status::Err_UnknownFormat;
if ( faceDimensionLocation != "" )
{
mNcFile.getDimension( faceDimensionLocation, &facesCount, &ncid );
if ( facesCount == faceDimension.at( 0 ) )
{
facesIndexDimensionId = faceDimensionId.at( 0 );
maxVerticesPerFaceDimensionId = faceDimensionId.at( 1 );
maxVerticesPerFace = faceDimension.at( 1 );
}
else
{
facesIndexDimensionId = faceDimensionId.at( 1 );
maxVerticesPerFaceDimensionId = faceDimensionId.at( 0 );
maxVerticesPerFace = faceDimension.at( 0 );
}
}
else
{
facesIndexDimensionId = faceDimensionId.at( 0 );
facesCount = faceDimension.at( 0 );
maxVerticesPerFaceDimensionId = faceDimensionId.at( 1 );
maxVerticesPerFace = faceDimension.at( 1 );
}
dims.setDimension( CFDimensions::Face2D, facesCount, facesIndexDimensionId );
dims.setDimension( CFDimensions::MaxVerticesInFace, maxVerticesPerFace, maxVerticesPerFaceDimensionId );
// number of edges in the mesh, not required for UGRID format
const std::string mesh2dEdge = mNcFile.getAttrStr( mMesh2dName, "edge_dimension" );
mNcFile.getDimension( mesh2dEdge, &count, &ncid );
dims.setDimension( CFDimensions::Face2DEdge, count, ncid );
if ( mNcFile.hasDimension( mesh2dEdge ) )
{
mNcFile.getDimension( mesh2dEdge, &count, &ncid );
dims.setDimension( CFDimensions::Face2DEdge, count, ncid );
}
else
{
dims.setDimension( CFDimensions::Face2DEdge, 0, -1 );
}
// max_face_nodes_dimension is usually something like max_nMesh2D_face_nodes
// maximum number of vertices in faces
const std::string mesh2dMaxNodesInFace = mNcFile.getAttrStr( mMesh2dName, "max_face_nodes_dimension" );
mNcFile.getDimension( mesh2dMaxNodesInFace, &count, &ncid );
dims.setDimension( CFDimensions::MaxVerticesInFace, count, ncid );
// Time
mNcFile.getDimension( "time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );
// Time not required for UGRID format
if ( mNcFile.hasDimension( "time" ) )
{
mNcFile.getDimension( "time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );
}
else
{
dims.setDimension( CFDimensions::Time, 0, -1 );
}
return dims;
}
@ -138,7 +209,9 @@ void MDAL::DriverUgrid::populateFaces( MDAL::Faces &faces )
const std::string mesh2dFaceNodeConnectivity = mNcFile.getAttrStr( mMesh2dName, "face_node_connectivity" );
size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
int fill_val = mNcFile.getAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" );
int fill_val = -1;
if ( mNcFile.hasAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" ) )
fill_val = mNcFile.getAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" );
int start_index = mNcFile.getAttrInt( mesh2dFaceNodeConnectivity, "start_index" );
std::vector<int> face_nodes_conn = mNcFile.readIntArr( mesh2dFaceNodeConnectivity, faceCount * verticesInFace );
@ -179,11 +252,16 @@ std::string MDAL::DriverUgrid::getCoordinateSystemVariableName()
std::string coordinate_system_variable;
// first try to get the coordinate system variable from grid definition
if ( mNcFile.hasArr( nodeZVariableName() ) )
std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile.getAttrStr( mMesh2dName, "node_coordinates" ), ' ' );
if ( nodeVariablesName.size() > 1 )
{
coordinate_system_variable = mNcFile.getAttrStr( nodeZVariableName(), "grid_mapping" );
if ( mNcFile.hasArr( nodeVariablesName[0] ) )
{
coordinate_system_variable = mNcFile.getAttrStr( nodeVariablesName[0], "grid_mapping" );
}
}
// if automatic discovery fails, try to check some hardcoded common variables that store projection
if ( coordinate_system_variable.empty() )
{
@ -281,16 +359,18 @@ void MDAL::DriverUgrid::parseNetCDFVariableMetadata( int varid, const std::strin
}
else
{
if ( MDAL::contains( long_name, ", x-component" ) )
if ( MDAL::contains( long_name, ", x-component" ) || MDAL::contains( long_name, "u component of " ) )
{
*is_vector = true;
name = MDAL::replace( long_name, ", x-component", "" );
name = MDAL::replace( name, "u component of ", "" );
}
else if ( MDAL::contains( long_name, ", y-component" ) )
else if ( MDAL::contains( long_name, ", y-component" ) || MDAL::contains( long_name, "v component of " ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( long_name, ", y-component", "" );
name = MDAL::replace( name, "v component of ", "" );
}
else
{

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

@ -247,7 +247,8 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXmdf::readXmdfGroupAsDatasetGrou
bool isVector = dimValues.size() == 3;
std::vector<double> times = dsTimes.readArrayDouble();
std::string timeUnit = rootGroup.attribute( "TimeUnits" ).readString();
convertTimeDataToHours( times, timeUnit );
// all fine, set group and return
group = std::make_shared<MDAL::DatasetGroup>(
name(),
@ -257,6 +258,7 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXmdf::readXmdfGroupAsDatasetGrou
);
group->setIsScalar( !isVector );
group->setIsOnVertices( true );
group->setMetadata( "TIMEUNITS", timeUnit );
// lazy loading of min and max of the dataset group
std::vector<float> mins = dsMins.readArray();
@ -279,3 +281,16 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXmdf::readXmdfGroupAsDatasetGrou
return group;
}
void MDAL::DriverXmdf::convertTimeDataToHours( std::vector<double> &times, const std::string &originalTimeDataUnit )
{
if ( originalTimeDataUnit != "Hours" )
{
for ( size_t i = 0; i < times.size(); i++ )
{
if ( originalTimeDataUnit == "Seconds" ) { times[i] /= 3600.0; }
else if ( originalTimeDataUnit == "Minutes" ) { times[i] /= 60.0; }
else if ( originalTimeDataUnit == "Days" ) { times[i] *= 24; }
}
}
}

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

@ -100,6 +100,8 @@ namespace MDAL
const std::string &nameSuffix,
size_t vertexCount,
size_t faceCount );
void convertTimeDataToHours( std::vector<double> &times, const std::string &originalTimeDataUnit );
};
} // namespace MDAL

58
external/mdal/mdal.cpp vendored
View File

@ -22,7 +22,7 @@ static MDAL_Status sLastStatus;
const char *MDAL_Version()
{
return "0.3.3";
return "0.4.0";
}
MDAL_Status MDAL_LastStatus()
@ -92,6 +92,18 @@ bool MDAL_DR_writeDatasetsCapability( DriverH driver )
return d->hasCapability( MDAL::Capability::WriteDatasets );
}
bool MDAL_DR_SaveMeshCapability( DriverH driver )
{
if ( !driver )
{
sLastStatus = MDAL_Status::Err_MissingDriver;
return false;
}
MDAL::Driver *d = static_cast< MDAL::Driver * >( driver );
return d->hasCapability( MDAL::Capability::SaveMesh );
}
const char *MDAL_DR_longName( DriverH driver )
{
if ( !driver )
@ -143,6 +155,39 @@ MeshH MDAL_LoadMesh( const char *meshFile )
return static_cast< MeshH >( MDAL::DriverManager::instance().load( filename, &sLastStatus ).release() );
}
void MDAL_SaveMesh( MeshH mesh, const char *meshFile, const char *driver )
{
if ( !meshFile )
{
sLastStatus = MDAL_Status::Err_FileNotFound;
return;
}
std::string driverName( driver );
auto d = MDAL::DriverManager::instance().driver( driver );
if ( !d )
{
sLastStatus = MDAL_Status::Err_MissingDriver;
return;
}
if ( !d->hasCapability( MDAL::Capability::SaveMesh ) )
{
sLastStatus = MDAL_Status::Err_MissingDriverCapability;
return;
}
if ( d->faceVerticesMaximumCount() < MDAL_M_faceVerticesMaximumCount( mesh ) )
{
sLastStatus = MDAL_Status::Err_IncompatibleMesh;
return;
}
std::string filename( meshFile );
MDAL::DriverManager::instance().save( static_cast< MDAL::Mesh * >( mesh ), filename, driverName, &sLastStatus );
}
void MDAL_CloseMesh( MeshH mesh )
{
@ -686,6 +731,16 @@ void MDAL_G_closeEditMode( DatasetGroupH group )
}
}
const char *MDAL_G_referenceTime( DatasetGroupH group )
{
if ( !group )
{
sLastStatus = MDAL_Status::Err_IncompatibleDataset;
return EMPTY_STR;
}
MDAL::DatasetGroup *g = static_cast< MDAL::DatasetGroup * >( group );
return _return_str( g->referenceTime() );
}
void MDAL_G_setMetadata( DatasetGroupH group, const char *key, const char *val )
{
@ -867,3 +922,4 @@ void MDAL_D_minimumMaximum( DatasetH dataset, double *min, double *max )
*min = stats.minimum;
*max = stats.maximum;
}

10
external/mdal/mdal_data_model.cpp vendored
View File

@ -155,6 +155,16 @@ void MDAL::DatasetGroup::setStatistics( const Statistics &statistics )
mStatistics = statistics;
}
std::string MDAL::DatasetGroup::referenceTime() const
{
return mReferenceTime;
}
void MDAL::DatasetGroup::setReferenceTime( const std::string &referenceTime )
{
mReferenceTime = referenceTime;
}
MDAL::Mesh *MDAL::DatasetGroup::mesh() const
{
return mParent;

4
external/mdal/mdal_data_model.hpp vendored
View File

@ -109,6 +109,9 @@ namespace MDAL
Statistics statistics() const;
void setStatistics( const Statistics &statistics );
std::string referenceTime() const;
void setReferenceTime( const std::string &referenceTime );
Mesh *mesh() const;
bool isInEditMode() const;
@ -124,6 +127,7 @@ namespace MDAL
bool mIsOnVertices = true;
std::string mUri; // file/uri from where it came
Statistics mStatistics;
std::string mReferenceTime;
};
typedef std::vector<std::shared_ptr<DatasetGroup>> DatasetGroups;

12
external/mdal/mdal_driver_manager.cpp vendored
View File

@ -9,6 +9,7 @@
#include "frmts/mdal_ascii_dat.hpp"
#include "frmts/mdal_binary_dat.hpp"
#include "frmts/mdal_selafin.hpp"
#include "frmts/mdal_esri_tin.hpp"
#include "mdal_utils.hpp"
#ifdef HAVE_HDF5
@ -92,6 +93,15 @@ void MDAL::DriverManager::loadDatasets( Mesh *mesh, const std::string &datasetFi
*status = MDAL_Status::Err_UnknownFormat;
}
void MDAL::DriverManager::save( MDAL::Mesh *mesh, const std::string &uri, const std::string &driverName, MDAL_Status *status ) const
{
auto selectedDriver = driver( driverName );
std::unique_ptr<Driver> drv( selectedDriver->create() );
drv->save( uri, mesh, status );
}
size_t MDAL::DriverManager::driversCount() const
{
return mDrivers.size();
@ -124,6 +134,7 @@ MDAL::DriverManager::DriverManager()
// MESH DRIVERS
mDrivers.push_back( std::make_shared<MDAL::Driver2dm>() );
mDrivers.push_back( std::make_shared<MDAL::DriverSelafin>() );
mDrivers.push_back( std::make_shared<MDAL::DriverEsriTin>() );
#ifdef HAVE_HDF5
mDrivers.push_back( std::make_shared<MDAL::DriverFlo2D>() );
@ -155,3 +166,4 @@ MDAL::DriverManager::DriverManager()
mDrivers.push_back( std::make_shared<MDAL::DriverXdmf>() );
#endif
}

2
external/mdal/mdal_driver_manager.hpp vendored
View File

@ -32,6 +32,8 @@ namespace MDAL
std::unique_ptr< Mesh > load( const std::string &meshFile, MDAL_Status *status ) const;
void loadDatasets( Mesh *mesh, const std::string &datasetFile, MDAL_Status *status ) const;
void save( Mesh *mesh, const std::string &uri, const std::string &driver, MDAL_Status *status ) const;
size_t driversCount() const;
std::shared_ptr<MDAL::Driver> driver( const std::string &driverName ) const;
std::shared_ptr<MDAL::Driver> driver( size_t index ) const;

40
external/mdal/mdal_utils.cpp vendored
View File

@ -586,3 +586,43 @@ bool MDAL::isNativeLittleEndian()
int n = 1;
return ( *( char * )&n == 1 );
}
std::string MDAL::coordinateToString( double coordinate, int precision )
{
std::ostringstream oss;
oss.setf( std::ios::fixed );
if ( fabs( coordinate ) > 180 )
oss.precision( precision ); //seems to not be a geographic coordinate, so 'precision' digits after the digital point
else
oss.precision( 6 + precision ); //could be a geographic coordinate, so 'precision'+6 digits after the digital point
oss << coordinate;
std::string returnString = oss.str();
returnString.back();
//remove unnecessary '0' or '.'
if ( returnString.size() > 0 )
{
while ( '0' == returnString.back() )
{
returnString.pop_back();
}
if ( '.' == returnString.back() )
returnString.pop_back();
}
return returnString;
}
std::string MDAL::doubleToString( double value, int precision )
{
std::ostringstream oss;
oss.precision( precision );
oss << value;
return oss.str();
}

29
external/mdal/mdal_utils.hpp vendored
View File

@ -14,6 +14,9 @@
#include <vector>
#include <stddef.h>
#include <limits>
#include <sstream>
#include <fstream>
#include <algorithm>
#include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp"
@ -66,6 +69,15 @@ namespace MDAL
double toDouble( const std::string &str );
bool toBool( const std::string &str );
//! Returns the string with a adapted format to coordinate
//! precision is the number of digits after the digital point if fabs(value)>180 (seems to not be a geographic coordinate)
//! precision+6 is the number of digits after the digital point if fabs(value)<=180 (could be a geographic coordinate)
std::string coordinateToString( double coordinate, int precision = 2 );
//! Returns a string with scientific format
//! precision is the number of signifiant digits
std::string doubleToString( double value, int precision = 6 );
/**
* Splits by deliminer and skips empty parts.
* Faster than version with std::string
@ -112,10 +124,25 @@ namespace MDAL
// mesh & datasets
//! Adds bed elevatiom dataset group to mesh
void addBedElevationDatasetGroup( MDAL::Mesh *mesh, const Vertices &vertices );
//! Adds face scalar dataset group to mesh
//! Adds altitude dataset group to mesh
void addFaceScalarDatasetGroup( MDAL::Mesh *mesh, const std::vector<double> &values, const std::string &name );
//! Loop through all faces and activate those which has all 4 values on vertices valid
void activateFaces( MDAL::MemoryMesh *mesh, std::shared_ptr<MemoryDataset> dataset );
//! function used to read all of type of value. Option to change the endianness is provided
template<typename T>
bool readValue( T &value, std::ifstream &in, bool changeEndianness = false )
{
char *const p = reinterpret_cast<char *>( &value );
if ( !in.read( p, sizeof( T ) ) )
return false;
if ( changeEndianness )
std::reverse( p, p + sizeof( T ) );
return true;
}
} // namespace MDAL
#endif //MDAL_UTILS_HPP

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

@ -44,6 +44,7 @@ IF (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ascii_dat.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_binary_dat.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_selafin.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_esri_tin.cpp
)
SET(MDAL_LIB_HDRS
@ -57,6 +58,7 @@ IF (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ascii_dat.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_binary_dat.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_selafin.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_esri_tin.hpp
)
IF(HDF5_FOUND)

8
tests/src/core/testqgsmeshlayer.cpp
View File

@ -227,7 +227,7 @@ void TestQgsMeshLayer::test_read_vertex_scalar_dataset()
QgsMeshDatasetIndex ds( 1, i );
const QgsMeshDatasetGroupMetadata meta = dp->datasetGroupMetadata( ds );
if ( meta.extraOptions().count() == 2 )
if ( dp->name() == QStringLiteral( "mesh_memory" ) )
{
QCOMPARE( meta.extraOptions()["description"], QString( "Vertex Scalar Dataset" ) );
QCOMPARE( meta.extraOptions()["meta2"], QString( "best dataset" ) );
@ -270,7 +270,7 @@ void TestQgsMeshLayer::test_read_vertex_vector_dataset()
QgsMeshDatasetIndex ds( 2, i );
const QgsMeshDatasetGroupMetadata meta = dp->datasetGroupMetadata( ds );
if ( meta.extraOptions().count() == 2 )
if ( dp->name() == QStringLiteral( "mesh_memory" ) )
QCOMPARE( meta.extraOptions()["description"], QString( "Vertex Vector Dataset" ) );
QCOMPARE( meta.name(), QString( "VertexVectorDataset" ) );
QVERIFY( !meta.isScalar() );
@ -310,7 +310,7 @@ void TestQgsMeshLayer::test_read_face_scalar_dataset()
QgsMeshDatasetIndex ds( 3, i );
const QgsMeshDatasetGroupMetadata meta = dp->datasetGroupMetadata( ds );
if ( meta.extraOptions().count() == 2 )
if ( dp->name() == QStringLiteral( "mesh_memory" ) )
QCOMPARE( meta.extraOptions()["description"], QString( "Face Scalar Dataset" ) );
QCOMPARE( meta.name(), QString( "FaceScalarDataset" ) );
QVERIFY( meta.isScalar() );
@ -348,7 +348,7 @@ void TestQgsMeshLayer::test_read_face_vector_dataset()
QgsMeshDatasetIndex ds( 4, i );
const QgsMeshDatasetGroupMetadata meta = dp->datasetGroupMetadata( ds );
if ( meta.extraOptions().count() == 2 )
if ( dp->name() == QStringLiteral( "mesh_memory" ) )
QCOMPARE( meta.extraOptions()["description"], QString( "Face Vector Dataset" ) );
QCOMPARE( meta.name(), QString( "FaceVectorDataset" ) );
QVERIFY( !meta.isScalar() );