Skip to content

Commit

Permalink
correctly read netcdf arrays with double/float/byte type
Browse files Browse the repository at this point in the history
  • Loading branch information
@PeterPetrik
PeterPetrik committed Nov 29, 2019
1 parent c579ebd commit 707a0b9
Show file tree
Hide file tree
Showing 4 changed files with 155 additions and 26 deletions.
69 changes: 66 additions & 3 deletions mdal/frmts/mdal_netcdf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
#include <vector>
#include <assert.h>
#include <netcdf.h>
#include <cmath>

#include "mdal_netcdf.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
Expand Down Expand Up @@ -82,7 +84,32 @@ std::vector<double> NetCDFFile::readDoubleArr( const std::string &name, size_t d
int arr_id;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
std::vector<double> arr_val( dim );
if ( nc_get_var_double( mNcid, arr_id, arr_val.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;

nc_type typep;
if ( nc_inq_varid( mNcid, name.c_str(), &arr_id ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_inq_vartype( mNcid, arr_id, &typep ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;

if ( typep == NC_FLOAT )
{
std::vector<float> arr_val_f( dim );
if ( nc_get_var_float( mNcid, arr_id, arr_val_f.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
for ( size_t i = 0; i < dim; ++i )
{
const float val = arr_val_f[i];
if ( std::isnan( val ) )
arr_val[i] = std::numeric_limits<double>::quiet_NaN();
else
arr_val[i] = static_cast<double>( val );
}
}
else if ( typep == NC_DOUBLE )
{
if ( nc_get_var_double( mNcid, arr_id, arr_val.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
else
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}

Expand All @@ -97,8 +124,44 @@ std::vector<double> NetCDFFile::readDoubleArr( int arr_id,
const std::vector<ptrdiff_t> stridep = {1, 1};

std::vector<double> arr_val( count_dim1 * count_dim2 );
int res = nc_get_vars_double( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val.data() );
if ( res != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;

nc_type typep;
if ( nc_inq_vartype( mNcid, arr_id, &typep ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;

if ( typep == NC_FLOAT )
{
std::vector<float> arr_val_f( count_dim1 * count_dim2 );
if ( nc_get_vars_float( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val_f.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
for ( size_t i = 0; i < count_dim1 * count_dim2; ++i )
{
const float val = arr_val_f[i];
if ( std::isnan( val ) )
arr_val[i] = std::numeric_limits<double>::quiet_NaN();
else
arr_val[i] = static_cast<double>( val );
}
}
else if ( typep == NC_BYTE )
{
std::vector<unsigned char> arr_val_b( count_dim1 * count_dim2 );
if ( nc_get_vars_ubyte( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val_b.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
for ( size_t i = 0; i < count_dim1 * count_dim2; ++i )
{
const unsigned char val = arr_val_b[i];
if ( val == 129 )
arr_val[i] = std::numeric_limits<double>::quiet_NaN();
else
arr_val[i] = double( int( val ) );
}
}
else if ( typep == NC_DOUBLE )
{
if ( nc_get_vars_double( mNcid, arr_id, startp.data(), countp.data(), stridep.data(), arr_val.data() ) != NC_NOERR ) throw MDAL_Status::Err_UnknownFormat;
}
else
{
throw MDAL_Status::Err_UnknownFormat;
}
return arr_val;
}

Expand Down
39 changes: 22 additions & 17 deletions mdal/frmts/mdal_tuflowfv.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,11 @@ size_t MDAL::TuflowFVDataset3D::faceToVolumeData( size_t indexStart, size_t coun
indexStart,
copyValues
);

// indexed from 1 in FV, from 0 in MDAL
for ( auto &element : vals )
element -= 1;

memcpy( buffer, vals.data(), copyValues * sizeof( int ) );
return copyValues;
}
Expand Down Expand Up @@ -234,6 +239,7 @@ void MDAL::DriverTuflowFV::populateFaces( MDAL::Faces &faces )
{
assert( faces.empty() );
size_t faceCount = mDimensions.size( CFDimensions::Face2D );
size_t vertexCount = mDimensions.size( CFDimensions::Vertex2D );
faces.resize( faceCount );

// Parse 2D Mesh
Expand All @@ -249,8 +255,9 @@ void MDAL::DriverTuflowFV::populateFaces( MDAL::Faces &faces )
for ( size_t j = 0; j < nVertices; ++j )
{
size_t idx = verticesInFace * i + j;
int val = face_nodes_conn[idx];
idxs.push_back( static_cast<size_t>( val ) );
size_t val = static_cast<size_t>( face_nodes_conn[idx] - 1 ); //indexed from 1
assert( val < vertexCount );
idxs.push_back( val );
}
faces[i] = idxs;
}
Expand Down Expand Up @@ -325,24 +332,22 @@ void MDAL::DriverTuflowFV::parseNetCDFVariableMetadata( int varid, const std::st
std::string long_name = mNcFile->getAttrStr( "long_name", varid );
if ( long_name.empty() )
{
if ( MDAL::endsWith( variableName, "_x" ) )
{
*is_vector = true;
name = MDAL::replace( variableName, "_x", "" );
}
else if ( MDAL::endsWith( variableName, "_y" ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( variableName, "_y", "" );
}
else
{
name = variableName;
}
name = variableName;
}
else
{
if ( MDAL::startsWith( long_name, "maximum value of " ) )
long_name = MDAL::replace( long_name, "maximum value of ", "" ) + "/Maximums";

if ( MDAL::startsWith( long_name, "minimum value of " ) )
long_name = MDAL::replace( long_name, "minimum value of ", "" ) + "/Minimums";

if ( MDAL::startsWith( long_name, "time at maximum value of " ) )
long_name = MDAL::replace( long_name, "time at maximum value of ", "" ) + "/Time at Maximums";

if ( MDAL::startsWith( long_name, "time at minimum value of " ) )
long_name = MDAL::replace( long_name, "time at minimum value of ", "" ) + "/Time at Minimums";

if ( MDAL::startsWith( long_name, "x_" ) )
{
*is_vector = true;
Expand Down
71 changes: 66 additions & 5 deletions tests/test_tuflowfv.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,13 +34,13 @@ TEST( MeshTuflowFVTest, TrapSteady053D )
int f_v_count = getFaceVerticesCountAt( m, 1 );
EXPECT_EQ( 4, f_v_count ); //quad
int f_v = getFaceVerticesIndexAt( m, 1, 0 );
EXPECT_EQ( 5, f_v );
EXPECT_EQ( 4, f_v );
f_v = getFaceVerticesIndexAt( m, 1, 1 );
EXPECT_EQ( 2, f_v );
EXPECT_EQ( 1, f_v );
f_v = getFaceVerticesIndexAt( m, 1, 2 );
EXPECT_EQ( 4, f_v );
EXPECT_EQ( 3, f_v );
f_v = getFaceVerticesIndexAt( m, 1, 3 );
EXPECT_EQ( 8, f_v );
EXPECT_EQ( 7, f_v );

ASSERT_EQ( 5, MDAL_M_datasetGroupCount( m ) );

Expand Down Expand Up @@ -135,7 +135,7 @@ TEST( MeshTuflowFVTest, TrapSteady053D )
EXPECT_DOUBLE_EQ( 0, getValueY( ds, 0 ) );

int index3d = get3DFrom2D( ds, 155 );
EXPECT_EQ( 1551, index3d );
EXPECT_EQ( 1550, index3d );

int levelsCount = getLevelsCount3D( ds, 0 );
EXPECT_EQ( 10, levelsCount );
Expand Down Expand Up @@ -272,7 +272,68 @@ TEST( MeshTuflowFVTest, TrapSteady053D )

// Close mesh
MDAL_CloseMesh( m );
}

TEST( MeshTuflowFVTest, TrapSteady053DWithMaxes )
{
std::string path = test_file( "/tuflowfv/withMaxes/trap_steady_05_3D.nc" );
MeshH m = MDAL_LoadMesh( path.c_str() );
EXPECT_NE( m, nullptr );
MDAL_Status s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

// ///////////
// Vertices
// ///////////
int v_count = MDAL_M_vertexCount( m );
EXPECT_EQ( v_count, 369 );

// ///////////
// Faces
// ///////////
int f_count = MDAL_M_faceCount( m );
EXPECT_EQ( 320, f_count );

ASSERT_EQ( 21, MDAL_M_datasetGroupCount( m ) );

// /////////////////////////////////
// Dataset: maximums dataset
// /////////////////////////////////
{
DatasetGroupH g = MDAL_M_datasetGroup( m, 2 );
ASSERT_NE( g, nullptr );

int meta_count = MDAL_G_metadataCount( g );
ASSERT_EQ( 1, meta_count );

const char *name = MDAL_G_name( g );
EXPECT_EQ( std::string( "temperature/Maximums" ), std::string( name ) );

bool scalar = MDAL_G_hasScalarData( g );
EXPECT_EQ( true, scalar );

MDAL_DataLocation dataLocation = MDAL_G_dataLocation( g );
EXPECT_EQ( dataLocation, MDAL_DataLocation::DataOnVolumes3D );

ASSERT_EQ( 1, MDAL_G_datasetCount( g ) );
DatasetH ds = MDAL_G_dataset( g, 0 );
ASSERT_NE( ds, nullptr );

double min, max;
MDAL_D_minimumMaximum( ds, &min, &max );
EXPECT_DOUBLE_EQ( 0, min );
EXPECT_DOUBLE_EQ( 0, max );

MDAL_G_minimumMaximum( g, &min, &max );
EXPECT_DOUBLE_EQ( 0, min );
EXPECT_DOUBLE_EQ( 0, max );

double time = MDAL_D_time( ds );
EXPECT_DOUBLE_EQ( 0, time );
}

// Close mesh
MDAL_CloseMesh( m );
}

int main( int argc, char **argv )
Expand Down
2 changes: 1 addition & 1 deletion tests/test_ugrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,7 +324,7 @@ TEST( MeshUgridTest, DFlow12RivierGridClm )
// this dataset contains only discrete set of results
// Mesh2D_s1:flag_values = "1 2 3 4 5 6 7 8 9 10 11 12" ;
MDAL_G_minimumMaximum( g, &min, &max );
EXPECT_DOUBLE_EQ( -2.0422003887246905e-301, min );
EXPECT_DOUBLE_EQ( 2, min );
EXPECT_DOUBLE_EQ( 12, max );

double time = MDAL_D_time( ds );
Expand Down

0 comments on commit 707a0b9

Please sign in to comment.