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

update to MDAL 0.3.2

Browse Source
This commit is contained in:
Peter Petrik 2019-05-03 10:33:13 +02:00
parent e213bde9c1
commit 435b594bbc
20 changed files with 1209 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -20,24 +20,24 @@ MDAL::Driver3Di *MDAL::Driver3Di::create()
return new Driver3Di();
}
MDAL::CFDimensions MDAL::Driver3Di::populateDimensions( const NetCDFFile &ncFile )
MDAL::CFDimensions MDAL::Driver3Di::populateDimensions( )
{
CFDimensions dims;
size_t count;
int ncid;
// 2D Mesh
ncFile.getDimension( "nMesh2D_nodes", &count, &ncid );
mNcFile.getDimension( "nMesh2D_nodes", &count, &ncid );
dims.setDimension( CFDimensions::Face2D, count, ncid );
ncFile.getDimension( "nCorner_Nodes", &count, &ncid );
mNcFile.getDimension( "nCorner_Nodes", &count, &ncid );
dims.setDimension( CFDimensions::MaxVerticesInFace, count, ncid );
// Vertices count is populated later in populateFacesAndVertices
// it is not known from the array dimensions
// Time
ncFile.getDimension( "time", &count, &ncid );
mNcFile.getDimension( "time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );
return dims;
@ -111,7 +111,7 @@ void MDAL::Driver3Di::populateFacesAndVertices( Vertices &vertices, Faces &faces
mDimensions.setDimension( CFDimensions::Vertex2D, vertices.size() );
}
void MDAL::Driver3Di::addBedElevation( MDAL::Mesh *mesh )
void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )
{
assert( mesh );
if ( 0 == mesh->facesCount() )
@ -189,12 +189,6 @@ std::set<std::string> MDAL::Driver3Di::ignoreNetCDFVariables()
return ignore_variables;
}
std::string MDAL::Driver3Di::nameSuffix( MDAL::CFDimensions::Type type )
{
MDAL_UNUSED( type );
return "";
}
void MDAL::Driver3Di::parseNetCDFVariableMetadata( int varid, const std::string &variableName, std::string &name, bool *is_vector, bool *is_x )
{
*is_vector = false;

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

@ -44,12 +44,11 @@ namespace MDAL
Driver3Di *create() override;
private:
CFDimensions populateDimensions( const NetCDFFile &ncFile ) override;
CFDimensions populateDimensions( ) override;
void populateFacesAndVertices( Vertices &vertices, Faces &faces ) override;
void addBedElevation( Mesh *mesh ) override;
void addBedElevation( MemoryMesh *mesh ) override;
std::string getCoordinateSystemVariableName() override;
std::set<std::string> ignoreNetCDFVariables() override;
std::string nameSuffix( CFDimensions::Type type ) override;
void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) override;

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

@ -79,7 +79,6 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
continue;
size_t arr_size = mDimensions.size( mDimensions.type( dimid ) ) * nTimesteps;
std::string suffix = nameSuffix( mDimensions.type( dimid ) );
// Get name, if it is vector and if it is x or y
std::string name;
@ -87,8 +86,6 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
bool is_x = false;
parseNetCDFVariableMetadata( varid, variable_name, name, &is_vector, &is_x );
if ( !suffix.empty() )
name = name + " (" + suffix + ")";
// Add it to the map
auto it = dsinfo_map.find( name );
@ -136,12 +133,12 @@ static void populate_vals( bool is_vector, double *vals, size_t i,
{
if ( is_vector )
{
vals[2 * i] = MDAL::safeValue( vals_x[idx], fill_val_x );
vals[2 * i + 1] = 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 );
}
else
{
vals[i] = MDAL::safeValue( vals_x[idx], fill_val_x );
vals[i] = MDAL::safeValue( vals_x[idx], fill_val_x );
}
}
@ -190,7 +187,7 @@ void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
// read X data
double fill_val_x = mNcFile.getFillValue( dsi.ncid_x );
std::vector<double> vals_x( dsi.arr_size );
std::vector<double> vals_x( dsi.arr_size, std::numeric_limits<double>::quiet_NaN() );
if ( nc_get_var_double( mNcFile.handle(), dsi.ncid_x, vals_x.data() ) ) CF_THROW_ERR;
// read Y data if vector
@ -199,7 +196,7 @@ void MDAL::DriverCF::addDatasetGroups( MDAL::Mesh *mesh, const std::vector<doubl
if ( dsi.is_vector )
{
fill_val_y = mNcFile.getFillValue( dsi.ncid_y );
vals_y.resize( dsi.arr_size );
vals_y.resize( dsi.arr_size, std::numeric_limits<double>::quiet_NaN() );
if ( nc_get_var_double( mNcFile.handle(), dsi.ncid_y, vals_y.data() ) ) CF_THROW_ERR;
}
@ -254,7 +251,8 @@ bool MDAL::DriverCF::canRead( const std::string &uri )
{
NetCDFFile ncFile;
ncFile.openFile( uri );
populateDimensions( ncFile );
mNcFile = ncFile;
populateDimensions( );
}
catch ( MDAL_Status )
{
@ -318,7 +316,7 @@ std::unique_ptr< MDAL::Mesh > MDAL::DriverCF::load( const std::string &fileName,
mNcFile.openFile( mFileName );
// Parse dimensions
mDimensions = populateDimensions( mNcFile );
mDimensions = populateDimensions( );
// Create mMesh
Faces faces;

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

@ -63,6 +63,7 @@ namespace MDAL
//! NetCDF Climate and Forecast (CF) Metadata Conventions
//! http://cfconventions.org
//! and http://ugrid-conventions.github.io/ugrid-conventions/
class DriverCF: public Driver
{
public:
@ -74,12 +75,11 @@ namespace MDAL
std::unique_ptr< Mesh > load( const std::string &fileName, MDAL_Status *status ) override;
protected:
virtual CFDimensions populateDimensions( const NetCDFFile &ncFile ) = 0;
virtual CFDimensions populateDimensions( ) = 0;
virtual void populateFacesAndVertices( Vertices &vertices, Faces &faces ) = 0;
virtual void addBedElevation( MDAL::Mesh *mesh ) = 0;
virtual void addBedElevation( MDAL::MemoryMesh *mesh ) = 0;
virtual std::string getCoordinateSystemVariableName() = 0;
virtual std::set<std::string> ignoreNetCDFVariables() = 0;
virtual std::string nameSuffix( CFDimensions::Type type ) = 0;
virtual void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) = 0;

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

@ -109,6 +109,7 @@ void MDAL::DriverFlo2D::parseCADPTSFile( const std::string &datFileName, std::ve
// CADPTS.DAT - COORDINATES OF CELL CENTERS (ELEM NUM, X, Y)
while ( std::getline( cadptsStream, line ) )
{
line = MDAL::rtrim( line );
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
if ( lineParts.size() != 3 )
{
@ -140,6 +141,7 @@ void MDAL::DriverFlo2D::parseFPLAINFile( std::vector<double> &elevations,
while ( std::getline( fplainStream, line ) )
{
line = MDAL::rtrim( line );
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
if ( lineParts.size() != 7 )
{
@ -220,6 +222,7 @@ void MDAL::DriverFlo2D::parseTIMDEPFile( const std::string &datFileName, const s
while ( std::getline( inStream, line ) )
{
line = MDAL::rtrim( line );
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
if ( lineParts.size() == 1 )
{
@ -301,6 +304,7 @@ void MDAL::DriverFlo2D::parseDEPTHFile( const std::string &datFileName, const st
// DEPTH.OUT - COORDINATES (ELEM NUM, X, Y, MAX DEPTH)
while ( std::getline( depthStream, line ) )
{
line = MDAL::rtrim( line );
if ( vertex_idx == nVertices ) throw MDAL_Status::Err_IncompatibleMesh;
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
@ -348,6 +352,7 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
{
if ( vertex_idx == nVertices ) throw MDAL_Status::Err_IncompatibleMesh;
line = MDAL::rtrim( line );
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
if ( lineParts.size() != 4 )
{
@ -378,6 +383,7 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
{
if ( vertex_idx == nVertices ) throw MDAL_Status::Err_IncompatibleMesh;
line = MDAL::rtrim( line );
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
if ( lineParts.size() != 4 )
{

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

@ -396,47 +396,6 @@ void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
}
}
void MDAL::DriverGdal::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 );
for ( size_t i = 0; i < 4; ++i )
{
const size_t vertexIndex = elem[i];
if ( isScalar )
{
double val = values[vertexIndex];
if ( std::isnan( val ) )
{
active[idx] = 0; //NOT ACTIVE
break;
}
}
else
{
double x = values[2 * vertexIndex];
double y = values[2 * vertexIndex + 1];
if ( std::isnan( x ) || std::isnan( y ) )
{
active[idx] = 0; //NOT ACTIVE
break;
}
}
}
}
}
void MDAL::DriverGdal::addDatasetGroups()
{
// Add dataset to mMesh
@ -465,7 +424,7 @@ void MDAL::DriverGdal::addDatasetGroups()
{
addDataToOutput( raster_bands[i], dataset, is_vector, i == 0 );
}
activateFaces( dataset );
MDAL::activateFaces( mMesh.get(), dataset );
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->datasets.push_back( dataset );
}

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

@ -90,7 +90,6 @@ namespace MDAL
metadata_hash parseMetadata( GDALMajorObjectH gdalBand, const char *pszDomain = nullptr );
void addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x );
bool addSrcProj();
void activateFaces( std::shared_ptr<MemoryDataset> tos );
void addDatasetGroups();
void createMesh();
void parseRasterBands( const GdalDataset *cfGDALDataset );

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

@ -103,13 +103,13 @@ std::string NetCDFFile::getAttrStr( const std::string &name, const std::string &
return getAttrStr( attr_name, arr_id );
}
std::string NetCDFFile::getAttrStr( const std::string &name, int varid ) const
std::string NetCDFFile::getAttrStr( const std::string &attr_name, int varid ) const
{
assert( mNcid != 0 );
size_t attlen = 0;
if ( nc_inq_attlen( mNcid, varid, name.c_str(), &attlen ) )
if ( nc_inq_attlen( mNcid, varid, attr_name.c_str(), &attlen ) )
{
// attribute is missing
return std::string();
@ -117,7 +117,7 @@ std::string NetCDFFile::getAttrStr( const std::string &name, int varid ) const
char *string_attr = static_cast<char *>( malloc( attlen + 1 ) );
if ( nc_get_att_text( mNcid, varid, name.c_str(), string_attr ) ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_get_att_text( mNcid, varid, attr_name.c_str(), string_attr ) ) throw MDAL_Status::Err_UnknownFormat;
string_attr[attlen] = '\0';
std::string res( string_attr );
@ -139,6 +139,7 @@ double NetCDFFile::getAttrDouble( int varid, const std::string &attr_name ) cons
return res;
}
int NetCDFFile::getVarId( const std::string &name )
{
int ncid_val;
@ -153,3 +154,9 @@ void NetCDFFile::getDimension( const std::string &name, size_t *val, int *ncid_v
if ( nc_inq_dimid( mNcid, name.c_str(), ncid_val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_dimlen( mNcid, *ncid_val, val ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
bool NetCDFFile::hasDimension( const std::string &name ) const
{
int ncid_val;
return nc_inq_dimid( mNcid, name.c_str(), &ncid_val ) == NC_NOERR;
}

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

@ -28,11 +28,19 @@ class NetCDFFile
int getAttrInt( const std::string &name, const std::string &attr_name ) const;
double getAttrDouble( int varid, const std::string &attr_name ) const;
/**
* Get string attribute
* \param name name of the variable
* \param attr_name name of the attribute of the variable
* \returns empty string if attribute is missing, else attribute value
*/
std::string getAttrStr( const std::string &name, const std::string &attr_name ) const;
std::string getAttrStr( const std::string &name, int varid ) const;
std::string getAttrStr( const std::string &attr_name, int varid ) const;
double getFillValue( int varid ) const;
int getVarId( const std::string &name );
void getDimension( const std::string &name, size_t *val, int *ncid_val ) const;
bool hasDimension( const std::string &name ) const;
private:
int mNcid; // C handle to the file
};

606
external/mdal/frmts/mdal_selafin.cpp vendored Normal file
View File

@ -0,0 +1,606 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#include <stdio.h>
#include <string.h>
#include <stddef.h>
#include <iosfwd>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cassert>
#include <memory>
#include <algorithm>
#include "mdal_selafin.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
#include <math.h>
MDAL::SerafinStreamReader::SerafinStreamReader() = default;
void MDAL::SerafinStreamReader::initialize( const std::string &fileName )
{
mFileName = fileName;
if ( !MDAL::fileExists( mFileName ) )
{
throw MDAL_Status::Err_FileNotFound;
}
mIn = std::ifstream( mFileName, std::ifstream::in | std::ifstream::binary );
if ( !mIn )
throw MDAL_Status::Err_FileNotFound; // Couldn't open the file
// get length of file:
mIn.seekg( 0, mIn.end );
mFileSize = mIn.tellg();
mIn.seekg( 0, mIn.beg );
mStreamInFloatPrecision = getStreamPrecision();
mIsNativeLittleEndian = MDAL::isNativeLittleEndian();
}
bool MDAL::SerafinStreamReader::getStreamPrecision( )
{
ignore_array_length( );
ignore( 72 );
std::string varType = read_string_without_length( 8 );
bool ret;
if ( varType == "SERAFIN" )
{
ret = true;
}
else if ( varType == "SERAFIND" )
{
ret = false;
}
else
{
throw MDAL_Status::Err_UnknownFormat;
}
ignore_array_length( );
return ret;
}
std::string MDAL::SerafinStreamReader::read_string( size_t len )
{
size_t length = read_sizet();
if ( length != len ) throw MDAL_Status::Err_UnknownFormat;
std::string ret = read_string_without_length( len );
ignore_array_length();
return ret;
}
std::vector<double> MDAL::SerafinStreamReader::read_double_arr( size_t len )
{
size_t length = read_sizet();
if ( mStreamInFloatPrecision )
{
if ( length != len * 4 ) throw MDAL_Status::Err_UnknownFormat;
}
else
{
if ( length != len * 8 ) throw MDAL_Status::Err_UnknownFormat;
}
std::vector<double> ret( len );
for ( size_t i = 0; i < len; ++i )
{
ret[i] = read_double();
}
ignore_array_length();
return ret;
}
std::vector<int> MDAL::SerafinStreamReader::read_int_arr( size_t len )
{
size_t length = read_sizet();
if ( length != len * 4 ) throw MDAL_Status::Err_UnknownFormat;
std::vector<int> ret( len );
for ( size_t i = 0; i < len; ++i )
{
ret[i] = read_int();
}
ignore_array_length();
return ret;
}
std::vector<size_t> MDAL::SerafinStreamReader::read_size_t_arr( size_t len )
{
size_t length = read_sizet();
if ( length != len * 4 ) throw MDAL_Status::Err_UnknownFormat;
std::vector<size_t> ret( len );
for ( size_t i = 0; i < len; ++i )
{
ret[i] = read_sizet();
}
ignore_array_length();
return ret;
}
std::string MDAL::SerafinStreamReader::read_string_without_length( size_t len )
{
std::vector<char> ptr( len );
mIn.read( ptr.data(), static_cast<int>( len ) );
if ( !mIn )
throw MDAL_Status::Err_UnknownFormat;
size_t str_length = 0;
for ( size_t i = len; i > 0; --i )
{
if ( ptr[i - 1] != 0x20 )
{
str_length = i;
break;
}
}
std::string ret( ptr.data(), str_length );
return ret;
}
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 );
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 );
}
return ret;
}
int MDAL::SerafinStreamReader::read_int( )
{
unsigned char data[4];
if ( mIn.read( reinterpret_cast< char * >( &data ), 4 ) )
if ( !mIn )
throw MDAL_Status::Err_UnknownFormat;
if ( mIsNativeLittleEndian )
{
std::reverse( reinterpret_cast< char * >( &data ), reinterpret_cast< char * >( &data ) + 4 );
}
int var;
memcpy( reinterpret_cast< char * >( &var ),
reinterpret_cast< char * >( &data ),
4 );
return var;
}
size_t MDAL::SerafinStreamReader::read_sizet()
{
int var = read_int( );
return static_cast<size_t>( var );
}
size_t MDAL::SerafinStreamReader::remainingBytes()
{
return static_cast<size_t>( mFileSize - mIn.tellg() );
}
void MDAL::SerafinStreamReader::ignore( int len )
{
mIn.ignore( len );
if ( !mIn )
throw MDAL_Status::Err_UnknownFormat;
}
void MDAL::SerafinStreamReader::ignore_array_length( )
{
ignore( 4 );
}
// //////////////////////////////
// DRIVER
// //////////////////////////////
MDAL::DriverSelafin::DriverSelafin():
Driver( "SELAFIN",
"Selafin File",
"*.slf",
Capability::ReadMesh
)
{
}
MDAL::DriverSelafin *MDAL::DriverSelafin::create()
{
return new DriverSelafin();
}
MDAL::DriverSelafin::~DriverSelafin() = default;
void MDAL::DriverSelafin::parseFile( std::vector<std::string> &var_names,
double *xOrigin,
double *yOrigin,
size_t *nElem,
size_t *nPoint,
size_t *nPointsPerElem,
std::vector<size_t> &ikle,
std::vector<double> &x,
std::vector<double> &y,
std::vector<timestep_map> &data )
{
/* 1 record containing the title of the study (72 characters) and a 8 characters
string indicating the type of format (SERAFIN or SERAFIND)
*/
mReader.initialize( mFileName );
/* 1 record containing the two integers NBV(1) and NBV(2) (number of linear
and quadratic variables, NBV(2) with the value of 0 for Telemac, as
quadratic values are not saved so far), numbered from 1 in docs
*/
std::vector<size_t> nbv = mReader.read_size_t_arr( 2 );
/* NBV(1) records containing the names and units of each variab
le (over 32 characters)
*/
for ( size_t i = 0; i < nbv[0]; ++i )
{
var_names.push_back( mReader.read_string( 32 ) );
}
/* 1 record containing the integers table IPARAM (10 integers, of which only
the 6 are currently being used), numbered from 1 in docs
- if IPARAM (3) != 0: the value corresponds to the x-coordinate of the
origin of the mesh,
- if IPARAM (4) != 0: the value corresponds to the y-coordinate of the
origin of the mesh,
- if IPARAM (7) != 0: the value corresponds to the number of planes
on the vertical (3D computation),
- if IPARAM (8) != 0: the value corresponds to the number of
boundary points (in parallel),
- if IPARAM (9) != 0: the value corresponds to the number of interface
points (in parallel),
- if IPARAM(8) or IPARAM(9) != 0: the array IPOBO below is replaced
by the array KNOLG (total initial number of points). All the other
numbers are local to the sub-domain, including IKLE
*/
std::vector<int> params = mReader.read_int_arr( 10 );
*xOrigin = static_cast<double>( params[2] );
*yOrigin = static_cast<double>( params[3] );
if ( params[6] != 0 )
{
// would need additional parsing
throw MDAL_Status::Err_MissingDriver;
}
/*
if IPARAM (10) = 1: a record containing the computation starting date
*/
if ( params[9] == 1 )
{
std::vector<int> datetime = mReader.read_int_arr( 6 );
MDAL_UNUSED( datetime )
}
/* 1 record containing the integers NELEM,NPOIN,NDP,1 (number of
elements, number of points, number of points per element and the value 1)
*/
std::vector<size_t> numbers = mReader.read_size_t_arr( 4 );
*nElem = numbers[0];
*nPoint = numbers[1];
*nPointsPerElem = numbers[2];
/* 1 record containing table IKLE (integer array
of dimension (NDP,NELEM)
which is the connectivity table.
Attention: in TELEMAC-2D, the dimensions of this array are (NELEM,NDP))
*/
ikle = mReader.read_size_t_arr( ( *nElem ) * ( *nPointsPerElem ) );
for ( size_t i = 0; i < ikle.size(); ++i )
{
-- ikle[i]; //numbered from 1
}
/* 1 record containing table IPOBO (integer array of dimension NPOIN); the
value of one element is 0 for an internal point, and
gives the numbering of boundary points for the others
*/
std::vector<int> iPointBoundary = mReader.read_int_arr( *nPoint );
MDAL_UNUSED( iPointBoundary );
/* 1 record containing table X (real array of dimension NPOIN containing the
abscissae of the points)
*/
x = mReader.read_double_arr( *nPoint );
/* 1 record containing table Y (real array of dimension NPOIN containing the
abscissae of the points)
*/
y = mReader.read_double_arr( *nPoint );
/* Next, for each time step, the following are found:
- 1 record containing time T (real),
- NBV(1)+NBV(2) records containing the results tables for each variable at time
*/
data.resize( var_names.size() );
size_t nTimesteps = mReader.remainingBytes() / ( 12 + ( 4 + ( *nPoint ) * 4 + 4 ) * var_names.size() );
for ( size_t nT = 0; nT < nTimesteps; ++nT )
{
std::vector<double> times = mReader.read_double_arr( 1 );
double time = times[0];
for ( size_t i = 0; i < var_names.size(); ++i )
{
timestep_map &datait = data[i];
std::vector<double> datavals = mReader.read_double_arr( *nPoint );
datait[time] = datavals;
}
}
}
void MDAL::DriverSelafin::createMesh(
double xOrigin,
double yOrigin,
size_t nElems,
size_t nPoints,
size_t nPointsPerElem,
std::vector<size_t> &ikle,
std::vector<double> &x,
std::vector<double> &y )
{
Vertices nodes( nPoints );
Vertex *nodesPtr = nodes.data();
for ( size_t n = 0; n < nPoints; ++n, ++nodesPtr )
{
nodesPtr->x = xOrigin + x[n];
nodesPtr->y = yOrigin + y[n];
}
Faces elements( nElems );
for ( size_t e = 0; e < nElems; ++e )
{
if ( nPointsPerElem != 3 )
{
throw MDAL_Status::Err_IncompatibleMesh; //we can add it, but it is uncommon for this format
}
// elemPtr->setId(e);
elements[e].resize( 3 );
for ( size_t p = 0; p < 3; p++ )
{
size_t val = ikle[e * 3 + p];
if ( val > nPoints - 1 )
{
elements[e][p] = 0;
}
else
{
elements[e][p] = val;
}
}
}
mMesh.reset(
new MemoryMesh(
"SELAFIN",
nodes.size(),
elements.size(),
3, //Triangles
computeExtent( nodes ),
mFileName
)
);
mMesh->faces = elements;
mMesh->vertices = nodes;
}
void MDAL::DriverSelafin::addData( const std::vector<std::string> &var_names, const std::vector<timestep_map> &data, size_t nPoints )
{
for ( size_t nName = 0; nName < var_names.size(); ++nName )
{
std::string var_name( var_names[nName] );
var_name = MDAL::toLower( MDAL::trim( var_name ) );
var_name = MDAL::replace( var_name, "/", "" ); // slash is represented as sub-dataset group but we do not have such subdatasets groups in Selafin
bool is_vector = false;
bool is_x = true;
if ( MDAL::contains( var_name, "velocity u" ) || MDAL::contains( var_name, "along x" ) )
{
is_vector = true;
var_name = MDAL::replace( var_name, "velocity u", "velocity" );
var_name = MDAL::replace( var_name, "along x", "" );
}
else if ( MDAL::contains( var_name, "velocity v" ) || MDAL::contains( var_name, "along y" ) )
{
is_vector = true;
is_x = false;
var_name = MDAL::replace( var_name, "velocity v", "velocity" );
var_name = MDAL::replace( var_name, "along y", "" );
}
std::shared_ptr<DatasetGroup> group = mMesh->group( var_name );
if ( !group )
{
group = std::make_shared< DatasetGroup >(
mMesh->driverName(),
mMesh.get(),
mMesh->uri(),
var_name
);
group->setIsOnVertices( true );
group->setIsScalar( !is_vector );
mMesh->datasetGroups.push_back( group );
}
size_t i = 0;
for ( timestep_map::const_iterator it = data[nName].begin(); it != data[nName].end(); ++it, ++i )
{
std::shared_ptr<MDAL::MemoryDataset> dataset;
if ( group->datasets.size() > i )
{
dataset = std::dynamic_pointer_cast<MDAL::MemoryDataset>( group->datasets[i] );
}
else
{
dataset = std::make_shared< MemoryDataset >( group.get() );
dataset->setTime( it->first );
group->datasets.push_back( dataset );
}
double *values = dataset->values();
for ( size_t nP = 0; nP < nPoints; nP++ )
{
double val = it->second.at( nP );
if ( MDAL::equals( val, 0 ) )
{
val = std::numeric_limits<double>::quiet_NaN();
}
if ( is_vector )
{
if ( is_x )
{
values[2 * nP] = val;
}
else
{
values[2 * nP + 1] = val;
}
}
else
{
values[nP] = val;
}
}
}
}
// now activate faces and calculate statistics
for ( auto group : mMesh->datasetGroups )
{
for ( auto dataset : group->datasets )
{
std::shared_ptr<MDAL::MemoryDataset> dts = std::dynamic_pointer_cast<MDAL::MemoryDataset>( dataset );
MDAL::activateFaces( mMesh.get(), dts );
MDAL::Statistics stats = MDAL::calculateStatistics( dataset );
dataset->setStatistics( stats );
}
MDAL::Statistics stats = MDAL::calculateStatistics( group );
group->setStatistics( stats );
}
}
bool MDAL::DriverSelafin::canRead( const std::string &uri )
{
if ( !MDAL::fileExists( uri ) ) return false;
std::ifstream in( uri, std::ifstream::in | std::ifstream::binary );
if ( !in ) return false;
// The first four bytes of the file should contain the values (in hexadecimal): 00 00 00 50.
// This actually indicates the start of a string of length 80 in the file.
// At position 84 in the file, the eight next bytes should read (in hexadecimal): 00 00 00 50 00 00 00 04.
unsigned char data[ 92 ];
in.read( reinterpret_cast< char * >( &data ), 92 );
if ( !in )
return false;
if ( data[0] != 0 || data[1] != 0 ||
data[2] != 0 || data[3] != 0x50 )
return false;
if ( data[84 + 0] != 0 || data[84 + 1] != 0 ||
data[84 + 2] != 0 || data[84 + 3] != 0x50 ||
data[84 + 4] != 0 || data[84 + 5] != 0 ||
data[84 + 6] != 0 || data[84 + 7] != 8 )
return false;
return true;
}
std::unique_ptr<MDAL::Mesh> MDAL::DriverSelafin::load( const std::string &meshFile, MDAL_Status *status )
{
if ( status ) *status = MDAL_Status::None;
mFileName = meshFile;
mMesh.reset();
std::vector<std::string> var_names;
double xOrigin;
double yOrigin;
size_t nElems;
size_t nPoints;
size_t nPointsPerElem;
std::vector<size_t> ikle;
std::vector<double> x;
std::vector<double> y;
std::vector<timestep_map> data;
try
{
parseFile( var_names,
&xOrigin,
&yOrigin,
&nElems,
&nPoints,
&nPointsPerElem,
ikle,
x,
y,
data );
createMesh( xOrigin,
yOrigin,
nElems,
nPoints,
nPointsPerElem,
ikle,
x,
y );
addData( var_names, data, nPoints );
}
catch ( MDAL_Status error )
{
if ( status ) *status = ( error );
mMesh.reset();
}
return std::unique_ptr<Mesh>( mMesh.release() );
}

99
external/mdal/frmts/mdal_selafin.hpp vendored Normal file
View File

@ -0,0 +1,99 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#ifndef MDAL_SELAFIN_HPP
#define MDAL_SELAFIN_HPP
#include <string>
#include <memory>
#include <map>
#include <iostream>
#include <fstream>
#include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp"
#include "mdal.h"
#include "mdal_driver.hpp"
namespace MDAL
{
class SerafinStreamReader
{
public:
SerafinStreamReader();
void initialize( const std::string &fileName );
std::string read_string( size_t len );
std::vector<double> read_double_arr( size_t len );
std::vector<int> read_int_arr( size_t len );
std::vector<size_t> read_size_t_arr( size_t len );
double read_double( );
int read_int( );
size_t read_sizet( );
size_t remainingBytes();
private:
void ignore_array_length( );
std::string read_string_without_length( size_t len );
void ignore( int len );
bool getStreamPrecision();
std::string mFileName;
bool mStreamInFloatPrecision = true;
bool mIsNativeLittleEndian = true;
long long mFileSize = -1;
std::ifstream mIn;
};
/**
* Serafin format (also called Selafin)
*
* Binary format for triangular mesh with datasets defined on vertices
* http://www.opentelemac.org/downloads/Archive/v6p0/telemac2d_user_manual_v6p0.pdf Appendix 3
* https://www.gdal.org/drv_selafin.html
*/
class DriverSelafin: public Driver
{
public:
DriverSelafin();
~DriverSelafin() override;
DriverSelafin *create() override;
bool canRead( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile, MDAL_Status *status ) override;
private:
typedef std::map<double, std::vector<double> > timestep_map; //TIME (sorted), nodeVal
void createMesh( double xOrigin,
double yOrigin,
size_t nElems,
size_t nPoints,
size_t nPointsPerElem,
std::vector<size_t> &ikle,
std::vector<double> &x,
std::vector<double> &y );
void addData( const std::vector<std::string> &var_names, const std::vector<timestep_map> &data, size_t nPoints );
void parseFile( std::vector<std::string> &var_names,
double *xOrigin,
double *yOrigin,
size_t *nElem,
size_t *nPoint,
size_t *nPointsPerElem,
std::vector<size_t> &ikle,
std::vector<double> &x,
std::vector<double> &y,
std::vector<timestep_map> &data );
bool getStreamPrecision( std::ifstream &in );
std::unique_ptr< MDAL::MemoryMesh > mMesh;
std::string mFileName;
SerafinStreamReader mReader;
};
} // namespace MDAL
#endif //MDAL_SELAFIN_HPP

323
external/mdal/frmts/mdal_ugrid.cpp vendored Normal file
View File

@ -0,0 +1,323 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#include "mdal_ugrid.hpp"
#include "mdal_utils.hpp"
#include <netcdf.h>
#include <assert.h>
MDAL::DriverUgrid::DriverUgrid()
: DriverCF(
"Ugrid",
"UGRID Results",
"*.nc" )
{
}
MDAL::DriverUgrid *MDAL::DriverUgrid::create()
{
return new DriverUgrid();
}
std::string MDAL::DriverUgrid::findMeshName( int dimension, bool optional ) const
{
const std::vector<std::string> variables = mNcFile.readArrNames();
for ( const std::string &varName : variables )
{
const std::string cf_role = mNcFile.getAttrStr( varName, "cf_role" );
if ( cf_role == "mesh_topology" )
{
int topology_dimension = mNcFile.getAttrInt( varName, "topology_dimension" );
if ( topology_dimension == dimension )
{
return varName;
}
}
}
if ( optional )
return "";
else
throw MDAL_Status::Err_UnknownFormat;
}
std::string MDAL::DriverUgrid::nodeZVariableName() const
{
// looks like mesh attributes does not have node_z array name
// reference
return mMesh2dName + "_node_z";
}
MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
{
CFDimensions dims;
size_t count;
int ncid;
mMesh1dName = findMeshName( 1, true ); // optional, may not be present
mMesh2dName = findMeshName( 2, false ); // force
// 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 );
// 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 );
// edge_dimension is usually something like nMesh2D_edge
// number of edges in the mesh
const std::string mesh2dEdge = mNcFile.getAttrStr( mMesh2dName, "edge_dimension" );
mNcFile.getDimension( mesh2dEdge, &count, &ncid );
dims.setDimension( CFDimensions::Face2DEdge, count, ncid );
// 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 );
return dims;
}
void MDAL::DriverUgrid::populateFacesAndVertices( Vertices &vertices, Faces &faces )
{
populateVertices( vertices );
populateFaces( faces );
}
void MDAL::DriverUgrid::populateVertices( MDAL::Vertices &vertices )
{
assert( vertices.empty() );
size_t vertexCount = mDimensions.size( CFDimensions::Vertex2D );
vertices.resize( vertexCount );
Vertex *vertexPtr = vertices.data();
// Parse 2D Mesh
// node_coordinates should be something like Mesh2D_node_x Mesh2D_node_y
std::string verticesXName, verticesYName;
parse2VariablesFromAttribute( mMesh2dName, "node_coordinates", verticesXName, verticesYName, false );
const std::vector<double> vertices2D_x = mNcFile.readDoubleArr( verticesXName, vertexCount );
const std::vector<double> vertices2D_y = mNcFile.readDoubleArr( verticesYName, vertexCount );
std::vector<double> vertices2D_z;
if ( mNcFile.hasArr( nodeZVariableName() ) )
{
vertices2D_z = mNcFile.readDoubleArr( nodeZVariableName(), vertexCount );
}
for ( size_t i = 0; i < vertexCount; ++i, ++vertexPtr )
{
vertexPtr->x = vertices2D_x[i];
vertexPtr->y = vertices2D_y[i];
if ( !vertices2D_z.empty() )
vertexPtr->z = vertices2D_z[i];
}
}
void MDAL::DriverUgrid::populateFaces( MDAL::Faces &faces )
{
assert( faces.empty() );
size_t faceCount = mDimensions.size( CFDimensions::Face2D );
faces.resize( faceCount );
// Parse 2D Mesh
// face_node_connectivity is usually something like Mesh2D_face_nodes
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 start_index = mNcFile.getAttrInt( mesh2dFaceNodeConnectivity, "start_index" );
std::vector<int> face_nodes_conn = mNcFile.readIntArr( mesh2dFaceNodeConnectivity, faceCount * verticesInFace );
for ( size_t i = 0; i < faceCount; ++i )
{
size_t nVertices = verticesInFace;
std::vector<size_t> idxs;
for ( size_t j = 0; j < verticesInFace; ++j )
{
size_t idx = verticesInFace * i + j;
int val = face_nodes_conn[idx];
if ( fill_val == val )
{
// found fill val
nVertices = j;
assert( nVertices > 1 );
break;
}
else
{
idxs.push_back( static_cast<size_t>( val - start_index ) );
}
}
faces[i] = idxs;
}
}
void MDAL::DriverUgrid::addBedElevation( MDAL::MemoryMesh *mesh )
{
if ( mNcFile.hasArr( nodeZVariableName() ) ) MDAL::addBedElevationDatasetGroup( mesh, mesh->vertices );
}
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() ) )
{
coordinate_system_variable = mNcFile.getAttrStr( nodeZVariableName(), "grid_mapping" );
}
// if automatic discovery fails, try to check some hardcoded common variables that store projection
if ( coordinate_system_variable.empty() )
{
if ( mNcFile.hasArr( "projected_coordinate_system" ) )
coordinate_system_variable = "projected_coordinate_system";
else if ( mNcFile.hasArr( "wgs84" ) )
coordinate_system_variable = "wgs84";
}
// return, may be empty
return coordinate_system_variable;
}
std::set<std::string> MDAL::DriverUgrid::ignoreNetCDFVariables()
{
std::set<std::string> ignore_variables;
ignore_variables.insert( "projected_coordinate_system" );
ignore_variables.insert( "time" );
ignore_variables.insert( "timestep" );
std::vector<std::string> meshes;
if ( mNcFile.hasArr( mMesh1dName ) )
meshes.push_back( mMesh1dName );
meshes.push_back( mMesh2dName );
for ( const std::string &mesh : meshes )
{
std::string xName, yName;
ignore_variables.insert( mesh );
parse2VariablesFromAttribute( mesh, "node_coordinates", xName, yName, true );
ignore_variables.insert( xName );
ignore_variables.insert( yName );
ignore_variables.insert( mNcFile.getAttrStr( mesh, "edge_node_connectivity" ) );
parse2VariablesFromAttribute( mesh, "edge_coordinates", xName, yName, true );
if ( !xName.empty() )
{
ignore_variables.insert( xName );
ignore_variables.insert( mNcFile.getAttrStr( xName, "bounds" ) );
}
if ( !yName.empty() )
{
ignore_variables.insert( yName );
ignore_variables.insert( mNcFile.getAttrStr( yName, "bounds" ) );
}
ignore_variables.insert( mNcFile.getAttrStr( mesh, "face_node_connectivity" ) );
parse2VariablesFromAttribute( mesh, "face_coordinates", xName, yName, true );
if ( !xName.empty() )
{
ignore_variables.insert( xName );
ignore_variables.insert( mNcFile.getAttrStr( xName, "bounds" ) );
}
if ( !yName.empty() )
{
ignore_variables.insert( yName );
ignore_variables.insert( mNcFile.getAttrStr( yName, "bounds" ) );
}
ignore_variables.insert( mNcFile.getAttrStr( mesh, "edge_face_connectivity" ) );
}
return ignore_variables;
}
void MDAL::DriverUgrid::parseNetCDFVariableMetadata( int varid, const std::string &variableName, std::string &name, bool *is_vector, bool *is_x )
{
*is_vector = false;
*is_x = true;
std::string long_name = mNcFile.getAttrStr( "long_name", varid );
if ( long_name.empty() )
{
std::string standard_name = mNcFile.getAttrStr( "standard_name", varid );
if ( standard_name.empty() )
{
name = variableName;
}
else
{
if ( MDAL::contains( standard_name, "_x_" ) )
{
*is_vector = true;
name = MDAL::replace( standard_name, "_x_", "" );
}
else if ( MDAL::contains( standard_name, "_y_" ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( standard_name, "_y_", "" );
}
else
{
name = standard_name;
}
}
}
else
{
if ( MDAL::contains( long_name, ", x-component" ) )
{
*is_vector = true;
name = MDAL::replace( long_name, ", x-component", "" );
}
else if ( MDAL::contains( long_name, ", y-component" ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( long_name, ", y-component", "" );
}
else
{
name = long_name;
}
}
}
void MDAL::DriverUgrid::parse2VariablesFromAttribute( const std::string &name, const std::string &attr_name,
std::string &var1, std::string &var2, bool optional ) const
{
const std::string mesh2dNodeCoordinates = mNcFile.getAttrStr( name, attr_name );
const std::vector<std::string> chunks = MDAL::split( mesh2dNodeCoordinates, ' ' );
if ( chunks.size() != 2 )
{
if ( optional )
{
var1 = "";
var2 = "";
}
else
throw MDAL_Status::Err_UnknownFormat;
}
else
{
var1 = chunks[0];
var2 = chunks[1];
}
}

52
external/mdal/frmts/mdal_ugrid.hpp vendored Normal file
View File

@ -0,0 +1,52 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
*/
#ifndef MDAL_UGRID_HPP
#define MDAL_UGRID_HPP
#include <map>
#include <string>
#include <stddef.h>
#include "mdal_cf.hpp"
#include "mdal_driver.hpp"
namespace MDAL
{
/**
* Driver of UGRID file format.
*
* The result UGRID NetCDF file is strictly based on CF-conventions 1.6
*/
class DriverUgrid: public DriverCF
{
public:
DriverUgrid();
~DriverUgrid() override = default;
DriverUgrid *create() override;
private:
CFDimensions populateDimensions( ) override;
void populateFacesAndVertices( Vertices &vertices, Faces &faces ) override;
void populateVertices( Vertices &vertices );
void populateFaces( Faces &faces );
void addBedElevation( MemoryMesh *mesh ) override;
std::string getCoordinateSystemVariableName() override;
std::set<std::string> ignoreNetCDFVariables() override;
void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) override;
void parse2VariablesFromAttribute( const std::string &name, const std::string &attr_name,
std::string &var1, std::string &var2,
bool optional ) const;
std::string findMeshName( int dimension, bool optional ) const;
std::string mMesh2dName;
std::string mMesh1dName;
std::string nodeZVariableName() const;
};
} // namespace MDAL
#endif // MDAL_UGRID_HPP

2
external/mdal/mdal.cpp vendored
View File

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

30
external/mdal/mdal_data_model.cpp vendored
View File

@ -224,6 +224,16 @@ std::string MDAL::Mesh::driverName() const
MDAL::Mesh::~Mesh() = default;
std::shared_ptr<MDAL::DatasetGroup> MDAL::Mesh::group( const std::string &name )
{
for ( auto grp : datasetGroups )
{
if ( grp->name() == name )
return grp;
}
return std::shared_ptr<MDAL::DatasetGroup>();
}
void MDAL::Mesh::setSourceCrs( const std::string &str )
{
mCrs = MDAL::trim( str );
@ -239,26 +249,6 @@ void MDAL::Mesh::setSourceCrsFromEPSG( int code )
setSourceCrs( std::string( "EPSG:" ) + std::to_string( code ) );
}
void MDAL::Mesh::setExtent( const BBox &extent )
{
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;

8
external/mdal/mdal_data_model.hpp vendored
View File

@ -164,16 +164,14 @@ namespace MDAL
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;
//! Find a dataset group by name
std::shared_ptr<DatasetGroup> group( const std::string &name );
size_t verticesCount() const;
size_t facesCount() const;
std::string uri() const;

4
external/mdal/mdal_driver_manager.cpp vendored
View File

@ -8,6 +8,7 @@
#include "frmts/mdal_2dm.hpp"
#include "frmts/mdal_ascii_dat.hpp"
#include "frmts/mdal_binary_dat.hpp"
#include "frmts/mdal_selafin.hpp"
#include "mdal_utils.hpp"
#ifdef HAVE_HDF5
@ -21,6 +22,7 @@
#endif
#ifdef HAVE_NETCDF
#include "frmts/mdal_ugrid.hpp"
#include "frmts/mdal_3di.hpp"
#include "frmts/mdal_sww.hpp"
#endif
@ -121,6 +123,7 @@ MDAL::DriverManager::DriverManager()
{
// MESH DRIVERS
mDrivers.push_back( std::make_shared<MDAL::Driver2dm>() );
mDrivers.push_back( std::make_shared<MDAL::DriverSelafin>() );
#ifdef HAVE_HDF5
mDrivers.push_back( std::make_shared<MDAL::DriverFlo2D>() );
@ -130,6 +133,7 @@ MDAL::DriverManager::DriverManager()
#ifdef HAVE_NETCDF
mDrivers.push_back( std::make_shared<MDAL::Driver3Di>() );
mDrivers.push_back( std::make_shared<MDAL::DriverSWW>() );
mDrivers.push_back( std::make_shared<MDAL::DriverUgrid>() );
#endif
#if defined HAVE_GDAL && defined HAVE_NETCDF

65
external/mdal/mdal_utils.cpp vendored
View File

@ -323,14 +323,16 @@ bool MDAL::equals( double val1, double val2, double eps )
double MDAL::safeValue( double val, double nodata, double eps )
{
if ( equals( val, nodata, eps ) )
{
return std::numeric_limits<double>::quiet_NaN();
}
else
{
if ( std::isnan( val ) )
return val;
}
if ( std::isnan( nodata ) )
return val;
if ( equals( val, nodata, eps ) )
return std::numeric_limits<double>::quiet_NaN();
return val;
}
double MDAL::parseTimeUnits( const std::string &units )
@ -541,3 +543,52 @@ void MDAL::addFaceScalarDatasetGroup( MDAL::Mesh *mesh,
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
}
void MDAL::activateFaces( MDAL::MemoryMesh *mesh, std::shared_ptr<MemoryDataset> dataset )
{
// only for data on vertices
if ( !dataset->group()->isOnVertices() )
return;
bool isScalar = dataset->group()->isScalar();
// Activate only Faces that do all Vertex's outputs with some data
int *active = dataset->active();
const double *values = dataset->constValues();
const size_t nFaces = mesh->facesCount();
for ( unsigned int idx = 0; idx < nFaces; ++idx )
{
Face elem = mesh->faces.at( idx );
for ( size_t i = 0; i < elem.size(); ++i )
{
const size_t vertexIndex = elem[i];
if ( isScalar )
{
double val = values[vertexIndex];
if ( std::isnan( val ) )
{
active[idx] = 0; //NOT ACTIVE
break;
}
}
else
{
double x = values[2 * vertexIndex];
double y = values[2 * vertexIndex + 1];
if ( std::isnan( x ) || std::isnan( y ) )
{
active[idx] = 0; //NOT ACTIVE
break;
}
}
}
}
}
bool MDAL::isNativeLittleEndian()
{
// https://stackoverflow.com/a/4181991/2838364
int n = 1;
return ( *( char * )&n == 1 );
}

5
external/mdal/mdal_utils.hpp vendored
View File

@ -24,6 +24,9 @@
namespace MDAL
{
// endianness
bool isNativeLittleEndian();
// numbers
bool equals( double val1, double val2, double eps = std::numeric_limits<double>::epsilon() );
@ -112,6 +115,8 @@ namespace MDAL
void addBedElevationDatasetGroup( MDAL::Mesh *mesh, const Vertices &vertices );
//! Adds face scalar 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 );
} // namespace MDAL
#endif //MDAL_UTILS_HPP

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

@ -42,6 +42,7 @@ IF (WITH_INTERNAL_MDAL)
${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
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_selafin.cpp
)
SET(MDAL_LIB_HDRS
@ -54,6 +55,7 @@ IF (WITH_INTERNAL_MDAL)
${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
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_selafin.hpp
)
IF(HDF5_FOUND)
@ -91,12 +93,14 @@ IF (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_3di.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_netcdf.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_sww.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ugrid.cpp
)
SET(MDAL_LIB_HDRS ${MDAL_LIB_HDRS}
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_cf.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_3di.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_netcdf.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_sww.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ugrid.hpp
)
SET (HAVE_NETCDF TRUE)
ENDIF(NETCDF_FOUND)