Skip to content

Commit

Permalink
Update 3Di driver to expose groundwater and surface water as individu…
Browse files Browse the repository at this point in the history 
…al meshes (#477)

* 3Di: expose groundwater and surface water as individual meshes

* Update docstring / cleanup
  • Loading branch information
@uclaros
uclaros authored Mar 20, 2024
1 parent 0cc2799 commit fa0421a
Show file tree
Hide file tree
Showing 7 changed files with 786 additions and 26 deletions.
285 changes: 276 additions & 9 deletions mdal/frmts/mdal_3di.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,11 @@ std::string MDAL::Driver3Di::buildUri( const std::string &meshFile )

populate2DMeshDimensions( dims );
if ( dims.size( CFDimensions::Face ) > 0 )
{
meshNames.push_back( "Mesh2D" );
meshNames.push_back( "Mesh2D_groundwater" );
meshNames.push_back( "Mesh2D_surface_water" );
}

if ( !meshNames.size() )
{
Expand Down Expand Up @@ -107,6 +111,7 @@ MDAL::CFDimensions MDAL::Driver3Di::populateDimensions( )

void MDAL::Driver3Di::populateElements( Vertices &vertices, Edges &edges, Faces &faces )
{
mRequestedMeshFaceIds.clear();
if ( mRequestedMeshName == "Mesh1D" )
populateMesh1DElements( vertices, edges );
else
Expand All @@ -117,11 +122,44 @@ void MDAL::Driver3Di::populateMesh2DElements( MDAL::Vertices &vertices, MDAL::Fa
{
assert( vertices.empty() );
size_t faceCount = mDimensions.size( CFDimensions::Face );
faces.resize( faceCount );
size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
size_t arrsize = faceCount * verticesInFace;
std::map<std::string, size_t> xyToVertex2DId;

if ( mRequestedMeshName == "Mesh2D_groundwater" ||
mRequestedMeshName == "Mesh2D_surface_water" )
{
const int ncidNodeType = mNcFile->getVarId( "Mesh2DNode_type" );
std::vector<int> nodeTypes( arrsize );
if ( nc_get_var_int( mNcFile->handle(), ncidNodeType, nodeTypes.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );

char w_id, wb_id;
if ( mRequestedMeshName == "Mesh2D_groundwater" )
{
if ( nc_get_att( mNcFile->handle(), ncidNodeType, "groundwater_2d", &w_id ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
if ( nc_get_att( mNcFile->handle(), ncidNodeType, "groundwater_boundary_2d", &wb_id ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
}
else if ( mRequestedMeshName == "Mesh2D_surface_water" )
{
if ( nc_get_att( mNcFile->handle(), ncidNodeType, "surface_water_2d", &w_id ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
if ( nc_get_att( mNcFile->handle(), ncidNodeType, "open_water_boundary_2d", &wb_id ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
}

for ( size_t nodeId = 0; nodeId < arrsize; ++nodeId )
{
const int nodeType = nodeTypes.at( nodeId );
if ( nodeType == w_id - '0' || nodeType == wb_id - '0' )
mRequestedMeshFaceIds.push_back( nodeId );
}
faces.reserve( mRequestedMeshFaceIds.size() );
vertices.reserve( mRequestedMeshFaceIds.size() * verticesInFace );
}
else
{
faces.reserve( faceCount );
vertices.reserve( faceCount * verticesInFace );
}

// X coordinate
int ncidX = mNcFile->getVarId( "Mesh2DContour_x" );
double fillX = mNcFile->getFillValue( ncidX );
Expand All @@ -138,6 +176,11 @@ void MDAL::Driver3Di::populateMesh2DElements( MDAL::Vertices &vertices, MDAL::Fa
// are used in multiple faces
for ( size_t faceId = 0; faceId < faceCount; ++faceId )
{
// If requesting a sub-mesh (groundwater / surface water), only use its face ids
if ( !mRequestedMeshFaceIds.empty() &&
std::find( mRequestedMeshFaceIds.cbegin(), mRequestedMeshFaceIds.cend(), faceId ) == mRequestedMeshFaceIds.cend() )
continue;

Face face;

for ( size_t faceVertexId = 0; faceVertexId < verticesInFace; ++faceVertexId )
Expand Down Expand Up @@ -173,7 +216,7 @@ void MDAL::Driver3Di::populateMesh2DElements( MDAL::Vertices &vertices, MDAL::Fa

}

faces[faceId] = face;
faces.push_back( face );
}

// Only now we have number of vertices, since we identified vertices that
Expand All @@ -187,12 +230,12 @@ void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )
if ( 0 == mesh->facesCount() )
return;

size_t faceCount = mesh->facesCount();
const size_t allFaceCount = mDimensions.size( CFDimensions::Face );

// read Z coordinate of 3di computation nodes centers
int ncidZ = mNcFile->getVarId( "Mesh2DFace_zcc" );
double fillZ = mNcFile->getFillValue( ncidZ );
std::vector<double> coordZ( faceCount );
const int ncidZ = mNcFile->getVarId( "Mesh2DFace_zcc" );
const double fillZ = mNcFile->getFillValue( ncidZ );
std::vector<double> coordZ( allFaceCount );
if ( nc_get_var_double( mNcFile->handle(), ncidZ, coordZ.data() ) )
return; //error reading the array

Expand All @@ -209,13 +252,27 @@ void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )

std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
dataset->setTime( MDAL::RelativeTimestamp() );
for ( size_t i = 0; i < faceCount; ++i )

if ( mRequestedMeshFaceIds.empty() )
{
dataset->setScalarValue( i, MDAL::safeValue( coordZ[i], fillZ ) );
for ( size_t i = 0; i < allFaceCount; ++i )
{
dataset->setScalarValue( i, MDAL::safeValue( coordZ[i], fillZ ) );
}
}
else
{
int i = 0;
for ( auto id : mRequestedMeshFaceIds )
{
dataset->setScalarValue( i, MDAL::safeValue( coordZ[id], fillZ ) );
++i;
}
}

dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
group->setStatistics( MDAL::calculateStatistics( group ) );
group->datasets.push_back( dataset );
group->setStatistics( MDAL::calculateStatistics( group ) );
mesh->datasetGroups.push_back( group );
}

Expand All @@ -229,6 +286,27 @@ std::string MDAL::Driver3Di::getTimeVariableName() const
return "time";
}

std::shared_ptr<MDAL::Dataset> MDAL::Driver3Di::create2DDataset( std::shared_ptr<MDAL::DatasetGroup> group, size_t ts, const MDAL::CFDatasetGroupInfo &dsi, double fill_val_x, double fill_val_y )
{
std::shared_ptr<MDAL::CF3DiDataset2D> dataset = std::make_shared<MDAL::CF3DiDataset2D>(
group.get(),
fill_val_x,
fill_val_y,
dsi.ncid_x,
dsi.ncid_y,
dsi.classification_x,
dsi.classification_y,
dsi.timeLocation,
dsi.nTimesteps,
dsi.nValues,
ts,
mNcFile,
mRequestedMeshFaceIds
);
dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
return std::move( dataset );
}

std::set<std::string> MDAL::Driver3Di::ignoreNetCDFVariables()
{
std::set<std::string> ignore_variables;
Expand Down Expand Up @@ -456,3 +534,192 @@ void MDAL::Driver3Di::parse1DConnection( const std::vector<int> &nodesId,
}
}
}

//////////////////////////////////////////////////////////////////////////////////////
MDAL::CF3DiDataset2D::CF3DiDataset2D( MDAL::DatasetGroup *parent,
double fill_val_x,
double fill_val_y,
int ncid_x,
int ncid_y,
Classification classification_x,
Classification classification_y,
CFDatasetGroupInfo::TimeLocation timeLocation,
size_t timesteps,
size_t values,
size_t ts,
std::shared_ptr<NetCDFFile> ncFile,
std::vector<size_t> mask )
: CFDataset2D( parent,
fill_val_x,
fill_val_y,
ncid_x,
ncid_y,
classification_x,
classification_y,
timeLocation,
timesteps,
values,
ts,
ncFile )
, mRequestedMeshFaceIds( mask )
{
}

MDAL::CF3DiDataset2D::~CF3DiDataset2D() = default;

size_t MDAL::CF3DiDataset2D::scalarData( size_t indexStart, size_t count, double *buffer )
{
// Use the basic CF implementation if we don't have sub-mesh face ids
if ( mRequestedMeshFaceIds.empty() )
return CFDataset2D::scalarData( indexStart, count, buffer );

assert( group()->isScalar() ); //checked in C API interface
if ( ( count < 1 ) || ( indexStart >= mRequestedMeshFaceIds.size() ) )
return 0;
if ( mTs >= mTimesteps )
return 0;

// This is the amount of data we should return. It may be limited by the size of mRequestedMeshFaceIds.
const size_t dataCount = indexStart + count < mRequestedMeshFaceIds.size() ? count : mRequestedMeshFaceIds.size() - indexStart;
// These are the first and last data indexes we need to featch from the file.
const size_t requestStart = mRequestedMeshFaceIds[ indexStart ];
const size_t requestEnd = indexStart + count < mRequestedMeshFaceIds.size() ? mRequestedMeshFaceIds[ indexStart + count ] : mRequestedMeshFaceIds.back();
// This is the amound of data we need to fetch from the file.
// It may be larger than count, since we might also fetch values that are not in mRequestedMeshFaceIds
const size_t copyValues = requestEnd - requestStart + 1;
std::vector<double> values_x;

if ( mTimeLocation == CFDatasetGroupInfo::NoTimeDimension )
{
values_x = mNcFile->readDoubleArr(
mNcidX,
requestStart,
copyValues
);
}
else
{
const bool timeFirstDim = mTimeLocation == CFDatasetGroupInfo::TimeDimensionFirst;
const size_t start_dim1 = timeFirstDim ? mTs : requestStart;
const size_t start_dim2 = timeFirstDim ? requestStart : mTs;
const size_t count_dim1 = timeFirstDim ? 1 : copyValues;
const size_t count_dim2 = timeFirstDim ? copyValues : 1;

values_x = mNcFile->readDoubleArr(
mNcidX,
start_dim1,
start_dim2,
count_dim1,
count_dim2
);
}

for ( size_t i = 0; i < dataCount; ++i )
{
const size_t idx = mRequestedMeshFaceIds[indexStart + i] - requestStart;
populate_scalar_vals( buffer,
i,
values_x,
idx,
mFillValX );
}
return dataCount;
}

size_t MDAL::CF3DiDataset2D::vectorData( size_t indexStart, size_t count, double *buffer )
{
// Use the basic CF implementation if we don't have sub-mesh face ids
if ( mRequestedMeshFaceIds.empty() )
return CFDataset2D::vectorData( indexStart, count, buffer );

assert( !group()->isScalar() ); //checked in C API interface
if ( ( count < 1 ) || ( indexStart >= mValues ) )
return 0;

if ( mTs >= mTimesteps )
return 0;

// This is the amount of data we should return. It may be limited by the size of mRequestedMeshFaceIds.
const size_t dataCount = indexStart + count < mRequestedMeshFaceIds.size() ? count : mRequestedMeshFaceIds.size() - indexStart;
// These are the first and last data indexes we need to featch from the file.
const size_t requestStart = mRequestedMeshFaceIds[ indexStart ];
const size_t requestEnd = indexStart + count < mRequestedMeshFaceIds.size() ? mRequestedMeshFaceIds[ indexStart + count ] : mRequestedMeshFaceIds.back();
// This is the amound of data we need to fetch from the file.
// It may be larger than count, since we might also fetch values that are not in mRequestedMeshFaceIds
const size_t copyValues = requestEnd - requestStart + 1;
std::vector<double> values_x;
std::vector<double> values_y;

if ( mTimeLocation == CFDatasetGroupInfo::NoTimeDimension )
{
values_x = mNcFile->readDoubleArr(
mNcidX,
requestStart,
copyValues
);

values_y = mNcFile->readDoubleArr(
mNcidX,
requestStart,
copyValues
);
}
else
{
const bool timeFirstDim = mTimeLocation == CFDatasetGroupInfo::TimeDimensionFirst;
const size_t start_dim1 = timeFirstDim ? mTs : requestStart;
const size_t start_dim2 = timeFirstDim ? requestStart : mTs;
const size_t count_dim1 = timeFirstDim ? 1 : copyValues;
const size_t count_dim2 = timeFirstDim ? copyValues : 1;

values_x = mNcFile->readDoubleArr(
mNcidX,
start_dim1,
start_dim2,
count_dim1,
count_dim2
);
values_y = mNcFile->readDoubleArr(
mNcidY,
start_dim1,
start_dim2,
count_dim1,
count_dim2
);
}

//if values component are classified convert from index to value
if ( !mClassificationX.empty() )
{
fromClassificationToValue( mClassificationX, values_x, 1 );
}

if ( !mClassificationY.empty() )
{
fromClassificationToValue( mClassificationY, values_y, 1 );
}

for ( size_t i = 0; i < dataCount; ++i )
{
const size_t idx = mRequestedMeshFaceIds[indexStart + i] - requestStart;
if ( group()->isPolar() )
populate_polar_vector_vals( buffer,
i,
values_x,
values_y,
idx,
mFillValX,
mFillValY,
group()->referenceAngles() );
else
populate_vector_vals( buffer,
i,
values_x,
values_y,
idx,
mFillValX,
mFillValY );
}

return dataCount;
}
Loading

0 comments on commit fa0421a

Please sign in to comment.