Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update 3Di driver to expose groundwater and surface water as individual meshes #477

Merged
merged 2 commits into from
Mar 20, 2024
+786 −26
Merged

Update 3Di driver to expose groundwater and surface water as individual meshes #477

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c0ff5c3
3Di: expose groundwater and surface water as individual meshes
uclaros Mar 7, 2024
04a97c6
Update docstring / cleanup
uclaros Mar 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
285 changes: 276 additions & 9 deletions mdal/frmts/mdal_3di.cpp
View file
Open in desktop
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
Loading